SpectralFunc1d.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 
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 //____________________________________________________________________________
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  LOG("SpectralFunc1", pDEBUG) << "Loading coonfiguration for SpectralFunc1d";
160 
161  this->CleanUp();
162 
163  // Load spectral function data.
164  // Hopefully analytical expressions will be available soon.
165  // Currently I have spectral functions for C12 and Fe56 only.
166  //
167  string data_dir =
168  string(gSystem->Getenv("GENIE")) + string("/data/evgen/nucl/spectral_functions/");
169 
170  string c12_sf1dk_file = data_dir + "benhar-sf1dk-12c.data";
171  string fe56_sf1dk_file = data_dir + "benhar-sf1dk-56fe.data";
172  string c12_sf1dw_file = data_dir + "benhar-sf1dw-12c.data";
173  string fe56_sf1dw_file = data_dir + "benhar-sf1dw-56fe.data";
174 
175  Spline * spl = 0;
176 
177  spl = new Spline(c12_sf1dk_file);
178  fSFk.insert(map<int, Spline*>::value_type(6,spl));
179  spl = new Spline(fe56_sf1dk_file);
180  fSFk.insert(map<int, Spline*>::value_type(26,spl));
181 
182  spl = new Spline(c12_sf1dw_file);
183  fSFw.insert(map<int, Spline*>::value_type(6,spl));
184  spl = new Spline(fe56_sf1dw_file);
185  fSFw.insert(map<int, Spline*>::value_type(26,spl));
186 
187  // scan for max.
188  map<int, Spline*>::const_iterator spliter;
189  int n = 100;
190  double pmin = 0.000;
191  double dp = 0.010;
192  for(spliter = fSFk.begin(); spliter != fSFk.end(); ++spliter) {
193  double prob_max = 0;
194  int Z = spliter->first;
195  spl = spliter->second;
196  for(int i=0; i<n; i++) {
197  double p = pmin + i*dp;
198  prob_max = TMath::Max(prob_max, spl->Evaluate(p));
199  }
200  fMaxProb.insert(map<int,double>::value_type(Z,prob_max));
201  }
202 
203  // Check whether to use the same removal energies as in the FG model or
204  // to use the average removal energy for the selected fermi momentum
205  // (computed from the spectral function itself)
206  GetParam( "UseRFGRemovalE", fUseRFGRemovalE ) ;
207 
208  // Check whether to use the same momentum cutoff as in the FG model
209  GetParam("UseRFGMomentumCutoff", fUseRFGMomentumCutoff ) ;
210 
211  //Get the momentum cutoff
212  GetParam( "RFG-MomentumCutOff", fPCutOff ) ;
213 
214  // Removal energies as used in the FG model
215  // Load removal energy for specific nuclei from either the algorithm's
216  // configuration file or the UserPhysicsOptions file.
217  // If none is used use Wapstra's semi-empirical formula.
218  for(int Z=1; Z<140; Z++) {
219  for(int A=Z; A<3*Z; A++) {
220  ostringstream key;
221  int pdgc = pdg::IonPdgCode(A,Z);
222  key << "RFG-NucRemovalE@Pdg=" << pdgc;
223  RgKey rgkey = key.str();
224  double eb ;
225  if ( GetParam( rgkey, eb, false ) ) {
226  eb = TMath::Max(eb, 0.);
227  LOG("BodekRitchie", pNOTICE)
228  << "Nucleus: " << pdgc << " -> using Eb = " << eb << " GeV";
229  fNucRmvE.insert(map<int,double>::value_type(Z,eb));
230  }
231  }
232  }
233 }
234 //____________________________________________________________________________
236 {
237  map<int, Spline*>::iterator spliter;
238 
239  for(spliter = fSFk.begin(); spliter != fSFk.end(); ++spliter) {
240  Spline * spl = spliter->second;
241  if(spl) delete spl;
242  }
243  for(spliter = fSFw.begin(); spliter != fSFw.end(); ++spliter) {
244  Spline * spl = spliter->second;
245  if(spl) delete spl;
246  }
247  fSFk.clear();
248  fSFw.clear();
249  fNucRmvE.clear();
250  fMaxProb.clear();
251 }
252 //____________________________________________________________________________
253 
const double kPi
void Configure(const Registry &config)
Basic constants.
const XML_Char * target
Definition: expat.h:268
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
bool GenerateNucleon(const Target &t) const
const char * p
Definition: xmltok.h:285
A numeric analysis tool class for interpolating 1-D functions.
Definition: Spline.h:47
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:42
Definition: config.py:1
double Evaluate(double x) const
Definition: Spline.cxx:362
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
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:97
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
int Z(void) const
Definition: Target.h:69
#define pINFO
Definition: Messenger.h:63
Misc GENIE control constants.
#define pWARN
Definition: Messenger.h:61
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
static const double A
Definition: Units.h:82
double Prob(double p, double w, const Target &t) const
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
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:81
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:69
#define pNOTICE
Definition: Messenger.h:62
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
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...
Root of GENIE utility namespaces.
#define pDEBUG
Definition: Messenger.h:64