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:

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	fXSecScale
	external xsec scaling factor More...

const QvalueShifter *	fQvalueShifter
	Optional algorithm to retrieve the qvalue shift for a given target. More...

double	fhbarc
	hbarc in GeVfm 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)

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

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

Physical Review C 70, 055503 (2004)

Author: Joe Johnston, University of Pittsburgh Steven Dytman, University of Pittsburgh

April 2016

Definition at line 48 of file NievesQELCCPXSec.h.

Constructor & Destructor Documentation

NievesQELCCPXSec::NievesQELCCPXSec ( )

Definition at line 53 of file NievesQELCCPXSec.cxx.

                                    :
 XSecAlgorithmI("genie::NievesQELCCPXSec")
 {
 
 }

NievesQELCCPXSec::NievesQELCCPXSec ( string config )

Definition at line 59 of file NievesQELCCPXSec.cxx.

                                                 :
 XSecAlgorithmI("genie::NievesQELCCPXSec", config)
 {
 
 }

NievesQELCCPXSec::~NievesQELCCPXSec ( )

virtual

Definition at line 65 of file NievesQELCCPXSec.cxx.

 {
 
  }

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 512 of file NievesQELCCPXSec.cxx.

 {
   if ( tgtIsNucleus && !assumeFreeNucleon ) {
     double dq = qTildeP4.Vect().Mag();
     double dq2 = TMath::Power(dq,2);
     double q2 = 1 * qTildeP4.Mag2();
     //Terms for polarization coefficients CN,CT, and CL
     double hbarc2 = TMath::Power(fhbarc,2);
     double c0 = 0.380/fhbarc;//Constant for CN in natural units
 
     //Density gives the nuclear density, normalized to 1
     //Input radius r must be in fm
     double rhop = nuclear::Density(r,A)*Z;
     double rhon = nuclear::Density(r,A)*N;
     double rho = rhop + rhon;
     double rho0 = A*nuclear::Density(0,A);
 
     double fPrime = (0.33*rho/rho0+0.45*(1-rho/rho0))*c0;
 
     // Get Fermi momenta
     double kF1, kF2;
     if(fLFG){
       if(hitNucIsProton){
         kF1 = TMath::Power(3*kPi2*rhop, 1.0/3.0) *fhbarc;
         kF2 = TMath::Power(3*kPi2*rhon, 1.0/3.0) *fhbarc;
       }else{
         kF1 = TMath::Power(3*kPi2*rhon, 1.0/3.0) *fhbarc;
         kF2 = TMath::Power(3*kPi2*rhop, 1.0/3.0) *fhbarc;
       }
     }else{
       if(hitNucIsProton){
         kF1 = fKFTable->FindClosestKF(tgt_pdgc, kPdgProton);
         kF2 = fKFTable->FindClosestKF(tgt_pdgc, kPdgNeutron);
       }else{
         kF1 = fKFTable->FindClosestKF(tgt_pdgc, kPdgNeutron);
         kF2 = fKFTable->FindClosestKF(tgt_pdgc, kPdgProton);
       }
     }
 
     double kF = TMath::Power(1.5*kPi2*rho, 1.0/3.0) *fhbarc;
 
     std::complex<double> imU(relLindhardIm(qTildeP4.E(),dq,kF1,kF2,
                                            M,is_neutrino,t0,r00));
 
     imaginaryU = imag(imU);
 
     std::complex<double> relLin(0,0),udel(0,0);
 
     // By comparison with Nieves' fortran code
     if(imaginaryU < 0.){
       relLin = relLindhard(qTildeP4.E(),dq,kF,M,is_neutrino,imU);
       udel = deltaLindhard(qTildeP4.E(),dq,rho,kF);
     }
     std::complex<double> relLinTot(relLin + udel);
 
   /* CRho = 2
      DeltaRho = 2500 MeV, (2.5 GeV)^2 = 6.25 GeV^2
      mRho = 770 MeV, (0.770 GeV)^2 = 0.5929 GeV^2
      g' = 0.63 */
     double Vt = 0.08*4*kPi/kPionMass2 *
       (2* TMath::Power((6.25-0.5929)/(6.25-q2),2)*dq2/(q2-0.5929) + 0.63);
   /* f^2/4/Pi = 0.08
      DeltaSubPi = 1200 MeV, (1.2 GeV)^2 = 1.44 GeV^2
      g' = 0.63 */
     double Vl = 0.08*4*kPi/kPionMass2 *
       (TMath::Power((1.44-kPionMass2)/(1.44-q2),2)*dq2/(q2-kPionMass2)+0.63);
 
     CN = 1.0/TMath::Power(abs(1.0-fPrime*relLin/hbarc2),2);
 
     CT = 1.0/TMath::Power(abs(1.0-relLinTot*Vt),2);
     CL = 1.0/TMath::Power(abs(1.0-relLinTot*Vl),2);
   }
   else {
     //Polarization Coefficients: all equal to 1.0 for free nucleon
     CN = 1.0;
     CT = 1.0;
     CL = 1.0;
     imaginaryU = 0.0;
   }
 }

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

private

Definition at line 1216 of file NievesQELCCPXSec.cxx.

         {
   Interaction * interaction = new Interaction(*in); // copy in
 
   // Set input values here
   double ein = 0.2;
   double ctl = 0.5;
   double rmaxfrac = 0.25;
 
   bool carbon = false; // true -> C12, false -> Pb208
 
   if(fRPA)
     fTensorsOutFile = "gen.RPA";
   else
     fTensorsOutFile = "gen.noRPA";
 
   // Calculate radius
   bool klave;
   double rp,ap,rn,an;
   if(carbon){
     klave = true;
     rp = 1.692;
     ap = 1.082;
     rn = 1.692;
     an = 1.082;
   }else{
     // Pb208
     klave = false;
     rp = 6.624;
     ap = 0.549;
     rn = 6.890;
     an = 0.549;
   }
   double rmax;
   if(!klave)
     rmax = TMath::Max(rp,rn) + 9.25*TMath::Max(ap,an);
   else
     rmax = TMath::Sqrt(20.0)*TMath::Max(rp,rn);
   double r = rmax *  rmaxfrac;
 
   // Relevant objects and parameters
   //const Kinematics &   kinematics = interaction -> Kine();
   const InitialState & init_state = interaction -> InitState();
   const Target & target = init_state.Tgt();
 
   // Parameters required for LmunuAnumu
   double M  = target.HitNucMass();
   double ml = interaction->FSPrimLepton()->Mass();
   bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
 
   // Iterate over lepton energy (which then affects q, which is passed to
   // LmunuAnumu using in and out NucleonMom
   double delta = (ein-0.025)/100.0;
   for(int it=0;it<100;it++){
     double tmu = it*delta;
     double eout = ml + tmu;
     double pout = TMath::Sqrt(eout*eout-ml*ml);
 
     double pin = TMath::Sqrt(ein*ein); // Assume massless neutrinos
 
     double q0 = ein-eout;
     double dq = TMath::Sqrt(pin*pin+pout*pout-2.0*ctl*pin*pout);
     double q2 = q0*q0-dq*dq;
     interaction->KinePtr()->SetQ2(-q2);
 
     // When this method is called, inNucleonMomOnShell is unused.
     // I can thus provide the calculated values using a null vector for
     // inNucleonMomOnShell. I also need to put qTildeP4 in the z direction, as
     // Nieves does in his paper.
     TLorentzVector qTildeP4(0, 0, dq, q0);
     TLorentzVector inNucleonMomOnShell(0,0,0,0);
 
     // neutrinoMom and leptonMom only directly affect the leptonic tensor, which
     // we are not calculating now. Use them to transfer q.
     TLorentzVector neutrinoMom(0,0,pout+dq,eout+q0);
     TLorentzVector leptonMom(0,0,pout,eout);
 
     if(fCoulomb){ // Use same steps as in XSec()
       // Coulomb potential
       double Vc = vcr(& target, r);
       fVc = Vc;
 
       // Outgoing lepton energy and momentum including coulomb potential
       int sign = is_neutrino ? 1 : -1;
       double El = leptonMom.E();
       double ElLocal = El - sign*Vc;
       if(ElLocal - ml <= 0.0){
         LOG("Nieves",pINFO) << "Event should be rejected. Coulomb effects "
                             << "push kinematics below threshold";
         return;
       }
       double plLocal = TMath::Sqrt(ElLocal*ElLocal-ml*ml);
 
       // Correction factor
       double coulombFactor= plLocal*ElLocal/leptonMom.Vect().Mag()/El;
       fCoulombFactor = coulombFactor; // Store and print
     }
 
     // TODO: apply Coulomb correction to 3-momentum transfer dq
 
     fFormFactors.Calculate(interaction);
     LmunuAnumu(neutrinoMom, inNucleonMomOnShell, leptonMom, qTildeP4,
                M, is_neutrino, target, false);
   }
   return;
 } // END TESTING CODE

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 382 of file NievesQELCCPXSec.cxx.

 {
   Algorithm::Configure(config);
   this->LoadConfig();
 }

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 388 of file NievesQELCCPXSec.cxx.

 {
   Algorithm::Configure(config);
   this->LoadConfig();
 }

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

private

Definition at line 738 of file NievesQELCCPXSec.cxx.

 {
   double q_zero = q0/fhbarc;
   double q_mod = dq/fhbarc;
   double k_fermi = kF/fhbarc;
   //Divide by hbarc in order to use natural units (rho is already in the correct units)
 
   //m = 939/197.3, md = 1232/197.3, mpi = 139/197.3
   double m = 4.7592;
   double md = 6.2433;
   double mpi = 0.7045;
 
   double fdel_f = 2.13;
   double wr = md-m;
   double gamma = 0;
   double gammap = 0;
 
   double q_zero2 =  TMath::Power(q_zero,  2);
   double q_mod2 =   TMath::Power(q_mod,   2);
   double k_fermi2 = TMath::Power(k_fermi, 2);
 
   double m2 =       TMath::Power(m,       2);
   double m4 =       TMath::Power(m,       4);
   double mpi2 =     TMath::Power(mpi,     2);
   double mpi4 =     TMath::Power(mpi,     4);
 
   double fdel_f2 =  TMath::Power(fdel_f,  2);
 
   //For the current code q_zero is always real
   //If q_zero can have an imaginary part then only the real part is used
   //until z and zp are calculated
 
   double s = m2+q_zero2-q_mod2+
     2.0*q_zero *TMath::Sqrt(m2+3.0/5.0*k_fermi2);
 
   if(s>TMath::Power(m+mpi,2)){
     double srot = TMath::Sqrt(s);
     double qcm = TMath::Sqrt(TMath::Power(s,2)+mpi4+m4-2.0*(s*mpi2+s*m2+
                 mpi2*m2)) /(2.0*srot);
     gamma = 1.0/3.0 * 1.0/(4.0*kPi) * fdel_f2*
      TMath::Power(qcm,3)/srot*(m+TMath::Sqrt(m2+TMath::Power(qcm,2)))/mpi2;
   }
   double sp = m2+q_zero2-q_mod2-
     2.0*q_zero *TMath::Sqrt(m2+3.0/5.0*k_fermi2);
 
 
   if(sp > TMath::Power(m+mpi,2)){
     double srotp = TMath::Sqrt(sp);
     double qcmp = TMath::Sqrt(TMath::Power(sp,2)+mpi4+m4-2.0*(sp*mpi2+sp*m2+
                  mpi2*m2))/(2.0*srotp);
     gammap = 1.0/3.0 * 1.0/(4.0*kPi) * fdel_f2*
       TMath::Power(qcmp,3)/srotp*(m+TMath::Sqrt(m2+TMath::Power(qcmp,2)))/mpi2;
   }
   //}//End if statement
   const std::complex<double> iNum(0,1.0);
 
   std::complex<double> z(md/(q_mod*k_fermi)*(q_zero-q_mod2/(2.0*md)
                          -wr +iNum*gamma/2.0));
   std::complex<double> zp(md/(q_mod*k_fermi)*(-q_zero-q_mod2/(2.0*md)
                           -wr +iNum*gammap/2.0));
 
   std::complex<double> pzeta(0.0);
   if(abs(z) > 50.0){
     pzeta = 2.0/(3.0*z)+2.0/(15.0*z*z*z);
   }else if(abs(z) < TMath::Power(10.0,-2)){
     pzeta = 2.0*z-2.0/3.0*z*z*z-iNum*kPi/2.0*(1.0-z*z);
   }else{
     pzeta = z + (1.0-z*z) * log((z+1.0)/(z-1.0))/2.0;
   }
 
   std::complex<double> pzetap(0);
   if(abs(zp) > 50.0){
     pzetap = 2.0/(3.0*zp)+2.0/(15.0*zp*zp*zp);
   }else if(abs(zp) < TMath::Power(10.0,-2)){
     pzetap = 2.0*zp-2.0/3.0*zp*zp*zp-iNum*kPi/2.0*(1.0-zp*zp);
   }else{
     pzetap = zp+ (1.0-zp*zp) * log((zp+1.0)/(zp-1.0))/2.0;
   }
 
   //Multiply by hbarc^2 to give answer in units of GeV^2
   return 2.0/3.0 * rho * md/(q_mod*k_fermi) * (pzeta +pzetap) * fdel_f2 *
     TMath::Power(fhbarc,2);
 }

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 342 of file NievesQELCCPXSec.cxx.

 {
   // If we're using the new spline generation method (which integrates
   // over the kPSQELEvGen phase space used by QELEventGenerator) then
   // let the cross section integrator do all of the work. It's smart
   // enough to handle free nucleon vs. nuclear targets, different
   // nuclear models (including the local Fermi gas model), etc.
   if ( fXSecIntegrator->Id().Name() == "genie::NewQELXSec" ) {
     return fXSecIntegrator->Integrate(this, in);
   }
   else {
     LOG("Nieves", pFATAL) << "Splines for the Nieves CCQE model must be"
       << " generated using genie::NewQELXSec";
     std::exit(1);
   }
 }

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

private

Definition at line 884 of file NievesQELCCPXSec.cxx.

                                                  {
   int copy[4] = {input[0],input[1],input[2],input[3]};
   int permutations = 0;
   int temp;
 
   for(int i=0;i<4;i++){
     for(int j=i+1;j<4;j++){
       //If any two elements are equal return 0
       if(input[i] == input[j])
         return 0;
       //If a larger number is before a smaller one, use permutations
       //(exchanges of two adjacent elements) to move the smaller element
       //so it is before the larger element, eg 2341->2314->2134->1234
       if(copy[i]>copy[j]){
         temp = copy[j];
         for(int k=j;k>i;k--){
           copy[k] = copy[k-1];
           permutations++;
         }
         copy[i] = temp;
       }
     }
   }
 
   if(permutations % 2 == 0){
     return 1;
   }else{
     return -1;
   }
 }

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 918 of file NievesQELCCPXSec.cxx.

 {
   double r = target.HitNucPosition();
   bool tgtIsNucleus = target.IsNucleus();
   int tgt_pdgc = target.Pdg();
   int A = target.A();
   int Z = target.Z();
   int N = target.N();
   bool hitNucIsProton = pdg::IsProton( target.HitNucPdg() );
 
   const double k[4] = {neutrinoMom.E(),neutrinoMom.Px(),neutrinoMom.Py(),neutrinoMom.Pz()};
   const double kPrime[4] = {leptonMom.E(),leptonMom.Px(),
                             leptonMom.Py(),leptonMom.Pz()};
 
   double q2 = qTildeP4.Mag2();
 
   const double q[4] = {qTildeP4.E(),qTildeP4.Px(),qTildeP4.Py(),qTildeP4.Pz()};
   double q0 = q[0];
   double dq = TMath::Sqrt(TMath::Power(q[1],2)+
                           TMath::Power(q[2],2)+TMath::Power(q[3],2));
 
   int sign = (is_neutrino) ? 1 : -1;
 
   // Get the QEL form factors (were calculated before this method was called)
   double F1V   = 0.5*fFormFactors.F1V();
   double xiF2V = 0.5*fFormFactors.xiF2V();
   double FA    = -fFormFactors.FA();
   // According to Nieves' paper, Fp = 2.0*M*FA/(kPionMass2-q2), but Llewelyn-
   // Smith uses Fp = 2.0*M^2*FA/(kPionMass2-q2), so I divide by M
   // This gives units of GeV^-1
   double Fp    = -1.0/M*fFormFactors.Fp();
 
 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
   LOG("Nieves", pDEBUG) << "\n" << fFormFactors;
 #endif
 
   // Calculate auxiliary parameters
   double M2      = TMath::Power(M,     2);
   double FA2     = TMath::Power(FA,    2);
   double Fp2     = TMath::Power(Fp,    2);
   double F1V2    = TMath::Power(F1V,   2);
   double xiF2V2  = TMath::Power(xiF2V, 2);
   double q02     = TMath::Power(q[0],  2);
   double dq2     = TMath::Power(dq,    2);
   double q4      = TMath::Power(q2,    2);
 
   double t0,r00;
   double CN=1.,CT=1.,CL=1.,imU=0;
   CNCTCLimUcalc(qTildeP4, M, r, is_neutrino, tgtIsNucleus,
     tgt_pdgc, A, Z, N, hitNucIsProton, CN, CT, CL, imU,
     t0, r00, assumeFreeNucleon);
 
   if ( imU > kASmallNum )
     return 0.;
 
   if ( !fRPA || assumeFreeNucleon ) {
     CN = 1.0;
     CT = 1.0;
     CL = 1.0;
   }
 
   double tulin[4] = {0.,0.,0.,0.};
   double rulin[4][4] = { {0.,0.,0.,0.},
                          {0.,0.,0.,0.},
                          {0.,0.,0.,0.},
                          {0.,0.,0.,0.} };
 
   // TESTING CODE:
   if(fCompareNievesTensors){
     // Use average values for initial momentum to calculate A, as given
     // in Appendix B of Nieves' paper. T gives average values of components
     // of p, and R gives the average value of two components multiplied
     // together
     double t3 = (0.5*q2 + q0*t0)/dq; // Average pz
 
     // Vector of p
 
     tulin[0] = t0;
     tulin[3] = t3;
 
     // R is a 4x4 matrix, with R[mu][nu] is the average
     // value of p[mu]*p[nu]
     double aR = r00-M2;
     double bR = (q4+4.0*r00*q02+4.0*q2*q0*t0)/(4.0*dq2);
     double gamma = (aR-bR)/2.0;
     double delta = (-aR+3.0*bR)/2.0/dq2;
 
     double r03 = (0.5*q2*t0 + q0*r00)/dq; // Average E(p)*pz
 
     rulin[0][0] = r00;
     rulin[0][3] = r03;
     rulin[1][1] = gamma;
     rulin[2][2] = gamma;
     rulin[3][0] = r03;
     rulin[3][3] = gamma+delta*dq2; // END TESTING CODE
   }
   else {
     // For normal code execulation, tulin is the initial nucleon momentum
     tulin[0] = inNucleonMomOnShell.E();
     tulin[1] = inNucleonMomOnShell.Px();
     tulin[2] = inNucleonMomOnShell.Py();
     tulin[3] = inNucleonMomOnShell.Pz();
 
     for(int i=0; i<4; i++)
       for(int j=0; j<4; j++)
         rulin[i][j] = tulin[i]*tulin[j];
   }
 
   //Additional constants and variables
   const int g[4][4] = {{1,0,0,0},{0,-1,0,0},{0,0,-1,0},{0,0,0,-1}};
   const std::complex<double> iNum(0,1);
   int leviCivitaIndexArray[4];
   double imaginaryPart = 0;
 
   std::complex<double> sum(0.0,0.0);
 
   double kPrimek = k[0]*kPrime[0]-k[1]*kPrime[1]-k[2]*kPrime[2]-k[3]*kPrime[3];
 
   std::complex<double> Lmunu(0.0,0.0),Lnumu(0.0,0.0);
   std::complex<double> Amunu(0.0,0.0),Anumu(0.0,0.0);
 
   // Calculate LmunuAnumu by iterating over mu and nu
   // In each iteration, add LmunuAnumu to sum if mu=nu, and add
   // LmunuAnumu + LnumuAmunu if mu != nu, since we start nu at mu
   double axx=0.,azz=0.,a0z=0.,a00=0.,axy=0.;
   for(int mu=0;mu<4;mu++){
     for(int nu=mu;nu<4;nu++){
       imaginaryPart = 0;
       if(mu == nu){
         //if mu==nu then levi-civita = 0, so imaginary part = 0
         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;
       }else{
         //if mu!=nu, then g[mu][nu] = 0
         //This same leviCivitaIndex array can be used in the else portion when
         //calculating Anumu
         leviCivitaIndexArray[0] = mu;
         leviCivitaIndexArray[1] = nu;
         for(int a=0;a<4;a++){
           for(int b=0;b<4;b++){
             leviCivitaIndexArray[2] = a;
             leviCivitaIndexArray[3] = b;
             imaginaryPart += - leviCivita(leviCivitaIndexArray)*kPrime[a]*k[b];
           }
         }
         //real(Lmunu) is symmetric, and imag(Lmunu) is antisymmetric
         //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);
         Lmunu = g[mu][mu]*kPrime[mu]*g[nu][nu]*k[nu]+g[nu][nu]*kPrime[nu]*g[mu][mu]*k[mu] + iNum*imaginaryPart;
         Lnumu = g[nu][nu]*kPrime[nu]*g[mu][mu]*k[mu]+g[mu][mu]*kPrime[mu]*g[nu][nu]*k[nu ]- iNum*imaginaryPart;
       } // End Lmunu calculation
 
       if(mu ==0 && nu == 0){
         Amunu = 16.0*F1V2*(2.0*rulin[0][0]*CN+2.0*q[0]*tulin[0]+q2/2.0)+
           2.0*q2*xiF2V2*
           (4.0-4.0*rulin[0][0]/M2-4.0*q[0]*tulin[0]/M2-q02*(4.0/q2+1.0/M2)) +
           4.0*FA2*(2.0*rulin[0][0]+2.0*q[0]*tulin[0]+(q2/2.0-2.0*M2))-
           (2.0*CL*Fp2*q2+8.0*FA*Fp*CL*M)*q02-16.0*F1V*xiF2V*(-q2+q02)*CN;
         a00 = real(Amunu); // TESTING CODE
         sum += Lmunu*Amunu;
       }else if(mu == 0 && nu == 3){
         Amunu = 16.0*F1V2*((2.0*rulin[0][3]+tulin[0]*dq)*CN+tulin[3]*q[0])+
           2.0*q2*xiF2V2*
           (-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))+
           4.0*FA2*((2.0*rulin[0][3]+dq*tulin[0])*CL+q[0]*tulin[3])-
           (2.0*CL*Fp2*q2+8.0*FA*Fp*CL*M)*dq*q[0]-
           16.0*F1V*xiF2V*dq*q[0];
         a0z= real(Amunu); // TESTING CODE
         Anumu = Amunu;
         sum += Lmunu*Anumu + Lnumu*Amunu;
       }else if(mu == 3 && nu == 3){
         Amunu = 16.0*F1V2*(2.0*rulin[3][3]+2.0*dq*tulin[3]-q2/2.0)+
           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))+
           4.0*FA2*(2.0*rulin[3][3]+2.0*dq*tulin[3]-(q2/2.0-2.0*CL*M2))-
           (2.0*CL*Fp2*q2+8.0*FA*Fp*CL*M)*dq2-
           16.0*F1V*xiF2V*(q2+dq2);
         azz = real(Amunu); // TESTING CODE
         sum += Lmunu*Amunu;
       }else if(mu ==1 && nu == 1){
         Amunu = 16.0*F1V2*(2.0*rulin[1][1]-q2/2.0)+
           2.0*q2*xiF2V2*(-4.0*CT-4.0*rulin[1][1]/M2) +
           4.0*FA2*(2.0*rulin[1][1]-(q2/2.0-2.0*CT*M2))-
           16.0*F1V*xiF2V*CT*q2;
         axx = real(Amunu); // TESTING CODE
         sum += Lmunu*Amunu;
       }else if(mu == 2 && nu == 2){
         // Ayy not explicitly listed in paper. This is included so rotating the
         // coordinates of k and k' about the z-axis does not change the xsec.
         Amunu = 16.0*F1V2*(2.0*rulin[2][2]-q2/2.0)+
           2.0*q2*xiF2V2*(-4.0*CT-4.0*rulin[2][2]/M2) +
           4.0*FA2*(2.0*rulin[2][2]-(q2/2.0-2.0*CT*M2))-
           16.0*F1V*xiF2V*CT*q2;
         sum += Lmunu*Amunu;
       }else if(mu ==1 && nu == 2){
         Amunu = sign*16.0*iNum*FA*(xiF2V+F1V)*(-dq*tulin[0]*CT + q[0]*tulin[3]);
         Anumu = -Amunu; // Im(A) is antisymmetric
         axy = imag(Amunu); // TESTING CODE
         sum += Lmunu*Anumu+Lnumu*Amunu;
       }
       // All other terms will be 0 because the initial nucleus is at rest and
       // qTilde is in the z direction
 
     } // End loop over nu
   } // End loop over mu
 
   // TESTING CODE
   if(fCompareNievesTensors){
     // get tmu
     double tmugev = leptonMom.E() - leptonMom.Mag();
     // Print Q2, form factors, and tensor elts
     std::ofstream ffstream;
     ffstream.open(fTensorsOutFile, std::ios_base::app);
     if(q0 > 0){
       ffstream << -q2 << "\t" << q[0] << "\t" << dq
                << "\t" << axx << "\t" << azz << "\t" << a0z
                << "\t" << a00 << "\t" << axy << "\t"
                << CT << "\t" << CL << "\t" << CN << "\t"
                << tmugev << "\t" << imU << "\t"
                << F1V << "\t" << xiF2V << "\t"
                << FA << "\t" << Fp << "\t"
                << tulin[0] << "\t"<< tulin[1] << "\t"
                << tulin[2] << "\t"<< tulin[3] << "\t"
                << rulin[0][0]<< "\t"<< rulin[0][1]<< "\t"
                << rulin[0][2]<< "\t"<< rulin[0][3]<< "\t"
                << rulin[1][0]<< "\t"<< rulin[1][1]<< "\t"
                << rulin[1][2]<< "\t"<< rulin[1][3]<< "\t"
                << rulin[2][0]<< "\t"<< rulin[2][1]<< "\t"
                << rulin[2][2]<< "\t"<< rulin[2][3]<< "\t"
                << rulin[3][0]<< "\t"<< rulin[3][1]<< "\t"
                << rulin[3][2]<< "\t"<< rulin[3][3]<< "\t"
                << fVc << "\t" << fCoulombFactor << "\t";
         ffstream << "\n";
     }
     ffstream.close();
   }
   // END TESTING CODE
 
   // Since the real parts of A and L are both symmetric and the imaginary
   // parts are antisymmetric, the contraction should be real
   if ( imag(sum) > kASmallNum )
     LOG("Nieves",pWARN) << "Imaginary part of tensor contraction is nonzero "
                         << "in QEL XSec, real(sum) = " << real(sum)
                         << "imag(sum) = " << imag(sum);
 
   return real(sum);
 }

void NievesQELCCPXSec::LoadConfig ( void )

private

Definition at line 394 of file NievesQELCCPXSec.cxx.

 {
   bool good_config = true ;
   double thc;
   GetParam( "CabibboAngle", thc ) ;
   fCos8c2 = TMath::Power(TMath::Cos(thc), 2);
 
   // Cross section scaling factor
   GetParam( "QEL-CC-XSecScale", fXSecScale ) ;
 
   // hbarc for unit conversion, GeV*fm
   fhbarc = kLightSpeed*kPlankConstant/genie::units::fermi;
 
   // load QEL form factors model
   fFormFactorsModel = dynamic_cast<const QELFormFactorsModelI *> (
     this->SubAlg("FormFactorsAlg") );
   assert(fFormFactorsModel);
   fFormFactors.SetModel( fFormFactorsModel ); // <-- attach algorithm
 
   // load XSec Integrator
   fXSecIntegrator = dynamic_cast<const XSecIntegratorI*>(
     this->SubAlg("XSec-Integrator") );
   assert(fXSecIntegrator);
 
   // Load settings for RPA and Coulomb effects
 
   // RPA corrections will not affect a free nucleon
   GetParamDef("RPA", fRPA, true ) ;
 
   // Coulomb Correction- adds a correction factor, and alters outgoing lepton
   // 3-momentum magnitude (but not direction)
   // Correction only becomes sizeable near threshold and/or for heavy nuclei
   GetParamDef( "Coulomb", fCoulomb, true ) ;
 
   LOG("Nieves", pNOTICE) << "RPA=" << fRPA << ", useCoulomb=" << fCoulomb;
 
   // Get nuclear model for use in Integral()
   RgKey nuclkey = "IntegralNuclearModel";
   fNuclModel = dynamic_cast<const NuclearModelI *> (this->SubAlg(nuclkey));
   assert(fNuclModel);
 
   // Check if the model is a local Fermi gas
   fLFG = fNuclModel->ModelType(Target()) == kNucmLocalFermiGas;
 
   if(!fLFG){
     // get the Fermi momentum table for relativistic Fermi gas
     GetParam( "FermiMomentumTable", fKFTableName ) ;
 
     fKFTable = 0;
     FermiMomentumTablePool * kftp = FermiMomentumTablePool::Instance();
     fKFTable = kftp->GetTable( fKFTableName );
     assert( fKFTable );
   }
 
   // TESTING CODE
   GetParamDef( "PrintDebugData", fCompareNievesTensors, false ) ;
   // END TESTING CODE
 
   // Nuclear radius parameter (R = R0*A^(1/3)) to use when computing
   // the maximum radius to use to integrate the Coulomb potential
   GetParam("NUCL-R0", fR0) ; // fm
 
   std::string temp_mode;
   GetParamDef( "RmaxMode", temp_mode, std::string("VertexGenerator") ) ;
 
   // Translate the string setting the Rmax mode to the appropriate
   // enum value, or complain if one couldn't be found
   if ( temp_mode == "VertexGenerator" ) {
     fCoulombRmaxMode = kMatchVertexGeneratorRmax;
   }
   else if ( temp_mode == "Nieves" ) {
     fCoulombRmaxMode = kMatchNieves;
   }
   else {
     LOG("Nieves", pFATAL) << "Unrecognized setting \"" << temp_mode
       << "\" requested for the RmaxMode parameter in the"
       << " configuration for NievesQELCCPXSec";
     gAbortingInErr = true;
     std::exit(1);
   }
 
   // Method to use to calculate the binding energy of the initial hit nucleon when
   // generating splines
   std::string temp_binding_mode;
   GetParamDef( "IntegralNucleonBindingMode", temp_binding_mode, std::string("UseNuclearModel") );
   fIntegralNucleonBindingMode = genie::utils::StringToQELBindingMode( temp_binding_mode );
 
   // Cutoff energy for integrating over nucleon momentum distribution (above this
   // lab-frame probe energy, the effects of Fermi motion and binding energy
   // are taken to be negligible for computing the total cross section)
   GetParamDef("IntegralNuclearInfluenceCutoffEnergy", fEnergyCutOff, 2.5 ) ;
 
   // Get PauliBlocker for possible use in XSec()
   fPauliBlocker = dynamic_cast<const PauliBlocker*>( this->SubAlg("PauliBlockerAlg") );
   assert( fPauliBlocker );
 
   // Decide whether or not it should be used in XSec()
   GetParamDef( "DoPauliBlocking", fDoPauliBlocking, true );
 
   // Read optional QvalueShifter:
   fQvalueShifter = nullptr;
   if( GetConfig().Exists("QvalueShifterAlg") ) {
     fQvalueShifter = dynamic_cast<const QvalueShifter *> ( this->SubAlg("QvalueShifterAlg") );
     if( !fQvalueShifter ) {
       good_config = false ;
       LOG("NievesQELCCPXSec", pERROR) << "The required QvalueShifterAlg does not exist. AlgID is : " << SubAlg("QvalueShifterAlg")->Id() ;
     }
   }
 
   if( ! good_config ) {
     LOG("NievesQELCCPXSec", pERROR) << "Configuration has failed.";
     exit(78) ;
   }
 
   // Scaling factor for the Coulomb potential
   GetParamDef( "CoulombScale", fCoulombScale, 1.0 );
 }

std::complex< double > NievesQELCCPXSec::relLindhard	(	double	q0gev,
		double	dqgev,
		double	kFgev,
		double	M,
		bool	isNeutrino,
		std::complex< double >	relLindIm
	)		const

private

Definition at line 655 of file NievesQELCCPXSec.cxx.

 {
   double q0 = q0gev/fhbarc;
   double qm = dqgev/fhbarc;
   double kf = kFgev/fhbarc;
   double m = M/fhbarc;
 
   if(q0>qm){
     LOG("Nieves", pWARN) << "relLindhard() failed";
     return 0.0;
   }
 
   std::complex<double> RealLinRel(ruLinRelX(q0,qm,kf,m)+ruLinRelX(-q0,qm,kf,m));
   double t0,r00;
   std::complex<double> ImLinRel(relLindhardIm(q0gev,dqgev,kFgev,kFgev,M,isNeutrino,t0,r00));
   //Units of GeV^2
   return(RealLinRel*TMath::Power(fhbarc,2) + 2.0*ImLinRel);
 }

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 598 of file NievesQELCCPXSec.cxx.

 {
   double M2 = TMath::Power(M,2);
   double EF1,EF2;
   if(isNeutrino){
     EF1 = TMath::Sqrt(M2+TMath::Power(kFn,2)); //EFn
     EF2 = TMath::Sqrt(M2+TMath::Power(kFp,2)); //EFp
   }else{
     EF1 = TMath::Sqrt(M2+TMath::Power(kFp,2)); //EFp
     EF2 = TMath::Sqrt(M2+TMath::Power(kFn,2)); //EFn
   }
 
   double q2 = TMath::Power(q0,2) - TMath::Power(dq,2);
   double a = (-q0+dq*TMath::Sqrt(1-4.0*M2/q2))/2.0;
   double epsRP = TMath::Max(TMath::Max(M,EF2-q0),a);
 
   // Other theta functions for q are handled by nuclear suppression
   // That is, q0>0 and -q2>0 are always handled, and q0>EF2-EF1 is
   // handled if pauli blocking is on, because otherwise the final
   // nucleon would be below the fermi sea
   //if(fNievesSuppression && !interaction->TestBit(kIAssumeFreeNucleon )
   //&& !EF1-epsRP<0){
   //LOG("Nieves", pINFO) << "Average value of E(p) above Fermi sea";
   //return 0;
   //}else{
   t0 = 0.5*(EF1+epsRP);
   r00 = (TMath::Power(EF1,2)+TMath::Power(epsRP,2)+EF1*epsRP)/3.0;
   std::complex<double> result(0.0,-M2/2.0/kPi/dq*(EF1-epsRP));
   return result;
   //}
 }

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

private

Definition at line 678 of file NievesQELCCPXSec.cxx.

 {
   double q02 = TMath::Power(q0, 2);
   double qm2 = TMath::Power(qm, 2);
   double kf2 = TMath::Power(kf, 2);
   double m2  = TMath::Power(m,  2);
   double m4  = TMath::Power(m,  4);
 
   double ef = TMath::Sqrt(m2+kf2);
   double q2 = q02-qm2;
   double q4 = TMath::Power(q2,2);
   double ds = TMath::Sqrt(1.0-4.0*m2/q2);
   double L1 = log((kf+ef)/m);
   std::complex<double> uL2(
        TMath::Log(TMath::Abs(
                     (ef + q0 - TMath::Sqrt(m2+TMath::Power(kf-qm,2)))/
                     (ef + q0 - TMath::Sqrt(m2 + TMath::Power(kf + qm,2))))) +
        TMath::Log(TMath::Abs(
                     (ef + q0 + TMath::Sqrt(m2 + TMath::Power(kf - qm,2)))/
                     (ef + q0 + TMath::Sqrt(m2 + TMath::Power(kf + qm,2))))));
 
   std::complex<double> uL3(
        TMath::Log(TMath::Abs((TMath::Power(2*kf + q0*ds,2)-qm2)/
                              (TMath::Power(2*kf - q0*ds,2)-qm2))) +
        TMath::Log(TMath::Abs((TMath::Power(kf-ef*ds,2) - (4*m4*qm2)/q4)/
                              (TMath::Power(kf+ef*ds,2) - (4*m4*qm2)/q4))));
 
   std::complex<double> RlinrelX(-L1/(16.0*kPi2)+
                                 uL2*(2.0*ef+q0)/(32.0*kPi2*qm)-
                                 uL3*ds/(64.0*kPi2));
 
   return RlinrelX*16.0*m2;
 }

bool NievesQELCCPXSec::ValidProcess ( const Interaction * i ) const

virtual

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 359 of file NievesQELCCPXSec.cxx.

 {
   if(interaction->TestBit(kISkipProcessChk)) return true;
 
   const InitialState & init_state = interaction->InitState();
   const ProcessInfo &  proc_info  = interaction->ProcInfo();
 
   if(!proc_info.IsQuasiElastic()) return false;
 
   int  nuc = init_state.Tgt().HitNucPdg();
   int  nu  = init_state.ProbePdg();
 
   bool isP   = pdg::IsProton(nuc);
   bool isN   = pdg::IsNeutron(nuc);
   bool isnu  = pdg::IsNeutrino(nu);
   bool isnub = pdg::IsAntiNeutrino(nu);
 
   bool prcok = proc_info.IsWeakCC() && ((isP&&isnub) || (isN&&isnu));
   if(!prcok) return false;
 
   return true;
 }

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

private

Definition at line 825 of file NievesQELCCPXSec.cxx.

                                                                      {
   if(target->IsNucleus()){
     int A = target->A();
     int Z = target->Z();
     double Rmax = 0.;
 
     if ( fCoulombRmaxMode == kMatchNieves ) {
       // Rmax calculated using formula from Nieves' fortran code and default
       // charge and neutron matter density parameters from NuclearUtils.cxx
       if (A > 20) {
         double c = TMath::Power(A,0.35), z = 0.54;
         Rmax = c + 9.25*z;
       }
       else {
         // c = 1.75 for A <= 20
         Rmax = TMath::Sqrt(20.0)*1.75;
       }
     }
     else if ( fCoulombRmaxMode == kMatchVertexGeneratorRmax ) {
       // TODO: This solution is fragile. If the formula used by VertexGenerator
       // changes, then this one will need to change too. Switch to using
       // a common function to get Rmax for both.
       Rmax = 3. * fR0 * std::pow(A, 1./3.);
     }
     else {
       LOG("Nieves", pFATAL) << "Unrecognized setting for fCoulombRmaxMode encountered"
         << " in NievesQELCCPXSec::vcr()";
       gAbortingInErr = true;
       std::exit(1);
     }
 
     if(Rcurr >= Rmax){
       LOG("Nieves",pNOTICE) << "Radius greater than maximum radius for coulomb corrections."
                           << " Integrating to max radius.";
       Rcurr = Rmax;
     }
 
     ROOT::Math::IBaseFunctionOneDim * func = new
       utils::gsl::wrap::NievesQELvcrIntegrand(Rcurr,A,Z);
     ROOT::Math::IntegrationOneDim::Type ig_type =
       utils::gsl::Integration1DimTypeFromString("adaptive");
 
     double abstol = 1; // We mostly care about relative tolerance;
     double reltol = 1E-4;
     int nmaxeval = 100000;
     ROOT::Math::Integrator ig(*func,ig_type,abstol,reltol,nmaxeval);
     double result = ig.Integral(0,Rmax);
     delete func;
 
     // Multiply by Z to normalize densities to number of protons
     // Multiply by hbarc to put result in GeV instead of fm
     // Multiply by an extra configurable scaling factor that defaults to unity
     return -kAem*4*kPi*result*fhbarc*fCoulombScale;
   }else{
     // If target is not a nucleus the potential will be 0
     return 0.0;
   }
 }

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.

 {
   /*// TESTING CODE:
   // The first time this method is called, output tensor elements and other
   // kinmeatics variables for various kinematics. This can the be compared
   // to Nieves' fortran code for validation purposes
   if(fCompareNievesTensors){
     LOG("Nieves",pNOTICE) << "Printing tensor elements for specific "
                           << "kinematics for testing purposes";
     CompareNievesTensors(interaction);
     fCompareNievesTensors = false;
     exit(0);
   }
   // END TESTING CODE*/
 
 
   if ( !this->ValidProcess   (interaction) ) return 0.;
   if ( !this->ValidKinematics(interaction) ) return 0.;
 
   // Get kinematics and init-state parameters
   const Kinematics &   kinematics = interaction -> Kine();
   const InitialState & init_state = interaction -> InitState();
   const Target & target = init_state.Tgt();
 
   // HitNucMass() looks up the PDGLibrary (on-shell) value for the initial
   // struck nucleon
   double mNi = target.HitNucMass();
 
   // Hadronic matrix element for CC neutrino interactions should really use
   // the "nucleon mass," i.e., the mean of the proton and neutrino masses.
   // This expression would also work for NC and EM scattering (since the
   // initial and final on-shell nucleon masses would be the same)
   double mNucleon = ( mNi + interaction->RecoilNucleon()->Mass() ) / 2.;
 
   // Create a copy of the struck nucleon 4-momentum that is forced
   // to be on-shell (this will be needed later for the tensor contraction,
   // in which the nucleon is treated in this way)
   double inNucleonOnShellEnergy = std::sqrt( std::pow(mNi, 2)
     + std::pow(target.HitNucP4().P(), 2) );
 
   // The Nieves CCQE model follows the de Forest prescription: free nucleon
   // (i.e., on-shell) form factors and spinors are used, but an effective
   // value of the 4-momentum transfer "qTilde" is used when computing the
   // contraction of the hadronic tensor. See comments in the
   // FullDifferentialXSec() method of LwlynSmithQELCCPXSec for more details.
   TLorentzVector inNucleonMomOnShell( target.HitNucP4().Vect(),
     inNucleonOnShellEnergy );
 
   // Get the four kinematic vectors and caluclate GFactor
   // Create copies of all kinematics, so they can be rotated
   // and boosted to the nucleon rest frame (because the tensor
   // constraction below only applies for the initial nucleon
   // at rest and q in the z direction)
 
   // All 4-momenta should already be stored, with the hit nucleon off-shell
   // as appropriate
   TLorentzVector* tempNeutrino = init_state.GetProbeP4(kRfLab);
   TLorentzVector neutrinoMom = *tempNeutrino;
   delete tempNeutrino;
   TLorentzVector inNucleonMom = target.HitNucP4();
   TLorentzVector leptonMom = kinematics.FSLeptonP4();
   TLorentzVector outNucleonMom = kinematics.HadSystP4();
 
   // Apply Pauli blocking if enabled
   if ( fDoPauliBlocking && target.IsNucleus() && !interaction->TestBit(kIAssumeFreeNucleon) ) {
     int final_nucleon_pdg = interaction->RecoilNucleonPdg();
     double kF = fPauliBlocker->GetFermiMomentum(target, final_nucleon_pdg,
       target.HitNucPosition());
     double pNf = outNucleonMom.P();
     if ( pNf < kF ) return 0.;
   }
 
   // Use the lab kinematics to calculate the Gfactor, in order to make
   // the XSec differential in initial nucleon momentum and energy
   // Divide by 4.0 because Nieves' conventions for the leptonic and hadronic
   // tensor contraction differ from LwlynSmith by a factor of 4
   double Gfactor = kGF2*fCos8c2 / (8.0*kPi*kPi*inNucleonOnShellEnergy
     *neutrinoMom.E()*outNucleonMom.E()*leptonMom.E()) / 4.0;
 
   // Calculate Coulomb corrections
   double ml = interaction->FSPrimLepton()->Mass();
   double ml2 = TMath::Power(ml, 2);
   double coulombFactor = 1.0;
   double plLocal = leptonMom.P();
 
   bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
   double r = target.HitNucPosition();
 
   if ( fCoulomb ) {
     // Coulomb potential
     double Vc = vcr(& target, r);
 
     // Outgoing lepton energy and momentum including Coulomb potential
     int sign = is_neutrino ? 1 : -1;
     double El = leptonMom.E();
     double pl = leptonMom.P();
     double ElLocal = El - sign*Vc;
 
     if ( ElLocal - ml <= 0. ) {
       LOG("Nieves", pDEBUG) << "Event should be rejected. Coulomb effects"
         << " push kinematics below threshold. Returning xsec = 0.0";
       return 0.0;
     }
 
     // The Coulomb correction factor blows up as pl -> 0. To guard against
     // unphysically huge corrections here, require that the lepton kinetic energy
     // (at infinity) is larger than the magnitude of the Coulomb potential
     // (should be around a few MeV)
     double KEl = El - ml;
     if ( KEl <= std::abs(Vc) ) {
       LOG("Nieves", pDEBUG) << "Outgoing lepton has a very small kinetic energy."
         << " Protecting against near-singularities in the Coulomb correction"
         << " factor by returning xsec = 0.0";
       return 0.0;
     }
 
     // Local value of the lepton 3-momentum magnitude for the Coulomb
     // correction
     plLocal = TMath::Sqrt( ElLocal*ElLocal - ml2 );
 
     // Correction factor
     coulombFactor= (plLocal * ElLocal) / (pl * El);
 
   }
 
   // When computing the contraction of the leptonic and hadronic tensors,
   // we need to use an effective value of the 4-momentum transfer q.
   // The energy transfer (q0) needs to be modified to account for the binding
   // energy of the struck nucleon, while the 3-momentum transfer needs to
   // be corrected for Coulomb effects.
   //
   // See the original Valencia model paper:
   // https://journals.aps.org/prc/abstract/10.1103/PhysRevC.70.055503
 
   double q0Tilde = outNucleonMom.E() - inNucleonMomOnShell.E();
 
   // Shift the q0Tilde if required:
   if( fQvalueShifter ) q0Tilde += q0Tilde * fQvalueShifter->Shift(*interaction) ;
 
   // If binding energy effects pull us into an unphysical region, return
   // zero for the differential cross section
   if ( q0Tilde <= 0. && target.IsNucleus() && !interaction->TestBit(kIAssumeFreeNucleon) ) return 0.;
 
   // Note that we're working in the lab frame (i.e., the rest frame
   // of the target nucleus). We can therefore use Nieves' explicit
   // form of the Amunu tensor if we rotate the 3-momenta so that
   // qTilde is in the +z direction
   TVector3 neutrinoMom3 = neutrinoMom.Vect();
   TVector3 leptonMom3 = leptonMom.Vect();
 
   TVector3 inNucleonMom3 = inNucleonMom.Vect();
   TVector3 outNucleonMom3 = outNucleonMom.Vect();
 
   // If Coulomb corrections are being used, adjust the lepton 3-momentum used
   // to get q3VecTilde so that its magnitude matches the local
   // Coulomb-corrected value calculated earlier. Note that, although the
   // treatment of Coulomb corrections by Nieves et al. doesn't change the
   // direction of the lepton 3-momentum, it *does* change the direction of the
   // 3-momentum transfer, and so the correction should be applied *before*
   // rotating coordinates into a frame where q3VecTilde lies along the positive
   // z axis.
   TVector3 leptonMomCoulomb3 = (! fCoulomb ) ? leptonMom3
     : plLocal * leptonMom3 * (1. / leptonMom3.Mag());
   TVector3 q3VecTilde = neutrinoMom3 - leptonMomCoulomb3;
 
   // Find the rotation angle needed to put q3VecTilde along z
   TVector3 zvec(0.0, 0.0, 1.0);
   TVector3 rot = ( q3VecTilde.Cross(zvec) ).Unit(); // Vector to rotate about
   // Angle between the z direction and q
   double angle = zvec.Angle( q3VecTilde );
 
   // Handle the edge case where q3VecTilde is along -z, so the
   // cross product above vanishes
   if ( q3VecTilde.Perp() == 0. && q3VecTilde.Z() < 0. ) {
     rot = TVector3(0., 1., 0.);
     angle = kPi;
   }
 
   // Rotate if the rotation vector is not 0
   if ( rot.Mag() >= kASmallNum ) {
 
     neutrinoMom3.Rotate(angle,rot);
     neutrinoMom.SetVect(neutrinoMom3);
 
     leptonMom3.Rotate(angle,rot);
     leptonMom.SetVect(leptonMom3);
 
     inNucleonMom3.Rotate(angle,rot);
     inNucleonMom.SetVect(inNucleonMom3);
     inNucleonMomOnShell.SetVect(inNucleonMom3);
 
     outNucleonMom3.Rotate(angle,rot);
     outNucleonMom.SetVect(outNucleonMom3);
 
   }
 
   // Calculate q and qTilde
   TLorentzVector qP4 = neutrinoMom - leptonMom;
   TLorentzVector qTildeP4(0., 0., q3VecTilde.Mag(), q0Tilde);
 
   double Q2 = -1. * qP4.Mag2();
   double Q2tilde = -1. * qTildeP4.Mag2();
 
   // Store Q2tilde in the interaction so that we get the correct
   // values of the form factors (according to the de Forest prescription)
   interaction->KinePtr()->SetQ2(Q2tilde);
 
   double q2 = -Q2tilde;
 
   // Check that q2 < 0 (accounting for rounding errors)
   if ( q2 >= kASmallNum ) {
     LOG("Nieves", pWARN) << "q2 >= 0, returning xsec = 0.0";
     return 0.0;
   }
 
   // Calculate form factors
   fFormFactors.Calculate( interaction );
 
   // Now that the form factors have been calculated, store Q2
   // in the event instead of Q2tilde
   interaction->KinePtr()->SetQ2( Q2 );
 
   // Do the contraction of the leptonic and hadronic tensors. See the
   // RPA-corrected expressions for the hadronic tensor elements in appendix A
   // of Phys. Rev. C 70, 055503 (2004). Note that the on-shell 4-momentum of
   // the initial struck nucleon should be used in the calculation, as well as
   // the effective 4-momentum transfer q tilde (corrected for the nucleon
   // binding energy and Coulomb effects)
   double LmunuAnumuResult = LmunuAnumu(neutrinoMom, inNucleonMomOnShell,
     leptonMom, qTildeP4, mNucleon, is_neutrino, target,
     interaction->TestBit( kIAssumeFreeNucleon ));
 
   // Calculate xsec
   double xsec = Gfactor*coulombFactor*LmunuAnumuResult;
 
   // Apply the factor that arises from elimination of the energy-conserving
   // delta function
   xsec *= genie::utils::EnergyDeltaFunctionSolutionQEL( *interaction );
 
   // Apply given scaling factor
   xsec *= fXSecScale;
 
   LOG("Nieves",pDEBUG) << "TESTING: RPA=" << fRPA
                        << ", Coulomb=" << fCoulomb
                        << ", q2 = " << q2 << ", xsec = " << xsec;
 
   //----- The algorithm computes dxsec/dQ2 or kPSQELEvGen
   //      Check whether variable tranformation is needed
   if ( kps != kPSQELEvGen ) {
 
     // Compute the appropriate Jacobian for transformation to the requested
     // phase space
     double J = utils::kinematics::Jacobian(interaction, kPSQELEvGen, kps);
 
 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
     LOG("Nieves", pDEBUG)
      << "Jacobian for transformation to: "
                   << KinePhaseSpace::AsString(kps) << ", J = " << J;
 #endif
     xsec *= J;
   }
 
   // Number of scattering centers in the target
   int nucpdgc = target.HitNucPdg();
   int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
 
   xsec *= NNucl; // nuclear xsec
 
   return xsec;
 }

Member Data Documentation

bool genie::NievesQELCCPXSec::fCompareNievesTensors

mutableprivate

print tensors

Definition at line 154 of file NievesQELCCPXSec.h.

double genie::NievesQELCCPXSec::fCos8c2

private

cos^2(cabibbo angle)

Definition at line 71 of file NievesQELCCPXSec.h.

bool genie::NievesQELCCPXSec::fCoulomb

private

use Coulomb corrections

Definition at line 80 of file NievesQELCCPXSec.h.

double genie::NievesQELCCPXSec::fCoulombFactor

mutableprivate

Definition at line 156 of file NievesQELCCPXSec.h.

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 112 of file NievesQELCCPXSec.h.

double genie::NievesQELCCPXSec::fCoulombScale

private

Scaling factor for the Coulomb potential.

Definition at line 108 of file NievesQELCCPXSec.h.

bool genie::NievesQELCCPXSec::fDoPauliBlocking

private

Whether to apply Pauli blocking in XSec()

Definition at line 98 of file NievesQELCCPXSec.h.

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 95 of file NievesQELCCPXSec.h.

QELFormFactors genie::NievesQELCCPXSec::fFormFactors

mutableprivate

Definition at line 68 of file NievesQELCCPXSec.h.

const QELFormFactorsModelI* genie::NievesQELCCPXSec::fFormFactorsModel

private

Definition at line 69 of file NievesQELCCPXSec.h.

double genie::NievesQELCCPXSec::fhbarc

private

hbar*c in GeV*fm

Definition at line 76 of file NievesQELCCPXSec.h.

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 91 of file NievesQELCCPXSec.h.

const FermiMomentumTable* genie::NievesQELCCPXSec::fKFTable

private

Definition at line 86 of file NievesQELCCPXSec.h.

string genie::NievesQELCCPXSec::fKFTableName

private

Definition at line 87 of file NievesQELCCPXSec.h.

bool genie::NievesQELCCPXSec::fLFG

private

Definition at line 85 of file NievesQELCCPXSec.h.

const NuclearModelI* genie::NievesQELCCPXSec::fNuclModel

private

Nuclear Model for integration.

Definition at line 82 of file NievesQELCCPXSec.h.

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

private

The PauliBlocker instance to use to apply that correction.

Definition at line 100 of file NievesQELCCPXSec.h.

const QvalueShifter* genie::NievesQELCCPXSec::fQvalueShifter

private

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

Definition at line 74 of file NievesQELCCPXSec.h.

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 105 of file NievesQELCCPXSec.h.

bool genie::NievesQELCCPXSec::fRPA

mutableprivate

use RPA corrections

Definition at line 79 of file NievesQELCCPXSec.h.

TString genie::NievesQELCCPXSec::fTensorsOutFile

mutableprivate

file to print tensors to

Definition at line 155 of file NievesQELCCPXSec.h.

double genie::NievesQELCCPXSec::fVc

mutableprivate

Definition at line 156 of file NievesQELCCPXSec.h.

const XSecIntegratorI* genie::NievesQELCCPXSec::fXSecIntegrator

private

Definition at line 70 of file NievesQELCCPXSec.h.

double genie::NievesQELCCPXSec::fXSecScale

private

external xsec scaling factor

Definition at line 73 of file NievesQELCCPXSec.h.

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

Generator/src/Physics/QuasiElastic/XSection/NievesQELCCPXSec.h
Generator/src/Physics/QuasiElastic/XSection/NievesQELCCPXSec.cxx

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation