Public Member Functions | Private Member Functions | Private Attributes | List of all members
genie::LwlynSmithQELCCPXSec Class Reference

Computes neutrino-nucleon(nucleus) QELCC differential cross section Is a concrete implementation of the XSecAlgorithmI interface.
. More...

#include <LwlynSmithQELCCPXSec.h>

Inheritance diagram for genie::LwlynSmithQELCCPXSec:
genie::XSecAlgorithmI genie::Algorithm

Public Member Functions

 LwlynSmithQELCCPXSec ()
 
 LwlynSmithQELCCPXSec (string config)
 
virtual ~LwlynSmithQELCCPXSec ()
 
double XSec (const Interaction *i, KinePhaseSpace_t k) const
 Compute the cross section for the input interaction. More...
 
double Integral (const Interaction *i) const
 
bool ValidProcess (const Interaction *i) const
 Can this cross section algorithm handle the input process? More...
 
void Configure (const Registry &config)
 
void Configure (string param_set)
 
- Public Member Functions inherited from genie::XSecAlgorithmI
virtual ~XSecAlgorithmI ()
 
virtual bool ValidKinematics (const Interaction *i) const
 Is the input kinematical point a physically allowed one? More...
 
- Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
 
virtual void FindConfig (void)
 
virtual const RegistryGetConfig (void) const
 
RegistryGetOwnedConfig (void)
 
virtual const AlgIdId (void) const
 Get algorithm ID. More...
 
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status. More...
 
virtual bool AllowReconfig (void) const
 
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm. More...
 
virtual void SetId (const AlgId &id)
 Set algorithm ID. More...
 
virtual void SetId (string name, string config)
 
const AlgorithmSubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

Private Member Functions

double FullDifferentialXSec (const Interaction *i) const
 
void LoadConfig (void)
 

Private Attributes

QELFormFactors fFormFactors
 
const QELFormFactorsModelIfFormFactorsModel
 
const XSecIntegratorIfXSecIntegrator
 
double fCos8c2
 cos^2(cabibbo angle) More...
 
double fXSecScale
 external xsec scaling factor More...
 
const NuclearModelIfNuclModel
 
bool fLFG
 If the nuclear model is lfg alway average over nucleons. More...
 
bool fDoAvgOverNucleonMomentum
 Average cross section over hit nucleon monentum? More...
 
double fEnergyCutOff
 
QELEvGen_BindingMode_t fIntegralNucleonBindingMode
 
bool fDoPauliBlocking
 Whether to apply Pauli blocking in FullDifferentialXSec. More...
 
const genie::PauliBlockerfPauliBlocker
 The PauliBlocker instance to use to apply that correction. More...
 

Additional Inherited Members

- Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
 
static string BuildParamVectSizeKey (const std::string &comm_name)
 
- Protected Member Functions inherited from genie::XSecAlgorithmI
 XSecAlgorithmI ()
 
 XSecAlgorithmI (string name)
 
 XSecAlgorithmI (string name, string config)
 
- Protected Member Functions inherited from genie::Algorithm
 Algorithm ()
 
 Algorithm (string name)
 
 Algorithm (string name, string config)
 
void Initialize (void)
 
void DeleteConfig (void)
 
void DeleteSubstructure (void)
 
RegistryExtractLocalConfig (const Registry &in) const
 
RegistryExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...
 
template<class T >
bool GetParam (const RgKey &name, T &p, bool is_top_call=true) const
 
template<class T >
bool GetParamDef (const RgKey &name, T &p, const T &def) const
 
template<class T >
int GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
 Handle to load vectors of parameters. More...
 
int GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
 
int AddTopRegistry (Registry *rp, bool owns=true)
 add registry with top priority, also update ownership More...
 
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership More...
 
int MergeTopRegistry (const Registry &r)
 
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships. More...
 
- Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
 
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...) More...
 
AlgId fID
 algorithm name and configuration set More...
 
vector< Registry * > fConfVect
 
vector< boolfOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

Computes neutrino-nucleon(nucleus) QELCC differential cross section Is a concrete implementation of the XSecAlgorithmI interface.
.

C.H.Llewellyn Smith, Physics Reports (Sect. C of Physics Letters) 3, No. 5 (1972) 261-379

Author
Costas Andreopoulos <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

May 05, 2004

Copyright (c) 2003-2020, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 37 of file LwlynSmithQELCCPXSec.h.

Constructor & Destructor Documentation

LwlynSmithQELCCPXSec::LwlynSmithQELCCPXSec ( )

Definition at line 41 of file LwlynSmithQELCCPXSec.cxx.

41  :
42 XSecAlgorithmI("genie::LwlynSmithQELCCPXSec")
43 {
44 
45 }
LwlynSmithQELCCPXSec::LwlynSmithQELCCPXSec ( string  config)

Definition at line 47 of file LwlynSmithQELCCPXSec.cxx.

47  :
48 XSecAlgorithmI("genie::LwlynSmithQELCCPXSec", config)
49 {
50 
51 }
static Config * config
Definition: config.cpp:1054
LwlynSmithQELCCPXSec::~LwlynSmithQELCCPXSec ( )
virtual

Definition at line 53 of file LwlynSmithQELCCPXSec.cxx.

54 {
55 
56 }

Member Function Documentation

void LwlynSmithQELCCPXSec::Configure ( const Registry config)
virtual

Configure the algorithm with an external registry The registry is merged with the top level registry if it is owned, Otherwise a copy of it is added with the highest priority

Reimplemented from genie::Algorithm.

Definition at line 402 of file LwlynSmithQELCCPXSec.cxx.

403 {
404  Algorithm::Configure(config);
405  this->LoadConfig();
406 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void LwlynSmithQELCCPXSec::Configure ( string  config)
virtual

Configure the algorithm from the AlgoConfigPool based on param_set string given in input An algorithm contains a vector of registries coming from different xml configuration files, which are loaded according a very precise prioriy This methods will load a number registries in order of priority: 1) "Tunable" parameter set from CommonParametes. This is loaded with the highest prioriry and it is designed to be used for tuning procedure Usage not expected from the user. 2) For every string defined in "CommonParame" the corresponding parameter set will be loaded from CommonParameter.xml 3) parameter set specified by the config string and defined in the xml file of the algorithm 4) if config is not "Default" also the Default parameter set from the same xml file will be loaded Effectively this avoids the repetion of a parameter when it is not changed in the requested configuration

Reimplemented from genie::Algorithm.

Definition at line 408 of file LwlynSmithQELCCPXSec.cxx.

409 {
411  this->LoadConfig();
412 }
static Config * config
Definition: config.cpp:1054
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
double LwlynSmithQELCCPXSec::FullDifferentialXSec ( const Interaction i) const
private

Definition at line 164 of file LwlynSmithQELCCPXSec.cxx.

166 {
167  // First we need access to all of the particles in the interaction
168  // The particles were stored in the lab frame
169  const Kinematics& kinematics = interaction -> Kine();
170  const InitialState& init_state = interaction -> InitState();
171 
172  const Target& tgt = init_state.Tgt();
173 
174  const TLorentzVector leptonMom = kinematics.FSLeptonP4();
175  const TLorentzVector outNucleonMom = kinematics.HadSystP4();
176 
177  // Apply Pauli blocking if enabled
178  if ( fDoPauliBlocking && tgt.IsNucleus() && !interaction->TestBit(kIAssumeFreeNucleon) ) {
179  int final_nucleon_pdg = interaction->RecoilNucleonPdg();
180  double kF = fPauliBlocker->GetFermiMomentum(tgt, final_nucleon_pdg,
181  tgt.HitNucPosition());
182  double pNf = outNucleonMom.P();
183  if ( pNf < kF ) return 0.;
184  }
185 
186  // Note that GetProbeP4 defaults to returning the probe 4-momentum in the
187  // struck nucleon rest frame, so we have to explicitly ask for the lab frame
188  // here
189  TLorentzVector * tempNeutrino = init_state.GetProbeP4(kRfLab);
190  TLorentzVector neutrinoMom = *tempNeutrino;
191  delete tempNeutrino;
192  TLorentzVector * inNucleonMom = init_state.TgtPtr()->HitNucP4Ptr();
193 
194  // *** CALCULATION OF "q" and "qTilde" ***
195  // According to the de Forest prescription for handling the off-shell
196  // initial struck nucleon, the cross section calculation should proceed
197  // as if for a free nucleon, except that an effective value of the 4-momentum
198  // transfer qTilde should be used in which the difference between the on-
199  // and off-shell energies of the hit nucleon has been subtracted from the
200  // energy transfer q0.
201 
202  // HitNucMass() looks up the PDGLibrary (on-shell) value for the initial
203  // struck nucleon
204  double mNi = init_state.Tgt().HitNucMass();
205 
206  // Hadronic matrix element for CC neutrino interactions should really use
207  // the "nucleon mass," i.e., the mean of the proton and neutrino masses.
208  // This expression would also work for NC and EM scattering (since the
209  // initial and final on-shell nucleon masses would be the same)
210  double mNucleon = ( mNi + interaction->RecoilNucleon()->Mass() ) / 2.;
211 
212  // Ordinary 4-momentum transfer
213  TLorentzVector qP4 = neutrinoMom - leptonMom;
214 
215  // Initial struck nucleon 4-momentum in which it is forced to be on-shell
216  double inNucleonOnShellEnergy = std::sqrt( std::pow(mNi, 2)
217  + std::pow(inNucleonMom->P(), 2) );
218  TLorentzVector inNucleonMomOnShell( inNucleonMom->Vect(), inNucleonOnShellEnergy );
219 
220  // Effective 4-momentum transfer (according to the deForest prescription) for
221  // use in computing the hadronic tensor
222  TLorentzVector qTildeP4 = outNucleonMom - inNucleonMomOnShell;
223 
224  double Q2 = -1. * qP4.Mag2();
225  double Q2tilde = -1. * qTildeP4.Mag2();
226 
227  // If the binding energy correction causes an unphysical value
228  // of q0Tilde or Q2tilde, just return 0.
229  if ( qTildeP4.E() <= 0. && init_state.Tgt().IsNucleus() &&
230  !interaction->TestBit(kIAssumeFreeNucleon) ) return 0.;
231  if ( Q2tilde <= 0. ) return 0.;
232 
233  // Store Q2tilde in the kinematic variable representing Q2.
234  // This will ensure that the form factors are calculated correctly
235  // using the de Forest prescription (Q2tilde instead of Q2).
236  interaction->KinePtr()->SetQ2(Q2tilde);
237 
238  // Calculate the QEL form factors
240 
241  double F1V = fFormFactors.F1V();
242  double xiF2V = fFormFactors.xiF2V();
243  double FA = fFormFactors.FA();
244  double Fp = fFormFactors.Fp();
245 
246  // Restore Q2 in the interaction's kinematic variables
247  // now that the form factors have been computed
248  interaction->KinePtr()->SetQ2( Q2 );
249 
250  // Overall factor in the differential cross section
251  double Gfactor = kGF2*fCos8c2 / ( 8. * kPi * kPi * inNucleonOnShellEnergy
252  * neutrinoMom.E() * outNucleonMom.E() * leptonMom.E() );
253 
254  // Now, we can calculate the cross section
255  double tau = Q2tilde / (4 * std::pow(mNucleon, 2));
256  double h1 = FA*FA*(1 + tau) + tau*(F1V + xiF2V)*(F1V + xiF2V);
257  double h2 = FA*FA + F1V*F1V + tau*xiF2V*xiF2V;
258  double h3 = 2.0 * FA * (F1V + xiF2V);
259  double h4 = 0.25 * xiF2V*xiF2V *(1-tau) + 0.5*F1V*xiF2V + FA*Fp - tau*Fp*Fp;
260 
261  bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
262  int sign = (is_neutrino) ? -1 : 1;
263  double l1 = 2*neutrinoMom.Dot(leptonMom)*std::pow(mNucleon, 2);
264  double l2 = 2*(neutrinoMom.Dot(inNucleonMomOnShell)) * (inNucleonMomOnShell.Dot(leptonMom)) - neutrinoMom.Dot(leptonMom)*std::pow(mNucleon, 2);
265  double l3 = (neutrinoMom.Dot(inNucleonMomOnShell) * qTildeP4.Dot(leptonMom)) - (neutrinoMom.Dot(qTildeP4) * leptonMom.Dot(inNucleonMomOnShell));
266  l3 *= sign;
267  double l4 = neutrinoMom.Dot(leptonMom) * qTildeP4.Dot(qTildeP4) - 2*neutrinoMom.Dot(qTildeP4)*leptonMom.Dot(qTildeP4);
268  double l5 = neutrinoMom.Dot(inNucleonMomOnShell) * leptonMom.Dot(qTildeP4) + leptonMom.Dot(inNucleonMomOnShell)*neutrinoMom.Dot(qTildeP4) - neutrinoMom.Dot(leptonMom)*inNucleonMomOnShell.Dot(qTildeP4);
269 
270  double LH = 2 *(l1*h1 + l2*h2 + l3*h3 + l4*h4 + l5*h2);
271 
272  double xsec = Gfactor * LH;
273 
274  // Apply the factor that arises from elimination of the energy-conserving
275  // delta function
277 
278  // Apply given scaling factor
279  xsec *= fXSecScale;
280 
281  // Number of scattering centers in the target
282  const Target & target = init_state.Tgt();
283  int nucpdgc = target.HitNucPdg();
284  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
285 
286 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
287  LOG("LwlynSmith", pDEBUG)
288  << "Nuclear suppression factor R(Q2) = " << R << ", NNucl = " << NNucl;
289 #endif
290 
291  xsec *= NNucl; // nuclear xsec
292 
293  return xsec;
294 
295 }
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
int HitNucPdg(void) const
Definition: Target.cxx:304
double HitNucMass(void) const
Definition: Target.cxx:233
bool IsNucleus(void) const
Definition: Target.cxx:272
constexpr T pow(T x)
Definition: pow.h:72
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:66
double fXSecScale
external xsec scaling factor
double GetFermiMomentum(const Target &tgt, int pdg_Nf, double radius=0.0) const
Get the Fermi momentum needed to check Pauli blocking.
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:333
const TLorentzVector & FSLeptonP4(void) const
Definition: Kinematics.h:65
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double EnergyDeltaFunctionSolutionQEL(const Interaction &inter)
Definition: QELUtils.cxx:50
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
int ProbePdg(void) const
Definition: InitialState.h:64
int Z(void) const
Definition: Target.h:68
bool fDoPauliBlocking
Whether to apply Pauli blocking in FullDifferentialXSec.
double xiF2V(void) const
Get the computed form factor xi*F2V.
void Calculate(const Interaction *interaction)
Compute the form factors for the input interaction using the attached model.
int sign(double val)
Definition: UtilFunc.cxx:104
int N(void) const
Definition: Target.h:69
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
double HitNucPosition(void) const
Definition: Target.h:89
Target * TgtPtr(void) const
Definition: InitialState.h:67
double fCos8c2
cos^2(cabibbo angle)
double F1V(void) const
Get the computed form factor F1V.
const Target & Tgt(void) const
Definition: InitialState.h:66
static const double kGF2
Definition: Constants.h:59
double Fp(void) const
Get the computed form factor Fp.
const genie::PauliBlocker * fPauliBlocker
The PauliBlocker instance to use to apply that correction.
static const double kPi
Definition: Constants.h:37
double FA(void) const
Get the computed form factor FA.
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
double LwlynSmithQELCCPXSec::Integral ( const Interaction i) const
virtual

Integrate the model over the kinematic phase space available to the input interaction (kinematical cuts can be included)

Implements genie::XSecAlgorithmI.

Definition at line 297 of file LwlynSmithQELCCPXSec.cxx.

298 {
299  // If we're using the new spline generation method (which integrates
300  // over the kPSQELEvGen phase space used by QELEventGenerator) then
301  // let the cross section integrator do all of the work. It's smart
302  // enough to handle free nucleon vs. nuclear targets, different
303  // nuclear models (including the local Fermi gas model), etc.
304  // TODO: think about doing this in a better way
305  if ( fXSecIntegrator->Id().Name() == "genie::NewQELXSec" ) {
306  return fXSecIntegrator->Integrate(this, in);
307  }
308 
309  // Otherwise, use the old integration method (kept for use with
310  // the historical default G18_00x series of tunes)
311  bool nuclear_target = in->InitState().Tgt().IsNucleus();
312  if(!nuclear_target || !fDoAvgOverNucleonMomentum) {
313  return fXSecIntegrator->Integrate(this,in);
314  }
315 
316  double E = in->InitState().ProbeE(kRfHitNucRest);
317  if(fLFG || E < fEnergyCutOff) {
318  // clone the input interaction so as to tweak the
319  // hit nucleon 4-momentum in the averaging loop
320  Interaction in_curr(*in);
321 
322  // hit target
323  Target * tgt = in_curr.InitState().TgtPtr();
324 
325  // get nuclear masses (init & final state nucleus)
326  int nucleon_pdgc = tgt->HitNucPdg();
327  bool is_p = pdg::IsProton(nucleon_pdgc);
328  int Zi = tgt->Z();
329  int Ai = tgt->A();
330  int Zf = (is_p) ? Zi-1 : Zi;
331  int Af = Ai-1;
332  PDGLibrary * pdglib = PDGLibrary::Instance();
333  TParticlePDG * nucl_i = pdglib->Find( pdg::IonPdgCode(Ai, Zi) );
334  TParticlePDG * nucl_f = pdglib->Find( pdg::IonPdgCode(Af, Zf) );
335  if(!nucl_f) {
336  LOG("LwlynSmith", pFATAL)
337  << "Unknwown nuclear target! No target with code: "
338  << pdg::IonPdgCode(Af, Zf) << " in PDGLibrary!";
339  exit(1);
340  }
341  double Mi = nucl_i -> Mass(); // initial nucleus mass
342  double Mf = nucl_f -> Mass(); // remnant nucleus mass
343 
344  // throw nucleons with fermi momenta and binding energies
345  // generated according to the current nuclear model for the
346  // input target and average the cross section
347  double xsec_sum = 0.;
348  const int nnuc = 2000;
349  // VertexGenerator for generating a position before generating
350  // each nucleon
351  VertexGenerator * vg = new VertexGenerator();
352  vg->Configure("Default");
353  for(int inuc=0;inuc<nnuc;inuc++){
354  // Generate a position in the nucleus
355  TVector3 nucpos = vg->GenerateVertex(&in_curr,tgt->A());
356  tgt->SetHitNucPosition(nucpos.Mag());
357 
358  // Generate a nucleon
359  fNuclModel->GenerateNucleon(*tgt, nucpos.Mag());
360  TVector3 p3N = fNuclModel->Momentum3();
361  double EN = Mi - TMath::Sqrt(p3N.Mag2() + Mf*Mf);
362  TLorentzVector* p4N = tgt->HitNucP4Ptr();
363  p4N->SetPx (p3N.Px());
364  p4N->SetPy (p3N.Py());
365  p4N->SetPz (p3N.Pz());
366  p4N->SetE (EN);
367 
368  double xsec = fXSecIntegrator->Integrate(this,&in_curr);
369  xsec_sum += xsec;
370  }
371  double xsec_avg = xsec_sum / nnuc;
372  delete vg;
373  return xsec_avg;
374  }else{
375  return fXSecIntegrator->Integrate(this,in);
376  }
377 }
TVector3 GenerateVertex(const Interaction *in, double A) const
int HitNucPdg(void) const
Definition: Target.cxx:304
bool fLFG
If the nuclear model is lfg alway average over nucleons.
int A(void) const
Definition: Target.h:70
#define pFATAL
Definition: Messenger.h:56
const NuclearModelI * fNuclModel
double Mass(Resonance_t res)
resonance mass (GeV)
void SetHitNucPosition(double r)
Definition: Target.cxx:210
bool fDoAvgOverNucleonMomentum
Average cross section over hit nucleon monentum?
const TVector3 & Momentum3(void) const
Definition: NuclearModelI.h:75
Summary information for an interaction.
Definition: Interaction.h:56
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:333
const XSecIntegratorI * fXSecIntegrator
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
string Name(void) const
Definition: AlgId.h:44
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
int Z(void) const
Definition: Target.h:68
void Configure(const Registry &config)
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:97
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:32
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:68
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
virtual bool GenerateNucleon(const Target &) const =0
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
void LwlynSmithQELCCPXSec::LoadConfig ( void  )
private

Definition at line 414 of file LwlynSmithQELCCPXSec.cxx.

415 {
416  // Cross section scaling factor
417  GetParamDef( "QEL-CC-XSecScale", fXSecScale, 1. ) ;
418 
419  double thc ;
420  GetParam( "CabibboAngle", thc ) ;
421  fCos8c2 = TMath::Power(TMath::Cos(thc), 2);
422 
423  // load QEL form factors model
424  fFormFactorsModel = dynamic_cast<const QELFormFactorsModelI *> (
425  this->SubAlg("FormFactorsAlg"));
426  assert(fFormFactorsModel);
427  fFormFactors.SetModel(fFormFactorsModel); // <-- attach algorithm
428 
429  // load XSec Integrator
431  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
432  assert(fXSecIntegrator);
433 
434  // Get nuclear model for use in Integral()
435  RgKey nuclkey = "IntegralNuclearModel";
436  fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
437  assert(fNuclModel);
438 
440 
441  bool average_over_nuc_mom ;
442  GetParamDef( "IntegralAverageOverNucleonMomentum", average_over_nuc_mom, false ) ;
443  // Always average over initial nucleons if the nuclear model is LFG
444  fDoAvgOverNucleonMomentum = fLFG || average_over_nuc_mom ;
445 
446  fEnergyCutOff = 0.;
447 
448  // Get averaging cutoff energy
449  GetParamDef("IntegralNuclearInfluenceCutoffEnergy", fEnergyCutOff, 2.0 ) ;
450 
451  // Method to use to calculate the binding energy of the initial hit nucleon when
452  // generating splines
453  std::string temp_binding_mode;
454  GetParamDef( "IntegralNucleonBindingMode", temp_binding_mode, std::string("UseNuclearModel") );
456 
457  // Get PauliBlocker for possible use in FullDifferentialXSec()
458  GetParamDef( "IntegralNucleonBindingMode", temp_binding_mode, std::string("UseNuclearModel") );
459  fPauliBlocker = dynamic_cast<const PauliBlocker*>( this->SubAlg("PauliBlockerAlg") );
460  assert( fPauliBlocker );
461 
462  // Decide whether or not it should be used in FullDifferentialXSec
463  GetParamDef( "DoPauliBlocking", fDoPauliBlocking, true );
464 }
Cross Section Integrator Interface.
std::string string
Definition: nybbler.cc:12
bool fLFG
If the nuclear model is lfg alway average over nucleons.
void SetModel(const QELFormFactorsModelI *model)
Attach an algorithm.
Examines whether the generated event should be Pauli blocked. Is a concerete implementation of the Ev...
Definition: PauliBlocker.h:36
const NuclearModelI * fNuclModel
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
double fXSecScale
external xsec scaling factor
bool fDoAvgOverNucleonMomentum
Average cross section over hit nucleon monentum?
QELEvGen_BindingMode_t fIntegralNucleonBindingMode
virtual NuclearModel_t ModelType(const Target &) const =0
const QELFormFactorsModelI * fFormFactorsModel
const XSecIntegratorI * fXSecIntegrator
Pure abstract base class. Defines the QELFormFactorsModelI interface to be implemented by any algorit...
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
bool fDoPauliBlocking
Whether to apply Pauli blocking in FullDifferentialXSec.
string RgKey
QELEvGen_BindingMode_t StringToQELBindingMode(const std::string &mode_str)
Definition: QELUtils.cxx:194
double fCos8c2
cos^2(cabibbo angle)
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 genie::PauliBlocker * fPauliBlocker
The PauliBlocker instance to use to apply that correction.
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
bool LwlynSmithQELCCPXSec::ValidProcess ( const Interaction i) const
virtual

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 379 of file LwlynSmithQELCCPXSec.cxx.

380 {
381  if(interaction->TestBit(kISkipProcessChk)) return true;
382 
383  const InitialState & init_state = interaction->InitState();
384  const ProcessInfo & proc_info = interaction->ProcInfo();
385 
386  if(!proc_info.IsQuasiElastic()) return false;
387 
388  int nuc = init_state.Tgt().HitNucPdg();
389  int nu = init_state.ProbePdg();
390 
391  bool isP = pdg::IsProton(nuc);
392  bool isN = pdg::IsNeutron(nuc);
393  bool isnu = pdg::IsNeutrino(nu);
394  bool isnub = pdg::IsAntiNeutrino(nu);
395 
396  bool prcok = proc_info.IsWeakCC() && ((isP&&isnub) || (isN&&isnu));
397  if(!prcok) return false;
398 
399  return true;
400 }
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
int HitNucPdg(void) const
Definition: Target.cxx:304
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:338
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:333
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:115
int ProbePdg(void) const
Definition: InitialState.h:64
const Target & Tgt(void) const
Definition: InitialState.h:66
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48
double LwlynSmithQELCCPXSec::XSec ( const Interaction i,
KinePhaseSpace_t  k 
) const
virtual

Compute the cross section for the input interaction.

Implements genie::XSecAlgorithmI.

Definition at line 58 of file LwlynSmithQELCCPXSec.cxx.

60 {
61  if(! this -> ValidProcess (interaction) ) {LOG("LwlynSmith",pWARN) << "not a valid process"; return 0.;}
62  if(! this -> ValidKinematics (interaction) ) {LOG("LwlynSmith",pWARN) << "not valid kinematics"; return 0.;}
63 
64  // If computing the full differential cross section, then all four momentum
65  // four-vectors (probe, hit nucleon, final lepton, and final nucleon) should
66  // have been set already, with the hit nucleon off-shell as appropriate.
67  if (kps == kPSQELEvGen) {
68  return this->FullDifferentialXSec(interaction);
69  }
70 
71  // Get kinematics & init-state parameters
72  const Kinematics & kinematics = interaction -> Kine();
73  const InitialState & init_state = interaction -> InitState();
74  const Target & target = init_state.Tgt();
75 
76  double E = init_state.ProbeE(kRfHitNucRest);
77  double E2 = TMath::Power(E,2);
78  double ml = interaction->FSPrimLepton()->Mass();
79  double M = target.HitNucMass();
80  double q2 = kinematics.q2();
81 
82  // One of the xsec terms changes sign for antineutrinos
83  bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
84  int sign = (is_neutrino) ? -1 : 1;
85 
86  // Calculate the QEL form factors
88 
89  double F1V = fFormFactors.F1V();
90  double xiF2V = fFormFactors.xiF2V();
91  double FA = fFormFactors.FA();
92  double Fp = fFormFactors.Fp();
93 
94 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
95  LOG("LwlynSmith", pDEBUG) << "\n" << fFormFactors;
96 #endif
97 
98  // Calculate auxiliary parameters
99  double ml2 = TMath::Power(ml, 2);
100  double M2 = TMath::Power(M, 2);
101  double M4 = TMath::Power(M2, 2);
102  double FA2 = TMath::Power(FA, 2);
103  double Fp2 = TMath::Power(Fp, 2);
104  double F1V2 = TMath::Power(F1V, 2);
105  double xiF2V2 = TMath::Power(xiF2V, 2);
106  double Gfactor = M2*kGF2*fCos8c2 / (8*kPi*E2);
107  double s_u = 4*E*M + q2 - ml2;
108  double q2_M2 = q2/M2;
109 
110  // Compute free nucleon differential cross section
111  double A = (0.25*(ml2-q2)/M2) * (
112  (4-q2_M2)*FA2 - (4+q2_M2)*F1V2 - q2_M2*xiF2V2*(1+0.25*q2_M2)
113  -4*q2_M2*F1V*xiF2V - (ml2/M2)*(
114  (F1V2+xiF2V2+2*F1V*xiF2V)+(FA2+4*Fp2+4*FA*Fp)+(q2_M2-4)*Fp2));
115  double B = -1 * q2_M2 * FA*(F1V+xiF2V);
116  double C = 0.25*(FA2 + F1V2 - 0.25*q2_M2*xiF2V2);
117 
118  double xsec = Gfactor * (A + sign*B*s_u/M2 + C*s_u*s_u/M4);
119 
120  // Apply given scaling factor
121  xsec *= fXSecScale;
122 
123 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
124  LOG("LwlynSmith", pDEBUG)
125  << "dXSec[QEL]/dQ2 [FreeN](E = "<< E << ", Q2 = "<< -q2 << ") = "<< xsec;
126  LOG("LwlynSmith", pDEBUG)
127  << "A(Q2) = " << A << ", B(Q2) = " << B << ", C(Q2) = " << C;
128 #endif
129 
130  //----- The algorithm computes dxsec/dQ2
131  // Check whether variable tranformation is needed
132  if(kps!=kPSQ2fE) {
134 
135 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
136  LOG("LwlynSmith", pDEBUG)
137  << "Jacobian for transformation to: "
138  << KinePhaseSpace::AsString(kps) << ", J = " << J;
139 #endif
140  xsec *= J;
141  }
142 
143  //----- if requested return the free nucleon xsec even for input nuclear tgt
144  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
145 
146  //----- compute nuclear suppression factor
147  // (R(Q2) is adapted from NeuGEN - see comments therein)
148  double R = nuclear::NuclQELXSecSuppression("Default", 0.5, interaction);
149 
150  //----- number of scattering centers in the target
151  int nucpdgc = target.HitNucPdg();
152  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
153 
154 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
155  LOG("LwlynSmith", pDEBUG)
156  << "Nuclear suppression factor R(Q2) = " << R << ", NNucl = " << NNucl;
157 #endif
158 
159  xsec *= (R*NNucl); // nuclear xsec
160 
161  return xsec;
162 }
double FullDifferentialXSec(const Interaction *i) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
int HitNucPdg(void) const
Definition: Target.cxx:304
double HitNucMass(void) const
Definition: Target.cxx:233
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double fXSecScale
external xsec scaling factor
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:141
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:333
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static string AsString(KinePhaseSpace_t kps)
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
int ProbePdg(void) const
Definition: InitialState.h:64
int Z(void) const
Definition: Target.h:68
double xiF2V(void) const
Get the computed form factor xi*F2V.
#define pWARN
Definition: Messenger.h:60
void Calculate(const Interaction *interaction)
Compute the form factors for the input interaction using the attached model.
int sign(double val)
Definition: UtilFunc.cxx:104
int N(void) const
Definition: Target.h:69
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
double fCos8c2
cos^2(cabibbo angle)
Definition: 018_def.c:13
#define A
Definition: memgrp.cpp:38
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
double F1V(void) const
Get the computed form factor F1V.
const Target & Tgt(void) const
Definition: InitialState.h:66
static const double kGF2
Definition: Constants.h:59
double Fp(void) const
Get the computed form factor Fp.
double ProbeE(RefFrame_t rf) const
static const double kPi
Definition: Constants.h:37
double FA(void) const
Get the computed form factor FA.
double NuclQELXSecSuppression(string kftable, double pmax, const Interaction *in)
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63

Member Data Documentation

double genie::LwlynSmithQELCCPXSec::fCos8c2
private

cos^2(cabibbo angle)

Definition at line 62 of file LwlynSmithQELCCPXSec.h.

bool genie::LwlynSmithQELCCPXSec::fDoAvgOverNucleonMomentum
private

Average cross section over hit nucleon monentum?

Definition at line 69 of file LwlynSmithQELCCPXSec.h.

bool genie::LwlynSmithQELCCPXSec::fDoPauliBlocking
private

Whether to apply Pauli blocking in FullDifferentialXSec.

Definition at line 78 of file LwlynSmithQELCCPXSec.h.

double genie::LwlynSmithQELCCPXSec::fEnergyCutOff
private

Average only for energies below this cutoff defining the region where nuclear modeling details do matter

Definition at line 70 of file LwlynSmithQELCCPXSec.h.

QELFormFactors genie::LwlynSmithQELCCPXSec::fFormFactors
mutableprivate

Definition at line 59 of file LwlynSmithQELCCPXSec.h.

const QELFormFactorsModelI* genie::LwlynSmithQELCCPXSec::fFormFactorsModel
private

Definition at line 60 of file LwlynSmithQELCCPXSec.h.

QELEvGen_BindingMode_t genie::LwlynSmithQELCCPXSec::fIntegralNucleonBindingMode
private

Enum specifying the method to use when calculating the binding energy of the initial hit nucleon during spline generation

Definition at line 75 of file LwlynSmithQELCCPXSec.h.

bool genie::LwlynSmithQELCCPXSec::fLFG
private

If the nuclear model is lfg alway average over nucleons.

Definition at line 68 of file LwlynSmithQELCCPXSec.h.

const NuclearModelI* genie::LwlynSmithQELCCPXSec::fNuclModel
private

Definition at line 67 of file LwlynSmithQELCCPXSec.h.

const genie::PauliBlocker* genie::LwlynSmithQELCCPXSec::fPauliBlocker
private

The PauliBlocker instance to use to apply that correction.

Definition at line 80 of file LwlynSmithQELCCPXSec.h.

const XSecIntegratorI* genie::LwlynSmithQELCCPXSec::fXSecIntegrator
private

Definition at line 61 of file LwlynSmithQELCCPXSec.h.

double genie::LwlynSmithQELCCPXSec::fXSecScale
private

external xsec scaling factor

Definition at line 64 of file LwlynSmithQELCCPXSec.h.


The documentation for this class was generated from the following files: