Contains auxiliary functions for Smith-Moniz model.
. More...

#include <SmithMonizUtils.h>

Inheritance diagram for genie::SmithMonizUtils:

Collaboration diagram for genie::SmithMonizUtils:

Classes
class	Func1D

class	Func1Dpar

class	Functor1D

Public Member Functions
	SmithMonizUtils ()

	SmithMonizUtils (string config)

virtual	~SmithMonizUtils ()

void	SetInteraction (const Interaction *i)

double	GetBindingEnergy (void) const

double	GetFermiMomentum (void) const

double	GetTheta_k (double v, double qv) const

double	GetTheta_p (double pv, double v, double qv, double &E_p) const

double	E_nu_thr_SM (void) const

Range1D_t	Q2QES_SM_lim (void) const

Range1D_t	vQES_SM_lim (double Q2) const

Range1D_t	kFQES_SM_lim (double nu, double Q2) const

double	PhaseSpaceVolume (KinePhaseSpace_t ps) const

double	vQES_SM_min (double Q2min, double Q2max) const

double	vQES_SM_max (double Q2min, double Q2max) const

void	Configure (const Registry &config)

void	Configure (string config)

double	QEL_EnuMin_SM (double E_nu) const

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

Static Public Member Functions
static double	rho (double P_Fermi, double T_Fermi, double p)

Static Public Member Functions inherited from genie::Algorithm
static string	BuildParamVectKey (const std::string &comm_name, unsigned int i)

static string	BuildParamVectSizeKey (const std::string &comm_name)

static string	BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)

static string	BuildParamMatRowSizeKey (const std::string &comm_name)

static string	BuildParamMatColSizeKey (const std::string &comm_name)

Private Member Functions
void	LoadConfig (void)

double	Q2lim1_SM (double Q2, double Enu) const

double	Q2lim2_SM (double Q2) const

void	DMINFC (Functor1D &F, double A, double B, double EPS, double DELTA, double &X, double &Y, bool &LLM) const

double	vlim1_SM (double Q2) const

double	vlim2_SM (double Q2) const

Private Attributes
map< int, double >	fNucRmvE

string	fKFTable

bool	fUseParametrization

const Interaction *	fInteraction

double	E_nu
	Neutrino energy (GeV) More...

double	m_lep
	Mass of final charged lepton (GeV) More...

double	mm_lep
	Squared mass of final charged lepton (GeV) More...

double	m_ini
	Mass of initial hadron or hadron system (GeV) More...

double	mm_ini
	Sqared mass of initial hadron or hadron system (GeV) More...

double	m_fin
	Mass of final hadron or hadron system (GeV) More...

double	mm_fin
	Squared mass of final hadron or hadron system (GeV) More...

double	m_tar
	Mass of target nucleus (GeV) More...

double	mm_tar
	Squared mass of target nucleus (GeV) More...

double	m_rnu
	Mass of residual nucleus (GeV) More...

double	mm_rnu
	Squared mass of residual nucleus (GeV) More...

double	P_Fermi
	Maximum value of Fermi momentum of target nucleon (GeV) More...

double	E_BIN
	Binding energy (GeV) More...

Additional Inherited Members
Protected Member Functions inherited from genie::Algorithm
	Algorithm ()

	Algorithm (string name)

	Algorithm (string name, string config)

void	Initialize (void)

void	DeleteConfig (void)

void	DeleteSubstructure (void)

Registry *	ExtractLocalConfig (const Registry &in) const

Registry *	ExtractLowerConfig (const Registry &in, const string &alg_key) const
	Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key. More...

template<class T >
bool	GetParam (const RgKey &name, T &p, bool is_top_call=true) const

template<class T >
bool	GetParamDef (const RgKey &name, T &p, const T &def) const

template<class T >
int	GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
	Handle to load vectors of parameters. More...

int	GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const

template<class T >
int	GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
	Handle to load matrix of parameters. More...

template<class T >
int	GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const

int	GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const

int	AddTopRegistry (Registry *rp, bool owns=true)
	add registry with top priority, also update ownership More...

int	AddLowRegistry (Registry *rp, bool owns=true)
	add registry with lowest priority, also update ownership More...

int	MergeTopRegistry (const Registry &r)

int	AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
	Add registries with top priority, also udated Ownerships. More...

Protected Attributes inherited from genie::Algorithm
bool	fAllowReconfig

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

AlgId	fID
	algorithm name and configuration set More...

vector< Registry * >	fConfVect

vector< bool >	fOwnerships
	ownership for every registry in fConfVect More...

AlgStatus_t	fStatus
	algorithm execution status More...

AlgMap *	fOwnedSubAlgMp
	local pool for owned sub-algs (taken out of the factory pool) More...

Detailed Description

Contains auxiliary functions for Smith-Moniz model.
.

References:: [1] R.A.Smith and E.J.Moniz, Nuclear Physics B43, (1972) 605-622
[2] K.S. Kuzmin, V.V. Lyubushkin, V.A.Naumov, Eur. Phys. J. C54, (2008) 517-538
[3] K.S.Kuzmin, V.A.Naumov, V.V.Lyubushkin, I.D.Kakorin, Kinematic formulas for GENIE implementation of Smith-Monitz model:
living document for internal use (http://theor.jinr.ru/NeutrinoOscillations/Papers/GENIE_formulas.pdf)

Author: Igor Kakorin kakor.nosp@m.in@j.nosp@m.inr.r.nosp@m.u Joint Institute for Nuclear Research
adapted from fortran code provided by:
Konstantin Kuzmin kkuzm.nosp@m.in@t.nosp@m.heor..nosp@m.jinr.nosp@m..ru Joint Institute for Nuclear Research
Vladimir Lyubushkin Joint Institute for Nuclear Research
Vadim Naumov vnaum.nosp@m.ov@t.nosp@m.heor..nosp@m.jinr.nosp@m..ru Joint Institute for Nuclear Research
based on code of:
Costas Andreopoulos const.nosp@m.anti.nosp@m.nos.a.nosp@m.ndre.nosp@m.opoul.nosp@m.os@c.nosp@m.ern.c.nosp@m.h University of Liverpool

Created:: May 05, 2017

License:: Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 53 of file SmithMonizUtils.h.

Constructor & Destructor Documentation

SmithMonizUtils::SmithMonizUtils ( )

Definition at line 46 of file SmithMonizUtils.cxx.

                                  :
 Algorithm("genie::SmithMonizUtils")
 {
 
 }

SmithMonizUtils::SmithMonizUtils ( string config )

Definition at line 52 of file SmithMonizUtils.cxx.

                                               :
 Algorithm("genie::SmithMonizUtils", config)
 {
 
 }

SmithMonizUtils::~SmithMonizUtils ( )

virtual

Definition at line 58 of file SmithMonizUtils.cxx.

 {
 
 }

Member Function Documentation

void SmithMonizUtils::Configure ( const Registry & config )

virtual

methods overloading the Algorithm() interface implementation to build the fragmentation function from configuration data

Reimplemented from genie::Algorithm.

Definition at line 63 of file SmithMonizUtils.cxx.

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

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

void SmithMonizUtils::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 69 of file SmithMonizUtils.cxx.

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

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

void SmithMonizUtils::DMINFC	(	Functor1D &	F,
		double	A,
		double	B,
		double	EPS,
		double	DELTA,
		double &	X,
		double &	Y,
		bool &	LLM
	)		const

private

Definition at line 548 of file SmithMonizUtils.cxx.

References genie::units::A, and genie::utils::kinematics::W().

Referenced by Q2QES_SM_lim(), QEL_EnuMin_SM(), vQES_SM_max(), and vQES_SM_min().

 {
   const double W5 = TMath::Sqrt(5);
   const double HV = (3-W5)/2;
   const double HW = (W5-1)/2;
   const double R1 = 1.0;
   const double HF = R1/2;
 
   int N = -1;
   if(A!=B) N = TMath::Nint(2.08*TMath::Log(TMath::Abs((A-B)/EPS)));
   double C = A;
   double D = B;
   if(A>B)
   {
      C = B;
      D = A;
   }
   bool LLT = true;
   bool LGE = true;
   double V, FV, W, FW, H;
   while(true)
   {
      H = D-C;
      if(N<0)
      {
        X = HF*(C+D);
        Y = F(X);
        LLM = TMath::Abs(X-A)>DELTA && TMath::Abs(X-B)>DELTA;
        return;
      }
      if(LLT)
      {
        V = C+HV*H;
        FV = F(V);
      }
      if(LGE)
      {
         W = C+HW*H;
         FW = F(W);
      }
      if(FV<FW)
      {
         LLT = true;
         LGE = false;
         D = W;
         W = V;
         FW = FV;
      }
      else
      {
         LLT = false;
         LGE = true;
         C = V;
         V = W;
         FV = FW;
      }
      N = N-1;
   }
 }

double SmithMonizUtils::E_nu_thr_SM ( void ) const

Definition at line 210 of file SmithMonizUtils.cxx.

References E_BIN, m_fin, m_lep, m_rnu, m_tar, mm_ini, mm_tar, P_Fermi, and QEL_EnuMin_SM().

Referenced by genie::SmithMonizQELCCXSec::Integrate().

 {
 
   Func1D<SmithMonizUtils> QEL_EnuMin_SM_(*this, &SmithMonizUtils::QEL_EnuMin_SM);
 
   // maximum of function call
   const int MFC = 10000;
   const double EPSABS = 0.;
   const double EPSREL = 1.0e-08;
   
   // Energy threshold of scattering on nucleus (Eq. (1) of Ref. 3)
   double E_min = ((m_lep + m_rnu + m_fin)*(m_lep + m_rnu + m_fin) - mm_tar)/2/m_tar;
   
   // Energy threshold of scattering on bound nucleon (Eq. (2) of Ref. 3)
   double E_min2 = ((m_lep + m_fin)*(m_lep + m_fin)-mm_ini-E_BIN*E_BIN+2*E_BIN*TMath::Sqrt(mm_ini+P_Fermi*P_Fermi))/2/(TMath::Sqrt(mm_ini+P_Fermi*P_Fermi)-E_BIN+P_Fermi);
   
   E_min = TMath::Max(E_min, E_min2);
   
   // if nu_1>nu_max at E_min then we try to find energy in range (E_min, 50.) where nu_1=nu_max and put E_min equal to it.
   // nu_1   -- minimal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3), 
   // nu_max -- maximal energy transfer on nucleus (see Eq. (9) of Ref.3)
   if (QEL_EnuMin_SM(E_min) > 0)
   {
     // C++ analog of fortran function Enu_rf= DZEROX(E_min,Enu_2,EPS,MFC,QEL_EnuMin_SM,1)
     ROOT::Math::RootFinder rfgb(ROOT::Math::RootFinder::kGSL_BRENT);
     //convergence is reached using tolerance = 2 *( epsrel * abs(x) + epsabs)
     if ( rfgb.Solve(QEL_EnuMin_SM_, E_min, 50., MFC, EPSABS, EPSREL) )
     {
        E_min = rfgb.Root();
     }
   }
   E_min = TMath::Max(E_min, 0.);
   return E_min;
 
 }

double SmithMonizUtils::GetBindingEnergy ( void ) const

Definition at line 629 of file SmithMonizUtils.cxx.

References E_BIN.

Referenced by genie::SmithMonizQELCCPXSec::d3sQES_dQ2dvdkF_SM(), and genie::QELEventGeneratorSM::ProcessEventRecord().

 {
   return E_BIN;
 }

double SmithMonizUtils::GetFermiMomentum ( void ) const

Definition at line 634 of file SmithMonizUtils.cxx.

References P_Fermi.

Referenced by genie::QELEventGeneratorSM::ComputeMaxXSec(), genie::SmithMonizQELCCPXSec::d3sQES_dQ2dvdkF_SM(), and genie::QELEventGeneratorSM::ProcessEventRecord().

 {
   return P_Fermi;
 }

double SmithMonizUtils::GetTheta_k	(	double	v,
		double	qv
	)		const

Definition at line 639 of file SmithMonizUtils.cxx.

References E_nu, and mm_lep.

 {
   return TMath::ACos((v + (mm_lep-v*v+qv*qv)/2/E_nu)/qv);
 }

double SmithMonizUtils::GetTheta_p	(	double	pv,
		double	v,
		double	qv,
		double &	E_p
	)		const

Definition at line 644 of file SmithMonizUtils.cxx.

References E_BIN, mm_fin, and mm_ini.

 {
   E_p = TMath::Sqrt(mm_ini+pv*pv)-E_BIN;
   if (pv != 0)
     return TMath::ACos(((v-E_BIN)*(2*E_p+v+E_BIN)-qv*qv+mm_ini-mm_fin)/(2*pv*qv));
   else
     return 0;
 }

Range1D_t SmithMonizUtils::kFQES_SM_lim	(	double	nu,
		double	Q2
	)		const

Definition at line 529 of file SmithMonizUtils.cxx.

References E_BIN, m_ini, mm_fin, mm_ini, and P_Fermi.

Referenced by genie::SmithMonizQELCCPXSec::d2sQES_dQ2dv_SM(), PhaseSpaceVolume(), and genie::QELEventGeneratorSM::ProcessEventRecord().

 {
   double qv = TMath::Sqrt(nu*nu+Q2);
   double c_f = (nu-E_BIN)/qv;
   double d_f = (E_BIN*E_BIN-2*nu*E_BIN-Q2+mm_ini-mm_fin)/(2*qv*m_ini);
   // minimal energy of initial nucleon (see Eq. (13) of Ref.3)
   double Ef_min= TMath::Max(m_ini*(c_f*d_f+TMath::Sqrt(1.0-c_f*c_f+d_f*d_f))/(1.0-c_f*c_f), TMath::Sqrt(P_Fermi*P_Fermi+mm_ini)-nu);
   double kF_min= P_Fermi!=0?TMath::Sqrt(TMath::Max(Ef_min*Ef_min-mm_ini,0.0)):0.;
   double kF_max= P_Fermi;
   Range1D_t R;
   if (kF_min<=kF_max)
     R = Range1D_t(kF_min,kF_max);
   else
     R = Range1D_t(0.5*(kF_min+kF_max),0.5*(kF_min+kF_max));
   return R;
 
 }

void SmithMonizUtils::LoadConfig ( void )

private

Definition at line 75 of file SmithMonizUtils.cxx.

References fKFTable, fNucRmvE, fUseParametrization, genie::Algorithm::GetConfig(), genie::Registry::GetItemMap(), genie::Algorithm::GetParam(), and genie::pdg::IonPdgCodeToZ().

Referenced by Configure().

 {
 
 
   GetParam( "FermiMomentumTable", fKFTable);
   GetParam( "RFG-UseParametrization", fUseParametrization);
 
 
   // load removal energy for specific nuclei from either the algorithm's
   // configuration file or the UserPhysicsOptions file.
   // if none is used use Wapstra's semi-empirical formula.
   const std::string keyStart = "RFG-NucRemovalE@Pdg=";
 
   RgIMap entries = GetConfig().GetItemMap();
 
   for(RgIMap::const_iterator it = entries.begin(); it != entries.end(); ++it)
   {
     const std::string& key = it->first;
     int pdg = 0;
     int Z = 0;
     if (0 == key.compare(0, keyStart.size(), keyStart.c_str()))
     {
       pdg = atoi(key.c_str() + keyStart.size());
       Z = pdg::IonPdgCodeToZ(pdg);
     }
     if (0 != pdg && 0 != Z)
     {
       ostringstream key_ss ;
       key_ss << keyStart << pdg;
       RgKey rgkey   = key_ss.str();
       double eb;
       GetParam( rgkey, eb ) ;
       eb = TMath::Max(eb, 0.);
       fNucRmvE.insert(map<int,double>::value_type(Z,eb));
     }
   }
 
 
 }

double SmithMonizUtils::PhaseSpaceVolume ( KinePhaseSpace_t ps ) const

Definition at line 653 of file SmithMonizUtils.cxx.

References kFQES_SM_lim(), genie::kPSQ2fE, genie::kPSQ2vpfE, genie::Range1D_t::max, genie::Range1D_t::min, genie::utils::kinematics::Q2(), Q2QES_SM_lim(), and vQES_SM_lim().

Referenced by genie::QELEventGeneratorSM::ProcessEventRecord().

 {
    double vol = 0.0;
    if (ps == kPSQ2fE)
    {
      Range1D_t rQ2 = Q2QES_SM_lim();
      vol = rQ2.max - rQ2.min;
      return vol;
    }
    else if (ps == kPSQ2vpfE)
    {
      const int kNQ2 = 101;
      const int kNv  = 101;
      Range1D_t rQ2 = Q2QES_SM_lim();
      double dQ2 = (rQ2.max-rQ2.min)/(kNQ2-1);
      for(int iq2=0; iq2<kNQ2-1; iq2++)
      {
        double Q2 = rQ2.min + iq2*dQ2;
        Range1D_t rv  = vQES_SM_lim(Q2);
        double dv = (rv.max-rv.min)/(kNv-1);
        for(int iv=0; iv<kNv-1; iv++)
        {
          double v = rv.min + iv*dv;
          Range1D_t rkF  = kFQES_SM_lim(Q2,v);
          double dkF = (rkF.max-rkF.min);
          vol += (dQ2*dv*dkF);
        }
      }
      return vol;
    }
    else
      return 0;
 }

double SmithMonizUtils::Q2lim1_SM	(	double	Q2,
		double	Enu
	)		const

private

Definition at line 273 of file SmithMonizUtils.cxx.

References a, genie::units::b, E_BIN, mm_fin, mm_ini, mm_lep, and P_Fermi.

Referenced by Q2QES_SM_lim(), and QEL_EnuMin_SM().

 {
   // maximal energy transfer (see Eq. (9) of Ref.3) 
   double nu_max = Enu*Q2/(Q2+mm_lep)-(Q2+mm_lep)/4/Enu;
   
   double E = sqrt(P_Fermi*P_Fermi+mm_ini);
   double b = (E-E_BIN)*(E-E_BIN)-P_Fermi*P_Fermi;
   double a = 0.5*(Q2+mm_fin-b);
   // minimal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3), 
   double nu_1 = (a*(E-E_BIN)-P_Fermi*TMath::Sqrt(a*a+Q2*b))/b;
   return nu_1-nu_max;
 
 }

double SmithMonizUtils::Q2lim2_SM ( double Q2 ) const

private

Definition at line 288 of file SmithMonizUtils.cxx.

References a, genie::units::b, E_BIN, m_fin, m_rnu, m_tar, mm_fin, mm_ini, mm_tar, and P_Fermi.

Referenced by Q2QES_SM_lim().

 {
   // minimal energy transfer for scattering on nucleus (see Eq. (7) of Ref.3)
   double nu_min = ((m_rnu+m_fin)*(m_rnu+m_fin)+Q2-mm_tar)/(2*m_tar);
   
   double E = sqrt(P_Fermi*P_Fermi+mm_ini);
   double b = (E-E_BIN)*(E-E_BIN)-P_Fermi*P_Fermi;
   double a = (Q2+mm_fin-b)*0.5;
   // maximal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3)
   double nu_2  = (a*(E-E_BIN)+P_Fermi*TMath::Sqrt(a*a+Q2*b))/b;
   return nu_min-nu_2;
 
 }

Range1D_t SmithMonizUtils::Q2QES_SM_lim ( void ) const

Definition at line 303 of file SmithMonizUtils.cxx.

References DMINFC(), E_BIN, E_nu, epsilon, LOG, m_fin, m_ini, m_lep, m_rnu, m_tar, mm_fin, mm_ini, mm_lep, mm_tar, P_Fermi, pFATAL, pWARN, Q2lim1_SM(), Q2lim2_SM(), and genie::units::s.

Referenced by genie::QELEventGeneratorSM::ComputeMaxXSec(), genie::utils::gsl::d2Xsec_dQ2dv::d2Xsec_dQ2dv(), PhaseSpaceVolume(), and genie::QELEventGeneratorSM::ProcessEventRecord().

 {
 
   // here we find Q2_min and Q2_max such that
   // 0. Q2_min>=0
   // 1. nu_1(Q2_min)<=nu_max(Q2_min)
   // 2. nu_1(Q2_max)<=nu_max(Q2_max)
   // 3. nu_min(Q2_min)<=nu_2(Q2_min)
   // 4. Q2_min>=the value of minimal Q2 defined for scattering on nucleus 
   // 5. Q2_max<=the value of maximal Q2 defined for scattering on nucleus (see Eqs. (3) of Ref.3)
   // nu_1 --   minimal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3), 
   // nu_max -- maximal energy transfer on nucleus (see Eq. (9) of Ref.3),
   // nu_2 --   maximal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3), 
   // nu_min -- minimal energy transfer on nucleus (see Eq. (7) of Ref.3)
   
   // maximum of function call
   const int MFC = 1000;
   const double EPS  = 1.0e-08;
   const double Delta= 1.0e-14;
   const double EPSABS = 0;
   const double EPSREL = 1.0e-08;
   const double Precision = std::numeric_limits<double>::epsilon();
   // if the nucleus mass is less than 4 then this is a special case
   if (E_BIN == 0 && P_Fermi == 0)
   {
     double s = 2*E_nu*m_ini+mm_ini;
     // minimal W2 for scattering on nucleus (see Eq. (6) of Ref.3)
     double W2 = mm_fin;
     // neutrino energy in CMS (see Eq. (4) of Ref.3)
     double E_nu_CM = (s-mm_ini)/2/TMath::Sqrt(s);
     // final lepton energy and momentum at W2_min (see Eqs. (5) of Ref.3)
     double E_l_CM  = (s+mm_lep-W2)/2/TMath::Sqrt(s);
     double P_l_CM  = E_l_CM>m_lep?TMath::Sqrt(E_l_CM*E_l_CM-mm_lep):Precision;
     // minimal and maximal allowed Q2 for scattering on nucleus (see Eqs. (3) of Ref.3)
     double Q2_min = 2*E_nu_CM*(E_l_CM-P_l_CM)-mm_lep;
     double Q2_max = 2*E_nu_CM*(E_l_CM+P_l_CM)-mm_lep;
     Q2_min= TMath::Max(Q2_min,0.0);
     Range1D_t R(Q2_min,Q2_max);
     return R;
   }
   double s = 2*E_nu*m_tar+mm_tar;
   // minimal W2 for scattering on nucleus (see Eq. (6) of Ref.3)
   double W2 = (m_rnu+m_fin)*(m_rnu+m_fin);
   // neutrino energy in CMS (see Eq. (4) of Ref.3)
   double E_nu_CM = (s-mm_tar)/2/TMath::Sqrt(s);
   // final lepton energy and momentum at W2_min (see Eqs. (5) of Ref.3)
   double E_l_CM  = (s+mm_lep-W2)/2/TMath::Sqrt(s);
   double P_l_CM  = E_l_CM>m_lep?TMath::Sqrt(E_l_CM*E_l_CM-mm_lep):Precision;
   // minimal and maximal allowed Q2 for scattering on nucleus (see Eqs. (3) of Ref.3)
   double Q2_min = 2*E_nu_CM*(E_l_CM-P_l_CM)-mm_lep;
   double Q2_max = 2*E_nu_CM*(E_l_CM+P_l_CM)-mm_lep;
   double F_MIN, Q2_0;
   bool LLM;
   // C++ analog of fortran function DMINFC(Q2lim1_SM,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM)
   Func1Dpar<SmithMonizUtils> Q2lim1_SM_(*this, &SmithMonizUtils::Q2lim1_SM, E_nu);
   // if minimum of nu_1-nu_max>0 then exit with error, because it's impossible
   DMINFC(Q2lim1_SM_,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM);
   if (F_MIN>0)
   {
       LOG("SmithMoniz", pFATAL)
       << "No overlapped area for energy " << E_nu << "\n" <<
       "Q2_min=" << Q2_min << " Q2_max=" << Q2_max << "\n" <<
       "Q2_0=" << Q2_0 << " F_MIN=" << F_MIN;
       exit(1);
   }
   // at Q2_0 here we have: nu_1(Q2_0)-nu_max(Q2_0)<0 
   // if nu_1(Q2_min)-nu_max(Q2_min)>0 we find corrected Q2_min_cor>Q2_min where nu_1(Q2_min_cor)-nu_max(Q2_min_cor)=0
   // (it is always possible because of above conditions) 
   if (Q2lim1_SM(Q2_min, E_nu)>0)
   {
     //C++ analog of fortran function Q2_RF=DZEROX(Q2_min,Q2_0,EPS,MFC,Q2lim1_SM,1)
     ROOT::Math::RootFinder rfgb(ROOT::Math::RootFinder::kGSL_BRENT);
     // convergence is reached using  tolerance = 2 *( epsrel * abs(x) + epsabs)
     if (rfgb.Solve(Q2lim1_SM_, Q2_min, Q2_0, MFC, EPSABS, EPSREL))
     {
       Q2_min= rfgb.Root();
     }
   }
   // if nu_1(Q2_max)-nu_max(Q2_max)>0 we find Q2_max_cor<Q2_max where nu_1(Q2_max_cor)-nu_max(Q2_max_cor)=0
   // (it is always possible because of above conditions) 
   if(Q2lim1_SM(Q2_max, E_nu)>0)
   {
      // C++ analog of fortran function Q2_RF=DZEROX(Q2_0,Q2_max,Eps,MFC,Q2lim1_SM,1)
      ROOT::Math::RootFinder rfgb(ROOT::Math::RootFinder::kGSL_BRENT);
      //convergence is reached using  tolerance = 2 *( epsrel * abs(x) + epsabs)
      if (rfgb.Solve(Q2lim1_SM_, Q2_0, Q2_max, MFC, EPSABS, EPSREL))
      {
        Q2_max= rfgb.Root();
      }
   }
   Func1D<SmithMonizUtils> Q2lim2_SM_(*this, &SmithMonizUtils::Q2lim2_SM);
   // if nu_min(Q2_min)-nu_2(Q2_min)>0 and nu_min(Q2_max)-nu_2(Q2_max)>0 then set Q2_min=Q2_max (it makes xsec equal to zero).
   if (Q2lim2_SM(Q2_min)>0)
   {
      if(Q2lim2_SM(Q2_max)>0)
      {
            LOG("SmithMoniz", pWARN) << "The RFG model is not applicable! The cross section is set zero!";
            Q2_min = Q2_max;
      }
      // here we have nu_min(Q2_min)-nu_2(Q2_min)>0 and nu_min(Q2_max)-nu_2(Q2_max)<0 or vice versa
      // so we always can find Q2_min_cor where nu_min(Q2_min_cor)-nu_2(Q2_min_cor)=0
      else
      {
        // C++ analog of fortran function Q2_RF = DZEROX(Q2_min,Q2_max,Eps,MFC,Q2lim2_SM,1)
        ROOT::Math::RootFinder rfgb(ROOT::Math::RootFinder::kGSL_BRENT);
        // convergence is reached using  tolerance = 2 *( epsrel * abs(x) + epsabs)
        if (rfgb.Solve(Q2lim2_SM_, Q2_min,Q2_max, MFC, EPSABS, EPSREL))
        {
           Q2_min= rfgb.Root();
        }
      }
   }
   Q2_min = TMath::Max(Q2_min,0.0);
 
   Range1D_t R(Q2_min,Q2_max);
   return R;
 
 }

double SmithMonizUtils::QEL_EnuMin_SM ( double E_nu ) const

Definition at line 247 of file SmithMonizUtils.cxx.

References DMINFC(), epsilon, m_fin, m_lep, m_rnu, m_tar, mm_lep, mm_tar, Q2lim1_SM(), and genie::units::s.

Referenced by E_nu_thr_SM().

 {
   // return the minimum of nu_1-nu_max as function of Q2 in range ( Q2_lim1(Enu), Q2_lim2(Enu) )
   const double EPS  = 1.0e-06;
   const double Delta= 1.0e-14;
   const double Precision = std::numeric_limits<double>::epsilon();
   double s = 2*Enu*m_tar+mm_tar;
   double W2 = (m_rnu+m_fin)*(m_rnu+m_fin);
   // neutrino energy in CMS (see Eq. (4) of Ref.3)
   double E_nu_CM = (s-mm_tar)/2/TMath::Sqrt(s);
   // final lepton energy and momentum at W2_min (see Eqs. (5) and (6) of Ref.3)
   double E_l_CM  = (s+mm_lep-W2)/2/TMath::Sqrt(s);
   double P_l_CM  = E_l_CM>m_lep?TMath::Sqrt(E_l_CM*E_l_CM-mm_lep):Precision;
   // minimal and maximal allowed Q2 for scattering on nucleus (see Eqs. (3) of Ref.3)
   double Q2_lim1 = 2*E_nu_CM*(E_l_CM-P_l_CM)-mm_lep;
   double Q2_lim2 = 2*E_nu_CM*(E_l_CM+P_l_CM)-mm_lep;
   // C++ analog of fortran function DMINFC(Q2lim1_SM,Q2_lim1,Q2_lim2,EPS,Delta,Q2_0,F_MIN,LLM)
   Func1Dpar<SmithMonizUtils> Q2lim1_SM_(*this, &SmithMonizUtils::Q2lim1_SM, Enu);
   double Q2_0,F_MIN;
   bool LLM;
   // find minimum of nu_1-nu_max as function of Q2 in range (Q2_lim1,Q2_lim2)
   DMINFC(Q2lim1_SM_,Q2_lim1,Q2_lim2,EPS,Delta,Q2_0,F_MIN,LLM);
   return F_MIN;
 }

double SmithMonizUtils::rho	(	double	P_Fermi,
		double	T_Fermi,
		double	p
	)

static

Definition at line 609 of file SmithMonizUtils.cxx.

Referenced by genie::SmithMonizQELCCPXSec::d3sQES_dQ2dvdkF_SM().

 {
 
   if (T_Fermi==0)
   {
     //Pure Fermi gaz with T_Fermi=0
     if(p<=P_Fermi)
       return 1.0;
     else
       return 0.0;
   }
   else
   {
       //Fermi-Dirac distribution
       return 1.0/(1.0 + TMath::Exp(-(P_Fermi-p)/T_Fermi));
   }
 
 
 }

void SmithMonizUtils::SetInteraction ( const Interaction * i )

Definition at line 116 of file SmithMonizUtils.cxx.

References genie::Target::A(), genie::units::A, genie::utils::nuclear::BindEnergyPerNucleon(), genie::utils::nuclear::BindEnergyPerNucleonParametrization(), E_BIN, E_nu, genie::utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization(), genie::PDGLibrary::Find(), genie::FermiMomentumTable::FindClosestKF(), fInteraction, fKFTable, fNucRmvE, genie::Interaction::FSPrimLepton(), fUseParametrization, genie::FermiMomentumTablePool::GetTable(), genie::Target::HitNucIsSet(), genie::Target::HitNucMass(), genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::FermiMomentumTablePool::Instance(), genie::PDGLibrary::Instance(), genie::pdg::IonPdgCode(), genie::pdg::IsProton(), genie::kRfLab, LOG, m_fin, m_ini, m_lep, m_rnu, m_tar, genie::utils::res::Mass(), genie::Target::Mass(), mm_fin, mm_ini, mm_lep, mm_rnu, mm_tar, P_Fermi, pFATAL, genie::InitialState::ProbeE(), genie::pdg::SwitchProtonNeutron(), genie::InitialState::Tgt(), and genie::Target::Z().

Referenced by genie::SmithMonizQELCCPXSec::d2sQES_dQ2dv_SM(), genie::utils::gsl::d2Xsec_dQ2dv::d2Xsec_dQ2dv(), genie::SmithMonizQELCCPXSec::d3sQES_dQ2dvdkF_SM(), genie::utils::gsl::dv_dQ2_E::dv_dQ2_E(), genie::SmithMonizQELCCXSec::Integrate(), and genie::QELEventGeneratorSM::ProcessEventRecord().

 {
 
   fInteraction = interaction;
   // get kinematics & init-state parameters
   const InitialState & init_state = interaction -> InitState();
   const Target & target = init_state.Tgt();
   PDGLibrary * pdglib = PDGLibrary::Instance();
   
   // neutrino energy (GeV)
   E_nu = interaction->InitState().ProbeE(kRfLab); 
 
   assert(target.HitNucIsSet());
   // get lepton&nuclear masses (init & final state nucleus)
   
   // mass of final charged lepton (GeV)
   m_lep = interaction->FSPrimLepton()->Mass();
   mm_lep     = TMath::Power(m_lep,    2);
   int nucl_pdg_ini = target.HitNucPdg();
   m_ini  = target.HitNucMass();
   mm_ini = TMath::Power(m_ini,    2);
   int nucl_pdg_fin = genie::pdg::SwitchProtonNeutron(nucl_pdg_ini);
   TParticlePDG * nucl_fin = pdglib->Find( nucl_pdg_fin );
   // mass of final hadron or hadron system (GeV)
   m_fin  = nucl_fin -> Mass();
   mm_fin = TMath::Power(m_fin,    2);
   // mass of target nucleus (GeV)
   m_tar = target.Mass();
   mm_tar = TMath::Power(m_tar,    2);
 
   // RFG is not applied for A<4
   if (target.A()<4)
   {
     E_BIN = P_Fermi = m_rnu = mm_rnu = 0;
     return;
   }
 
   bool is_p = pdg::IsProton(nucl_pdg_ini);
   int Zi = target.Z();
   int Ai = target.A();
   int Zf = (is_p) ? Zi-1 : Zi;
   int Af = Ai-1;
   TParticlePDG * nucl_f = pdglib->Find( pdg::IonPdgCode(Af, Zf) );
   if(!nucl_f)
   {
     LOG("SmithMoniz", pFATAL)
     << "Unknwown nuclear target! No target with code: "
     << pdg::IonPdgCode(Af, Zf) << " in PDGLibrary!";
     exit(1);
   }
   // mass of residual nucleus (GeV)
   m_rnu = nucl_f -> Mass();
   mm_rnu = TMath::Power(m_rnu, 2);
 
   int Z = target.Z();
   int A = target.A();
   int N = A-Z;
 
 
   // maximum value of Fermi momentum of target nucleon (GeV)
   if (A < 6 || !fUseParametrization)
   {
      // look up the Fermi momentum for this Target
      FermiMomentumTablePool * kftp = FermiMomentumTablePool::Instance();
      const FermiMomentumTable * kft = kftp->GetTable(fKFTable);
      P_Fermi = kft->FindClosestKF(pdg::IonPdgCode(A, Z), nucl_pdg_ini);
   }
   else
   {
     // define the Fermi momentum for this Target
     //
     P_Fermi = utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization(target);
     // correct the Fermi momentum for the struck nucleon
     if(is_p) P_Fermi *= TMath::Power( 2.*Z/A, 1./3);
     else
            P_Fermi *= TMath::Power( 2.*N/A, 1./3);
   }
 
   // neutrino binding energy (GeV)
   if (target.A() < 6 || !fUseParametrization)
   {
      map<int,double>::const_iterator it = fNucRmvE.find(Z);
      if(it != fNucRmvE.end()) E_BIN = it->second;
        else E_BIN = utils::nuclear::BindEnergyPerNucleon(target);
   }
   else
     E_BIN = utils::nuclear::BindEnergyPerNucleonParametrization(target);
 
 
 
 
 }

double SmithMonizUtils::vlim1_SM ( double Q2 ) const

private

Definition at line 499 of file SmithMonizUtils.cxx.

References a, genie::units::b, E_BIN, m_fin, m_rnu, m_tar, mm_fin, mm_ini, mm_tar, and P_Fermi.

Referenced by vQES_SM_min().

 {
   // minimal energy transfer for scattering on nucleus (see Eq. (7) of Ref.3)
   double nu_min = ((m_rnu+m_fin)*(m_rnu+m_fin)+Q2-mm_tar)/(2*m_tar);
   
   double E = sqrt(P_Fermi*P_Fermi+mm_ini);
   double b = (E-E_BIN)*(E-E_BIN)-P_Fermi*P_Fermi;
   double a = (Q2+mm_fin-b)*0.5;
   // minimal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3)
   double nu_1  = (a*(E-E_BIN)-P_Fermi*TMath::Sqrt(a*a+Q2*b))/b;
   nu_min= TMath::Max(nu_min,nu_1);
   return nu_min;
 }

double SmithMonizUtils::vlim2_SM ( double Q2 ) const

private

Definition at line 514 of file SmithMonizUtils.cxx.

References a, genie::units::b, E_BIN, E_nu, mm_fin, mm_ini, mm_lep, and P_Fermi.

Referenced by vQES_SM_max().

 {
   // maximal energy transfer (see Eq. (9) of Ref.3) 
   double nu_max = E_nu*Q2/(Q2+mm_lep)-(Q2+mm_lep)/4/E_nu;
   
   double E = sqrt(P_Fermi*P_Fermi+mm_ini);
   double b = (E-E_BIN)*(E-E_BIN)-P_Fermi*P_Fermi;
   double a = 0.5*(Q2+mm_fin-b);
   // maximal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3)
   double nu_2 = (a*(E-E_BIN)+P_Fermi*TMath::Sqrt(a*a+Q2*b))/b;
   nu_max= TMath::Min(nu_max,nu_2);
   return -nu_max;
 }

Range1D_t SmithMonizUtils::vQES_SM_lim ( double Q2 ) const

Definition at line 423 of file SmithMonizUtils.cxx.

References a, genie::units::b, E_BIN, E_nu, m_fin, m_rnu, m_tar, mm_fin, mm_ini, mm_lep, mm_tar, and P_Fermi.

Referenced by genie::QELEventGeneratorSM::ComputeMaxXSec(), genie::utils::gsl::dv_dQ2_E::DoEval(), PhaseSpaceVolume(), genie::QELEventGeneratorSM::ProcessEventRecord(), vQES_SM_max(), and vQES_SM_min().

 {
 
   // minimal energy transfer for scattering on nucleus (see Eq. (7) of Ref.3)
   double nu_min= ((m_rnu+m_fin)*(m_rnu+m_fin)+Q2-mm_tar)/2/m_tar;
   
   // if the target is nucleon then nu_min=nu_max=(m_fin^2+Q^2-m_ini^2)/(2*m_ini)
   if (E_BIN == 0 && P_Fermi == 0)
     return Range1D_t(nu_min, nu_min);
   
   // maximal energy transfer (see Eq. (9) of Ref.3)
   double nu_max = E_nu*Q2/(Q2+mm_lep)-(Q2+mm_lep)/4/E_nu;
   
   // now we find limits for bound nucleon
   double E = TMath::Sqrt(P_Fermi*P_Fermi+mm_ini);
   double b = (E-E_BIN)*(E-E_BIN)-P_Fermi*P_Fermi;
   double a = (Q2+mm_fin-b)*0.5;
   double tmp1 = a*(E-E_BIN);
   double tmp2  = P_Fermi*TMath::Sqrt(a*a+Q2*b);
   // minimal and maximal energy transfer for bound nucleon (see Eqs. (11) of Ref. 3)
   double nu_1 = (tmp1-tmp2)/b;
   double nu_2 = (tmp1+tmp2)/b;
   // for minimal energy transfer we take maximum of corresponding values on nucleus and bound nucleon
   nu_min= TMath::Max(nu_min,nu_1);
   // for maximal energy transfer we take minimum of corresponding values on nucleus and bound nucleon
   nu_max= TMath::Min(nu_max,nu_2);
   
   if (nu_min<=nu_max)
     return Range1D_t(nu_min,nu_max);
   else 
   // to avoid machine precision errors
     return Range1D_t(0.5*(nu_min+nu_max),0.5*(nu_min+nu_max));
 
 }

double SmithMonizUtils::vQES_SM_max	(	double	Q2min,
		double	Q2max
	)		const

Definition at line 479 of file SmithMonizUtils.cxx.

References DMINFC(), E_BIN, genie::Range1D_t::min, P_Fermi, vlim2_SM(), and vQES_SM_lim().

Referenced by genie::QELEventGeneratorSM::ComputeMaxXSec().

 {
   // maximum of function call
   const double EPS  = 1.0e-08;
   const double Delta= 1.0e-14;
      
   if (E_BIN == 0 && P_Fermi == 0)
     return vQES_SM_lim(Q2_max).min;
   
   double F_MIN, Q2_0;
   bool LLM;
   // C++ analog of fortran function DMINFC(Q2lim1_SM,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM)
   Func1D<SmithMonizUtils> vlim2_SM_(*this, &SmithMonizUtils::vlim2_SM);
   DMINFC(vlim2_SM_,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM);
   
   return -F_MIN;
   
 }

double SmithMonizUtils::vQES_SM_min	(	double	Q2min,
		double	Q2max
	)		const

Definition at line 459 of file SmithMonizUtils.cxx.

References DMINFC(), E_BIN, genie::Range1D_t::min, P_Fermi, vlim1_SM(), and vQES_SM_lim().

Referenced by genie::QELEventGeneratorSM::ComputeMaxXSec().

 {
   // maximum of function call
   const double EPS  = 1.0e-08;
   const double Delta= 1.0e-14;
    
   if (E_BIN == 0 && P_Fermi == 0)
     return vQES_SM_lim(Q2_min).min;
   
   double F_MIN, Q2_0;
   bool LLM;
   // C++ analog of fortran function DMINFC(Q2lim1_SM,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM)
   Func1D<SmithMonizUtils> vlim1_SM_(*this, &SmithMonizUtils::vlim1_SM);
   DMINFC(vlim1_SM_,Q2_min,Q2_max,EPS,Delta,Q2_0,F_MIN,LLM);
   
   return F_MIN;
   
 }