genie::utils::nuclear Namespace Reference

Simple nuclear physics empirical formulas (densities, radii, ...) and empirical nuclear corrections. More...

Functions
double	BindEnergy (const Target &target)
double	BindEnergy (int nucA, int nucZ)
double	BindEnergyPerNucleon (const Target &target)
double	BindEnergyLastNucleon (const Target &target)
double	Radius (int A, double Ro=constants::kNucRo)
double	NuclQELXSecSuppression (string kftable, double pmax, const Interaction *in)
double	RQEFG_generic (double q2, double Mn, double kFi, double kFf, double pmax)
double	FmI1 (double alpha, double beta, double a, double b, double kFi, double kFf, double q)
double	FmI2 (double alpha, double beta, double a, double b, double kFi, double kFf, double q)
double	FmArea (double alpha, double beta, double kf, double pmax)
double	DISNuclFactor (double x, int A)
double	Density (double r, int A, double ring=0.)
double	DensityGaus (double r, double ap, double alf, double ring=0.)
double	DensityWoodsSaxon (double r, double c, double z, double ring=0.)
double	BindEnergyPerNucleonParametrization (const Target &target)
double	FermiMomentumForIsoscalarNucleonParametrization (const Target &target)

Detailed Description

Simple nuclear physics empirical formulas (densities, radii, ...) and empirical nuclear corrections.

Author

Costas Andreopoulos <c.andreopoulos \at cern.ch> University of Liverpool

Created:\n May 06, 2004

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

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

Function Documentation

BindEnergy() [1/2]

double genie::utils::nuclear::BindEnergy

(

const Target &

target

)

Definition at line 59 of file NuclearUtils.cxx.

{
// Compute the average binding energy (in GeV) using the semi-empirical
// formula from Wapstra (Handbuch der Physik, XXXVIII/1)
 
  if(!target.IsNucleus()) return 0;
 
  int Z = (int) target.Z();
  int A = (int) target.A();
 
  return BindEnergy(A,Z);
}

References genie::Target::A(), BindEnergy(), genie::Target::IsNucleus(), and genie::Target::Z().

Referenced by BindEnergy(), BindEnergyLastNucleon(), and BindEnergyPerNucleon().

BindEnergy() [2/2]

double genie::utils::nuclear::BindEnergy	(	int	nucA,
		int	nucZ )

Definition at line 72 of file NuclearUtils.cxx.

{
// Compute the average binding energy (in GeV) using the semi-empirical
// formula from Wapstra (Handbuch der Physik, XXXVIII/1)
 
  if (nucA<=0 || nucZ<=0) return 0;
  
  if (nucZ == 2)
  {
   if (nucA == 3)
     return 7.5e-3;
   if (nucA == 4)
     return 28.4e-3;
   return 0.;
  }
 
  double a =  15.835;
  double b =  18.33;
  double s =  23.20;
  double d =   0.714;
 
  double delta = 0;                          /*E-O*/
  if (nucZ%2==1 && nucA%2==1) delta =  11.2; /*O-O*/
  if (nucZ%2==0 && nucA%2==0) delta = -11.2; /*E-E*/
 
  double N = (double) (nucA-nucZ);
  double Z = (double) nucZ;
  double A = (double) nucA;
 
  double BE =  a * A
             - b * TMath::Power(A,0.667)
             - s * TMath::Power(N-Z,2.0)/A
             - d * TMath::Power(Z,2.0)/TMath::Power(A,0.333)
             - delta / TMath::Sqrt(A); // MeV
 
  return (1e-3 * BE); // GeV
}

References a, and e.

BindEnergyLastNucleon()

double genie::utils::nuclear::BindEnergyLastNucleon

(

const Target &

target

)

Definition at line 119 of file NuclearUtils.cxx.

{
// Compute the binding for the most loose nucleon (in GeV)
 
  if(!target.IsNucleus()) return 0;
 
  //-- temporarily, return the binding energy per nucleon rather than the
  //   separation energy of the last nucleon
   return (utils::nuclear::BindEnergy(target) / target.A());
}

References genie::Target::A(), BindEnergy(), and genie::Target::IsNucleus().

BindEnergyPerNucleon()

double genie::utils::nuclear::BindEnergyPerNucleon

(

const Target &

target

)

Definition at line 110 of file NuclearUtils.cxx.

{
// Compute the average binding energy per nucleon (in GeV)
 
  if(!target.IsNucleus()) return 0;
 
  return (utils::nuclear::BindEnergy(target) / target.A());
}

References genie::Target::A(), BindEnergy(), and genie::Target::IsNucleus().

Referenced by genie::FGMBodekRitchie::GenerateNucleon(), genie::LocalFGM::GenerateNucleon(), genie::SpectralFunc1d::GenerateNucleon(), and genie::SmithMonizUtils::SetInteraction().

BindEnergyPerNucleonParametrization()

double genie::utils::nuclear::BindEnergyPerNucleonParametrization

(

const Target &

target

)

Definition at line 493 of file NuclearUtils.cxx.

{
// Compute the average binding energy per nucleon (in GeV)
if(!target.IsNucleus()) return 0;
double x = TMath::Power(target.A(),1/3.0) / target.Z();
return (0.05042772591+x*(-0.11377355795+x*(0.15159890400-0.08825307197*x)));
}

References genie::Target::A(), genie::Target::IsNucleus(), and genie::Target::Z().

Referenced by genie::FGMBodekRitchie::GenerateNucleon(), and genie::SmithMonizUtils::SetInteraction().

Density()

double genie::utils::nuclear::Density	(	double	r,
		int	A,
		double	ring = 0. )

Definition at line 396 of file NuclearUtils.cxx.

{
// [by S.Dytman]
//
  if(A>20) {
    double c = 1., z = 1.;
 
    if      (A ==  27) { c = 3.07; z = 0.52; }  // aluminum
    else if (A ==  28) { c = 3.07; z = 0.54; }  // silicon
    else if (A ==  40) { c = 3.53; z = 0.54; }  // argon
    else if (A ==  56) { c = 4.10; z = 0.56; }  // iron
    else if (A == 208) { c = 6.62; z = 0.55; }  // lead
    else { 
       c = TMath::Power(A,0.35); z = 0.54; 
    } //others
 
    LOG("Nuclear",pINFO)
        << "r= " << r << ", ring= " << ring;
    double rho = DensityWoodsSaxon(r,c,z,ring);
    return rho;
  }
  else if (A>4) {
    double ap = 1., alf = 1.;
 
    if      (A ==  7) { ap = 1.77; alf = 0.327; } // lithium
    else if (A == 12) { ap = 1.69; alf = 1.08;  } // carbon
    else if (A == 14) { ap = 1.76; alf = 1.23;  } // nitrogen
    else if (A == 16) { ap = 1.83; alf = 1.54;  } // oxygen
    else  { 
      ap=1.75; alf=-0.4+.12*A; 
    }  //others- alf=0.08 if A=4
 
    double rho = DensityGaus(r,ap,alf,ring);
    return rho;
  }
  else {
    // helium
    double ap = 1.9/TMath::Sqrt(2.);  
    double alf=0.;    
    double rho = DensityGaus(r,ap,alf,ring);
    return rho;
  }
 
  return 0;
}

References DensityGaus(), DensityWoodsSaxon(), LOG, and pINFO.

Referenced by CheckVertexDistribution(), genie::NievesQELCCPXSec::CNCTCLimUcalc(), genie::utils::gsl::wrap::NievesQELvcrIntegrand::DoEval(), genie::utils::gsl::wrap::NuclDensityMomentIntegrand::DoEval(), genie::NucleonDecayPrimaryVtxGenerator::GenerateDecayedNucleonPosition(), genie::NNBarOscPrimaryVtxGenerator::GenerateOscillatingNeutronPosition(), genie::VertexGenerator::GenerateVertex(), genie::PauliBlocker::GetFermiMomentum(), genie::LocalFGM::LocalFermiMomentum(), genie::utils::intranuke2018::MeanFreePath(), genie::utils::intranuke2025::MeanFreePath(), genie::utils::intranuke::MeanFreePath(), genie::utils::intranuke2018::MeanFreePath_Delta(), genie::utils::intranuke2025::MeanFreePath_Delta(), genie::utils::intranuke::MeanFreePath_Delta(), and NuclQELXSecSuppression().

DensityGaus()

double genie::utils::nuclear::DensityGaus	(	double	r,
		double	ap,
		double	alf,
		double	ring = 0. )

Definition at line 442 of file NuclearUtils.cxx.

{
// [adapted from neugen3 density_gaus.F written by S.Dytman]
//
// Modified harmonic osc density distribution. 
// Norm gives normalization to 1
//
// input  : radial distance in nucleus [units: fm]
// output : nuclear density            [units: fm^-3]
 
  ring = TMath::Min(ring, 0.3*a);
 
  double aeval = a + ring;
  double norm  = 1./((5.568 + alf*8.353)*TMath::Power(a,3.));  //0.0132;
  double b     = TMath::Power(r/aeval, 2.);
  double dens  = norm * (1. + alf*b) * TMath::Exp(-b);
 
  LOG("Nuclear", pINFO) 
        << "r = " << r << ", norm = " << norm << ", dens = " << dens 
        << ", aeval= " << aeval;
 
  return dens;
}

References a, LOG, and pINFO.

Referenced by Density().

DensityWoodsSaxon()

double genie::utils::nuclear::DensityWoodsSaxon	(	double	r,
		double	c,
		double	z,
		double	ring = 0. )

Definition at line 467 of file NuclearUtils.cxx.

{
// [adapted from neugen3 density_ws.F written by S.Dytman]
//
// Woods-Saxon desity distribution. Norn gives normalization to 1
//
// input  : radial distance in nucleus [units: fm]
// output : nuclear density            [units: fm^-3]
 
  LOG("Nuclear",pINFO)
      << "c= " << c << ", z= " << z << ",ring= " << ring;
 
  ring = TMath::Min(ring, 0.75*c);
 
  double ceval = c + ring;
  double norm  = (3./(4.*kPi*TMath::Power(c,3)))*1./(1.+TMath::Power((kPi*z/c),2));
  double dens  = norm / (1 + TMath::Exp((r-ceval)/z));
 
  LOG("Nuclear", pINFO) 
     << "r = " << r << ", norm = " << norm 
     << ", dens = " << dens << " , ceval= " << ceval;
 
  return dens;
}

References genie::constants::kPi, LOG, and pINFO.

Referenced by Density().

DISNuclFactor()

double genie::utils::nuclear::DISNuclFactor	(	double	x,
		int	A )

Definition at line 371 of file NuclearUtils.cxx.

{
// Adapted from NeuGEN's nuc_factor(). Kept original comments from Hugh.
 
  double xv     = TMath::Min(0.75, x);
  double xv2    = xv  * xv;
  double xv3    = xv2 * xv;
  double xv4    = xv3 * xv;
  double xv5    = xv4 * xv;
  double xvp    = TMath::Power(xv, 14.417);
  double expaxv = TMath::Exp(-21.94*xv);
 
  double f = 1.;
 
  // first factor goes from free nucleons to deuterium
  if(A >= 2) {
    f= 0.985*(1.+0.422*xv - 2.745*xv2 + 7.570*xv3 - 10.335*xv4 + 5.422*xv5);
  }
  // 2nd factor goes from deuterium to iso-scalar iron
  if(A > 2) {
    f *= (1.096 - 0.364*xv - 0.278*expaxv + 2.722*xvp);
  }
  return f;
}

Referenced by genie::QPMDISStrucFuncBase::NuclMod(), and genie::QPMDMDISStrucFuncBase::NuclMod().

FermiMomentumForIsoscalarNucleonParametrization()

double genie::utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization

(

const Target &

target

)

Definition at line 501 of file NuclearUtils.cxx.

{
// Compute Fermi momentum for isoscalar nucleon (in GeV)
if(!target.IsNucleus()) return 0;
double x = 1.0 / target.A();
return (0.27+x*(-1.12887857491+x*(9.72670908033-39.53095724456*x)));
}

References genie::Target::A(), and genie::Target::IsNucleus().

Referenced by genie::FGMBodekRitchie::FermiMomentum(), genie::SmithMonizUtils::SetInteraction(), genie::BSKLNBaseRESPXSec2014::XSec(), genie::MKSPPPXSec2020::XSec(), and genie::ReinSehgalRESPXSec::XSec().

FmArea()

double genie::utils::nuclear::FmArea	(	double	alpha,
		double	beta,
		double	kf,
		double	pmax )

Definition at line 361 of file NuclearUtils.cxx.

{
// Adapted from NeuGEN's fm_area() used in r_factor()
 
  double kf3 = TMath::Power(kf,3.);
  double sum = 4.*kPi* (alpha * kf3/3. + beta*(1./kf - 1./pmax));
  return sum;
}

References genie::constants::kPi.

Referenced by RQEFG_generic().

FmI1()

double genie::utils::nuclear::FmI1	(	double	alpha,
		double	beta,
		double	a,
		double	b,
		double	kFi,
		double	kFf,
		double	q )

Definition at line 298 of file NuclearUtils.cxx.

{
// Adapted from NeuGEN's fm_integral1() used in r_factor()
 
  double f=0;
 
  double q2   = TMath::Power(q,  2);
  double a2   = TMath::Power(a,  2);
  double b2   = TMath::Power(b,  2);
  double kFi2 = TMath::Power(kFi,2);
  double kFf2 = TMath::Power(kFf,2);
 
  if(kFi < a) {
     double lg = TMath::Log(b/a);
 
     f = -beta * (kFf2-q2)/(4.*q) * (1./a2 - 1./b2) + beta/(2.*q) * lg;
 
  } else if (kFi > b) {
     double a4 = TMath::Power(a2,2);
     double b4 = TMath::Power(b2,2);
 
     f = - (kFf2-q2) * alpha/(4.*q) * (b2-a2) + alpha/(8.*q) * (b4-a4);
 
  } else {
     double a4   = TMath::Power(a2,2);
     double kFi4 = TMath::Power(kFi2,2);
     double lg   = TMath::Log(b/kFi);
 
     f = - alpha*(kFf2-q2)/(4.*q)*(kFi2-a2) + alpha/(8.*q)*(kFi4 - a4)
         - beta*(kFf2-q2)/(4.*q)*(1./kFi2 - 1./b2) + beta/(2.*q)*lg;
  }
 
  double integral2 = FmI2(alpha,beta,a,b,kFi,kFf,q);
  double integral1 = integral2 + f;
  
  return integral1;
}

References a, and FmI2().

Referenced by RQEFG_generic().

FmI2()

double genie::utils::nuclear::FmI2	(	double	alpha,
		double	beta,
		double	a,
		double	b,
		double	kFi,
		double	kFf,
		double	q )

Definition at line 337 of file NuclearUtils.cxx.

{
// Adapted from NeuGEN's fm_integral2() used in r_factor()
 
  double integral2 = 0;
 
  if(kFi < a) {
     integral2 = beta * (1./a - 1./b);
 
  } else if(kFi > b) {
     double a3 = TMath::Power(a,3);
     double b3 = TMath::Power(b,3);
     integral2 = alpha/3. * (b3 - a3);
 
  } else {
     double a3   = TMath::Power(a,3);
     double kFi3 = TMath::Power(kFi,3);
 
     integral2 = alpha/3. * (kFi3 - a3) + beta * (1./kFi - 1./b);
  }
  return integral2;
}

References a.

Referenced by FmI1(), and RQEFG_generic().

NuclQELXSecSuppression()

double genie::utils::nuclear::NuclQELXSecSuppression	(	string	kftable,
		double	pmax,
		const Interaction *	in )

Definition at line 140 of file NuclearUtils.cxx.

{
  const InitialState & init_state = interaction->InitState();
  const ProcessInfo &  proc_info  = interaction->ProcInfo();
  const Target &       target     = init_state.Tgt();
 
  if (!target.IsNucleus()) {
     return 1.0; // no suppression for free nucleon tgt
  }
 
  //
  // special case: deuterium   
  //
 
  if (target.A() == 2) {
    NuclearData * nucldat = NuclearData::Instance();
    return nucldat->DeuteriumSuppressionFactor(interaction->Kine().Q2());
  }
 
  //
  // general case
  //
 
  int target_pdgc         = target.Pdg();
  int struck_nucleon_pdgc = target.HitNucPdg();
  int final_nucleon_pdgc  = struck_nucleon_pdgc;
 
  if(proc_info.IsWeakCC()) {
     final_nucleon_pdgc = pdg::SwitchProtonNeutron(struck_nucleon_pdgc);
  }
 
  // Check if an LFG model should be used for Fermi momentum
  // Create a nuclear model object to check the model type
  AlgConfigPool * confp = AlgConfigPool::Instance();
  const Registry * gc = confp->GlobalParameterList();
  RgKey nuclkey = "NuclearModel";
  RgAlg nuclalg = gc->GetAlg(nuclkey);
  AlgFactory * algf = AlgFactory::Instance();
  const genie::NuclearModelI* nuclModel =
    dynamic_cast<const genie::NuclearModelI*>(
                             algf->GetAlgorithm(nuclalg.name,nuclalg.config));
  // Check if the model is a local Fermi gas
  bool lfg = (nuclModel && nuclModel->ModelType(Target()) == kNucmLocalFermiGas);
 
  double kFi, kFf;
  if(lfg){
    double hbarc = kLightSpeed*kPlankConstant/genie::units::fermi;
    Target* tgt = interaction->InitStatePtr()->TgtPtr();
    double radius = tgt->HitNucPosition();
 
    int A = tgt->A();
    // kFi
    bool is_p_i = pdg::IsProton(struck_nucleon_pdgc);
    double numNuci = (is_p_i) ? (double)tgt->Z():(double)tgt->N();
    kFi = TMath::Power(3*kPi2*numNuci*
                     genie::utils::nuclear::Density(radius,A),1.0/3.0) *hbarc;
    // kFi
    bool is_p_f = pdg::IsProton(final_nucleon_pdgc);
    double numNucf = (is_p_f) ? (double)tgt->Z():(double)tgt->N();
    kFf = TMath::Power(3*kPi2*numNucf*
                     genie::utils::nuclear::Density(radius,A),1.0/3.0) *hbarc;
  }else{
    // get the requested Fermi momentum table 
    FermiMomentumTablePool * kftp = FermiMomentumTablePool::Instance();
    const FermiMomentumTable * kft  = kftp->GetTable(kftable);
    
    // Fermi momentum for initial, final nucleons
    kFi = kft->FindClosestKF(target_pdgc, struck_nucleon_pdgc);
    kFf = (struck_nucleon_pdgc==final_nucleon_pdgc) ? kFi : 
          kft->FindClosestKF(target_pdgc, final_nucleon_pdgc );
  }
  
  double Mn = target.HitNucP4Ptr()->M(); // can be off m/shell
 
  const Kinematics & kine = interaction->Kine();
  double q2 = kine.q2();
 
  double R = RQEFG_generic(q2, Mn, kFi, kFf, pmax);
  return R;
}

References genie::Target::A(), RgAlg::config, Density(), genie::NuclearData::DeuteriumSuppressionFactor(), genie::units::fermi, genie::FermiMomentumTable::FindClosestKF(), genie::Registry::GetAlg(), genie::AlgFactory::GetAlgorithm(), genie::FermiMomentumTablePool::GetTable(), genie::AlgConfigPool::GlobalParameterList(), genie::Target::HitNucP4Ptr(), genie::Target::HitNucPdg(), genie::Target::HitNucPosition(), genie::Interaction::InitState(), genie::Interaction::InitStatePtr(), genie::AlgConfigPool::Instance(), genie::AlgFactory::Instance(), genie::FermiMomentumTablePool::Instance(), genie::NuclearData::Instance(), genie::Target::IsNucleus(), genie::pdg::IsProton(), genie::ProcessInfo::IsWeakCC(), genie::Interaction::Kine(), genie::constants::kLightSpeed, genie::kNucmLocalFermiGas, genie::constants::kPi2, genie::constants::kPlankConstant, genie::NuclearModelI::ModelType(), genie::Target::N(), RgAlg::name, genie::Target::Pdg(), genie::Interaction::ProcInfo(), genie::Kinematics::Q2(), genie::Kinematics::q2(), RQEFG_generic(), genie::pdg::SwitchProtonNeutron(), genie::InitialState::Tgt(), genie::InitialState::TgtPtr(), and genie::Target::Z().

Referenced by genie::AhrensDMELPXSec::XSec(), genie::AhrensNCELPXSec::XSec(), genie::LwlynSmithQELCCPXSec::XSec(), and genie::RosenbluthPXSec::XSec().

Radius()

double genie::utils::nuclear::Radius	(	int	A,
		double	Ro = constants::kNucRo )

Definition at line 130 of file NuclearUtils.cxx.

{
// Compute the nuclear radius (in fm)
 
  if(A<1) return 0;
 
  double Rn = Ro*TMath::Power(1.*A, 0.3333333);
  return Rn;
}

Referenced by genie::masterclass::MCTruthDisplay::DrawDiagram(), genie::PattonCEvNSPXSec::NuclearDensityMoment(), and genie::INukeDeltaPropg::ProcessEventRecord().

RQEFG_generic()

double genie::utils::nuclear::RQEFG_generic	(	double	q2,
		double	Mn,
		double	kFi,
		double	kFf,
		double	pmax )

Definition at line 222 of file NuclearUtils.cxx.

{
// Computes the nuclear suppression of the QEL differential cross section 
//
// Inputs:
//  - q2   : momentum transfer (< 0)
//  - Mn   : hit nucleon mass (nucleon can be off the mass shell)
//  - kFi  : Fermi momentum, initial state (hit) nucleon @ nuclear target
//  - kFf  : Fermi momentum, final state (recoil) nucleon @ nuclear target
//  - pmax : A Fermi momentum cutoff
//
// A direct adaptation of NeuGEN's qelnuc() and r_factor()
//
// Hugh's comments from the original code: 
// "This routine is based on an analytic calculation of the rejection factor 
//  in the Fermi Gas model using the form for the fermi momentum distribution 
//  given in the  Bodek and Ritchie paper.  [P.R.  D23 (1981) 1070]
//  R is the ratio of the differential cross section from the nuclear material 
//  specified by (kf,pf) to the differential cross section for a free nucleon".
//  (kf,pf = Fermi Gas model Fermi momentum for initial,final nucleons)
//
 
  double R = 1.;
 
  // Compute magnitude of the 3-momentum transfer to the nucleon
  double Mn2    = Mn*Mn;
  double magq2  = q2 * (0.25*q2/Mn2 - 1.);
  double q      = TMath::Sqrt(TMath::Max(0.,magq2));
 
  double kfa   = kFi * 2./kPi;
  double kfa2  = TMath::Power(kfa,2);
  double kFi4  = TMath::Power(kFi,4);
  double rkf   = 1./(1. - kFi/4.);
  double alpha = 1. - 6.*kfa2;
  double beta  = 2. * rkf * kfa2 * kFi4;
 
  double fm_area = FmArea(alpha,beta,kFi,pmax);
 
  if (q <= kFf) {
 
     double p1   = kFf - q;
     double p2   = kFf + q;
     double fmi1 = FmI1(alpha,beta,p1,p2,  kFi,kFf,q);
     double fmi2 = FmI2(alpha,beta,p2,pmax,kFi,kFf,q);
 
     R = 2*kPi * (fmi1 + 2*fmi2) / fm_area;
 
  } else if (q > kFf && q <= (pmax-kFf)) {
 
     double p1    = q - kFf;
     double p2    = q + kFf;
     double fmi1  = FmI1(alpha,beta, p1, p2,   kFi, kFf,q);
     double fmi2a = FmI2(alpha,beta, 0., p1,   kFi, kFf,q);
     double fmi2b = FmI2(alpha,beta, p2, pmax, kFi, kFf,q);
 
     R = 2*kPi * (fmi1 + 2*(fmi2a+fmi2b)) / fm_area;
 
  } else if (q > (pmax-kFf) && q <= (pmax+kFf)) {
 
     double p1   = q - kFf;
     double fmi2 = FmI2(alpha,beta, 0.,p1,  kFi,kFf,q);
     double fmi1 = FmI1(alpha,beta, p1,pmax,kFi,kFf,q);
 
     R = 2*kPi * (2.*fmi2 + fmi1) / fm_area;
 
  } else if (q > (pmax+kFf)) {
     R = 1.; 
 
  } else {
     LOG("Nuclear", pFATAL) << "Illegal input q = " << q;
     exit(1);
  }
  return R;
}

References FmArea(), FmI1(), FmI2(), genie::constants::kPi, LOG, and pFATAL.

Referenced by main(), and NuclQELXSecSuppression().