FGMBodekRitchie.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 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 
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 //____________________________________________________________________________
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 //____________________________________________________________________________
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  double Z = (double) target.Z();
145  double N = (double) target.N();
146  double A = (double) target.A();
147  double KF;
148  if (A<6 || !fUseParametrization)
149  {
150  //-- look up the Fermi momentum for this Target
152  const FermiMomentumTable * kft = kftp->GetTable(fKFTable);
153  KF = kft->FindClosestKF(target_pdgc, nucleon_pdgc);
154  LOG("BodekRitchie", pNOTICE) << "KF = " << KF;
155  }
156  else
157  {
158  //-- define the Fermi momentum for this Target
159  //
161  //-- correct the Fermi momentum for the struck nucleon
162  assert(target.HitNucIsSet());
163  bool is_p = pdg::IsProton(nucleon_pdgc);
164  if(is_p) KF *= TMath::Power( 2*Z/A, 1./3.);
165  else
166  KF *= TMath::Power( 2*N/A, 1./3.);
167  LOG("BodekRitchie", pINFO) << "Corrected KF = " << KF;
168  }
169 
170 
171 
172 
173  double a = 2.0;
174  double C = 4. * kPi * TMath::Power(KF,3) / 3.;
175  double R = 1. / (1.- KF/fPMax);
176 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
177  LOG("BodekRitchie", pDEBUG) << "a = " << a;
178  LOG("BodekRitchie", pDEBUG) << "C = " << C;
179  LOG("BodekRitchie", pDEBUG) << "R = " << R;
180 #endif
181 
182  //-- create the probability distribution
183 
184  int npbins = (int) (1000*fPMax);
185  TH1D * prob = new TH1D("", "", npbins, 0, fPMax);
186  prob->SetDirectory(0);
187 
188  double dp = fPMax / (npbins-1);
189  double iC = (C>0) ? 1./C : 0.;
190  double kfa_pi_2 = TMath::Power(KF*a/kPi,2);
191 
192  for(int i = 0; i < npbins; i++) {
193  double p = i * dp;
194  double p2 = TMath::Power(p,2);
195 
196  // calculate |phi(p)|^2
197  double phi2 = 0;
198  if (p <= KF)
199  phi2 = iC * (1. - 6.*kfa_pi_2);
200  else if ( p > KF && p < fPCutOff)
201  phi2 = iC * (2*R*kfa_pi_2*TMath::Power(KF/p,4.));
202 
203  // calculate probability density : dProbability/dp
204  double dP_dp = 4*kPi * p2 * phi2;
205 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
206  LOG("BodekRitchie", pDEBUG) << "p = " << p << ", dP/dp = " << dP_dp;
207 #endif
208  prob->Fill(p, dP_dp);
209  }
210 
211  //-- normalize the probability distribution
212  prob->Scale( 1.0 / prob->Integral("width") );
213 
214  //-- store
215  fProbDistroMap.insert(
216  map<string, TH1D*>::value_type(target.AsString(),prob));
217 
218  return prob;
219 }
220 //____________________________________________________________________________
222 {
223  Algorithm::Configure(config);
224  this->LoadConfig();
225 }
226 //____________________________________________________________________________
227 void FGMBodekRitchie::Configure(string param_set)
228 {
229  Algorithm::Configure(param_set);
230  this->LoadConfig();
231 }
232 //____________________________________________________________________________
234 {
235  this->GetParam( "FermiMomentumTable", fKFTable);
236 
237 
238  // Default value 4.0 from original paper by A. Bodek and J. L. Ritchie. Phys. Rev. D 23, 1070
239  this->GetParamDef("MomentumMax", fPMax, 4.0);
240  this->GetParam("RFG-MomentumCutOff", fPCutOff);
241  this->GetParam("RFG-UseParametrization", fUseParametrization);
242 
243  assert(fPMax > 0 && fPCutOff > 0 && fPCutOff < fPMax);
244 
245  // Load removal energy for specific nuclei from either the algorithm's
246  // configuration file or the UserPhysicsOptions file.
247  // If none is used use Wapstra's semi-empirical formula.
248  const std::string keyStart = "RFG-NucRemovalE@Pdg=";
249 
250  RgIMap entries = GetConfig().GetItemMap();
251 
252  for(RgIMap::const_iterator it = entries.begin(); it != entries.end(); ++it){
253  const std::string& key = it->first;
254  int pdg = 0;
255  int Z = 0;
256  if (0 == key.compare(0, keyStart.size(), keyStart.c_str())) {
257  pdg = atoi(key.c_str() + keyStart.size());
258  Z = pdg::IonPdgCodeToZ(pdg);
259  }
260  if (0 != pdg && 0 != Z) {
261  ostringstream key_ss ;
262  key_ss << keyStart << pdg;
263  RgKey rgkey = key_ss.str();
264  double eb ;
265  GetParam( rgkey, eb ) ;
266  eb = TMath::Max(eb, 0.);
267  LOG("BodekRitchie", pINFO)
268  << "Nucleus: " << pdg << " -> using Eb = " << eb << " GeV";
269  fNucRmvE.insert(map<int,double>::value_type(Z,eb));
270  }
271  }
272 
273  LOG("BodekRitchie", pDEBUG)
274  << "Finished LoadConfig";
275 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
276  for (map<int,double>::iterator it = fNucRmvE.begin(); it != fNucRmvE.end(); ++it) {
277  LOG("BodekRitchie", pDEBUG)
278  << "Z = " << (*it).first << "; eb = " << (*it).second;
279  }
280 #endif
281 }
282 //____________________________________________________________________________
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
int A(void) const
Definition: Target.h:71
const char * p
Definition: xmltok.h:285
map< int, double > fNucRmvE
static FermiMomentumTablePool * Instance(void)
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.
Definition: config.py:1
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:254
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
Singleton class to load & serve tables of Fermi momentum constants.
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 Configure(const Registry &config)
const RgIMap & GetItemMap(void) const
Definition: Registry.h:162
const FermiMomentumTable * GetTable(string name)
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:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
#define R(x)
int Z(void) const
Definition: Target.h:69
bool GenerateNucleon(const Target &t) const
#define pINFO
Definition: Messenger.h:63
TH1D * ProbDistro(const Target &t) const
float bin[41]
Definition: plottest35.C:14
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double Prob(double mom, double w, const Target &t) const
double FermiMomentumForIsoscalarNucleonParametrization(const Target &target)
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
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)
map< string, TH1D * > fProbDistroMap
assert(nhit_max >=nhit_nbins)
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:53
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
#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
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Root of GENIE utility namespaces.
#define pDEBUG
Definition: Messenger.h:64
map< RgKey, RegistryItemI * > RgIMap
Definition: Registry.h:46