NuclearUtils.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
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 //___________________________________________________________________________
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 }
100 //___________________________________________________________________________
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 }
109 //___________________________________________________________________________
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 }
120 //___________________________________________________________________________
121 double genie::utils::nuclear::Radius(int A, double Ro)
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 }
130 //___________________________________________________________________________
132  string kftable, double pmax, const Interaction * interaction)
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
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
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 }
212 //___________________________________________________________________________
214  double q2, double Mn, double kFi, double kFf, double pmax)
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 }
288 //___________________________________________________________________________
289 double genie::utils::nuclear::FmI1(double alpha, double beta,
290  double a, double b, double kFi, double kFf, double q)
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 }
327 //___________________________________________________________________________
328 double genie::utils::nuclear::FmI2(double alpha, double beta,
329  double a, double b, double kFi, double /*kFf*/, double /*q*/)
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 }
351 //___________________________________________________________________________
353  double alpha, double beta, double kf, double pmax)
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 }
361 //___________________________________________________________________________
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 }
386 //___________________________________________________________________________
387 double genie::utils::nuclear::Density(double r, int A, double ring)
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 }
432 //___________________________________________________________________________
434  double r, double a, double alf, double ring)
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 }
457 //___________________________________________________________________________
459  double r, double c, double z, double ring)
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 }
483 //___________________________________________________________________________
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 }
491 //___________________________________________________________________________
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 }
499 //___________________________________________________________________________
500 
double DensityGaus(double r, double ap, double alf, double ring=0.)
const int pmax
Definition: cellShifts.C:13
Basic constants.
TH1F * a3
Definition: f2_nu.C:606
static NuclearData * Instance(void)
Definition: NuclearData.cxx:54
bool IsWeakCC(void) const
const XML_Char * target
Definition: expat.h:268
TH1F * a2
Definition: f2_nu.C:545
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
double FmArea(double alpha, double beta, double kf, double pmax)
int HitNucPdg(void) const
Definition: Target.cxx:321
double delta
Definition: runWimpSim.h:98
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)
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
#define pFATAL
Definition: Messenger.h:57
bool IsNucleus(void) const
Definition: Target.cxx:289
int Pdg(void) const
Definition: Target.h:72
static FermiMomentumTablePool * Instance(void)
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.
Double_t beta
virtual NuclearModel_t ModelType(const Target &) const =0
Float_t Z
Definition: plot.C:38
double Radius(int A, double Ro=constants::kNucRo)
Double_t q2[12][num]
Definition: f2_nu.C:137
double RQEFG_generic(double q2, double Mn, double kFi, double kFf, double pmax)
const XML_Char * s
Definition: expat.h:262
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:151
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
Singleton class to load & serve tables of Fermi momentum constants.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
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
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:41
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:86
Float_t d
Definition: plot.C:236
#define R(x)
int Z(void) const
Definition: Target.h:69
Double_t radius
#define pINFO
Definition: Messenger.h:63
double FmI1(double alpha, double beta, double a, double b, double kFi, double kFf, double q)
z
Definition: test.py:28
Registry * GlobalParameterList(void) const
static const double kLightSpeed
Definition: Constants.h:32
double BindEnergyLastNucleon(const Target &target)
double DensityWoodsSaxon(double r, double c, double z, double ring=0.)
double FermiMomentumForIsoscalarNucleonParametrization(const Target &target)
static const double A
Definition: Units.h:82
int N(void) const
Definition: Target.h:70
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
static AlgFactory * Instance()
Definition: AlgFactory.cxx:75
Float_t norm
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
double HitNucPosition(void) const
Definition: Target.h:90
Target * TgtPtr(void) const
Definition: InitialState.h:68
TH1F * a4
Definition: f2_nu.C:666
exit(0)
const hit & b
Definition: hits.cxx:21
double DeuteriumSuppressionFactor(double Q2)
Definition: NuclearData.cxx:65
TRandom3 r(0)
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:135
const Target & Tgt(void) const
Definition: InitialState.h:67
The GENIE Algorithm Factory.
Definition: AlgFactory.h:40
Double_t sum
Definition: plot.C:31
Float_t e
Definition: plot.C:35
double DISNuclFactor(double x, int A)
static const double kPlankConstant
Definition: Constants.h:33
static const double kPi
Definition: Constants.h:38
static const double kPi2
Definition: Constants.h:39
string config
RgAlg GetAlg(RgKey key) const
Definition: Registry.cxx:503
static AlgConfigPool * Instance()
double NuclQELXSecSuppression(string kftable, double pmax, const Interaction *in)
Initial State information.
Definition: InitialState.h:49