GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
genie::NievesQELCCPXSec Class Reference

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

#include <NievesQELCCPXSec.h>

Inheritance diagram for genie::NievesQELCCPXSec:
Inheritance graph
[legend]
Collaboration diagram for genie::NievesQELCCPXSec:
Collaboration graph
[legend]

Public Member Functions

 NievesQELCCPXSec ()
 
 NievesQELCCPXSec (string config)
 
virtual ~NievesQELCCPXSec ()
 
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 Registry & GetConfig (void) const
 
Registry * GetOwnedConfig (void)
 
virtual const AlgId & Id (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 Algorithm * SubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

Private Member Functions

void LoadConfig (void)
 
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
 
std::complex< double > relLindhardIm (double q0gev, double dqgev, double kFngev, double kFpgev, double M, bool isNeutrino, double &t0, double &r00) const
 
std::complex< double > relLindhard (double q0gev, double dqgev, double kFgev, double M, bool isNeutrino, std::complex< double > relLindIm) const
 
std::complex< double > ruLinRelX (double q0, double qm, double kf, double m) const
 
std::complex< double > deltaLindhard (double q0gev, double dqgev, double rho, double kFgev) const
 
double vcr (const Target *target, double r) const
 
int leviCivita (int input[]) 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
 
void CompareNievesTensors (const Interaction *i) const
 

Private Attributes

QELFormFactors fFormFactors
 
const QELFormFactorsModelI * fFormFactorsModel
 
const XSecIntegratorI * fXSecIntegrator
 
double fCos8c2
 cos^2(cabibbo angle) More...
 
double fXSecCCScale
 external xsec scaling factor for CC More...
 
double fXSecNCScale
 external xsec scaling factor for NC More...
 
double fXSecEMScale
 external xsec scaling factor for EM More...
 
const QvalueShifter * fQvalueShifter
 Optional algorithm to retrieve the qvalue shift for a given target. More...
 
double fhbarc
 hbar*c in GeV*fm More...
 
bool fRPA
 use RPA corrections More...
 
bool fCoulomb
 use Coulomb corrections More...
 
const NuclearModelI * fNuclModel
 Nuclear Model for integration. More...
 
bool fLFG
 
const FermiMomentumTable * fKFTable
 
string fKFTableName
 
QELEvGen_BindingMode_t fIntegralNucleonBindingMode
 
double fEnergyCutOff
 
bool fDoPauliBlocking
 Whether to apply Pauli blocking in XSec() More...
 
const genie::PauliBlocker * fPauliBlocker
 The PauliBlocker instance to use to apply that correction. More...
 
double fR0
 
double fCoulombScale
 Scaling factor for the Coulomb potential. More...
 
Nieves_Coulomb_Rmax_t fCoulombRmaxMode
 
bool fCompareNievesTensors
 print tensors More...
 
TString fTensorsOutFile
 file to print tensors to More...
 
double fVc
 
double fCoulombFactor
 

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)
 
static string BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)
 
static string BuildParamMatRowSizeKey (const std::string &comm_name)
 
static string BuildParamMatColSizeKey (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)
 
Registry * ExtractLocalConfig (const Registry &in) const
 
Registry * ExtractLowerConfig (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
 
template<class T >
int GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
 Handle to load matrix of parameters. More...
 
template<class T >
int GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const
 
int GetParamMatKeys (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< bool > fOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMap * fOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

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

References:
Physical Review C 70, 055503 (2004)
Author
Joe Johnston, University of Pittsburgh Steven Dytman, University of Pittsburgh
Created:
April 2016
License:
Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 48 of file NievesQELCCPXSec.h.

Constructor & Destructor Documentation

NievesQELCCPXSec::NievesQELCCPXSec ( )

Definition at line 53 of file NievesQELCCPXSec.cxx.

53  :
54 XSecAlgorithmI("genie::NievesQELCCPXSec")
55 {
56 
57 }
NievesQELCCPXSec::NievesQELCCPXSec ( string  config)

Definition at line 59 of file NievesQELCCPXSec.cxx.

59  :
60 XSecAlgorithmI("genie::NievesQELCCPXSec", config)
61 {
62 
63 }
NievesQELCCPXSec::~NievesQELCCPXSec ( )
virtual

Definition at line 65 of file NievesQELCCPXSec.cxx.

66 {
67 
68  }

Member Function Documentation

void NievesQELCCPXSec::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
private

Definition at line 519 of file NievesQELCCPXSec.cxx.

References deltaLindhard(), genie::utils::prem::Density(), fhbarc, genie::FermiMomentumTable::FindClosestKF(), fKFTable, fLFG, genie::kPdgNeutron, genie::kPdgProton, genie::constants::kPi, genie::constants::kPi2, genie::constants::kPionMass2, relLindhard(), and relLindhardIm().

Referenced by LmunuAnumu().

523 {
524  if ( tgtIsNucleus && !assumeFreeNucleon ) {
525  double dq = qTildeP4.Vect().Mag();
526  double dq2 = TMath::Power(dq,2);
527  double q2 = 1 * qTildeP4.Mag2();
528  //Terms for polarization coefficients CN,CT, and CL
529  double hbarc2 = TMath::Power(fhbarc,2);
530  double c0 = 0.380/fhbarc;//Constant for CN in natural units
531 
532  //Density gives the nuclear density, normalized to 1
533  //Input radius r must be in fm
534  double rhop = nuclear::Density(r,A)*Z;
535  double rhon = nuclear::Density(r,A)*N;
536  double rho = rhop + rhon;
537  double rho0 = A*nuclear::Density(0,A);
538 
539  double fPrime = (0.33*rho/rho0+0.45*(1-rho/rho0))*c0;
540 
541  // Get Fermi momenta
542  double kF1, kF2;
543  if(fLFG){
544  if(hitNucIsProton){
545  kF1 = TMath::Power(3*kPi2*rhop, 1.0/3.0) *fhbarc;
546  kF2 = TMath::Power(3*kPi2*rhon, 1.0/3.0) *fhbarc;
547  }else{
548  kF1 = TMath::Power(3*kPi2*rhon, 1.0/3.0) *fhbarc;
549  kF2 = TMath::Power(3*kPi2*rhop, 1.0/3.0) *fhbarc;
550  }
551  }else{
552  if(hitNucIsProton){
553  kF1 = fKFTable->FindClosestKF(tgt_pdgc, kPdgProton);
554  kF2 = fKFTable->FindClosestKF(tgt_pdgc, kPdgNeutron);
555  }else{
556  kF1 = fKFTable->FindClosestKF(tgt_pdgc, kPdgNeutron);
557  kF2 = fKFTable->FindClosestKF(tgt_pdgc, kPdgProton);
558  }
559  }
560 
561  double kF = TMath::Power(1.5*kPi2*rho, 1.0/3.0) *fhbarc;
562 
563  std::complex<double> imU(relLindhardIm(qTildeP4.E(),dq,kF1,kF2,
564  M,is_neutrino,t0,r00));
565 
566  imaginaryU = imag(imU);
567 
568  std::complex<double> relLin(0,0),udel(0,0);
569 
570  // By comparison with Nieves' fortran code
571  if(imaginaryU < 0.){
572  relLin = relLindhard(qTildeP4.E(),dq,kF,M,is_neutrino,imU);
573  udel = deltaLindhard(qTildeP4.E(),dq,rho,kF);
574  }
575  std::complex<double> relLinTot(relLin + udel);
576 
577  /* CRho = 2
578  DeltaRho = 2500 MeV, (2.5 GeV)^2 = 6.25 GeV^2
579  mRho = 770 MeV, (0.770 GeV)^2 = 0.5929 GeV^2
580  g' = 0.63 */
581  double Vt = 0.08*4*kPi/kPionMass2 *
582  (2* TMath::Power((6.25-0.5929)/(6.25-q2),2)*dq2/(q2-0.5929) + 0.63);
583  /* f^2/4/Pi = 0.08
584  DeltaSubPi = 1200 MeV, (1.2 GeV)^2 = 1.44 GeV^2
585  g' = 0.63 */
586  double Vl = 0.08*4*kPi/kPionMass2 *
587  (TMath::Power((1.44-kPionMass2)/(1.44-q2),2)*dq2/(q2-kPionMass2)+0.63);
588 
589  CN = 1.0/TMath::Power(abs(1.0-fPrime*relLin/hbarc2),2);
590 
591  CT = 1.0/TMath::Power(abs(1.0-relLinTot*Vt),2);
592  CL = 1.0/TMath::Power(abs(1.0-relLinTot*Vl),2);
593  }
594  else {
595  //Polarization Coefficients: all equal to 1.0 for free nucleon
596  CN = 1.0;
597  CT = 1.0;
598  CL = 1.0;
599  imaginaryU = 0.0;
600  }
601 }
genie::NievesQELCCPXSec::fKFTable
const FermiMomentumTable * fKFTable
Definition: NievesQELCCPXSec.h:88
genie::NievesQELCCPXSec::deltaLindhard
std::complex< double > deltaLindhard(double q0gev, double dqgev, double rho, double kFgev) const
Definition: NievesQELCCPXSec.cxx:745
genie::NievesQELCCPXSec::fhbarc
double fhbarc
hbar*c in GeV*fm
Definition: NievesQELCCPXSec.h:78
genie::units::A
static constexpr double A
Definition: Units.h:74
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
genie::FermiMomentumTable::FindClosestKF
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
Definition: FermiMomentumTable.cxx:45
genie::NievesQELCCPXSec::relLindhard
std::complex< double > relLindhard(double q0gev, double dqgev, double kFgev, double M, bool isNeutrino, std::complex< double > relLindIm) const
Definition: NievesQELCCPXSec.cxx:662
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
genie::utils::prem::Density
double Density(double r)
Definition: PREM.cxx:18
genie::NievesQELCCPXSec::relLindhardIm
std::complex< double > relLindhardIm(double q0gev, double dqgev, double kFngev, double kFpgev, double M, bool isNeutrino, double &t0, double &r00) const
Definition: NievesQELCCPXSec.cxx:605
void NievesQELCCPXSec::CompareNievesTensors ( const Interaction *  i) const
private

Definition at line 1223 of file NievesQELCCPXSec.cxx.

References genie::Interaction::FSPrimLepton(), genie::Target::HitNucMass(), genie::pdg::IsNeutrino(), genie::Interaction::KinePtr(), LOG, pINFO, genie::InitialState::ProbePdg(), genie::Kinematics::SetQ2(), and genie::InitialState::Tgt().

1224  {
1225  Interaction * interaction = new Interaction(*in); // copy in
1226 
1227  // Set input values here
1228  double ein = 0.2;
1229  double ctl = 0.5;
1230  double rmaxfrac = 0.25;
1231 
1232  bool carbon = false; // true -> C12, false -> Pb208
1233 
1234  if(fRPA)
1235  fTensorsOutFile = "gen.RPA";
1236  else
1237  fTensorsOutFile = "gen.noRPA";
1238 
1239  // Calculate radius
1240  bool klave;
1241  double rp,ap,rn,an;
1242  if(carbon){
1243  klave = true;
1244  rp = 1.692;
1245  ap = 1.082;
1246  rn = 1.692;
1247  an = 1.082;
1248  }else{
1249  // Pb208
1250  klave = false;
1251  rp = 6.624;
1252  ap = 0.549;
1253  rn = 6.890;
1254  an = 0.549;
1255  }
1256  double rmax;
1257  if(!klave)
1258  rmax = TMath::Max(rp,rn) + 9.25*TMath::Max(ap,an);
1259  else
1260  rmax = TMath::Sqrt(20.0)*TMath::Max(rp,rn);
1261  double r = rmax * rmaxfrac;
1262 
1263  // Relevant objects and parameters
1264  //const Kinematics & kinematics = interaction -> Kine();
1265  const InitialState & init_state = interaction -> InitState();
1266  const Target & target = init_state.Tgt();
1267 
1268  // Parameters required for LmunuAnumu
1269  double M = target.HitNucMass();
1270  double ml = interaction->FSPrimLepton()->Mass();
1271  bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
1272 
1273  // Iterate over lepton energy (which then affects q, which is passed to
1274  // LmunuAnumu using in and out NucleonMom
1275  double delta = (ein-0.025)/100.0;
1276  for(int it=0;it<100;it++){
1277  double tmu = it*delta;
1278  double eout = ml + tmu;
1279  double pout = TMath::Sqrt(eout*eout-ml*ml);
1280 
1281  double pin = TMath::Sqrt(ein*ein); // Assume massless neutrinos
1282 
1283  double q0 = ein-eout;
1284  double dq = TMath::Sqrt(pin*pin+pout*pout-2.0*ctl*pin*pout);
1285  double q2 = q0*q0-dq*dq;
1286  interaction->KinePtr()->SetQ2(-q2);
1287 
1288  // When this method is called, inNucleonMomOnShell is unused.
1289  // I can thus provide the calculated values using a null vector for
1290  // inNucleonMomOnShell. I also need to put qTildeP4 in the z direction, as
1291  // Nieves does in his paper.
1292  TLorentzVector qTildeP4(0, 0, dq, q0);
1293  TLorentzVector inNucleonMomOnShell(0,0,0,0);
1294 
1295  // neutrinoMom and leptonMom only directly affect the leptonic tensor, which
1296  // we are not calculating now. Use them to transfer q.
1297  TLorentzVector neutrinoMom(0,0,pout+dq,eout+q0);
1298  TLorentzVector leptonMom(0,0,pout,eout);
1299 
1300  if(fCoulomb){ // Use same steps as in XSec()
1301  // Coulomb potential
1302  double Vc = vcr(& target, r);
1303  fVc = Vc;
1304 
1305  // Outgoing lepton energy and momentum including coulomb potential
1306  int sign = is_neutrino ? 1 : -1;
1307  double El = leptonMom.E();
1308  double ElLocal = El - sign*Vc;
1309  if(ElLocal - ml <= 0.0){
1310  LOG("Nieves",pINFO) << "Event should be rejected. Coulomb effects "
1311  << "push kinematics below threshold";
1312  return;
1313  }
1314  double plLocal = TMath::Sqrt(ElLocal*ElLocal-ml*ml);
1315 
1316  // Correction factor
1317  double coulombFactor= plLocal*ElLocal/leptonMom.Vect().Mag()/El;
1318  fCoulombFactor = coulombFactor; // Store and print
1319  }
1320 
1321  // TODO: apply Coulomb correction to 3-momentum transfer dq
1322 
1323  fFormFactors.Calculate(interaction);
1324  LmunuAnumu(neutrinoMom, inNucleonMomOnShell, leptonMom, qTildeP4,
1325  M, is_neutrino, target, false);
1326  }
1327  return;
1328 } // END TESTING CODE
genie::NievesQELCCPXSec::fRPA
bool fRPA
use RPA corrections
Definition: NievesQELCCPXSec.h:81
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
genie::Kinematics::SetQ2
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:255
genie::Interaction::KinePtr
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
genie::Target::HitNucMass
double HitNucMass(void) const
Definition: Target.cxx:233
genie::NievesQELCCPXSec::fFormFactors
QELFormFactors fFormFactors
Definition: NievesQELCCPXSec.h:68
genie::NievesQELCCPXSec::fCoulomb
bool fCoulomb
use Coulomb corrections
Definition: NievesQELCCPXSec.h:82
genie::Interaction
Summary information for an interaction.
Definition: Interaction.h:56
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
genie::InitialState::ProbePdg
int ProbePdg(void) const
Definition: InitialState.h:64
pINFO
#define pINFO
Definition: Messenger.h:62
genie::Interaction::FSPrimLepton
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
Definition: Interaction.cxx:126
genie::QELFormFactors::Calculate
void Calculate(const Interaction *interaction)
Compute the form factors for the input interaction using the attached model.
Definition: QELFormFactors.cxx:49
genie::NievesQELCCPXSec::vcr
double vcr(const Target *target, double r) const
Definition: NievesQELCCPXSec.cxx:832
genie::NievesQELCCPXSec::fTensorsOutFile
TString fTensorsOutFile
file to print tensors to
Definition: NievesQELCCPXSec.h:157
genie::InitialState::Tgt
const Target & Tgt(void) const
Definition: InitialState.h:66
genie::NievesQELCCPXSec::LmunuAnumu
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
Definition: NievesQELCCPXSec.cxx:925
genie::InitialState
Initial State information.
Definition: InitialState.h:48
void NievesQELCCPXSec::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 388 of file NievesQELCCPXSec.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

389 {
390  Algorithm::Configure(config);
391  this->LoadConfig();
392 }
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void NievesQELCCPXSec::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 394 of file NievesQELCCPXSec.cxx.

References genie::Algorithm::Configure(), and LoadConfig().

395 {
396  Algorithm::Configure(config);
397  this->LoadConfig();
398 }
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
std::complex< double > NievesQELCCPXSec::deltaLindhard ( double  q0gev,
double  dqgev,
double  rho,
double  kFgev 
) const
private

Definition at line 745 of file NievesQELCCPXSec.cxx.

References fhbarc, genie::constants::kPi, genie::units::m, genie::units::m2, and genie::units::s.

Referenced by CNCTCLimUcalc().

747 {
748  double q_zero = q0/fhbarc;
749  double q_mod = dq/fhbarc;
750  double k_fermi = kF/fhbarc;
751  //Divide by hbarc in order to use natural units (rho is already in the correct units)
752 
753  //m = 939/197.3, md = 1232/197.3, mpi = 139/197.3
754  double m = 4.7592;
755  double md = 6.2433;
756  double mpi = 0.7045;
757 
758  double fdel_f = 2.13;
759  double wr = md-m;
760  double gamma = 0;
761  double gammap = 0;
762 
763  double q_zero2 = TMath::Power(q_zero, 2);
764  double q_mod2 = TMath::Power(q_mod, 2);
765  double k_fermi2 = TMath::Power(k_fermi, 2);
766 
767  double m2 = TMath::Power(m, 2);
768  double m4 = TMath::Power(m, 4);
769  double mpi2 = TMath::Power(mpi, 2);
770  double mpi4 = TMath::Power(mpi, 4);
771 
772  double fdel_f2 = TMath::Power(fdel_f, 2);
773 
774  //For the current code q_zero is always real
775  //If q_zero can have an imaginary part then only the real part is used
776  //until z and zp are calculated
777 
778  double s = m2+q_zero2-q_mod2+
779  2.0*q_zero *TMath::Sqrt(m2+3.0/5.0*k_fermi2);
780 
781  if(s>TMath::Power(m+mpi,2)){
782  double srot = TMath::Sqrt(s);
783  double qcm = TMath::Sqrt(TMath::Power(s,2)+mpi4+m4-2.0*(s*mpi2+s*m2+
784  mpi2*m2)) /(2.0*srot);
785  gamma = 1.0/3.0 * 1.0/(4.0*kPi) * fdel_f2*
786  TMath::Power(qcm,3)/srot*(m+TMath::Sqrt(m2+TMath::Power(qcm,2)))/mpi2;
787  }
788  double sp = m2+q_zero2-q_mod2-
789  2.0*q_zero *TMath::Sqrt(m2+3.0/5.0*k_fermi2);
790 
791 
792  if(sp > TMath::Power(m+mpi,2)){
793  double srotp = TMath::Sqrt(sp);
794  double qcmp = TMath::Sqrt(TMath::Power(sp,2)+mpi4+m4-2.0*(sp*mpi2+sp*m2+
795  mpi2*m2))/(2.0*srotp);
796  gammap = 1.0/3.0 * 1.0/(4.0*kPi) * fdel_f2*
797  TMath::Power(qcmp,3)/srotp*(m+TMath::Sqrt(m2+TMath::Power(qcmp,2)))/mpi2;
798  }
799  //}//End if statement
800  const std::complex<double> iNum(0,1.0);
801 
802  std::complex<double> z(md/(q_mod*k_fermi)*(q_zero-q_mod2/(2.0*md)
803  -wr +iNum*gamma/2.0));
804  std::complex<double> zp(md/(q_mod*k_fermi)*(-q_zero-q_mod2/(2.0*md)
805  -wr +iNum*gammap/2.0));
806 
807  std::complex<double> pzeta(0.0);
808  if(abs(z) > 50.0){
809  pzeta = 2.0/(3.0*z)+2.0/(15.0*z*z*z);
810  }else if(abs(z) < TMath::Power(10.0,-2)){
811  pzeta = 2.0*z-2.0/3.0*z*z*z-iNum*kPi/2.0*(1.0-z*z);
812  }else{
813  pzeta = z + (1.0-z*z) * log((z+1.0)/(z-1.0))/2.0;
814  }
815 
816  std::complex<double> pzetap(0);
817  if(abs(zp) > 50.0){
818  pzetap = 2.0/(3.0*zp)+2.0/(15.0*zp*zp*zp);
819  }else if(abs(zp) < TMath::Power(10.0,-2)){
820  pzetap = 2.0*zp-2.0/3.0*zp*zp*zp-iNum*kPi/2.0*(1.0-zp*zp);
821  }else{
822  pzetap = zp+ (1.0-zp*zp) * log((zp+1.0)/(zp-1.0))/2.0;
823  }
824 
825  //Multiply by hbarc^2 to give answer in units of GeV^2
826  return 2.0/3.0 * rho * md/(q_mod*k_fermi) * (pzeta +pzetap) * fdel_f2 *
827  TMath::Power(fhbarc,2);
828 }
genie::units::s
static constexpr double s
Definition: Units.h:95
genie::NievesQELCCPXSec::fhbarc
double fhbarc
hbar*c in GeV*fm
Definition: NievesQELCCPXSec.h:78
genie::units::m2
static constexpr double m2
Definition: Units.h:72
genie::units::m
static constexpr double m
Definition: Units.h:71
double NievesQELCCPXSec::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 348 of file NievesQELCCPXSec.cxx.

References fXSecIntegrator, genie::Algorithm::Id(), genie::XSecIntegratorI::Integrate(), LOG, genie::AlgId::Name(), and pFATAL.

349 {
350  // If we're using the new spline generation method (which integrates
351  // over the kPSQELEvGen phase space used by QELEventGenerator) then
352  // let the cross section integrator do all of the work. It's smart
353  // enough to handle free nucleon vs. nuclear targets, different
354  // nuclear models (including the local Fermi gas model), etc.
355  if ( fXSecIntegrator->Id().Name() == "genie::NewQELXSec" ) {
356  return fXSecIntegrator->Integrate(this, in);
357  }
358  else {
359  LOG("Nieves", pFATAL) << "Splines for the Nieves CCQE model must be"
360  << " generated using genie::NewQELXSec";
361  std::exit(1);
362  }
363 }
pFATAL
#define pFATAL
Definition: Messenger.h:56
genie::NievesQELCCPXSec::fXSecIntegrator
const XSecIntegratorI * fXSecIntegrator
Definition: NievesQELCCPXSec.h:70
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::AlgId::Name
string Name(void) const
Definition: AlgId.h:44
genie::Algorithm::Id
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
genie::XSecIntegratorI::Integrate
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
int NievesQELCCPXSec::leviCivita ( int  input[]) const
private

Definition at line 891 of file NievesQELCCPXSec.cxx.

Referenced by LmunuAnumu().

891  {
892  int copy[4] = {input[0],input[1],input[2],input[3]};
893  int permutations = 0;
894  int temp;
895 
896  for(int i=0;i<4;i++){
897  for(int j=i+1;j<4;j++){
898  //If any two elements are equal return 0
899  if(input[i] == input[j])
900  return 0;
901  //If a larger number is before a smaller one, use permutations
902  //(exchanges of two adjacent elements) to move the smaller element
903  //so it is before the larger element, eg 2341->2314->2134->1234
904  if(copy[i]>copy[j]){
905  temp = copy[j];
906  for(int k=j;k>i;k--){
907  copy[k] = copy[k-1];
908  permutations++;
909  }
910  copy[i] = temp;
911  }
912  }
913  }
914 
915  if(permutations % 2 == 0){
916  return 1;
917  }else{
918  return -1;
919  }
920 }
double NievesQELCCPXSec::LmunuAnumu ( const TLorentzVector  neutrinoMom,
const TLorentzVector  inNucleonMom,
const TLorentzVector  leptonMom,
const TLorentzVector  outNucleonMom,
double  M,
bool  is_neutrino,
const Target &  target,
bool  assumeFreeNucleon 
) const
private

Definition at line 925 of file NievesQELCCPXSec.cxx.

References genie::Target::A(), genie::units::A, a, genie::units::b, CNCTCLimUcalc(), genie::QELFormFactors::F1V(), genie::QELFormFactors::FA(), fCompareNievesTensors, fCoulombFactor, fFormFactors, genie::QELFormFactors::Fp(), fRPA, fTensorsOutFile, fVc, genie::units::g, genie::Target::HitNucPdg(), genie::Target::HitNucPosition(), genie::Target::IsNucleus(), genie::pdg::IsProton(), genie::controls::kASmallNum, leviCivita(), LOG, genie::Target::N(), pDEBUG, genie::Target::Pdg(), pWARN, genie::QELFormFactors::xiF2V(), and genie::Target::Z().

Referenced by XSec().

929 {
930  double r = target.HitNucPosition();
931  bool tgtIsNucleus = target.IsNucleus();
932  int tgt_pdgc = target.Pdg();
933  int A = target.A();
934  int Z = target.Z();
935  int N = target.N();
936  bool hitNucIsProton = pdg::IsProton( target.HitNucPdg() );
937 
938  const double k[4] = {neutrinoMom.E(),neutrinoMom.Px(),neutrinoMom.Py(),neutrinoMom.Pz()};
939  const double kPrime[4] = {leptonMom.E(),leptonMom.Px(),
940  leptonMom.Py(),leptonMom.Pz()};
941 
942  double q2 = qTildeP4.Mag2();
943 
944  const double q[4] = {qTildeP4.E(),qTildeP4.Px(),qTildeP4.Py(),qTildeP4.Pz()};
945  double q0 = q[0];
946  double dq = TMath::Sqrt(TMath::Power(q[1],2)+
947  TMath::Power(q[2],2)+TMath::Power(q[3],2));
948 
949  int sign = (is_neutrino) ? 1 : -1;
950 
951  // Get the QEL form factors (were calculated before this method was called)
952  double F1V = 0.5*fFormFactors.F1V();
953  double xiF2V = 0.5*fFormFactors.xiF2V();
954  double FA = -fFormFactors.FA();
955  // According to Nieves' paper, Fp = 2.0*M*FA/(kPionMass2-q2), but Llewelyn-
956  // Smith uses Fp = 2.0*M^2*FA/(kPionMass2-q2), so I divide by M
957  // This gives units of GeV^-1
958  double Fp = -1.0/M*fFormFactors.Fp();
959 
960 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
961  LOG("Nieves", pDEBUG) << "\n" << fFormFactors;
962 #endif
963 
964  // Calculate auxiliary parameters
965  double M2 = TMath::Power(M, 2);
966  double FA2 = TMath::Power(FA, 2);
967  double Fp2 = TMath::Power(Fp, 2);
968  double F1V2 = TMath::Power(F1V, 2);
969  double xiF2V2 = TMath::Power(xiF2V, 2);
970  double q02 = TMath::Power(q[0], 2);
971  double dq2 = TMath::Power(dq, 2);
972  double q4 = TMath::Power(q2, 2);
973 
974  double t0,r00;
975  double CN=1.,CT=1.,CL=1.,imU=0;
976  CNCTCLimUcalc(qTildeP4, M, r, is_neutrino, tgtIsNucleus,
977  tgt_pdgc, A, Z, N, hitNucIsProton, CN, CT, CL, imU,
978  t0, r00, assumeFreeNucleon);
979 
980  if ( imU > kASmallNum )
981  return 0.;
982 
983  if ( !fRPA || assumeFreeNucleon ) {
984  CN = 1.0;
985  CT = 1.0;
986  CL = 1.0;
987  }
988 
989  double tulin[4] = {0.,0.,0.,0.};
990  double rulin[4][4] = { {0.,0.,0.,0.},
991  {0.,0.,0.,0.},
992  {0.,0.,0.,0.},
993  {0.,0.,0.,0.} };
994 
995  // TESTING CODE:
996  if(fCompareNievesTensors){
997  // Use average values for initial momentum to calculate A, as given
998  // in Appendix B of Nieves' paper. T gives average values of components
999  // of p, and R gives the average value of two components multiplied
1000  // together
1001  double t3 = (0.5*q2 + q0*t0)/dq; // Average pz
1002 
1003  // Vector of p
1004 
1005  tulin[0] = t0;
1006  tulin[3] = t3;
1007 
1008  // R is a 4x4 matrix, with R[mu][nu] is the average
1009  // value of p[mu]*p[nu]
1010  double aR = r00-M2;
1011  double bR = (q4+4.0*r00*q02+4.0*q2*q0*t0)/(4.0*dq2);
1012  double gamma = (aR-bR)/2.0;
1013  double delta = (-aR+3.0*bR)/2.0/dq2;
1014 
1015  double r03 = (0.5*q2*t0 + q0*r00)/dq; // Average E(p)*pz
1016 
1017  rulin[0][0] = r00;
1018  rulin[0][3] = r03;
1019  rulin[1][1] = gamma;
1020  rulin[2][2] = gamma;
1021  rulin[3][0] = r03;
1022  rulin[3][3] = gamma+delta*dq2; // END TESTING CODE
1023  }
1024  else {
1025  // For normal code execulation, tulin is the initial nucleon momentum
1026  tulin[0] = inNucleonMomOnShell.E();
1027  tulin[1] = inNucleonMomOnShell.Px();
1028  tulin[2] = inNucleonMomOnShell.Py();
1029  tulin[3] = inNucleonMomOnShell.Pz();
1030 
1031  for(int i=0; i<4; i++)
1032  for(int j=0; j<4; j++)
1033  rulin[i][j] = tulin[i]*tulin[j];
1034  }
1035 
1036  //Additional constants and variables
1037  const int g[4][4] = {{1,0,0,0},{0,-1,0,0},{0,0,-1,0},{0,0,0,-1}};
1038  const std::complex<double> iNum(0,1);
1039  int leviCivitaIndexArray[4];
1040  double imaginaryPart = 0;
1041 
1042  std::complex<double> sum(0.0,0.0);
1043 
1044  double kPrimek = k[0]*kPrime[0]-k[1]*kPrime[1]-k[2]*kPrime[2]-k[3]*kPrime[3];
1045 
1046  std::complex<double> Lmunu(0.0,0.0),Lnumu(0.0,0.0);
1047  std::complex<double> Amunu(0.0,0.0),Anumu(0.0,0.0);
1048 
1049  // Calculate LmunuAnumu by iterating over mu and nu
1050  // In each iteration, add LmunuAnumu to sum if mu=nu, and add
1051  // LmunuAnumu + LnumuAmunu if mu != nu, since we start nu at mu
1052  double axx=0.,azz=0.,a0z=0.,a00=0.,axy=0.;
1053  for(int mu=0;mu<4;mu++){
1054  for(int nu=mu;nu<4;nu++){
1055  imaginaryPart = 0;
1056  if(mu == nu){
1057  //if mu==nu then levi-civita = 0, so imaginary part = 0
1058  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;
1059  }else{
1060  //if mu!=nu, then g[mu][nu] = 0
1061  //This same leviCivitaIndex array can be used in the else portion when
1062  //calculating Anumu
1063  leviCivitaIndexArray[0] = mu;
1064  leviCivitaIndexArray[1] = nu;
1065  for(int a=0;a<4;a++){
1066  for(int b=0;b<4;b++){
1067  leviCivitaIndexArray[2] = a;
1068  leviCivitaIndexArray[3] = b;
1069  imaginaryPart += - leviCivita(leviCivitaIndexArray)*kPrime[a]*k[b];
1070  }
1071  }
1072  //real(Lmunu) is symmetric, and imag(Lmunu) is antisymmetric
1073  //std::complex<double> num(g[mu][mu]*kPrime[mu]*g[nu][nu]*k[nu]+g[nu][nu]*kPrime[nu]*g[mu][mu]*k[mu],imaginaryPart);
1074  Lmunu = g[mu][mu]*kPrime[mu]*g[nu][nu]*k[nu]+g[nu][nu]*kPrime[nu]*g[mu][mu]*k[mu] + iNum*imaginaryPart;
1075  Lnumu = g[nu][nu]*kPrime[nu]*g[mu][mu]*k[mu]+g[mu][mu]*kPrime[mu]*g[nu][nu]*k[nu ]- iNum*imaginaryPart;
1076  } // End Lmunu calculation
1077 
1078  if(mu ==0 && nu == 0){
1079  Amunu = 16.0*F1V2*(2.0*rulin[0][0]*CN+2.0*q[0]*tulin[0]+q2/2.0)+
1080  2.0*q2*xiF2V2*
1081  (4.0-4.0*rulin[0][0]/M2-4.0*q[0]*tulin[0]/M2-q02*(4.0/q2+1.0/M2)) +
1082  4.0*FA2*(2.0*rulin[0][0]+2.0*q[0]*tulin[0]+(q2/2.0-2.0*M2))-
1083  (2.0*CL*Fp2*q2+8.0*FA*Fp*CL*M)*q02-16.0*F1V*xiF2V*(-q2+q02)*CN;
1084  a00 = real(Amunu); // TESTING CODE
1085  sum += Lmunu*Amunu;
1086  }else if(mu == 0 && nu == 3){
1087  Amunu = 16.0*F1V2*((2.0*rulin[0][3]+tulin[0]*dq)*CN+tulin[3]*q[0])+
1088  2.0*q2*xiF2V2*
1089  (-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))+
1090  4.0*FA2*((2.0*rulin[0][3]+dq*tulin[0])*CL+q[0]*tulin[3])-
1091  (2.0*CL*Fp2*q2+8.0*FA*Fp*CL*M)*dq*q[0]-
1092  16.0*F1V*xiF2V*dq*q[0];
1093  a0z= real(Amunu); // TESTING CODE
1094  Anumu = Amunu;
1095  sum += Lmunu*Anumu + Lnumu*Amunu;
1096  }else if(mu == 3 && nu == 3){
1097  Amunu = 16.0*F1V2*(2.0*rulin[3][3]+2.0*dq*tulin[3]-q2/2.0)+
1098  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))+
1099  4.0*FA2*(2.0*rulin[3][3]+2.0*dq*tulin[3]-(q2/2.0-2.0*CL*M2))-
1100  (2.0*CL*Fp2*q2+8.0*FA*Fp*CL*M)*dq2-
1101  16.0*F1V*xiF2V*(q2+dq2);
1102  azz = real(Amunu); // TESTING CODE
1103  sum += Lmunu*Amunu;
1104  }else if(mu ==1 && nu == 1){
1105  Amunu = 16.0*F1V2*(2.0*rulin[1][1]-q2/2.0)+
1106  2.0*q2*xiF2V2*(-4.0*CT-4.0*rulin[1][1]/M2) +
1107  4.0*FA2*(2.0*rulin[1][1]-(q2/2.0-2.0*CT*M2))-
1108  16.0*F1V*xiF2V*CT*q2;
1109  axx = real(Amunu); // TESTING CODE
1110  sum += Lmunu*Amunu;
1111  }else if(mu == 2 && nu == 2){
1112  // Ayy not explicitly listed in paper. This is included so rotating the
1113  // coordinates of k and k' about the z-axis does not change the xsec.
1114  Amunu = 16.0*F1V2*(2.0*rulin[2][2]-q2/2.0)+
1115  2.0*q2*xiF2V2*(-4.0*CT-4.0*rulin[2][2]/M2) +
1116  4.0*FA2*(2.0*rulin[2][2]-(q2/2.0-2.0*CT*M2))-
1117  16.0*F1V*xiF2V*CT*q2;
1118  sum += Lmunu*Amunu;
1119  }else if(mu ==1 && nu == 2){
1120  Amunu = sign*16.0*iNum*FA*(xiF2V+F1V)*(-dq*tulin[0]*CT + q[0]*tulin[3]);
1121  Anumu = -Amunu; // Im(A) is antisymmetric
1122  axy = imag(Amunu); // TESTING CODE
1123  sum += Lmunu*Anumu+Lnumu*Amunu;
1124  }
1125  // All other terms will be 0 because the initial nucleus is at rest and
1126  // qTilde is in the z direction
1127 
1128  } // End loop over nu
1129  } // End loop over mu
1130 
1131  // TESTING CODE
1132  if(fCompareNievesTensors){
1133  // get tmu
1134  double tmugev = leptonMom.E() - leptonMom.Mag();
1135  // Print Q2, form factors, and tensor elts
1136  std::ofstream ffstream;
1137  ffstream.open(fTensorsOutFile, std::ios_base::app);
1138  if(q0 > 0){
1139  ffstream << -q2 << "\t" << q[0] << "\t" << dq
1140  << "\t" << axx << "\t" << azz << "\t" << a0z
1141  << "\t" << a00 << "\t" << axy << "\t"
1142  << CT << "\t" << CL << "\t" << CN << "\t"
1143  << tmugev << "\t" << imU << "\t"
1144  << F1V << "\t" << xiF2V << "\t"
1145  << FA << "\t" << Fp << "\t"
1146  << tulin[0] << "\t"<< tulin[1] << "\t"
1147  << tulin[2] << "\t"<< tulin[3] << "\t"
1148  << rulin[0][0]<< "\t"<< rulin[0][1]<< "\t"
1149  << rulin[0][2]<< "\t"<< rulin[0][3]<< "\t"
1150  << rulin[1][0]<< "\t"<< rulin[1][1]<< "\t"
1151  << rulin[1][2]<< "\t"<< rulin[1][3]<< "\t"
1152  << rulin[2][0]<< "\t"<< rulin[2][1]<< "\t"
1153  << rulin[2][2]<< "\t"<< rulin[2][3]<< "\t"
1154  << rulin[3][0]<< "\t"<< rulin[3][1]<< "\t"
1155  << rulin[3][2]<< "\t"<< rulin[3][3]<< "\t"
1156  << fVc << "\t" << fCoulombFactor << "\t";
1157  ffstream << "\n";
1158  }
1159  ffstream.close();
1160  }
1161  // END TESTING CODE
1162 
1163  // Since the real parts of A and L are both symmetric and the imaginary
1164  // parts are antisymmetric, the contraction should be real
1165  if ( imag(sum) > kASmallNum )
1166  LOG("Nieves",pWARN) << "Imaginary part of tensor contraction is nonzero "
1167  << "in QEL XSec, real(sum) = " << real(sum)
1168  << "imag(sum) = " << imag(sum);
1169 
1170  return real(sum);
1171 }
genie::NievesQELCCPXSec::fRPA
bool fRPA
use RPA corrections
Definition: NievesQELCCPXSec.h:81
genie::units::g
static constexpr double g
Definition: Units.h:144
genie::Target::HitNucPdg
int HitNucPdg(void) const
Definition: Target.cxx:304
genie::Target::A
int A(void) const
Definition: Target.h:70
genie::Target::IsNucleus
bool IsNucleus(void) const
Definition: Target.cxx:272
genie::NievesQELCCPXSec::CNCTCLimUcalc
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
Definition: NievesQELCCPXSec.cxx:519
genie::Target::Pdg
int Pdg(void) const
Definition: Target.h:71
genie::NievesQELCCPXSec::fFormFactors
QELFormFactors fFormFactors
Definition: NievesQELCCPXSec.h:68
genie::units::b
static constexpr double b
Definition: Units.h:78
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::units::A
static constexpr double A
Definition: Units.h:74
genie::NievesQELCCPXSec::fCompareNievesTensors
bool fCompareNievesTensors
print tensors
Definition: NievesQELCCPXSec.h:156
a
const double a
Definition: gUpMuFluxGen.cxx:164
genie::Target::Z
int Z(void) const
Definition: Target.h:68
genie::controls::kASmallNum
static const double kASmallNum
Definition: Controls.h:40
genie::QELFormFactors::xiF2V
double xiF2V(void) const
Get the computed form factor xi*F2V.
Definition: QELFormFactors.h:56
pWARN
#define pWARN
Definition: Messenger.h:60
genie::NievesQELCCPXSec::fTensorsOutFile
TString fTensorsOutFile
file to print tensors to
Definition: NievesQELCCPXSec.h:157
genie::Target::N
int N(void) const
Definition: Target.h:69
genie::Target::HitNucPosition
double HitNucPosition(void) const
Definition: Target.h:89
genie::NievesQELCCPXSec::leviCivita
int leviCivita(int input[]) const
Definition: NievesQELCCPXSec.cxx:891
genie::QELFormFactors::F1V
double F1V(void) const
Get the computed form factor F1V.
Definition: QELFormFactors.h:53
genie::QELFormFactors::Fp
double Fp(void) const
Get the computed form factor Fp.
Definition: QELFormFactors.h:62
genie::QELFormFactors::FA
double FA(void) const
Get the computed form factor FA.
Definition: QELFormFactors.h:59
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
void NievesQELCCPXSec::LoadConfig ( void  )
private

Definition at line 400 of file NievesQELCCPXSec.cxx.

References fCompareNievesTensors, fCos8c2, fCoulomb, fCoulombRmaxMode, fCoulombScale, fDoPauliBlocking, fEnergyCutOff, genie::units::fermi, fFormFactors, fFormFactorsModel, fhbarc, fIntegralNucleonBindingMode, fKFTable, fKFTableName, fLFG, fNuclModel, fPauliBlocker, fQvalueShifter, fR0, fRPA, fXSecCCScale, fXSecIntegrator, fXSecNCScale, genie::gAbortingInErr, genie::Algorithm::GetConfig(), genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), genie::FermiMomentumTablePool::GetTable(), genie::Algorithm::Id(), genie::FermiMomentumTablePool::Instance(), genie::constants::kLightSpeed, genie::kMatchNieves, genie::kMatchVertexGeneratorRmax, genie::kNucmLocalFermiGas, genie::constants::kPlankConstant, LOG, genie::NuclearModelI::ModelType(), pERROR, pFATAL, pNOTICE, genie::QELFormFactors::SetModel(), genie::utils::StringToQELBindingMode(), and genie::Algorithm::SubAlg().

Referenced by Configure().

401 {
402  bool good_config = true ;
403  double thc;
404  GetParam( "CabibboAngle", thc ) ;
405  fCos8c2 = TMath::Power(TMath::Cos(thc), 2);
406 
407  // Cross section scaling factor
408  GetParam( "QEL-CC-XSecScale", fXSecCCScale ) ;
409  GetParam( "QEL-NC-XSecScale", fXSecNCScale ) ;
410 
411  // hbarc for unit conversion, GeV*fm
412  fhbarc = kLightSpeed*kPlankConstant/genie::units::fermi;
413 
414  // load QEL form factors model
415  fFormFactorsModel = dynamic_cast<const QELFormFactorsModelI *> (
416  this->SubAlg("FormFactorsAlg") );
417  assert(fFormFactorsModel);
418  fFormFactors.SetModel( fFormFactorsModel ); // <-- attach algorithm
419 
420  // load XSec Integrator
421  fXSecIntegrator = dynamic_cast<const XSecIntegratorI*>(
422  this->SubAlg("XSec-Integrator") );
423  assert(fXSecIntegrator);
424 
425  // Load settings for RPA and Coulomb effects
426 
427  // RPA corrections will not affect a free nucleon
428  GetParamDef("RPA", fRPA, true ) ;
429 
430  // Coulomb Correction- adds a correction factor, and alters outgoing lepton
431  // 3-momentum magnitude (but not direction)
432  // Correction only becomes sizeable near threshold and/or for heavy nuclei
433  GetParamDef( "Coulomb", fCoulomb, true ) ;
434 
435  LOG("Nieves", pNOTICE) << "RPA=" << fRPA << ", useCoulomb=" << fCoulomb;
436 
437  // Get nuclear model for use in Integral()
438  RgKey nuclkey = "IntegralNuclearModel";
439  fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
440  assert(fNuclModel);
441 
442  // Check if the model is a local Fermi gas
443  fLFG = fNuclModel->ModelType(Target()) == kNucmLocalFermiGas;
444 
445  if(!fLFG){
446  // get the Fermi momentum table for relativistic Fermi gas
447  GetParam( "FermiMomentumTable", fKFTableName ) ;
448 
449  fKFTable = 0;
450  FermiMomentumTablePool * kftp = FermiMomentumTablePool::Instance();
451  fKFTable = kftp->GetTable( fKFTableName );
452  assert( fKFTable );
453  }
454 
455  // TESTING CODE
456  GetParamDef( "PrintDebugData", fCompareNievesTensors, false ) ;
457  // END TESTING CODE
458 
459  // Nuclear radius parameter (R = R0*A^(1/3)) to use when computing
460  // the maximum radius to use to integrate the Coulomb potential
461  GetParam("NUCL-R0", fR0) ; // fm
462 
463  std::string temp_mode;
464  GetParamDef( "RmaxMode", temp_mode, std::string("VertexGenerator") ) ;
465 
466  // Translate the string setting the Rmax mode to the appropriate
467  // enum value, or complain if one couldn't be found
468  if ( temp_mode == "VertexGenerator" ) {
469  fCoulombRmaxMode = kMatchVertexGeneratorRmax;
470  }
471  else if ( temp_mode == "Nieves" ) {
472  fCoulombRmaxMode = kMatchNieves;
473  }
474  else {
475  LOG("Nieves", pFATAL) << "Unrecognized setting \"" << temp_mode
476  << "\" requested for the RmaxMode parameter in the"
477  << " configuration for NievesQELCCPXSec";
478  gAbortingInErr = true;
479  std::exit(1);
480  }
481 
482  // Method to use to calculate the binding energy of the initial hit nucleon when
483  // generating splines
484  std::string temp_binding_mode;
485  GetParamDef( "IntegralNucleonBindingMode", temp_binding_mode, std::string("UseNuclearModel") );
486  fIntegralNucleonBindingMode = genie::utils::StringToQELBindingMode( temp_binding_mode );
487 
488  // Cutoff energy for integrating over nucleon momentum distribution (above this
489  // lab-frame probe energy, the effects of Fermi motion and binding energy
490  // are taken to be negligible for computing the total cross section)
491  GetParamDef("IntegralNuclearInfluenceCutoffEnergy", fEnergyCutOff, 2.5 ) ;
492 
493  // Get PauliBlocker for possible use in XSec()
494  fPauliBlocker = dynamic_cast<const PauliBlocker*>( this->SubAlg("PauliBlockerAlg") );
495  assert( fPauliBlocker );
496 
497  // Decide whether or not it should be used in XSec()
498  GetParamDef( "DoPauliBlocking", fDoPauliBlocking, true );
499 
500  // Read optional QvalueShifter:
501  fQvalueShifter = nullptr;
502  if( GetConfig().Exists("QvalueShifterAlg") ) {
503  fQvalueShifter = dynamic_cast<const QvalueShifter *> ( this->SubAlg("QvalueShifterAlg") );
504  if( !fQvalueShifter ) {
505  good_config = false ;
506  LOG("NievesQELCCPXSec", pERROR) << "The required QvalueShifterAlg does not exist. AlgID is : " << SubAlg("QvalueShifterAlg")->Id() ;
507  }
508  }
509 
510  if( ! good_config ) {
511  LOG("NievesQELCCPXSec", pERROR) << "Configuration has failed.";
512  exit(78) ;
513  }
514 
515  // Scaling factor for the Coulomb potential
516  GetParamDef( "CoulombScale", fCoulombScale, 1.0 );
517 }
genie::NievesQELCCPXSec::fQvalueShifter
const QvalueShifter * fQvalueShifter
Optional algorithm to retrieve the qvalue shift for a given target.
Definition: NievesQELCCPXSec.h:76
genie::NievesQELCCPXSec::fRPA
bool fRPA
use RPA corrections
Definition: NievesQELCCPXSec.h:81
genie::NievesQELCCPXSec::fXSecCCScale
double fXSecCCScale
external xsec scaling factor for CC
Definition: NievesQELCCPXSec.h:73
pERROR
#define pERROR
Definition: Messenger.h:59
genie::XSecIntegratorI
Cross Section Integrator Interface.
Definition: XSecIntegratorI.h:29
genie::NievesQELCCPXSec::fIntegralNucleonBindingMode
QELEvGen_BindingMode_t fIntegralNucleonBindingMode
Definition: NievesQELCCPXSec.h:93
genie::NievesQELCCPXSec::fKFTable
const FermiMomentumTable * fKFTable
Definition: NievesQELCCPXSec.h:88
genie::QELFormFactors::SetModel
void SetModel(const QELFormFactorsModelI *model)
Attach an algorithm.
Definition: QELFormFactors.cxx:43
pFATAL
#define pFATAL
Definition: Messenger.h:56
genie::FermiMomentumTablePool::Instance
static FermiMomentumTablePool * Instance(void)
Definition: FermiMomentumTablePool.cxx:62
genie::PauliBlocker
Examines whether the generated event should be Pauli blocked. Is a concerete implementation of the Ev...
Definition: PauliBlocker.h:36
genie::NievesQELCCPXSec::fXSecIntegrator
const XSecIntegratorI * fXSecIntegrator
Definition: NievesQELCCPXSec.h:70
genie::NuclearModelI
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
genie::NievesQELCCPXSec::fFormFactors
QELFormFactors fFormFactors
Definition: NievesQELCCPXSec.h:68
genie::NievesQELCCPXSec::fCoulombScale
double fCoulombScale
Scaling factor for the Coulomb potential.
Definition: NievesQELCCPXSec.h:110
genie::NievesQELCCPXSec::fCoulomb
bool fCoulomb
use Coulomb corrections
Definition: NievesQELCCPXSec.h:82
genie::NievesQELCCPXSec::fCoulombRmaxMode
Nieves_Coulomb_Rmax_t fCoulombRmaxMode
Definition: NievesQELCCPXSec.h:114
genie::NuclearModelI::ModelType
virtual NuclearModel_t ModelType(const Target &) const =0
genie::Algorithm::GetConfig
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:246
genie::NievesQELCCPXSec::fFormFactorsModel
const QELFormFactorsModelI * fFormFactorsModel
Definition: NievesQELCCPXSec.h:69
genie::NievesQELCCPXSec::fhbarc
double fhbarc
hbar*c in GeV*fm
Definition: NievesQELCCPXSec.h:78
genie::NievesQELCCPXSec::fCos8c2
double fCos8c2
cos^2(cabibbo angle)
Definition: NievesQELCCPXSec.h:71
genie::NievesQELCCPXSec::fNuclModel
const NuclearModelI * fNuclModel
Nuclear Model for integration.
Definition: NievesQELCCPXSec.h:84
genie::FermiMomentumTablePool
Singleton class to load &amp; serve tables of Fermi momentum constants.
Definition: FermiMomentumTablePool.h:34
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::QELFormFactorsModelI
Pure abstract base class. Defines the QELFormFactorsModelI interface to be implemented by any algorit...
Definition: QELFormFactorsModelI.h:28
genie::FermiMomentumTablePool::GetTable
const FermiMomentumTable * GetTable(string name)
Definition: FermiMomentumTablePool.cxx:76
genie::NievesQELCCPXSec::fCompareNievesTensors
bool fCompareNievesTensors
print tensors
Definition: NievesQELCCPXSec.h:156
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
genie::NievesQELCCPXSec::fDoPauliBlocking
bool fDoPauliBlocking
Whether to apply Pauli blocking in XSec()
Definition: NievesQELCCPXSec.h:100
genie::NievesQELCCPXSec::fXSecNCScale
double fXSecNCScale
external xsec scaling factor for NC
Definition: NievesQELCCPXSec.h:74
genie::Algorithm::Id
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
RgKey
string RgKey
Definition: RegistryItemTypeDef.h:33
genie::utils::StringToQELBindingMode
QELEvGen_BindingMode_t StringToQELBindingMode(const std::string &mode_str)
Definition: QELUtils.cxx:194
genie::units::fermi
static constexpr double fermi
Definition: Units.h:55
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::Algorithm::GetParamDef
bool GetParamDef(const RgKey &name, T &p, const T &def) const
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
genie::gAbortingInErr
bool gAbortingInErr
Definition: Messenger.cxx:34
genie::NievesQELCCPXSec::fPauliBlocker
const genie::PauliBlocker * fPauliBlocker
The PauliBlocker instance to use to apply that correction.
Definition: NievesQELCCPXSec.h:102
genie::constants::kPlankConstant
static const double kPlankConstant
Definition: Framework/Conventions/Constants.h:32
genie::Algorithm::SubAlg
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
std::complex< double > NievesQELCCPXSec::relLindhard ( double  q0gev,
double  dqgev,
double  kFgev,
double  M,
bool  isNeutrino,
std::complex< double >  relLindIm 
) const
private

Definition at line 662 of file NievesQELCCPXSec.cxx.

References fhbarc, LOG, genie::units::m, pWARN, relLindhardIm(), and ruLinRelX().

Referenced by CNCTCLimUcalc().

666 {
667  double q0 = q0gev/fhbarc;
668  double qm = dqgev/fhbarc;
669  double kf = kFgev/fhbarc;
670  double m = M/fhbarc;
671 
672  if(q0>qm){
673  LOG("Nieves", pWARN) << "relLindhard() failed";
674  return 0.0;
675  }
676 
677  std::complex<double> RealLinRel(ruLinRelX(q0,qm,kf,m)+ruLinRelX(-q0,qm,kf,m));
678  double t0,r00;
679  std::complex<double> ImLinRel(relLindhardIm(q0gev,dqgev,kFgev,kFgev,M,isNeutrino,t0,r00));
680  //Units of GeV^2
681  return(RealLinRel*TMath::Power(fhbarc,2) + 2.0*ImLinRel);
682 }
genie::NievesQELCCPXSec::ruLinRelX
std::complex< double > ruLinRelX(double q0, double qm, double kf, double m) const
Definition: NievesQELCCPXSec.cxx:685
genie::NievesQELCCPXSec::fhbarc
double fhbarc
hbar*c in GeV*fm
Definition: NievesQELCCPXSec.h:78
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
pWARN
#define pWARN
Definition: Messenger.h:60
genie::units::m
static constexpr double m
Definition: Units.h:71
genie::NievesQELCCPXSec::relLindhardIm
std::complex< double > relLindhardIm(double q0gev, double dqgev, double kFngev, double kFpgev, double M, bool isNeutrino, double &t0, double &r00) const
Definition: NievesQELCCPXSec.cxx:605
std::complex< double > NievesQELCCPXSec::relLindhardIm ( double  q0gev,
double  dqgev,
double  kFngev,
double  kFpgev,
double  M,
bool  isNeutrino,
double &  t0,
double &  r00 
) const
private

Definition at line 605 of file NievesQELCCPXSec.cxx.

References a, and genie::constants::kPi.

Referenced by CNCTCLimUcalc(), and relLindhard().

611 {
612  double M2 = TMath::Power(M,2);
613  double EF1,EF2;
614  if(isNeutrino){
615  EF1 = TMath::Sqrt(M2+TMath::Power(kFn,2)); //EFn
616  EF2 = TMath::Sqrt(M2+TMath::Power(kFp,2)); //EFp
617  }else{
618  EF1 = TMath::Sqrt(M2+TMath::Power(kFp,2)); //EFp
619  EF2 = TMath::Sqrt(M2+TMath::Power(kFn,2)); //EFn
620  }
621 
622  double q2 = TMath::Power(q0,2) - TMath::Power(dq,2);
623  double a = (-q0+dq*TMath::Sqrt(1-4.0*M2/q2))/2.0;
624  double epsRP = TMath::Max(TMath::Max(M,EF2-q0),a);
625 
626  // Other theta functions for q are handled by nuclear suppression
627  // That is, q0>0 and -q2>0 are always handled, and q0>EF2-EF1 is
628  // handled if pauli blocking is on, because otherwise the final
629  // nucleon would be below the fermi sea
630  //if(fNievesSuppression && !interaction->TestBit(kIAssumeFreeNucleon )
631  //&& !EF1-epsRP<0){
632  //LOG("Nieves", pINFO) << "Average value of E(p) above Fermi sea";
633  //return 0;
634  //}else{
635  t0 = 0.5*(EF1+epsRP);
636  r00 = (TMath::Power(EF1,2)+TMath::Power(epsRP,2)+EF1*epsRP)/3.0;
637  std::complex<double> result(0.0,-M2/2.0/kPi/dq*(EF1-epsRP));
638  return result;
639  //}
640 }
a
const double a
Definition: gUpMuFluxGen.cxx:164
std::complex< double > NievesQELCCPXSec::ruLinRelX ( double  q0,
double  qm,
double  kf,
double  m 
) const
private

Definition at line 685 of file NievesQELCCPXSec.cxx.

References genie::constants::kPi2, and genie::units::m2.

Referenced by relLindhard().

687 {
688  double q02 = TMath::Power(q0, 2);
689  double qm2 = TMath::Power(qm, 2);
690  double kf2 = TMath::Power(kf, 2);
691  double m2 = TMath::Power(m, 2);
692  double m4 = TMath::Power(m, 4);
693 
694  double ef = TMath::Sqrt(m2+kf2);
695  double q2 = q02-qm2;
696  double q4 = TMath::Power(q2,2);
697  double ds = TMath::Sqrt(1.0-4.0*m2/q2);
698  double L1 = log((kf+ef)/m);
699  std::complex<double> uL2(
700  TMath::Log(TMath::Abs(
701  (ef + q0 - TMath::Sqrt(m2+TMath::Power(kf-qm,2)))/
702  (ef + q0 - TMath::Sqrt(m2 + TMath::Power(kf + qm,2))))) +
703  TMath::Log(TMath::Abs(
704  (ef + q0 + TMath::Sqrt(m2 + TMath::Power(kf - qm,2)))/
705  (ef + q0 + TMath::Sqrt(m2 + TMath::Power(kf + qm,2))))));
706 
707  std::complex<double> uL3(
708  TMath::Log(TMath::Abs((TMath::Power(2*kf + q0*ds,2)-qm2)/
709  (TMath::Power(2*kf - q0*ds,2)-qm2))) +
710  TMath::Log(TMath::Abs((TMath::Power(kf-ef*ds,2) - (4*m4*qm2)/q4)/
711  (TMath::Power(kf+ef*ds,2) - (4*m4*qm2)/q4))));
712 
713  std::complex<double> RlinrelX(-L1/(16.0*kPi2)+
714  uL2*(2.0*ef+q0)/(32.0*kPi2*qm)-
715  uL3*ds/(64.0*kPi2));
716 
717  return RlinrelX*16.0*m2;
718 }
genie::units::m2
static constexpr double m2
Definition: Units.h:72
genie::units::m
static constexpr double m
Definition: Units.h:71
bool NievesQELCCPXSec::ValidProcess ( const Interaction *  i) const
virtual

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 365 of file NievesQELCCPXSec.cxx.

References genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::pdg::IsAntiNeutrino(), genie::pdg::IsNeutrino(), genie::pdg::IsNeutron(), genie::pdg::IsProton(), genie::ProcessInfo::IsQuasiElastic(), genie::ProcessInfo::IsWeakCC(), genie::kISkipProcessChk, genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), and genie::InitialState::Tgt().

Referenced by XSec().

366 {
367  if(interaction->TestBit(kISkipProcessChk)) return true;
368 
369  const InitialState & init_state = interaction->InitState();
370  const ProcessInfo & proc_info = interaction->ProcInfo();
371 
372  if(!proc_info.IsQuasiElastic()) return false;
373 
374  int nuc = init_state.Tgt().HitNucPdg();
375  int nu = init_state.ProbePdg();
376 
377  bool isP = pdg::IsProton(nuc);
378  bool isN = pdg::IsNeutron(nuc);
379  bool isnu = pdg::IsNeutrino(nu);
380  bool isnub = pdg::IsAntiNeutrino(nu);
381 
382  bool prcok = proc_info.IsWeakCC() && ((isP&&isnub) || (isN&&isnu));
383  if(!prcok) return false;
384 
385  return true;
386 }
genie::ProcessInfo::IsWeakCC
bool IsWeakCC(void) const
Definition: ProcessInfo.cxx:203
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
genie::Target::HitNucPdg
int HitNucPdg(void) const
Definition: Target.cxx:304
genie::ProcessInfo::IsQuasiElastic
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
genie::pdg::IsNeutron
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
genie::ProcessInfo
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
genie::pdg::IsAntiNeutrino
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
genie::InitialState::ProbePdg
int ProbePdg(void) const
Definition: InitialState.h:64
genie::InitialState::Tgt
const Target & Tgt(void) const
Definition: InitialState.h:66
genie::kISkipProcessChk
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
genie::InitialState
Initial State information.
Definition: InitialState.h:48
double NievesQELCCPXSec::vcr ( const Target *  target,
double  r 
) const
private

Definition at line 832 of file NievesQELCCPXSec.cxx.

References genie::Target::A(), genie::units::A, fCoulombRmaxMode, fCoulombScale, fhbarc, fR0, func(), genie::gAbortingInErr, genie::utils::gsl::Integration1DimTypeFromString(), genie::Target::IsNucleus(), genie::constants::kAem, genie::kMatchNieves, genie::kMatchVertexGeneratorRmax, genie::constants::kPi, LOG, pFATAL, pNOTICE, and genie::Target::Z().

Referenced by XSec().

832  {
833  if(target->IsNucleus()){
834  int A = target->A();
835  int Z = target->Z();
836  double Rmax = 0.;
837 
838  if ( fCoulombRmaxMode == kMatchNieves ) {
839  // Rmax calculated using formula from Nieves' fortran code and default
840  // charge and neutron matter density parameters from NuclearUtils.cxx
841  if (A > 20) {
842  double c = TMath::Power(A,0.35), z = 0.54;
843  Rmax = c + 9.25*z;
844  }
845  else {
846  // c = 1.75 for A <= 20
847  Rmax = TMath::Sqrt(20.0)*1.75;
848  }
849  }
850  else if ( fCoulombRmaxMode == kMatchVertexGeneratorRmax ) {
851  // TODO: This solution is fragile. If the formula used by VertexGenerator
852  // changes, then this one will need to change too. Switch to using
853  // a common function to get Rmax for both.
854  Rmax = 3. * fR0 * std::pow(A, 1./3.);
855  }
856  else {
857  LOG("Nieves", pFATAL) << "Unrecognized setting for fCoulombRmaxMode encountered"
858  << " in NievesQELCCPXSec::vcr()";
859  gAbortingInErr = true;
860  std::exit(1);
861  }
862 
863  if(Rcurr >= Rmax){
864  LOG("Nieves",pNOTICE) << "Radius greater than maximum radius for coulomb corrections."
865  << " Integrating to max radius.";
866  Rcurr = Rmax;
867  }
868 
869  ROOT::Math::IBaseFunctionOneDim * func = new
870  utils::gsl::wrap::NievesQELvcrIntegrand(Rcurr,A,Z);
871  ROOT::Math::IntegrationOneDim::Type ig_type =
872  utils::gsl::Integration1DimTypeFromString("adaptive");
873 
874  double abstol = 1; // We mostly care about relative tolerance;
875  double reltol = 1E-4;
876  int nmaxeval = 100000;
877  ROOT::Math::Integrator ig(*func,ig_type,abstol,reltol,nmaxeval);
878  double result = ig.Integral(0,Rmax);
879  delete func;
880 
881  // Multiply by Z to normalize densities to number of protons
882  // Multiply by hbarc to put result in GeV instead of fm
883  // Multiply by an extra configurable scaling factor that defaults to unity
884  return -kAem*4*kPi*result*fhbarc*fCoulombScale;
885  }else{
886  // If target is not a nucleus the potential will be 0
887  return 0.0;
888  }
889 }
genie::utils::gsl::Integration1DimTypeFromString
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
Definition: GSLUtils.cxx:23
genie::Target::A
int A(void) const
Definition: Target.h:70
pFATAL
#define pFATAL
Definition: Messenger.h:56
genie::Target::IsNucleus
bool IsNucleus(void) const
Definition: Target.cxx:272
genie::NievesQELCCPXSec::fCoulombScale
double fCoulombScale
Scaling factor for the Coulomb potential.
Definition: NievesQELCCPXSec.h:110
genie::NievesQELCCPXSec::fCoulombRmaxMode
Nieves_Coulomb_Rmax_t fCoulombRmaxMode
Definition: NievesQELCCPXSec.h:114
genie::NievesQELCCPXSec::fhbarc
double fhbarc
hbar*c in GeV*fm
Definition: NievesQELCCPXSec.h:78
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::units::A
static constexpr double A
Definition: Units.h:74
func
double func(double x, double y)
Definition: gtestBLI2DUnifGrid.cxx:81
genie::Target::Z
int Z(void) const
Definition: Target.h:68
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::gAbortingInErr
bool gAbortingInErr
Definition: Messenger.cxx:34
double NievesQELCCPXSec::XSec ( const Interaction *  i,
KinePhaseSpace_t  k 
) const
virtual

Compute the cross section for the input interaction.

Implements genie::XSecAlgorithmI.

Definition at line 70 of file NievesQELCCPXSec.cxx.

References genie::KinePhaseSpace::AsString(), genie::QELFormFactors::Calculate(), genie::utils::EnergyDeltaFunctionSolutionQEL(), fCos8c2, fCoulomb, fDoPauliBlocking, fFormFactors, fPauliBlocker, fQvalueShifter, fRPA, genie::Kinematics::FSLeptonP4(), genie::Interaction::FSPrimLepton(), fXSecCCScale, fXSecNCScale, genie::PauliBlocker::GetFermiMomentum(), genie::InitialState::GetProbeP4(), genie::Kinematics::HadSystP4(), genie::Target::HitNucMass(), genie::Target::HitNucP4(), genie::Target::HitNucPdg(), genie::Target::HitNucPosition(), genie::pdg::IsNeutrino(), genie::Target::IsNucleus(), genie::pdg::IsProton(), genie::ProcessInfo::IsWeakCC(), genie::ProcessInfo::IsWeakNC(), genie::utils::mec::J(), genie::utils::kinematics::Jacobian(), genie::controls::kASmallNum, genie::constants::kGF2, genie::kIAssumeFreeNucleon, genie::Interaction::KinePtr(), genie::constants::kPi, genie::kPSQELEvGen, genie::kRfLab, LmunuAnumu(), LOG, genie::Target::N(), pDEBUG, genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), pWARN, genie::utils::kinematics::Q2(), genie::Interaction::RecoilNucleon(), genie::Interaction::RecoilNucleonPdg(), genie::Kinematics::SetQ2(), genie::QvalueShifter::Shift(), genie::InitialState::Tgt(), genie::XSecAlgorithmI::ValidKinematics(), ValidProcess(), vcr(), and genie::Target::Z().

72 {
73  /*// TESTING CODE:
74  // The first time this method is called, output tensor elements and other
75  // kinmeatics variables for various kinematics. This can the be compared
76  // to Nieves' fortran code for validation purposes
77  if(fCompareNievesTensors){
78  LOG("Nieves",pNOTICE) << "Printing tensor elements for specific "
79  << "kinematics for testing purposes";
80  CompareNievesTensors(interaction);
81  fCompareNievesTensors = false;
82  exit(0);
83  }
84  // END TESTING CODE*/
85 
86 
87  if ( !this->ValidProcess (interaction) ) return 0.;
88  if ( !this->ValidKinematics(interaction) ) return 0.;
89 
90  // Get kinematics and init-state parameters
91  const Kinematics & kinematics = interaction -> Kine();
92  const InitialState & init_state = interaction -> InitState();
93  const Target & target = init_state.Tgt();
94 
95  // HitNucMass() looks up the PDGLibrary (on-shell) value for the initial
96  // struck nucleon
97  double mNi = target.HitNucMass();
98 
99  // Hadronic matrix element for CC neutrino interactions should really use
100  // the "nucleon mass," i.e., the mean of the proton and neutrino masses.
101  // This expression would also work for NC and EM scattering (since the
102  // initial and final on-shell nucleon masses would be the same)
103  double mNucleon = ( mNi + interaction->RecoilNucleon()->Mass() ) / 2.;
104 
105  // Create a copy of the struck nucleon 4-momentum that is forced
106  // to be on-shell (this will be needed later for the tensor contraction,
107  // in which the nucleon is treated in this way)
108  double inNucleonOnShellEnergy = std::sqrt( std::pow(mNi, 2)
109  + std::pow(target.HitNucP4().P(), 2) );
110 
111  // The Nieves CCQE model follows the de Forest prescription: free nucleon
112  // (i.e., on-shell) form factors and spinors are used, but an effective
113  // value of the 4-momentum transfer "qTilde" is used when computing the
114  // contraction of the hadronic tensor. See comments in the
115  // FullDifferentialXSec() method of LwlynSmithQELCCPXSec for more details.
116  TLorentzVector inNucleonMomOnShell( target.HitNucP4().Vect(),
117  inNucleonOnShellEnergy );
118 
119  // Get the four kinematic vectors and caluclate GFactor
120  // Create copies of all kinematics, so they can be rotated
121  // and boosted to the nucleon rest frame (because the tensor
122  // constraction below only applies for the initial nucleon
123  // at rest and q in the z direction)
124 
125  // All 4-momenta should already be stored, with the hit nucleon off-shell
126  // as appropriate
127  TLorentzVector* tempNeutrino = init_state.GetProbeP4(kRfLab);
128  TLorentzVector neutrinoMom = *tempNeutrino;
129  delete tempNeutrino;
130  TLorentzVector inNucleonMom = target.HitNucP4();
131  TLorentzVector leptonMom = kinematics.FSLeptonP4();
132  TLorentzVector outNucleonMom = kinematics.HadSystP4();
133 
134  // Apply Pauli blocking if enabled
135  if ( fDoPauliBlocking && target.IsNucleus() && !interaction->TestBit(kIAssumeFreeNucleon) ) {
136  int final_nucleon_pdg = interaction->RecoilNucleonPdg();
137  double kF = fPauliBlocker->GetFermiMomentum(target, final_nucleon_pdg,
138  target.HitNucPosition());
139  double pNf = outNucleonMom.P();
140  if ( pNf < kF ) return 0.;
141  }
142 
143  // Use the lab kinematics to calculate the Gfactor, in order to make
144  // the XSec differential in initial nucleon momentum and energy
145  // Divide by 4.0 because Nieves' conventions for the leptonic and hadronic
146  // tensor contraction differ from LwlynSmith by a factor of 4
147  double Gfactor = kGF2*fCos8c2 / (8.0*kPi*kPi*inNucleonOnShellEnergy
148  *neutrinoMom.E()*outNucleonMom.E()*leptonMom.E()) / 4.0;
149 
150  // Calculate Coulomb corrections
151  double ml = interaction->FSPrimLepton()->Mass();
152  double ml2 = TMath::Power(ml, 2);
153  double coulombFactor = 1.0;
154  double plLocal = leptonMom.P();
155 
156  bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
157  double r = target.HitNucPosition();
158 
159  if ( fCoulomb ) {
160  // Coulomb potential
161  double Vc = vcr(& target, r);
162 
163  // Outgoing lepton energy and momentum including Coulomb potential
164  int sign = is_neutrino ? 1 : -1;
165  double El = leptonMom.E();
166  double pl = leptonMom.P();
167  double ElLocal = El - sign*Vc;
168 
169  if ( ElLocal - ml <= 0. ) {
170  LOG("Nieves", pDEBUG) << "Event should be rejected. Coulomb effects"
171  << " push kinematics below threshold. Returning xsec = 0.0";
172  return 0.0;
173  }
174 
175  // The Coulomb correction factor blows up as pl -> 0. To guard against
176  // unphysically huge corrections here, require that the lepton kinetic energy
177  // (at infinity) is larger than the magnitude of the Coulomb potential
178  // (should be around a few MeV)
179  double KEl = El - ml;
180  if ( KEl <= std::abs(Vc) ) {
181  LOG("Nieves", pDEBUG) << "Outgoing lepton has a very small kinetic energy."
182  << " Protecting against near-singularities in the Coulomb correction"
183  << " factor by returning xsec = 0.0";
184  return 0.0;
185  }
186 
187  // Local value of the lepton 3-momentum magnitude for the Coulomb
188  // correction
189  plLocal = TMath::Sqrt( ElLocal*ElLocal - ml2 );
190 
191  // Correction factor
192  coulombFactor= (plLocal * ElLocal) / (pl * El);
193 
194  }
195 
196  // When computing the contraction of the leptonic and hadronic tensors,
197  // we need to use an effective value of the 4-momentum transfer q.
198  // The energy transfer (q0) needs to be modified to account for the binding
199  // energy of the struck nucleon, while the 3-momentum transfer needs to
200  // be corrected for Coulomb effects.
201  //
202  // See the original Valencia model paper:
203  // https://journals.aps.org/prc/abstract/10.1103/PhysRevC.70.055503
204 
205  double q0Tilde = outNucleonMom.E() - inNucleonMomOnShell.E();
206 
207  // Shift the q0Tilde if required:
208  if( fQvalueShifter ) q0Tilde += q0Tilde * fQvalueShifter->Shift(*interaction) ;
209 
210  // If binding energy effects pull us into an unphysical region, return
211  // zero for the differential cross section
212  if ( q0Tilde <= 0. && target.IsNucleus() && !interaction->TestBit(kIAssumeFreeNucleon) ) return 0.;
213 
214  // Note that we're working in the lab frame (i.e., the rest frame
215  // of the target nucleus). We can therefore use Nieves' explicit
216  // form of the Amunu tensor if we rotate the 3-momenta so that
217  // qTilde is in the +z direction
218  TVector3 neutrinoMom3 = neutrinoMom.Vect();
219  TVector3 leptonMom3 = leptonMom.Vect();
220 
221  TVector3 inNucleonMom3 = inNucleonMom.Vect();
222  TVector3 outNucleonMom3 = outNucleonMom.Vect();
223 
224  // If Coulomb corrections are being used, adjust the lepton 3-momentum used
225  // to get q3VecTilde so that its magnitude matches the local
226  // Coulomb-corrected value calculated earlier. Note that, although the
227  // treatment of Coulomb corrections by Nieves et al. doesn't change the
228  // direction of the lepton 3-momentum, it *does* change the direction of the
229  // 3-momentum transfer, and so the correction should be applied *before*
230  // rotating coordinates into a frame where q3VecTilde lies along the positive
231  // z axis.
232  TVector3 leptonMomCoulomb3 = (! fCoulomb ) ? leptonMom3
233  : plLocal * leptonMom3 * (1. / leptonMom3.Mag());
234  TVector3 q3VecTilde = neutrinoMom3 - leptonMomCoulomb3;
235 
236  // Find the rotation angle needed to put q3VecTilde along z
237  TVector3 zvec(0.0, 0.0, 1.0);
238  TVector3 rot = ( q3VecTilde.Cross(zvec) ).Unit(); // Vector to rotate about
239  // Angle between the z direction and q
240  double angle = zvec.Angle( q3VecTilde );
241 
242  // Handle the edge case where q3VecTilde is along -z, so the
243  // cross product above vanishes
244  if ( q3VecTilde.Perp() == 0. && q3VecTilde.Z() < 0. ) {
245  rot = TVector3(0., 1., 0.);
246  angle = kPi;
247  }
248 
249  // Rotate if the rotation vector is not 0
250  if ( rot.Mag() >= kASmallNum ) {
251 
252  neutrinoMom3.Rotate(angle,rot);
253  neutrinoMom.SetVect(neutrinoMom3);
254 
255  leptonMom3.Rotate(angle,rot);
256  leptonMom.SetVect(leptonMom3);
257 
258  inNucleonMom3.Rotate(angle,rot);
259  inNucleonMom.SetVect(inNucleonMom3);
260  inNucleonMomOnShell.SetVect(inNucleonMom3);
261 
262  outNucleonMom3.Rotate(angle,rot);
263  outNucleonMom.SetVect(outNucleonMom3);
264 
265  }
266 
267  // Calculate q and qTilde
268  TLorentzVector qP4 = neutrinoMom - leptonMom;
269  TLorentzVector qTildeP4(0., 0., q3VecTilde.Mag(), q0Tilde);
270 
271  double Q2 = -1. * qP4.Mag2();
272  double Q2tilde = -1. * qTildeP4.Mag2();
273 
274  // Store Q2tilde in the interaction so that we get the correct
275  // values of the form factors (according to the de Forest prescription)
276  interaction->KinePtr()->SetQ2(Q2tilde);
277 
278  double q2 = -Q2tilde;
279 
280  // Check that q2 < 0 (accounting for rounding errors)
281  if ( q2 >= kASmallNum ) {
282  LOG("Nieves", pWARN) << "q2 >= 0, returning xsec = 0.0";
283  return 0.0;
284  }
285 
286  // Calculate form factors
287  fFormFactors.Calculate( interaction );
288 
289  // Now that the form factors have been calculated, store Q2
290  // in the event instead of Q2tilde
291  interaction->KinePtr()->SetQ2( Q2 );
292 
293  // Do the contraction of the leptonic and hadronic tensors. See the
294  // RPA-corrected expressions for the hadronic tensor elements in appendix A
295  // of Phys. Rev. C 70, 055503 (2004). Note that the on-shell 4-momentum of
296  // the initial struck nucleon should be used in the calculation, as well as
297  // the effective 4-momentum transfer q tilde (corrected for the nucleon
298  // binding energy and Coulomb effects)
299  double LmunuAnumuResult = LmunuAnumu(neutrinoMom, inNucleonMomOnShell,
300  leptonMom, qTildeP4, mNucleon, is_neutrino, target,
301  interaction->TestBit( kIAssumeFreeNucleon ));
302 
303  // Calculate xsec
304  double xsec = Gfactor*coulombFactor*LmunuAnumuResult;
305 
306  // Apply the factor that arises from elimination of the energy-conserving
307  // delta function
308  xsec *= genie::utils::EnergyDeltaFunctionSolutionQEL( *interaction );
309 
310  // Apply given scaling factor
311  double xsec_scale = 1 ;
312  const ProcessInfo& proc_info = interaction->ProcInfo();
313 
314  if( proc_info.IsWeakCC() ) xsec_scale = fXSecCCScale;
315  else if( proc_info.IsWeakNC() ) xsec_scale = fXSecNCScale;
316 
317  xsec *= xsec_scale ;
318 
319  LOG("Nieves",pDEBUG) << "TESTING: RPA=" << fRPA
320  << ", Coulomb=" << fCoulomb
321  << ", q2 = " << q2 << ", xsec = " << xsec;
322 
323  //----- The algorithm computes dxsec/dQ2 or kPSQELEvGen
324  // Check whether variable tranformation is needed
325  if ( kps != kPSQELEvGen ) {
326 
327  // Compute the appropriate Jacobian for transformation to the requested
328  // phase space
329  double J = utils::kinematics::Jacobian(interaction, kPSQELEvGen, kps);
330 
331 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
332  LOG("Nieves", pDEBUG)
333  << "Jacobian for transformation to: "
334  << KinePhaseSpace::AsString(kps) << ", J = " << J;
335 #endif
336  xsec *= J;
337  }
338 
339  // Number of scattering centers in the target
340  int nucpdgc = target.HitNucPdg();
341  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
342 
343  xsec *= NNucl; // nuclear xsec
344 
345  return xsec;
346 }
genie::NievesQELCCPXSec::fQvalueShifter
const QvalueShifter * fQvalueShifter
Optional algorithm to retrieve the qvalue shift for a given target.
Definition: NievesQELCCPXSec.h:76
genie::ProcessInfo::IsWeakCC
bool IsWeakCC(void) const
Definition: ProcessInfo.cxx:203
genie::NievesQELCCPXSec::fRPA
bool fRPA
use RPA corrections
Definition: NievesQELCCPXSec.h:81
genie::NievesQELCCPXSec::fXSecCCScale
double fXSecCCScale
external xsec scaling factor for CC
Definition: NievesQELCCPXSec.h:73
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
genie::utils::mec::J
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
genie::utils::kinematics::Q2
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
genie::Target::HitNucPdg
int HitNucPdg(void) const
Definition: Target.cxx:304
genie::Target::HitNucMass
double HitNucMass(void) const
Definition: Target.cxx:233
genie::Target::IsNucleus
bool IsNucleus(void) const
Definition: Target.cxx:272
genie::Kinematics
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
genie::Kinematics::HadSystP4
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:66
genie::NievesQELCCPXSec::fFormFactors
QELFormFactors fFormFactors
Definition: NievesQELCCPXSec.h:68
genie::NievesQELCCPXSec::fCoulomb
bool fCoulomb
use Coulomb corrections
Definition: NievesQELCCPXSec.h:82
genie::NievesQELCCPXSec::fCos8c2
double fCos8c2
cos^2(cabibbo angle)
Definition: NievesQELCCPXSec.h:71
genie::Target::HitNucP4
const TLorentzVector & HitNucP4(void) const
Definition: Target.h:91
genie::PauliBlocker::GetFermiMomentum
double GetFermiMomentum(const Target &tgt, int pdg_Nf, double radius=0.0) const
Get the Fermi momentum needed to check Pauli blocking.
Definition: PauliBlocker.cxx:168
genie::XSecAlgorithmI::ValidKinematics
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
Definition: XSecAlgorithmI.cxx:40
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
genie::ProcessInfo::IsWeakNC
bool IsWeakNC(void) const
Definition: ProcessInfo.cxx:208
genie::Kinematics::FSLeptonP4
const TLorentzVector & FSLeptonP4(void) const
Definition: Kinematics.h:65
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::utils::EnergyDeltaFunctionSolutionQEL
double EnergyDeltaFunctionSolutionQEL(const Interaction &inter)
Definition: QELUtils.cxx:50
genie::KinePhaseSpace::AsString
static string AsString(KinePhaseSpace_t kps)
Definition: KinePhaseSpace.h:85
genie::ProcessInfo
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
genie::InitialState::ProbePdg
int ProbePdg(void) const
Definition: InitialState.h:64
genie::NievesQELCCPXSec::fDoPauliBlocking
bool fDoPauliBlocking
Whether to apply Pauli blocking in XSec()
Definition: NievesQELCCPXSec.h:100
genie::NievesQELCCPXSec::fXSecNCScale
double fXSecNCScale
external xsec scaling factor for NC
Definition: NievesQELCCPXSec.h:74
genie::Target::Z
int Z(void) const
Definition: Target.h:68
genie::controls::kASmallNum
static const double kASmallNum
Definition: Controls.h:40
genie::QvalueShifter::Shift
virtual double Shift(const Target &target) const
Definition: QvalueShifter.cxx:51
pWARN
#define pWARN
Definition: Messenger.h:60
genie::QELFormFactors::Calculate
void Calculate(const Interaction *interaction)
Compute the form factors for the input interaction using the attached model.
Definition: QELFormFactors.cxx:49
genie::NievesQELCCPXSec::vcr
double vcr(const Target *target, double r) const
Definition: NievesQELCCPXSec.cxx:832
genie::Target::N
int N(void) const
Definition: Target.h:69
genie::kIAssumeFreeNucleon
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
genie::Target::HitNucPosition
double HitNucPosition(void) const
Definition: Target.h:89
genie::utils::kinematics::Jacobian
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
genie::NievesQELCCPXSec::ValidProcess
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Definition: NievesQELCCPXSec.cxx:365
genie::InitialState::Tgt
const Target & Tgt(void) const
Definition: InitialState.h:66
genie::NievesQELCCPXSec::LmunuAnumu
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
Definition: NievesQELCCPXSec.cxx:925
genie::NievesQELCCPXSec::fPauliBlocker
const genie::PauliBlocker * fPauliBlocker
The PauliBlocker instance to use to apply that correction.
Definition: NievesQELCCPXSec.h:102
genie::InitialState::GetProbeP4
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
Definition: InitialState.cxx:325
genie::InitialState
Initial State information.
Definition: InitialState.h:48
pDEBUG
#define pDEBUG
Definition: Messenger.h:63

Member Data Documentation

bool genie::NievesQELCCPXSec::fCompareNievesTensors
mutableprivate

print tensors

Definition at line 156 of file NievesQELCCPXSec.h.

Referenced by LmunuAnumu(), and LoadConfig().

double genie::NievesQELCCPXSec::fCos8c2
private

cos^2(cabibbo angle)

Definition at line 71 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

bool genie::NievesQELCCPXSec::fCoulomb
private

use Coulomb corrections

Definition at line 82 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::NievesQELCCPXSec::fCoulombFactor
mutableprivate

Definition at line 158 of file NievesQELCCPXSec.h.

Referenced by LmunuAnumu().

Nieves_Coulomb_Rmax_t genie::NievesQELCCPXSec::fCoulombRmaxMode
private

Enum variable describing which method of computing Rmax should be used for integrating the Coulomb potential

Definition at line 114 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and vcr().

double genie::NievesQELCCPXSec::fCoulombScale
private

Scaling factor for the Coulomb potential.

Definition at line 110 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and vcr().

bool genie::NievesQELCCPXSec::fDoPauliBlocking
private

Whether to apply Pauli blocking in XSec()

Definition at line 100 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::NievesQELCCPXSec::fEnergyCutOff
private

Cutoff lab-frame probe energy above which the effects of Fermi motion and binding energy are ignored when computing the total cross section

Definition at line 97 of file NievesQELCCPXSec.h.

Referenced by LoadConfig().

QELFormFactors genie::NievesQELCCPXSec::fFormFactors
mutableprivate

Definition at line 68 of file NievesQELCCPXSec.h.

Referenced by LmunuAnumu(), LoadConfig(), and XSec().

const QELFormFactorsModelI* genie::NievesQELCCPXSec::fFormFactorsModel
private

Definition at line 69 of file NievesQELCCPXSec.h.

Referenced by LoadConfig().

double genie::NievesQELCCPXSec::fhbarc
private

hbar*c in GeV*fm

Definition at line 78 of file NievesQELCCPXSec.h.

Referenced by CNCTCLimUcalc(), deltaLindhard(), LoadConfig(), relLindhard(), and vcr().

QELEvGen_BindingMode_t genie::NievesQELCCPXSec::fIntegralNucleonBindingMode
private

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

Definition at line 93 of file NievesQELCCPXSec.h.

Referenced by LoadConfig().

const FermiMomentumTable* genie::NievesQELCCPXSec::fKFTable
private

Definition at line 88 of file NievesQELCCPXSec.h.

Referenced by CNCTCLimUcalc(), and LoadConfig().

string genie::NievesQELCCPXSec::fKFTableName
private

Definition at line 89 of file NievesQELCCPXSec.h.

Referenced by LoadConfig().

bool genie::NievesQELCCPXSec::fLFG
private

Definition at line 87 of file NievesQELCCPXSec.h.

Referenced by CNCTCLimUcalc(), and LoadConfig().

const NuclearModelI* genie::NievesQELCCPXSec::fNuclModel
private

Nuclear Model for integration.

Definition at line 84 of file NievesQELCCPXSec.h.

Referenced by LoadConfig().

const genie::PauliBlocker* genie::NievesQELCCPXSec::fPauliBlocker
private

The PauliBlocker instance to use to apply that correction.

Definition at line 102 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

const QvalueShifter* genie::NievesQELCCPXSec::fQvalueShifter
private

Optional algorithm to retrieve the qvalue shift for a given target.

Definition at line 76 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::NievesQELCCPXSec::fR0
private

Nuclear radius parameter r = R0*A^(1/3) used to compute the maximum radius for integration of the Coulomb potential when matching the VertexGenerator method

Definition at line 107 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and vcr().

bool genie::NievesQELCCPXSec::fRPA
mutableprivate

use RPA corrections

Definition at line 81 of file NievesQELCCPXSec.h.

Referenced by LmunuAnumu(), LoadConfig(), and XSec().

TString genie::NievesQELCCPXSec::fTensorsOutFile
mutableprivate

file to print tensors to

Definition at line 157 of file NievesQELCCPXSec.h.

Referenced by LmunuAnumu().

double genie::NievesQELCCPXSec::fVc
mutableprivate

Definition at line 158 of file NievesQELCCPXSec.h.

Referenced by LmunuAnumu().

double genie::NievesQELCCPXSec::fXSecCCScale
private

external xsec scaling factor for CC

Definition at line 73 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::NievesQELCCPXSec::fXSecEMScale
private

external xsec scaling factor for EM

Definition at line 75 of file NievesQELCCPXSec.h.

const XSecIntegratorI* genie::NievesQELCCPXSec::fXSecIntegrator
private

Definition at line 70 of file NievesQELCCPXSec.h.

Referenced by Integral(), and LoadConfig().

double genie::NievesQELCCPXSec::fXSecNCScale
private

external xsec scaling factor for NC

Definition at line 74 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

The documentation for this class was generated from the following files:
Generated on Wed Jul 17 2024 14:47:01 for GENIEGenerator by   doxygen 1.8.5