LocalFGM.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: Joe Johnston, University of Pittsburgh (Advisor Steven Dytman)
8 
9  For the class documentation see the corresponding header file.
10 
11 */
12 //____________________________________________________________________________
13 
14 #include <sstream>
15 #include <cstdlib>
16 #include <TMath.h>
17 
28 
29 using std::ostringstream;
30 using namespace genie;
31 using namespace genie::constants;
32 using namespace genie::utils;
33 
34 //____________________________________________________________________________
36 NuclearModelI("genie::LocalFGM")
37 {
38 
39 }
40 //____________________________________________________________________________
42 NuclearModelI("genie::LocalFGM", config)
43 {
44 
45 }
46 //____________________________________________________________________________
48 {
49 
50 }
51 //____________________________________________________________________________
53  double hitNucleonRadius) const
54 {
55  assert(target.HitNucIsSet());
56 
58  fCurrMomentum.SetXYZ(0,0,0);
59 
60  //-- set fermi momentum vector
61  //
62  TH1D * prob = this->ProbDistro(target,hitNucleonRadius);
63  if(!prob) {
64  LOG("LocalFGM", pNOTICE)
65  << "Null nucleon momentum probability distribution";
66  exit(1);
67  }
68  double p = prob->GetRandom();
69  delete prob;
70  LOG("LocalFGM", pINFO) << "|p,nucleon| = " << p;
71 
73 
74  double costheta = -1. + 2. * rnd->RndGen().Rndm();
75  double sintheta = TMath::Sqrt(1.-costheta*costheta);
76  double fi = 2 * kPi * rnd->RndGen().Rndm();
77  double cosfi = TMath::Cos(fi);
78  double sinfi = TMath::Sin(fi);
79 
80  double px = p*sintheta*cosfi;
81  double py = p*sintheta*sinfi;
82  double pz = p*costheta;
83 
84  fCurrMomentum.SetXYZ(px,py,pz);
85 
86  //-- set removal energy
87  //
88  int Z = target.Z();
89  map<int,double>::const_iterator it = fNucRmvE.find(Z);
90  if(it != fNucRmvE.end()) fCurrRemovalEnergy = it->second;
91  else fCurrRemovalEnergy = nuclear::BindEnergyPerNucleon(target);
92 
93  return true;
94 }
95 //____________________________________________________________________________
96 double LocalFGM::Prob(double p, double w, const Target & target,
97  double hitNucleonRadius) const
98 {
99  if(w<0) {
100  TH1D * prob = this->ProbDistro(target, hitNucleonRadius);
101  int bin = prob->FindBin(p);
102  double y = prob->GetBinContent(bin);
103  double dx = prob->GetBinWidth(bin);
104  double pr = y*dx;
105  delete prob;
106  return pr;
107  }
108  return 1;
109 }
110 //____________________________________________________________________________
111 // *** The TH1D object must be deleted after it is used ***
112 TH1D * LocalFGM::ProbDistro(const Target & target, double r) const
113 {
114  LOG("LocalFGM", pNOTICE)
115  << "Computing P = f(p_nucleon) for: " << target.AsString()
116  << ", Nucleon Radius = " << r;
117  LOG("LocalFGM", pNOTICE)
118  << ", P(max) = " << fPMax;
119 
120  //-- get information for the nuclear target
121  int nucleon_pdgc = target.HitNucPdg();
122  assert(pdg::IsProton(nucleon_pdgc) || pdg::IsNeutron(nucleon_pdgc));
123  int A = target.A();
124 
125  assert(target.HitNucIsSet());
126  bool is_p = pdg::IsProton(nucleon_pdgc);
127  double numNuc = (is_p) ? (double)target.Z():(double)target.N();
128 
129  // Calculate Fermi Momentum using Local FG equations
131  double KF= TMath::Power(3*kPi2*numNuc*genie::utils::nuclear::Density(r,A),
132  1.0/3.0) *hbarc;
133 
134  LOG("LocalFGM",pNOTICE) << "KF = " << KF;
135 
136  double a = 2.0;
137  double C = 4. * kPi * TMath::Power(KF,3) / 3.;
138 
139  // Do not include nucleon correlation tail
140  //double R = 1. / (1.- KF/4.);
141 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
142  LOG("LocalFGM", pDEBUG) << "a = " << a;
143  LOG("LocalFGM", pDEBUG) << "C = " << C;
144  //LOG("LocalFGM", pDEBUG) << "R = " << R;
145 #endif
146 
147  //-- create the probability distribution
148 
149  int npbins = (int) (1000*fPMax);
150  TH1D * prob = new TH1D("", "", npbins, 0, fPMax);
151  prob->SetDirectory(0);
152 
153  double dp = fPMax / (npbins-1);
154  double iC = (C>0) ? 1./C : 0.;
155  double kfa_pi_2 = TMath::Power(KF*a/kPi,2);
156 
157  for(int i = 0; i < npbins; i++) {
158  double p = i * dp;
159  double p2 = TMath::Power(p,2);
160 
161  // calculate |phi(p)|^2
162  double phi2 = 0;
163  if (p <= KF)
164  phi2 = iC * (1. - 6.*kfa_pi_2);
165 
166  // Do not include nucleon correlation tail
167  //else if ( p > KF && p < fPCutOff)
168  // phi2 = iC * (2*R*kfa_pi_2*TMath::Power(KF/p,4.));
169 
170  // calculate probability density : dProbability/dp
171  double dP_dp = 4*kPi * p2 * phi2;
172 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
173  LOG("LocalFGM", pDEBUG) << "p = " << p << ", dP/dp = " << dP_dp;
174 #endif
175  prob->Fill(p, dP_dp);
176  }
177 
178  //-- normalize the probability distribution
179  prob->Scale( 1.0 / prob->Integral("width") );
180 
181  return prob;
182 }
183 //____________________________________________________________________________
185 {
186  Algorithm::Configure(config);
187  this->LoadConfig();
188 }
189 //____________________________________________________________________________
190 void LocalFGM::Configure(string param_set)
191 {
192  Algorithm::Configure(param_set);
193  this->LoadConfig();
194 }
195 //____________________________________________________________________________
197 {
198  this->GetParamDef("MomentumMax", fPMax, 1.0);
199  assert(fPMax > 0);
200 
201  // Load removal energy for specific nuclei from either the algorithm's
202  // configuration file or the UserPhysicsOptions file.
203  // If none is used use Wapstra's semi-empirical formula.
204  //
205 
206  for(int Z=1; Z<140; Z++) {
207  for(int A=Z; A<3*Z; A++) {
208  ostringstream key ;
209  int pdgc = pdg::IonPdgCode(A,Z);
210  key << "RFG-NucRemovalE@Pdg=" << pdgc;
211  RgKey rgkey = key.str();
212  double eb ;
213  if ( GetParam( rgkey, eb, false ) ) {
214  eb = TMath::Max(eb, 0.);
215  LOG("LocalFGM", pINFO) << "Nucleus: " << pdgc << " -> using Eb = " << eb << " GeV";
216  fNucRmvE.insert(map<int,double>::value_type(Z,eb));
217  }
218  }
219  }
220 }
221 //____________________________________________________________________________
const double kPi
Basic constants.
string AsString(void) const
Definition: Target.cxx:400
const XML_Char * target
Definition: expat.h:268
set< int >::iterator it
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
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
const char * p
Definition: xmltok.h:285
double Prob(double p, double w, const Target &t, double hitNucleonRadius) const
Definition: LocalFGM.cxx:96
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
Definition: config.py:1
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
Float_t Z
Definition: plot.C:38
const double C
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:304
double BindEnergyPerNucleon(const Target &target)
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
map< int, double > fNucRmvE
Definition: LocalFGM.h:70
double dx[NP][NC]
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void LoadConfig(void)
Definition: LocalFGM.cxx:196
const double a
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
bool GenerateNucleon(const Target &t, double hitNucleonRadius) const
Definition: LocalFGM.cxx:52
int Z(void) const
Definition: Target.h:69
#define pINFO
Definition: Messenger.h:63
double fPMax
Definition: LocalFGM.h:72
float bin[41]
Definition: plottest35.C:14
static const double kLightSpeed
Definition: Constants.h:32
TH1D * ProbDistro(const Target &t, double r) const
Definition: LocalFGM.cxx:112
static const double A
Definition: Units.h:82
int N(void) const
Definition: Target.h:70
bool HitNucIsSet(void) const
Definition: Target.cxx:300
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
virtual ~LocalFGM()
Definition: LocalFGM.cxx:47
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:81
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
exit(0)
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
void Configure(const Registry &config)
Definition: LocalFGM.cxx:184
#define pNOTICE
Definition: Messenger.h:62
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
Float_t w
Definition: plot.C:20
static const double kPlankConstant
Definition: Constants.h:33
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
static const double kPi2
Definition: Constants.h:39
Root of GENIE utility namespaces.
#define pDEBUG
Definition: Messenger.h:64