34 #define INCLXX_IN_GEANT4_MODE 1
51 const G4double xrat=ekin*oneOverThreshold;
69 s11pzHC(-2.228000000000294018,8.7560000000005723725,-0.61000000000023239325,-5.4139999999999780324,3.3338333333333348023,-0.75835000000000022049,0.060623611111111114688),
70 s01ppHC(2.0570000000126518344,-6.029000000012135826,36.768500000002462784,-45.275666666666553533,25.112666666666611953,-7.2174166666666639187,1.0478875000000000275,-0.060804365079365080846),
71 s01pzHC(0.18030000000000441851,7.8700999999999953598,-4.0548999999999990425,0.555199999999999959),
72 s11pmHC(0.20590000000000031866,3.3450999999999993936,-1.4401999999999997825,0.17076666666666664973),
73 s12pmHC(-0.77235999999999901328,4.2626599999999991117,-1.9008899999999997323,0.30192266666666663379,-0.012270833333333331986),
74 s12ppHC(-0.75724999999999975664,2.0934399999999998565,-0.3803099999999999814),
75 s12zzHC(-0.89599999999996965072,7.882999999999978632,-7.1049999999999961928,1.884333333333333089),
76 s02pzHC(-1.0579999999999967036,11.113999999999994089,-8.5259999999999990196,2.0051666666666666525),
77 s02pmHC(2.4009000000012553286,-7.7680000000013376183,20.619000000000433505,-16.429666666666723928,5.2525708333333363472,-0.58969166666666670206),
78 s12mzHC(-0.21858699999999976269,1.9148999999999999722,-0.31727500000000001065,-0.027695000000000000486)
101 return 34.*std::pow(plab/0.4, (-2.104));
103 else if (plab < 0.800) {
104 return 23.5+1000.*std::pow(plab-0.7, 4);
106 else if (plab <= 2.0) {
107 return 1250./(50.+plab)-4.*std::pow(plab-1.3, 2);
110 return 77./(plab+1.5);
125 return 6.3555*std::exp(-3.2481*alp-0.377*alp*alp);
127 else if (plab < 0.851) {
128 return 33.+196.*std::pow(std::fabs(plab-0.95),2.5);
130 else if (plab <= 2.0) {
131 return 31./std::sqrt(plab);
134 return 77./(plab+1.5);
139 return 34.*std::pow(plab/0.4, (-2.104));
141 else if (plab < 0.8067) {
142 return 23.5+1000.*std::pow(plab-0.7, 4);
144 else if (plab <= 2.0) {
145 return 1250./(50.+plab)-4.*std::pow(plab-1.3, 2);
147 else if (plab <= 3.0956) {
148 return 77./(plab+1.5);
152 return 11.2+25.5*std::pow(plab, -1.12)+0.151*std::pow(alp, 2)-1.62*alp;
184 return 6.3555*std::exp(-3.2481*alp-0.377*std::pow(alp, 2));
186 else if (plab < 1.0) {
187 return 33.+196.*std::sqrt(std::pow(std::fabs(plab-0.95),5));
189 else if (plab < 1.924) {
190 return 24.2+8.9*plab;
194 return 48.9-33.7*std::pow(plab, -3.08)+0.619*std::pow(alp, 2)-5.12*alp;
199 return 34.*std::pow(plab/0.4, (-2.104));
201 else if (plab < 0.8734) {
202 return 23.5+1000.*std::pow(plab-0.7, 4);
204 else if (plab < 1.5) {
205 return 23.5+24.6/(1.+std::exp(-10.*(plab-1.2)));
207 else if (plab < 3.0044) {
208 return 41.+60.*(plab-0.9)*std::exp(-1.2*plab);
212 return 45.6+219.*std::pow(plab, -4.23)+0.41*std::pow(alp, 2)-3.41*alp;
223 if(s>=4074595.287720512986) {
230 if(s>=4074595.287720512986) {
237 if (sincl < 0.) sincl = 0.;
264 if ((iso != 0) && (plab < 2.1989)) {
265 snnpit = xsiso -
NNTwoPi(ener, iso, xsiso);
266 if (snnpit < 1.e-8) snnpit=0.;
269 else if ((iso == 0) && (plab < 1.7369)) {
271 if (snnpit < 1.e-8) snnpit=0.;
277 s11pz=55.185/std::pow((0.1412*plab+5),2);
279 else if (plab > 13.9) {
281 s11pz=6.67-13.3*std::pow(plab, -6.18)+0.456*alp*alp-3.29*alp;
283 else if (plab >= 0.7765) {
288 if (plab >= 0.79624) {
295 if (snnpit1 < 1.e-8) snnpit1=0.;
302 s01pz=15289.4/std::pow((11.573*plab+5),2);
304 else if (plab >= 0.777) {
310 s11pm=46.68/std::pow((0.2231*plab+5),2);
312 else if (plab >= 0.788) {
318 snnpit=2*(s01pz+2*s11pm)-snnpit1;
319 if (snnpit < 1.e-8) snnpit=0.;
346 if (iso==0 && plab<3.3) {
348 if (snn2pit < 1.e-8) snn2pit=0.;
357 else if (plab >= 1.3817) {
363 s12pp=141.505/std::pow((-0.1016*plab-7),2);
365 else if (plab >= 1.5739) {
372 s12zz=97.355/std::pow((1.1579*plab+5),2);
374 else if (plab >= 1.72207) {
380 s02pz=178.082/std::pow((0.2014*plab+5),2);
382 else if (plab >= 1.5656) {
389 snn2pit=s12pm+s12pp+s12zz+s02pz;
390 if (snn2pit < 1.e-8) snn2pit=0.;
396 s02pm=135.826/std::pow(plab,2);
398 else if (plab >= 1.21925) {
403 if (plab >= 1.29269) {
409 snn2pit=3*(s02pm+0.5*s12mz-0.5*s02pz-s12zz);
410 if (snn2pit < 1.e-8) snn2pit=0.;
424 return 46.72/std::pow((plab - 5.8821),2);
427 snn3pit=xsiso-xs1pi-xs2pi;
428 if (snn3pit < 1.e-8) snn3pit=0.;
435 return 5592.92/std::pow((plab+14.9764),2);
437 else if (plab > 2.1989){
438 snn3pit=xsiso-xs1pi-xs2pi;
439 if (snn3pit < 1.e-8) snn3pit=0.;
484 return NNTwoPi(ener, 2, xsiso2);
506 return NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2);
511 return 0.5*(
NNThreePi(ener, 0, xsiso0, xs1pi0, xs2pi0)+
NNThreePi(ener, 2, xsiso2, xs1pi2, xs2pi2));
520 return ((sigma>0.) ? sigma : 0.);
531 return NNOnePi(particle1, particle2);
533 return NNTwoPi(particle1, particle2);
537 return NNFourPi(particle1, particle2);
548 q2=(y-std::pow(1076.0, 2))*(y-std::pow(800.0, 2))/(4.0*y);
553 sdel=326.5/(std::pow((x-1215.0-ramass)*2.0/110.0,2)+1.0);
554 return sdel*f3*(1.0-5.0*ramass/1215.0);
561 return -2.33730e-06*std::pow(x, 3)+1.13819e-02*std::pow(x,2)
562 -1.83993e+01*x+9893.4;
563 }
else if (x <= 2150.0) {
564 return 1.13531e-06*std::pow(x, 3)-6.91694e-03*std::pow(x, 2)
565 +1.39907e+01*x-9360.76;
567 return -3.18087*std::log(x)+52.9784;
578 q2=(y-std::pow(1076.0, 2))*(y-std::pow(800.0, 2))/(4.0*y);
583 sdel=326.5/(std::pow((x-1215.0-ramass)*2.0/110.0,2)+1.0);
584 return sdel*f3*(1.0-5.0*ramass/1215.0)/3.;
591 return 0.00120683*(x-1372.52)*(x-1372.52)+26.2058;
592 }
else if(x <= 1578.0) {
593 return 1.15873e-05*x*x+49965.6/((x-1519.59)*(x-1519.59)+2372.55);
594 }
else if(x <= 2028.4) {
595 return 34.0248+43262.2/((x-1681.65)*(x-1681.65)+1689.35);
596 }
else if(x <= 7500.0) {
597 return 3.3e-7*(x-7500.0)*(x-7500.0)+24.5;
606 return NNTot(p1, p2);
617 return inelastic +
elastic(p1, p2);
638 if(pLab>212677. || pLab<296.367)
643 const G4int cg = 4 + ind2t3*ipit3;
654 xpipp=17.965*std::pow(p1, 5.4606);
656 xpipp=24.3-12.3*std::pow(p1, -1.91)+0.324*p2*p2-2.44*
p2;
666 xpimp=26.6-7.18*std::pow(p1, -1.86)+0.327*p2*p2-2.81*
p2;
673 return 0.5*(xpipp+xpimp);
680 if(x>20000.)
return 0.0;
689 }
else if(particle2->
isPion()) {
695 const G4double q2=(y-1076.0*1076.0)*(y-800.0*800.0)/y/4.0;
699 const G4double q3 = std::pow(std::sqrt(q2),3);
701 G4double sdelResult = 326.5/(std::pow((x-1215.0-ramass)*2.0/(110.0-ramass), 2)+1.0);
702 sdelResult = sdelResult*(1.0-5.0*ramass/1215.0);
703 const G4int cg = 4 + ind2t3*ipit3;
704 sdelResult = sdelResult*f3*cg/6.0;
726 }
else if(particle2->
isPion()) {
734 if((ind2t3 == 1 && ipit3 == 2) || (ind2t3 == -1 && ipit3 == -2))
736 else if((ind2t3 == 1 && ipit3 == -2) || (ind2t3 == -1 && ipit3 == 2))
748 if(isospin==4 || isospin==-4)
return 0.0;
762 if(Ecm <= 938.3 + deltaMass) {
766 if(Ecm < 938.3 + deltaMass + 2.0) {
767 Ecm = 938.3 + deltaMass + 2.0;
785 G4double result = 0.5 * x * y * sDelta;
790 result *= 3.*(32.0 + isospin * isospin * (deltaIsospin * deltaIsospin - 5))/64.0;
791 result /= 1.0 + 0.25 * (isospin * isospin);
826 return 5.5e-6 * x/(7.7 + x);
828 return (5.34 + 0.67*(x - 2.0)) * 1.0e-6;
832 G4double b = (7.16 - 1.63*x) * 1.0e-6;
833 return b/(1.0 + std::exp(-(x - 0.45)/0.05));
834 }
else if(pl < 1100.0) {
835 return (9.87 - 4.88 * x) * 1.0e-6;
837 return (3.68 + 0.76*x) * 1.0e-6;
850 return piNOnePi(particle1,particle2);
852 return piNTwoPi(particle1,particle2);
880 const G4int cg = 4 + ind2*ipi;
890 tamp6=
piNIne(particle1, particle2);
892 tamp6=0.204+18.2*std::pow(p1, -1.72)+6.33*std::pow(p1, -1.13);
899 tamp2=
piNIne(particle1, particle2);
901 tamp2=9.04*std::pow(p1, -1.17)+18.*std::pow(p1, -1.21);
902 if (tamp2 < 0.0) tamp2=0;
910 if (s1pin > inelastic)
941 const G4int cg = 4 + ind2*ipi;
951 tamp6=
piNIne(particle1, particle2)-
piNOnePi(particle1, particle2);
953 tamp6=1.59+25.5*std::pow(p1, -1.04);
961 tamp2=
piNIne(particle1, particle2)-
piNOnePi(particle1, particle2);
963 tamp2=2.457794117647+18.066176470588*std::pow(p1, -0.92);
969 const G4double s2pin=0.5*(tamp6+tamp2);
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.
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.
static const G4double s12mzOOT
One over threshold for s12mz.
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.
G4bool isDelta() const
Is it a Delta?
std::string print() const
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.
virtual G4double piNTwoPi(Particle const *const p1, Particle const *const p2)
Cross section for Two (more) pion production - piN entrance channel.
const HornerC4 s01pzHC
Horner coefficients for s01pz.
CrossSectionsMultiPions()
G4bool nucleon(G4int ityp)
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.
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
Get the particle type.
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.
static const G4double s11pzOOT
One over threshold for s11pz.
G4bool isNucleon() const
Is this a nucleon?
const HornerC4 s12mzHC
Horner coefficients for s12mz.
static const G4double s02pmOOT
One over threshold for s02pm.
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)
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
virtual G4double piNOnePi(Particle const *const p1, Particle const *const p2)
Cross section for One (more) pion production - piN entrance channel.
Cross sections used in INCL Multipions.
const HornerC6 s02pmHC
Horner coefficients for s02pm.
const HornerC5 s12pmHC
Horner coefficients for s12pm.