49 using namespace genie;
74 if(!
this ->
ValidProcess (interaction) ) {
LOG(
"LwlynSmith",
pWARN) <<
"not a valid process";
return 0.;}
86 const InitialState & init_state = interaction -> InitState();
90 double E2 = TMath::Power(E,2);
93 double q2 = kinematics.
q2();
97 int sign = (is_neutrino) ? -1 : 1;
107 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 112 double ml2 = TMath::Power(ml, 2);
113 double M2 = TMath::Power(M, 2);
114 double M4 = TMath::Power(M2, 2);
115 double FA2 = TMath::Power(FA, 2);
116 double Fp2 = TMath::Power(Fp, 2);
117 double F1V2 = TMath::Power(F1V, 2);
118 double xiF2V2 = TMath::Power(xiF2V, 2);
120 double s_u = 4*E*M + q2 - ml2;
121 double q2_M2 = q2/M2;
124 double A = (0.25*(ml2-
q2)/M2) * (
125 (4-q2_M2)*FA2 - (4+q2_M2)*F1V2 - q2_M2*xiF2V2*(1+0.25*q2_M2)
126 -4*q2_M2*F1V*xiF2V - (ml2/M2)*(
127 (F1V2+xiF2V2+2*F1V*xiF2V)+(FA2+4*Fp2+4*FA*Fp)+(q2_M2-4)*Fp2));
128 double B = -1 * q2_M2 * FA*(F1V+xiF2V);
129 double C = 0.25*(FA2 + F1V2 - 0.25*q2_M2*xiF2V2);
131 double xsec = Gfactor * (A + sign*B*s_u/M2 + C*s_u*s_u/M4);
136 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 138 <<
"dXSec[QEL]/dQ2 [FreeN](E = "<< E <<
", Q2 = "<< -q2 <<
") = "<<
xsec;
140 <<
"A(Q2) = " << A <<
", B(Q2) = " << B <<
", C(Q2) = " <<
C;
148 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 150 <<
"Jacobian for transformation to: " 167 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 169 <<
"Nuclear suppression factor R(Q2) = " << R <<
", NNucl = " << NNucl;
183 const InitialState& init_state = interaction -> InitState();
187 const TLorentzVector leptonMom = kinematics.
FSLeptonP4();
188 const TLorentzVector outNucleonMom = kinematics.
HadSystP4();
195 double pNf = outNucleonMom.P();
196 if ( pNf < kF )
return 0.;
203 TLorentzVector neutrinoMom = *tempNeutrino;
223 double mNucleon = ( mNi + interaction->
RecoilNucleon()->Mass() ) / 2.;
226 TLorentzVector qP4 = neutrinoMom - leptonMom;
231 TLorentzVector inNucleonMomOnShell( inNucleonMom->Vect(), inNucleonOnShellEnergy );
235 TLorentzVector qTildeP4 = outNucleonMom - inNucleonMomOnShell;
237 double Q2 = -1. * qP4.Mag2();
238 double Q2tilde = -1. * qTildeP4.Mag2();
242 if ( qTildeP4.E() <= 0. && init_state.
Tgt().
IsNucleus() &&
244 if ( Q2tilde <= 0. )
return 0.;
265 * neutrinoMom.E() * outNucleonMom.E() * leptonMom.E() );
268 double tau = Q2tilde / (4 *
std::pow(mNucleon, 2));
269 double h1 = FA*FA*(1 + tau) + tau*(F1V + xiF2V)*(F1V + xiF2V);
270 double h2 = FA*FA + F1V*F1V + tau*xiF2V*xiF2V;
271 double h3 = 2.0 * FA * (F1V + xiF2V);
272 double h4 = 0.25 * xiF2V*xiF2V *(1-tau) + 0.5*F1V*xiF2V + FA*Fp - tau*Fp*Fp;
275 int sign = (is_neutrino) ? -1 : 1;
276 double l1 = 2*neutrinoMom.Dot(leptonMom)*
std::pow(mNucleon, 2);
277 double l2 = 2*(neutrinoMom.Dot(inNucleonMomOnShell)) * (inNucleonMomOnShell.Dot(leptonMom)) - neutrinoMom.Dot(leptonMom)*
std::pow(mNucleon, 2);
278 double l3 = (neutrinoMom.Dot(inNucleonMomOnShell) * qTildeP4.Dot(leptonMom)) - (neutrinoMom.Dot(qTildeP4) * leptonMom.Dot(inNucleonMomOnShell));
280 double l4 = neutrinoMom.Dot(leptonMom) * qTildeP4.Dot(qTildeP4) - 2*neutrinoMom.Dot(qTildeP4)*leptonMom.Dot(qTildeP4);
281 double l5 = neutrinoMom.Dot(inNucleonMomOnShell) * leptonMom.Dot(qTildeP4) + leptonMom.Dot(inNucleonMomOnShell)*neutrinoMom.Dot(qTildeP4) - neutrinoMom.Dot(leptonMom)*inNucleonMomOnShell.Dot(qTildeP4);
283 double LH = 2 *(l1*h1 + l2*h2 + l3*h3 + l4*h4 + l5*
h2);
285 double xsec = Gfactor * LH;
299 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 301 <<
"Nuclear suppression factor R(Q2) = " <<
R <<
", NNucl = " << NNucl;
343 int Zf = (is_p) ? Zi-1 : Zi;
350 <<
"Unknwown nuclear target! No target with code: " 354 double Mi = nucl_i ->
Mass();
355 double Mf = nucl_f ->
Mass();
360 double xsec_sum = 0.;
361 const int nnuc = 2000;
366 for(
int inuc=0;inuc<nnuc;inuc++){
376 p4N->SetPx (p3N.Px());
377 p4N->SetPy (p3N.Py());
378 p4N->SetPz (p3N.Pz());
384 double xsec_avg = xsec_sum / nnuc;
409 bool prcok = proc_info.
IsWeakCC() && ((isP&&isnub) || (isN&&isnu));
410 if(!prcok)
return false;
434 fCos8c2 = TMath::Power(TMath::Cos(thc), 2);
438 this->
SubAlg(
"FormFactorsAlg"));
448 RgKey nuclkey =
"IntegralNuclearModel";
454 bool average_over_nuc_mom ;
455 GetParamDef(
"IntegralAverageOverNucleonMomentum", average_over_nuc_mom,
false ) ;
double FullDifferentialXSec(const Interaction *i) const
Cross Section Calculation Interface.
TVector3 GenerateVertex(const Interaction *in, double A) const
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
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
TParticlePDG * RecoilNucleon(void) const
recoil nucleon
bool fLFG
If the nuclear model is lfg alway average over nucleons.
int RecoilNucleonPdg(void) const
recoil nucleon pdg
bool IsQuasiElastic(void) const
double HitNucMass(void) const
bool IsNucleus(void) const
Generated/set kinematical variables for an event.
Examines whether the generated event should be Pauli blocked. Is a concerete implementation of the Ev...
const NuclearModelI * fNuclModel
const TLorentzVector & HadSystP4(void) const
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
double Mass(Resonance_t res)
resonance mass (GeV)
void SetHitNucPosition(double r)
double fXSecScale
external xsec scaling factor
enum genie::EKinePhaseSpace KinePhaseSpace_t
bool fDoAvgOverNucleonMomentum
Average cross section over hit nucleon monentum?
double Integral(const Interaction *i) const
const TVector3 & Momentum3(void) const
QELEvGen_BindingMode_t fIntegralNucleonBindingMode
virtual NuclearModel_t ModelType(const Target &) const =0
Summary information for an interaction.
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?
double q2(bool selected=false) const
const QELFormFactorsModelI * fFormFactorsModel
const XSecIntegratorI * fXSecIntegrator
const TLorentzVector & FSLeptonP4(void) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
virtual ~LwlynSmithQELCCPXSec()
double EnergyDeltaFunctionSolutionQEL(const Interaction &inter)
static string AsString(KinePhaseSpace_t kps)
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
bool IsAntiNeutrino(int pdgc)
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 FullDifferentialXSec.
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
QELFormFactors fFormFactors
void Configure(const Registry &config)
virtual const AlgId & Id(void) const
Get algorithm ID.
TLorentzVector * HitNucP4Ptr(void) const
static PDGLibrary * Instance(void)
Singleton class to load & serve a TDatabasePDG.
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
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
double HitNucPosition(void) const
QELEvGen_BindingMode_t StringToQELBindingMode(const std::string &mode_str)
Target * TgtPtr(void) const
double fCos8c2
cos^2(cabibbo angle)
TParticlePDG * Find(int pdgc)
int IonPdgCode(int A, int Z)
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)
virtual bool GenerateNucleon(const Target &) const =0
const InitialState & InitState(void) const
const ProcessInfo & ProcInfo(void) 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
double ProbeE(RefFrame_t rf) const
const genie::PauliBlocker * fPauliBlocker
The PauliBlocker instance to use to apply that correction.
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Root of GENIE utility namespaces.
void Configure(const Registry &config)
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
double NuclQELXSecSuppression(string kftable, double pmax, const Interaction *in)
Initial State information.
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
const Algorithm * SubAlg(const RgKey ®istry_key) const