34#define INCLXX_IN_GEANT4_MODE 1
53 struct BystrickyEvaluator {
57 const G4double xrat=ekin*oneOverThreshold;
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)
114 return inelastic +
elastic(p1, p2);
156 else if (oldXS3Pi != 0.) {
157 newXS3Pi=oldXS3Pi-xsEta-xsOmega;
158 if (newXS3Pi < 1.e-09)
159 newXS2Pi=oldXS2Pi-(xsEta+xsOmega-oldXS3Pi);
164 newXS2Pi=oldXS2Pi-xsEta-xsOmega;
165 if (newXS2Pi < 1.e-09)
171 if (oldXS4Pi != 0.) {
172 newXS4Pi=oldXS4Pi-xsEta-xsOmega;
173 if (newXS4Pi < 1.e-09)
174 newXS3Pi=oldXS3Pi-(xsEta+xsOmega-oldXS4Pi);
179 newXS3Pi=oldXS3Pi-xsEta-xsOmega;
180 if (newXS3Pi < 1.e-09)
186 newXS4Pi=oldXS4Pi-xsEta-xsOmega;
187 if (newXS4Pi < 1.e-09)
208 else return 0.5 * sigma;
210 else if (isoin == 1) {
212 else return 0.5 * sigma;
232 else return 0.5 * sigma;
234 else if (isoin == 1) {
236 else return 0.5 * sigma;
243#if defined(NDEBUG) || defined(INCLXX_IN_GEANT4_MODE)
265 if (particle1->
isEta()) {
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;
292 massnucleon=nucleon->getMass();
298 if (sigma < 0.) sigma=0.;
314 if (particle1->
isEta()) {
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;
335 if (sigma < 0.) sigma = 0.;
351 if (particle1->
isEta()) {
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;
373 if (sigma < 0.) sigma = 0.;
399 sigma = 20. + 4.0/pLab;
427 sigma = 5.4 + 10.*std::exp(-0.6*pLab);
453 massomega=particle1->
getMass();
454 massnucleon=particle2->
getMass();
457 massomega=particle2->
getMass();
458 massnucleon=particle1->
getMass();
469 if (sigma >
omegaNInelastic(particle1, particle2) || (pLab_omega < 200.)) {
492#if defined(NDEBUG) || defined(INCLXX_IN_GEANT4_MODE)
513 if (particle1->
isPion()) {
515 massnucleon=particle2->
getMass();
519 massnucleon=particle1->
getMass();
528 if (ECM < 1486.5) sigma=0.;
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;
535 else if (ECM < 1670.)
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;
539 else if (ECM < 1714.)
541 sigma = 0.000003737765*std::pow(ECM,2) - 0.005664062*ECM;
543 else sigma=1.47*std::pow(plab, -1.68);
563 if (ECM < 1486.5) sigma=0.;
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;
570 else if (ECM < 1670.)
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;
574 else if (ECM < 1714.)
576 sigma = 0.000003737765*std::pow(ECM,2) - 0.005664062*ECM;
578 else sigma=1.47*std::pow(plab, -1.68);
598 if (particle1->
isPion()) {
600 massnucleon=particle2->
getMass();
604 massnucleon=particle1->
getMass();
610 if (plab < param) sigma=0.;
611 else sigma=13.76*(plab-param)/(std::pow(plab, 3.33) - 1.07);
632 if (plab < param) sigma=0.;
633 else sigma=13.76*(plab-param)/(std::pow(plab, 3.33)-1.07);
648 sNNEta = 2.5*std::pow((x-1.),1.47)*std::pow(x,-1.25)*1000.;
650 else if(Ecm >= 2.6) {
651 sNNEta = -327.29*Ecm*Ecm*Ecm + 2870.*Ecm*Ecm - 7229.3*Ecm + 5273.3;
658 if (sNNEta < 1.e-9) sNNEta = 0.;
667 else if (Ecm >= 2.6) {
668 sNNEta1 = sNNEta*std::exp(-(-5.53151576/Ecm+0.8850425));
670 else if (Ecm >= 2.525) {
671 sNNEta1 = -4433.586*Ecm*Ecm*Ecm*Ecm + 56581.54*Ecm*Ecm*Ecm - 270212.6*Ecm*Ecm + 571650.6*Ecm - 451091.6;
674 sNNEta1 = 17570.217219*Ecm*Ecm - 84910.985402*Ecm + 102585.55847;
677 sNNEta2 = -10220.89518466*Ecm*Ecm+51227.30841724*Ecm-64097.96025731;
678 if (sNNEta2 < 0.) sNNEta2=0.;
680 sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;
685 if (sNNEta < 1.e-9 || Ecm < Mn+Mp+Meta) sNNEta = 0.;
712 sNNEta = -13.008*Ecm*Ecm + 84.531*Ecm + 36.234;
714 else if(Ecm>=2.725) {
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;
717 else if(Ecm>=2.575) {
718 sNNEta = -2640.3*Ecm*Ecm + 14692*Ecm - 20225;
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;
738 if (sNNEta < 1.e-9 || Ecm < Thr0) sNNEta = 0.;
747 else if (Ecm >= 3.5) {
748 sNNEta1 = -1916.2*Ecm*Ecm*Ecm + 21556.0*Ecm*Ecm - 80828.0*Ecm + 101200.0;
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;
754 sNNEta1 = 17570.217219*Ecm*Ecm - 84910.985402*Ecm + 102585.55847;
757 sNNEta2 = -10220.89518466*Ecm*Ecm+51227.30841724*Ecm-64097.96025731;
758 if (sNNEta2 < 0.) sNNEta2=0.;
760 sNNEta = 2*(sNNEta1+sNNEta2)-sNNEta;
761 if (sNNEta < 1.e-9 || Ecm < Thr0) sNNEta = 0.;
789 sNNOmega = 2.5*std::pow(x-1, 1.47)*std::pow(x, -1.11) ;
792 sNNOmega = (568.5254*Ecm*Ecm - 2694.045*Ecm + 3106.247)/1000.;
799 if (sNNOmega < 1.e-9) sNNOmega = 0.;
805 sNNOmega1 = 3.*sNNOmega;
807 sNNOmega = 2.*sNNOmega1-sNNOmega;
809 if (sNNOmega < 1.e-9) sNNOmega = 0.;
838 sNNOmega = 330.*(Ecm-sthroot)/(1.05+std::pow((Ecm-sthroot),2));
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 ;
861 if (sNNOmega < 1.e-9 || Ecm < Thr0) sNNOmega = 0.;
864 return sNNOmega/1000.;
867 sNNOmega1 = 3*sNNOmega;
869 sNNOmega = 2*sNNOmega1-sNNOmega;
870 if (sNNOmega < 1.e-9 || Ecm < Thr0) sNNOmega = 0.;
872 return sNNOmega/1000.;
910 if (oldXS4Pi != 0. || oldXS3Pi != 0.)
912 else if (oldXS2Pi != 0.) {
913 newXS2Pi=oldXS2Pi-xsEtaOmega;
915 newXS1Pi=oldXS1Pi-(xsEtaOmega-oldXS2Pi);
920 newXS1Pi=oldXS1Pi-xsEtaOmega;
926 else if (oldXS3Pi != 0.) {
927 newXS3Pi=oldXS3Pi-xsEtaOmega;
929 newXS2Pi=oldXS2Pi-(xsEtaOmega-oldXS3Pi);
934 newXS2Pi=oldXS2Pi-xsEtaOmega;
941 if (oldXS4Pi != 0.) {
942 newXS4Pi=oldXS4Pi-xsEtaOmega;
944 newXS3Pi=oldXS3Pi-(xsEtaOmega-oldXS4Pi);
949 newXS3Pi=oldXS3Pi-xsEtaOmega;
956 newXS4Pi=oldXS4Pi-xsEtaOmega;
975 if (ener < 2018.563)
return 0.;
985 if (ener < 2018.563)
return 0.;
1002 if (ener < 2018.563)
return 0.;
1022 if (ener < 2018.563)
return 0.;
1045 if (ener < 2018.563)
return 0.;
1053 if (xsinelas <= 1.e-9)
return 0.;
1058 return ((sigma>1.e-9) ? sigma : 0.);
1069 if (ener < 2018.563)
return 0.;
1076 if (xsinelas <= 1.e-9)
return 0.;
1096 if (ener < 2018.563)
return 0.;
1102 if (xsinelas <= 1.e-9)
return 0.;
1119 if (ener < 2018.563)
return 0.;
1129 if (ener < 2018.563)
return 0.;
1146 if (ener < 2018.563)
return 0.;
1165 if (ener < 2018.563)
return 0.;
1191 if (ener < 2018.563)
return 0.;
1199 if (xsinelas <= 1.e-9)
return 0.;
1204 return ((sigma>1.e-9) ? sigma : 0.);
1218 if (ener < 2018.563)
return 0.;
1225 if (xsinelas <= 1.e-9)
return 0.;
1248 if (ener < 2018.563)
return 0.;
1254 if (xsinelas <= 1.e-9)
return 0.;
Multipion and mesonic Resonances cross sections.
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.
CrossSectionsMultiPionsAndResonances()
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?
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)