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 <costas.andreopoulos stfc.ac.uk> University of Liverpool & STFC Rutherford Appleton Lab

May 06, 2004

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

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, makeTrainCVSamples::int, genie::Target::IsNucleus(), Z, 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:71
bool IsNucleus(void) const
Definition: Target.cxx:289
Float_t Z
Definition: plot.C:38
double BindEnergy(const Target &target)
int Z(void) const
Definition: Target.h:69
static const double A
Definition: Units.h:82
double genie::utils::nuclear::BindEnergy ( int  nucA,
int  nucZ 
)

Definition at line 72 of file NuclearUtils.cxx.

References genie::units::A, a, b, d, delta, e, ana::Sqrt(), and Z.

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  double a = 15.835;
80  double b = 18.33;
81  double s = 23.20;
82  double d = 0.714;
83 
84  double delta = 0; /*E-O*/
85  if (nucZ%2==1 && nucA%2==1) delta = 11.2; /*O-O*/
86  if (nucZ%2==0 && nucA%2==0) delta = -11.2; /*E-E*/
87 
88  double N = (double) (nucA-nucZ);
89  double Z = (double) nucZ;
90  double A = (double) nucA;
91 
92  double BE = a * A
93  - b * TMath::Power(A,0.667)
94  - s * TMath::Power(N-Z,2.0)/A
95  - d * TMath::Power(Z,2.0)/TMath::Power(A,0.333)
96  - delta / TMath::Sqrt(A); // MeV
97 
98  return (1e-3 * BE); // GeV
99 }
double delta
Definition: runWimpSim.h:98
Float_t Z
Definition: plot.C:38
const XML_Char * s
Definition: expat.h:262
const double a
Float_t d
Definition: plot.C:236
static const double A
Definition: Units.h:82
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
const hit & b
Definition: hits.cxx:21
Float_t e
Definition: plot.C:35
double genie::utils::nuclear::BindEnergyLastNucleon ( const Target target)

Definition at line 110 of file NuclearUtils.cxx.

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

111 {
112 // Compute the binding for the most loose nucleon (in GeV)
113 
114  if(!target.IsNucleus()) return 0;
115 
116  //-- temporarily, return the binding energy per nucleon rather than the
117  // separation energy of the last nucleon
118  return (utils::nuclear::BindEnergy(target) / target.A());
119 }
int A(void) const
Definition: Target.h:71
bool IsNucleus(void) const
Definition: Target.cxx:289
double BindEnergy(const Target &target)
double genie::utils::nuclear::BindEnergyPerNucleon ( const Target target)

Definition at line 101 of file NuclearUtils.cxx.

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

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

102 {
103 // Compute the average binding energy per nucleon (in GeV)
104 
105  if(!target.IsNucleus()) return 0;
106 
107  return (utils::nuclear::BindEnergy(target) / target.A());
108 }
int A(void) const
Definition: Target.h:71
bool IsNucleus(void) const
Definition: Target.cxx:289
double BindEnergy(const Target &target)
double genie::utils::nuclear::BindEnergyPerNucleonParametrization ( const Target target)

Definition at line 484 of file NuclearUtils.cxx.

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

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

485 {
486 // Compute the average binding energy per nucleon (in GeV)
487 if(!target.IsNucleus()) return 0;
488 double x = TMath::Power(target.A(),1/3.0) / target.Z();
489 return (0.05042772591+x*(-0.11377355795+x*(0.15159890400-0.08825307197*x)));
490 }
int A(void) const
Definition: Target.h:71
bool IsNucleus(void) const
Definition: Target.cxx:289
int Z(void) const
Definition: Target.h:69
double genie::utils::nuclear::Density ( double  r,
int  A,
double  ring = 0. 
)

Definition at line 387 of file NuclearUtils.cxx.

References genie::units::A, make_syst_table_plots::c, DensityGaus(), DensityWoodsSaxon(), LOG, pINFO, ana::Sqrt(), and test::z.

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

388 {
389 // [by S.Dytman]
390 //
391  if(A>20) {
392  double c = 1., z = 1.;
393 
394  if (A == 27) { c = 3.07; z = 0.52; } // aluminum
395  else if (A == 28) { c = 3.07; z = 0.54; } // silicon
396  else if (A == 40) { c = 3.53; z = 0.54; } // argon
397  else if (A == 56) { c = 4.10; z = 0.56; } // iron
398  else if (A == 208) { c = 6.62; z = 0.55; } // lead
399  else {
400  c = TMath::Power(A,0.35); z = 0.54;
401  } //others
402 
403  LOG("Nuclear",pINFO)
404  << "r= " << r << ", ring= " << ring;
405  double rho = DensityWoodsSaxon(r,c,z,ring);
406  return rho;
407  }
408  else if (A>4) {
409  double ap = 1., alf = 1.;
410 
411  if (A == 7) { ap = 1.77; alf = 0.327; } // lithium
412  else if (A == 12) { ap = 1.69; alf = 1.08; } // carbon
413  else if (A == 14) { ap = 1.76; alf = 1.23; } // nitrogen
414  else if (A == 16) { ap = 1.83; alf = 1.54; } // oxygen
415  else {
416  ap=1.75; alf=-0.4+.12*A;
417  } //others- alf=0.08 if A=4
418 
419  double rho = DensityGaus(r,ap,alf,ring);
420  return rho;
421  }
422  else {
423  // helium
424  double ap = 1.9/TMath::Sqrt(2.);
425  double alf=0.;
426  double rho = DensityGaus(r,ap,alf,ring);
427  return rho;
428  }
429 
430  return 0;
431 }
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:97
#define pINFO
Definition: Messenger.h:63
z
Definition: test.py:28
double DensityWoodsSaxon(double r, double c, double z, double ring=0.)
static const double A
Definition: Units.h:82
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
TRandom3 r(0)
double genie::utils::nuclear::DensityGaus ( double  r,
double  ap,
double  alf,
double  ring = 0. 
)

Definition at line 433 of file NuclearUtils.cxx.

References b, LOG, norm, and pINFO.

Referenced by Density().

435 {
436 // [adapted from neugen3 density_gaus.F written by S.Dytman]
437 //
438 // Modified harmonic osc density distribution.
439 // Norm gives normalization to 1
440 //
441 // input : radial distance in nucleus [units: fm]
442 // output : nuclear density [units: fm^-3]
443 
444  ring = TMath::Min(ring, 0.3*a);
445 
446  double aeval = a + ring;
447  double norm = 1./((5.568 + alf*8.353)*TMath::Power(a,3.)); //0.0132;
448  double b = TMath::Power(r/aeval, 2.);
449  double dens = norm * (1. + alf*b) * TMath::Exp(-b);
450 
451  LOG("Nuclear", pINFO)
452  << "r = " << r << ", norm = " << norm << ", dens = " << dens
453  << ", aeval= " << aeval;
454 
455  return dens;
456 }
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const double a
#define pINFO
Definition: Messenger.h:63
Float_t norm
const hit & b
Definition: hits.cxx:21
TRandom3 r(0)
double genie::utils::nuclear::DensityWoodsSaxon ( double  r,
double  c,
double  z,
double  ring = 0. 
)

Definition at line 458 of file NuclearUtils.cxx.

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

Referenced by Density().

460 {
461 // [adapted from neugen3 density_ws.F written by S.Dytman]
462 //
463 // Woods-Saxon desity distribution. Norn gives normalization to 1
464 //
465 // input : radial distance in nucleus [units: fm]
466 // output : nuclear density [units: fm^-3]
467 
468  LOG("Nuclear",pINFO)
469  << "c= " << c << ", z= " << z << ",ring= " << ring;
470 
471  ring = TMath::Min(ring, 0.75*c);
472 
473  double ceval = c + ring;
474  double norm = (3./(4.*kPi*TMath::Power(c,3)))*1./(1.+TMath::Power((kPi*z/c),2));
475  double dens = norm / (1 + TMath::Exp((r-ceval)/z));
476 
477  LOG("Nuclear", pINFO)
478  << "r = " << r << ", norm = " << norm
479  << ", dens = " << dens << " , ceval= " << ceval;
480 
481  return dens;
482 }
const double kPi
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
#define pINFO
Definition: Messenger.h:63
z
Definition: test.py:28
Float_t norm
TRandom3 r(0)
double genie::utils::nuclear::DISNuclFactor ( double  x,
int  A 
)

Definition at line 362 of file NuclearUtils.cxx.

References MakeMiniprodValidationCuts::f.

Referenced by genie::QPMDISStrucFuncBase::NuclMod().

363 {
364 // Adapted from NeuGEN's nuc_factor(). Kept original comments from Hugh.
365 
366  double xv = TMath::Min(0.75, x);
367  double xv2 = xv * xv;
368  double xv3 = xv2 * xv;
369  double xv4 = xv3 * xv;
370  double xv5 = xv4 * xv;
371  double xvp = TMath::Power(xv, 14.417);
372  double expaxv = TMath::Exp(-21.94*xv);
373 
374  double f = 1.;
375 
376  // first factor goes from free nucleons to deuterium
377  if(A >= 2) {
378  f= 0.985*(1.+0.422*xv - 2.745*xv2 + 7.570*xv3 - 10.335*xv4 + 5.422*xv5);
379  }
380  // 2nd factor goes from deuterium to iso-scalar iron
381  if(A > 2) {
382  f *= (1.096 - 0.364*xv - 0.278*expaxv + 2.722*xvp);
383  }
384  return f;
385 }
static const double A
Definition: Units.h:82
double genie::utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization ( const Target target)

Definition at line 492 of file NuclearUtils.cxx.

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

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

493 {
494 // Compute Fermi momentum for isoscalar nucleon (in GeV)
495 if(!target.IsNucleus()) return 0;
496 double x = 1.0 / target.A();
497 return (0.27+x*(-1.12887857491+x*(9.72670908033-39.53095724456*x)));
498 }
int A(void) const
Definition: Target.h:71
bool IsNucleus(void) const
Definition: Target.cxx:289
double genie::utils::nuclear::FmArea ( double  alpha,
double  beta,
double  kf,
double  pmax 
)

Definition at line 352 of file NuclearUtils.cxx.

References genie::constants::kPi, pmax, and sum.

Referenced by RQEFG_generic().

354 {
355 // Adapted from NeuGEN's fm_area() used in r_factor()
356 
357  double kf3 = TMath::Power(kf,3.);
358  double sum = 4.*kPi* (alpha * kf3/3. + beta*(1./kf - 1./pmax));
359  return sum;
360 }
const int pmax
Definition: cellShifts.C:13
const double kPi
Double_t beta
Double_t sum
Definition: plot.C:31
double genie::utils::nuclear::FmI1 ( double  alpha,
double  beta,
double  a,
double  b,
double  kFi,
double  kFf,
double  q 
)

Definition at line 289 of file NuclearUtils.cxx.

References a2, a4, MakeMiniprodValidationCuts::f, FmI2(), and q2.

Referenced by RQEFG_generic().

291 {
292 // Adapted from NeuGEN's fm_integral1() used in r_factor()
293 
294  double f=0;
295 
296  double q2 = TMath::Power(q, 2);
297  double a2 = TMath::Power(a, 2);
298  double b2 = TMath::Power(b, 2);
299  double kFi2 = TMath::Power(kFi,2);
300  double kFf2 = TMath::Power(kFf,2);
301 
302  if(kFi < a) {
303  double lg = TMath::Log(b/a);
304 
305  f = -beta * (kFf2-q2)/(4.*q) * (1./a2 - 1./b2) + beta/(2.*q) * lg;
306 
307  } else if (kFi > b) {
308  double a4 = TMath::Power(a2,2);
309  double b4 = TMath::Power(b2,2);
310 
311  f = - (kFf2-q2) * alpha/(4.*q) * (b2-a2) + alpha/(8.*q) * (b4-a4);
312 
313  } else {
314  double a4 = TMath::Power(a2,2);
315  double kFi4 = TMath::Power(kFi2,2);
316  double lg = TMath::Log(b/kFi);
317 
318  f = - alpha*(kFf2-q2)/(4.*q)*(kFi2-a2) + alpha/(8.*q)*(kFi4 - a4)
319  - beta*(kFf2-q2)/(4.*q)*(1./kFi2 - 1./b2) + beta/(2.*q)*lg;
320  }
321 
322  double integral2 = FmI2(alpha,beta,a,b,kFi,kFf,q);
323  double integral1 = integral2 + f;
324 
325  return integral1;
326 }
TH1F * a2
Definition: f2_nu.C:545
double FmI2(double alpha, double beta, double a, double b, double kFi, double kFf, double q)
Double_t beta
Double_t q2[12][num]
Definition: f2_nu.C:137
const double a
TH1F * a4
Definition: f2_nu.C:666
const hit & b
Definition: hits.cxx:21
double genie::utils::nuclear::FmI2 ( double  alpha,
double  beta,
double  a,
double  b,
double  kFi,
double  kFf,
double  q 
)

Definition at line 328 of file NuclearUtils.cxx.

References a3, and b.

Referenced by FmI1(), and RQEFG_generic().

330 {
331 // Adapted from NeuGEN's fm_integral2() used in r_factor()
332 
333  double integral2 = 0;
334 
335  if(kFi < a) {
336  integral2 = beta * (1./a - 1./b);
337 
338  } else if(kFi > b) {
339  double a3 = TMath::Power(a,3);
340  double b3 = TMath::Power(b,3);
341  integral2 = alpha/3. * (b3 - a3);
342 
343  } else {
344  double a3 = TMath::Power(a,3);
345  double kFi3 = TMath::Power(kFi,3);
346 
347  integral2 = alpha/3. * (kFi3 - a3) + beta * (1./kFi - 1./b);
348  }
349  return integral2;
350 }
TH1F * a3
Definition: f2_nu.C:606
Double_t beta
const double a
const hit & b
Definition: hits.cxx:21
double genie::utils::nuclear::NuclQELXSecSuppression ( string  kftable,
double  pmax,
const Interaction in 
)

Definition at line 131 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(), python.hepunit::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(), q2, R, radius, 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().

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

Definition at line 121 of file NuclearUtils.cxx.

Referenced by genie::masterclass::MCTruthDisplay::DrawDiagram(), genie::INukeDeltaPropg::ProcessEventRecord(), and jmshower::JMShower::~JMShower().

122 {
123 // Compute the nuclear radius (in fm)
124 
125  if(A<1) return 0;
126 
127  double Rn = Ro*TMath::Power(1.*A, 0.3333333);
128  return Rn;
129 }
static const double A
Definition: Units.h:82
double genie::utils::nuclear::RQEFG_generic ( double  q2,
double  Mn,
double  kFi,
double  kFf,
double  pmax 
)

Definition at line 213 of file NuclearUtils.cxx.

References beta, exit(), FmArea(), FmI1(), FmI2(), genie::constants::kPi, LOG, plot_validation_datamc::p1, plot_validation_datamc::p2, pFATAL, R, and ana::Sqrt().

Referenced by main(), and NuclQELXSecSuppression().

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