GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
SpectralFunc.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  @ May 01, 2012 - CA
14  Pick spectral function data from $GENIE/data/evgen/nucl/spectral_functions
15 */
16 //____________________________________________________________________________
17 
18 #include <TSystem.h>
19 #include <TNtupleD.h>
20 #include <TGraph2D.h>
21 
29 
30 using namespace genie;
31 using namespace genie::constants;
32 using namespace genie::controls;
33 
34 //____________________________________________________________________________
36 NuclearModelI("genie::SpectralFunc")
37 {
38  fSfFe56 = 0;
39  fSfC12 = 0;
40 }
41 //____________________________________________________________________________
43 NuclearModelI("genie::SpectralFunc", config)
44 {
45  fSfFe56 = 0;
46  fSfC12 = 0;
47 }
48 //____________________________________________________________________________
50 {
51 //if (fSfFe56) delete fSfFe56;
52 //if (fSfC12 ) delete fSfC12;
53 }
54 //____________________________________________________________________________
55 bool SpectralFunc::GenerateNucleon(const Target & target) const
56 {
57  TGraph2D * sf = this->SelectSpectralFunction(target);
58 
59  if(!sf) {
60  fCurrRemovalEnergy = 0.;
61  fCurrMomentum.SetXYZ(0.,0.,0.);
62  return false;
63  }
64 
65  double kmin = sf->GetXmin(); // momentum range
66  double kmax = sf->GetXmax();
67  double wmin = sf->GetYmin(); // removal energy range
68  double wmax = sf->GetYmax();
69  double probmax = sf->GetZmax(); // maximum probability
70  probmax *= 1.1;
71 
72  double dk = kmax - kmin;
73  double dw = wmax - wmin;
74 
75  LOG("SpectralFunc", pINFO) << "Momentum range = [" << kmin << ", " << kmax << "]";
76  LOG("SpectralFunc", pINFO) << "Rmv energy range = [" << wmin << ", " << wmax << "]";
77 
79 
80  unsigned int niter = 0;
81  while(1) {
82  if(niter > kRjMaxIterations) {
83  LOG("SpectralFunc", pWARN)
84  << "Couldn't generate a hit nucleon after " << niter << " iterations";
85  return false;
86  }
87  niter++;
88 
89  // random pair
90  double kc = kmin + dk * rnd->RndGen().Rndm();
91  double wc = wmin + dw * rnd->RndGen().Rndm();
92  LOG("SpectralFunc", pINFO) << "Trying p = " << kc << ", w = " << wc;
93 
94  // accept/reject
95  double prob = this->Prob(kc,wc, target);
96  double probg = probmax * rnd->RndGen().Rndm();
97  bool accept = (probg < prob);
98  if(!accept) continue;
99 
100  LOG("SpectralFunc", pINFO) << "|p,nucleon| = " << kc;
101  LOG("SpectralFunc", pINFO) << "|w,nucleon| = " << wc;
102 
103  // generate momentum components
104  double costheta = -1. + 2. * rnd->RndGen().Rndm();
105  double sintheta = TMath::Sqrt(1.-costheta*costheta);
106  double fi = 2 * kPi * rnd->RndGen().Rndm();
107  double cosfi = TMath::Cos(fi);
108  double sinfi = TMath::Sin(fi);
109 
110  double kx = kc*sintheta*cosfi;
111  double ky = kc*sintheta*sinfi;
112  double kz = kc*costheta;
113 
114  // set generated values
115  fCurrRemovalEnergy = wc;
116  fCurrMomentum.SetXYZ(kx,ky,kz);
117 
118  return true;
119  }
120  return false;
121 }
122 //____________________________________________________________________________
124  double p, double w, const Target & target) const
125 {
126  TGraph2D * sf = this->SelectSpectralFunction(target);
127  if(!sf) return 0;
128 
129  return sf->Interpolate(p,w);
130 }
131 //____________________________________________________________________________
132 void SpectralFunc::Configure(const Registry & config)
133 {
134  Algorithm::Configure(config);
135  this->LoadConfig();
136 }
137 //____________________________________________________________________________
138 void SpectralFunc::Configure(string param_set)
139 {
140  Algorithm::Configure(param_set);
141  this->LoadConfig();
142 }
143 //____________________________________________________________________________
145 {
146 
148 
149  LOG("SpectralFunc", pDEBUG) << "Loading Benhar et al. spectral functions";
150 
151  string data_dir =
152  string(gSystem->Getenv("GENIE")) +
153  string("/data/evgen/nucl/spectral_functions/");
154  string c12file = data_dir + "benhar-sf-12c.data";
155  string fe56file = data_dir + "benhar-sf-56fe.data";
156 
157  TNtupleD sfdata_fe56("sfdata_fe56","","k:e:prob");
158  TNtupleD sfdata_c12 ("sfdata_c12", "","k:e:prob");
159 
160  sfdata_fe56.ReadFile ( fe56file.c_str() );
161  sfdata_c12. ReadFile ( c12file.c_str () );
162 
163  LOG("SpectralFunc", pDEBUG) << "Loaded " << sfdata_fe56.GetEntries() << " Fe56 points";
164  LOG("SpectralFunc", pDEBUG) << "Loaded " << sfdata_c12.GetEntries() << " C12 points";
165 
166  if (fSfFe56) delete fSfFe56;
167  if (fSfC12 ) delete fSfC12;
168 
169  fSfFe56 = this->Convert2Graph(sfdata_fe56);
170  fSfC12 = this->Convert2Graph(sfdata_c12);
171 
172  fSfFe56->SetName("sf_fe56");
173  fSfC12 ->SetName("sf_c12");
174 }
175 //____________________________________________________________________________
176 TGraph2D * SpectralFunc::Convert2Graph(TNtupleD & sfdata) const
177 {
178  int np = sfdata.GetEntries();
179  TGraph2D * sfgraph = new TGraph2D(np);
180 
181  sfdata.Draw("k:e:prob","","GOFF");
182  assert(np==sfdata.GetSelectedRows());
183  double * k = sfdata.GetV1();
184  double * e = sfdata.GetV2();
185  double * p = sfdata.GetV3();
186 
187  for(int i=0; i<np; i++) {
188  double ki = k[i] * (units::MeV/units::GeV); // momentum
189  double ei = e[i] * (units::MeV/units::GeV); // removal energy
190  double pi = p[i] * TMath::Power(ki,2); // probabillity
191  sfgraph->SetPoint(i, ki,ei,pi);
192  }
193  return sfgraph;
194 }
195 //____________________________________________________________________________
196 TGraph2D * SpectralFunc::SelectSpectralFunction(const Target & t) const
197 {
198  TGraph2D * sf = 0;
199  int pdgc = t.Pdg();
200 
201  if (pdgc == kPdgTgtC12) sf = fSfC12;
202  else if (pdgc == kPdgTgtFe56) sf = fSfFe56;
203  else {
204  LOG("SpectralFunc", pERROR)
205  << "** The spectral function for target " << pdgc << " isn't available";
206  }
207  if(!sf) {
208  LOG("SpectralFunc", pERROR) << "** Null spectral function";
209  }
210  return sf;
211 }
212 //____________________________________________________________________________
TGraph2D * fSfFe56
Benhar&#39;s Fe56 SF.
Definition: SpectralFunc.h:62
#define pERROR
Definition: Messenger.h:59
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
TGraph2D * Convert2Graph(TNtupleD &data) const
int Pdg(void) const
Definition: Target.h:71
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:46
static constexpr double MeV
Definition: Units.h:129
double Prob(double p, double w, const Target &t) const
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
const double e
#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 GeV
Definition: Units.h:28
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
#define pINFO
Definition: Messenger.h:62
bool GenerateNucleon(const Target &t) const
#define pWARN
Definition: Messenger.h:60
virtual void LoadConfig()
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
void Configure(const Registry &config)
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:80
TGraph2D * fSfC12
Benhar&#39;s C12 SF.
Definition: SpectralFunc.h:63
const int kPdgTgtC12
Definition: PDGCodes.h:202
TGraph2D * SelectSpectralFunction(const Target &target) const
const int kPdgTgtFe56
Definition: PDGCodes.h:205
static const unsigned int kRjMaxIterations
Definition: Controls.h:26
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
#define pDEBUG
Definition: Messenger.h:63