GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SpectralFunc1d.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 
14 */
15 //____________________________________________________________________________
16 
17 #include <sstream>
18 
19 #include <TSystem.h>
20 
31 
32 using std::ostringstream;
33 
34 using namespace genie;
35 using namespace genie::constants;
36 using namespace genie::controls;
37 using namespace genie::utils;
38 
39 //____________________________________________________________________________
41 NuclearModelI("genie::SpectralFunc1d")
42 {
43 
44 }
45 //____________________________________________________________________________
47 NuclearModelI("genie::SpectralFunc1d", config)
48 {
49 
50 }
51 //____________________________________________________________________________
53 {
54  this->CleanUp();
55 }
56 //____________________________________________________________________________
57 bool SpectralFunc1d::GenerateNucleon(const Target & target) const
58 {
60  int Z = target.Z();
61 
62  map<int, double>::const_iterator dbl_it;
63  map<int, Spline*>::const_iterator spl_it;
64 
65  // Select fermi momentum from the integrated (over removal energies) s/f.
66  //
67  spl_it = fSFk.find(Z);
68  dbl_it = fMaxProb.find(Z);
69  if(spl_it == fSFk.end() || dbl_it == fMaxProb.end()) {
70  fCurrRemovalEnergy = 0.;
71  fCurrMomentum.SetXYZ(0.,0.,0.);
72  return false;
73  }
74 
75  double prob_max = dbl_it->second;
76  LOG("SpectralFunc1", pDEBUG) << "Max probability = " << prob_max;
77 
78  double p = 0;
79  unsigned int niter = 0;
80  while(1) {
81  if(niter > kRjMaxIterations) {
82  LOG("SpectralFunc1", pWARN)
83  << "Couldn't generate a hit nucleon after " << niter << " iterations";
84  return false;
85  }
86  niter++;
87 
88  if(fUseRFGMomentumCutoff) p = fPCutOff * rnd->RndGen().Rndm();
89  else p = rnd->RndGen().Rndm();
90 
91  double prob = spl_it->second->Evaluate(p);
92  double probg = prob_max * rnd->RndGen().Rndm();
93  LOG("SpectralFunc1", pDEBUG) << "Trying p = " << p << " / prob = " << prob;
94 
95  bool accept = (probg<prob);
96  if(accept) break;
97  }
98 
99  LOG("SpectralFunc1", pINFO) << "|p,nucleon| = " << p;
100 
101  double costheta = -1. + 2. * rnd->RndGen().Rndm();
102  double sintheta = TMath::Sqrt(1.-costheta*costheta);
103  double fi = 2 * kPi * rnd->RndGen().Rndm();
104  double cosfi = TMath::Cos(fi);
105  double sinfi = TMath::Sin(fi);
106 
107  double px = p*sintheta*cosfi;
108  double py = p*sintheta*sinfi;
109  double pz = p*costheta;
110 
111  fCurrMomentum.SetXYZ(px,py,pz);
112 
113  // Set removal energy
114  // Do it either in the same way as in the FG model or by using the average
115  // removal energy for the seleced pF as calculated from the s/f itself
116  //
117  if(fUseRFGRemovalE) {
118  dbl_it = fNucRmvE.find(Z);
119  if(dbl_it != fNucRmvE.end()) fCurrRemovalEnergy = dbl_it->second;
121  } else {
122  spl_it = fSFw.find(Z);
123  if(spl_it==fSFw.end()) {
124  fCurrRemovalEnergy = 0.;
125  fCurrMomentum.SetXYZ(0.,0.,0.);
126  return false;
127  } else fCurrRemovalEnergy = spl_it->second->Evaluate(p);
128  }
129 
130  return true;
131 }
132 //____________________________________________________________________________
134  double p, double /*w*/, const Target & target) const
135 {
136  if(fUseRFGMomentumCutoff) { if(p > fPCutOff) return 0; }
137 
138  int Z = target.Z();
139  map<int, Spline*>::const_iterator spl_it = fSFk.find(Z);
140 
141  if(spl_it == fSFk.end()) return 0;
142  else return TMath::Max(0.,spl_it->second->Evaluate(p));
143 }
144 //____________________________________________________________________________
146 {
147  Algorithm::Configure(config);
148  this->LoadConfig();
149 }
150 //____________________________________________________________________________
151 void SpectralFunc1d::Configure(string param_set)
152 {
153  Algorithm::Configure(param_set);
154  this->LoadConfig();
155 }
156 //____________________________________________________________________________
158 {
159 
161 
162  LOG("SpectralFunc1", pDEBUG) << "Loading coonfiguration for SpectralFunc1d";
163 
164  this->CleanUp();
165 
166  // Load spectral function data.
167  // Hopefully analytical expressions will be available soon.
168  // Currently I have spectral functions for C12 and Fe56 only.
169  //
170  string data_dir =
171  string(gSystem->Getenv("GENIE")) + string("/data/evgen/nucl/spectral_functions/");
172 
173  string c12_sf1dk_file = data_dir + "benhar-sf1dk-12c.data";
174  string fe56_sf1dk_file = data_dir + "benhar-sf1dk-56fe.data";
175  string c12_sf1dw_file = data_dir + "benhar-sf1dw-12c.data";
176  string fe56_sf1dw_file = data_dir + "benhar-sf1dw-56fe.data";
177 
178  Spline * spl = 0;
179 
180  spl = new Spline(c12_sf1dk_file);
181  fSFk.insert(map<int, Spline*>::value_type(6,spl));
182  spl = new Spline(fe56_sf1dk_file);
183  fSFk.insert(map<int, Spline*>::value_type(26,spl));
184 
185  spl = new Spline(c12_sf1dw_file);
186  fSFw.insert(map<int, Spline*>::value_type(6,spl));
187  spl = new Spline(fe56_sf1dw_file);
188  fSFw.insert(map<int, Spline*>::value_type(26,spl));
189 
190  // scan for max.
191  map<int, Spline*>::const_iterator spliter;
192  int n = 100;
193  double pmin = 0.000;
194  double dp = 0.010;
195  for(spliter = fSFk.begin(); spliter != fSFk.end(); ++spliter) {
196  double prob_max = 0;
197  int Z = spliter->first;
198  spl = spliter->second;
199  for(int i=0; i<n; i++) {
200  double p = pmin + i*dp;
201  prob_max = TMath::Max(prob_max, spl->Evaluate(p));
202  }
203  fMaxProb.insert(map<int,double>::value_type(Z,prob_max));
204  }
205 
206  // Check whether to use the same removal energies as in the FG model or
207  // to use the average removal energy for the selected fermi momentum
208  // (computed from the spectral function itself)
209  GetParam( "UseRFGRemovalE", fUseRFGRemovalE ) ;
210 
211  // Check whether to use the same momentum cutoff as in the FG model
212  GetParam("UseRFGMomentumCutoff", fUseRFGMomentumCutoff ) ;
213 
214  //Get the momentum cutoff
215  GetParam( "RFG-MomentumCutOff", fPCutOff ) ;
216 
217  // Removal energies as used in the FG model
218  // Load removal energy for specific nuclei from either the algorithm's
219  // configuration file or the UserPhysicsOptions file.
220  // If none is used use Wapstra's semi-empirical formula.
221  for(int Z=1; Z<140; Z++) {
222  for(int A=Z; A<3*Z; A++) {
223  ostringstream key;
224  int pdgc = pdg::IonPdgCode(A,Z);
225  key << "RFG-NucRemovalE@Pdg=" << pdgc;
226  RgKey rgkey = key.str();
227  double eb ;
228  if ( GetParam( rgkey, eb, false ) ) {
229  eb = TMath::Max(eb, 0.);
230  LOG("BodekRitchie", pNOTICE)
231  << "Nucleus: " << pdgc << " -> using Eb = " << eb << " GeV";
232  fNucRmvE.insert(map<int,double>::value_type(Z,eb));
233  }
234  }
235  }
236 }
237 //____________________________________________________________________________
239 {
240  map<int, Spline*>::iterator spliter;
241 
242  for(spliter = fSFk.begin(); spliter != fSFk.end(); ++spliter) {
243  Spline * spl = spliter->second;
244  if(spl) delete spl;
245  }
246  for(spliter = fSFw.begin(); spliter != fSFw.end(); ++spliter) {
247  Spline * spl = spliter->second;
248  if(spl) delete spl;
249  }
250  fSFk.clear();
251  fSFw.clear();
252  fNucRmvE.clear();
253  fMaxProb.clear();
254 }
255 //____________________________________________________________________________
256 
void Configure(const Registry &config)
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
bool GenerateNucleon(const Target &t) const
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:58
map< int, double > fNucRmvE
Removal energies as used in FG model.
map< int, Spline * > fSFk
All available spectral funcs integrated over removal energy.
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
double Evaluate(double x) const
Definition: Spline.cxx:363
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
double BindEnergyPerNucleon(const Target &target)
#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
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
#define pINFO
Definition: Messenger.h:62
#define pWARN
Definition: Messenger.h:60
virtual void LoadConfig()
double Prob(double p, double w, const Target &t) const
string RgKey
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
map< int, double > fMaxProb
Max SF(k) probability used in rejection method.
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:80
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
#define pNOTICE
Definition: Messenger.h:61
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
map< int, Spline * > fSFw
Average nucleon removal as a function of pF - computed from the spectral function.
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
#define pDEBUG
Definition: Messenger.h:63