GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
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().

60 {
61 // Compute the average binding energy (in GeV) using the semi-empirical
62 // formula from Wapstra (Handbuch der Physik, XXXVIII/1)
63 
64  if(!target.IsNucleus()) return 0;
65 
66  int Z = (int) target.Z();
67  int A = (int) target.A();
68 
69  return BindEnergy(A,Z);
70 }
int A(void) const
Definition: Target.h:70
bool IsNucleus(void) const
Definition: Target.cxx:272
double BindEnergy(const Target &target)
static constexpr double A
Definition: Units.h:74
int Z(void) const
Definition: Target.h:68
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.

73 {
74 // Compute the average binding energy (in GeV) using the semi-empirical
75 // formula from Wapstra (Handbuch der Physik, XXXVIII/1)
76 
77  if (nucA<=0 || nucZ<=0) return 0;
78 
79  if (nucZ == 2)
80  {
81  if (nucA == 3)
82  return 7.5e-3;
83  if (nucA == 4)
84  return 28.4e-3;
85  return 0.;
86  }
87 
88  double a = 15.835;
89  double b = 18.33;
90  double s = 23.20;
91  double d = 0.714;
92 
93  double delta = 0; /*E-O*/
94  if (nucZ%2==1 && nucA%2==1) delta = 11.2; /*O-O*/
95  if (nucZ%2==0 && nucA%2==0) delta = -11.2; /*E-E*/
96 
97  double N = (double) (nucA-nucZ);
98  double Z = (double) nucZ;
99  double A = (double) nucA;
100 
101  double BE = a * A
102  - b * TMath::Power(A,0.667)
103  - s * TMath::Power(N-Z,2.0)/A
104  - d * TMath::Power(Z,2.0)/TMath::Power(A,0.333)
105  - delta / TMath::Sqrt(A); // MeV
106 
107  return (1e-3 * BE); // GeV
108 }
static constexpr double s
Definition: Units.h:95
static constexpr double b
Definition: Units.h:78
const double e
static constexpr double A
Definition: Units.h:74
const double a
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().

120 {
121 // Compute the binding for the most loose nucleon (in GeV)
122 
123  if(!target.IsNucleus()) return 0;
124 
125  //-- temporarily, return the binding energy per nucleon rather than the
126  // separation energy of the last nucleon
127  return (utils::nuclear::BindEnergy(target) / target.A());
128 }
int A(void) const
Definition: Target.h:70
bool IsNucleus(void) const
Definition: Target.cxx:272
double BindEnergy(const Target &target)
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().

111 {
112 // Compute the average binding energy per nucleon (in GeV)
113 
114  if(!target.IsNucleus()) return 0;
115 
116  return (utils::nuclear::BindEnergy(target) / target.A());
117 }
int A(void) const
Definition: Target.h:70
bool IsNucleus(void) const
Definition: Target.cxx:272
double BindEnergy(const Target &target)
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().

494 {
495 // Compute the average binding energy per nucleon (in GeV)
496 if(!target.IsNucleus()) return 0;
497 double x = TMath::Power(target.A(),1/3.0) / target.Z();
498 return (0.05042772591+x*(-0.11377355795+x*(0.15159890400-0.08825307197*x)));
499 }
int A(void) const
Definition: Target.h:70
bool IsNucleus(void) const
Definition: Target.cxx:272
int Z(void) const
Definition: Target.h:68
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().

397 {
398 // [by S.Dytman]
399 //
400  if(A>20) {
401  double c = 1., z = 1.;
402 
403  if (A == 27) { c = 3.07; z = 0.52; } // aluminum
404  else if (A == 28) { c = 3.07; z = 0.54; } // silicon
405  else if (A == 40) { c = 3.53; z = 0.54; } // argon
406  else if (A == 56) { c = 4.10; z = 0.56; } // iron
407  else if (A == 208) { c = 6.62; z = 0.55; } // lead
408  else {
409  c = TMath::Power(A,0.35); z = 0.54;
410  } //others
411 
412  LOG("Nuclear",pINFO)
413  << "r= " << r << ", ring= " << ring;
414  double rho = DensityWoodsSaxon(r,c,z,ring);
415  return rho;
416  }
417  else if (A>4) {
418  double ap = 1., alf = 1.;
419 
420  if (A == 7) { ap = 1.77; alf = 0.327; } // lithium
421  else if (A == 12) { ap = 1.69; alf = 1.08; } // carbon
422  else if (A == 14) { ap = 1.76; alf = 1.23; } // nitrogen
423  else if (A == 16) { ap = 1.83; alf = 1.54; } // oxygen
424  else {
425  ap=1.75; alf=-0.4+.12*A;
426  } //others- alf=0.08 if A=4
427 
428  double rho = DensityGaus(r,ap,alf,ring);
429  return rho;
430  }
431  else {
432  // helium
433  double ap = 1.9/TMath::Sqrt(2.);
434  double alf=0.;
435  double rho = DensityGaus(r,ap,alf,ring);
436  return rho;
437  }
438 
439  return 0;
440 }
double DensityGaus(double r, double ap, double alf, double ring=0.)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static constexpr double A
Definition: Units.h:74
#define pINFO
Definition: Messenger.h:62
double DensityWoodsSaxon(double r, double c, double z, double ring=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().

444 {
445 // [adapted from neugen3 density_gaus.F written by S.Dytman]
446 //
447 // Modified harmonic osc density distribution.
448 // Norm gives normalization to 1
449 //
450 // input : radial distance in nucleus [units: fm]
451 // output : nuclear density [units: fm^-3]
452 
453  ring = TMath::Min(ring, 0.3*a);
454 
455  double aeval = a + ring;
456  double norm = 1./((5.568 + alf*8.353)*TMath::Power(a,3.)); //0.0132;
457  double b = TMath::Power(r/aeval, 2.);
458  double dens = norm * (1. + alf*b) * TMath::Exp(-b);
459 
460  LOG("Nuclear", pINFO)
461  << "r = " << r << ", norm = " << norm << ", dens = " << dens
462  << ", aeval= " << aeval;
463 
464  return dens;
465 }
static constexpr double b
Definition: Units.h:78
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const double a
#define pINFO
Definition: Messenger.h:62
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().

469 {
470 // [adapted from neugen3 density_ws.F written by S.Dytman]
471 //
472 // Woods-Saxon desity distribution. Norn gives normalization to 1
473 //
474 // input : radial distance in nucleus [units: fm]
475 // output : nuclear density [units: fm^-3]
476 
477  LOG("Nuclear",pINFO)
478  << "c= " << c << ", z= " << z << ",ring= " << ring;
479 
480  ring = TMath::Min(ring, 0.75*c);
481 
482  double ceval = c + ring;
483  double norm = (3./(4.*kPi*TMath::Power(c,3)))*1./(1.+TMath::Power((kPi*z/c),2));
484  double dens = norm / (1 + TMath::Exp((r-ceval)/z));
485 
486  LOG("Nuclear", pINFO)
487  << "r = " << r << ", norm = " << norm
488  << ", dens = " << dens << " , ceval= " << ceval;
489 
490  return dens;
491 }
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pINFO
Definition: Messenger.h:62
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().

372 {
373 // Adapted from NeuGEN's nuc_factor(). Kept original comments from Hugh.
374 
375  double xv = TMath::Min(0.75, x);
376  double xv2 = xv * xv;
377  double xv3 = xv2 * xv;
378  double xv4 = xv3 * xv;
379  double xv5 = xv4 * xv;
380  double xvp = TMath::Power(xv, 14.417);
381  double expaxv = TMath::Exp(-21.94*xv);
382 
383  double f = 1.;
384 
385  // first factor goes from free nucleons to deuterium
386  if(A >= 2) {
387  f= 0.985*(1.+0.422*xv - 2.745*xv2 + 7.570*xv3 - 10.335*xv4 + 5.422*xv5);
388  }
389  // 2nd factor goes from deuterium to iso-scalar iron
390  if(A > 2) {
391  f *= (1.096 - 0.364*xv - 0.278*expaxv + 2.722*xvp);
392  }
393  return f;
394 }
static constexpr double A
Definition: Units.h:74
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().

502 {
503 // Compute Fermi momentum for isoscalar nucleon (in GeV)
504 if(!target.IsNucleus()) return 0;
505 double x = 1.0 / target.A();
506 return (0.27+x*(-1.12887857491+x*(9.72670908033-39.53095724456*x)));
507 }
int A(void) const
Definition: Target.h:70
bool IsNucleus(void) const
Definition: Target.cxx:272
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().

363 {
364 // Adapted from NeuGEN's fm_area() used in r_factor()
365 
366  double kf3 = TMath::Power(kf,3.);
367  double sum = 4.*kPi* (alpha * kf3/3. + beta*(1./kf - 1./pmax));
368  return sum;
369 }
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().

300 {
301 // Adapted from NeuGEN's fm_integral1() used in r_factor()
302 
303  double f=0;
304 
305  double q2 = TMath::Power(q, 2);
306  double a2 = TMath::Power(a, 2);
307  double b2 = TMath::Power(b, 2);
308  double kFi2 = TMath::Power(kFi,2);
309  double kFf2 = TMath::Power(kFf,2);
310 
311  if(kFi < a) {
312  double lg = TMath::Log(b/a);
313 
314  f = -beta * (kFf2-q2)/(4.*q) * (1./a2 - 1./b2) + beta/(2.*q) * lg;
315 
316  } else if (kFi > b) {
317  double a4 = TMath::Power(a2,2);
318  double b4 = TMath::Power(b2,2);
319 
320  f = - (kFf2-q2) * alpha/(4.*q) * (b2-a2) + alpha/(8.*q) * (b4-a4);
321 
322  } else {
323  double a4 = TMath::Power(a2,2);
324  double kFi4 = TMath::Power(kFi2,2);
325  double lg = TMath::Log(b/kFi);
326 
327  f = - alpha*(kFf2-q2)/(4.*q)*(kFi2-a2) + alpha/(8.*q)*(kFi4 - a4)
328  - beta*(kFf2-q2)/(4.*q)*(1./kFi2 - 1./b2) + beta/(2.*q)*lg;
329  }
330 
331  double integral2 = FmI2(alpha,beta,a,b,kFi,kFf,q);
332  double integral1 = integral2 + f;
333 
334  return integral1;
335 }
double FmI2(double alpha, double beta, double a, double b, double kFi, double kFf, double q)
static constexpr double b
Definition: Units.h:78
const double a
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().

339 {
340 // Adapted from NeuGEN's fm_integral2() used in r_factor()
341 
342  double integral2 = 0;
343 
344  if(kFi < a) {
345  integral2 = beta * (1./a - 1./b);
346 
347  } else if(kFi > b) {
348  double a3 = TMath::Power(a,3);
349  double b3 = TMath::Power(b,3);
350  integral2 = alpha/3. * (b3 - a3);
351 
352  } else {
353  double a3 = TMath::Power(a,3);
354  double kFi3 = TMath::Power(kFi,3);
355 
356  integral2 = alpha/3. * (kFi3 - a3) + beta * (1./kFi - 1./b);
357  }
358  return integral2;
359 }
static constexpr double b
Definition: Units.h:78
const double a
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().

142 {
143  const InitialState & init_state = interaction->InitState();
144  const ProcessInfo & proc_info = interaction->ProcInfo();
145  const Target & target = init_state.Tgt();
146 
147  if (!target.IsNucleus()) {
148  return 1.0; // no suppression for free nucleon tgt
149  }
150 
151  //
152  // special case: deuterium
153  //
154 
155  if (target.A() == 2) {
156  NuclearData * nucldat = NuclearData::Instance();
157  return nucldat->DeuteriumSuppressionFactor(interaction->Kine().Q2());
158  }
159 
160  //
161  // general case
162  //
163 
164  int target_pdgc = target.Pdg();
165  int struck_nucleon_pdgc = target.HitNucPdg();
166  int final_nucleon_pdgc = struck_nucleon_pdgc;
167 
168  if(proc_info.IsWeakCC()) {
169  final_nucleon_pdgc = pdg::SwitchProtonNeutron(struck_nucleon_pdgc);
170  }
171 
172  // Check if an LFG model should be used for Fermi momentum
173  // Create a nuclear model object to check the model type
174  AlgConfigPool * confp = AlgConfigPool::Instance();
175  const Registry * gc = confp->GlobalParameterList();
176  RgKey nuclkey = "NuclearModel";
177  RgAlg nuclalg = gc->GetAlg(nuclkey);
178  AlgFactory * algf = AlgFactory::Instance();
179  const genie::NuclearModelI* nuclModel =
180  dynamic_cast<const genie::NuclearModelI*>(
181  algf->GetAlgorithm(nuclalg.name,nuclalg.config));
182  // Check if the model is a local Fermi gas
183  bool lfg = (nuclModel && nuclModel->ModelType(Target()) == kNucmLocalFermiGas);
184 
185  double kFi, kFf;
186  if(lfg){
188  Target* tgt = interaction->InitStatePtr()->TgtPtr();
189  double radius = tgt->HitNucPosition();
190 
191  int A = tgt->A();
192  // kFi
193  bool is_p_i = pdg::IsProton(struck_nucleon_pdgc);
194  double numNuci = (is_p_i) ? (double)tgt->Z():(double)tgt->N();
195  kFi = TMath::Power(3*kPi2*numNuci*
196  genie::utils::nuclear::Density(radius,A),1.0/3.0) *hbarc;
197  // kFi
198  bool is_p_f = pdg::IsProton(final_nucleon_pdgc);
199  double numNucf = (is_p_f) ? (double)tgt->Z():(double)tgt->N();
200  kFf = TMath::Power(3*kPi2*numNucf*
201  genie::utils::nuclear::Density(radius,A),1.0/3.0) *hbarc;
202  }else{
203  // get the requested Fermi momentum table
204  FermiMomentumTablePool * kftp = FermiMomentumTablePool::Instance();
205  const FermiMomentumTable * kft = kftp->GetTable(kftable);
206 
207  // Fermi momentum for initial, final nucleons
208  kFi = kft->FindClosestKF(target_pdgc, struck_nucleon_pdgc);
209  kFf = (struck_nucleon_pdgc==final_nucleon_pdgc) ? kFi :
210  kft->FindClosestKF(target_pdgc, final_nucleon_pdgc );
211  }
212 
213  double Mn = target.HitNucP4Ptr()->M(); // can be off m/shell
214 
215  const Kinematics & kine = interaction->Kine();
216  double q2 = kine.q2();
217 
218  double R = RQEFG_generic(q2, Mn, kFi, kFf, pmax);
219  return R;
220 }
bool IsWeakCC(void) const
static constexpr double hbarc
Definition: Units.h:34
int HitNucPdg(void) const
Definition: Target.cxx:304
double Density(double r, int A, double ring=0.)
int A(void) const
Definition: Target.h:70
A singleton class holding all configuration registries built while parsing all loaded XML configurati...
Definition: AlgConfigPool.h:40
bool IsNucleus(void) const
Definition: Target.cxx:272
int Pdg(void) const
Definition: Target.h:71
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
int SwitchProtonNeutron(int pdgc)
Definition: PDGUtils.cxx:356
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
A table of Fermi momentum constants.
virtual NuclearModel_t ModelType(const Target &) const =0
double RQEFG_generic(double q2, double Mn, double kFi, double kFf, double pmax)
double q2(bool selected=false) const
Definition: Kinematics.cxx:141
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
Singleton class to load &amp; serve tables of Fermi momentum constants.
static constexpr double A
Definition: Units.h:74
const FermiMomentumTable * GetTable(string name)
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:75
int Z(void) const
Definition: Target.h:68
Registry * GlobalParameterList(void) const
int N(void) const
Definition: Target.h:69
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
double HitNucPosition(void) const
Definition: Target.h:89
double DeuteriumSuppressionFactor(double Q2)
Definition: NuclearData.cxx:65
static constexpr double fermi
Definition: Units.h:55
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
const Target & Tgt(void) const
Definition: InitialState.h:66
The GENIE Algorithm Factory.
Definition: AlgFactory.h:39
static const double kPlankConstant
string config
RgAlg GetAlg(RgKey key) const
Definition: Registry.cxx:488
Initial State information.
Definition: InitialState.h:48
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().

131 {
132 // Compute the nuclear radius (in fm)
133 
134  if(A<1) return 0;
135 
136  double Rn = Ro*TMath::Power(1.*A, 0.3333333);
137  return Rn;
138 }
static constexpr double A
Definition: Units.h:74
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().

224 {
225 // Computes the nuclear suppression of the QEL differential cross section
226 //
227 // Inputs:
228 // - q2 : momentum transfer (< 0)
229 // - Mn : hit nucleon mass (nucleon can be off the mass shell)
230 // - kFi : Fermi momentum, initial state (hit) nucleon @ nuclear target
231 // - kFf : Fermi momentum, final state (recoil) nucleon @ nuclear target
232 // - pmax : A Fermi momentum cutoff
233 //
234 // A direct adaptation of NeuGEN's qelnuc() and r_factor()
235 //
236 // Hugh's comments from the original code:
237 // "This routine is based on an analytic calculation of the rejection factor
238 // in the Fermi Gas model using the form for the fermi momentum distribution
239 // given in the Bodek and Ritchie paper. [P.R. D23 (1981) 1070]
240 // R is the ratio of the differential cross section from the nuclear material
241 // specified by (kf,pf) to the differential cross section for a free nucleon".
242 // (kf,pf = Fermi Gas model Fermi momentum for initial,final nucleons)
243 //
244 
245  double R = 1.;
246 
247  // Compute magnitude of the 3-momentum transfer to the nucleon
248  double Mn2 = Mn*Mn;
249  double magq2 = q2 * (0.25*q2/Mn2 - 1.);
250  double q = TMath::Sqrt(TMath::Max(0.,magq2));
251 
252  double kfa = kFi * 2./kPi;
253  double kfa2 = TMath::Power(kfa,2);
254  double kFi4 = TMath::Power(kFi,4);
255  double rkf = 1./(1. - kFi/4.);
256  double alpha = 1. - 6.*kfa2;
257  double beta = 2. * rkf * kfa2 * kFi4;
258 
259  double fm_area = FmArea(alpha,beta,kFi,pmax);
260 
261  if (q <= kFf) {
262 
263  double p1 = kFf - q;
264  double p2 = kFf + q;
265  double fmi1 = FmI1(alpha,beta,p1,p2, kFi,kFf,q);
266  double fmi2 = FmI2(alpha,beta,p2,pmax,kFi,kFf,q);
267 
268  R = 2*kPi * (fmi1 + 2*fmi2) / fm_area;
269 
270  } else if (q > kFf && q <= (pmax-kFf)) {
271 
272  double p1 = q - kFf;
273  double p2 = q + kFf;
274  double fmi1 = FmI1(alpha,beta, p1, p2, kFi, kFf,q);
275  double fmi2a = FmI2(alpha,beta, 0., p1, kFi, kFf,q);
276  double fmi2b = FmI2(alpha,beta, p2, pmax, kFi, kFf,q);
277 
278  R = 2*kPi * (fmi1 + 2*(fmi2a+fmi2b)) / fm_area;
279 
280  } else if (q > (pmax-kFf) && q <= (pmax+kFf)) {
281 
282  double p1 = q - kFf;
283  double fmi2 = FmI2(alpha,beta, 0.,p1, kFi,kFf,q);
284  double fmi1 = FmI1(alpha,beta, p1,pmax,kFi,kFf,q);
285 
286  R = 2*kPi * (2.*fmi2 + fmi1) / fm_area;
287 
288  } else if (q > (pmax+kFf)) {
289  R = 1.;
290 
291  } else {
292  LOG("Nuclear", pFATAL) << "Illegal input q = " << q;
293  exit(1);
294  }
295  return R;
296 }
double FmArea(double alpha, double beta, double kf, double pmax)
double FmI2(double alpha, double beta, double a, double b, double kFi, double kFf, double q)
#define pFATAL
Definition: Messenger.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double FmI1(double alpha, double beta, double a, double b, double kFi, double kFf, double q)