| Geant4 Cross Reference |
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
38 /** \file G4INCLCrossSectionsMultiPionsAndResonances.cc
39 * \brief Multipion and mesonic Resonances cross sections
40 *
41 * \date 4th February 2014
42 * \author Jean-Christophe David
43 */
44
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
1265