GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NuclearUtils.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6 
7  Costas Andreopoulos <c.andreopoulos \at cern.ch>
8  University of Liverpool
9 
10  For documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13 
14  @ Sep 19, 2007 - CA
15  Copied here the nuclear density methods used by intranuke so that they
16  can be used by other modules as well (eg position vertex selection)
17  @ Nov 30, 2007 - CA
18  Added possibility to increase the nuclear size (intranuke may increase the
19  nuclear size by a const times the de Broglie wavelength of a transported
20  hadron). Density() methods have a new default argument (ring) which is 0
21  if not explicity set.
22  @ Nov 30, 2009 - CA
23  Overriding NuclQELXSecSuppression() calc for deuterium and using input
24  data from S.K.Singh, Nucl. Phys. B 36, 419 (1972) [data used by Hugh for
25  the Merenyi test]
26  @ May 18, 2010 - CA
27  Restructure NuclQELXSecSuppression() (Add RQEFG_generic()) to aid
28  reweighting.
29  @ Jul 15, 2010 - AM
30  Added BindEnergy(int nucA, int nucZ), used in Intranuke
31  @ Mar 18, 2016 - JJ (SD)
32  Check if a local Fermi gas model should be used when calculating the
33  Fermi momentum
34 */
35 //____________________________________________________________________________
36 
37 #include <cstdlib>
38 
39 #include <TMath.h>
40 
53 
54 using namespace genie;
55 using namespace genie::constants;
56 
57 
58 //____________________________________________________________________________
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 }
71 //___________________________________________________________________________
72 double genie::utils::nuclear::BindEnergy(int nucA, int nucZ)
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 }
109 //___________________________________________________________________________
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 }
118 //___________________________________________________________________________
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 }
129 //___________________________________________________________________________
130 double genie::utils::nuclear::Radius(int A, double Ro)
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 }
139 //___________________________________________________________________________
141  string kftable, double pmax, const Interaction * interaction)
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
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
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 }
221 //___________________________________________________________________________
223  double q2, double Mn, double kFi, double kFf, double pmax)
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 }
297 //___________________________________________________________________________
298 double genie::utils::nuclear::FmI1(double alpha, double beta,
299  double a, double b, double kFi, double kFf, double q)
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 }
336 //___________________________________________________________________________
337 double genie::utils::nuclear::FmI2(double alpha, double beta,
338  double a, double b, double kFi, double /*kFf*/, double /*q*/)
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 }
360 //___________________________________________________________________________
362  double alpha, double beta, double kf, double pmax)
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 }
370 //___________________________________________________________________________
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 }
395 //___________________________________________________________________________
396 double genie::utils::nuclear::Density(double r, int A, double ring)
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 }
441 //___________________________________________________________________________
443  double r, double a, double alf, double ring)
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 }
466 //___________________________________________________________________________
468  double r, double c, double z, double ring)
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 }
492 //___________________________________________________________________________
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 }
500 //___________________________________________________________________________
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 }
508 //___________________________________________________________________________
509 
double DensityGaus(double r, double ap, double alf, double ring=0.)
static NuclearData * Instance(void)
Definition: NuclearData.cxx:54
bool IsWeakCC(void) const
static constexpr double hbarc
Definition: Units.h:34
double FmArea(double alpha, double beta, double kf, double pmax)
int HitNucPdg(void) const
Definition: Target.cxx:304
double Density(double r, int A, double ring=0.)
double FmI2(double alpha, double beta, double a, double b, double kFi, double kFf, double q)
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
#define pFATAL
Definition: Messenger.h:56
bool IsNucleus(void) const
Definition: Target.cxx:272
int Pdg(void) const
Definition: Target.h:71
static constexpr double s
Definition: Units.h:95
static FermiMomentumTablePool * Instance(void)
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
static constexpr double b
Definition: Units.h:78
double Radius(int A, double Ro=constants::kNucRo)
double RQEFG_generic(double q2, double Mn, double kFi, double kFf, double pmax)
Summary information for an interaction.
Definition: Interaction.h:56
double BindEnergy(const Target &target)
double BindEnergyPerNucleon(const Target &target)
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.
const double e
#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
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
const double a
double BindEnergyPerNucleonParametrization(const Target &target)
const Kinematics & Kine(void) const
Definition: Interaction.h:71
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
#define pINFO
Definition: Messenger.h:62
double FmI1(double alpha, double beta, double a, double b, double kFi, double kFf, double q)
Registry * GlobalParameterList(void) const
double BindEnergyLastNucleon(const Target &target)
double DensityWoodsSaxon(double r, double c, double z, double ring=0.)
double FermiMomentumForIsoscalarNucleonParametrization(const Target &target)
int N(void) const
Definition: Target.h:69
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:247
static AlgFactory * Instance()
Definition: AlgFactory.cxx:64
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
double HitNucPosition(void) const
Definition: Target.h:89
Target * TgtPtr(void) const
Definition: InitialState.h:67
double DeuteriumSuppressionFactor(double Q2)
Definition: NuclearData.cxx:65
static constexpr double fermi
Definition: Units.h:55
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
const Target & Tgt(void) const
Definition: InitialState.h:66
The GENIE Algorithm Factory.
Definition: AlgFactory.h:39
double DISNuclFactor(double x, int A)
static const double kPlankConstant
string config
RgAlg GetAlg(RgKey key) const
Definition: Registry.cxx:488
static AlgConfigPool * Instance()
double NuclQELXSecSuppression(string kftable, double pmax, const Interaction *in)
Initial State information.
Definition: InitialState.h:48