genie::SmithMonizUtils Class Reference

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.
virtual AlgStatus_t	GetStatus (void) const
	Get algorithm status.
virtual bool	AllowReconfig (void) const
virtual AlgCmp_t	Compare (const Algorithm *alg) const
	Compare with input algorithm.
virtual void	SetId (const AlgId &id)
	Set algorithm ID.
virtual void	SetId (string name, string config)
const Algorithm *	SubAlg (const RgKey &registry_key) const
void	AdoptConfig (void)
void	AdoptSubstructure (void)
virtual void	Print (ostream &stream) const
	Print algorithm info.

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)
double	m_lep
	Mass of final charged lepton (GeV)
double	mm_lep
	Squared mass of final charged lepton (GeV)
double	m_ini
	Mass of initial hadron or hadron system (GeV)
double	mm_ini
	Sqared mass of initial hadron or hadron system (GeV)
double	m_fin
	Mass of final hadron or hadron system (GeV)
double	mm_fin
	Squared mass of final hadron or hadron system (GeV)
double	m_tar
	Mass of target nucleus (GeV)
double	mm_tar
	Squared mass of target nucleus (GeV)
double	m_rnu
	Mass of residual nucleus (GeV)
double	mm_rnu
	Squared mass of residual nucleus (GeV)
double	P_Fermi
	Maximum value of Fermi momentum of target nucleon (GeV)
double	E_BIN
	Binding energy (GeV)

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.
template<class T>
bool	GetParam (const RgKey &name, T &p, bool is_top_call=true) const
template<class T>
bool	GetParamDef (const RgKey &name, T &p, const T &def) const
template<class T>
int	GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
	Handle to load vectors of parameters.
int	GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
template<class T>
int	GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
	Handle to load matrix of parameters.
template<class T>
int	GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const
int	GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
int	AddTopRegistry (Registry *rp, bool owns=true)
	add registry with top priority, also update ownership
int	AddLowRegistry (Registry *rp, bool owns=true)
	add registry with lowest priority, also update ownership
int	MergeTopRegistry (const Registry &r)
int	AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
	Add registries with top priority, also udated Ownerships.
Protected Attributes inherited from genie::Algorithm
bool	fAllowReconfig
bool	fOwnsSubstruc
	true if it owns its substructure (sub-algs,...)
AlgId	fID
	algorithm name and configuration set
vector< Registry * >	fConfVect
vector< bool >	fOwnerships
	ownership for every registry in fConfVect
AlgStatus_t	fStatus
	algorithm execution status
AlgMap *	fOwnedSubAlgMp
	local pool for owned sub-algs (taken out of the factory pool)

Detailed Description

Contains auxiliary functions for Smith-Moniz model.
.

References:\n [1] R.A.Smith and E.J.Moniz, Nuclear Physics B43, (1972) 605-622 \n

[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:\n May 05, 2017

License:\n Copyright (c) 2003-2025, The GENIE Collaboration

For the full text of the license visit http://copyright.genie-mc.org

Definition at line 53 of file SmithMonizUtils.h.

Constructor & Destructor Documentation

SmithMonizUtils() [1/2]

SmithMonizUtils::SmithMonizUtils

(

)

Definition at line 46 of file SmithMonizUtils.cxx.

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

References genie::Algorithm::Algorithm().

SmithMonizUtils() [2/2]

SmithMonizUtils::SmithMonizUtils

(

string

config

)

Definition at line 52 of file SmithMonizUtils.cxx.

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

References genie::Algorithm::Algorithm().

~SmithMonizUtils()

SmithMonizUtils::~SmithMonizUtils ( )

virtual

Definition at line 58 of file SmithMonizUtils.cxx.

{
 
}

Member Function Documentation

Configure() [1/2]

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.

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

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

Configure() [2/2]

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.

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

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

DMINFC()

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.

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

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

E_nu_thr_SM()

double SmithMonizUtils::E_nu_thr_SM

(

void

)

const

Definition at line 210 of file SmithMonizUtils.cxx.

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

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

GetBindingEnergy()

double SmithMonizUtils::GetBindingEnergy

(

void

)

const

Definition at line 629 of file SmithMonizUtils.cxx.

{
  return E_BIN;
}

References E_BIN.

GetFermiMomentum()

double SmithMonizUtils::GetFermiMomentum

(

void

)

const

Definition at line 634 of file SmithMonizUtils.cxx.

{
  return P_Fermi;
}

References P_Fermi.

GetTheta_k()

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

Definition at line 639 of file SmithMonizUtils.cxx.

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

References E_nu, and mm_lep.

GetTheta_p()

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

Definition at line 644 of file SmithMonizUtils.cxx.

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

References E_BIN, mm_fin, and mm_ini.

kFQES_SM_lim()

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

Definition at line 529 of file SmithMonizUtils.cxx.

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

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

Referenced by PhaseSpaceVolume().

LoadConfig()

void SmithMonizUtils::LoadConfig ( void )

private

Definition at line 75 of file SmithMonizUtils.cxx.

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

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

Referenced by Configure(), and Configure().

PhaseSpaceVolume()

double SmithMonizUtils::PhaseSpaceVolume

(

KinePhaseSpace_t

ps

)

const

Definition at line 653 of file SmithMonizUtils.cxx.

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

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

Q2lim1_SM()

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

private

Definition at line 273 of file SmithMonizUtils.cxx.

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

References a, E_BIN, mm_fin, mm_ini, mm_lep, and P_Fermi.

Referenced by Q2QES_SM_lim(), and QEL_EnuMin_SM().

Q2lim2_SM()

double SmithMonizUtils::Q2lim2_SM ( double Q2 ) const

private

Definition at line 288 of file SmithMonizUtils.cxx.

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

References a, E_BIN, m_fin, m_rnu, m_tar, mm_fin, mm_ini, mm_tar, and P_Fermi.

Referenced by Q2QES_SM_lim().

Q2QES_SM_lim()

Range1D_t SmithMonizUtils::Q2QES_SM_lim

(

void

)

const

Definition at line 303 of file SmithMonizUtils.cxx.

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

References DMINFC(), E_BIN, E_nu, 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(), and Q2lim2_SM().

Referenced by PhaseSpaceVolume().

QEL_EnuMin_SM()

double SmithMonizUtils::QEL_EnuMin_SM

(

double

E_nu

)

const

Definition at line 247 of file SmithMonizUtils.cxx.

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

References DMINFC(), m_fin, m_lep, m_rnu, m_tar, mm_lep, mm_tar, and Q2lim1_SM().

Referenced by E_nu_thr_SM().

rho()

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

static

Definition at line 609 of file SmithMonizUtils.cxx.

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

References P_Fermi.

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

SetInteraction()

void SmithMonizUtils::SetInteraction

(

const Interaction *

i

)

Definition at line 116 of file SmithMonizUtils.cxx.

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

References genie::Target::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::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().

vlim1_SM()

double SmithMonizUtils::vlim1_SM ( double Q2 ) const

private

Definition at line 499 of file SmithMonizUtils.cxx.

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

References a, E_BIN, m_fin, m_rnu, m_tar, mm_fin, mm_ini, mm_tar, and P_Fermi.

Referenced by vQES_SM_min().

vlim2_SM()

double SmithMonizUtils::vlim2_SM ( double Q2 ) const

private

Definition at line 514 of file SmithMonizUtils.cxx.

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

References a, E_BIN, E_nu, mm_fin, mm_ini, mm_lep, and P_Fermi.

Referenced by vQES_SM_max().

vQES_SM_lim()

Range1D_t SmithMonizUtils::vQES_SM_lim

(

double

Q2

)

const

Definition at line 423 of file SmithMonizUtils.cxx.

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

References a, E_BIN, E_nu, m_fin, m_rnu, m_tar, mm_fin, mm_ini, mm_lep, mm_tar, and P_Fermi.

Referenced by PhaseSpaceVolume(), vQES_SM_max(), and vQES_SM_min().

vQES_SM_max()

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

Definition at line 479 of file SmithMonizUtils.cxx.

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

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

vQES_SM_min()

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

Definition at line 459 of file SmithMonizUtils.cxx.

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

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

Member Data Documentation

E_BIN

double genie::SmithMonizUtils::E_BIN

private

Binding energy (GeV)

Definition at line 142 of file SmithMonizUtils.h.

Referenced by E_nu_thr_SM(), GetBindingEnergy(), GetTheta_p(), kFQES_SM_lim(), Q2lim1_SM(), Q2lim2_SM(), Q2QES_SM_lim(), SetInteraction(), vlim1_SM(), vlim2_SM(), vQES_SM_lim(), vQES_SM_max(), and vQES_SM_min().

E_nu

double genie::SmithMonizUtils::E_nu

private

Neutrino energy (GeV)

Definition at line 130 of file SmithMonizUtils.h.

Referenced by GetTheta_k(), Q2QES_SM_lim(), SetInteraction(), vlim2_SM(), and vQES_SM_lim().

fInteraction

const Interaction* genie::SmithMonizUtils::fInteraction

private

Definition at line 125 of file SmithMonizUtils.h.

Referenced by SetInteraction().

fKFTable

string genie::SmithMonizUtils::fKFTable

private

Definition at line 122 of file SmithMonizUtils.h.

Referenced by LoadConfig(), and SetInteraction().

fNucRmvE

map<int, double> genie::SmithMonizUtils::fNucRmvE

private

Definition at line 121 of file SmithMonizUtils.h.

Referenced by LoadConfig(), and SetInteraction().

fUseParametrization

bool genie::SmithMonizUtils::fUseParametrization

private

Definition at line 123 of file SmithMonizUtils.h.

Referenced by LoadConfig(), and SetInteraction().

m_fin

double genie::SmithMonizUtils::m_fin

private

Mass of final hadron or hadron system (GeV)

Definition at line 135 of file SmithMonizUtils.h.

Referenced by E_nu_thr_SM(), Q2lim2_SM(), Q2QES_SM_lim(), QEL_EnuMin_SM(), SetInteraction(), vlim1_SM(), and vQES_SM_lim().

m_ini

double genie::SmithMonizUtils::m_ini

private

Mass of initial hadron or hadron system (GeV)

Definition at line 133 of file SmithMonizUtils.h.

Referenced by kFQES_SM_lim(), Q2QES_SM_lim(), and SetInteraction().

m_lep

double genie::SmithMonizUtils::m_lep

private

Mass of final charged lepton (GeV)

Definition at line 131 of file SmithMonizUtils.h.

Referenced by E_nu_thr_SM(), Q2QES_SM_lim(), QEL_EnuMin_SM(), and SetInteraction().

m_rnu

double genie::SmithMonizUtils::m_rnu

private

Mass of residual nucleus (GeV)

Definition at line 139 of file SmithMonizUtils.h.

Referenced by E_nu_thr_SM(), Q2lim2_SM(), Q2QES_SM_lim(), QEL_EnuMin_SM(), SetInteraction(), vlim1_SM(), and vQES_SM_lim().

m_tar

double genie::SmithMonizUtils::m_tar

private

Mass of target nucleus (GeV)

Definition at line 137 of file SmithMonizUtils.h.

Referenced by E_nu_thr_SM(), Q2lim2_SM(), Q2QES_SM_lim(), QEL_EnuMin_SM(), SetInteraction(), vlim1_SM(), and vQES_SM_lim().

mm_fin

double genie::SmithMonizUtils::mm_fin

private

Squared mass of final hadron or hadron system (GeV)

Definition at line 136 of file SmithMonizUtils.h.

Referenced by GetTheta_p(), kFQES_SM_lim(), Q2lim1_SM(), Q2lim2_SM(), Q2QES_SM_lim(), SetInteraction(), vlim1_SM(), vlim2_SM(), and vQES_SM_lim().

mm_ini

double genie::SmithMonizUtils::mm_ini

private

Sqared mass of initial hadron or hadron system (GeV)

Definition at line 134 of file SmithMonizUtils.h.

Referenced by E_nu_thr_SM(), GetTheta_p(), kFQES_SM_lim(), Q2lim1_SM(), Q2lim2_SM(), Q2QES_SM_lim(), SetInteraction(), vlim1_SM(), vlim2_SM(), and vQES_SM_lim().

mm_lep

double genie::SmithMonizUtils::mm_lep

private

Squared mass of final charged lepton (GeV)

Definition at line 132 of file SmithMonizUtils.h.

Referenced by GetTheta_k(), Q2lim1_SM(), Q2QES_SM_lim(), QEL_EnuMin_SM(), SetInteraction(), vlim2_SM(), and vQES_SM_lim().

mm_rnu

double genie::SmithMonizUtils::mm_rnu

private

Squared mass of residual nucleus (GeV)

Definition at line 140 of file SmithMonizUtils.h.

Referenced by SetInteraction().

mm_tar

double genie::SmithMonizUtils::mm_tar

private

Squared mass of target nucleus (GeV)

Definition at line 138 of file SmithMonizUtils.h.

Referenced by E_nu_thr_SM(), Q2lim2_SM(), Q2QES_SM_lim(), QEL_EnuMin_SM(), SetInteraction(), vlim1_SM(), and vQES_SM_lim().

P_Fermi

double genie::SmithMonizUtils::P_Fermi

private

Maximum value of Fermi momentum of target nucleon (GeV)

Definition at line 141 of file SmithMonizUtils.h.

Referenced by E_nu_thr_SM(), GetFermiMomentum(), kFQES_SM_lim(), Q2lim1_SM(), Q2lim2_SM(), Q2QES_SM_lim(), rho(), SetInteraction(), vlim1_SM(), vlim2_SM(), vQES_SM_lim(), vQES_SM_max(), and vQES_SM_min().

---

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

• SmithMonizUtils.h
• SmithMonizUtils.cxx