GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCylindTH1Flux.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  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <cassert>
12 #include <algorithm>
13 
14 #include <TH1D.h>
15 #include <TF1.h>
16 #include <TVector3.h>
17 
26 
29 
30 using namespace genie;
31 using namespace genie::constants;
32 using namespace genie::flux;
33 
34 //____________________________________________________________________________
36 {
37  this->Initialize();
38 }
39 //___________________________________________________________________________
40 GCylindTH1Flux::~GCylindTH1Flux()
41 {
42  this->CleanUp();
43 }
44 //___________________________________________________________________________
46 {
47  //-- Reset previously generated neutrino code / 4-p / 4-x
48  this->ResetSelection();
49 
50  //-- Generate an energy from the 'combined' spectrum histogram
51  // and compute the momentum vector
52  double Ev = (double) fTotSpectrum->GetRandom();
53 
54  TVector3 p3(*fDirVec); // momentum along the neutrino direction
55  p3.SetMag(Ev); // with |p|=Ev
56  // EDIT: Check if we're running with DM beam
58  double Md = PDGLibrary::Instance()->Find(kPdgDarkMatter)->Mass();
59  p3.SetMag(TMath::Sqrt(Ev*Ev - Md*Md));
60  }
61 
62  fgP4.SetPxPyPzE(p3.Px(), p3.Py(), p3.Pz(), Ev);
63 
64  //-- Select a neutrino species from the flux fractions at the
65  // selected energy
66  fgPdgC = (*fPdgCList)[this->SelectNeutrino(Ev)];
67 
68  //-- Compute neutrino 4-x
69 
70  if(fRt <= 0) {
71  fgX4.SetXYZT(0.,0.,0.,0.);
72  }
73  else {
74  // Create a vector (vec) that points to a random position at a disk
75  // of radius Rt passing through the origin, perpendicular to the
76  // input direction.
77  TVector3 vec0(*fDirVec); // vector along the input direction
78  TVector3 vec = vec0.Orthogonal(); // orthogonal vector
79 
80  double psi = this->GeneratePhi(); // rndm angle [0,2pi]
81  double Rt = this->GenerateRt(); // rndm R [0,Rtransverse]
82 
83  vec.Rotate(psi,vec0); // rotate around original vector
84  vec.SetMag(Rt); // set new norm
85 
86  // Set the neutrino position as beam_spot + vec
87  double x = fBeamSpot->X() + vec.X();
88  double y = fBeamSpot->Y() + vec.Y();
89  double z = fBeamSpot->Z() + vec.Z();
90 
91  fgX4.SetXYZT(x,y,z,0.);
92  }
93 
94  LOG("Flux", pINFO) << "Generated neutrino pdg-code: " << fgPdgC;
95  LOG("Flux", pINFO)
96  << "Generated neutrino p4: " << utils::print::P4AsShortString(&fgP4);
97  LOG("Flux", pINFO)
98  << "Generated neutrino x4: " << utils::print::X4AsString(&fgX4);
99 
100  return true;
101 }
102 //___________________________________________________________________________
103 void GCylindTH1Flux::Clear(Option_t * opt)
104 {
105 // Dummy clear method needed to conform to GFluxI interface
106 //
107  LOG("Flux", pERROR) <<
108  "No Clear(Option_t * opt) method implemented for opt: "<< opt;
109 }
110 //___________________________________________________________________________
111 void GCylindTH1Flux::GenerateWeighted(bool gen_weighted)
112 {
113 // Dummy implementation needed to conform to GFluxI interface
114 //
115  LOG("Flux", pERROR) <<
116  "No GenerateWeighted(bool gen_weighted) method implemented for " <<
117  "gen_weighted: " << gen_weighted;
118 }
119 //___________________________________________________________________________
121 {
122  LOG("Flux", pNOTICE) << "Initializing GCylindTH1Flux driver";
123 
124  fMaxEv = 0;
125  fPdgCList = new PDGCodeList;
126  fTotSpectrum = 0;
127  fDirVec = 0;
128  fBeamSpot = 0;
129  fRt =-1;
130  fRtDep = 0;
131 
132  this->ResetSelection();
133  this->SetRtDependence("x");
134  //eg, other example: this->SetRtDependence("pow(x,2)");
135 }
136 //___________________________________________________________________________
138 {
139 // initializing running neutrino pdg-code, 4-position, 4-momentum
140  fgPdgC = 0;
141  fgP4.SetPxPyPzE (0.,0.,0.,0.);
142  fgX4.SetXYZT (0.,0.,0.,0.);
143 }
144 //___________________________________________________________________________
146 {
147  LOG("Flux", pNOTICE) << "Cleaning up...";
148 
149  if (fDirVec ) delete fDirVec;
150  if (fBeamSpot ) delete fBeamSpot;
151  if (fPdgCList ) delete fPdgCList;
152  if (fTotSpectrum) delete fTotSpectrum;
153  if (fRtDep ) delete fRtDep;
154 
155  unsigned int nspectra = fSpectrum.size();
156  for(unsigned int i = 0; i < nspectra; i++) {
157  TH1D * spectrum = fSpectrum[i];
158  delete spectrum;
159  spectrum = 0;
160  }
161 }
162 //___________________________________________________________________________
163 void GCylindTH1Flux::SetNuDirection(const TVector3 & direction)
164 {
165  if(fDirVec) delete fDirVec;
166  fDirVec = new TVector3(direction);
167 }
168 //___________________________________________________________________________
169 void GCylindTH1Flux::SetBeamSpot(const TVector3 & spot)
170 {
171  if(fBeamSpot) delete fBeamSpot;
172  fBeamSpot = new TVector3(spot);
173 }
174 //___________________________________________________________________________
176 {
177  LOG ("Flux", pNOTICE) << "Setting R[transverse] = " << Rt;
178  fRt = Rt;
179 
180  if(fRtDep) fRtDep->SetRange(0,Rt);
181 }
182 //___________________________________________________________________________
183 void GCylindTH1Flux::AddEnergySpectrum(int nu_pdgc, TH1D * spectrum)
184 {
185  LOG("Flux", pNOTICE) << "Adding flux spectrum for pdg = " << nu_pdgc;
186 
187  fPdgCList->push_back(nu_pdgc);
188 
189  bool accepted = (count(fPdgCList->begin(),fPdgCList->end(),nu_pdgc) == 1);
190  if(!accepted) {
191  LOG ("Flux", pWARN)
192  << "The pdg-code isn't recognized and the spectrum was ignored";
193  } else {
194  fSpectrum.push_back(spectrum);
195 
196  int nb = spectrum->GetNbinsX();
197  Axis_t max = spectrum->GetBinLowEdge(nb)+spectrum->GetBinWidth(nb);
198  fMaxEv = TMath::Max(fMaxEv, (double)max);
199 
200  LOG("Flux", pNOTICE)
201  << "Updating maximum energy of flux particles to: " << fMaxEv;
202 
203  this->AddAllFluxes(); // update combined flux
204  }
205 }
206 //___________________________________________________________________________
208 {
209 // Set the (functional form of) Rt dependence as string, eg "x*x+sin(x)"
210 // You do not need to set this method. The default behaviour is to generate
211 // flux neutrinos uniformly over the area of the cylinder's cross section.
212 
213  if(fRtDep) delete fRtDep;
214 
215  fRtDep = new TF1("rdep", rdep.c_str(), 0,fRt);
216 }
217 //___________________________________________________________________________
219 {
220  LOG("Flux", pNOTICE) << "Computing combined flux";
221 
222  if(fTotSpectrum) delete fTotSpectrum;
223 
224  vector<TH1D *>::const_iterator spectrum_iter;
225 
226  unsigned int inu=0;
227  for(spectrum_iter = fSpectrum.begin();
228  spectrum_iter != fSpectrum.end(); ++spectrum_iter) {
229  TH1D * spectrum = *spectrum_iter;
230 
231  if(inu==0) { fTotSpectrum = new TH1D(*spectrum); }
232  else { fTotSpectrum->Add(spectrum); }
233  inu++;
234  }
235 }
236 //___________________________________________________________________________
238 {
239  const unsigned int n = fPdgCList->size();
240  double fraction[n];
241 
242  vector<TH1D *>::const_iterator spectrum_iter;
243 
244  unsigned int inu=0;
245  for(spectrum_iter = fSpectrum.begin();
246  spectrum_iter != fSpectrum.end(); ++spectrum_iter) {
247  TH1D * spectrum = *spectrum_iter;
248  fraction[inu++] = spectrum->GetBinContent(spectrum->FindBin(Ev));
249  }
250 
251  double sum = 0;
252  for(inu = 0; inu < n; inu++) {
253  sum += fraction[inu];
254  fraction[inu] = sum;
255  LOG("Flux", pDEBUG) << "SUM-FRACTION(0->" << inu <<") = " << sum;
256  }
257 
258  RandomGen * rnd = RandomGen::Instance();
259  double R = sum * rnd->RndFlux().Rndm();
260 
261  LOG("Flux", pDEBUG) << "R e [0,SUM] = " << R;
262 
263  for(inu = 0; inu < n; inu++) {if ( R < fraction[inu] ) return inu;}
264 
265  LOG("Flux", pERROR) << "Could not select a neutrino species";
266  assert(false);
267 
268  return -1;
269 }
270 //___________________________________________________________________________
271 double GCylindTH1Flux::GeneratePhi(void) const
272 {
273  RandomGen * rnd = RandomGen::Instance();
274  double phi = 2.*kPi * rnd->RndFlux().Rndm(); // [0,2pi]
275  return phi;
276 }
277 //___________________________________________________________________________
278 double GCylindTH1Flux::GenerateRt(void) const
279 {
280  double Rt = fRtDep->GetRandom(); // rndm R [0,Rtransverse]
281  return Rt;
282 }
283 //___________________________________________________________________________
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:45
static constexpr double nb
Definition: Units.h:81
#define pERROR
Definition: Messenger.h:59
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
PDGCodeList * fPdgCList
list of neutrino pdg-codes
bool GenerateNext(void)
generate the next flux neutrino (return false in err)
const int kPdgDarkMatter
Definition: PDGCodes.h:218
void SetRtDependence(string rdep)
TLorentzVector fgP4
running generated nu 4-momentum
bool ExistsInPDGCodeList(int pdg_code) const
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:29
TVector3 * fBeamSpot
beam spot position
TLorentzVector fgX4
running generated nu 4-position
vector< TH1D * > fSpectrum
flux = f(Ev), 1/neutrino species
A list of PDG codes.
Definition: PDGCodeList.h:32
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:96
const int kPdgAntiDarkMatter
Definition: PDGCodes.h:219
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
void Initialize(void)
#define pINFO
Definition: Messenger.h:62
void SetTransverseRadius(double Rt)
void SetNuDirection(const TVector3 &direction)
#define pWARN
Definition: Messenger.h:60
double GeneratePhi(void) const
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
TH1D * fTotSpectrum
combined flux = f(Ev)
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)
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
#define pNOTICE
Definition: Messenger.h:61
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:71
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:57
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
#define FLUXDRIVERREG4(_nsa, _nsb, _name, _fqname)
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
double fMaxEv
maximum energy
#define pDEBUG
Definition: Messenger.h:63