NievesQELCCPXSec.h
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \class genie::NievesQELCCPXSec
5 
6 \brief Computes neutrino-nucleon(nucleus) QELCC differential cross section
7  with RPA corrections
8  Is a concrete implementation of the XSecAlgorithmI interface. \n
9 
10 \ref Physical Review C 70, 055503 (2004)
11 
12 \author Joe Johnston, University of Pittsburgh
13  Steven Dytman, University of Pittsburgh
14 
15 \created April 2016
16 
17 \cpright Copyright (c) 2003-2019, The GENIE Collaboration
18  For the full text of the license visit http://copyright.genie-mc.org
19  or see $GENIE/LICENSE
20 */
21 //____________________________________________________________________________
22 
23 #ifndef _NIEVES_QELCC_CROSS_SECTION_H_
24 #define _NIEVES_QELCC_CROSS_SECTION_H_
25 
29 #include <complex>
30 #include <Math/IFunction.h>
34 
35 namespace genie {
36 
37 typedef enum EQELRmax {
38  // Use the same maximum radius as VertexGenerator (3*R0*A^(1/3))
40 
41  // Use the method for calculting Rmax from Nieves' Fortran code
44 
46 class XSecIntegratorI;
47 
49 
50 public:
52  NievesQELCCPXSec(string config);
53  virtual ~NievesQELCCPXSec();
54 
55  // XSecAlgorithmI interface implementation
56  double XSec (const Interaction * i, KinePhaseSpace_t k) const;
57  double Integral (const Interaction * i) const;
58  bool ValidProcess (const Interaction * i) const;
59 
60  // Override the Algorithm::Configure methods to load configuration
61  // data to private data members
62  void Configure (const Registry & config);
63  void Configure (string param_set);
64 
65 private:
66  void LoadConfig (void);
67 
71  double fCos8c2; ///< cos^2(cabibbo angle)
72 
73  double fXSecScale; ///< external xsec scaling factor
74 
75  double fhbarc; ///< hbar*c in GeV*fm
76 
77  // mutable for testing purposes only!
78  mutable bool fRPA; ///< use RPA corrections
79  bool fCoulomb; ///< use Coulomb corrections
80 
81  const NuclearModelI* fNuclModel; ///< Nuclear Model for integration
82  // Detect whether the nuclear model is local Fermi gas, and store
83  // the relativistic Fermi momentum table if not
84  bool fLFG;
86  string fKFTableName;
87 
88  /// Enum specifying the method to use when calculating the binding energy of
89  /// the initial hit nucleon during spline generation
91 
92  /// Cutoff lab-frame probe energy above which the effects of Fermi motion and
93  /// binding energy are ignored when computing the total cross section
94  double fEnergyCutOff;
95 
96  /// Whether to apply Pauli blocking in XSec()
98  /// The PauliBlocker instance to use to apply that correction
100 
101  /// Nuclear radius parameter r = R0*A^(1/3) used to compute the
102  /// maximum radius for integration of the Coulomb potential
103  /// when matching the VertexGenerator method
104  double fR0;
105 
106  /// Enum variable describing which method of computing Rmax should be used
107  /// for integrating the Coulomb potential
108  Nieves_Coulomb_Rmax_t fCoulombRmaxMode;
109 
110  //Functions needed to calculate XSec:
111 
112  // Calculates values of CN, CT, CL, and imU, and stores them in the provided
113  // variables. If target is not a nucleus, then CN, CN, and CL are all 1.0.
114  // r must be in units of fm.
115  void CNCTCLimUcalc(TLorentzVector qTildeP4, double M, double r,
116  bool is_neutrino, bool tgtIsNucleus, int tgt_pdgc, int A, int Z, int N,
117  bool hitNucIsProton, double & CN, double & CT, double & CL, double & imU,
118  double & t0, double & r00, bool assumeFreeNucleon) const;
119 
120  //Equations to calculate the relativistic Lindhard function for Amunu
121  std::complex<double> relLindhardIm(double q0gev, double dqgev,
122  double kFngev, double kFpgev,
123  double M, bool isNeutrino,
124  double & t0, double & r00) const;
125  std::complex<double> relLindhard(double q0gev, double dqgev,
126  double kFgev, double M,
127  bool isNeutrino,
128  std::complex<double> relLindIm) const;
129  std::complex<double> ruLinRelX(double q0, double qm,
130  double kf, double m) const;
131  std::complex<double> deltaLindhard(double q0gev, double dqgev,
132  double rho, double kFgev) const;
133 
134  // Potential for coulomb correction
135  double vcr(const Target * target, double r) const;
136 
137  //input must be length 4. Returns 1 if input is an even permutation of 0123,
138  //-1 if input is an odd permutation of 0123, and 0 if any two elements
139  //are equal
140  int leviCivita(int input[]) const;
141 
142  double LmunuAnumu(const TLorentzVector neutrinoMom,
143  const TLorentzVector inNucleonMom, const TLorentzVector leptonMom,
144  const TLorentzVector outNucleonMom, double M, bool is_neutrino,
145  const Target& target, bool assumeFreeNucleon) const;
146 
147  // NOTE: THE FOLLOWING CODE IS FOR TESTING PURPOSES ONLY
148  // Used to print tensor elements and various inputs for comparison to Nieves'
149  // fortran code
150  mutable bool fCompareNievesTensors; ///< print tensors
151  mutable TString fTensorsOutFile; ///< file to print tensors to
152  mutable double fVc,fCoulombFactor;
153  void CompareNievesTensors(const Interaction* i) const;
154  // END TESTING CODE
155 };
156 } // genie namespace
157 
158 //____________________________________________________________________________
159 /*!
160 \class genie::utils::gsl::wrap::NievesQELIntegrand
161 
162 \brief Auxiliary scalar function for integration over the nuclear density
163  when calculaing the Coulomb correction in the Nieves QEL xsec model
164 
165 \author Joe Johnston, University of Pittsburgh
166  Steven Dytman, University of Pittsburgh
167 
168 \created June 03, 2016
169 */
170 //____________________________________________________________________________
171 
172 namespace genie {
173  namespace utils {
174  namespace gsl {
175  namespace wrap {
176 
177  class NievesQELvcrIntegrand : public ROOT::Math::IBaseFunctionOneDim
178  {
179  public:
180  NievesQELvcrIntegrand(double Rcurr, int A, int Z);
182  // ROOT::Math::IBaseFunctionOneDim interface
183  unsigned int NDim (void) const;
184  double DoEval (double rin) const;
185  ROOT::Math::IBaseFunctionOneDim * Clone (void) const;
186  private:
187  double fRcurr;
188  double fA;
189  double fZ;
190  };
191 
192  } // wrap namespace
193  } // gsl namespace
194  } // utils namespace
195 } // genie namespace
196 
197 
198 #endif
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
bool fRPA
use RPA corrections
A class holding Quasi Elastic (QEL) Form Factors.
const XML_Char * target
Definition: expat.h:268
void Configure(const Registry &config)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
enum genie::EQELRmax Nieves_Coulomb_Rmax_t
QELEvGen_BindingMode_t fIntegralNucleonBindingMode
const FermiMomentumTable * fKFTable
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
Examines whether the generated event should be Pauli blocked. Is a concerete implementation of the Ev...
Definition: PauliBlocker.h:36
const XSecIntegratorI * fXSecIntegrator
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
A table of Fermi momentum constants.
Definition: config.py:1
QELFormFactors fFormFactors
bool fCoulomb
use Coulomb corrections
enum genie::EKinePhaseSpace KinePhaseSpace_t
enum genie::EQELEvGenBindingMode QELEvGen_BindingMode_t
Nieves_Coulomb_Rmax_t fCoulombRmaxMode
std::complex< double > deltaLindhard(double q0gev, double dqgev, double rho, double kFgev) const
const QELFormFactorsModelI * fFormFactorsModel
double fhbarc
hbar*c in GeV*fm
Float_t Z
Definition: plot.C:38
double fCos8c2
cos^2(cabibbo angle)
Summary information for an interaction.
Definition: Interaction.h:56
const NuclearModelI * fNuclModel
Nuclear Model for integration.
Pure abstract base class. Defines the QELFormFactorsModelI interface to be implemented by any algorit...
bool fCompareNievesTensors
print tensors
double Integral(const Interaction *i) const
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
bool fDoPauliBlocking
Whether to apply Pauli blocking in XSec()
double fXSecScale
external xsec scaling factor
Computes neutrino-nucleon(nucleus) QELCC differential cross section with RPA corrections Is a concret...
double vcr(const Target *target, double r) const
TString fTensorsOutFile
file to print tensors to
static const double A
Definition: Units.h:82
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
TRandom3 r(0)
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?
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.
std::complex< double > relLindhardIm(double q0gev, double dqgev, double kFngev, double kFpgev, double M, bool isNeutrino, double &t0, double &r00) const