ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLCrossSectionsMultiPionsAndResonances.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLCrossSectionsMultiPionsAndResonances.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
45 #include "G4INCLCrossSectionsMultiPionsAndResonances.hh"
46 #include "G4INCLKinematicsUtils.hh"
47 #include "G4INCLParticleTable.hh"
48 // #include <cassert>
49 
50 namespace G4INCL {
51 
52  template<G4int N>
53  struct BystrickyEvaluator {
54  static G4double eval(const G4double pLab, const G4double oneOverThreshold, HornerCoefficients<N> const &coeffs) {
55  const G4double pMeV = pLab*1E3;
56  const G4double ekin=std::sqrt(ParticleTable::effectiveNucleonMass2+pMeV*pMeV)-ParticleTable::effectiveNucleonMass;
57  const G4double xrat=ekin*oneOverThreshold;
58  const G4double x=std::log(xrat);
59  return HornerEvaluator<N>::eval(x, coeffs) * x * std::exp(-0.5*x);
60  }
61  };
62 
63  const G4int CrossSectionsMultiPionsAndResonances::nMaxPiNN = 4;
64  const G4int CrossSectionsMultiPionsAndResonances::nMaxPiPiN = 4;
65 
66  const G4double CrossSectionsMultiPionsAndResonances::s11pzOOT = 0.0035761542037692665889;
67  const G4double CrossSectionsMultiPionsAndResonances::s01ppOOT = 0.003421025623481919853;
68  const G4double CrossSectionsMultiPionsAndResonances::s01pzOOT = 0.0035739814152966403123;
69  const G4double CrossSectionsMultiPionsAndResonances::s11pmOOT = 0.0034855350296270480281;
70  const G4double CrossSectionsMultiPionsAndResonances::s12pmOOT = 0.0016672224074691565119;
71  const G4double CrossSectionsMultiPionsAndResonances::s12ppOOT = 0.0016507643038726931312;
72  const G4double CrossSectionsMultiPionsAndResonances::s12zzOOT = 0.0011111111111111111111;
73  const G4double CrossSectionsMultiPionsAndResonances::s02pzOOT = 0.00125;
74  const G4double CrossSectionsMultiPionsAndResonances::s02pmOOT = 0.0016661112962345883443;
75  const G4double CrossSectionsMultiPionsAndResonances::s12mzOOT = 0.0017047391749062392793;
76 
77  CrossSectionsMultiPionsAndResonances::CrossSectionsMultiPionsAndResonances() :
78  s11pzHC(-2.228000000000294018,8.7560000000005723725,-0.61000000000023239325,-5.4139999999999780324,3.3338333333333348023,-0.75835000000000022049,0.060623611111111114688),
79  s01ppHC(2.0570000000126518344,-6.029000000012135826,36.768500000002462784,-45.275666666666553533,25.112666666666611953,-7.2174166666666639187,1.0478875000000000275,-0.060804365079365080846),
80  s01pzHC(0.18030000000000441851,7.8700999999999953598,-4.0548999999999990425,0.555199999999999959),
81  s11pmHC(0.20590000000000031866,3.3450999999999993936,-1.4401999999999997825,0.17076666666666664973),
82  s12pmHC(-0.77235999999999901328,4.2626599999999991117,-1.9008899999999997323,0.30192266666666663379,-0.012270833333333331986),
83  s12ppHC(-0.75724999999999975664,2.0934399999999998565,-0.3803099999999999814),
84  s12zzHC(-0.89599999999996965072,7.882999999999978632,-7.1049999999999961928,1.884333333333333089),
85  s02pzHC(-1.0579999999999967036,11.113999999999994089,-8.5259999999999990196,2.0051666666666666525),
86  s02pmHC(2.4009000000012553286,-7.7680000000013376183,20.619000000000433505,-16.429666666666723928,5.2525708333333363472,-0.58969166666666670206),
87  s12mzHC(-0.21858699999999976269,1.9148999999999999722,-0.31727500000000001065,-0.027695000000000000486)
88  {
89  }
90 
91  G4double CrossSectionsMultiPionsAndResonances::total(Particle const * const p1, Particle const * const p2) {
92  G4double inelastic;
93  if(p1->isNucleon() && p2->isNucleon()) {
94  return CrossSectionsMultiPions::NNTot(p1, p2);
95  } else if((p1->isNucleon() && p2->isDelta()) ||
96  (p1->isDelta() && p2->isNucleon())) {
97  inelastic = CrossSectionsMultiPions::NDeltaToNN(p1, p2);
98  } else if((p1->isNucleon() && p2->isPion()) ||
99  (p1->isPion() && p2->isNucleon())) {
100  return CrossSectionsMultiPions::piNTot(p1,p2);
101  } else if((p1->isNucleon() && p2->isEta()) ||
102  (p1->isEta() && p2->isNucleon())) {
103  inelastic = etaNToPiN(p1,p2) + etaNToPiPiN(p1,p2);
104  } else if((p1->isNucleon() && p2->isOmega()) ||
105  (p1->isOmega() && p2->isNucleon())) {
106  inelastic = omegaNInelastic(p1,p2);
107  } else if((p1->isNucleon() && p2->isEtaPrime()) ||
108  (p1->isEtaPrime() && p2->isNucleon())) {
109  inelastic = etaPrimeNToPiN(p1,p2);
110  } else {
111  inelastic = 0.;
112  }
113 
114  return inelastic + elastic(p1, p2);
115  }
116 
117 
118  G4double CrossSectionsMultiPionsAndResonances::elastic(Particle const * const p1, Particle const * const p2) {
119  if((p1->isNucleon()||p1->isDelta()) && (p2->isNucleon()||p2->isDelta())){
120  return CrossSectionsMultiPions::elastic(p1, p2);
121  }
122  else if ((p1->isNucleon() && p2->isPion()) || (p2->isNucleon() && p1->isPion())){
123  return CrossSectionsMultiPions::elastic(p1, p2);
124  }
125  else if ((p1->isNucleon() && p2->isEta()) || (p2->isNucleon() && p1->isEta())){
126  return etaNElastic(p1, p2);
127  }
128  else if ((p1->isNucleon() && p2->isOmega()) || (p2->isNucleon() && p1->isOmega())){
129  return omegaNElastic(p1, p2);
130  }
131  else {
132  return 0.0;
133  }
134  }
135 
136 
137  G4double CrossSectionsMultiPionsAndResonances::piNToxPiN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
138  //
139  // pion-Nucleon producing xpi pions cross sections (corrected due to eta and omega)
140  //
141 // assert(xpi>1 && xpi<=nMaxPiPiN);
142 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
143 
144  const G4double oldXS2Pi=CrossSectionsMultiPions::piNToxPiN(2,particle1, particle2);
145  const G4double oldXS3Pi=CrossSectionsMultiPions::piNToxPiN(3,particle1, particle2);
146  const G4double oldXS4Pi=CrossSectionsMultiPions::piNToxPiN(4,particle1, particle2);
147  const G4double xsEta=piNToEtaN(particle1, particle2);
148  const G4double xsOmega=piNToOmegaN(particle1, particle2);
149  G4double newXS2Pi=0.;
150  G4double newXS3Pi=0.;
151  G4double newXS4Pi=0.;
152 
153  if (xpi == 2) {
154  if (oldXS4Pi != 0.)
155  newXS2Pi=oldXS2Pi;
156  else if (oldXS3Pi != 0.) {
157  newXS3Pi=oldXS3Pi-xsEta-xsOmega;
158  if (newXS3Pi < 1.e-09)
159  newXS2Pi=oldXS2Pi-(xsEta+xsOmega-oldXS3Pi);
160  else
161  newXS2Pi=oldXS2Pi;
162  }
163  else {
164  newXS2Pi=oldXS2Pi-xsEta-xsOmega;
165  if (newXS2Pi < 1.e-09)
166  newXS2Pi=0.;
167  }
168  return newXS2Pi;
169  }
170  else if (xpi == 3) {
171  if (oldXS4Pi != 0.) {
172  newXS4Pi=oldXS4Pi-xsEta-xsOmega;
173  if (newXS4Pi < 1.e-09)
174  newXS3Pi=oldXS3Pi-(xsEta+xsOmega-oldXS4Pi);
175  else
176  newXS3Pi=oldXS3Pi;
177  }
178  else {
179  newXS3Pi=oldXS3Pi-xsEta-xsOmega;
180  if (newXS3Pi < 1.e-09)
181  newXS3Pi=0.;
182  }
183  return newXS3Pi;
184  }
185  else if (xpi == 4) {
186  newXS4Pi=oldXS4Pi-xsEta-xsOmega;
187  if (newXS4Pi < 1.e-09)
188  newXS4Pi=0.;
189  return newXS4Pi;
190  }
191  else // should never reach this point
192  return 0.;
193  }
194 
195  G4double CrossSectionsMultiPionsAndResonances::piNToEtaN(Particle const * const particle1, Particle const * const particle2) {
196  //
197  // Pion-Nucleon producing Eta cross sections
198  //
199 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
200 
201  G4double sigma;
202  sigma=piMinuspToEtaN(particle1,particle2);
203 
204  const G4int isoin = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
205 
206  if (isoin == -1) {
207  if ((particle1->getType()) == Proton || (particle2->getType()) == Proton) return sigma;
208  else return 0.5 * sigma;
209  }
210  else if (isoin == 1) {
211  if ((particle1->getType()) == Neutron || (particle2->getType()) == Neutron) return sigma;
212  else return 0.5 * sigma;
213  }
214  else return 0. ; // should never return 0. (?) // pi+ p and pi- n return 0.
215 
216 // return sigma;
217  }
218 
219  G4double CrossSectionsMultiPionsAndResonances::piNToOmegaN(Particle const * const particle1, Particle const * const particle2) {
220  //
221  // Pion-Nucleon producing Omega cross sections
222  //
223 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
224 
225  G4double sigma;
226  sigma=piMinuspToOmegaN(particle1,particle2);
227 
228  const G4int isoin = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
229 
230  if (isoin == -1) {
231  if ((particle1->getType()) == Proton || (particle2->getType()) == Proton) return sigma;
232  else return 0.5 * sigma;
233  }
234  else if (isoin == 1) {
235  if ((particle1->getType()) == Neutron || (particle2->getType()) == Neutron) return sigma;
236  else return 0.5 * sigma;
237  }
238  else return 0. ; // should never return 0. (?) // pi+ p and pi- n return 0.
239 
240 // return sigma;
241  }
242 
243 #if defined(NDEBUG) || defined(INCLXX_IN_GEANT4_MODE)
244  G4double CrossSectionsMultiPionsAndResonances::piNToEtaPrimeN(Particle const * const /*particle1*/, Particle const * const /*particle2*/) {
245 #else
246  G4double CrossSectionsMultiPionsAndResonances::piNToEtaPrimeN(Particle const * const particle1, Particle const * const particle2) {
247 #endif
248  //
249  // Pion-Nucleon producing EtaPrime cross sections
250  //
251 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
252 
253  return 0.;
254  }
255 
256  G4double CrossSectionsMultiPionsAndResonances::etaNToPiN(Particle const * const particle1, Particle const * const particle2) {
257  //
258  // Eta-Nucleon producing Pion cross sections
259  //
260 // assert((particle1->isNucleon() && particle2->isEta()) || (particle1->isEta() && particle2->isNucleon()));
261 
262  const Particle *eta;
263  const Particle *nucleon;
264 
265  if (particle1->isEta()) {
266  eta = particle1;
267  nucleon = particle2;
268  }
269  else {
270  eta = particle2;
271  nucleon = particle1;
272  }
273 
274  const G4double pLab = KinematicsUtils::momentumInLab(eta, nucleon);
275  G4double sigma=0.;
276 
277  if (pLab <= 574.)
278  sigma= 1.511147E-13*std::pow(pLab,6)- 3.603636E-10*std::pow(pLab,5)+ 3.443487E-07*std::pow(pLab,4)- 1.681980E-04*std::pow(pLab,3)+ 4.437913E-02*std::pow(pLab,2)- 6.172108E+00*pLab+ 4.031449E+02;
279  else if (pLab <= 850.)
280  sigma= -8.00018E-14*std::pow(pLab,6)+ 3.50041E-10*std::pow(pLab,5)- 6.33891E-07*std::pow(pLab,4)+ 6.07658E-04*std::pow(pLab,3)- 3.24936E-01*std::pow(pLab,2)+ 9.18098E+01*pLab- 1.06943E+04;
281  else if (pLab <= 1300.)
282  sigma= 6.56364E-09*std::pow(pLab,3)- 2.07653E-05*std::pow(pLab,2)+ 1.84148E-02*pLab- 1.70427E+00;
283  else {
284  G4double ECM=KinematicsUtils::totalEnergyInCM(eta, nucleon);
285  G4double massPiZero=ParticleTable::getINCLMass(PiZero);
286  G4double massPiMinus=ParticleTable::getINCLMass(PiMinus);
287  G4double massProton=ParticleTable::getINCLMass(Proton);
288  G4double masseta;
289  G4double massnucleon;
290  G4double pCM_eta;
291  masseta=eta->getMass();
292  massnucleon=nucleon->getMass();
293  pCM_eta=KinematicsUtils::momentumInCM(ECM, masseta, massnucleon);
294  G4double pCM_PiZero=KinematicsUtils::momentumInCM(ECM, massPiZero, massProton);
295  G4double pCM_PiMinus=KinematicsUtils::momentumInCM(ECM, massPiMinus, massProton); // = pCM_PiPlus (because massPiMinus = massPiPlus)
296  sigma = (piMinuspToEtaN(ECM)/2.) * std::pow((pCM_PiZero/pCM_eta), 2) + piMinuspToEtaN(ECM) * std::pow((pCM_PiMinus/pCM_eta), 2);
297  }
298  if (sigma < 0.) sigma=0.;
299 
300  return sigma;
301  }
302 
303  G4double CrossSectionsMultiPionsAndResonances::etaNToPiPiN(Particle const * const particle1, Particle const * const particle2) {
304  //
305  // Eta-Nucleon producing Two Pions cross sections
306  //
307 // assert((particle1->isNucleon() && particle2->isEta()) || (particle1->isEta() && particle2->isNucleon()));
308 
309  G4double sigma=0.;
310 
311  const Particle *eta;
312  const Particle *nucleon;
313 
314  if (particle1->isEta()) {
315  eta = particle1;
316  nucleon = particle2;
317  }
318  else {
319  eta = particle2;
320  nucleon = particle1;
321  }
322 
323  const G4double pLab = KinematicsUtils::momentumInLab(eta, nucleon);
324 
325  if (pLab < 450.)
326  sigma = 2.01854221E-13*std::pow(pLab,6) - 3.49750459E-10*std::pow(pLab,5) + 2.46011585E-07*std::pow(pLab,4) - 9.01422901E-05*std::pow(pLab,3) + 1.83382964E-02*std::pow(pLab,2) - 2.03113098E+00*pLab + 1.10358550E+02;
327  else if (pLab < 600.)
328  sigma = 2.01854221E-13*std::pow(450.,6) - 3.49750459E-10*std::pow(450.,5) + 2.46011585E-07*std::pow(450.,4) - 9.01422901E-05*std::pow(450.,3) + 1.83382964E-02*std::pow(450.,2) - 2.03113098E+00*450. + 1.10358550E+02;
329  else if (pLab <= 1300.)
330  sigma = -6.32793049e-16*std::pow(pLab,6) + 3.95985900e-12*std::pow(pLab,5) - 1.01727714e-8*std::pow(pLab,4) +
331  1.37055547e-05*std::pow(pLab,3) - 1.01830486e-02*std::pow(pLab,2) + 3.93492126*pLab - 609.447145;
332  else
333  sigma = etaNToPiN(particle1,particle2);
334 
335  if (sigma < 0.) sigma = 0.;
336  return sigma; // Parameterization from the ANL-Osaka DCC model [PRC88(2013)035209] - eta p --> "pi+pi0 n" + "pi0 pi0 p" total XS
337  }
338 
339 
340  G4double CrossSectionsMultiPionsAndResonances::etaNElastic(Particle const * const particle1, Particle const * const particle2) {
341  //
342  // Eta-Nucleon elastic cross sections
343  //
344 // assert((particle1->isNucleon() && particle2->isEta()) || (particle1->isEta() && particle2->isNucleon()));
345 
346  G4double sigma=0.;
347 
348  const Particle *eta;
349  const Particle *nucleon;
350 
351  if (particle1->isEta()) {
352  eta = particle1;
353  nucleon = particle2;
354  }
355  else {
356  eta = particle2;
357  nucleon = particle1;
358  }
359 
360  const G4double pLab = KinematicsUtils::momentumInLab(eta, nucleon);
361 
362  if (pLab < 700.)
363  sigma = 3.6838e-15*std::pow(pLab,6) - 9.7815e-12*std::pow(pLab,5) + 9.7914e-9*std::pow(pLab,4) -
364  4.3222e-06*std::pow(pLab,3) + 7.9188e-04*std::pow(pLab,2) - 1.8379e-01*pLab + 84.965;
365  else if (pLab < 1400.)
366  sigma = 3.562630e-16*std::pow(pLab,6) - 2.384766e-12*std::pow(pLab,5) + 6.601312e-9*std::pow(pLab,4) -
367  9.667078e-06*std::pow(pLab,3) + 7.894845e-03*std::pow(pLab,2) - 3.409200*pLab + 609.8501;
368  else if (pLab < 2025.)
369  sigma = -1.041950E-03*pLab + 2.110529E+00;
370  else
371  sigma=0.;
372 
373  if (sigma < 0.) sigma = 0.;
374  return sigma; // Parameterization from the ANL-Osaka DCC model [PRC88(2013)035209]
375  }
376 
377  G4double CrossSectionsMultiPionsAndResonances::omegaNInelastic(Particle const * const particle1, Particle const * const particle2) {
378  //
379  // Omega-Nucleon inelastic cross sections
380  //
381 // assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
382 
383  G4double sigma=0.;
384 
385  const Particle *omega;
386  const Particle *nucleon;
387 
388  if (particle1->isOmega()) {
389  omega = particle1;
390  nucleon = particle2;
391  }
392  else {
393  omega = particle2;
394  nucleon = particle1;
395  }
396 
397  const G4double pLab = KinematicsUtils::momentumInLab(omega, nucleon)/1000.; // GeV/c
398 
399  sigma = 20. + 4.0/pLab; // Eq.(24) in G.I. Lykasov et al., EPJA 6, 71-81 (1999)
400 
401  return sigma;
402  }
403 
404 
405  G4double CrossSectionsMultiPionsAndResonances::omegaNElastic(Particle const * const particle1, Particle const * const particle2) {
406  //
407  // Omega-Nucleon elastic cross sections
408  //
409 // assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
410 
411  G4double sigma=0.;
412 
413  const Particle *omega;
414  const Particle *nucleon;
415 
416  if (particle1->isOmega()) {
417  omega = particle1;
418  nucleon = particle2;
419  }
420  else {
421  omega = particle2;
422  nucleon = particle1;
423  }
424 
425  const G4double pLab = KinematicsUtils::momentumInLab(omega, nucleon)/1000.; // GeV/c
426 
427  sigma = 5.4 + 10.*std::exp(-0.6*pLab); // Eq.(21) in G.I. Lykasov et al., EPJA 6, 71-81 (1999)
428 
429  return sigma;
430  }
431 
432 
433  G4double CrossSectionsMultiPionsAndResonances::omegaNToPiN(Particle const * const particle1, Particle const * const particle2) {
434  //
435  // Omega-Nucleon producing Pion cross sections
436  //
437 // assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
438 
439  G4double ECM=KinematicsUtils::totalEnergyInCM(particle1, particle2);
440 
441  G4double massPiZero=ParticleTable::getINCLMass(PiZero);
442  G4double massPiMinus=ParticleTable::getINCLMass(PiMinus);
443  G4double massProton=ParticleTable::getINCLMass(Proton);
444 
445  G4double massomega;
446  G4double massnucleon;
447  G4double pCM_omega;
448  G4double pLab_omega;
449 
450  G4double sigma=0.;
451 
452  if (particle1->isOmega()) {
453  massomega=particle1->getMass();
454  massnucleon=particle2->getMass();
455  }
456  else {
457  massomega=particle2->getMass();
458  massnucleon=particle1->getMass();
459  }
460  pCM_omega=KinematicsUtils::momentumInCM(ECM, massomega, massnucleon);
461  pLab_omega=KinematicsUtils::momentumInLab(ECM*ECM, massomega, massnucleon);
462 
463  G4double pCM_PiZero=KinematicsUtils::momentumInCM(ECM, massPiZero, massProton);
464  G4double pCM_PiMinus=KinematicsUtils::momentumInCM(ECM, massPiMinus, massProton); // = pCM_PiPlus (because massPiMinus = massPiPlus)
465 
466  sigma = (piMinuspToOmegaN(ECM)/2.) * std::pow((pCM_PiZero/pCM_omega), 2)
467  + piMinuspToOmegaN(ECM) * std::pow((pCM_PiMinus/pCM_omega), 2);
468 
469  if (sigma > omegaNInelastic(particle1, particle2) || (pLab_omega < 200.)) {
470 // if (sigma > omegaNInelastic(particle1, particle2)) {
471  sigma = omegaNInelastic(particle1, particle2);
472  }
473 
474  return sigma;
475  }
476 
477 
478  G4double CrossSectionsMultiPionsAndResonances::omegaNToPiPiN(Particle const * const particle1, Particle const * const particle2) {
479  //
480  // Omega-Nucleon producing 2 PionS cross sections
481  //
482 // assert((particle1->isNucleon() && particle2->isOmega()) || (particle1->isOmega() && particle2->isNucleon()));
483 
484  G4double sigma=0.;
485 
486  sigma = omegaNInelastic(particle1,particle2) - omegaNToPiN(particle1,particle2) ;
487 
488  return sigma;
489  }
490 
491 
492 #if defined(NDEBUG) || defined(INCLXX_IN_GEANT4_MODE)
493  G4double CrossSectionsMultiPionsAndResonances::etaPrimeNToPiN(Particle const * const /*particle1*/, Particle const * const /*particle2*/) {
494 #else
495  G4double CrossSectionsMultiPionsAndResonances::etaPrimeNToPiN(Particle const * const particle1, Particle const * const particle2) {
496 #endif
497  //
498  // EtaPrime-Nucleon producing Pion cross sections
499  //
500 // assert((particle1->isNucleon() && particle2->isEtaPrime()) || (particle1->isEtaPrime() && particle2->isNucleon()));
501 
502  return 0.;
503  }
504 
505  G4double CrossSectionsMultiPionsAndResonances::piMinuspToEtaN(Particle const * const particle1, Particle const * const particle2) {
506  //
507  // Pion-Nucleon producing Eta cross sections
508  //
509 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
510 
511  G4double masspion;
512  G4double massnucleon;
513  if (particle1->isPion()) {
514  masspion=particle1->getMass();
515  massnucleon=particle2->getMass();
516  }
517  else {
518  masspion=particle2->getMass();
519  massnucleon=particle1->getMass();
520  }
521 
522  G4double ECM=KinematicsUtils::totalEnergyInCM(particle1, particle2);
523  G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.; // GeV/c
524 
525  G4double sigma;
526 
527 // new parameterization (JCD) - end of february 2016
528  if (ECM < 1486.5) sigma=0.;
529  else
530  {
531  if (ECM < 1535.)
532  {
533  sigma = -0.0000003689197974814*std::pow(ECM,4) + 0.002260193900097*std::pow(ECM,3) - 5.193105877187*std::pow(ECM,2) + 5303.505273919*ECM - 2031265.900648;
534  }
535  else if (ECM < 1670.)
536  {
537  sigma = -0.0000000337986446*std::pow(ECM,4) + 0.000218279989*std::pow(ECM,3) - 0.528276144*std::pow(ECM,2) + 567.828367*ECM - 228709.42;
538  }
539  else if (ECM < 1714.)
540  {
541  sigma = 0.000003737765*std::pow(ECM,2) - 0.005664062*ECM;
542  }
543  else sigma=1.47*std::pow(plab, -1.68);
544  }
545 //
546 
547  return sigma;
548  }
549 
550  G4double CrossSectionsMultiPionsAndResonances::piMinuspToEtaN(const G4double ECM) {
551  //
552  // Pion-Nucleon producing Eta cross sections
553  //
554 
555  const G4double masspion = ParticleTable::getRealMass(PiMinus);
556  const G4double massnucleon = ParticleTable::getRealMass(Proton);
557 
558  G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.; // GeV/c
559 
560  G4double sigma;
561 
562 // new parameterization (JCD) - end of february 2016
563  if (ECM < 1486.5) sigma=0.;
564  else
565  {
566  if (ECM < 1535.)
567  {
568  sigma = -0.0000003689197974814*std::pow(ECM,4) + 0.002260193900097*std::pow(ECM,3) - 5.193105877187*std::pow(ECM,2) + 5303.505273919*ECM - 2031265.900648;
569  }
570  else if (ECM < 1670.)
571  {
572  sigma = -0.0000000337986446*std::pow(ECM,4) + 0.000218279989*std::pow(ECM,3) - 0.528276144*std::pow(ECM,2) + 567.828367*ECM - 228709.42;
573  }
574  else if (ECM < 1714.)
575  {
576  sigma = 0.000003737765*std::pow(ECM,2) - 0.005664062*ECM;
577  }
578  else sigma=1.47*std::pow(plab, -1.68);
579  }
580 
581  return sigma;
582  }
583 
584  G4double CrossSectionsMultiPionsAndResonances::piMinuspToOmegaN(Particle const * const particle1, Particle const * const particle2) {
585  //
586  // Pion-Nucleon producing Omega cross sections
587  //
588 // assert((particle1->isNucleon() && particle2->isPion()) || (particle1->isPion() && particle2->isNucleon()));
589 //jcd to be removed
590 // return 0.;
591 //jcd
592 
593 // G4double param=1.095 ; // Deneye (Thesis)
594  G4double param=1.0903 ; // JCD (threshold taken into account)
595 
596  G4double masspion;
597  G4double massnucleon;
598  if (particle1->isPion()) {
599  masspion=particle1->getMass();
600  massnucleon=particle2->getMass();
601  }
602  else {
603  masspion=particle2->getMass();
604  massnucleon=particle1->getMass();
605  }
606  G4double ECM=KinematicsUtils::totalEnergyInCM(particle1, particle2);
607  G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.; // GeV/c
608 
609  G4double sigma;
610  if (plab < param) sigma=0.;
611  else sigma=13.76*(plab-param)/(std::pow(plab, 3.33) - 1.07); // Phys. Rev. C 41, 1701–1718 (1990)
612 
613  return sigma;
614 }
615  G4double CrossSectionsMultiPionsAndResonances::piMinuspToOmegaN(const G4double ECM) {
616  //
617  // Pion-Nucleon producing Omega cross sections
618  //
619 //jcd to be removed
620 // return 0.;
621 //jcd
622 
623 // G4double param=1.095 ; // Deneye (Thesis)
624  G4double param=1.0903 ; // JCD (threshold taken into account)
625 
626  const G4double masspion = ParticleTable::getRealMass(PiMinus);
627  const G4double massnucleon = ParticleTable::getRealMass(Proton);
628 
629  G4double plab=KinematicsUtils::momentumInLab(ECM*ECM, masspion, massnucleon)/1000.; // GeV/c
630 
631  G4double sigma;
632  if (plab < param) sigma=0.;
633  else sigma=13.76*(plab-param)/(std::pow(plab, 3.33)-1.07);
634 
635  return sigma;
636  }
637 
638  G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaIso(const G4double ener, const G4int iso) {
639 
640  const G4double Ecm=0.001*ener;
641  G4double sNNEta; // pp->pp+eta(+X)
642  G4double sNNEta1; // np->np+eta(+X)
643  G4double sNNEta2; // np->d+eta (d will be considered as np - How far is this right?)
644  G4double x=Ecm*Ecm/5.88;
645 
646 //jcd
647  if (Ecm >= 3.05) {
648  sNNEta = 2.5*std::pow((x-1.),1.47)*std::pow(x,-1.25)*1000.;
649  }
650  else if(Ecm >= 2.6) {
651  sNNEta = -327.29*Ecm*Ecm*Ecm + 2870.*Ecm*Ecm - 7229.3*Ecm + 5273.3;
652  if (sNNEta <= NNToNNEtaExcluIso(ener, 2)*1000.) sNNEta = NNToNNEtaExcluIso(ener, 2)*1000.;
653  }
654  else {
655  sNNEta = NNToNNEtaExcluIso(ener, 2)*1000.;
656  }
657 //jcd
658  if (sNNEta < 1.e-9) sNNEta = 0.;
659 
660  if (iso != 0) {
661  return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
662  }
663 
664  if(Ecm >= 6.25) {
665  sNNEta1 = sNNEta;
666  }
667  else if (Ecm >= 2.6) {
668  sNNEta1 = sNNEta*std::exp(-(-5.53151576/Ecm+0.8850425));
669  }
670  else if (Ecm >= 2.525) { // = exclusive pn
671  sNNEta1 = -4433.586*Ecm*Ecm*Ecm*Ecm + 56581.54*Ecm*Ecm*Ecm - 270212.6*Ecm*Ecm + 571650.6*Ecm - 451091.6;
672  }
673  else { // = exclusive pn
674  sNNEta1 = 17570.217219*Ecm*Ecm - 84910.985402*Ecm + 102585.55847;
675  }
676 
677  sNNEta2 = -10220.89518466*Ecm*Ecm+51227.30841724*Ecm-64097.96025731;
678  if (sNNEta2 < 0.) sNNEta2=0.;
679 
680  sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;
681 
682  G4double Mn=ParticleTable::getRealMass(Neutron)/1000.;
683  G4double Mp=ParticleTable::getRealMass(Proton)/1000.;
684  G4double Meta=ParticleTable::getRealMass(Eta)/1000.;
685  if (sNNEta < 1.e-9 || Ecm < Mn+Mp+Meta) sNNEta = 0.;
686 
687  return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
688  }
689 
690 
691  G4double CrossSectionsMultiPionsAndResonances::NNToNNEta(Particle const * const particle1, Particle const * const particle2) {
692 
693  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
694  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
695 
696  if (iso != 0) {
697  return NNToNNEtaIso(ener, iso);
698  }
699  else {
700  return 0.5*(NNToNNEtaIso(ener, 0)+NNToNNEtaIso(ener, 2));
701  }
702  }
703 
704  G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaExcluIso(const G4double ener, const G4int iso) {
705 
706  const G4double Ecm=0.001*ener;
707  G4double sNNEta; // pp->pp+eta
708  G4double sNNEta1; // np->np+eta
709  G4double sNNEta2; // np->d+eta (d wil be considered as np - How far is this right?)
710 
711  if(Ecm>=3.875) { // By hand (JCD)
712  sNNEta = -13.008*Ecm*Ecm + 84.531*Ecm + 36.234;
713  }
714  else if(Ecm>=2.725) { // By hand (JCD)
715  sNNEta = -913.2809*std::pow(Ecm,5) + 15564.27*std::pow(Ecm,4) - 105054.9*std::pow(Ecm,3) + 351294.2*std::pow(Ecm,2) - 582413.9*Ecm + 383474.7;
716  }
717  else if(Ecm>=2.575) { // By hand (JCD)
718  sNNEta = -2640.3*Ecm*Ecm + 14692*Ecm - 20225;
719  }
720  else {
721  sNNEta = -147043.497285*std::pow(Ecm,4) + 1487222.5438123*std::pow(Ecm,3) - 5634399.900744*std::pow(Ecm,2) + 9477290.199378*Ecm - 5972174.353438;
722  }
723 
724  G4double Mn=ParticleTable::getRealMass(Neutron)/1000.;
725  G4double Mp=ParticleTable::getRealMass(Proton)/1000.;
726  G4double Meta=ParticleTable::getRealMass(Eta)/1000.;
727  G4double Thr0=0.;
728  if (iso > 0) {
729  Thr0=2.*Mp+Meta;
730  }
731  else if (iso < 0) {
732  Thr0=2.*Mn+Meta;
733  }
734  else {
735  Thr0=Mn+Mp+Meta;
736  }
737 
738  if (sNNEta < 1.e-9 || Ecm < Thr0) sNNEta = 0.; // Thr0: Ecm threshold
739 
740  if (iso != 0) {
741  return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
742  }
743 
744  if(Ecm>=3.9) {
745  sNNEta1 = sNNEta;
746  }
747  else if (Ecm >= 3.5) {
748  sNNEta1 = -1916.2*Ecm*Ecm*Ecm + 21556.0*Ecm*Ecm - 80828.0*Ecm + 101200.0;
749  }
750  else if (Ecm >= 2.525) {
751  sNNEta1 = -4433.586*Ecm*Ecm*Ecm*Ecm + 56581.54*Ecm*Ecm*Ecm - 270212.6*Ecm*Ecm + 571650.6*Ecm - 451091.6;
752  }
753  else {
754  sNNEta1 = 17570.217219*Ecm*Ecm - 84910.985402*Ecm + 102585.55847;
755  }
756 
757  sNNEta2 = -10220.89518466*Ecm*Ecm+51227.30841724*Ecm-64097.96025731;
758  if (sNNEta2 < 0.) sNNEta2=0.;
759 
760  sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;
761  if (sNNEta < 1.e-9 || Ecm < Thr0) sNNEta = 0.; // Thr0: Ecm threshold
762 
763  return sNNEta/1000.; // parameterization in microbarn (not millibarn)!
764 
765  }
766 
767  G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaExclu(Particle const * const particle1, Particle const * const particle2) {
768 
769  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
770  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
771 
772  if (iso != 0) {
773  return NNToNNEtaExcluIso(ener, iso);
774  }
775  else {
776  return 0.5*(NNToNNEtaExcluIso(ener, 0)+NNToNNEtaExcluIso(ener, 2));
777  }
778  }
779 
780 
781  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaIso(const G4double ener, const G4int iso) {
782 
783  const G4double Ecm=0.001*ener;
784  G4double sNNOmega; // pp->pp+eta(+X)
785  G4double sNNOmega1; // np->np+eta(+X)
786  G4double x=Ecm*Ecm/7.06;
787 
788  if(Ecm>4.0) {
789  sNNOmega = 2.5*std::pow(x-1, 1.47)*std::pow(x, -1.11) ;
790  }
791  else if(Ecm>2.802) { // 2802 MeV -> threshold to open inclusive (based on multipion threshold and omega mass)
792  sNNOmega = (568.5254*Ecm*Ecm - 2694.045*Ecm + 3106.247)/1000.;
793  if (sNNOmega <= NNToNNOmegaExcluIso(ener, 2)) sNNOmega = NNToNNOmegaExcluIso(ener, 2);
794  }
795  else {
796  sNNOmega = NNToNNOmegaExcluIso(ener, 2);
797  }
798 
799  if (sNNOmega < 1.e-9) sNNOmega = 0.;
800 
801  if (iso != 0) {
802  return sNNOmega;
803  }
804 
805  sNNOmega1 = 3.*sNNOmega; // 3.0: ratio pn/pp (5 from theory; 2 from experiments)
806 
807  sNNOmega = 2.*sNNOmega1-sNNOmega;
808 
809  if (sNNOmega < 1.e-9) sNNOmega = 0.;
810 
811  return sNNOmega;
812  }
813 
814 
815  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmega(Particle const * const particle1, Particle const * const particle2) {
816 
817  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
818  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
819 //jcd to be removed
820 // return 0.;
821 //jcd
822  if (iso != 0) {
823  return NNToNNOmegaIso(ener, iso);
824  }
825  else {
826  return 0.5*(NNToNNOmegaIso(ener, 0)+NNToNNOmegaIso(ener, 2));
827  }
828  }
829 
830  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaExcluIso(const G4double ener, const G4int iso) {
831 
832  const G4double Ecm=0.001*ener;
833  G4double sNNOmega; // pp->pp+eta
834  G4double sNNOmega1; // np->np+eta
835  G4double sthroot=std::sqrt(7.06);
836 
837  if(Ecm>=3.0744) { // By hand (JCD)
838  sNNOmega = 330.*(Ecm-sthroot)/(1.05+std::pow((Ecm-sthroot),2));
839  }
840  else if(Ecm>=2.65854){
841  sNNOmega = -1208.09757*std::pow(Ecm,3) + 10773.3322*std::pow(Ecm,2) - 31661.0223*Ecm + 30728.7241 ;
842  }
843  else {
844  sNNOmega = 0. ;
845  }
846 
847  G4double Mn=ParticleTable::getRealMass(Neutron)/1000.;
848  G4double Mp=ParticleTable::getRealMass(Proton)/1000.;
849  G4double Momega=ParticleTable::getRealMass(Omega)/1000.;
850  G4double Thr0=0.;
851  if (iso > 0) {
852  Thr0=2.*Mp+Momega;
853  }
854  else if (iso < 0) {
855  Thr0=2.*Mn+Momega;
856  }
857  else {
858  Thr0=Mn+Mp+Momega;
859  }
860 
861  if (sNNOmega < 1.e-9 || Ecm < Thr0) sNNOmega = 0.; // Thr0: Ecm threshold
862 
863  if (iso != 0) {
864  return sNNOmega/1000.; // parameterization in microbarn (not millibarn)!
865  }
866 
867  sNNOmega1 = 3*sNNOmega; // 3.0: ratio pn/pp
868 
869  sNNOmega = 2*sNNOmega1-sNNOmega;
870  if (sNNOmega < 1.e-9 || Ecm < Thr0) sNNOmega = 0.;
871 
872  return sNNOmega/1000.; // parameterization in microbarn (not millibarn)!
873  }
874 
875  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaExclu(Particle const * const particle1, Particle const * const particle2) {
876 //jcd to be removed
877 // return 0.;
878 //jcd
879 
880  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2);
881  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
882 
883  if (iso != 0) {
884  return NNToNNOmegaExcluIso(ener, iso);
885  }
886  else {
887  return 0.5*(NNToNNOmegaExcluIso(ener, 0)+NNToNNOmegaExcluIso(ener, 2));
888  }
889  }
890 
891 
892  G4double CrossSectionsMultiPionsAndResonances::NNToxPiNN(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
893  //
894  // Nucleon-Nucleon producing xpi pions cross sections
895  //
896 // assert(xpi>0 && xpi<=nMaxPiNN);
897 // assert(particle1->isNucleon() && particle2->isNucleon());
898 
899  G4double oldXS1Pi=CrossSectionsMultiPions::NNToxPiNN(1,particle1, particle2);
900  G4double oldXS2Pi=CrossSectionsMultiPions::NNToxPiNN(2,particle1, particle2);
901  G4double oldXS3Pi=CrossSectionsMultiPions::NNToxPiNN(3,particle1, particle2);
902  G4double oldXS4Pi=CrossSectionsMultiPions::NNToxPiNN(4,particle1, particle2);
903  G4double xsEtaOmega=NNToNNEta(particle1, particle2)+NNToNNOmega(particle1, particle2);
904  G4double newXS1Pi=0.;
905  G4double newXS2Pi=0.;
906  G4double newXS3Pi=0.;
907  G4double newXS4Pi=0.;
908 
909  if (xpi == 1) {
910  if (oldXS4Pi != 0. || oldXS3Pi != 0.)
911  newXS1Pi=oldXS1Pi;
912  else if (oldXS2Pi != 0.) {
913  newXS2Pi=oldXS2Pi-xsEtaOmega;
914  if (newXS2Pi < 0.)
915  newXS1Pi=oldXS1Pi-(xsEtaOmega-oldXS2Pi);
916  else
917  newXS1Pi=oldXS1Pi;
918  }
919  else
920  newXS1Pi=oldXS1Pi-xsEtaOmega;
921  return newXS1Pi;
922  }
923  else if (xpi == 2) {
924  if (oldXS4Pi != 0.)
925  newXS2Pi=oldXS2Pi;
926  else if (oldXS3Pi != 0.) {
927  newXS3Pi=oldXS3Pi-xsEtaOmega;
928  if (newXS3Pi < 0.)
929  newXS2Pi=oldXS2Pi-(xsEtaOmega-oldXS3Pi);
930  else
931  newXS2Pi=oldXS2Pi;
932  }
933  else {
934  newXS2Pi=oldXS2Pi-xsEtaOmega;
935  if (newXS2Pi < 0.)
936  newXS2Pi=0.;
937  }
938  return newXS2Pi;
939  }
940  else if (xpi == 3) {
941  if (oldXS4Pi != 0.) {
942  newXS4Pi=oldXS4Pi-xsEtaOmega;
943  if (newXS4Pi < 0.)
944  newXS3Pi=oldXS3Pi-(xsEtaOmega-oldXS4Pi);
945  else
946  newXS3Pi=oldXS3Pi;
947  }
948  else {
949  newXS3Pi=oldXS3Pi-xsEtaOmega;
950  if (newXS3Pi < 0.)
951  newXS3Pi=0.;
952  }
953  return newXS3Pi;
954  }
955  else if (xpi == 4) {
956  newXS4Pi=oldXS4Pi-xsEtaOmega;
957  if (newXS4Pi < 0.)
958  newXS4Pi=0.;
959  return newXS4Pi;
960  }
961 
962  else // should never reach this point
963  return 0.;
964  }
965 
966 
967  G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaOnePi(Particle const * const particle1, Particle const * const particle2) {
968  // Cross section for nucleon-nucleon producing one eta and one pion
969 
970  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
971  if (iso!=0)
972  return 0.;
973 
974  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta (= 2705.55 - 2018.563; 4074595.287720512986=2018.563*2018.563)
975  if (ener < 2018.563) return 0.;
976 
977  const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
978  const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
979 
980  return 0.25*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
981  }
982 
983  G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaOnePiOrDelta(Particle const * const particle1, Particle const * const particle2) {
984  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
985  if (ener < 2018.563) return 0.;
986  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
987 
988  const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
989  if (iso != 0)
990  return CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
991  else {
992  const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
993  return 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
994  }
995  }
996 
997  G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaTwoPi(Particle const * const particle1, Particle const * const particle2) {
998  //
999  // Nucleon-Nucleon producing one eta and two pions
1000  //
1001  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
1002  if (ener < 2018.563) return 0.;
1003  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1004 
1005 
1006  const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1007  if (iso != 0) {
1008  return CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
1009  }
1010  else {
1011  const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1012  return 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
1013  }
1014  }
1015 
1016  G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaThreePi(Particle const * const particle1, Particle const * const particle2) {
1017  //
1018  // Nucleon-Nucleon producing one eta and three pions
1019  //
1020 
1021  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
1022  if (ener < 2018.563) return 0.;
1023  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1024 
1025 
1026  const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1027  const G4double xs1pi2=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2);
1028  const G4double xs2pi2=CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
1029  if (iso != 0)
1030  return CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2);
1031  else {
1032  const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1033  const G4double xs1pi0=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0);
1034  const G4double xs2pi0=CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0);
1035  return 0.5*(CrossSectionsMultiPions::NNThreePi(ener, 0, xsiso0, xs1pi0, xs2pi0)+ CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2));
1036  }
1037  }
1038 
1039  G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaFourPi(Particle const * const particle1, Particle const * const particle2) {
1040  //
1041  // Nucleon-Nucleon producing one eta and four pions
1042  //
1043 
1044  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
1045  if (ener < 2018.563) return 0.;
1046  const G4double s = ener*ener;
1047  const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1048  G4double xsinelas;
1049  if (i!=0)
1050  xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1051  else
1052  xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1053  if (xsinelas <= 1.e-9) return 0.;
1054  G4double ratio=(NNToNNEta(particle1, particle2)-NNToNNEtaExclu(particle1, particle2))/xsinelas;
1055  if(s<6.25E6)
1056  return 0.;
1057  const G4double sigma = NNToNNEta(particle1, particle2) - NNToNNEtaExclu(particle1, particle2) - ratio*(NNToNNEtaOnePiOrDelta(particle1, particle2) + NNToNNEtaTwoPi(particle1, particle2) + NNToNNEtaThreePi(particle1, particle2));
1058  return ((sigma>1.e-9) ? sigma : 0.);
1059  }
1060 
1061  G4double CrossSectionsMultiPionsAndResonances::NNToNNEtaxPi(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
1062  //
1063  // Nucleon-Nucleon producing one eta and xpi pions
1064  //
1065 // assert(xpi>0 && xpi<=nMaxPiNN);
1066 // assert(particle1->isNucleon() && particle2->isNucleon());
1067 
1068  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
1069  if (ener < 2018.563) return 0.;
1070  const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1071  G4double xsinelas;
1072  if (i!=0)
1073  xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1074  else
1075  xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1076  if (xsinelas <= 1.e-9) return 0.;
1077  G4double ratio=(NNToNNEta(particle1, particle2)-NNToNNEtaExclu(particle1, particle2))/xsinelas;
1078 
1079  if (xpi == 1)
1080  return NNToNNEtaOnePi(particle1, particle2)*ratio;
1081  else if (xpi == 2)
1082  return NNToNNEtaTwoPi(particle1, particle2)*ratio;
1083  else if (xpi == 3)
1084  return NNToNNEtaThreePi(particle1, particle2)*ratio;
1085  else if (xpi == 4)
1086  return NNToNNEtaFourPi(particle1, particle2);
1087  else // should never reach this point
1088  return 0.;
1089  }
1090 
1091 
1092  G4double CrossSectionsMultiPionsAndResonances::NNToNDeltaEta(Particle const * const p1, Particle const * const p2) {
1093 // assert(p1->isNucleon() && p2->isNucleon());
1094  const G4int i = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1095  const G4double ener=KinematicsUtils::totalEnergyInCM(p1, p2) - 581.437; // 581.437 MeV translation to open pion production in NNEta
1096  if (ener < 2018.563) return 0.;
1097  G4double xsinelas;
1098  if (i!=0)
1099  xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1100  else
1101  xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1102  if (xsinelas <= 1.e-9) return 0.;
1103  G4double ratio=(NNToNNEta(p1, p2)-NNToNNEtaExclu(p1, p2))/xsinelas;
1104  G4double sigma = NNToNNEtaOnePiOrDelta(p1, p2)*ratio;
1105  if(i==0)
1106  sigma *= 0.5;
1107  return sigma;
1108  }
1109 
1110 
1111  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaOnePi(Particle const * const particle1, Particle const * const particle2) {
1112  // Cross section for nucleon-nucleon producing one omega and one pion
1113 
1114  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1115  if (iso!=0)
1116  return 0.;
1117 
1118  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega (= 2802. - 2018.563; 4074595.287720512986=2018.563*2018.563)
1119  if (ener < 2018.563) return 0.;
1120 
1121  const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1122  const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1123 
1124  return 0.25*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
1125  }
1126 
1127  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaOnePiOrDelta(Particle const * const particle1, Particle const * const particle2) {
1128  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
1129  if (ener < 2018.563) return 0.;
1130  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1131 
1132  const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1133  if (iso != 0)
1134  return CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
1135  else {
1136  const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1137  return 0.5*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
1138  }
1139  }
1140 
1141  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaTwoPi(Particle const * const particle1, Particle const * const particle2) {
1142  //
1143  // Nucleon-Nucleon producing one omega and two pions
1144  //
1145  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
1146  if (ener < 2018.563) return 0.;
1147  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1148 
1149  const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1150  if (iso != 0) {
1151  return CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
1152  }
1153  else {
1154  const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1155  return 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
1156  }
1157  }
1158 
1159  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaThreePi(Particle const * const particle1, Particle const * const particle2) {
1160  //
1161  // Nucleon-Nucleon producing one omega and three pions
1162  //
1163 
1164  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
1165  if (ener < 2018.563) return 0.;
1166  const G4int iso=ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1167 
1168 
1169  const G4double xsiso2=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1170  const G4double xs1pi2=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2);
1171  const G4double xs2pi2=CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
1172  if (iso != 0)
1173  return CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2);
1174  else {
1175  const G4double xsiso0=CrossSectionsMultiPions::NNInelasticIso(ener, 0);
1176  const G4double xs1pi0=CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0);
1177  const G4double xs2pi0=CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0);
1178  return 0.5*(CrossSectionsMultiPions::NNThreePi(ener, 0, xsiso0, xs1pi0, xs2pi0)+ CrossSectionsMultiPions::NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2));
1179  }
1180  }
1181 
1182  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaFourPi(Particle const * const particle1, Particle const * const particle2) {
1183  //
1184  // Nucleon-Nucleon producing one omega and four pions
1185  //
1186 //jcd to be removed
1187 // return 0.;
1188 //jcd
1189 
1190  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
1191  if (ener < 2018.563) return 0.;
1192  const G4double s = ener*ener;
1193  const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1194  G4double xsinelas;
1195  if (i!=0)
1196  xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1197  else
1198  xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1199  if (xsinelas <= 1.e-9) return 0.;
1200  G4double ratio=(NNToNNOmega(particle1, particle2)-NNToNNOmegaExclu(particle1, particle2))/xsinelas;
1201  if(s<6.25E6)
1202  return 0.;
1203  const G4double sigma = NNToNNOmega(particle1, particle2) - NNToNNOmegaExclu(particle1, particle2) - ratio*(NNToNNOmegaOnePiOrDelta(particle1, particle2) + NNToNNOmegaTwoPi(particle1, particle2) + NNToNNOmegaThreePi(particle1, particle2));
1204  return ((sigma>1.e-9) ? sigma : 0.);
1205  }
1206 
1207  G4double CrossSectionsMultiPionsAndResonances::NNToNNOmegaxPi(const G4int xpi, Particle const * const particle1, Particle const * const particle2) {
1208  //
1209  // Nucleon-Nucleon producing one omega and xpi pions
1210  //
1211 // assert(xpi>0 && xpi<=nMaxPiNN);
1212 // assert(particle1->isNucleon() && particle2->isNucleon());
1213 //jcd to be removed
1214 // return 0.;
1215 //jcd
1216 
1217  const G4double ener=KinematicsUtils::totalEnergyInCM(particle1, particle2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
1218  if (ener < 2018.563) return 0.;
1219  const G4int i = ParticleTable::getIsospin(particle1->getType()) + ParticleTable::getIsospin(particle2->getType());
1220  G4double xsinelas;
1221  if (i!=0)
1222  xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1223  else
1224  xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1225  if (xsinelas <= 1.e-9) return 0.;
1226  G4double ratio=(NNToNNOmega(particle1, particle2)-NNToNNOmegaExclu(particle1, particle2))/xsinelas;
1227 
1228  if (xpi == 1)
1229  return NNToNNOmegaOnePi(particle1, particle2)*ratio;
1230  else if (xpi == 2)
1231  return NNToNNOmegaTwoPi(particle1, particle2)*ratio;
1232  else if (xpi == 3)
1233  return NNToNNOmegaThreePi(particle1, particle2)*ratio;
1234  else if (xpi == 4)
1235  return NNToNNOmegaFourPi(particle1, particle2);
1236  else // should never reach this point
1237  return 0.;
1238  }
1239 
1240 
1241  G4double CrossSectionsMultiPionsAndResonances::NNToNDeltaOmega(Particle const * const p1, Particle const * const p2) {
1242 // assert(p1->isNucleon() && p2->isNucleon());
1243 //jcd to be removed
1244 // return 0.;
1245 //jcd
1246  const G4int i = ParticleTable::getIsospin(p1->getType()) + ParticleTable::getIsospin(p2->getType());
1247  const G4double ener=KinematicsUtils::totalEnergyInCM(p1, p2) - 783.437; // 783.437 MeV translation to open pion production in NNOmega
1248  if (ener < 2018.563) return 0.;
1249  G4double xsinelas;
1250  if (i!=0)
1251  xsinelas=CrossSectionsMultiPions::NNInelasticIso(ener, 2);
1252  else
1253  xsinelas=0.5*(CrossSectionsMultiPions::NNInelasticIso(ener, 0)+CrossSectionsMultiPions::NNInelasticIso(ener, 2));
1254  if (xsinelas <= 1.e-9) return 0.;
1255  G4double ratio=(NNToNNOmega(p1, p2)-NNToNNOmegaExclu(p1, p2))/xsinelas;
1256  G4double sigma = NNToNNOmegaOnePiOrDelta(p1, p2)*ratio;
1257  if(i==0)
1258  sigma *= 0.5;
1259  return sigma;
1260  }
1261 
1262 
1263 } // namespace G4INCL
1264