40 #include <TLorentzVector.h> 71 using std::ostringstream;
72 using namespace genie;
79 int pdgc,
const TLorentzVector & x4,
const TLorentzVector &
p4,
80 double A,
double Z,
double nRpi,
double nRnuc,
const bool useOset,
const bool altOset,
const bool xsecNNCorr,
string INukeMode)
94 bool is_kaon = pdgc ==
kPdgKP;
97 if(!is_pion && !is_nucleon && !is_kaon && !is_gamma)
return 0.;
108 double ring = (momentum>0) ? 1.240/momentum : 0;
110 if(A<=20) { ring /= 2.; }
117 if(INukeMode==
"hN2018")
119 if (is_pion ) { ring *= nRpi; }
120 else if (is_nucleon ) { ring *= nRnuc; }
121 else if (is_gamma || is_kaon || useOset) { ring = 0.;}
125 if (is_pion || is_kaon ) { ring *= nRpi; }
126 else if (is_nucleon ) { ring *= nRnuc; }
127 else if (is_gamma ) { ring = 0.; }
131 double rnow = x4.Vect().Mag();
137 double ke = (p4.Energy() - p4.M()) /
units::MeV;
144 double ppcnt = (double) Z/ (
double)
A;
147 if (is_pion and (INukeMode ==
"hN2018") and useOset and ke < 350.0)
150 { sigtot = fHadroData2018 -> XSecPipp_Tot() -> Evaluate(ke)*ppcnt;
151 sigtot+= fHadroData2018 -> XSecPipn_Tot() -> Evaluate(ke)*(1-ppcnt);}
153 { sigtot = fHadroData2018 -> XSecPi0p_Tot() -> Evaluate(ke)*ppcnt;
154 sigtot+= fHadroData2018 -> XSecPi0n_Tot() -> Evaluate(ke)*(1-ppcnt);}
156 { sigtot = fHadroData2018 -> XSecPipn_Tot() -> Evaluate(ke)*ppcnt;
157 sigtot+= fHadroData2018 -> XSecPipp_Tot() -> Evaluate(ke)*(1-ppcnt);}
160 sigtot = fHadroData2018 -> XSecPp_Tot() -> Evaluate(ke)*ppcnt;
165 double R0 = 1.25 * TMath::Power(A,1./3.) + 2.0 * 0.65;
166 double Mp = pLib->
Find(2212)->Mass();
167 double M = pLib->
Find(pdgc)->Mass();
170 if (Z*hc/137./x4.Vect().Mag() >
E)
173 double Bc = z*Z*hc/137./
R0;
175 double f = TMath::ACos(TMath::Power(x,0.5)) - TMath::Power(x*(1-x),0.5);
176 double B = 0.63*z*Z*TMath::Power((M/Mp)/E,0.5);
177 double Pc = TMath::Exp(-B*f);
180 sigtot+= fHadroData2018 -> XSecPn_Tot() -> Evaluate(ke)*(1-ppcnt);
182 double E0 = TMath::Power(A,0.2)*12.;
183 if (INukeMode==
"hN2018"){
if(ke<E0){sigtot=0.0;}}
188 sigtot = fHadroData2018 -> XSecPn_Tot() -> Evaluate(ke)*ppcnt;
189 sigtot+= fHadroData2018 -> XSecNn_Tot() -> Evaluate(ke)*(1-ppcnt);
190 double E0 = TMath::Power(A,0.2)*12.;
191 if (INukeMode==
"hN2018"){
if(ke<E0){sigtot=0.0;}}
195 { sigtot = fHadroData2018 -> XSecKpN_Tot() -> Evaluate(ke);
199 { sigtot = fHadroData2018 -> XSecGamp_fs() -> Evaluate(ke)*ppcnt;
200 sigtot+= fHadroData2018 -> XSecGamn_fs() -> Evaluate(ke)*(1-ppcnt);}
212 if (xsecNNCorr and is_nucleon)
217 if(sigtot<1
E-6){sigtot=1
E-6;}
220 double lamda = 1. / (rho * sigtot);
223 if( ! TMath::Finite(lamda) ) {
236 int pdgc,
const TLorentzVector & x4,
const TLorentzVector &
p4,
double A)
251 if(!is_deltapp)
return 0.;
254 double rnow = x4.Vect().Mag();
259 double ke = (p4.Energy() - p4.M()) /
units::MeV;
260 ke = TMath::Max(1500., ke);
261 ke = TMath::Min( 0., ke);
266 else if (ke<1000) sig=40;
273 double lamda = 1. / (rho * sig);
276 if( ! TMath::Finite(lamda) ) {
284 int pdgc,
const TLorentzVector & x4,
const TLorentzVector &
p4,
double A,
double Z,
285 double mfp_scale_factor,
double nRpi,
double nRnuc,
double NR,
double R0)
304 double R = NR * R0 * TMath::Power(A, 1./3.);
306 TVector3 dr3 = p4.Vect().Unit();
307 TLorentzVector dr4(dr3,0);
310 <<
"Calculating survival probability for hadron with PDG code = " << pdgc
311 <<
" and momentum = " << p4.P() <<
" GeV";
313 <<
"mfp scale = " << mfp_scale_factor
314 <<
", nRpi = " << nRpi <<
", nRnuc = " << nRnuc <<
", NR = " << NR
315 <<
", R0 = " << R0 <<
" fm";
317 TLorentzVector x4_curr(x4);
320 double rnow = x4_curr.Vect().Mag();
321 if (rnow > (R+step))
break;
323 x4_curr += (step*dr4);
324 rnow = x4_curr.Vect().Mag();
327 double mfp_twk = mfp * mfp_scale_factor;
329 double dprob = (mfp_twk>0) ? TMath::Exp(-step/mfp_twk) : 0.;
339 LOG(
"INukeUtils",
pDEBUG) <<
"Psurv = " << prob;
345 const TLorentzVector & x4,
const TLorentzVector &
p4,
346 double A,
double NR,
double R0)
352 double R = NR * R0 * TMath::Power(A, 1./3.);
355 TVector3 dr3 = p4.Vect().Unit();
356 TLorentzVector dr4(dr3,0);
358 TLorentzVector x4_curr(x4);
362 double rnow = x4_curr.Vect().Mag();
363 x4_curr += (step*dr4);
365 rnow = x4_curr.Vect().Mag();
372 int pdgc,
const TLorentzVector & x4,
const TLorentzVector &
p4,
373 double A,
double Z,
double NR,
double R0)
382 double R = NR * R0 * TMath::Power(A, 1./3.);
385 TVector3 dr3 = p4.Vect().Unit();
386 TLorentzVector dr4(dr3,0);
388 TLorentzVector x4_curr(x4);
393 double rnow = x4_curr.Vect().Mag();
394 x4_curr += (step*dr4);
396 rnow = x4_curr.Vect().Mag();
399 d_mfp += (step/lambda);
417 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 419 <<
"Stepping particle [" << p->
Name() <<
"] by dr = " << step <<
" fm";
423 TVector3 dr = p->
P4()->Vect().Unit();
426 TLorentzVector dx4(dr,dt);
427 TLorentzVector x4new = *(p->
X4()) + dx4;
429 if(nuclear_radius > 0.) {
434 double r = x4new.Vect().Mag();
435 double rmax = nuclear_radius+
epsilon;
438 <<
"Particle was stepped too far away (r = " << r <<
" fm)";
440 <<
"Placing it " << epsilon
441 <<
" fm outside the nucleus (r' = " << rmax <<
" fm)";
447 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 464 int &RemnA,
int &RemnZ, TLorentzVector &RemnP4,
469 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 471 <<
"PreEquilibrium() is invoked for a : " << p->
Name()
472 <<
" whose kinetic energy is : " << p->
KinE();
479 bool allow_dup =
true;
482 double ppcnt = (double) RemnZ / (
double) RemnA;
491 if(rnd->
RndFsi().Rndm()<ppcnt)
500 ppcnt = (double) RemnZ / (
double) RemnA;
531 if(success)
LOG(
"INukeUtils2018",
pINFO) <<
"Successful phase space decay for pre-equilibrium nucleus FSI event";
535 exception.
SetReason(
"Phase space generation of pre-equilibrium nucleus final state failed, details above");
555 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 557 <<
"Particle at: " << p_loc;
564 int f_loc = p_loc + 1;
584 for(
unsigned int j=0;
j<descendants->size();
j++)
586 loc = (*descendants)[
j];
615 int &RemnA,
int &RemnZ, TLorentzVector &RemnP4,
620 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 622 <<
"Equilibrium() is invoked for a : " << p->
Name()
623 <<
" whose kinetic energy is : " << p->
KinE();
630 bool allow_dup =
true;
634 double ppcnt = (double) RemnZ / (
double) RemnA;
644 if(rnd->
RndFsi().Rndm()<ppcnt)
653 ppcnt = (double) RemnZ / (
double) RemnA;
656 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 658 <<
"Remnant nucleus (A,Z) = (" << RemnA <<
", " << RemnZ <<
")";
688 if (success)
LOG(
"INukeUtils",
pINFO) <<
"successful equilibrium interaction";
692 exception.
SetReason(
"Phase space generation of compound nucleus final state failed, details above");
701 GHepRecord* ev,
int pcode,
int tcode,
int scode,
int s2code,
double C3CM,
715 double E3L, P3L, E4L, P4L;
716 TVector3 tP1L, tPtot, tbeta, tbetadir, tTrans, tVect;
717 TVector3 tP1zCM, tP2zCM, tP3L, tP4L;
729 M1 = pLib->
Find(pcode)->Mass();
731 M3 = pLib->
Find(scode)->Mass();
732 M4 = pLib->
Find(s2code)->Mass();
735 TLorentzVector t4P1L = *p->
P4();
736 TLorentzVector t4P2L = *t->
P4();
739 double bindE = 0.025;
742 LOG(
"TwoBodyCollision",
pNOTICE) <<
"M1 = " << t4P1L.M() <<
" , M2 = " << t4P2L.M();
743 LOG(
"TwoBodyCollision",
pNOTICE) <<
"E1 = " << t4P1L.E() <<
" , E2 = " << t4P2L.E();
745 if ( (p->
Energy()-p->
Mass()) < bindE ){bindE = 0.;}
748 if((pcode==2112||pcode==2212)&&(t4P1L.E()-M1)<.1) bindE = 0.0;
749 if((pcode==211||pcode==-211||pcode==111)&&(t4P1L.E()-M1)<.08) bindE = 0.0;
750 if((pcode==321)&&(t4P1L.E()-M1)<.1) bindE = 0.0;
753 TLorentzVector t4P3L, t4P4L;
756 P3L = t4P3L.Vect().Mag();
757 P4L = t4P4L.Vect().Mag();
762 <<
"TwoBodyKinematics fails: C3CM, P3 = " << C3CM <<
" " 763 << P3L <<
" " << E3L <<
"\n" <<
" P4 = " 764 << P4L <<
" " << E4L ;
769 P3L = t4P3L.Vect().Mag();
770 P4L = t4P4L.Vect().Mag();
774 <<
"C3CM, P3 = " << C3CM <<
" " 775 << P3L <<
" " << E3L <<
"\n" <<
" P4 = " 776 << P4L <<
" " << E4L ;
779 if(!(TMath::Finite(P3L)) || P3L<.001)
782 <<
"Particle 3 momentum small or non-finite: " << P3L
783 <<
"\n" <<
"--> Assigning .001 as new momentum";
787 if(!(TMath::Finite(P4L)) || P4L<.001)
790 <<
"Particle 4 momentum small or non-finite: " << P4L
791 <<
"\n" <<
"--> Assigning .001 as new momentum";
812 <<
"t4P2L= " << t4P2L.E() <<
" " << t4P2L.Z()
813 <<
" RemnP4= " << RemnP4.E() <<
" " << RemnP4.Z() ;
839 LOG(
"INukeUtils",
pINFO) <<
"Successful 2 body collision";
845 double M3,
double M4, TLorentzVector t4P1L, TLorentzVector t4P2L,
846 TLorentzVector &t4P3L, TLorentzVector &t4P4L,
double C3CM, TLorentzVector &RemnP4,
double bindE)
866 double E1L, E2L, P1L, P2L, E3L, P3L;
870 double E1CM, E2CM, E3CM, P3CM;
873 double theta1, theta2, theta5, P1zL, P2zL, P1tL, P2tL;
874 TVector3 tbeta, tbetadir, tTrans, tVect;
875 TVector3 tP1zCM, tP2zCM, vP3L;
876 TLorentzVector t4P1buf, t4P2buf, t4Ptot;
885 if (C3CM < -1. || C3CM > 1.)
return false;
899 t4Ptot = t4P1buf + t4P2buf;
901 LOG(
"INukeUtils",
pINFO) <<
"M1 "<<t4P1L.M()<<
", M2 "<<t4P2L.M();
902 LOG(
"INukeUtils",
pINFO) <<
"E1L "<<E1L<<
", E1CM "<<E1CM;
903 LOG(
"INukeUtils",
pINFO) <<
"bindE = " << bindE;
913 LOG(
"INukeUtils",
pNOTICE) <<
"TwoBodyKinematics Failed: Forbidden by binding energy";
914 LOG(
"INukeUtils",
pNOTICE) <<
"E1L, E2L, M3, M4 : "<< E1L <<
", "<< E2L <<
", "<< M3 <<
", "<< M4;
915 t4P3L.SetPxPyPzE(0,0,0,0);
916 t4P4L.SetPxPyPzE(0,0,0,0);
922 tbeta = t4Ptot.Vect() * (1.0 / (E1L + E2L));
923 tbetadir = tbeta.Unit();
928 theta1 = t4P1buf.Angle(tbeta);
929 theta2 = t4P2buf.Angle(tbeta);
930 P1zL = P1L*TMath::Cos(theta1);
931 P2zL = P2L*TMath::Cos(theta2);
935 if(TMath::Abs((tVect - tbetadir).
Mag())<.01) tVect.SetXYZ(0,1,0);
936 theta5 = tVect.Angle(tbetadir);
937 tTrans = (tVect - TMath::Cos(theta5)*tbetadir).
Unit();
940 E1CM = gm*E1L - gm*beta*P1zL;
941 tP1zCM = gm*P1zL*tbetadir - gm*tbeta*E1L;
942 E2CM = gm*E2L - gm*beta*P2zL;
943 tP2zCM = gm*P2zL*tbetadir - gm*tbeta*E2L;
947 LOG(
"INukeUtils",
pINFO) <<
"M1 "<<t4P1L.M()<<
", M2 "<<t4P2L.M();
948 LOG(
"INukeUtils",
pINFO) <<
"E1L "<<E1L<<
", E1CM "<<E1CM;
949 LOG(
"INukeUtils",
pINFO) <<
"P1zL "<<P1zL<<
", P1zCM "<<tP1zCM.Mag()<<
", P1tL "<<P1tL;
950 LOG(
"INukeUtils",
pINFO) <<
"E2L "<<E2L<<
", E2CM "<<E2CM;
951 LOG(
"INukeUtils",
pINFO) <<
"P2zL "<<P2zL<<
", P2zCM "<<tP2zCM.Mag()<<
", P2tL "<<P2tL;
952 LOG(
"INukeUtils",
pINFO) <<
"C3CM "<<C3CM;
955 E3CM = (Et*Et + M3*M3 - M4*M4) / (2.0 * Et);
958 if(E3CM*E3CM - M3*M3<0 || E3CM<0 || Et<0)
960 if (Et<0)
LOG(
"INukeUtils",
pNOTICE) <<
"TwoBodyKinematics Failed: Total energy is negative";
961 if (E3CM<M3)
LOG(
"INukeUtils",
pNOTICE) <<
"TwoBodyKinematics Failed: Scattered Particle 3 CM energy is too small";
962 if (E3CM*E3CM - M3*M3<0)
LOG(
"INukeUtils",
pNOTICE) <<
"TwoBodyKinematics Failed: Scattered Particle 3 CM momentum is nonreal";
963 t4P3L.SetPxPyPzE(0,0,0,0);
964 t4P4L.SetPxPyPzE(0,0,0,0);
970 P3zL = gm*beta*E3CM + gm*P3CM*C3CM;
978 double E4CM = Et-E3CM;
979 double P4zL = gm*beta*E4CM - gm*P3CM*C3CM;
980 double P4tL = -1.*P3tL;
984 LOG(
"INukeUtils",
pINFO) <<
"M3 "<< M3 <<
", M4 "<< M4;
985 LOG(
"INukeUtils",
pINFO) <<
"E3L "<<E3L<<
", E3CM "<<E3CM;
986 LOG(
"INukeUtils",
pINFO) <<
"P3zL "<<P3zL<<
", P3tL "<<P3tL;
987 LOG(
"INukeUtils",
pINFO) <<
"C3L "<<P3zL/P3L;
989 LOG(
"INukeUtils",
pINFO) <<
"E4L "<<E4L<<
", E4CM "<<E4CM;
990 LOG(
"INukeUtils",
pINFO) <<
"P4zL "<<P4zL<<
", P4tL "<<P4tL;
991 LOG(
"INukeUtils",
pINFO) <<
"P4L "<<P4L;
992 LOG(
"INukeUtils",
pINFO) <<
"C4L "<<P4zL/P4L;
994 double echeck = E1L + E2L - (E3L + E4L);
995 double pzcheck = P1zL+ P2zL - (P3zL + P4zL);
996 double ptcheck = P1tL+ P2tL - (P3tL + P4tL);
998 LOG(
"INukeUtils",
pINFO) <<
"Check 4-momentum conservation - Energy "<<echeck<<
", z momentum "<<pzcheck <<
", transverse momentum " << ptcheck ;
1003 if(!(TMath::Finite(P3L)) || P3L<.001)
1006 <<
"Particle 3 momentum small or non-finite: " << P3L
1007 <<
"\n" <<
"--> Assigning .001 as new momentum";
1017 vP3L = P3zL*tbetadir + P3tL*tTrans;
1018 vP3L.Rotate(PHI3,tbetadir);
1020 t4P3L.SetVect(vP3L);
1023 t4P4L = t4P1buf + t4P2buf - t4P3L;
1024 t4P4L-= TLorentzVector(0,0,0,bindE);
1031 if(t4P4L.Mag2()<0 || t4P4L.E()<0)
1033 LOG(
"INukeUtils",
pNOTICE)<<
"TwoBodyKinematics Failed: Target mass or energy is negative";
1034 t4P3L.SetPxPyPzE(0,0,0,0);
1035 t4P4L.SetPxPyPzE(0,0,0,0);
1039 if (bindE!=0) RemnP4 += TLorentzVector(0,0,0,bindE);
1045 bool DoFermi,
double FermiFac,
double FermiMomentum,
const NuclearModelI* Nuclmodel)
1057 double M1, M2, M3, M4, M5;
1058 double P1L, P2L, P3L, P4L, P5L;
1059 double E1L, E2L, E3L, E4L, E5L;
1060 double E1CM, E2CM, P3tL;
1061 double PizL, PitL, PiL, EiL;
1062 double EiCM, P4CM2, E4CM2, E5CM2, P3CM, E3CM;
1063 double beta,
gm, beta2, gm2;
1064 double P3zL, P4zL, P4tL, P5zL, P5tL;
1065 double Et, M, theta1, theta2;
1067 double theta3, theta4, phi3, phi4, theta5;
1068 TVector3 tP2L, tP1L, tPtot, tbeta, tbetadir, tTrans, tP4L, tP5L;
1069 TVector3 tP1zCM, tP2zCM, tP3L, tPiL, tbeta2, tbetadir2, tVect, tTrans2;
1078 M2 = pLib->
Find(tcode)->Mass();
1089 target.SetHitNucPdg(tcode);
1091 tP2L = FermiFac * Nuclmodel->
Momentum3();
1097 tP2L.SetXYZ(0.0, 0.0, 0.0);
1106 tP1L = p->
P4()->Vect();
1107 tPtot = tP1L + tP2L;
1109 tbeta = tPtot * (1.0 / (E1L + E2L));
1110 tbetadir = tbeta.Unit();
1114 theta1 = tP1L.Angle(tbeta);
1115 theta2 = tP2L.Angle(tbeta);
1116 P1zL = P1L*TMath::Cos(theta1);
1117 P2zL = P2L*TMath::Cos(theta2);
1118 tVect.SetXYZ(1,0,0);
1119 if(TMath::Abs((tVect - tbetadir).
Mag())<.01) tVect.SetXYZ(0,1,0);
1120 theta5 = tVect.Angle(tbetadir);
1121 tTrans = (tVect - TMath::Cos(theta5)*tbetadir).
Unit();
1123 E1CM = gm*E1L - gm*beta*P1zL;
1124 tP1zCM = gm*P1zL*tbetadir - gm*tbeta*E1L;
1125 E2CM = gm*E2L - gm*beta*P2zL;
1126 tP2zCM = gm*P2zL*tbetadir - gm*tbeta*E2L;
1128 M = (rnd->
RndFsi().Rndm()*(Et - M3 - M4 - M5)) + (M4 + M5);
1129 E3CM = (Et*Et + M3*M3 - M*M)/(2*Et);
1131 if(E3CM*E3CM - M3*M3<0)
1134 <<
"PionProduction P3 has non-real momentum - retry kinematics";
1135 LOG(
"INukeUtils",
pNOTICE) <<
"Energy, masses of 3 fs particales:" 1136 << E3CM <<
" " << M3 <<
" " <<
" " << M4 <<
" " << M5;
1138 exception.
SetReason(
"PionProduction particle 3 has non-real momentum");
1149 P3zL = gm*beta*E3CM + gm*P3CM*TMath::Cos(theta3);
1150 P3tL = P3CM*TMath::Sin(theta3);
1151 PizL = gm*beta*EiCM - gm*P3CM*TMath::Cos(theta3);
1152 PitL = -P3CM*TMath::Sin(theta3);
1160 if(!(TMath::Finite(P3L)) || P3L < .001)
1163 <<
"Particle 3 " << M3 <<
" momentum small or non-finite: " << P3L
1164 <<
"\n" <<
"--> Assigning .001 as new momentum";
1171 tP3L = P3zL*tbetadir + P3tL*tTrans;
1172 tPiL = PizL*tbetadir + PitL*tTrans;
1173 tP3L.Rotate(phi3,tbetadir);
1174 tPiL.Rotate(phi3,tbetadir);
1177 tbeta2 = tPiL * (1.0 / EiL);
1178 tbetadir2 = tbeta2.Unit();
1179 beta2 = tbeta2.Mag();
1182 E4CM2 = (M*M + M4*M4 - M5*M5) / (2*M);
1186 tVect.SetXYZ(1,0,0);
1187 if(TMath::Abs((tVect - tbetadir2).
Mag())<.01) tVect.SetXYZ(0,1,0);
1188 theta5 = tVect.Angle(tbetadir2);
1189 tTrans2 = (tVect - TMath::Cos(theta5)*tbetadir2).
Unit();
1191 P4zL = gm2*beta2*E4CM2 + gm2*P4CM2*TMath::Cos(theta4);
1192 P4tL = P4CM2*TMath::Sin(theta4);
1193 P5zL = gm2*beta2*E5CM2 - gm2*P4CM2*TMath::Cos(theta4);
1202 if(!(TMath::Finite(P4L)) || P4L < .001)
1205 <<
"Particle 4 " << M4 <<
" momentum small or non-finite: " << P4L
1206 <<
"\n" <<
"--> Assigning .001 as new momentum";
1212 if(!(TMath::Finite(P5L)) || P5L < .001)
1215 <<
"Particle 5 " << M5 <<
" momentum small or non-finite: " << P5L
1216 <<
"\n" <<
"--> Assigning .001 as new momentum";
1223 tP4L = P4zL*tbetadir2 + P4tL*tTrans2;
1224 tP5L = P5zL*tbetadir2 + P5tL*tTrans2;
1225 tP4L.Rotate(phi4,tbetadir2);
1226 tP5L.Rotate(phi4,tbetadir2);
1232 <<
"PionProduction fails because of Pauli blocking - retry kinematics";
1234 exception.
SetReason(
"PionProduction final state not determined");
1248 TLorentzVector(tP3L,E3L);
1249 TLorentzVector(tP4L,E4L);
1250 TLorentzVector(tP5L,E5L);
1274 TLorentzVector &RemnP4,
bool DoFermi,
double FermiFac,
double FermiMomentum,
const NuclearModelI* Nuclmodel)
1295 int pcode = p->
Pdg();
1297 int p1code = p->
Pdg();
1299 int p3code = 0, p4code = 0, p5code = 0;
1315 double kine = 1000*p->
KinE();
1323 double xsec2pi0n = TMath::Max(0.,TMath::Exp(8.82 - (573.2/(kine-107.3))));
1326 double xsecpippimn = TMath::Max(0.,TMath::Exp(11.06 - (985.9/(kine-88.2))));
1329 double xsecpimpi0p = TMath::Max(0.,TMath::Exp(9.58 - (1229.4/(kine-60.5))));
1330 double totpimp = xsec2pi0n + xsecpippimn + xsecpimpi0p;
1336 double xsec2pipn = TMath::Max(0.,TMath::Exp(5.64 - (222.6/(kine-150.))));
1339 double xsecpippi0p = TMath::Max(0.,TMath::Exp(7.95 - (852.6/(kine-77.8))));
1340 double totpipp = xsec2pipn + xsecpippi0p;
1342 if (totpimp<=0 && totpipp<=0) {
1343 LOG(
"INukeUtils",
pNOTICE) <<
"InelasticHN called below threshold energy";
1349 double xsecp, xsecn;
1351 case kPdgPi0: xsecp = 0.5 * (totpimp + totpipp); xsecn = xsecp;
break;
1352 case kPdgPiP: xsecp = totpipp; xsecn = totpimp;
break;
1353 case kPdgPiM: xsecp = totpimp; xsecn = totpipp;
break;
1355 LOG(
"INukeUtils",
pWARN) <<
"InelasticHN cannot handle probe: " 1358 exception.
SetReason(
"PionProduction final state not determined");
1367 xsecn *= RemnA-RemnZ;
1371 double rand = rnd->
RndFsi().Rndm() * (xsecp + xsecn);
1373 { rand /= RemnZ; ptarg =
true;}
1375 { rand -= xsecp; rand /= RemnA-RemnZ; ptarg =
false;}
1380 if (((ptarg==
true)&&(p1code==
kPdgPiP))
1381 || ((ptarg==
false)&&(p1code==
kPdgPiM)))
1383 if (rand < xsec2pipn)
1395 else if (((ptarg==
false)&&(p1code==
kPdgPiP))
1396 || ((ptarg==
true)&&(p1code==
kPdgPiM)))
1398 if (rand < xsec2pi0n)
1404 else if (rand < (xsec2pi0n + xsecpippimn))
1419 rand = rnd->
RndFsi().Rndm();
1420 if (rand < 191./270.)
1426 else if (rand < 7./135.)
1443 exception.
SetReason(
"PionProduction final state not determined");
1451 double tote = p->
Energy();
1452 double pMass = pLib->
Find(2212)->Mass();
1453 double nMass = pLib->
Find(2112)->Mass();
1454 double etapp2ppPi0 =
1456 double etapp2pnPip =
1458 pMass+nMass,pLib->
Find(211)->Mass());
1459 double etapn2nnPip =
1461 double etapn2ppPim =
1464 if ((etapp2ppPi0<=0.)&&(etapp2pnPip<=0.)&&(etapn2nnPip<=0.)&&(etapn2ppPim<=0.)) {
1465 LOG(
"INukeUtils",
pNOTICE) <<
"PionProduction() called below threshold energy";
1467 exception.
SetReason(
"PionProduction final state not possible - below threshold");
1473 double xsecppPi0=0,xsecpnPiP=0,xsecnnPiP=0,xsecppPiM=0;
1475 xsecppPi0 = 4511*(1-.91*TMath::Exp(-TMath::Power((etapp2ppPi0-.705),2)));
1476 xsecppPi0 *= (1-TMath::Exp(-TMath::Power((.556*etapp2ppPi0),3.5)));
1477 xsecppPi0 *= (1-TMath::Exp(-56.897/(etapp2ppPi0*etapp2ppPi0)));
1478 xsecppPi0 = TMath::Max(0.,xsecppPi0);}
1481 xsecpnPiP = 18840*(1-.808*TMath::Exp(-TMath::Power((etapp2pnPip-.371),2)));
1482 xsecpnPiP *= (1-TMath::Exp(-TMath::Power((.568*etapp2pnPip),3.2)));
1483 xsecpnPiP *= (1-TMath::Exp(-39.818/(etapp2pnPip*etapp2pnPip)));
1484 xsecpnPiP = TMath::Max(0.,xsecpnPiP);}
1487 xsecnnPiP = 7670*(1-.479*TMath::Exp(-TMath::Power((etapn2nnPip-.947),2)));
1488 xsecnnPiP *= (1-TMath::Exp(-TMath::Power((.35*etapn2nnPip),3.2)));
1489 xsecnnPiP *= (1-TMath::Exp(-71.53/(etapn2nnPip*etapn2nnPip)));
1490 xsecnnPiP = TMath::Max(0.,xsecnnPiP);}
1493 xsecppPiM = 7670*(1-.479*TMath::Exp(-TMath::Power((etapn2ppPim-.947),2)));
1494 xsecppPiM *= (1-TMath::Exp(-TMath::Power((.35*etapn2ppPim),3.2)));
1495 xsecppPiM *= (1-TMath::Exp(-71.53/(etapn2ppPim*etapn2ppPim)));
1496 xsecppPiM = TMath::Max(0.,xsecppPiM);}
1499 double sigma10 = TMath::Max(0.,xsecpnPiP - xsecppPi0);
1500 double sigma01 = TMath::Max(0.,xsecppPiM + xsecnnPiP - xsecppPi0);
1502 double xsecpnPi0 = .5*(sigma10 + sigma01);
1503 xsecpnPi0 = TMath::Max(xsecpnPi0,0.);
1505 LOG(
"INukeUtils",
pDEBUG) <<
'\n' <<
"Cross section values: "<<
'\n' 1506 << xsecppPi0 <<
" PP pi0" <<
'\n' 1507 << xsecpnPiP <<
" PN pi+" <<
'\n' 1508 << xsecnnPiP <<
" NN pi+" <<
'\n' 1509 << xsecpnPi0 <<
" PN pi0";
1511 double xsecp=0,xsecn=0;
1513 case kPdgProton: xsecp=xsecppPi0+xsecpnPiP; xsecn=xsecppPiM+xsecnnPiP+xsecpnPi0;
break;
1514 case kPdgNeutron: xsecp=xsecppPiM+xsecnnPiP+xsecpnPi0; xsecn=xsecppPi0+xsecpnPiP;
break;
1516 LOG(
"INukeUtils",
pWARN) <<
"InelasticHN cannot handle probe: " 1525 xsecn *= RemnA-RemnZ;
1529 double rand = rnd->
RndFsi().Rndm() * (xsecp + xsecn);
1531 { rand /= RemnZ; ptarg =
true;}
1533 { rand -= xsecp; rand /= RemnA-RemnZ; ptarg =
false;}
1567 <<
"Unable to handle probe (=" << p1code <<
") in InelasticHN()";
1574 LOG(
"INukeUtils",
pNOTICE) <<
"PionProduction() failed : no nucleons to produce pions";
1582 LOG(
"INukeUtils",
pNOTICE) <<
"PionProduction() failed : too few protons in nucleus";
1584 exception.
SetReason(
"PionProduction fails - too few protons available");
1610 LOG(
"INukeUtils",
pDEBUG) <<
"Remnant (A,Z) = (" <<RemnA<<
','<<RemnZ<<
')';
1612 RemnP4 -= *s1->
P4() + *s2->
P4() + *s3->
P4() - *p->
P4();
1617 exception.
SetReason(
"PionProduction final state not determined");
1624 double Mtwopart,
double Mpi)
1634 double Scm = Minc*Minc + Mtarg*Mtarg + 2*Mtarg*nrg;
1636 eta= TMath::Power(Scm,2) + TMath::Power(Mtwopart,4) + TMath::Power(Mpi,4);
1637 eta-= 2*TMath::Power(Mtwopart*Mpi,2);
1638 eta-= 2*Scm*TMath::Power(Mtwopart,2);
1639 eta-= 2*Scm*TMath::Power(Mpi,2);
1640 eta = TMath::Power(eta/Scm,1./2.);
1643 eta=TMath::Max(eta,0.);
1657 LOG(
"INukeUtils",
pINFO) <<
"*** Performing a Phase Space Decay";
1660 LOG(
"INukeUtils",
pINFO) <<
"probe mass: M = " << p->
Mass();
1664 ostringstream state_sstream;
1665 state_sstream <<
"( ";
1666 vector<int>::const_iterator pdg_iter;
1668 double * mass =
new double[pdgv.size()];
1669 double mass_sum = 0;
1670 for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
1671 int pdgc = *pdg_iter;
1676 state_sstream << nm <<
" ";
1678 state_sstream <<
")";
1680 TLorentzVector * pd = p->
GetP4();
1687 double availE = pd->Energy() + mass_sum;
1688 if(is_nuc||is_kaon) availE -= p->
Mass();
1692 <<
"size, mass_sum, availE, pd mass, energy = " << pdgv.size() <<
" " 1693 << mass_sum <<
" " << availE <<
" " << p->
Mass() <<
" " << p->
Energy() ;
1696 double dE = mass_sum;
1697 if(is_nuc||is_kaon) dE -= p->
Mass();
1698 TLorentzVector premnsub(0,0,0,dE);
1702 <<
"Final state = " << state_sstream.str() <<
" has N = " << pdgv.size()
1703 <<
" particles / total mass = " << mass_sum;
1708 TGenPhaseSpace GenPhaseSpace;
1709 bool permitted = GenPhaseSpace.SetDecay(*pd, pdgv.size(), mass);
1712 <<
" *** Phase space decay is not permitted \n" 1713 <<
" Total particle mass = " << mass_sum <<
"\n" 1730 for(
int k=0; k<200; k++) {
1731 double w = GenPhaseSpace.Generate();
1732 wmax = TMath::Max(wmax,w);
1737 <<
"Max phase space gen. weight @ current hadronic interaction: " << wmax;
1744 bool accept_decay=
false;
1745 unsigned int itry=0;
1747 while(!accept_decay)
1754 <<
"Couldn't generate an unweighted phase space decay after " 1755 << itry <<
" attempts";
1761 double w = GenPhaseSpace.Generate();
1762 double gw = wmax * rnd->
RndFsi().Rndm();
1766 <<
"Decay weight = " << w <<
" > max decay weight = " << wmax;
1769 LOG(
"INukeUtils",
pNOTICE) <<
"Decay weight = " << w <<
" / R = " << gw;
1770 accept_decay = (gw<=
w);
1778 LOG(
"INukeUtils",
pNOTICE) <<
"mother index = " << mom;
1782 TLorentzVector * v4 = p->
GetX4();
1784 double checkpx = p->
Px();
1785 double checkpy = p->
Py();
1786 double checkpz = p->
Pz();
1787 double checkE = p->
E();
1791 for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) {
1794 int pdgc = *pdg_iter;
1798 TLorentzVector * p4fin = GenPhaseSpace.GetDecay(i++);
1805 double En = p4fin->Energy();
1812 double dE_leftover = TMath::Min(NucRmvE, KE);
1815 double pmag_old = p4fin->P();
1816 double pmag_new =
TMath::Sqrt(TMath::Max(0.,En*En-M*M));
1817 double scale = pmag_new / pmag_old;
1818 double pxn = scale * p4fin->Px();
1819 double pyn = scale * p4fin->Py();
1820 double pzn = scale * p4fin->Pz();
1822 TLorentzVector p4n(pxn,pyn,pzn,En);
1833 if (p4n.Vect().Mag()>=0.001)
1835 GHepParticle new_particle(pdgc, ist_pi, mom,-1,-1,-1, p4n, *v4);
1843 LOG(
"INukeUtils",
pINFO)<<
"Momentum too small; assigning 0.001 as new momentum";
1845 double phi = 2*
kPi*rnd->
RndFsi().Rndm();
1846 double omega = 2*rnd->
RndFsi().Rndm();
1850 p4n.SetPxPyPzE(0.001,0,0,E4n);
1851 p4n.Rotate(TMath::ACos(1-omega),TVector3(0,0,1));
1852 p4n.Rotate(phi,TVector3(1,0,0));
1854 RemnP4 -= (p4n - TLorentzVector(0,0,0,M));
1856 GHepParticle new_particle(pdgc, ist, mom,-1,-1,-1, p4n, *v4);
1862 GHepParticle new_particle(pdgc, ist, mom,-1,-1,-1, p4n, *v4);
1868 double dpx = (1-
scale)*p4fin->Px();
1869 double dpy = (1-
scale)*p4fin->Py();
1870 double dpz = (1-
scale)*p4fin->Pz();
1871 TLorentzVector premnadd(dpx,dpy,dpz,dE_leftover);
1875 LOG(
"INukeUtils",
pNOTICE) <<
"check conservation: Px = " << checkpx <<
" Py = " << checkpy
1876 <<
" Pz = " << checkpz <<
" E = " << checkE;
1887 const double &pionKineticEnergy,
1890 const double &protonFraction,
1891 const bool &isTableChosen
1897 if (iNukeOset == NULL)
1902 static const std::string dataDir = (gSystem->Getenv(
"GINUKEHADRONDATA")) ?
1903 string(gSystem->Getenv(
"GINUKEHADRONDATA")) :
1904 string(gSystem->Getenv(
"GENIE")) +
1905 string(
"/data/evgen/intranuke/");
1907 static const std::string dataFile = dataDir +
"tot_xsec/" 1908 "intranuke-xsection-pi+N-Oset.dat";
1917 iNukeOset->
setupOset (density, pionKineticEnergy, pionPDG, protonFraction);
double MeanFreePath_Delta(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A)
Mean free path (Delta++ test)
const int kPdgP33m1232_DeltaPP
An exception thrown by SimulateHadronState for kinematics problems. TwoBodyCollision/Kinematics used ...
void SetFirstMother(int m)
virtual GHepParticle * Particle(int position) const
void PreEquilibrium(GHepRecord *ev, GHepParticle *p, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, const NuclearModelI *Nuclmodel, double NucRmvE, EINukeMode mode=kIMdHN)
virtual int ParticlePosition(int pdg, GHepStatus_t i, int start=0) const
TLorentzVector * GetX4(void) const
TRandom3 & RndFsi(void) const
rnd number generator used by intranuclear cascade monte carlos
static INukeNucleonCorr * getInstance()
get single instance of INukeNucleonCorr; create if necessary
THE MAIN GENIE PROJECT NAMESPACE
double E(void) const
Get energy.
static RandomGen * Instance()
Access instance.
const TLorentzVector * P4(void) const
bool PhaseSpaceDecay(GHepRecord *ev, GHepParticle *p, const PDGCodeList &pdgv, TLorentzVector &RemnP4, double NucRmvE, EINukeMode mode=kIMdHA)
general phase space decay method
double Density(double r, int A, double ring=0.)
Table-based implementation of Oset model.
enum genie::EGHepStatus GHepStatus_t
void StepParticle(GHepParticle *p, double step, double nuclear_radius=-1.)
Step particle.
bool ThreeBodyKinematics(GHepRecord *ev, GHepParticle *p, int tcode, GHepParticle *s1, GHepParticle *s2, GHepParticle *s3, bool DoFermi=false, double FermiFac=0, double FermiMomentum=0, const NuclearModelI *Nuclmodel=(const NuclearModelI *) 0)
bool PionProduction(GHepRecord *ev, GHepParticle *p, GHepParticle *s1, GHepParticle *s2, GHepParticle *s3, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, double FermiMomentum, const NuclearModelI *Nuclmodel)
static constexpr Double_t nm
double ProbSurvival(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A, double Z, double mfp_scale_factor=1.0, double nRpi=0.5, double nRnuc=1.0, double NR=3, double R0=1.4)
Hadron survival probability.
string P4AsString(const TLorentzVector *p)
static const unsigned int kMaxUnweightDecayIterations
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
::xsd::cxx::tree::exception< char > exception
double Mass(Resonance_t res)
resonance mass (GeV)
void SetMomentum(const TLorentzVector &p4)
double Mass(void) const
Mass that corresponds to the PDG code.
A singleton holding random number generator classes. All random number generation in GENIE should tak...
double Pz(void) const
Get Pz.
const TVector3 & Momentum3(void) const
double Energy(void) const
Get energy.
double Px(void) const
Get Px.
bool CompareMomentum(const GHepParticle *p) const
int LastMother(void) const
double Dist2Exit(const TLorentzVector &x4, const TLorentzVector &p4, double A, double NR=3, double R0=1.4)
Distance to exit.
int FirstMother(void) const
string Name(void) const
Name that corresponds to the PDG code.
const int kPdgCompNuclCluster
void SetPosition(const TLorentzVector &v4)
double getTotalCrossSection() const
return total = (qel+cex+abs) cross section
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
void SetReason(string reason)
TLorentzVector * GetP4(void) const
bool ComparePdgCodes(const GHepParticle *p) const
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
bool TwoBodyKinematics(double M3, double M4, TLorentzVector tP1L, TLorentzVector tP2L, TLorentzVector &tP3L, TLorentzVector &tP4L, double C3CM, TLorentzVector &RemnP4, double bindE=0)
bool TwoBodyCollision(GHepRecord *ev, int pcode, int tcode, int scode, int s2code, double C3CM, GHepParticle *p, GHepParticle *t, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, EINukeMode mode=kIMdHA)
Intranuke utility functions.
virtual GHepParticle * TargetNucleus(void) const
virtual vector< int > * GetStableDescendants(int position) const
double CalculateEta(double Minc, double ke, double Mtarg, double Mtwopart, double Mpi)
double MeanFreePath(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A, double Z, double nRpi=0.5, double nRnuc=1.0, const bool useOset=false, const bool altOset=false, const bool xsecNNCorr=false, string INukeMode="XX2018")
Mean free path (pions, nucleons)
void Equilibrium(GHepRecord *ev, GHepParticle *p, int &RemnA, int &RemnZ, TLorentzVector &RemnP4, bool DoFermi, double FermiFac, const NuclearModelI *Nuclmodel, double NucRmvE, EINukeMode mode=kIMdHN)
Misc GENIE control constants.
virtual void setupOset(const double &density, const double &pionTk, const int &pionPDG, const double &protonFraction)=0
use to set up Oset class (assign pion Tk, nuclear density etc)
static double fMinKinEnergy
bool CompareStatusCodes(const GHepParticle *p) const
void SetRemovalEnergy(double Erm)
double KinE(bool mass_from_pdg=false) const
Get kinetic energy.
static PDGLibrary * Instance(void)
static INukeHadroData2018 * Instance(void)
void SetLastMother(int m)
Singleton class to load & serve a TDatabasePDG.
Var Sqrt(const Var &v)
Use to take sqrt of a var.
const TLorentzVector * X4(void) const
void SetStatus(GHepStatus_t s)
bool IsNeutronOrProton(int pdgc)
TParticlePDG * Find(int pdgc)
assert(nhit_max >=nhit_nbins)
virtual bool GenerateNucleon(const Target &) const =0
virtual void AddParticle(const GHepParticle &p)
static double fMaxKinEnergyHN
string X4AsString(const TLorentzVector *x)
double sigmaTotalOset(const double &pionKineticEnergy, const double &density, const int &pionPDG, const double &protonFraction, const bool &isTableChosen=true)
GENIE's GHEP MC event record.
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
double Dist2ExitMFP(int pdgc, const TLorentzVector &x4, const TLorentzVector &p4, double A, double Z, double NR=3, double R0=1.4)
Distance to exit.
string Vec3AsString(const TVector3 *vec)
Root of GENIE utility namespaces.
void push_back(int pdg_code)
double Py(void) const
Get Py.