GENIEGenerator
Loading...
Searching...
No Matches
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:
[legend]
Collaboration diagram for genie::NievesQELCCPXSec:
[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.
double Integral (const Interaction *i) const
bool ValidProcess (const Interaction *i) const
 Can this cross section algorithm handle the input process?
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?
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.
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status.
virtual bool AllowReconfig (void) const
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm.
virtual void SetId (const AlgId &id)
 Set algorithm ID.
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.

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)
double fXSecCCScale
 external xsec scaling factor for CC
double fXSecNCScale
 external xsec scaling factor for NC
double fXSecEMScale
 external xsec scaling factor for EM
const QvalueShifter * fQvalueShifter
 Optional algorithm to retrieve the qvalue shift for a given target.
double fhbarc
 hbar*c in GeV*fm
bool fRPA
 use RPA corrections
bool fCoulomb
 use Coulomb corrections
const NuclearModelI * fNuclModel
 Nuclear Model for integration.
bool fLFG
const FermiMomentumTable * fKFTable
string fKFTableName
QELEvGen_BindingMode_t fIntegralNucleonBindingMode
double fEnergyCutOff
bool fDoPauliBlocking
 Whether to apply Pauli blocking in XSec()
const genie::PauliBlocker * fPauliBlocker
 The PauliBlocker instance to use to apply that correction.
double fR0
double fCoulombScale
 Scaling factor for the Coulomb potential.
Nieves_Coulomb_Rmax_t fCoulombRmaxMode
bool fCompareNievesTensors
 print tensors
TString fTensorsOutFile
 file to print tensors to
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.
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.
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.
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
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership
int MergeTopRegistry (const Registry &r)
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships.
Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...)
AlgId fID
 algorithm name and configuration set
vector< Registry * > fConfVect
vector< bool > fOwnerships
 ownership for every registry in fConfVect
AlgStatus_t fStatus
 algorithm execution status
AlgMap * fOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool)

Detailed Description

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

References:\n Physical Review C 70, 055503 (2004)
Author
Joe Johnston, University of Pittsburgh Steven Dytman, University of Pittsburgh
Created:\n April 2016
License:\n Copyright (c) 2003-2025, 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() [1/2]

NievesQELCCPXSec::NievesQELCCPXSec ( )

Definition at line 53 of file NievesQELCCPXSec.cxx.

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

References genie::XSecAlgorithmI::XSecAlgorithmI().

◆ NievesQELCCPXSec() [2/2]

NievesQELCCPXSec::NievesQELCCPXSec ( string config)

Definition at line 59 of file NievesQELCCPXSec.cxx.

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

References genie::XSecAlgorithmI::XSecAlgorithmI().

◆ ~NievesQELCCPXSec()

NievesQELCCPXSec::~NievesQELCCPXSec ( )
virtual

Definition at line 65 of file NievesQELCCPXSec.cxx.

66{
67
68 }

Member Function Documentation

◆ CNCTCLimUcalc()

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.

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::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::NievesQELCCPXSec::fhbarc
double fhbarc
hbar*c in GeV*fm
Definition NievesQELCCPXSec.h:78
genie::NievesQELCCPXSec::deltaLindhard
std::complex< double > deltaLindhard(double q0gev, double dqgev, double rho, double kFgev) const
Definition NievesQELCCPXSec.cxx:745
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
genie::NievesQELCCPXSec::fKFTable
const FermiMomentumTable * fKFTable
Definition NievesQELCCPXSec.h:88
genie::utils::nuclear::Density
double Density(double r, int A, double ring=0.)
Definition NuclearUtils.cxx:396
genie::kPdgProton
const int kPdgProton
Definition PDGCodes.h:81
genie::kPdgNeutron
const int kPdgNeutron
Definition PDGCodes.h:83

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

Referenced by LmunuAnumu().

◆ CompareNievesTensors()

void NievesQELCCPXSec::CompareNievesTensors ( const Interaction * i) const
private

Definition at line 1223 of file NievesQELCCPXSec.cxx.

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
pINFO
#define pINFO
Definition Messenger.h:62
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
genie::InitialState::Tgt
const Target & Tgt(void) const
Definition InitialState.h:66
genie::InitialState::ProbePdg
int ProbePdg(void) const
Definition InitialState.h:64
genie::Interaction::FSPrimLepton
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
Definition Interaction.cxx:126
genie::Interaction::KinePtr
Kinematics * KinePtr(void) const
Definition Interaction.h:76
genie::Kinematics::SetQ2
void SetQ2(double Q2, bool selected=false)
Definition Kinematics.cxx:255
genie::NievesQELCCPXSec::fCoulomb
bool fCoulomb
use Coulomb corrections
Definition NievesQELCCPXSec.h:82
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::vcr
double vcr(const Target *target, double r) const
Definition NievesQELCCPXSec.cxx:832
genie::NievesQELCCPXSec::fRPA
bool fRPA
use RPA corrections
Definition NievesQELCCPXSec.h:81
genie::NievesQELCCPXSec::fTensorsOutFile
TString fTensorsOutFile
file to print tensors to
Definition NievesQELCCPXSec.h:157
genie::Target::HitNucMass
double HitNucMass(void) const
Definition Target.cxx:233
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition PDGUtils.cxx:110

References fCoulomb, fCoulombFactor, fFormFactors, fRPA, genie::Interaction::FSPrimLepton(), fTensorsOutFile, fVc, genie::Target::HitNucMass(), genie::pdg::IsNeutrino(), genie::Interaction::KinePtr(), LmunuAnumu(), LOG, pINFO, genie::InitialState::ProbePdg(), genie::Kinematics::SetQ2(), genie::InitialState::Tgt(), and vcr().

◆ Configure() [1/2]

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.

389{
390 Algorithm::Configure(config);
391 this->LoadConfig();
392}
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition Algorithm.cxx:62

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

◆ Configure() [2/2]

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.

395{
396 Algorithm::Configure(config);
397 this->LoadConfig();
398}

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

◆ deltaLindhard()

std::complex< double > NievesQELCCPXSec::deltaLindhard ( double q0gev,
double dqgev,
double rho,
double kFgev ) const
private

Definition at line 745 of file NievesQELCCPXSec.cxx.

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}

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

Referenced by CNCTCLimUcalc().

◆ Integral()

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.

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

References fXSecIntegrator, LOG, and pFATAL.

◆ leviCivita()

int NievesQELCCPXSec::leviCivita ( int input[]) const
private

Definition at line 891 of file NievesQELCCPXSec.cxx.

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}

Referenced by LmunuAnumu().

◆ LmunuAnumu()

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.

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}
pDEBUG
#define pDEBUG
Definition Messenger.h:63
pWARN
#define pWARN
Definition Messenger.h:60
genie::NievesQELCCPXSec::fCompareNievesTensors
bool fCompareNievesTensors
print tensors
Definition NievesQELCCPXSec.h:156
genie::NievesQELCCPXSec::leviCivita
int leviCivita(int input[]) const
Definition NievesQELCCPXSec.cxx:891
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::HitNucPdg
int HitNucPdg(void) const
Definition Target.cxx:304
genie::Target::N
int N(void) const
Definition Target.h:69
genie::Target::Z
int Z(void) const
Definition Target.h:68
genie::Target::A
int A(void) const
Definition Target.h:70
genie::Target::HitNucPosition
double HitNucPosition(void) const
Definition Target.h:89
genie::Target::Pdg
int Pdg(void) const
Definition Target.h:71
genie::Target::IsNucleus
bool IsNucleus(void) const
Definition Target.cxx:272
a
const double a
Definition gUpMuFluxGen.cxx:164
genie::controls::kASmallNum
static const double kASmallNum
Definition Controls.h:40
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition PDGUtils.cxx:336

References genie::Target::A(), a, CNCTCLimUcalc(), fCompareNievesTensors, fCoulombFactor, fFormFactors, fRPA, fTensorsOutFile, fVc, 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, and genie::Target::Z().

Referenced by CompareNievesTensors(), and XSec().

◆ LoadConfig()

void NievesQELCCPXSec::LoadConfig ( void )
private

Definition at line 400 of file NievesQELCCPXSec.cxx.

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}
pNOTICE
#define pNOTICE
Definition Messenger.h:61
pERROR
#define pERROR
Definition Messenger.h:59
RgKey
string RgKey
Definition RegistryItemTypeDef.h:33
genie::Algorithm::GetConfig
virtual const Registry & GetConfig(void) const
Definition Algorithm.cxx:246
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
genie::Algorithm::GetParamDef
bool GetParamDef(const RgKey &name, T &p, const T &def) const
genie::Algorithm::SubAlg
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition Algorithm.cxx:345
genie::Algorithm::Id
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition Algorithm.h:98
genie::FermiMomentumTablePool::GetTable
const FermiMomentumTable * GetTable(string name)
Definition FermiMomentumTablePool.cxx:76
genie::FermiMomentumTablePool::Instance
static FermiMomentumTablePool * Instance(void)
Definition FermiMomentumTablePool.cxx:62
genie::NievesQELCCPXSec::fCoulombRmaxMode
Nieves_Coulomb_Rmax_t fCoulombRmaxMode
Definition NievesQELCCPXSec.h:114
genie::NievesQELCCPXSec::fIntegralNucleonBindingMode
QELEvGen_BindingMode_t fIntegralNucleonBindingMode
Definition NievesQELCCPXSec.h:93
genie::NievesQELCCPXSec::fCoulombScale
double fCoulombScale
Scaling factor for the Coulomb potential.
Definition NievesQELCCPXSec.h:110
genie::NievesQELCCPXSec::fNuclModel
const NuclearModelI * fNuclModel
Nuclear Model for integration.
Definition NievesQELCCPXSec.h:84
genie::NievesQELCCPXSec::fXSecNCScale
double fXSecNCScale
external xsec scaling factor for NC
Definition NievesQELCCPXSec.h:74
genie::NievesQELCCPXSec::fDoPauliBlocking
bool fDoPauliBlocking
Whether to apply Pauli blocking in XSec()
Definition NievesQELCCPXSec.h:100
genie::NievesQELCCPXSec::fPauliBlocker
const genie::PauliBlocker * fPauliBlocker
The PauliBlocker instance to use to apply that correction.
Definition NievesQELCCPXSec.h:102
genie::NievesQELCCPXSec::fFormFactorsModel
const QELFormFactorsModelI * fFormFactorsModel
Definition NievesQELCCPXSec.h:69
genie::NievesQELCCPXSec::fXSecCCScale
double fXSecCCScale
external xsec scaling factor for CC
Definition NievesQELCCPXSec.h:73
genie::NievesQELCCPXSec::fQvalueShifter
const QvalueShifter * fQvalueShifter
Optional algorithm to retrieve the qvalue shift for a given target.
Definition NievesQELCCPXSec.h:76
genie::NievesQELCCPXSec::fCos8c2
double fCos8c2
cos^2(cabibbo angle)
Definition NievesQELCCPXSec.h:71
genie::units::fermi
static constexpr double fermi
Definition Units.h:55
genie::utils::StringToQELBindingMode
QELEvGen_BindingMode_t StringToQELBindingMode(const std::string &mode_str)
Definition QELUtils.cxx:194
genie::gAbortingInErr
bool gAbortingInErr
Definition Messenger.cxx:34
genie::kNucmLocalFermiGas
@ kNucmLocalFermiGas
Definition NuclearModel.h:32
genie::kMatchVertexGeneratorRmax
@ kMatchVertexGeneratorRmax
Definition NievesQELCCPXSec.h:39

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, pERROR, pFATAL, pNOTICE, genie::utils::StringToQELBindingMode(), and genie::Algorithm::SubAlg().

Referenced by Configure(), and Configure().

◆ relLindhard()

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.

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

References fhbarc, LOG, pWARN, relLindhardIm(), and ruLinRelX().

Referenced by CNCTCLimUcalc().

◆ relLindhardIm()

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.

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}

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

Referenced by CNCTCLimUcalc(), and relLindhard().

◆ ruLinRelX()

std::complex< double > NievesQELCCPXSec::ruLinRelX ( double q0,
double qm,
double kf,
double m ) const
private

Definition at line 685 of file NievesQELCCPXSec.cxx.

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}

References genie::constants::kPi2.

Referenced by relLindhard().

◆ ValidProcess()

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.

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::ProcessInfo::IsQuasiElastic
bool IsQuasiElastic(void) const
Definition ProcessInfo.cxx:69
genie::pdg::IsNeutron
bool IsNeutron(int pdgc)
Definition PDGUtils.cxx:341
genie::pdg::IsAntiNeutrino
bool IsAntiNeutrino(int pdgc)
Definition PDGUtils.cxx:118
genie::kISkipProcessChk
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition Interaction.h:47

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().

◆ vcr()

double NievesQELCCPXSec::vcr ( const Target * target,
double r ) const
private

Definition at line 832 of file NievesQELCCPXSec.cxx.

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}
func
double func(double x, double y)
Definition gtestBLI2DUnifGrid.cxx:81
genie::utils::gsl::Integration1DimTypeFromString
ROOT::Math::IntegrationOneDim::Type Integration1DimTypeFromString(string type)
Definition GSLUtils.cxx:23

References genie::Target::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 CompareNievesTensors(), and XSec().

◆ XSec()

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.

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 // Check that Q2tilde > 0 (accounting for rounding errors)
275 if ( Q2tilde <= kASmallNum ) {
276 LOG("Nieves", pWARN) << "Q2tilde <= 0, returning xsec = 0.0";
277 return 0.0;
278 }
279
280 // Store Q2tilde in the interaction so that we get the correct
281 // values of the form factors (according to the de Forest prescription)
282 interaction->KinePtr()->SetQ2(Q2tilde);
283
284 double q2 = -Q2tilde;
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::InitialState::GetProbeP4
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
Definition InitialState.cxx:325
genie::KinePhaseSpace::AsString
static string AsString(KinePhaseSpace_t kps)
Definition KinePhaseSpace.h:85
genie::Kinematics::HadSystP4
const TLorentzVector & HadSystP4(void) const
Definition Kinematics.h:66
genie::Kinematics::FSLeptonP4
const TLorentzVector & FSLeptonP4(void) const
Definition Kinematics.h:65
genie::NievesQELCCPXSec::ValidProcess
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Definition NievesQELCCPXSec.cxx:365
genie::ProcessInfo::IsWeakNC
bool IsWeakNC(void) const
Definition ProcessInfo.cxx:208
genie::Target::HitNucP4
const TLorentzVector & HitNucP4(void) const
Definition Target.h:91
genie::XSecAlgorithmI::ValidKinematics
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
Definition XSecAlgorithmI.cxx:40
genie::utils::kinematics::Q2
double Q2(const Interaction *const i)
Definition KineUtils.cxx:1077
genie::utils::kinematics::Jacobian
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition KineUtils.cxx:130
genie::utils::mec::J
double J(double q0, double q3, double Enu, double ml)
Definition MECUtils.cxx:147
genie::utils::res::Mass
double Mass(Resonance_t res)
resonance mass (GeV)
Definition BaryonResUtils.cxx:436
genie::utils::EnergyDeltaFunctionSolutionQEL
double EnergyDeltaFunctionSolutionQEL(const Interaction &inter)
Definition QELUtils.cxx:50
genie::kRfLab
@ kRfLab
Definition RefFrame.h:26
genie::kIAssumeFreeNucleon
const UInt_t kIAssumeFreeNucleon
Definition Interaction.h:49

References genie::KinePhaseSpace::AsString(), genie::utils::EnergyDeltaFunctionSolutionQEL(), fCos8c2, fCoulomb, fDoPauliBlocking, fFormFactors, fPauliBlocker, fQvalueShifter, fRPA, genie::Interaction::FSPrimLepton(), fXSecCCScale, fXSecNCScale, genie::InitialState::GetProbeP4(), 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::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::Interaction::RecoilNucleon(), genie::Interaction::RecoilNucleonPdg(), genie::Kinematics::SetQ2(), genie::InitialState::Tgt(), genie::XSecAlgorithmI::ValidKinematics(), ValidProcess(), vcr(), and genie::Target::Z().

Member Data Documentation

◆ fCompareNievesTensors

bool genie::NievesQELCCPXSec::fCompareNievesTensors
mutableprivate

print tensors

Definition at line 156 of file NievesQELCCPXSec.h.

Referenced by LmunuAnumu(), and LoadConfig().

◆ fCos8c2

double genie::NievesQELCCPXSec::fCos8c2
private

cos^2(cabibbo angle)

Definition at line 71 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

◆ fCoulomb

bool genie::NievesQELCCPXSec::fCoulomb
private

use Coulomb corrections

Definition at line 82 of file NievesQELCCPXSec.h.

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

◆ fCoulombFactor

double genie::NievesQELCCPXSec::fCoulombFactor
private

Definition at line 158 of file NievesQELCCPXSec.h.

Referenced by CompareNievesTensors(), and LmunuAnumu().

◆ fCoulombRmaxMode

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().

◆ fCoulombScale

double genie::NievesQELCCPXSec::fCoulombScale
private

Scaling factor for the Coulomb potential.

Definition at line 110 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and vcr().

◆ fDoPauliBlocking

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().

◆ fEnergyCutOff

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().

◆ fFormFactors

QELFormFactors genie::NievesQELCCPXSec::fFormFactors
mutableprivate

Definition at line 68 of file NievesQELCCPXSec.h.

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

◆ fFormFactorsModel

const QELFormFactorsModelI* genie::NievesQELCCPXSec::fFormFactorsModel
private

Definition at line 69 of file NievesQELCCPXSec.h.

Referenced by LoadConfig().

◆ fhbarc

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().

◆ fIntegralNucleonBindingMode

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().

◆ fKFTable

const FermiMomentumTable* genie::NievesQELCCPXSec::fKFTable
private

Definition at line 88 of file NievesQELCCPXSec.h.

Referenced by CNCTCLimUcalc(), and LoadConfig().

◆ fKFTableName

string genie::NievesQELCCPXSec::fKFTableName
private

Definition at line 89 of file NievesQELCCPXSec.h.

Referenced by LoadConfig().

◆ fLFG

bool genie::NievesQELCCPXSec::fLFG
private

Definition at line 87 of file NievesQELCCPXSec.h.

Referenced by CNCTCLimUcalc(), and LoadConfig().

◆ fNuclModel

const NuclearModelI* genie::NievesQELCCPXSec::fNuclModel
private

Nuclear Model for integration.

Definition at line 84 of file NievesQELCCPXSec.h.

Referenced by LoadConfig().

◆ fPauliBlocker

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().

◆ fQvalueShifter

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().

◆ fR0

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().

◆ fRPA

bool genie::NievesQELCCPXSec::fRPA
mutableprivate

use RPA corrections

Definition at line 81 of file NievesQELCCPXSec.h.

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

◆ fTensorsOutFile

TString genie::NievesQELCCPXSec::fTensorsOutFile
mutableprivate

file to print tensors to

Definition at line 157 of file NievesQELCCPXSec.h.

Referenced by CompareNievesTensors(), and LmunuAnumu().

◆ fVc

double genie::NievesQELCCPXSec::fVc
mutableprivate

Definition at line 158 of file NievesQELCCPXSec.h.

Referenced by CompareNievesTensors(), and LmunuAnumu().

◆ fXSecCCScale

double genie::NievesQELCCPXSec::fXSecCCScale
private

external xsec scaling factor for CC

Definition at line 73 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

◆ fXSecEMScale

double genie::NievesQELCCPXSec::fXSecEMScale
private

external xsec scaling factor for EM

Definition at line 75 of file NievesQELCCPXSec.h.

◆ fXSecIntegrator

const XSecIntegratorI* genie::NievesQELCCPXSec::fXSecIntegrator
private

Definition at line 70 of file NievesQELCCPXSec.h.

Referenced by Integral(), and LoadConfig().

◆ fXSecNCScale

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 for GENIEGenerator by doxygen 1.14.0