17 #include <TLorentzVector.h> 18 #include <Math/IFunction.h> 19 #include <Math/Integrator.h> 51 using namespace genie;
96 const InitialState & init_state = interaction -> InitState();
107 double mNucleon = ( mNi + interaction->
RecoilNucleon()->Mass() ) / 2.;
120 TLorentzVector inNucleonMomOnShell( target.
HitNucP4().Vect(),
121 inNucleonOnShellEnergy );
132 TLorentzVector neutrinoMom = *tempNeutrino;
134 TLorentzVector inNucleonMom = target.
HitNucP4();
135 TLorentzVector leptonMom = kinematics.
FSLeptonP4();
136 TLorentzVector outNucleonMom = kinematics.
HadSystP4();
143 double pNf = outNucleonMom.P();
144 if ( pNf < kF )
return 0.;
152 *neutrinoMom.E()*outNucleonMom.E()*leptonMom.E()) / 4.0;
156 double ml2 = TMath::Power(ml, 2);
157 double coulombFactor = 1.0;
158 double plLocal = leptonMom.P();
165 double Vc =
vcr(& target, r);
168 int sign = is_neutrino ? 1 : -1;
169 double El = leptonMom.E();
170 double pl = leptonMom.P();
171 double ElLocal = El - sign*Vc;
173 if ( ElLocal - ml <= 0. ) {
174 LOG(
"Nieves",
pDEBUG) <<
"Event should be rejected. Coulomb effects" 175 <<
" push kinematics below threshold. Returning xsec = 0.0";
183 double KEl = El - ml;
185 LOG(
"Nieves",
pDEBUG) <<
"Outgoing lepton has a very small kinetic energy." 186 <<
" Protecting against near-singularities in the Coulomb correction" 187 <<
" factor by returning xsec = 0.0";
196 coulombFactor= (plLocal * ElLocal) / (pl * El);
209 double q0Tilde = outNucleonMom.E() - inNucleonMomOnShell.E();
219 TVector3 neutrinoMom3 = neutrinoMom.Vect();
220 TVector3 leptonMom3 = leptonMom.Vect();
222 TVector3 inNucleonMom3 = inNucleonMom.Vect();
223 TVector3 outNucleonMom3 = outNucleonMom.Vect();
233 TVector3 leptonMomCoulomb3 = (!
fCoulomb ) ? leptonMom3
234 : plLocal * leptonMom3 * (1. / leptonMom3.Mag());
235 TVector3 q3VecTilde = neutrinoMom3 - leptonMomCoulomb3;
238 TVector3 zvec(0.0, 0.0, 1.0);
239 TVector3 rot = ( q3VecTilde.Cross(zvec) ).
Unit();
241 double angle = zvec.Angle( q3VecTilde );
245 if ( q3VecTilde.Perp() == 0. && q3VecTilde.Z() < 0. ) {
246 rot = TVector3(0., 1., 0.);
253 neutrinoMom3.Rotate(angle,rot);
254 neutrinoMom.SetVect(neutrinoMom3);
256 leptonMom3.Rotate(angle,rot);
257 leptonMom.SetVect(leptonMom3);
259 inNucleonMom3.Rotate(angle,rot);
260 inNucleonMom.SetVect(inNucleonMom3);
261 inNucleonMomOnShell.SetVect(inNucleonMom3);
263 outNucleonMom3.Rotate(angle,rot);
264 outNucleonMom.SetVect(outNucleonMom3);
269 TLorentzVector qP4 = neutrinoMom - leptonMom;
270 TLorentzVector qTildeP4(0., 0., q3VecTilde.Mag(), q0Tilde);
272 double Q2 = -1. * qP4.Mag2();
273 double Q2tilde = -1. * qTildeP4.Mag2();
279 double q2 = -Q2tilde;
283 LOG(
"Nieves",
pWARN) <<
"q2 >= 0, returning xsec = 0.0";
300 double LmunuAnumuResult =
LmunuAnumu(neutrinoMom, inNucleonMomOnShell,
301 leptonMom, qTildeP4, mNucleon, is_neutrino, target,
305 double xsec = Gfactor*coulombFactor*LmunuAnumuResult;
316 <<
", q2 = " << q2 <<
", xsec = " <<
xsec;
326 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 328 <<
"Jacobian for transformation to: " 354 LOG(
"Nieves",
pFATAL) <<
"Splines for the Nieves CCQE model must be" 355 <<
" generated using genie::NewQELXSec";
377 bool prcok = proc_info.
IsWeakCC() && ((isP&&isnub) || (isN&&isnu));
378 if(!prcok)
return false;
399 fCos8c2 = TMath::Power(TMath::Cos(thc), 2);
409 this->
SubAlg(
"FormFactorsAlg") );
415 this->
SubAlg(
"XSec-Integrator") );
431 RgKey nuclkey =
"IntegralNuclearModel";
461 if ( temp_mode ==
"VertexGenerator" ) {
464 else if ( temp_mode ==
"Nieves" ) {
468 LOG(
"Nieves",
pFATAL) <<
"Unrecognized setting \"" << temp_mode
469 <<
"\" requested for the RmaxMode parameter in the" 470 <<
" configuration for NievesQELCCPXSec";
495 double M,
double r,
bool is_neutrino,
bool tgtIsNucleus,
int tgt_pdgc,
496 int A,
int Z,
int N,
bool hitNucIsProton,
double & CN,
double & CT,
double & CL,
497 double & imaginaryU,
double &
t0,
double & r00,
bool assumeFreeNucleon)
const 499 if ( tgtIsNucleus && !assumeFreeNucleon ) {
500 double dq = qTildeP4.Vect().Mag();
501 double dq2 = TMath::Power(dq,2);
502 double q2 = 1 * qTildeP4.Mag2();
504 double hbarc2 = TMath::Power(
fhbarc,2);
511 double rho = rhop + rhon;
514 double fPrime = (0.33*rho/rho0+0.45*(1-rho/rho0))*c0;
520 kF1 = TMath::Power(3*
kPi2*rhop, 1.0/3.0) *
fhbarc;
521 kF2 = TMath::Power(3*
kPi2*rhon, 1.0/3.0) *
fhbarc;
523 kF1 = TMath::Power(3*
kPi2*rhon, 1.0/3.0) *
fhbarc;
524 kF2 = TMath::Power(3*
kPi2*rhop, 1.0/3.0) *
fhbarc;
536 double kF = TMath::Power(1.5*
kPi2*rho, 1.0/3.0) *
fhbarc;
538 std::complex<double> imU(
relLindhardIm(qTildeP4.E(),dq,kF1,kF2,
539 M,is_neutrino,
t0,r00));
541 imaginaryU = imag(imU);
543 std::complex<double> relLin(0,0),udel(0,0);
547 relLin =
relLindhard(qTildeP4.E(),dq,kF,M,is_neutrino,imU);
550 std::complex<double> relLinTot(relLin + udel);
557 (2* TMath::Power((6.25-0.5929)/(6.25-q2),2)*dq2/(q2-0.5929) + 0.63);
564 CN = 1.0/TMath::Power(
abs(1.0-fPrime*relLin/hbarc2),2);
566 CT = 1.0/TMath::Power(
abs(1.0-relLinTot*Vt),2);
567 CL = 1.0/TMath::Power(
abs(1.0-relLinTot*Vl),2);
581 double kFn,
double kFp,
587 double M2 = TMath::Power(M,2);
597 double q2 = TMath::Power(q0,2) - TMath::Power(dq,2);
599 double epsRP = TMath::Max(TMath::Max(M,EF2-q0),a);
610 t0 = 0.5*(EF1+epsRP);
611 r00 = (TMath::Power(EF1,2)+TMath::Power(epsRP,2)+EF1*epsRP)/3.0;
612 std::complex<double>
result(0.0,-M2/2.0/
kPi/dq*(EF1-epsRP));
638 double dqgev,
double kFgev,
double M,
640 std::complex<double> )
const 648 LOG(
"Nieves",
pWARN) <<
"relLindhard() failed";
654 std::complex<double> ImLinRel(
relLindhardIm(q0gev,dqgev,kFgev,kFgev,M,isNeutrino,t0,r00));
656 return(RealLinRel*TMath::Power(
fhbarc,2) + 2.0*ImLinRel);
661 double kf,
double m)
const 663 double q02 = TMath::Power(q0, 2);
664 double qm2 = TMath::Power(qm, 2);
665 double kf2 = TMath::Power(kf, 2);
666 double m2 = TMath::Power(m, 2);
667 double m4 = TMath::Power(m, 4);
671 double q4 = TMath::Power(q2,2);
673 double L1 =
log((kf+ef)/m);
674 std::complex<double> uL2(
675 TMath::Log(TMath::Abs(
677 (ef + q0 -
TMath::Sqrt(m2 + TMath::Power(kf + qm,2))))) +
678 TMath::Log(TMath::Abs(
679 (ef + q0 +
TMath::Sqrt(m2 + TMath::Power(kf - qm,2)))/
680 (ef + q0 +
TMath::Sqrt(m2 + TMath::Power(kf + qm,2))))));
682 std::complex<double> uL3(
683 TMath::Log(TMath::Abs((TMath::Power(2*kf + q0*ds,2)-qm2)/
684 (TMath::Power(2*kf - q0*ds,2)-qm2))) +
685 TMath::Log(TMath::Abs((TMath::Power(kf-ef*ds,2) - (4*m4*qm2)/q4)/
686 (TMath::Power(kf+ef*ds,2) - (4*m4*qm2)/q4))));
688 std::complex<double> RlinrelX(-L1/(16.0*
kPi2)+
689 uL2*(2.0*ef+q0)/(32.0*
kPi2*qm)-
692 return RlinrelX*16.0*
m2;
721 double dq,
double rho,
double kF)
const 723 double q_zero = q0/
fhbarc;
725 double k_fermi = kF/
fhbarc;
733 double fdel_f = 2.13;
738 double q_zero2 = TMath::Power(q_zero, 2);
739 double q_mod2 = TMath::Power(q_mod, 2);
740 double k_fermi2 = TMath::Power(k_fermi, 2);
742 double m2 = TMath::Power(m, 2);
743 double m4 = TMath::Power(m, 4);
744 double mpi2 = TMath::Power(mpi, 2);
745 double mpi4 = TMath::Power(mpi, 4);
747 double fdel_f2 = TMath::Power(fdel_f, 2);
753 double s = m2+q_zero2-q_mod2+
756 if(s>TMath::Power(m+mpi,2)){
758 double qcm =
TMath::Sqrt(TMath::Power(s,2)+mpi4+m4-2.0*(s*mpi2+s*m2+
759 mpi2*m2)) /(2.0*srot);
760 gamma = 1.0/3.0 * 1.0/(4.0*
kPi) * fdel_f2*
761 TMath::Power(qcm,3)/srot*(m+
TMath::Sqrt(m2+TMath::Power(qcm,2)))/mpi2;
763 double sp = m2+q_zero2-q_mod2-
767 if(sp > TMath::Power(m+mpi,2)){
769 double qcmp =
TMath::Sqrt(TMath::Power(sp,2)+mpi4+m4-2.0*(sp*mpi2+sp*m2+
770 mpi2*m2))/(2.0*srotp);
771 gammap = 1.0/3.0 * 1.0/(4.0*
kPi) * fdel_f2*
772 TMath::Power(qcmp,3)/srotp*(m+
TMath::Sqrt(m2+TMath::Power(qcmp,2)))/mpi2;
775 const std::complex<double> iNum(0,1.0);
777 std::complex<double>
z(md/(q_mod*k_fermi)*(q_zero-q_mod2/(2.0*md)
778 -wr +iNum*gamma/2.0));
779 std::complex<double> zp(md/(q_mod*k_fermi)*(-q_zero-q_mod2/(2.0*md)
780 -wr +iNum*gammap/2.0));
782 std::complex<double> pzeta(0.0);
784 pzeta = 2.0/(3.0*
z)+2.0/(15.0*z*z*z);
785 }
else if(
abs(z) < TMath::Power(10.0,-2)){
786 pzeta = 2.0*z-2.0/3.0*z*z*z-iNum*
kPi/2.0*(1.0-z*
z);
788 pzeta = z + (1.0-z*
z) *
log((z+1.0)/(z-1.0))/2.0;
791 std::complex<double> pzetap(0);
793 pzetap = 2.0/(3.0*zp)+2.0/(15.0*zp*zp*zp);
794 }
else if(
abs(zp) < TMath::Power(10.0,-2)){
795 pzetap = 2.0*zp-2.0/3.0*zp*zp*zp-iNum*
kPi/2.0*(1.0-zp*zp);
797 pzetap = zp+ (1.0-zp*zp) *
log((zp+1.0)/(zp-1.0))/2.0;
801 return 2.0/3.0 * rho * md/(q_mod*k_fermi) * (pzeta +pzetap) * fdel_f2 *
817 double c = TMath::Power(A,0.35),
z = 0.54;
832 LOG(
"Nieves",
pFATAL) <<
"Unrecognized setting for fCoulombRmaxMode encountered" 833 <<
" in NievesQELCCPXSec::vcr()";
839 LOG(
"Nieves",
pNOTICE) <<
"Radius greater than maximum radius for coulomb corrections." 840 <<
" Integrating to max radius.";
844 ROOT::Math::IBaseFunctionOneDim *
func =
new 850 double reltol = 1
E-4;
851 int nmaxeval = 100000;
852 ROOT::Math::Integrator ig(*func,ig_type,abstol,reltol,nmaxeval);
853 double result = ig.Integral(0,Rmax);
866 int copy[4] = {input[0],input[1],input[2],input[3]};
867 int permutations = 0;
870 for(
int i=0;
i<4;
i++){
871 for(
int j=
i+1;
j<4;
j++){
873 if(input[
i] == input[
j])
880 for(
int k=j;k>
i;k--){
889 if(permutations % 2 == 0){
900 const TLorentzVector inNucleonMomOnShell,
const TLorentzVector leptonMom,
901 const TLorentzVector qTildeP4,
double M,
bool is_neutrino,
906 int tgt_pdgc = target.
Pdg();
912 const double k[4] = {neutrinoMom.E(),neutrinoMom.Px(),neutrinoMom.Py(),neutrinoMom.Pz()};
913 const double kPrime[4] = {leptonMom.E(),leptonMom.Px(),
914 leptonMom.Py(),leptonMom.Pz()};
916 double q2 = qTildeP4.Mag2();
918 const double q[4] = {qTildeP4.E(),qTildeP4.Px(),qTildeP4.Py(),qTildeP4.Pz()};
921 TMath::Power(q[2],2)+TMath::Power(q[3],2));
923 int sign = (is_neutrino) ? 1 : -1;
934 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 939 double M2 = TMath::Power(M, 2);
940 double FA2 = TMath::Power(FA, 2);
941 double Fp2 = TMath::Power(Fp, 2);
942 double F1V2 = TMath::Power(F1V, 2);
943 double xiF2V2 = TMath::Power(xiF2V, 2);
944 double q02 = TMath::Power(q[0], 2);
945 double dq2 = TMath::Power(dq, 2);
946 double q4 = TMath::Power(q2, 2);
949 double CN=1.,CT=1.,CL=1.,imU=0;
951 tgt_pdgc, A, Z, N, hitNucIsProton, CN, CT, CL, imU,
952 t0, r00, assumeFreeNucleon);
957 if ( !
fRPA || assumeFreeNucleon ) {
963 double tulin[4] = {0.,0.,0.,0.};
964 double rulin[4][4] = { {0.,0.,0.,0.},
975 double t3 = (0.5*q2 + q0*
t0)/dq;
985 double bR = (q4+4.0*r00*q02+4.0*q2*q0*
t0)/(4.0*dq2);
986 double gamma = (aR-bR)/2.0;
987 double delta = (-aR+3.0*bR)/2.0/dq2;
989 double r03 = (0.5*q2*t0 + q0*r00)/dq;
996 rulin[3][3] = gamma+delta*dq2;
1000 tulin[0] = inNucleonMomOnShell.E();
1001 tulin[1] = inNucleonMomOnShell.Px();
1002 tulin[2] = inNucleonMomOnShell.Py();
1003 tulin[3] = inNucleonMomOnShell.Pz();
1005 for(
int i=0;
i<4;
i++)
1006 for(
int j=0;
j<4;
j++)
1007 rulin[
i][
j] = tulin[
i]*tulin[
j];
1011 const int g[4][4] = {{1,0,0,0},{0,-1,0,0},{0,0,-1,0},{0,0,0,-1}};
1012 const std::complex<double> iNum(0,1);
1013 int leviCivitaIndexArray[4];
1014 double imaginaryPart = 0;
1016 std::complex<double>
sum(0.0,0.0);
1018 double kPrimek = k[0]*kPrime[0]-k[1]*kPrime[1]-k[2]*kPrime[2]-k[3]*kPrime[3];
1020 std::complex<double> Lmunu(0.0,0.0),Lnumu(0.0,0.0);
1021 std::complex<double> Amunu(0.0,0.0),Anumu(0.0,0.0);
1026 double axx=0.,azz=0.,a0z=0.,a00=0.,axy=0.;
1027 for(
int mu=0;mu<4;mu++){
1028 for(
int nu=mu;nu<4;nu++){
1032 Lmunu = g[mu][mu]*kPrime[mu]*g[nu][nu]*k[nu]+g[nu][nu]*kPrime[nu]*g[mu][mu]*k[mu]-g[mu][nu]*kPrimek;
1037 leviCivitaIndexArray[0] = mu;
1038 leviCivitaIndexArray[1] = nu;
1039 for(
int a=0;
a<4;
a++){
1040 for(
int b=0;
b<4;
b++){
1041 leviCivitaIndexArray[2] =
a;
1042 leviCivitaIndexArray[3] =
b;
1043 imaginaryPart += -
leviCivita(leviCivitaIndexArray)*kPrime[
a]*k[
b];
1048 Lmunu = g[mu][mu]*kPrime[mu]*g[nu][nu]*k[nu]+g[nu][nu]*kPrime[nu]*g[mu][mu]*k[mu] + iNum*imaginaryPart;
1049 Lnumu = g[nu][nu]*kPrime[nu]*g[mu][mu]*k[mu]+g[mu][mu]*kPrime[mu]*g[nu][nu]*k[nu ]- iNum*imaginaryPart;
1052 if(mu ==0 && nu == 0){
1053 Amunu = 16.0*F1V2*(2.0*rulin[0][0]*CN+2.0*q[0]*tulin[0]+q2/2.0)+
1055 (4.0-4.0*rulin[0][0]/M2-4.0*q[0]*tulin[0]/M2-q02*(4.0/q2+1.0/M2)) +
1056 4.0*FA2*(2.0*rulin[0][0]+2.0*q[0]*tulin[0]+(q2/2.0-2.0*M2))-
1057 (2.0*CL*Fp2*q2+8.0*FA*Fp*CL*M)*q02-16.0*F1V*xiF2V*(-q2+q02)*CN;
1060 }
else if(mu == 0 && nu == 3){
1061 Amunu = 16.0*F1V2*((2.0*rulin[0][3]+tulin[0]*dq)*CN+tulin[3]*q[0])+
1063 (-4.0*rulin[0][3]/M2-2.0*(dq*tulin[0]+q[0]*tulin[3])/M2-dq*q[0]*(4.0/q2+1.0/M2))+
1064 4.0*FA2*((2.0*rulin[0][3]+dq*tulin[0])*CL+q[0]*tulin[3])-
1065 (2.0*CL*Fp2*q2+8.0*FA*Fp*CL*M)*dq*q[0]-
1066 16.0*F1V*xiF2V*dq*q[0];
1069 sum += Lmunu*Anumu + Lnumu*Amunu;
1070 }
else if(mu == 3 && nu == 3){
1071 Amunu = 16.0*F1V2*(2.0*rulin[3][3]+2.0*dq*tulin[3]-q2/2.0)+
1072 2.0*q2*xiF2V2*(-4.0-4.0*rulin[3][3]/M2-4.0*dq*tulin[3]/M2-dq2*(4.0/q2+1.0/M2))+
1073 4.0*FA2*(2.0*rulin[3][3]+2.0*dq*tulin[3]-(q2/2.0-2.0*CL*M2))-
1074 (2.0*CL*Fp2*q2+8.0*FA*Fp*CL*M)*dq2-
1075 16.0*F1V*xiF2V*(q2+dq2);
1078 }
else if(mu ==1 && nu == 1){
1079 Amunu = 16.0*F1V2*(2.0*rulin[1][1]-q2/2.0)+
1080 2.0*q2*xiF2V2*(-4.0*CT-4.0*rulin[1][1]/M2) +
1081 4.0*FA2*(2.0*rulin[1][1]-(q2/2.0-2.0*CT*M2))-
1082 16.0*F1V*xiF2V*CT*q2;
1085 }
else if(mu == 2 && nu == 2){
1088 Amunu = 16.0*F1V2*(2.0*rulin[2][2]-q2/2.0)+
1089 2.0*q2*xiF2V2*(-4.0*CT-4.0*rulin[2][2]/M2) +
1090 4.0*FA2*(2.0*rulin[2][2]-(q2/2.0-2.0*CT*M2))-
1091 16.0*F1V*xiF2V*CT*q2;
1093 }
else if(mu ==1 && nu == 2){
1094 Amunu = sign*16.0*iNum*FA*(xiF2V+F1V)*(-dq*tulin[0]*CT + q[0]*tulin[3]);
1097 sum += Lmunu*Anumu+Lnumu*Amunu;
1108 double tmugev = leptonMom.E() - leptonMom.Mag();
1110 std::ofstream ffstream;
1113 ffstream << -q2 <<
"\t" << q[0] <<
"\t" << dq
1114 <<
"\t" << axx <<
"\t" << azz <<
"\t" << a0z
1115 <<
"\t" << a00 <<
"\t" << axy <<
"\t" 1116 << CT <<
"\t" << CL <<
"\t" << CN <<
"\t" 1117 << tmugev <<
"\t" << imU <<
"\t" 1118 << F1V <<
"\t" << xiF2V <<
"\t" 1119 << FA <<
"\t" << Fp <<
"\t" 1120 << tulin[0] <<
"\t"<< tulin[1] <<
"\t" 1121 << tulin[2] <<
"\t"<< tulin[3] <<
"\t" 1122 << rulin[0][0]<<
"\t"<< rulin[0][1]<<
"\t" 1123 << rulin[0][2]<<
"\t"<< rulin[0][3]<<
"\t" 1124 << rulin[1][0]<<
"\t"<< rulin[1][1]<<
"\t" 1125 << rulin[1][2]<<
"\t"<< rulin[1][3]<<
"\t" 1126 << rulin[2][0]<<
"\t"<< rulin[2][1]<<
"\t" 1127 << rulin[2][2]<<
"\t"<< rulin[2][3]<<
"\t" 1128 << rulin[3][0]<<
"\t"<< rulin[3][1]<<
"\t" 1129 << rulin[3][2]<<
"\t"<< rulin[3][3]<<
"\t" 1140 LOG(
"Nieves",
pWARN) <<
"Imaginary part of tensor contraction is nonzero " 1141 <<
"in QEL XSec, real(sum) = " << real(sum)
1142 <<
"imag(sum) = " << imag(sum);
1151 double Rcurr,
int A,
int Z):
1152 ROOT::Math::IBaseFunctionOneDim()
1173 return rhop*rin*rin/
fRcurr;
1179 ROOT::Math::IBaseFunctionOneDim *
1204 double rmaxfrac = 0.25;
1206 bool carbon =
false;
1209 fTensorsOutFile =
"gen.RPA";
1211 fTensorsOutFile =
"gen.noRPA";
1232 rmax = TMath::Max(rp,rn) + 9.25*TMath::Max(ap,an);
1235 double r = rmax * rmaxfrac;
1239 const InitialState & init_state = interaction -> InitState();
1249 double delta = (ein-0.025)/100.0;
1250 for(
int it=0;
it<100;
it++){
1252 double eout = ml + tmu;
1257 double q0 = ein-eout;
1258 double dq =
TMath::Sqrt(pin*pin+pout*pout-2.0*ctl*pin*pout);
1259 double q2 = q0*q0-dq*dq;
1266 TLorentzVector qTildeP4(0, 0, dq, q0);
1267 TLorentzVector inNucleonMomOnShell(0,0,0,0);
1271 TLorentzVector neutrinoMom(0,0,pout+dq,eout+q0);
1272 TLorentzVector leptonMom(0,0,pout,eout);
1276 double Vc = vcr(& target, r);
1280 int sign = is_neutrino ? 1 : -1;
1281 double El = leptonMom.E();
1282 double ElLocal = El - sign*Vc;
1283 if(ElLocal - ml <= 0.0){
1284 LOG(
"Nieves",
pINFO) <<
"Event should be rejected. Coulomb effects " 1285 <<
"push kinematics below threshold";
1288 double plLocal =
TMath::Sqrt(ElLocal*ElLocal-ml*ml);
1291 double coulombFactor= plLocal*ElLocal/leptonMom.Vect().Mag()/El;
1292 fCoulombFactor = coulombFactor;
1297 fFormFactors.Calculate(interaction);
1298 LmunuAnumu(neutrinoMom, inNucleonMomOnShell, leptonMom, qTildeP4,
1299 M, is_neutrino, target,
false);
Cross Section Calculation Interface.
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
std::complex< double > ruLinRelX(double q0, double qm, double kf, double m) const
double DoEval(double rin) const
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
bool IsWeakCC(void) const
bool fRPA
use RPA corrections
bool IsNeutrino(int pdgc)
void Configure(const Registry &config)
double J(double q0, double q3, double Enu, double ml)
THE MAIN GENIE PROJECT NAMESPACE
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
int HitNucPdg(void) const
void SetQ2(double Q2, bool selected=false)
Kinematics * KinePtr(void) const
double Density(double r, int A, double ring=0.)
TParticlePDG * RecoilNucleon(void) const
recoil nucleon
QELEvGen_BindingMode_t fIntegralNucleonBindingMode
const FermiMomentumTable * fKFTable
int RecoilNucleonPdg(void) const
recoil nucleon pdg
static const double fermi
bool IsQuasiElastic(void) const
double HitNucMass(void) const
bool IsNucleus(void) const
void CNCTCLimUcalc(TLorentzVector qTildeP4, double M, double r, bool is_neutrino, bool tgtIsNucleus, int tgt_pdgc, int A, int Z, int N, bool hitNucIsProton, double &CN, double &CT, double &CL, double &imU, double &t0, double &r00, bool assumeFreeNucleon) const
static FermiMomentumTablePool * Instance(void)
Generated/set kinematical variables for an event.
Examines whether the generated event should be Pauli blocked. Is a concerete implementation of the Ev...
const XSecIntegratorI * fXSecIntegrator
const TLorentzVector & HadSystP4(void) const
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
QELFormFactors fFormFactors
bool fCoulomb
use Coulomb corrections
enum genie::EKinePhaseSpace KinePhaseSpace_t
Nieves_Coulomb_Rmax_t fCoulombRmaxMode
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
std::complex< double > deltaLindhard(double q0gev, double dqgev, double rho, double kFgev) const
virtual NuclearModel_t ModelType(const Target &) const =0
const QELFormFactorsModelI * fFormFactorsModel
double fhbarc
hbar*c in GeV*fm
double fCos8c2
cos^2(cabibbo angle)
Summary information for an interaction.
const TLorentzVector & HitNucP4(void) const
double GetFermiMomentum(const Target &tgt, int pdg_Nf, double radius=0.0) const
Get the Fermi momentum needed to check Pauli blocking.
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
const NuclearModelI * fNuclModel
Nuclear Model for integration.
const TLorentzVector & FSLeptonP4(void) const
Singleton class to load & serve tables of Fermi momentum constants.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
NievesQELvcrIntegrand(double Rcurr, int A, int Z)
const FermiMomentumTable * GetTable(string name)
double EnergyDeltaFunctionSolutionQEL(const Interaction &inter)
static string AsString(KinePhaseSpace_t kps)
bool fCompareNievesTensors
print tensors
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
bool IsAntiNeutrino(int pdgc)
double Integral(const Interaction *i) const
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
virtual void Configure(const Registry &config)
bool fDoPauliBlocking
Whether to apply Pauli blocking in XSec()
double func(double x, double y)
static const double kASmallNum
unsigned int NDim(void) const
double fXSecScale
external xsec scaling factor
Misc GENIE control constants.
A very simple service to remember what detector we're working in.
static const double kLightSpeed
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
double vcr(const Target *target, double r) const
static constexpr Double_t m2
TString fTensorsOutFile
file to print tensors to
virtual const AlgId & Id(void) const
Get algorithm ID.
virtual ~NievesQELCCPXSec()
Var Sqrt(const Var &v)
Use to take sqrt of a var.
A registry. Provides the container for algorithm configuration parameters.
const UInt_t kIAssumeFreeNucleon
double HitNucPosition(void) const
QELEvGen_BindingMode_t StringToQELBindingMode(const std::string &mode_str)
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
assert(nhit_max >=nhit_nbins)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
const InitialState & InitState(void) const
const ProcessInfo & ProcInfo(void) const
int leviCivita(int input[]) const
void CompareNievesTensors(const Interaction *i) const
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
std::complex< double > relLindhard(double q0gev, double dqgev, double kFgev, double M, bool isNeutrino, std::complex< double > relLindIm) const
double LmunuAnumu(const TLorentzVector neutrinoMom, const TLorentzVector inNucleonMom, const TLorentzVector leptonMom, const TLorentzVector outNucleonMom, double M, bool is_neutrino, const Target &target, bool assumeFreeNucleon) const
const genie::PauliBlocker * fPauliBlocker
The PauliBlocker instance to use to apply that correction.
static const double kPlankConstant
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Root of GENIE utility namespaces.
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
std::complex< double > relLindhardIm(double q0gev, double dqgev, double kFngev, double kFpgev, double M, bool isNeutrino, double &t0, double &r00) const
Initial State information.
static const double kPionMass2
const Algorithm * SubAlg(const RgKey ®istry_key) const