LwlynSmithQELCCPXSec.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Sep 19, 2009 - CA
14  Renamed QELPXSec -> LwlynSmithQELCCPXSec
15  @ Mar 18, 2016 - JJ (SD)
16  Moved code to average over initial nucleons from QELXSec to the Integral()
17  method here. For each nucleon, generate a struck nucleon position, then a
18  momentum, then integrate.
19  @ 2015 - AF
20  Added FullDifferentialXSec method to work with QELEventGenerator
21 */
22 //____________________________________________________________________________
23 
24 #include <TMath.h>
25 
38 
48 
49 using namespace genie;
50 using namespace genie::constants;
51 using namespace genie::utils;
52 
53 //____________________________________________________________________________
55 XSecAlgorithmI("genie::LwlynSmithQELCCPXSec")
56 {
57 
58 }
59 //____________________________________________________________________________
61 XSecAlgorithmI("genie::LwlynSmithQELCCPXSec", config)
62 {
63 
64 }
65 //____________________________________________________________________________
67 {
68 
69 }
70 //____________________________________________________________________________
72  const Interaction * interaction, KinePhaseSpace_t kps) const
73 {
74  if(! this -> ValidProcess (interaction) ) {LOG("LwlynSmith",pWARN) << "not a valid process"; return 0.;}
75  if(! this -> ValidKinematics (interaction) ) {LOG("LwlynSmith",pWARN) << "not valid kinematics"; return 0.;}
76 
77  // If computing the full differential cross section, then all four momentum
78  // four-vectors (probe, hit nucleon, final lepton, and final nucleon) should
79  // have been set already, with the hit nucleon off-shell as appropriate.
80  if (kps == kPSQELEvGen) {
81  return this->FullDifferentialXSec(interaction);
82  }
83 
84  // Get kinematics & init-state parameters
85  const Kinematics & kinematics = interaction -> Kine();
86  const InitialState & init_state = interaction -> InitState();
87  const Target & target = init_state.Tgt();
88 
89  double E = init_state.ProbeE(kRfHitNucRest);
90  double E2 = TMath::Power(E,2);
91  double ml = interaction->FSPrimLepton()->Mass();
92  double M = target.HitNucMass();
93  double q2 = kinematics.q2();
94 
95  // One of the xsec terms changes sign for antineutrinos
96  bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
97  int sign = (is_neutrino) ? -1 : 1;
98 
99  // Calculate the QEL form factors
100  fFormFactors.Calculate(interaction);
101 
102  double F1V = fFormFactors.F1V();
103  double xiF2V = fFormFactors.xiF2V();
104  double FA = fFormFactors.FA();
105  double Fp = fFormFactors.Fp();
106 
107 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
108  LOG("LwlynSmith", pDEBUG) << "\n" << fFormFactors;
109 #endif
110 
111  // Calculate auxiliary parameters
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);
119  double Gfactor = M2*kGF2*fCos8c2 / (8*kPi*E2);
120  double s_u = 4*E*M + q2 - ml2;
121  double q2_M2 = q2/M2;
122 
123  // Compute free nucleon differential cross section
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);
130 
131  double xsec = Gfactor * (A + sign*B*s_u/M2 + C*s_u*s_u/M4);
132 
133  // Apply given scaling factor
134  xsec *= fXSecScale;
135 
136 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
137  LOG("LwlynSmith", pDEBUG)
138  << "dXSec[QEL]/dQ2 [FreeN](E = "<< E << ", Q2 = "<< -q2 << ") = "<< xsec;
139  LOG("LwlynSmith", pDEBUG)
140  << "A(Q2) = " << A << ", B(Q2) = " << B << ", C(Q2) = " << C;
141 #endif
142 
143  //----- The algorithm computes dxsec/dQ2
144  // Check whether variable tranformation is needed
145  if(kps!=kPSQ2fE) {
146  double J = utils::kinematics::Jacobian(interaction,kPSQ2fE,kps);
147 
148 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
149  LOG("LwlynSmith", pDEBUG)
150  << "Jacobian for transformation to: "
151  << KinePhaseSpace::AsString(kps) << ", J = " << J;
152 #endif
153  xsec *= J;
154  }
155 
156  //----- if requested return the free nucleon xsec even for input nuclear tgt
157  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
158 
159  //----- compute nuclear suppression factor
160  // (R(Q2) is adapted from NeuGEN - see comments therein)
161  double R = nuclear::NuclQELXSecSuppression("Default", 0.5, interaction);
162 
163  //----- number of scattering centers in the target
164  int nucpdgc = target.HitNucPdg();
165  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
166 
167 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
168  LOG("LwlynSmith", pDEBUG)
169  << "Nuclear suppression factor R(Q2) = " << R << ", NNucl = " << NNucl;
170 #endif
171 
172  xsec *= (R*NNucl); // nuclear xsec
173 
174  return xsec;
175 }
176 //____________________________________________________________________________
178  const
179 {
180  // First we need access to all of the particles in the interaction
181  // The particles were stored in the lab frame
182  const Kinematics& kinematics = interaction -> Kine();
183  const InitialState& init_state = interaction -> InitState();
184 
185  const Target& tgt = init_state.Tgt();
186 
187  const TLorentzVector leptonMom = kinematics.FSLeptonP4();
188  const TLorentzVector outNucleonMom = kinematics.HadSystP4();
189 
190  // Apply Pauli blocking if enabled
191  if ( fDoPauliBlocking && tgt.IsNucleus() && !interaction->TestBit(kIAssumeFreeNucleon) ) {
192  int final_nucleon_pdg = interaction->RecoilNucleonPdg();
193  double kF = fPauliBlocker->GetFermiMomentum(tgt, final_nucleon_pdg,
194  tgt.HitNucPosition());
195  double pNf = outNucleonMom.P();
196  if ( pNf < kF ) return 0.;
197  }
198 
199  // Note that GetProbeP4 defaults to returning the probe 4-momentum in the
200  // struck nucleon rest frame, so we have to explicitly ask for the lab frame
201  // here
202  TLorentzVector * tempNeutrino = init_state.GetProbeP4(kRfLab);
203  TLorentzVector neutrinoMom = *tempNeutrino;
204  delete tempNeutrino;
205  TLorentzVector * inNucleonMom = init_state.TgtPtr()->HitNucP4Ptr();
206 
207  // *** CALCULATION OF "q" and "qTilde" ***
208  // According to the de Forest prescription for handling the off-shell
209  // initial struck nucleon, the cross section calculation should proceed
210  // as if for a free nucleon, except that an effective value of the 4-momentum
211  // transfer qTilde should be used in which the difference between the on-
212  // and off-shell energies of the hit nucleon has been subtracted from the
213  // energy transfer q0.
214 
215  // HitNucMass() looks up the PDGLibrary (on-shell) value for the initial
216  // struck nucleon
217  double mNi = init_state.Tgt().HitNucMass();
218 
219  // Hadronic matrix element for CC neutrino interactions should really use
220  // the "nucleon mass," i.e., the mean of the proton and neutrino masses.
221  // This expression would also work for NC and EM scattering (since the
222  // initial and final on-shell nucleon masses would be the same)
223  double mNucleon = ( mNi + interaction->RecoilNucleon()->Mass() ) / 2.;
224 
225  // Ordinary 4-momentum transfer
226  TLorentzVector qP4 = neutrinoMom - leptonMom;
227 
228  // Initial struck nucleon 4-momentum in which it is forced to be on-shell
229  double inNucleonOnShellEnergy = std::sqrt( std::pow(mNi, 2)
230  + std::pow(inNucleonMom->P(), 2) );
231  TLorentzVector inNucleonMomOnShell( inNucleonMom->Vect(), inNucleonOnShellEnergy );
232 
233  // Effective 4-momentum transfer (according to the deForest prescription) for
234  // use in computing the hadronic tensor
235  TLorentzVector qTildeP4 = outNucleonMom - inNucleonMomOnShell;
236 
237  double Q2 = -1. * qP4.Mag2();
238  double Q2tilde = -1. * qTildeP4.Mag2();
239 
240  // If the binding energy correction causes an unphysical value
241  // of q0Tilde or Q2tilde, just return 0.
242  if ( qTildeP4.E() <= 0. && init_state.Tgt().IsNucleus() &&
243  !interaction->TestBit(kIAssumeFreeNucleon) ) return 0.;
244  if ( Q2tilde <= 0. ) return 0.;
245 
246  // Store Q2tilde in the kinematic variable representing Q2.
247  // This will ensure that the form factors are calculated correctly
248  // using the de Forest prescription (Q2tilde instead of Q2).
249  interaction->KinePtr()->SetQ2(Q2tilde);
250 
251  // Calculate the QEL form factors
252  fFormFactors.Calculate(interaction);
253 
254  double F1V = fFormFactors.F1V();
255  double xiF2V = fFormFactors.xiF2V();
256  double FA = fFormFactors.FA();
257  double Fp = fFormFactors.Fp();
258 
259  // Restore Q2 in the interaction's kinematic variables
260  // now that the form factors have been computed
261  interaction->KinePtr()->SetQ2( Q2 );
262 
263  // Overall factor in the differential cross section
264  double Gfactor = kGF2*fCos8c2 / ( 8. * kPi * kPi * inNucleonOnShellEnergy
265  * neutrinoMom.E() * outNucleonMom.E() * leptonMom.E() );
266 
267  // Now, we can calculate the cross section
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;
273 
274  bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
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));
279  l3 *= sign;
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);
282 
283  double LH = 2 *(l1*h1 + l2*h2 + l3*h3 + l4*h4 + l5*h2);
284 
285  double xsec = Gfactor * LH;
286 
287  // Apply the factor that arises from elimination of the energy-conserving
288  // delta function
289  xsec *= genie::utils::EnergyDeltaFunctionSolutionQEL( *interaction );
290 
291  // Apply given scaling factor
292  xsec *= fXSecScale;
293 
294  // Number of scattering centers in the target
295  const Target & target = init_state.Tgt();
296  int nucpdgc = target.HitNucPdg();
297  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
298 
299 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
300  LOG("LwlynSmith", pDEBUG)
301  << "Nuclear suppression factor R(Q2) = " << R << ", NNucl = " << NNucl;
302 #endif
303 
304  xsec *= NNucl; // nuclear xsec
305 
306  return xsec;
307 
308 }
309 //____________________________________________________________________________
311 {
312  // If we're using the new spline generation method (which integrates
313  // over the kPSQELEvGen phase space used by QELEventGenerator) then
314  // let the cross section integrator do all of the work. It's smart
315  // enough to handle free nucleon vs. nuclear targets, different
316  // nuclear models (including the local Fermi gas model), etc.
317  // TODO: think about doing this in a better way
318  if ( fXSecIntegrator->Id().Name() == "genie::NewQELXSec" ) {
319  return fXSecIntegrator->Integrate(this, in);
320  }
321 
322  // Otherwise, use the old integration method (kept for use with
323  // the historical default G18_00x series of tunes)
324  bool nuclear_target = in->InitState().Tgt().IsNucleus();
325  if(!nuclear_target || !fDoAvgOverNucleonMomentum) {
326  return fXSecIntegrator->Integrate(this,in);
327  }
328 
329  double E = in->InitState().ProbeE(kRfHitNucRest);
330  if(fLFG || E < fEnergyCutOff) {
331  // clone the input interaction so as to tweak the
332  // hit nucleon 4-momentum in the averaging loop
333  Interaction in_curr(*in);
334 
335  // hit target
336  Target * tgt = in_curr.InitState().TgtPtr();
337 
338  // get nuclear masses (init & final state nucleus)
339  int nucleon_pdgc = tgt->HitNucPdg();
340  bool is_p = pdg::IsProton(nucleon_pdgc);
341  int Zi = tgt->Z();
342  int Ai = tgt->A();
343  int Zf = (is_p) ? Zi-1 : Zi;
344  int Af = Ai-1;
345  PDGLibrary * pdglib = PDGLibrary::Instance();
346  TParticlePDG * nucl_i = pdglib->Find( pdg::IonPdgCode(Ai, Zi) );
347  TParticlePDG * nucl_f = pdglib->Find( pdg::IonPdgCode(Af, Zf) );
348  if(!nucl_f) {
349  LOG("LwlynSmith", pFATAL)
350  << "Unknwown nuclear target! No target with code: "
351  << pdg::IonPdgCode(Af, Zf) << " in PDGLibrary!";
352  exit(1);
353  }
354  double Mi = nucl_i -> Mass(); // initial nucleus mass
355  double Mf = nucl_f -> Mass(); // remnant nucleus mass
356 
357  // throw nucleons with fermi momenta and binding energies
358  // generated according to the current nuclear model for the
359  // input target and average the cross section
360  double xsec_sum = 0.;
361  const int nnuc = 2000;
362  // VertexGenerator for generating a position before generating
363  // each nucleon
364  VertexGenerator * vg = new VertexGenerator();
365  vg->Configure("Default");
366  for(int inuc=0;inuc<nnuc;inuc++){
367  // Generate a position in the nucleus
368  TVector3 nucpos = vg->GenerateVertex(&in_curr,tgt->A());
369  tgt->SetHitNucPosition(nucpos.Mag());
370 
371  // Generate a nucleon
372  fNuclModel->GenerateNucleon(*tgt, nucpos.Mag());
373  TVector3 p3N = fNuclModel->Momentum3();
374  double EN = Mi - TMath::Sqrt(p3N.Mag2() + Mf*Mf);
375  TLorentzVector* p4N = tgt->HitNucP4Ptr();
376  p4N->SetPx (p3N.Px());
377  p4N->SetPy (p3N.Py());
378  p4N->SetPz (p3N.Pz());
379  p4N->SetE (EN);
380 
381  double xsec = fXSecIntegrator->Integrate(this,&in_curr);
382  xsec_sum += xsec;
383  }
384  double xsec_avg = xsec_sum / nnuc;
385  delete vg;
386  return xsec_avg;
387  }else{
388  return fXSecIntegrator->Integrate(this,in);
389  }
390 }
391 //____________________________________________________________________________
392 bool LwlynSmithQELCCPXSec::ValidProcess(const Interaction * interaction) const
393 {
394  if(interaction->TestBit(kISkipProcessChk)) return true;
395 
396  const InitialState & init_state = interaction->InitState();
397  const ProcessInfo & proc_info = interaction->ProcInfo();
398 
399  if(!proc_info.IsQuasiElastic()) return false;
400 
401  int nuc = init_state.Tgt().HitNucPdg();
402  int nu = init_state.ProbePdg();
403 
404  bool isP = pdg::IsProton(nuc);
405  bool isN = pdg::IsNeutron(nuc);
406  bool isnu = pdg::IsNeutrino(nu);
407  bool isnub = pdg::IsAntiNeutrino(nu);
408 
409  bool prcok = proc_info.IsWeakCC() && ((isP&&isnub) || (isN&&isnu));
410  if(!prcok) return false;
411 
412  return true;
413 }
414 //____________________________________________________________________________
416 {
417  Algorithm::Configure(config);
418  this->LoadConfig();
419 }
420 //____________________________________________________________________________
422 {
423  Algorithm::Configure(config);
424  this->LoadConfig();
425 }
426 //____________________________________________________________________________
428 {
429  // Cross section scaling factor
430  GetParamDef( "QEL-CC-XSecScale", fXSecScale, 1. ) ;
431 
432  double thc ;
433  GetParam( "CabibboAngle", thc ) ;
434  fCos8c2 = TMath::Power(TMath::Cos(thc), 2);
435 
436  // load QEL form factors model
437  fFormFactorsModel = dynamic_cast<const QELFormFactorsModelI *> (
438  this->SubAlg("FormFactorsAlg"));
440  fFormFactors.SetModel(fFormFactorsModel); // <-- attach algorithm
441 
442  // load XSec Integrator
444  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
446 
447  // Get nuclear model for use in Integral()
448  RgKey nuclkey = "IntegralNuclearModel";
449  fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
451 
453 
454  bool average_over_nuc_mom ;
455  GetParamDef( "IntegralAverageOverNucleonMomentum", average_over_nuc_mom, false ) ;
456  // Always average over initial nucleons if the nuclear model is LFG
457  fDoAvgOverNucleonMomentum = fLFG || average_over_nuc_mom ;
458 
459  fEnergyCutOff = 0.;
460 
461  // Get averaging cutoff energy
462  GetParamDef("IntegralNuclearInfluenceCutoffEnergy", fEnergyCutOff, 2.0 ) ;
463 
464  // Method to use to calculate the binding energy of the initial hit nucleon when
465  // generating splines
466  std::string temp_binding_mode;
467  GetParamDef( "IntegralNucleonBindingMode", temp_binding_mode, std::string("UseNuclearModel") );
469 
470  // Get PauliBlocker for possible use in FullDifferentialXSec()
471  GetParamDef( "IntegralNucleonBindingMode", temp_binding_mode, std::string("UseNuclearModel") );
472  fPauliBlocker = dynamic_cast<const PauliBlocker*>( this->SubAlg("PauliBlockerAlg") );
474 
475  // Decide whether or not it should be used in FullDifferentialXSec
476  GetParamDef( "DoPauliBlocking", fDoPauliBlocking, true );
477 }
double FullDifferentialXSec(const Interaction *i) const
Cross Section Calculation Interface.
const double kPi
TVector3 GenerateVertex(const Interaction *in, double A) const
Basic constants.
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
const XML_Char * target
Definition: expat.h:268
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
int HitNucPdg(void) const
Definition: Target.cxx:321
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:265
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
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
int A(void) const
Definition: Target.h:71
void SetModel(const QELFormFactorsModelI *model)
Attach an algorithm.
TH1F * h3
Definition: berger.C:36
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:67
double HitNucMass(void) const
Definition: Target.cxx:250
T sqrt(T number)
Definition: d0nt_math.hpp:156
#define pFATAL
Definition: Messenger.h:57
bool IsNucleus(void) const
Definition: Target.cxx:289
constexpr T pow(T x)
Definition: pow.h:75
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
Examines whether the generated event should be Pauli blocked. Is a concerete implementation of the Ev...
Definition: PauliBlocker.h:36
const NuclearModelI * fNuclModel
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:67
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
double Mass(Resonance_t res)
resonance mass (GeV)
void SetHitNucPosition(double r)
Definition: Target.cxx:227
double fXSecScale
external xsec scaling factor
Definition: config.py:1
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
Definition: NuclearModelI.h:67
QELEvGen_BindingMode_t fIntegralNucleonBindingMode
virtual NuclearModel_t ModelType(const Target &) const =0
Double_t q2[12][num]
Definition: f2_nu.C:137
const double C
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
Summary information for an interaction.
Definition: Interaction.h:56
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
Definition: Kinematics.cxx:151
const QELFormFactorsModelI * fFormFactorsModel
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
const XSecIntegratorI * fXSecIntegrator
const TLorentzVector & FSLeptonP4(void) const
Definition: Kinematics.h:66
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Pure abstract base class. Defines the QELFormFactorsModelI interface to be implemented by any algorit...
string Name(void) const
Definition: AlgId.h:45
Float_t E
Definition: plot.C:20
double EnergyDeltaFunctionSolutionQEL(const Interaction &inter)
Definition: QELUtils.cxx:51
static string AsString(KinePhaseSpace_t kps)
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
int ProbePdg(void) const
Definition: InitialState.h:65
#define R(x)
int Z(void) const
Definition: Target.h:69
bool fDoPauliBlocking
Whether to apply Pauli blocking in FullDifferentialXSec.
double xiF2V(void) const
Get the computed form factor xi*F2V.
Double_t xsec[nknots]
Definition: testXsec.C:47
#define pWARN
Definition: Messenger.h:61
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
TH1F * h2
Definition: plot.C:45
void Calculate(const Interaction *interaction)
Compute the form factors for the input interaction using the attached model.
void Configure(const Registry &config)
TH1F * h1
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
static const double A
Definition: Units.h:82
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
ifstream in
Definition: comparison.C:7
int N(void) const
Definition: Target.h:70
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
double HitNucPosition(void) const
Definition: Target.h:90
QELEvGen_BindingMode_t StringToQELBindingMode(const std::string &mode_str)
Definition: QELUtils.cxx:195
Target * TgtPtr(void) const
Definition: InitialState.h:68
double fCos8c2
cos^2(cabibbo angle)
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
exit(0)
assert(nhit_max >=nhit_nbins)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:128
virtual bool GenerateNucleon(const Target &) const =0
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double F1V(void) const
Get the computed form factor F1V.
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
Definition: InitialState.h:67
static const double kGF2
Definition: Constants.h:60
double Fp(void) const
Get the computed form factor Fp.
void kinematics()
Definition: kinematics.C:10
double ProbeE(RefFrame_t rf) const
const genie::PauliBlocker * fPauliBlocker
The PauliBlocker instance to use to apply that correction.
double FA(void) const
Get the computed form factor FA.
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
def sign(x)
Definition: canMan.py:197
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
double NuclQELXSecSuppression(string kftable, double pmax, const Interaction *in)
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:353