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

#include "/cvmfs/nova.opensciencegrid.org/externals/genie/v3_00_06_p01a/Linux64bit+3.10-2.17-e19-debug/GENIE-Generator/src/Physics/QuasiElastic/XSection/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)

virtual bool	ValidKinematics (const Interaction *i) const
	Is the input kinematical point a physically allowed one? More...

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

Protected Member Functions
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 >
bool	GetParamVect (const std::string &comm_name, std::vector< T > &v, unsigned int max, 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
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...

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

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

Nieves_Coulomb_Rmax_t	fCoulombRmaxMode

bool	fCompareNievesTensors
	print tensors More...

TString	fTensorsOutFile
	file to print tensors to More...

double	fVc

double	fCoulombFactor

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

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

NievesQELCCPXSec::NievesQELCCPXSec ( string config )

Definition at line 63 of file NievesQELCCPXSec.cxx.

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

NievesQELCCPXSec::~NievesQELCCPXSec ( )

virtual

Definition at line 69 of file NievesQELCCPXSec.cxx.

 {
 
  }

Member Function Documentation

int Algorithm::AddLowRegistry	(	Registry *	rp,
		bool	owns = `true`
	)

protectedinherited

add registry with lowest priority, also update ownership

Definition at line 601 of file Algorithm.cxx.

Referenced by genie::EventGenerator::Configure().

                                                        {
 
   fConfVect.push_back( rp ) ;
   fOwnerships.push_back( own ) ;
 
   if ( fConfig ) {
     delete fConfig ;
     fConfig = 0 ;
   }
 
   return fConfVect.size() ;
 
 }

int Algorithm::AddTopRegisties	(	const vector< Registry * > &	rs,
		bool	owns = `false`
	)

protectedinherited

Add registries with top priority, also udated Ownerships.

Definition at line 653 of file Algorithm.cxx.

                                                                        {
 
   fConfVect.insert( fConfVect.begin(), rs.begin(), rs.end() ) ;
   
   fOwnerships.insert( fOwnerships.begin(), rs.size(), own ) ;
   
   if ( fConfig ) {
     delete fConfig ;
     fConfig = 0 ;
   }
 
   return fConfVect.size() ;  
 
 }

int Algorithm::AddTopRegistry	(	Registry *	rp,
		bool	owns = `true`
	)

protectedinherited

add registry with top priority, also update ownership

Definition at line 585 of file Algorithm.cxx.

Referenced by genie::EventGeneratorListAssembler::AssembleGeneratorList().

                                                        {  
 
   fConfVect.insert( fConfVect.begin(), rp ) ;
   fOwnerships.insert( fOwnerships.begin(), own ) ;
 
   if ( fConfig ) {
     delete fConfig ;
     fConfig = 0 ;
   }
 
   return fConfVect.size() ;
 
 }

void Algorithm::AdoptConfig ( void )

inherited

Clone the configuration registry looked up from the configuration pool and take its ownership

Definition at line 394 of file Algorithm.cxx.

References Configure(), GetConfig(), LOG, and pNOTICE.

Referenced by genie::Algorithm::AllowReconfig().

                                 {
 
   LOG("Algorithm", pNOTICE)
     << this->Id().Key() << " is taking ownership of its configuration";
   
   // if(fOwnsConfig) {
   //   LOG("Algorithm", pWARN)
   //     << this->Id().Key() << " already owns its configuration!";
   //   return;
   // }
 
   this->Configure( GetConfig() );
 }

void Algorithm::AdoptSubstructure ( void )

inherited

Take ownership of the algorithms subtructure (sub-algorithms,...) by copying them from the AlgFactory pool to the local pool Also bring all the configuration variables to the top level config Registry. This can be used to group together a series of algorithms & their configurations and extract (a clone of) this group from the shared pools. Having a series of algorithms/configurations behaving as a monolithic block, with a single point of configuration (the top level) is to be used when bits & pieces of GENIE are used in isolation for data fitting or reweighting

Definition at line 408 of file Algorithm.cxx.

References genie::AlgFactory::AdoptAlgorithm(), genie::Algorithm::AdoptSubstructure(), GetConfig(), genie::AlgFactory::Instance(), genie::kRgAlg, LOG, pDEBUG, pNOTICE, and genie::RegistryItemI::TypeInfo().

Referenced by genie::Algorithm::AdoptSubstructure(), genie::Algorithm::AllowReconfig(), main(), and testReconfigInOwnedModules().

 {
 // Take ownership of the algorithms subtructure (sub-algorithms,..) by copying 
 // them from the AlgFactory pool to the local pool. Also bring all the 
 // configuration variables to the top level. See the header for more details.
 // A secial naming convention is required for configuration parameter keys 
 // for parameters belonging to sub-algorithms (or sub-algorithms of these 
 // sub-algorithms and so on...). 
 // The convention is: "sub-alg-key/sub-sub-alg-key/.../original name"
 // This is a recursive method: The AdoptSubtructure()of all sub-algorithms is
 // invoked.
 //
   LOG("Algorithm", pNOTICE) 
      << "Algorithm: " << this->Id().Key() << " is adopting its substructure";
 
 //  Registry deep_config;
 //  deep_config.UnLock();
 //  deep_config.SetName(this->Id().Key());
 
   //  deep_config.SetName(this->Id().Config() + ";D");
   //  fID.SetConfig(this->Id().Config() + ";D");
 
   if(fOwnsSubstruc) this->DeleteSubstructure();
 
   fOwnedSubAlgMp = new AlgMap;
   fOwnsSubstruc  = true;
 
   AlgFactory * algf = AlgFactory::Instance();
 
   const RgIMap & rgmap = GetConfig().GetItemMap();
 
   RgIMapConstIter iter = rgmap.begin();
   for( ; iter != rgmap.end(); ++iter) {
 
     RgKey reg_key = iter->first;
     RegistryItemI * ri = iter->second;
 
     if(ri->TypeInfo() == kRgAlg) {
         LOG("Algorithm", pDEBUG) 
                  << "Found sub-algorithm pointed to by " << reg_key;
         RgAlg reg_alg = fConfig->GetAlg(reg_key);
         AlgId id(reg_alg);
 
         LOG("Algorithm", pDEBUG) << "Adopting sub-algorithm = " << id.Key();
         Algorithm * subalg = algf->AdoptAlgorithm(id.Name(),id.Config());
         subalg->AdoptSubstructure();
 
         LOG("Algorithm", pDEBUG) << "Adding sub-algorithm to local pool";
         AlgMapPair key_alg_pair(reg_key, subalg);
         fOwnedSubAlgMp->insert(key_alg_pair);
 
     }
 
   }
 
 
   if ( fConfig ) {
     delete fConfig ;
     fConfig = 0 ;
   }
 
 }

virtual bool genie::Algorithm::AllowReconfig ( void ) const

inlinevirtualinherited

Allow reconfigration after initializaton? Algorithms may opt-out, if reconfiguration is not necessary, to improve event reweighting speed.

Definition at line 106 of file Algorithm.h.

Referenced by genie::AlgFactory::ForceReconfiguration().

106 { return fAllowReconfig; }

genie::Algorithm::fAllowReconfig

bool fAllowReconfig

Definition: Algorithm.h:153

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

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

Referenced by LmunuAnumu().

 {
   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;
   }
 }

AlgCmp_t Algorithm::Compare ( const Algorithm * alg ) const

virtualinherited

Compare with input algorithm.

Definition at line 294 of file Algorithm.cxx.

References genie::AlgId::Config(), genie::Algorithm::Id(), genie::kAlgCmpDiffAlg, genie::kAlgCmpDiffConfig, genie::kAlgCmpIdentical, genie::kAlgCmpUnknown, and genie::AlgId::Name().

Referenced by genie::Algorithm::AllowReconfig().

 {
 // Compares itself with the input algorithm
 
   string alg1    = this->Id().Name();
   string config1 = this->Id().Config();
   string alg2    = algo->Id().Name();
   string config2 = algo->Id().Config();
 
   if(alg1 == alg2)
   {
     if(config1 == config2) return kAlgCmpIdentical;
     else                   return kAlgCmpDiffConfig;
   }
   else return kAlgCmpDiffAlg;
 
   return kAlgCmpUnknown;
 }

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

private

Definition at line 1197 of file NievesQELCCPXSec.cxx.

References delta, genie::Interaction::FSPrimLepton(), genie::Target::HitNucMass(), genie::pdg::IsNeutrino(), it, genie::Interaction::KinePtr(), LOG, pINFO, genie::InitialState::ProbePdg(), q2, r(), genie::Kinematics::SetQ2(), canMan::sign(), ana::Sqrt(), and genie::InitialState::Tgt().

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

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

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

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

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

void Algorithm::DeleteConfig ( void )

protectedinherited

Definition at line 471 of file Algorithm.cxx.

References MECModelEnuComparisons::i.

Referenced by genie::Algorithm::AllowReconfig().

 {
   // there is nothing to delete if the configuration is not owned but is 
   // rather looked up from the configuration pool
   //
   
   for ( unsigned int i = 0 ; i < fConfVect.size() ; ++i ) {
     if ( fOwnerships[i] ) {
       delete fConfVect[i] ;
     }
   }
 
   fConfVect.clear() ;
   fOwnerships.clear() ;
 
   // delete owned configuration registry
 
   if(fConfig) {
     delete fConfig; 
     fConfig=0;
   }
 
 }

void Algorithm::DeleteSubstructure ( void )

protectedinherited

Definition at line 496 of file Algorithm.cxx.

Referenced by genie::Algorithm::AllowReconfig().

 {
   // there is nothing to delete if the sub-algorithms are not owned but rather
   // taken from the AlgFactory's pool
   //
   if(!fOwnsSubstruc) return;
 
   // delete local algorithm pool
   //
   AlgMapIter iter = fOwnedSubAlgMp->begin();
   for( ; iter != fOwnedSubAlgMp->end(); ++iter) {
     Algorithm * alg = iter->second;
     if(alg) {
       delete alg;
       alg=0;
     }
   }
   delete fOwnedSubAlgMp;
   fOwnedSubAlgMp = 0;
 }

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

private

Definition at line 720 of file NievesQELCCPXSec.cxx.

References abs(), fhbarc, kPi, test_ParserArtEvents::log, m, Munits::m2, fetch_tb_beamline_files::md, ana::Sqrt(), and test::z.

Referenced by CNCTCLimUcalc().

 {
   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);
 }

Registry * Algorithm::ExtractLocalConfig ( const Registry & in ) const

protectedinherited

Split an incoming configuration Registry into a block valid for this algorithm Ownership of the returned registry belongs to the algo

Definition at line 518 of file Algorithm.cxx.

References genie::RegistryItemI::Clone(), genie::Registry::GetItemMap(), genie::Registry::Name(), and confusionMatrixTree::out.

Referenced by genie::Algorithm::AllowReconfig().

                                                                     {
 
   const RgIMap & rgmap = in.GetItemMap();
   Registry * out = new Registry( in.Name(), false );
   
   for( RgIMapConstIter reg_iter = rgmap.begin(); 
        reg_iter != rgmap.end(); ++reg_iter ) {
 
     RgKey reg_key = reg_iter->first;
     if( reg_key.find( '/' ) != string::npos) continue; 
 
     // at this point
     // this key is referred to the local algorithm
     // it has to be copied in out;
       
     RegistryItemI * ri = reg_iter->second;
     RgIMapPair key_item_pair( reg_key, ri->Clone() );
     out -> Set(key_item_pair);
 
   }
 
   if ( out -> NEntries() <= 0 ) {
     delete out ;
     out = 0 ;
   }
 
   return out ;
 }

Registry * Algorithm::ExtractLowerConfig	(	const Registry &	in,
		const string &	alg_key
	)		const

protectedinherited

Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key.

Definition at line 549 of file Algorithm.cxx.

References genie::RegistryItemI::Clone(), genie::Registry::GetItemMap(), genie::Registry::Name(), and confusionMatrixTree::out.

Referenced by genie::Algorithm::AllowReconfig().

                                                                                             {
 
   const RgIMap & rgmap = in.GetItemMap();
   Registry * out = new Registry( in.Name(), false );
   
   for( RgIMapConstIter reg_iter = rgmap.begin(); 
        reg_iter != rgmap.end(); ++reg_iter ) {
     
     RgKey reg_key = reg_iter->first;
     if( reg_key.find(alg_key+"/") == string::npos) continue; 
 
     // at this point
     // this key is referred to the sub-algorithm
     // indicated by alg_key: it has to be copied in out;
       
     int new_key_start = reg_key.find_first_of('/')+1;
     RgKey new_reg_key = reg_key.substr( new_key_start, reg_key.length() );
   
     RegistryItemI * ri = reg_iter->second;
     RgIMapPair key_item_pair(new_reg_key, ri->Clone());
     out -> Set(key_item_pair);
 
   }
 
   if ( out -> NEntries() <= 0 ) {
     delete out ;
     out = 0 ;
   }
   
   return out ;
   
 }

void Algorithm::FindConfig ( void )

virtualinherited

Lookup configuration from the config pool Similar logic from void Configure(string)

Definition at line 135 of file Algorithm.cxx.

References gen_hdf5record::config, exit(), genie::AlgConfigPool::FindRegistry(), genie::Registry::GetItemMap(), genie::Registry::GetString(), MECModelEnuComparisons::i, genie::AlgConfigPool::Instance(), it, genie::Registry::ItemIsLocal(), parse_dependency_file_t::list, LOG, pDEBUG, pFATAL, time_estimates::pool, pWARN, moon_position_table_new3::second, genie::utils::str::Split(), string, and APDHVSetting::temp.

 {
 // Finds its configration Registry from the ConfigPool and gets a pointer to
 // it. If the Registry comes from the ConfigPool then the Algorithm does not
 // own its configuration (the ConfigPool singleton keeps the ownership and the
 // responsibility to -eventually- delete all the Registries it instantiates
 // by parsing the XML config files).
 
   DeleteConfig() ;
 
   AlgConfigPool * pool = AlgConfigPool::Instance();
 
   Registry * config = 0 ;
 
   // load the Default config if config is not the default
   if ( fID.Config() != "Default" ) {
     config = pool -> FindRegistry( fID.Name(), "Default" );
     if ( config ) {
       if ( config -> NEntries() > 0 ) {
         AddTopRegistry( config, false ) ;
         LOG("Algorithm", pDEBUG) << "\n" << *config;
       }
     }
   } 
 
   // Load the right config 
   config = pool->FindRegistry(this);
 
   if(!config)
     // notify & keep whatever config Registry was used before.
     LOG("Algorithm", pWARN)
        << "No Configuration available for "
        << this->Id().Key() << " at the ConfigPool";
   else {
     if ( config -> NEntries() > 0 ) {
       AddTopRegistry( config, false ) ;
       LOG("Algorithm", pDEBUG) << "\n" << config;
     }
   }
   
   const string common_key_root = "Common" ;
   std::map<string, string> common_lists;
 
   // Load Common Parameters if key that start with "Common" is found
   for ( unsigned int i = 0 ; i < fConfVect.size() ; ++i ) {
     const Registry & temp = * fConfVect[i] ;
     for ( RgIMapConstIter it = temp.GetItemMap().begin() ; it !=  temp.GetItemMap().end() ; ++it ) {
       
       // check if it is a "Common" entry
       if ( it -> first.find( common_key_root ) == 0 ) {
         // retrieve the type of the common entry
         std::string type = it -> first.substr(common_key_root.size() ) ;
         
         if ( temp.ItemIsLocal( it -> first ) ) {
           
           string temp_list = temp.GetString( it -> first ) ;
           if ( temp_list.length() > 0 ) {
             common_lists[type] = temp_list ;
           }
         }
       }
       
     }
     
   } // loop over the local registries
 
 
   for ( std::map<string, string>::const_iterator it = common_lists.begin() ;
         it != common_lists.end() ; ++it ) {
 
     vector<string> list = str::Split( it -> second , "," ) ;
 
     for ( unsigned int i = 0; i < list.size(); ++i ) {
 
       config = pool -> CommonList( it -> first, list[i] ) ;
 
       if ( ! config ) {
         LOG("Algorithm", pFATAL)
           << "No Commom parameters available for " << it -> first << " list "
           << list[i] << " at the ConfigPool";
 
             exit( 78 ) ;
       }
       else  {
             AddLowRegistry( config, false ) ;
             LOG("Algorithm", pDEBUG) << "Loading " 
                                      << it -> first << " registry " 
                                      << list[i] << " \n" << config;
       }
 
     }
 
   }
 
 
   // Load Tunable from CommonParameters 
   // only if the option is specified in RunOpt
   config = pool -> CommonList( "Param", "Tunable" ) ;
   if ( config ) {
     if ( config -> NEntries() > 0 ) {
       AddTopRegistry( config, false ) ;
       LOG("Algorithm", pDEBUG) << "Loading Tunable registry \n" << config;
     }
   }
   else {
     // notify & keep whatever config Registry was used before.
     LOG("Algorithm", pWARN)
       << "No Tunable parameter set available at the ConfigPool";
   }
 
   if ( fConfig ) {
     delete fConfig ;
     fConfig = 0 ;
   }
 
 }

const Registry & Algorithm::GetConfig ( void ) const

virtualinherited

Get configuration registry Evaluate the summary of the configuration and returns it The summary of a configuration is a merge of all the registries known to the algorithm (see Configure() methods) but every parameter is appearing only once and in case of repetitions, only the parameter from the registry with the highest prioriry is considered.

Definition at line 254 of file Algorithm.cxx.

References febshutoff_auto::end, genie::Algorithm::GetConfig(), MECModelEnuComparisons::i, LOG, pDEBUG, r(), and moon_position_table_new3::second.

Referenced by genie::EventGeneratorListAssembler::AssembleGeneratorList(), GetAlgorithms(), genie::Algorithm::GetConfig(), genie::GRV98LO::GRV98LO(), genie::NewQELXSec::Integrate(), genie::LHAPDF5::LHAPDF5(), genie::IBDXSecMap::LoadConfig(), genie::Decayer::LoadConfig(), genie::PythiaHadronization::LoadConfig(), genie::FGMBodekRitchie::LoadConfig(), genie::NuclearModelMap::LoadConfig(), genie::SmithMonizUtils::LoadConfig(), main(), genie::AlgFactory::Print(), TestPythiaTauDecays(), testReconfigInOwnedModules(), and genie::P33PaschosLalakulichPXSec::XSec().

                                                 {
 
   if ( fConfig ) return * fConfig ;
 
   const_cast<Algorithm*>( this ) -> fConfig = new Registry( fID.Key() + "_summary", false ) ;
   
   // loop and append 
   // understand the append mechanism
   for ( unsigned int i = 0 ; i < fConfVect.size(); ++i ) {
     fConfig -> Append( * fConfVect[i] ) ;
   } 
 
   if ( fOwnsSubstruc ) {
 
     for ( AlgMapConstIter iter = fOwnedSubAlgMp -> begin() ;
           iter != fOwnedSubAlgMp -> end() ; ++iter ) {
 
       Algorithm * subalg = iter -> second ;
 
       LOG("Algorithm", pDEBUG) << "Appending config from " << iter -> first << " -> " << subalg -> Id() ;
       const Registry & r = subalg->GetConfig();
       RgKey prefix = iter -> first + "/";
       fConfig -> Append(r,prefix);
 
     }
 
   } //if owned substructure
 
   return * fConfig ;
 }

Registry * Algorithm::GetOwnedConfig ( void )

inherited

Returns the pointer of the summary registry, see previous method Gives access to the summary so it could be changed. The usage of this method is deprecated as it is mantained only for back compatibility. If you need to add or chage a parter (or more), use the AddTopRegistry() instead

Definition at line 287 of file Algorithm.cxx.

References GetConfig().

Referenced by genie::TransverseEnhancementFFModel::LoadConfig(), and genie::EffectiveSF::LoadConfig().

 {
   
   GetConfig() ;
   return fConfig;
 }

template<class T >

bool genie::Algorithm::GetParam	(	const RgKey &	name,
		T &	p,
		bool	is_top_call = `true`
	)		const

protectedinherited

Ideal access to a parameter value from the vector of registries Returns true if the value is found and the parameters is set

template<class T >

bool genie::Algorithm::GetParamDef	(	const RgKey &	name,
		T &	p,
		const T &	def
	)		const

protectedinherited

Ideal access to a parameter value from the vector of registries, With default value. Returns true if the value is set from the registries, false if the value is the default

template<class T >

bool genie::Algorithm::GetParamVect	(	const std::string &	comm_name,
		std::vector< T > &	v,
		unsigned int	max,
		bool	is_top_call = `true`
	)		const

protectedinherited

Handle to load vectors of parameters It looks for different registry item with name comm_name0, comm_name1, etc...

virtual AlgStatus_t genie::Algorithm::GetStatus ( void ) const

inlinevirtualinherited

Get algorithm status.

Definition at line 101 of file Algorithm.h.

References genie::Algorithm::fStatus.

101 { return fStatus; }

genie::Algorithm::fStatus

AlgStatus_t fStatus

algorithm execution status

Definition: Algorithm.h:166

virtual const AlgId& genie::Algorithm::Id ( void ) const

inlinevirtualinherited

Get algorithm ID.

Definition at line 98 of file Algorithm.h.

References genie::Algorithm::fID.

Referenced by genie::KineGeneratorWithCache::AccessCacheBranch(), genie::QELEventGeneratorSM::AccessCacheBranch2(), genie::QELEventGeneratorSM::AccessCacheBranchDiffv(), genie::InteractionListAssembler::AssembleInteractionList(), genie::XSecAlgorithmMap::BuildMap(), genie::InteractionGeneratorMap::BuildMap(), genie::XSecSplineList::BuildSplineKey(), genie::DISXSec::CacheBranchName(), genie::ReinSehgalRESXSecWithCache::CacheBranchName(), genie::DMDISXSec::CacheBranchName(), genie::ReinSehgalRESXSecWithCacheFast::CacheBranchName(), genie::Algorithm::Compare(), genie::RESKinematicsGenerator::ComputeMaxXSec(), genie::COHElKinematicsGenerator::ComputeMaxXSec(), genie::SKKinematicsGenerator::ComputeMaxXSec(), genie::COHKinematicsGenerator::ComputeMaxXSec(), genie::Algorithm::Configure(), genie::GEVGDriver::CreateSplines(), genie::QPMDISPXSec::DISRESJoinSuppressionFactor(), genie::QPMDMDISPXSec::DMDISRESJoinSuppressionFactor(), genie::AlgConfigPool::FindRegistry(), genie::AlgFactory::ForceReconfiguration(), genie::GEVGDriver::GenerateEvent(), GetAlgorithms(), genie::LwlynSmithQELCCPXSec::Integral(), Integral(), genie::COHXSec::Integrate(), genie::QPMDISPXSec::LoadConfig(), genie::QPMDMDISPXSec::LoadConfig(), genie::EventGenerator::LoadConfig(), genie::EventGeneratorListAssembler::LoadGenerator(), main(), genie::COHKinematicsGenerator::MaxXSec_AlvarezRuso(), genie::XSecAlgorithmMap::Print(), genie::InteractionGeneratorMap::Print(), genie::AlgFactory::Print(), genie::COHHadronicSystemGenerator::ProcessEventRecord(), genie::COHPrimaryLeptonGenerator::ProcessEventRecord(), genie::COHKinematicsGenerator::ProcessEventRecord(), genie::MECGenerator::ProcessEventRecord(), genie::EventGenerator::ProcessEventRecord(), genie::KNOPythiaHadronization::SelectHadronizer(), TestPythiaTauDecays(), and genie::GEVGDriver::UseSplines().

98 { return fID; }

genie::Algorithm::fID

AlgId fID

algorithm name and configuration set

Definition: Algorithm.h:156

void Algorithm::Initialize ( void )

protectedinherited

Definition at line 343 of file Algorithm.cxx.

Referenced by genie::Algorithm::AllowReconfig().

 {
 // Algorithm initialization
 //
   fAllowReconfig  = true;
   fOwnsSubstruc   = false;
   fConfig         = 0;
   fOwnedSubAlgMp  = 0;
 }

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

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

 {
   // 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 865 of file NievesQELCCPXSec.cxx.

References MECModelEnuComparisons::i, calib::j, and APDHVSetting::temp.

Referenced by LmunuAnumu().

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

References genie::Target::A(), genie::units::A, a, demo::app, b, CNCTCLimUcalc(), delta, genie::QELFormFactors::F1V(), genie::QELFormFactors::FA(), fCompareNievesTensors, fCoulombFactor, fFormFactors, genie::QELFormFactors::Fp(), fRPA, fTensorsOutFile, fVc, MECModelEnuComparisons::g, genie::Target::HitNucPdg(), genie::Target::HitNucPosition(), MECModelEnuComparisons::i, genie::Target::IsNucleus(), genie::pdg::IsProton(), calib::j, genie::controls::kASmallNum, leviCivita(), LOG, genie::Target::N(), pDEBUG, genie::Target::Pdg(), pWARN, q2, r(), canMan::sign(), ana::Sqrt(), sum, getGoodRuns4SAM::t0, genie::QELFormFactors::xiF2V(), Z, and genie::Target::Z().

Referenced by XSec().

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

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

Referenced by Configure().

 {
   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 );
 }

int Algorithm::MergeTopRegistry ( const Registry & r )

protectedinherited

Merge with top level registry if first reg of the vector is owned Otherwise an owned copy is added as a top registry

Definition at line 618 of file Algorithm.cxx.

                                                       {
 
   if ( fOwnerships.empty() ) {
 
         // this algorithm is not configured right now, the incoming registry is the only configuration
     Registry * p = new Registry( r ) ;
         AddTopRegistry( p ) ;
 
         return 1 ;
   }
 
   if ( fOwnerships[0] ) {
         //the top registry is owned: it can be changed with no consequences for other algorithms
     fConfVect[0] -> Merge( r ) ;
   }
   else {
         // The top registry is not owned so it cannot be changed
         // The registry will be added with top priority
 
         Registry * p = new Registry( r ) ;
         AddTopRegistry( p ) ;
   }
 
   // The configuration has changed so the summary is not updated anymore and must be deleted
   if ( fConfig ) {
      delete fConfig ;
      fConfig = 0 ;
   }
 
   return fConfVect.size() ;
 }

void Algorithm::Print ( ostream & stream ) const

virtualinherited

Print algorithm info.

Definition at line 323 of file Algorithm.cxx.

References GetConfig(), and r().

Referenced by genie::Algorithm::AllowReconfig(), and genie::operator<<().

 {
   // print algorithm name & parameter-set
   stream << "\nAlgorithm Key: "  << this->fID.Key();
   stream << " - Owns Substruc: " << ((fOwnsSubstruc) ? "[true]" : "[false]");
 
   // print algorithm configuration
   const Registry & r = this->GetConfig();
   stream << r;
 
   if(fOwnsSubstruc) {
     AlgMapConstIter iter = fOwnedSubAlgMp->begin();
     for(; iter!=fOwnedSubAlgMp->end(); ++iter) {
       Algorithm * alg = iter->second;
       stream << "<Next algorithm is owned by : " << this->fID.Key() << ">";
       stream << *alg;
     }
   }
 }

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

private

Definition at line 637 of file NievesQELCCPXSec.cxx.

References fhbarc, LOG, m, pWARN, relLindhardIm(), ruLinRelX(), and getGoodRuns4SAM::t0.

Referenced by CNCTCLimUcalc().

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

References a, kPi, q2, fillBadChanDBTables::result, and ana::Sqrt().

Referenced by CNCTCLimUcalc(), and relLindhard().

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

References kf2, genie::constants::kPi2, test_ParserArtEvents::log, Munits::m2, q2, and ana::Sqrt().

Referenced by relLindhard().

 {
   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;
 }

void Algorithm::SetId ( const AlgId & id )

virtualinherited

Set algorithm ID.

Definition at line 313 of file Algorithm.cxx.

Referenced by genie::Algorithm::AllowReconfig().

 {
   fID.Copy(id);
 }

void Algorithm::SetId	(	string	name,
		string	config
	)

virtualinherited

Definition at line 318 of file Algorithm.cxx.

 {
   fID.SetId(name, config);
 }

const Algorithm * Algorithm::SubAlg ( const RgKey & registry_key ) const

inherited

Access the sub-algorithm pointed to by the input key, either from the local pool or from AlgFactory's pool

Definition at line 353 of file Algorithm.cxx.

References ana::assert(), genie::AlgFactory::GetAlgorithm(), genie::AlgFactory::Instance(), LOG, pERROR, and pINFO.

 {
 // Returns the sub-algorithm pointed to this algorithm's XML config file using
 // the the values of the key.
 // This method asserts the existence of these keys in the XML config.
 // Note: Since only 1 parameter is used, the key value should contain both the
 // algorithm name and its configuration set according to the usual scheme:
 // namespace::algorithm_name/configuration_set
 //
   LOG("Algorithm", pINFO)
     << "Fetching sub-alg within alg: " << this->Id().Key()
                                    << " pointed to by key: " << registry_key;
   
   //-- if the algorithm owns its substructure:
   //      return the sub-algorithm from the local pool
   //
   if(fOwnsSubstruc) {
     AlgMapConstIter iter = fOwnedSubAlgMp->find(registry_key);
     if(iter!=fOwnedSubAlgMp->end()) return iter->second;
     LOG("Algorithm", pERROR) 
       << "Owned sub-alg pointed to by key: " << registry_key 
       << " was not found within alg: " << this->Id().Key();
      return 0;
   }
 
   //-- if the algorithm does not own its substructure:
   //      return the sub-algorithm from the AlgFactory's pool
   RgAlg alg ;
   GetParam( registry_key, alg ) ;
 
   LOG("Algorithm", pINFO)
     << "Registry key: " << registry_key << " points to algorithm: " << alg;
   
   // retrieve the Algorithm object from the the Algorithm factory
   AlgFactory * algf = AlgFactory::Instance();
   const Algorithm * algbase = algf->GetAlgorithm(alg.name, alg.config);
   assert(algbase);
 
   return algbase;
 }

bool XSecAlgorithmI::ValidKinematics ( const Interaction * i ) const

virtualinherited

Is the input kinematical point a physically allowed one?

Reimplemented in genie::KovalenkoQELCharmPXSec, genie::PaisQELLambdaPXSec, genie::NuElectronPXSec, genie::KLVOxygenIBDPXSec, genie::IMDAnnihilationPXSec, genie::StrumiaVissaniIBDPXSec, genie::H3AMNuGammaPXSec, genie::GLRESPXSec, and genie::IBDXSecMap.

Definition at line 46 of file XSecAlgorithmI.cxx.

References genie::KPhaseSpace::IsAboveThreshold(), genie::KPhaseSpace::IsAllowed(), genie::kISkipKinematicChk, LOG, genie::Interaction::PhaseSpace(), and pINFO.

 {
 // can offer common implementation for all concrete x-section models because
 // the input interaction is aware of its kinematic limits
 
   if ( interaction->TestBit(kISkipKinematicChk) ) return true;
 
   const KPhaseSpace& kps = interaction->PhaseSpace();
 
   if ( ! kps.IsAboveThreshold() ) {
      LOG("XSecBase", pINFO)  << "*** Below energy threshold";
      return false;
   }
   if ( ! kps.IsAllowed() ) {
      LOG("XSecBase", pINFO)  << "*** Not in allowed kinematical space";
      return false;
   }
   return true;
 }

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

virtual

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 360 of file NievesQELCCPXSec.cxx.

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

Referenced by XSec().

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

References genie::Target::A(), genie::units::A, plot_validation_datamc::c, E, exit(), fCoulombRmaxMode, fhbarc, fR0, func(), genie::gAbortingInErr, genie::utils::gsl::Integration1DimTypeFromString(), genie::Target::IsNucleus(), genie::constants::kAem, genie::kMatchNieves, genie::kMatchVertexGeneratorRmax, kPi, LOG, pFATAL, pNOTICE, cet::pow(), fillBadChanDBTables::result, ana::Sqrt(), test::z, Z, and genie::Target::Z().

Referenced by XSec().

                                                                      {
   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
     return -kAem*4*kPi*result*fhbarc;
   }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 74 of file NievesQELCCPXSec.cxx.

References std::abs(), angle, genie::KinePhaseSpace::AsString(), genie::QELFormFactors::Calculate(), genie::utils::EnergyDeltaFunctionSolutionQEL(), fCos8c2, fCoulomb, fDoPauliBlocking, fFormFactors, fPauliBlocker, fRPA, genie::Kinematics::FSLeptonP4(), genie::Interaction::FSPrimLepton(), fXSecScale, genie::PauliBlocker::GetFermiMomentum(), genie::InitialState::GetProbeP4(), genie::Kinematics::HadSystP4(), genie::Target::HitNucMass(), genie::Target::HitNucP4(), genie::Target::HitNucPdg(), genie::Target::HitNucPosition(), genie::pdg::IsNeutrino(), genie::Target::IsNucleus(), genie::pdg::IsProton(), genie::utils::mec::J(), genie::utils::kinematics::Jacobian(), genie::controls::kASmallNum, genie::constants::kGF2, genie::kIAssumeFreeNucleon, kinematics(), genie::Interaction::KinePtr(), kPi, genie::kPSQELEvGen, genie::kRfLab, LmunuAnumu(), LOG, genie::Target::N(), pDEBUG, std_candles::pl, cet::pow(), genie::InitialState::ProbePdg(), pWARN, genie::utils::kinematics::Q2(), q2, r(), genie::Interaction::RecoilNucleon(), genie::Interaction::RecoilNucleonPdg(), genie::Kinematics::SetQ2(), canMan::sign(), std::sqrt(), ana::Sqrt(), genie::InitialState::Tgt(), Unit(), genie::XSecAlgorithmI::ValidKinematics(), ValidProcess(), vcr(), xsec, and genie::Target::Z().

 {
   /*// 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();
 
   // 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::Algorithm::fAllowReconfig

protectedinherited

Definition at line 153 of file Algorithm.h.

Referenced by genie::Algorithm::AllowReconfig(), genie::Decayer::Configure(), genie::UnstableParticleDecayer::Configure(), genie::LHAPDF5::Configure(), and genie::LHAPDF6::Configure().

bool genie::NievesQELCCPXSec::fCompareNievesTensors

mutableprivate

print tensors

Definition at line 150 of file NievesQELCCPXSec.h.

Referenced by LmunuAnumu(), and LoadConfig().

vector<Registry*> genie::Algorithm::fConfVect

protectedinherited

ideally these members should go private Registry will be access only through the GetParam method configurations registries from various sources the order of the vector is the precedence in case of repeated parameters position 0 -> Highest precedence

Definition at line 161 of file Algorithm.h.

double genie::NievesQELCCPXSec::fCos8c2

private

cos^2(cabibbo angle)

Definition at line 71 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

bool genie::NievesQELCCPXSec::fCoulomb

private

use Coulomb corrections

Definition at line 79 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::NievesQELCCPXSec::fCoulombFactor

mutableprivate

Definition at line 152 of file NievesQELCCPXSec.h.

Referenced by LmunuAnumu().

Nieves_Coulomb_Rmax_t genie::NievesQELCCPXSec::fCoulombRmaxMode

private

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

Definition at line 108 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and vcr().

bool genie::NievesQELCCPXSec::fDoPauliBlocking

private

Whether to apply Pauli blocking in XSec()

Definition at line 97 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::NievesQELCCPXSec::fEnergyCutOff

private

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

Definition at line 94 of file NievesQELCCPXSec.h.

Referenced by LoadConfig().

QELFormFactors genie::NievesQELCCPXSec::fFormFactors

mutableprivate

Definition at line 68 of file NievesQELCCPXSec.h.

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

const QELFormFactorsModelI* genie::NievesQELCCPXSec::fFormFactorsModel

private

Definition at line 69 of file NievesQELCCPXSec.h.

Referenced by LoadConfig().

double genie::NievesQELCCPXSec::fhbarc

private

hbar*c in GeV*fm

Definition at line 75 of file NievesQELCCPXSec.h.

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

AlgId genie::Algorithm::fID

protectedinherited

algorithm name and configuration set

Definition at line 156 of file Algorithm.h.

Referenced by genie::Algorithm::Id().

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

Referenced by LoadConfig().

const FermiMomentumTable* genie::NievesQELCCPXSec::fKFTable

private

Definition at line 85 of file NievesQELCCPXSec.h.

Referenced by CNCTCLimUcalc(), and LoadConfig().

string genie::NievesQELCCPXSec::fKFTableName

private

Definition at line 86 of file NievesQELCCPXSec.h.

Referenced by LoadConfig().

bool genie::NievesQELCCPXSec::fLFG

private

Definition at line 84 of file NievesQELCCPXSec.h.

Referenced by CNCTCLimUcalc(), and LoadConfig().

const NuclearModelI* genie::NievesQELCCPXSec::fNuclModel

private

Nuclear Model for integration.

Definition at line 81 of file NievesQELCCPXSec.h.

Referenced by LoadConfig().

AlgMap* genie::Algorithm::fOwnedSubAlgMp

protectedinherited

local pool for owned sub-algs (taken out of the factory pool)

Definition at line 167 of file Algorithm.h.

vector<bool> genie::Algorithm::fOwnerships

protectedinherited

ownership for every registry in fConfVect

Definition at line 164 of file Algorithm.h.

bool genie::Algorithm::fOwnsSubstruc

protectedinherited

true if it owns its substructure (sub-algs,...)

Definition at line 155 of file Algorithm.h.

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

private

The PauliBlocker instance to use to apply that correction.

Definition at line 99 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::NievesQELCCPXSec::fR0

private

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

Definition at line 104 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and vcr().

bool genie::NievesQELCCPXSec::fRPA

mutableprivate

use RPA corrections

Definition at line 78 of file NievesQELCCPXSec.h.

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

AlgStatus_t genie::Algorithm::fStatus

protectedinherited

algorithm execution status

Definition at line 166 of file Algorithm.h.

Referenced by genie::Algorithm::GetStatus().

TString genie::NievesQELCCPXSec::fTensorsOutFile

mutableprivate

file to print tensors to

Definition at line 151 of file NievesQELCCPXSec.h.

Referenced by LmunuAnumu().

double genie::NievesQELCCPXSec::fVc

mutableprivate

Definition at line 152 of file NievesQELCCPXSec.h.

Referenced by LmunuAnumu().

const XSecIntegratorI* genie::NievesQELCCPXSec::fXSecIntegrator

private

Definition at line 70 of file NievesQELCCPXSec.h.

Referenced by Integral(), and LoadConfig().

double genie::NievesQELCCPXSec::fXSecScale

private

external xsec scaling factor

Definition at line 73 of file NievesQELCCPXSec.h.

Referenced by LoadConfig(), and XSec().

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

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

Public Member Functions

Protected Member Functions

Protected Attributes

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation