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 cern.ch> University of Liverpool

Created:

May 06, 2004

License:

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

Function Documentation

double genie::utils::nuclear::BindEnergy

(

const Target &

target

)

Definition at line 59 of file NuclearUtils.cxx.

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

Referenced by BindEnergyLastNucleon(), and BindEnergyPerNucleon().

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

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

Definition at line 72 of file NuclearUtils.cxx.

References genie::units::A, a, genie::units::b, e, and genie::units::s.

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

double genie::utils::nuclear::BindEnergyLastNucleon

(

const Target &

target

)

Definition at line 119 of file NuclearUtils.cxx.

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

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

double genie::utils::nuclear::BindEnergyPerNucleon

(

const Target &

target

)

Definition at line 110 of file NuclearUtils.cxx.

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

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

{
// Compute the average binding energy per nucleon (in GeV)

  if(!target.IsNucleus()) return 0;

  return (utils::nuclear::BindEnergy(target) / target.A());
}

double genie::utils::nuclear::BindEnergyPerNucleonParametrization

(

const Target &

target

)

Definition at line 493 of file NuclearUtils.cxx.

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

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

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

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

Definition at line 396 of file NuclearUtils.cxx.

References genie::units::A, DensityGaus(), DensityWoodsSaxon(), LOG, and pINFO.

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

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

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

Definition at line 442 of file NuclearUtils.cxx.

References genie::units::b, LOG, and pINFO.

Referenced by Density().

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

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

Definition at line 467 of file NuclearUtils.cxx.

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

Referenced by Density().

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

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

Definition at line 371 of file NuclearUtils.cxx.

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

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

double genie::utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization

(

const Target &

target

)

Definition at line 501 of file NuclearUtils.cxx.

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

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

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

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

Definition at line 361 of file NuclearUtils.cxx.

References genie::constants::kPi.

Referenced by RQEFG_generic().

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

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.

References FmI2().

Referenced by RQEFG_generic().

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

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.

References genie::units::b.

Referenced by FmI1(), and RQEFG_generic().

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

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

Definition at line 140 of file NuclearUtils.cxx.

References genie::Target::A(), genie::units::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::units::hbarc, genie::Target::HitNucP4Ptr(), genie::Target::HitNucPdg(), genie::Target::HitNucPosition(), genie::Interaction::InitState(), genie::Interaction::InitStatePtr(), genie::NuclearData::Instance(), genie::FermiMomentumTablePool::Instance(), genie::AlgFactory::Instance(), genie::AlgConfigPool::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::AhrensNCELPXSec::XSec(), genie::RosenbluthPXSec::XSec(), genie::AhrensDMELPXSec::XSec(), and genie::LwlynSmithQELCCPXSec::XSec().

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

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

Definition at line 130 of file NuclearUtils.cxx.

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

{
// Compute the nuclear radius (in fm)

  if(A<1) return 0;

  double Rn = Ro*TMath::Power(1.*A, 0.3333333);
  return Rn;
}

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

Definition at line 222 of file NuclearUtils.cxx.

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

Referenced by main(), and NuclQELXSecSuppression().

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