SpectralFunc.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  @ 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 //____________________________________________________________________________
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 //____________________________________________________________________________
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  LOG("SpectralFunc", pDEBUG) << "Loading Benhar et al. spectral functions";
147 
148  string data_dir =
149  string(gSystem->Getenv("GENIE")) +
150  string("/data/evgen/nucl/spectral_functions/");
151  string c12file = data_dir + "benhar-sf-12c.data";
152  string fe56file = data_dir + "benhar-sf-56fe.data";
153 
154  TNtupleD sfdata_fe56("sfdata_fe56","","k:e:prob");
155  TNtupleD sfdata_c12 ("sfdata_c12", "","k:e:prob");
156 
157  sfdata_fe56.ReadFile ( fe56file.c_str() );
158  sfdata_c12. ReadFile ( c12file.c_str () );
159 
160  LOG("SpectralFunc", pDEBUG) << "Loaded " << sfdata_fe56.GetEntries() << " Fe56 points";
161  LOG("SpectralFunc", pDEBUG) << "Loaded " << sfdata_c12.GetEntries() << " C12 points";
162 
163  if (fSfFe56) delete fSfFe56;
164  if (fSfC12 ) delete fSfC12;
165 
166  fSfFe56 = this->Convert2Graph(sfdata_fe56);
167  fSfC12 = this->Convert2Graph(sfdata_c12);
168 
169  fSfFe56->SetName("sf_fe56");
170  fSfC12 ->SetName("sf_c12");
171 }
172 //____________________________________________________________________________
173 TGraph2D * SpectralFunc::Convert2Graph(TNtupleD & sfdata) const
174 {
175  int np = sfdata.GetEntries();
176  TGraph2D * sfgraph = new TGraph2D(np);
177 
178  sfdata.Draw("k:e:prob","","GOFF");
179  assert(np==sfdata.GetSelectedRows());
180  double * k = sfdata.GetV1();
181  double * e = sfdata.GetV2();
182  double * p = sfdata.GetV3();
183 
184  for(int i=0; i<np; i++) {
185  double ki = k[i] * (units::MeV/units::GeV); // momentum
186  double ei = e[i] * (units::MeV/units::GeV); // removal energy
187  double pi = p[i] * TMath::Power(ki,2); // probabillity
188  sfgraph->SetPoint(i, ki,ei,pi);
189  }
190  return sfgraph;
191 }
192 //____________________________________________________________________________
194 {
195  TGraph2D * sf = 0;
196  int pdgc = t.Pdg();
197 
198  if (pdgc == kPdgTgtC12) sf = fSfC12;
199  else if (pdgc == kPdgTgtFe56) sf = fSfFe56;
200  else {
201  LOG("SpectralFunc", pERROR)
202  << "** The spectral function for target " << pdgc << " isn't available";
203  }
204  if(!sf) {
205  LOG("SpectralFunc", pERROR) << "** Null spectral function";
206  }
207  return sf;
208 }
209 //____________________________________________________________________________
const double kPi
TGraph2D * fSfFe56
Benhar&#39;s Fe56 SF.
Definition: SpectralFunc.h:59
Basic constants.
const XML_Char * target
Definition: expat.h:268
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
TGraph2D * Convert2Graph(TNtupleD &data) const
const char * p
Definition: xmltok.h:285
static const double MeV
Definition: Units.h:130
int Pdg(void) const
Definition: Target.h:72
Pure abstract base class. Defines the NuclearModelI interface to be implemented by any physics model ...
Definition: NuclearModelI.h:42
Definition: config.py:1
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:30
#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
#define pINFO
Definition: Messenger.h:63
bool GenerateNucleon(const Target &t) const
Misc GENIE control constants.
#define pWARN
Definition: Messenger.h:61
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
void Configure(const Registry &config)
cosmicTree ReadFile("CosmicList.txt","isSA/I:isA17/I:HadEfrac/F:CCE/F:CCESA/F")
TRandom3 & RndGen(void) const
rnd number generator for generic usage
Definition: RandomGen.h:81
TGraph2D * fSfC12
Benhar&#39;s C12 SF.
Definition: SpectralFunc.h:60
assert(nhit_max >=nhit_nbins)
const int kPdgTgtC12
Definition: PDGCodes.h:179
TGraph2D * SelectSpectralFunction(const Target &target) const
const int kPdgTgtFe56
Definition: PDGCodes.h:182
static const unsigned int kRjMaxIterations
Definition: Controls.h:27
Float_t e
Definition: plot.C:35
static const double GeV
Definition: Units.h:29
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...
#define pDEBUG
Definition: Messenger.h:64
enum BeamMode string