GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
FGMBodekRitchie.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 the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Feb 07, 2008 - CA
14  Call SetDirectory(0) at the temp momentum distribution histogram to stop
15  it from being automatically written out at the event file.
16  @ Jun 18, 2008 - CA
17  Deallocate the momentum distribution histograms map at dtor
18 */
19 //____________________________________________________________________________
20 
21 #include <sstream>
22 #include <cstdlib>
23 #include <TMath.h>
24 
26 #include "Framework/Conventions/GBuild.h"
36 
37 using std::ostringstream;
38 using namespace genie;
39 using namespace genie::constants;
40 using namespace genie::utils;
41 
42 //____________________________________________________________________________
44 NuclearModelI("genie::FGMBodekRitchie")
45 {
46 
47 }
48 //____________________________________________________________________________
50 NuclearModelI("genie::FGMBodekRitchie", config)
51 {
52 
53 }
54 //____________________________________________________________________________
56 {
57  map<string, TH1D*>::iterator iter = fProbDistroMap.begin();
58  for( ; iter != fProbDistroMap.end(); ++iter) {
59  TH1D * hst = iter->second;
60  if(hst) {
61  delete hst;
62  hst=0;
63  }
64  }
65  fProbDistroMap.clear();
66 }
67 //____________________________________________________________________________
68 bool FGMBodekRitchie::GenerateNucleon(const Target & target) const
69 {
70  assert(target.HitNucIsSet());
71 
73  fCurrMomentum.SetXYZ(0,0,0);
74 
75  //-- set fermi momentum vector
76  //
77  TH1D * prob = this->ProbDistro(target);
78  if ( ! prob ) {
79  LOG("BodekRitchie", pNOTICE)
80  << "Null nucleon momentum probability distribution";
81  exit(1);
82  }
83  double p = prob->GetRandom();
84  LOG("BodekRitchie", pINFO) << "|p,nucleon| = " << p;
85 
87 
88  double costheta = -1. + 2. * rnd->RndGen().Rndm();
89  double sintheta = TMath::Sqrt(1.-costheta*costheta);
90  double fi = 2 * kPi * rnd->RndGen().Rndm();
91  double cosfi = TMath::Cos(fi);
92  double sinfi = TMath::Sin(fi);
93 
94  double px = p*sintheta*cosfi;
95  double py = p*sintheta*sinfi;
96  double pz = p*costheta;
97 
98  fCurrMomentum.SetXYZ(px,py,pz);
99 
100  //-- set removal energy
101  //
102  if ( target.A() < 6 || ! fUseParametrization )
103  {
104  int Z = target.Z();
105  map<int,double>::const_iterator it = fNucRmvE.find(Z);
106  if(it != fNucRmvE.end()) fCurrRemovalEnergy = it->second;
107  else fCurrRemovalEnergy = nuclear::BindEnergyPerNucleon(target);
108  }
109  else {
110  fCurrRemovalEnergy = nuclear::BindEnergyPerNucleonParametrization(target);
111  }
112 
113  return true;
114 }
115 //____________________________________________________________________________
116 double FGMBodekRitchie::Prob(double mom, double w, const Target & target) const
117 {
118  if(w<0) {
119  TH1D * prob_distr = this->ProbDistro(target);
120  int bin = prob_distr->FindBin(mom);
121  double y = prob_distr->GetBinContent(bin);
122  double dx = prob_distr->GetBinWidth(bin);
123  double prob = y*dx;
124  return prob;
125  }
126  return 1;
127 }
128 //____________________________________________________________________________
129 TH1D * FGMBodekRitchie::ProbDistro(const Target & target) const
130 {
131  //-- return stored /if already computed/
132  map<string, TH1D*>::iterator it = fProbDistroMap.find(target.AsString());
133  if(it != fProbDistroMap.end()) return it->second;
134 
135  LOG("BodekRitchie", pNOTICE)
136  << "Computing P = f(p_nucleon) for: " << target.AsString();
137  LOG("BodekRitchie", pNOTICE)
138  << "P(cut-off) = " << fPCutOff << ", P(max) = " << fPMax;
139 
140  //-- get information for the nuclear target
141  //int target_pdgc = pdg::IonPdgCode(target.A(), target.Z());
142  int nucleon_pdgc = target.HitNucPdg();
143  assert( pdg::IsProton(nucleon_pdgc) || pdg::IsNeutron(nucleon_pdgc) );
144 
145  double KF = FermiMomentum( target, nucleon_pdgc ) ;
146 
147  double a = 2.0;
148  double C = 4. * kPi * TMath::Power(KF,3) / 3.;
149  double R = 1. / (1.- KF/fPMax);
150 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
151  LOG("BodekRitchie", pDEBUG) << "a = " << a;
152  LOG("BodekRitchie", pDEBUG) << "C = " << C;
153  LOG("BodekRitchie", pDEBUG) << "R = " << R;
154 #endif
155 
156  //-- create the probability distribution
157 
158  int npbins = (int) (1000*fPMax);
159  TH1D * prob = new TH1D("", "", npbins, 0, fPMax);
160  prob->SetDirectory(0);
161 
162  double dp = fPMax / (npbins-1);
163  double iC = (C>0) ? 1./C : 0.;
164  double kfa_pi_2 = TMath::Power(KF*a/kPi,2);
165 
166  for(int i = 0; i < npbins; i++) {
167  double p = i * dp;
168  double p2 = TMath::Power(p,2);
169 
170  // calculate |phi(p)|^2
171  double phi2 = 0;
172  if (p <= KF)
173  phi2 = iC * (1. - 6.*kfa_pi_2);
174  else if ( p > KF && p < fPCutOff)
175  phi2 = iC * (2*R*kfa_pi_2*TMath::Power(KF/p,4.));
176 
177  // calculate probability density : dProbability/dp
178  double dP_dp = 4*kPi * p2 * phi2;
179 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
180  LOG("BodekRitchie", pDEBUG) << "p = " << p << ", dP/dp = " << dP_dp;
181 #endif
182  prob->Fill(p, dP_dp);
183  }
184 
185  //-- normalize the probability distribution
186  prob->Scale( 1.0 / prob->Integral("width") );
187 
188  //-- store
189  fProbDistroMap.insert(
190  map<string, TH1D*>::value_type(target.AsString(),prob));
191 
192  return prob;
193 }
194 //____________________________________________________________________________
195 double FGMBodekRitchie::FermiMomentum( const Target & t, int nucleon_pdg ) const {
196 
197  assert( pdg::IsProton(nucleon_pdg) || pdg::IsNeutron(nucleon_pdg) );
198 
199  if ( t.A()<6 || !fUseParametrization ) {
200  //-- look up the Fermi momentum for this Target from tables
201  return NuclearModelI::FermiMomentum( t, nucleon_pdg ) ;
202  }
203  else
204  {
205  //-- define the Fermi momentum for this Target
207 
208  //-- correct the Fermi momentum for the struck nucleon
209  assert(t.HitNucIsSet());
210 
211  double A = (double) t.A() ;
212  if( pdg::IsProton(nucleon_pdg) ) {
213  double Z = (double) t.Z() ;
214  kF *= TMath::Power( 2*Z/A, 1./3.);
215  }
216  else {
217  double N = (double) t.N() ;
218  kF *= TMath::Power( 2*N/A, 1./3.);
219  }
220 
221  LOG("BodekRitchie", pINFO) << "Corrected KF = " << kF;
222 
223  return kF ;
224 
225  }
226 
227 }
228 //____________________________________________________________________________
230 {
231  Algorithm::Configure(config);
232  this->LoadConfig();
233 }
234 //____________________________________________________________________________
235 void FGMBodekRitchie::Configure(string param_set)
236 {
237  Algorithm::Configure(param_set);
238  this->LoadConfig();
239 }
240 //____________________________________________________________________________
242 {
243 
245 
246  // Default value 4.0 from original paper by A. Bodek and J. L. Ritchie. Phys. Rev. D 23, 1070
247  this->GetParamDef("MomentumMax", fPMax, 4.0);
248  this->GetParam("RFG-MomentumCutOff", fPCutOff);
249  this->GetParam("RFG-UseParametrization", fUseParametrization);
250 
251  assert(fPMax > 0 && fPCutOff > 0 && fPCutOff < fPMax);
252 
253  // Load removal energy for specific nuclei from either the algorithm's
254  // configuration file or the UserPhysicsOptions file.
255  // If none is used use Wapstra's semi-empirical formula.
256  const std::string keyStart = "RFG-NucRemovalE@Pdg=";
257 
258  RgIMap entries = GetConfig().GetItemMap();
259 
260  for(RgIMap::const_iterator it = entries.begin(); it != entries.end(); ++it){
261  const std::string& key = it->first;
262  int pdg = 0;
263  int Z = 0;
264  if (0 == key.compare(0, keyStart.size(), keyStart.c_str())) {
265  pdg = atoi(key.c_str() + keyStart.size());
266  Z = pdg::IonPdgCodeToZ(pdg);
267  }
268  if (0 != pdg && 0 != Z) {
269  ostringstream key_ss ;
270  key_ss << keyStart << pdg;
271  RgKey rgkey = key_ss.str();
272  double eb ;
273  GetParam( rgkey, eb ) ;
274  eb = TMath::Max(eb, 0.);
275  LOG("BodekRitchie", pINFO)
276  << "Nucleus: " << pdg << " -> using Eb = " << eb << " GeV";
277  fNucRmvE.insert(map<int,double>::value_type(Z,eb));
278  }
279  }
280 
281  LOG("BodekRitchie", pDEBUG)
282  << "Finished LoadConfig";
283 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
284  for (map<int,double>::iterator it = fNucRmvE.begin(); it != fNucRmvE.end(); ++it) {
285  LOG("BodekRitchie", pDEBUG)
286  << "Z = " << (*it).first << "; eb = " << (*it).second;
287  }
288 #endif
289 }
290 //____________________________________________________________________________
string AsString(void) const
Definition: Target.cxx:383
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
int HitNucPdg(void) const
Definition: Target.cxx:304
int A(void) const
Definition: Target.h:70
map< int, double > fNucRmvE
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:246
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
virtual double FermiMomentum(const Target &, int nucleon_pdg) const
double BindEnergyPerNucleon(const Target &target)
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void Configure(const Registry &config)
const RgIMap & GetItemMap(void) const
Definition: Registry.h:161
static constexpr double A
Definition: Units.h:74
const double a
double BindEnergyPerNucleonParametrization(const Target &target)
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
int Z(void) const
Definition: Target.h:68
bool GenerateNucleon(const Target &t) const
#define pINFO
Definition: Messenger.h:62
TH1D * ProbDistro(const Target &t) const
virtual void LoadConfig()
double Prob(double mom, double w, const Target &t) const
double FermiMomentumForIsoscalarNucleonParametrization(const Target &target)
int N(void) const
Definition: Target.h:69
bool HitNucIsSet(void) const
Definition: Target.cxx:283
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:80
map< string, TH1D * > fProbDistroMap
virtual double FermiMomentum(const Target &t, int nucleon_pdg) const
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:55
#define pNOTICE
Definition: Messenger.h:61
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
#define pDEBUG
Definition: Messenger.h:63
map< RgKey, RegistryItemI * > RgIMap
Definition: Registry.h:45