Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4INCLCrossSectionsMultiPionsAndResonances.cc
Go to the documentation of this file.
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
48// #include <cassert>
49
50namespace 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;
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
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;
74 const G4double CrossSectionsMultiPionsAndResonances::s02pmOOT = 0.0016661112962345883443;
75 const G4double CrossSectionsMultiPionsAndResonances::s12mzOOT = 0.0017047391749062392793;
76
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
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())) {
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
119 if((p1->isNucleon()||p1->isDelta()) && (p2->isNucleon()||p2->isDelta())){
121 }
122 else if ((p1->isNucleon() && p2->isPion()) || (p2->isNucleon() && p1->isPion())){
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 {
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
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
551 //
552 // Pion-Nucleon producing Eta cross sections
553 //
554
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}
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
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
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
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
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
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
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
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
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
979
980 return 0.25*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
981 }
982
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
989 if (iso != 0)
990 return CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
991 else {
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
1007 if (iso != 0) {
1008 return CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
1009 }
1010 else {
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
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 {
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)
1051 else
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)
1074 else
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
1093// assert(p1->isNucleon() && p2->isNucleon());
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)
1100 else
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
1123
1124 return 0.25*(CrossSectionsMultiPions::NNOnePiOrDelta(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNOnePiOrDelta(ener, 2, xsiso2));
1125 }
1126
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
1133 if (iso != 0)
1134 return CrossSectionsMultiPions::NNOnePiOrDelta(ener, iso, xsiso2);
1135 else {
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
1150 if (iso != 0) {
1151 return CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2);
1152 }
1153 else {
1155 return 0.5*(CrossSectionsMultiPions::NNTwoPi(ener, 0, xsiso0)+ CrossSectionsMultiPions::NNTwoPi(ener, 2, xsiso2));
1156 }
1157 }
1158
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
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 {
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)
1197 else
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)
1223 else
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
1242// assert(p1->isNucleon() && p2->isNucleon());
1243//jcd to be removed
1244// return 0.;
1245//jcd
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)
1252 else
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
Multipion and mesonic Resonances cross sections.
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
virtual G4double NNToNNOmegaExcluIso(const G4double ener, const G4int iso)
Isotopic Cross section for Omega production (exclusive) - NN entrance channel.
virtual G4double piNToEtaN(Particle const *const p1, Particle const *const p2)
Cross sections for mesonic resonance production - piN Channel.
static const G4double s11pmOOT
One over threshold for s11pm.
virtual G4double etaNElastic(Particle const *const p1, Particle const *const p2)
Cross sections for mesonic resonance absorption on nucleon - elastic Channel.
virtual G4double omegaNToPiPiN(Particle const *const p1, Particle const *const p2)
Cross sections for omega-induced 2Pi emission on nucleon.
virtual G4double NNToNNEtaIso(const G4double ener, const G4int iso)
Cross section for One (more) pion production - piN entrance channel.
virtual G4double omegaNToPiN(Particle const *const p1, Particle const *const p2)
Cross section for OmegaN->PiN.
virtual G4double NNToNNOmegaExclu(Particle const *const particle1, Particle const *const particle2)
Cross section for Omega production (exclusive) - NN entrance channel.
virtual G4double NNToNDeltaEta(Particle const *const p1, Particle const *const p2)
Cross section for N-Delta-Eta production - NNEta Channel.
virtual G4double NNToNNOmegaTwoPi(Particle const *const part1, Particle const *const part2)
Cross section for direct 2-pion production - NNOmega channel.
virtual G4double NNToNNOmegaOnePiOrDelta(Particle const *const part1, Particle const *const part2)
Cross section for direct 1-pion production - NNOmega channel.
virtual G4double piNToxPiN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Cross section for X pion production - piN Channel (modified due to the mesonic resonances)
virtual G4double etaNToPiPiN(Particle const *const p1, Particle const *const p2)
Cross sections for mesonic resonance absorption on nucleon - pipiN Channel.
virtual G4double piNToOmegaN(Particle const *const p1, Particle const *const p2)
Cross section for PiN->OmegaN.
virtual G4double NNToNNOmegaThreePi(Particle const *const part1, Particle const *const part2)
Cross section for direct 3-pion production - NNOmega channel.
virtual G4double omegaNInelastic(Particle const *const p1, Particle const *const p2)
Cross sections for mesonic resonance absorption on nucleon - inelastic Channel.
static const G4int nMaxPiNN
Maximum number of outgoing pions in NN collisions.
virtual G4double NNToNNOmegaIso(const G4double ener, const G4int iso)
Isotopic Cross section for Omega production (inclusive) - NN entrance channel.
static const G4double s12ppOOT
One over threshold for s12pp.
virtual G4double NNToNNEta(Particle const *const particle1, Particle const *const particle2)
Cross section for Eta production (inclusive) - NN entrance channel.
static const G4double s11pzOOT
One over threshold for s11pz.
virtual G4double elastic(Particle const *const p1, Particle const *const p2)
new elastic particle-particle cross section
virtual G4double etaPrimeNToPiN(Particle const *const p1, Particle const *const p2)
Cross section for EtaPrimeN->PiN.
virtual G4double total(Particle const *const p1, Particle const *const p2)
new total particle-particle cross section
virtual G4double NNToNNEtaExcluIso(const G4double ener, const G4int iso)
Isotopic Cross section for Eta production (exclusive) - NN entrance channel.
virtual G4double NNToNNOmegaxPi(const G4int xpi, Particle const *const p1, Particle const *const p2)
Cross section for X pion production - NNOmega Channel.
virtual G4double NNToxPiNN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Cross section for X pion production - NN Channel.
virtual G4double NNToNNEtaExclu(Particle const *const particle1, Particle const *const particle2)
Cross section for Eta production (exclusive) - NN entrance channel.
virtual G4double NNToNNEtaOnePiOrDelta(Particle const *const part1, Particle const *const part2)
Cross section for direct 1-pion production - NNEta channel.
virtual G4double etaNToPiN(Particle const *const p1, Particle const *const p2)
Cross sections for mesonic resonance absorption on nucleon - piN Channel.
virtual G4double omegaNElastic(Particle const *const p1, Particle const *const p2)
static const G4double s01ppOOT
One over threshold for s01pp.
virtual G4double piNToEtaPrimeN(Particle const *const p1, Particle const *const p2)
Cross section for PiN->EtaPrimeN.
static const G4double s12zzOOT
One over threshold for s12zz.
static const G4double s02pzOOT
One over threshold for s02pz.
static const G4int nMaxPiPiN
Maximum number of outgoing pions in piN collisions.
static const G4double s12mzOOT
One over threshold for s12mz.
virtual G4double NNToNNOmegaFourPi(Particle const *const part1, Particle const *const part2)
Cross section for direct 4-pion production - NNOmega channel.
virtual G4double NNToNNEtaThreePi(Particle const *const part1, Particle const *const part2)
Cross section for direct 3-pion production - NNEta channel.
virtual G4double NNToNNEtaTwoPi(Particle const *const part1, Particle const *const part2)
Cross section for direct 2-pion production - NNEta channel.
static const G4double s02pmOOT
One over threshold for s02pm.
virtual G4double NNToNNOmega(Particle const *const particle1, Particle const *const particle2)
Cross section for Omega production (inclusive) - NN entrance channel.
virtual G4double NNToNNEtaFourPi(Particle const *const part1, Particle const *const part2)
Cross section for direct 4-pion production - NNEta channel.
virtual G4double NNToNNEtaOnePi(Particle const *const part1, Particle const *const part2)
Cross section for direct 1-pion production - NNEta channel.
static const G4double s12pmOOT
One over threshold for s12pm.
G4double piMinuspToEtaN(Particle const *const p1, Particle const *const p2)
Internal function for pion cross sections.
static const G4double s01pzOOT
One over threshold for s01pz.
virtual G4double NNToNNEtaxPi(const G4int xpi, Particle const *const p1, Particle const *const p2)
Cross section for X pion production - NNEta Channel.
virtual G4double NNToNDeltaOmega(Particle const *const p1, Particle const *const p2)
Cross section for N-Delta-Eta production - NNOmega Channel.
virtual G4double NNToNNOmegaOnePi(Particle const *const part1, Particle const *const part2)
Cross section for direct 1-pion production - NNOmega channel.
G4double piMinuspToOmegaN(Particle const *const p1, Particle const *const p2)
virtual G4double piNToxPiN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Cross section for X pion production - piN Channel.
G4double piNTot(Particle const *const p1, Particle const *const p2)
virtual G4double elastic(Particle const *const p1, Particle const *const p2)
Elastic particle-particle cross section.
G4double NNTot(Particle const *const part1, Particle const *const part2)
Internal implementation of the NN total cross section.
virtual G4double NDeltaToNN(Particle const *const p1, Particle const *const p2)
Cross section for NDelta->NN.
virtual G4double NNTwoPi(const G4double ener, const G4int iso, const G4double xsiso)
Cross section for direct 2-pion production - NN entrance channel.
virtual G4double NNOnePiOrDelta(const G4double ener, const G4int iso, const G4double xsiso)
Cross section for direct 1-pion production + delta production - NN entrance channel.
virtual G4double NNThreePi(const G4double ener, const G4int iso, const G4double xsiso, const G4double xs1pi, const G4double xs2pi)
Cross section for direct 3-pion production - NN entrance channel.
virtual G4double NNToxPiNN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Cross section for X pion production - NN Channel.
G4double NNInelasticIso(const G4double ener, const G4int iso)
Internal implementation of the isospin dependent NN reaction cross section.
G4bool isEtaPrime() const
Is this an etaprime?
G4bool isOmega() const
Is this an omega?
G4bool isEta() const
Is this an eta?
G4bool isPion() const
Is this a pion?
G4INCL::ParticleType getType() const
G4double getMass() const
Get the cached particle mass.
G4bool isDelta() const
Is it a Delta?
G4bool isNucleon() const
G4double totalEnergyInCM(Particle const *const p1, Particle const *const p2)
G4double momentumInLab(Particle const *const p1, Particle const *const p2)
gives the momentum in the lab frame of two particles.
G4double momentumInCM(Particle const *const p1, Particle const *const p2)
gives the momentum in the CM frame of two particles.
G4double getRealMass(const G4INCL::ParticleType t)
Get particle mass (in MeV/c^2)
G4double getINCLMass(const G4int A, const G4int Z, const G4int S)
Get INCL nuclear mass (in MeV/c^2)
const G4double effectiveNucleonMass2
G4int getIsospin(const ParticleType t)
Get the isospin of a particle.
const G4double effectiveNucleonMass
static G4double eval(const G4double pLab, const G4double oneOverThreshold, HornerCoefficients< N > const &coeffs)
static G4double eval(const G4double x, HornerCoefficients< N > const &coeffs)