GCylindTH1Flux.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 - July 04, 2005
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13 
14  @ Feb 22, 2011 - JD
15  Implemented dummy versions of the new GFluxI::Clear, GFluxI::Index and
16  GFluxI::GenerateWeighted methods needed for pre-generation of flux
17  interaction probabilities in GMCJDriver.
18 
19 */
20 //____________________________________________________________________________
21 
22 #include <cassert>
23 #include <algorithm>
24 
25 #include <TH1D.h>
26 #include <TF1.h>
27 #include <TVector3.h>
28 
35 
38 
39 using namespace genie;
40 using namespace genie::constants;
41 using namespace genie::flux;
42 
43 //____________________________________________________________________________
45 {
46  this->Initialize();
47 }
48 //___________________________________________________________________________
49 GCylindTH1Flux::~GCylindTH1Flux()
50 {
51  this->CleanUp();
52 }
53 //___________________________________________________________________________
55 {
56  //-- Reset previously generated neutrino code / 4-p / 4-x
57  this->ResetSelection();
58 
59  //-- Generate an energy from the 'combined' spectrum histogram
60  // and compute the momentum vector
61  double Ev = (double) fTotSpectrum->GetRandom();
62 
63  TVector3 p3(*fDirVec); // momentum along the neutrino direction
64  p3.SetMag(Ev); // with |p|=Ev
65 
66  fgP4.SetPxPyPzE(p3.Px(), p3.Py(), p3.Pz(), Ev);
67 
68  //-- Select a neutrino species from the flux fractions at the
69  // selected energy
70  fgPdgC = (*fPdgCList)[this->SelectNeutrino(Ev)];
71 
72  //-- Compute neutrino 4-x
73 
74  if(fRt <= 0) {
75  fgX4.SetXYZT(0.,0.,0.,0.);
76  }
77  else {
78  // Create a vector (vec) that points to a random position at a disk
79  // of radius Rt passing through the origin, perpendicular to the
80  // input direction.
81  TVector3 vec0(*fDirVec); // vector along the input direction
82  TVector3 vec = vec0.Orthogonal(); // orthogonal vector
83 
84  double psi = this->GeneratePhi(); // rndm angle [0,2pi]
85  double Rt = this->GenerateRt(); // rndm R [0,Rtransverse]
86 
87  vec.Rotate(psi,vec0); // rotate around original vector
88  vec.SetMag(Rt); // set new norm
89 
90  // Set the neutrino position as beam_spot + vec
91  double x = fBeamSpot->X() + vec.X();
92  double y = fBeamSpot->Y() + vec.Y();
93  double z = fBeamSpot->Z() + vec.Z();
94 
95  fgX4.SetXYZT(x,y,z,0.);
96  }
97 
98  LOG("Flux", pINFO) << "Generated neutrino pdg-code: " << fgPdgC;
99  LOG("Flux", pINFO)
100  << "Generated neutrino p4: " << utils::print::P4AsShortString(&fgP4);
101  LOG("Flux", pINFO)
102  << "Generated neutrino x4: " << utils::print::X4AsString(&fgX4);
103 
104  return true;
105 }
106 //___________________________________________________________________________
107 void GCylindTH1Flux::Clear(Option_t * opt)
108 {
109 // Dummy clear method needed to conform to GFluxI interface
110 //
111  LOG("Flux", pERROR) <<
112  "No Clear(Option_t * opt) method implemented for opt: "<< opt;
113 }
114 //___________________________________________________________________________
115 void GCylindTH1Flux::GenerateWeighted(bool gen_weighted)
116 {
117 // Dummy implementation needed to conform to GFluxI interface
118 //
119  LOG("Flux", pERROR) <<
120  "No GenerateWeighted(bool gen_weighted) method implemented for " <<
121  "gen_weighted: " << gen_weighted;
122 }
123 //___________________________________________________________________________
125 {
126  LOG("Flux", pNOTICE) << "Initializing GCylindTH1Flux driver";
127 
128  fMaxEv = 0;
129  fPdgCList = new PDGCodeList;
130  fTotSpectrum = 0;
131  fDirVec = 0;
132  fBeamSpot = 0;
133  fRt =-1;
134  fRtDep = 0;
135 
136  this->ResetSelection();
137  this->SetRtDependence("x");
138  //eg, other example: this->SetRtDependence("pow(x,2)");
139 }
140 //___________________________________________________________________________
142 {
143 // initializing running neutrino pdg-code, 4-position, 4-momentum
144  fgPdgC = 0;
145  fgP4.SetPxPyPzE (0.,0.,0.,0.);
146  fgX4.SetXYZT (0.,0.,0.,0.);
147 }
148 //___________________________________________________________________________
150 {
151  LOG("Flux", pNOTICE) << "Cleaning up...";
152 
153  if (fDirVec ) delete fDirVec;
154  if (fBeamSpot ) delete fBeamSpot;
155  if (fPdgCList ) delete fPdgCList;
156  if (fTotSpectrum) delete fTotSpectrum;
157  if (fRtDep ) delete fRtDep;
158 
159  unsigned int nspectra = fSpectrum.size();
160  for(unsigned int i = 0; i < nspectra; i++) {
161  TH1D * spectrum = fSpectrum[i];
162  delete spectrum;
163  spectrum = 0;
164  }
165 }
166 //___________________________________________________________________________
167 void GCylindTH1Flux::SetNuDirection(const TVector3 & direction)
168 {
169  if(fDirVec) delete fDirVec;
170  fDirVec = new TVector3(direction);
171 }
172 //___________________________________________________________________________
173 void GCylindTH1Flux::SetBeamSpot(const TVector3 & spot)
174 {
175  if(fBeamSpot) delete fBeamSpot;
176  fBeamSpot = new TVector3(spot);
177 }
178 //___________________________________________________________________________
180 {
181  LOG ("Flux", pNOTICE) << "Setting R[transverse] = " << Rt;
182  fRt = Rt;
183 
184  if(fRtDep) fRtDep->SetRange(0,Rt);
185 }
186 //___________________________________________________________________________
187 void GCylindTH1Flux::AddEnergySpectrum(int nu_pdgc, TH1D * spectrum)
188 {
189  LOG("Flux", pNOTICE) << "Adding flux spectrum for pdg = " << nu_pdgc;
190 
191  fPdgCList->push_back(nu_pdgc);
192 
193  bool accepted = (count(fPdgCList->begin(),fPdgCList->end(),nu_pdgc) == 1);
194  if(!accepted) {
195  LOG ("Flux", pWARN)
196  << "The pdg-code isn't recognized and the spectrum was ignored";
197  } else {
198  fSpectrum.push_back(spectrum);
199 
200  int nb = spectrum->GetNbinsX();
201  Axis_t max = spectrum->GetBinLowEdge(nb)+spectrum->GetBinWidth(nb);
202  fMaxEv = TMath::Max(fMaxEv, (double)max);
203 
204  LOG("Flux", pNOTICE)
205  << "Updating maximum energy of flux particles to: " << fMaxEv;
206 
207  this->AddAllFluxes(); // update combined flux
208  }
209 }
210 //___________________________________________________________________________
212 {
213 // Set the (functional form of) Rt dependence as string, eg "x*x+sin(x)"
214 // You do not need to set this method. The default behaviour is to generate
215 // flux neutrinos uniformly over the area of the cylinder's cross section.
216 
217  if(fRtDep) delete fRtDep;
218 
219  fRtDep = new TF1("rdep", rdep.c_str(), 0,fRt);
220 }
221 //___________________________________________________________________________
223 {
224  LOG("Flux", pNOTICE) << "Computing combined flux";
225 
226  if(fTotSpectrum) delete fTotSpectrum;
227 
228  vector<TH1D *>::const_iterator spectrum_iter;
229 
230  unsigned int inu=0;
231  for(spectrum_iter = fSpectrum.begin();
232  spectrum_iter != fSpectrum.end(); ++spectrum_iter) {
233  TH1D * spectrum = *spectrum_iter;
234 
235  if(inu==0) { fTotSpectrum = new TH1D(*spectrum); }
236  else { fTotSpectrum->Add(spectrum); }
237  inu++;
238  }
239 }
240 //___________________________________________________________________________
242 {
243  const unsigned int n = fPdgCList->size();
244  double fraction[n];
245 
246  vector<TH1D *>::const_iterator spectrum_iter;
247 
248  unsigned int inu=0;
249  for(spectrum_iter = fSpectrum.begin();
250  spectrum_iter != fSpectrum.end(); ++spectrum_iter) {
251  TH1D * spectrum = *spectrum_iter;
252  fraction[inu++] = spectrum->GetBinContent(spectrum->FindBin(Ev));
253  }
254 
255  double sum = 0;
256  for(inu = 0; inu < n; inu++) {
257  sum += fraction[inu];
258  fraction[inu] = sum;
259  LOG("Flux", pDEBUG) << "SUM-FRACTION(0->" << inu <<") = " << sum;
260  }
261 
263  double R = sum * rnd->RndFlux().Rndm();
264 
265  LOG("Flux", pDEBUG) << "R e [0,SUM] = " << R;
266 
267  for(inu = 0; inu < n; inu++) {if ( R < fraction[inu] ) return inu;}
268 
269  LOG("Flux", pERROR) << "Could not select a neutrino species";
270  assert(false);
271 
272  return -1;
273 }
274 //___________________________________________________________________________
275 double GCylindTH1Flux::GeneratePhi(void) const
276 {
278  double phi = 2.*kPi * rnd->RndFlux().Rndm(); // [0,2pi]
279  return phi;
280 }
281 //___________________________________________________________________________
282 double GCylindTH1Flux::GenerateRt(void) const
283 {
284  double Rt = fRtDep->GetRandom(); // rndm R [0,Rtransverse]
285  return Rt;
286 }
287 //___________________________________________________________________________
const double kPi
A class for generating concrete GFluxI derived classes based on the factory pattern. This code supplies a CPP macro which allows the classes to self-register and thus no modification of this class is needed in order to expand the list of classes it knows about.
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:52
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
PDGCodeList * fPdgCList
list of neutrino pdg-codes
bool GenerateNext(void)
generate the next flux neutrino (return false in err)
void SetRtDependence(string rdep)
TLorentzVector fgP4
running generated nu 4-momentum
void Clear(Option_t *opt)
reset state variables based on opt
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
TVector3 * fBeamSpot
beam spot position
void Initialize(Bool_t useTMVAStyle=kTRUE)
Definition: tmvaglob.C:184
Loaders::FluxType flux
TLorentzVector fgX4
running generated nu 4-position
vector< TH1D * > fSpectrum
flux = f(Ev), 1/neutrino species
A list of PDG codes.
Definition: PDGCodeList.h:33
TF1 * fRtDep
transverse radius dependence
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double GenerateRt(void) const
TVector3 * fDirVec
neutrino direction
int fgPdgC
running generated nu pdg-code
A generic GENIE flux driver. Generates a &#39;cylindrical&#39; neutrino beam along the input direction...
double fRt
transverse size of neutrino beam
#define R(x)
#define pINFO
Definition: Messenger.h:63
Eigen::VectorXd vec
void SetTransverseRadius(double Rt)
z
Definition: test.py:28
void SetNuDirection(const TVector3 &direction)
#define pWARN
Definition: Messenger.h:61
double GeneratePhi(void) const
TH1D * fTotSpectrum
combined flux = f(Ev)
assert(nhit_max >=nhit_nbins)
void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
void SetBeamSpot(const TVector3 &spot)
void AddEnergySpectrum(int nu_pdgc, TH1D *spectrum)
#define pNOTICE
Definition: Messenger.h:62
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:72
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:64
static const double nb
Definition: Units.h:89
Double_t sum
Definition: plot.C:31
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
#define FLUXDRIVERREG4(_nsa, _nsb, _name, _fqname)
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
double fMaxEv
maximum energy
#define pDEBUG
Definition: Messenger.h:64