34 #define INCLXX_IN_GEANT4_MODE 1
51 const G4double xrat=ekin*oneOverThreshold;
72 s11pzHC(-2.228000000000294018,8.7560000000005723725,-0.61000000000023239325,-5.4139999999999780324,3.3338333333333348023,-0.75835000000000022049,0.060623611111111114688),
73 s01ppHC(2.0570000000126518344,-6.029000000012135826,36.768500000002462784,-45.275666666666553533,25.112666666666611953,-7.2174166666666639187,1.0478875000000000275,-0.060804365079365080846),
74 s01pzHC(0.18030000000000441851,7.8700999999999953598,-4.0548999999999990425,0.555199999999999959),
75 s11pmHC(0.20590000000000031866,3.3450999999999993936,-1.4401999999999997825,0.17076666666666664973),
76 s12pmHC(-0.77235999999999901328,4.2626599999999991117,-1.9008899999999997323,0.30192266666666663379,-0.012270833333333331986),
77 s12ppHC(-0.75724999999999975664,2.0934399999999998565,-0.3803099999999999814),
78 s12zzHC(-0.89599999999996965072,7.882999999999978632,-7.1049999999999961928,1.884333333333333089),
79 s02pzHC(-1.0579999999999967036,11.113999999999994089,-8.5259999999999990196,2.0051666666666666525),
80 s02pmHC(2.4009000000012553286,-7.7680000000013376183,20.619000000000433505,-16.429666666666723928,5.2525708333333363472,-0.58969166666666670206),
81 s12mzHC(-0.21858699999999976269,1.9148999999999999722,-0.31727500000000001065,-0.027695000000000000486)
104 return 34.*std::pow(plab/0.4, (-2.104));
106 else if (plab < 0.800) {
107 return 23.5+1000.*std::pow(plab-0.7, 4);
109 else if (plab <= 2.0) {
110 return 1250./(50.+plab)-4.*std::pow(plab-1.3, 2);
113 return 77./(plab+1.5);
128 return 6.3555*std::exp(-3.2481*alp-0.377*alp*alp);
130 else if (plab < 0.851) {
131 return 33.+196.*std::pow(std::fabs(plab-0.95),2.5);
133 else if (plab <= 2.0) {
134 return 31./std::sqrt(plab);
137 return 77./(plab+1.5);
142 return 34.*std::pow(plab/0.4, (-2.104));
144 else if (plab < 0.8067) {
145 return 23.5+1000.*std::pow(plab-0.7, 4);
147 else if (plab <= 2.0) {
148 return 1250./(50.+plab)-4.*std::pow(plab-1.3, 2);
150 else if (plab <= 3.0956) {
151 return 77./(plab+1.5);
155 return 11.2+25.5*std::pow(plab, -1.12)+0.151*std::pow(alp, 2)-1.62*alp;
187 return 6.3555*std::exp(-3.2481*alp-0.377*std::pow(alp, 2));
189 else if (plab < 1.0) {
190 return 33.+196.*std::sqrt(std::pow(std::fabs(plab-0.95),5));
192 else if (plab < 1.924) {
193 return 24.2+8.9*plab;
197 return 48.9-33.7*std::pow(plab, -3.08)+0.619*std::pow(alp, 2)-5.12*alp;
202 return 34.*std::pow(plab/0.4, (-2.104));
204 else if (plab < 0.8734) {
205 return 23.5+1000.*std::pow(plab-0.7, 4);
207 else if (plab < 1.5) {
208 return 23.5+24.6/(1.+std::exp(-10.*(plab-1.2)));
210 else if (plab < 3.0044) {
211 return 41.+60.*(plab-0.9)*std::exp(-1.2*plab);
215 return 45.6+219.*std::pow(plab, -4.23)+0.41*std::pow(alp, 2)-3.41*alp;
226 if(s>=4074595.287720512986) {
233 if(s>=4074595.287720512986) {
240 if (sincl < 0.) sincl = 0.;
267 if ((iso != 0) && (plab < 2.1989)) {
268 snnpit = xsiso -
NNTwoPi(ener, iso, xsiso);
269 if (snnpit < 1.e-8) snnpit=0.;
272 else if ((iso == 0) && (plab < 1.7369)) {
274 if (snnpit < 1.e-8) snnpit=0.;
280 s11pz=55.185/std::pow((0.1412*plab+5),2);
282 else if (plab > 13.9) {
284 s11pz=6.67-13.3*std::pow(plab, -6.18)+0.456*alp*alp-3.29*alp;
286 else if (plab >= 0.7765) {
291 if (plab >= 0.79624) {
298 if (snnpit1 < 1.e-8) snnpit1=0.;
305 s01pz=15289.4/std::pow((11.573*plab+5),2);
307 else if (plab >= 0.777) {
313 s11pm=46.68/std::pow((0.2231*plab+5),2);
315 else if (plab >= 0.788) {
321 snnpit=2*(s01pz+2*s11pm)-snnpit1;
322 if (snnpit < 1.e-8) snnpit=0.;
349 if (iso==0 && plab<3.3) {
351 if (snn2pit < 1.e-8) snn2pit=0.;
360 else if (plab >= 1.3817) {
366 s12pp=141.505/std::pow((-0.1016*plab-7),2);
368 else if (plab >= 1.5739) {
375 s12zz=97.355/std::pow((1.1579*plab+5),2);
377 else if (plab >= 1.72207) {
383 s02pz=178.082/std::pow((0.2014*plab+5),2);
385 else if (plab >= 1.5656) {
392 snn2pit=s12pm+s12pp+s12zz+s02pz;
393 if (snn2pit < 1.e-8) snn2pit=0.;
399 s02pm=135.826/std::pow(plab,2);
401 else if (plab >= 1.21925) {
406 if (plab >= 1.29269) {
412 snn2pit=3*(s02pm+0.5*s12mz-0.5*s02pz-s12zz);
413 if (snn2pit < 1.e-8) snn2pit=0.;
427 return 46.72/std::pow((plab - 5.8821),2);
430 snn3pit=xsiso-xs1pi-xs2pi;
431 if (snn3pit < 1.e-8) snn3pit=0.;
438 return 5592.92/std::pow((plab+14.9764),2);
440 else if (plab > 2.1989){
441 snn3pit=xsiso-xs1pi-xs2pi;
442 if (snn3pit < 1.e-8) snn3pit=0.;
487 return NNTwoPi(ener, 2, xsiso2);
509 return NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2);
514 return 0.5*(
NNThreePi(ener, 0, xsiso0, xs1pi0, xs2pi0)+
NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2));
523 return ((sigma>1.e-9) ? sigma : 0.);
534 return NNOnePi(particle1, particle2);
536 return NNTwoPi(particle1, particle2);
540 return NNFourPi(particle1, particle2);
553 q2=(y-std::pow(1076.0, 2))*(y-std::pow(800.0, 2))/(4.0*y);
556 G4double f3=q3/(q3+std::pow(180.0, 3));
558 sdel=326.5/(std::pow((x-1215.0-ramass)*2.0/110.0,2)+1.0);
559 return sdel*f3*(1.0-5.0*ramass/1215.0);
566 return -2.33730e-06*std::pow(x, 3)+1.13819e-02*std::pow(x,2)
567 -1.83993e+01*x+9893.4;
568 }
else if (x <= 2150.0) {
569 return 1.13531e-06*std::pow(x, 3)-6.91694e-03*std::pow(x, 2)
570 +1.39907e+01*x-9360.76;
572 return -3.18087*std::log(x)+52.9784;
583 q2=(y-std::pow(1076.0, 2))*(y-std::pow(800.0, 2))/(4.0*y);
586 G4double f3=q3/(q3+std::pow(180.0, 3));
588 sdel=326.5/(std::pow((x-1215.0-ramass)*2.0/110.0,2)+1.0);
589 return sdel*f3*(1.0-5.0*ramass/1215.0)/3.;
596 return 0.00120683*(x-1372.52)*(x-1372.52)+26.2058;
597 }
else if(x <= 1578.0) {
598 return 1.15873e-05*x*x+49965.6/((x-1519.59)*(x-1519.59)+2372.55);
599 }
else if(x <= 2028.4) {
600 return 34.0248+43262.2/((x-1681.65)*(x-1681.65)+1689.35);
601 }
else if(x <= 7500.0) {
602 return 3.3e-7*(x-7500.0)*(x-7500.0)+24.5;
611 return NNTot(p1, p2);
622 return inelastic +
elastic(p1, p2);
643 if(pLab>212677. || pLab<296.367)
648 const G4int cg = 4 + ind2t3*ipit3;
670 return 0.5*(xpipp+xpimp);
677 if(x>20000.)
return 0.0;
686 }
else if(particle2->
isPion()) {
692 const G4double q2=(y-1076.0*1076.0)*(y-800.0*800.0)/y/4.0;
696 const G4double q3 = std::pow(std::sqrt(q2),3);
697 const G4double f3 = q3/(q3 + 5832000.);
698 G4double sdelResult = 326.5/(std::pow((x-1215.0-ramass)*2.0/(110.0-ramass), 2)+1.0);
699 sdelResult = sdelResult*(1.0-5.0*ramass/1215.0);
700 const G4int cg = 4 + ind2t3*ipit3;
701 sdelResult = sdelResult*f3*cg/6.0;
723 }
else if(particle2->
isPion()) {
731 if((ind2t3 == 1 && ipit3 == 2) || (ind2t3 == -1 && ipit3 == -2))
733 else if((ind2t3 == 1 && ipit3 == -2) || (ind2t3 == -1 && ipit3 == 2))
745 if(isospin==4 || isospin==-4)
return 0.0;
759 if(Ecm <= 938.3 + deltaMass) {
763 if(Ecm < 938.3 + deltaMass + 2.0) {
764 Ecm = 938.3 + deltaMass + 2.0;
787 result *= 3.*(32.0 + isospin * isospin * (deltaIsospin * deltaIsospin - 5))/64.0;
788 result /= 1.0 + 0.25 * (isospin * isospin);
825 return 5.5e-6 * x/(7.7 + x);
827 return (5.34 + 0.67*(x - 2.0)) * 1.0e-6;
831 G4double b = (7.16 - 1.63*x) * 1.0e-6;
832 return b/(1.0 + std::exp(-(x - 0.45)/0.05));
833 }
else if(pl < 1100.0) {
834 return (9.87 - 4.88 * x) * 1.0e-6;
836 return (3.68 + 0.76*x) * 1.0e-6;
851 if (OnePi < 1.e-09) OnePi = 0.;
856 if (TwoPi < 1.e-09) TwoPi = 0.;
861 if (piNThreePi < 1.e-09) piNThreePi = 0.;
887 const G4int cg = 4 + ind2*ipi;
903 if (tamp2 < 0.0) tamp2=0;
910 if (s1pin > inelastic)
941 const G4int cg = 4 + ind2*ipi;
961 const G4double s2pin=0.5*(tamp6+tamp2);
983 if(pLab>212677. || pLab<296.367)
997 xpipp=17.965*std::pow(p1, 5.4606);
999 xpipp=24.3-12.3*std::pow(p1, -1.91)+0.324*p2*p2-2.44*p2;
1011 nucleon = particle1;
1015 nucleon = particle2;
1022 if(pLab>212677. || pLab<296.367)
1038 xpimp=26.6-7.18*std::pow(p1, -1.86)+0.327*p2*p2-2.81*p2;
1051 nucleon = particle1;
1055 nucleon = particle2;
1077 tamp6=0.204+18.2*std::pow(p1, -1.72)+6.33*std::pow(p1, -1.13);
1088 nucleon = particle1;
1092 nucleon = particle2;
1114 tamp2=9.04*std::pow(p1, -1.17)+18.*std::pow(p1, -1.21);
1115 if (tamp2 < 0.0) tamp2=0;
1130 nucleon = particle1;
1134 nucleon = particle2;
1156 tamp6=1.59+25.5*std::pow(p1, -1.04);
1172 nucleon = particle1;
1176 nucleon = particle2;
1198 tamp2=2.457794117647+18.066176470588*std::pow(p1, -0.92);
G4double G4ParticleHPJENDLHEData::G4double result
virtual G4double NNToxPiNN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Cross section for X pion production - NN Channel.
static const G4double s01pzOOT
One over threshold for s01pz.
virtual G4double NNToNDelta(Particle const *const p1, Particle const *const p2)
Cross section for Delta production - NN 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 piNToEtaN(Particle const *const p1, Particle const *const p2)
Cross sections for mesonic resonance production - piN Channel.
G4double getMass() const
Get the cached particle mass.
G4double NNInelasticIso(const G4double ener, const G4int iso)
Internal implementation of the isospin dependent NN reaction cross section.
virtual G4double NNToNNOmega(Particle const *const particle1, Particle const *const particle2)
Cross section for Eta production - NN entrance channel.
static const G4double s12mzOOT
One over threshold for s12mz.
virtual G4double etaNToPiN(Particle const *const p1, Particle const *const p2)
Cross sections for mesonic resonance absorption on nucleon - piN Channel.
G4double piNIne(Particle const *const p1, Particle const *const p2)
virtual G4double NNTwoPi(const G4double ener, const G4int iso, const G4double xsiso)
Cross section for direct 2-pion production - NN entrance channel.
static const G4double s12ppOOT
One over threshold for s12pp.
G4double squareTotalEnergyInCM(Particle const *const p1, Particle const *const p2)
virtual G4double piNToDelta(Particle const *const p1, Particle const *const p2)
Cross section for Delta production - piN Channel.
G4double NNElasticFixed(const G4double s, const G4int i)
Internal implementation of the NN elastic cross section with fixed isospin.
const HornerC4 s12zzHC
Horner coefficients for s12zz.
virtual G4double etaPrimeNToPiN(Particle const *const p1, Particle const *const p2)
Cross section for EtaPrimeN->PiN.
G4bool isDelta() const
Is it a Delta?
std::string print() const
virtual G4double NNToNDeltaOmega(Particle const *const p1, Particle const *const p2)
Cross section for N-Delta-Eta production - NNEta Channel.
virtual G4double piNToxPiN(const G4int xpi, Particle const *const p1, Particle const *const p2)
Cross section for X pion production - piN Channel.
G4double spnPiPlusPHE(const G4double x)
Internal function for pion cross sections.
const HornerC4 s02pzHC
Horner coefficients for s02pz.
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.
G4double NNElastic(Particle const *const part1, Particle const *const part2)
Internal implementation of the NN elastic cross section.
const HornerC8 s01ppHC
Horner coefficients for s01pp.
virtual G4double NNOnePi(Particle const *const part1, Particle const *const part2)
Cross section for direct 1-pion production - NN entrance channel.
static const G4double s01ppOOT
One over threshold for s01pp.
G4double NNTot(Particle const *const part1, Particle const *const part2)
Internal implementation of the NN total cross section.
static const G4int nMaxPiNN
Maximum number of outgoing pions in NN collisions.
virtual G4double piNTwoPi(Particle const *const p1, Particle const *const p2)
Cross section for Two (more) pion production - piN entrance channel.
G4double piPluspOnePi(Particle const *const p1, Particle const *const p2)
virtual G4double piNToEtaPrimeN(Particle const *const p1, Particle const *const p2)
Cross section for PiN->EtaPrimeN.
virtual G4double omegaNToPiN(Particle const *const p1, Particle const *const p2)
Cross section for OmegaN->PiN.
const HornerC4 s01pzHC
Horner coefficients for s01pz.
virtual G4double NNToNNOmegaExclu(Particle const *const particle1, Particle const *const particle2)
Cross section for Eta production (exclusive) - NN entrance channel.
CrossSectionsMultiPions()
G4bool nucleon(G4int ityp)
virtual G4double NNToNNOmegaxPi(const G4int xpi, Particle const *const p1, Particle const *const p2)
Cross section for X pion production - NNEta Channel.
G4double piMinuspTwoPi(Particle const *const p1, Particle const *const p2)
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.
const HornerC3 s12ppHC
Horner coefficients for s12pp.
static const G4double s11pmOOT
One over threshold for s11pm.
virtual G4double NNFourPi(Particle const *const part1, Particle const *const part2)
Cross section for direct 4-pion production - NN entrance channel.
virtual G4double NDeltaToNN(Particle const *const p1, Particle const *const p2)
Cross section for NDelta->NN.
virtual G4double NNToNNEtaExclu(Particle const *const particle1, Particle const *const particle2)
Cross section for Eta production (exclusive) - NN entrance channel.
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.
G4double piMinuspOnePi(Particle const *const p1, Particle const *const p2)
const HornerC7 s11pzHC
Horner coefficients for s11pz.
G4int getIsospin(const ParticleType t)
Get the isospin of a particle.
const G4double effectiveNucleonMass2
static const G4double s02pzOOT
One over threshold for s02pz.
G4INCL::ParticleType getType() const
static G4double eval(const G4double pLab, const G4double oneOverThreshold, HornerCoefficients< N > const &coeffs)
static const G4double s12pmOOT
One over threshold for s12pm.
static const G4double s12zzOOT
One over threshold for s12zz.
G4double NNTotFixed(const G4double s, const G4int i)
Internal implementation of the NN total cross section with fixed isospin.
G4double piPluspTwoPi(Particle const *const p1, Particle const *const p2)
static const G4double s11pzOOT
One over threshold for s11pz.
const HornerC4 s12mzHC
Horner coefficients for s12mz.
static const G4double s02pmOOT
One over threshold for s02pm.
virtual G4double NNToNNEta(Particle const *const particle1, Particle const *const particle2)
Cross section for Eta production - NN entrance channel.
virtual G4double NNToNNEtaxPi(const G4int xpi, Particle const *const p1, Particle const *const p2)
Cross section for X pion production - NNEta Channel.
G4double piMinuspIne(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.
virtual G4double calculateNNAngularSlope(G4double energyCM, G4int iso)
Calculate the slope of the NN DDXS.
G4double totalEnergyInCM(Particle const *const p1, Particle const *const p2)
virtual G4double omegaNToPiPiN(Particle const *const p1, Particle const *const p2)
Cross section for OmegaN->PiPiN.
const HornerC4 s11pmHC
Horner coefficients for s11pm.
G4bool isPion() const
Is this a pion?
G4double spnPiMinusPHE(const G4double x)
Internal function for pion cross sections.
static G4double eval(const G4double x, HornerCoefficients< N > const &coeffs)
virtual G4double total(Particle const *const p1, Particle const *const p2)
Total (elastic+inelastic) particle-particle cross section.
const G4double effectiveNucleonMass
static const G4int nMaxPiPiN
Maximum number of outgoing pions in piN collisions.
virtual G4double piNOnePi(Particle const *const p1, Particle const *const p2)
Cross section for One (more) pion production - piN entrance channel.
virtual G4double NNToNDeltaEta(Particle const *const p1, Particle const *const p2)
Cross section for N-Delta-Eta production - NNEta Channel.
Cross sections used in INCL Multipions.
const HornerC6 s02pmHC
Horner coefficients for s02pm.
G4double piPluspIne(Particle const *const p1, Particle const *const p2)
const HornerC5 s12pmHC
Horner coefficients for s12pm.