GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GAtmoFlux.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 <iostream>
13 #include <fstream>
14 #include <cmath>
15 
16 #include <TH3D.h>
17 #include <TMath.h>
18 
20 #include "Tools/Flux/GAtmoFlux.h"
28 
29 using namespace genie;
30 using namespace genie::flux;
31 using namespace genie::constants;
32 
33 //____________________________________________________________________________
34 GAtmoFlux::GAtmoFlux()
35 {
36  fInitialized = 0;
37 }
38 //___________________________________________________________________________
39 GAtmoFlux::~GAtmoFlux()
40 {
41  this->CleanUp();
42 }
43 //___________________________________________________________________________
44 double GAtmoFlux::MaxEnergy(void)
45 {
46  return TMath::Min(fMaxEv, fMaxEvCut);
47 }
48 //___________________________________________________________________________
49 bool GAtmoFlux::GenerateNext(void)
50 {
51  while(1) {
52  // Attempt to generate next flux neutrino
53  bool nextok = this->GenerateNext_1try();
54  if(!nextok) continue;
55 
56  // Check generated neutrino energy against max energy.
57  // We may have to reject the current neutrino if a user-defined max
58  // energy cut restricts the available range of energies.
59  const TLorentzVector & p4 = this->Momentum();
60  double E = p4.Energy();
61  double Emin = this->MinEnergy();
62  double Emax = this->MaxEnergy();
63  double wght = this->Weight();
64 
65  bool accept = (E<=Emax && E>=Emin && wght>0);
66  if(accept) return true;
67  }
68  return false;
69 }
70 //___________________________________________________________________________
71 bool GAtmoFlux::GenerateNext_1try(void)
72 {
73  // Must have run intitialization
74  assert(fInitialized);
75 
76  // Reset previously generated neutrino code / 4-p / 4-x
77  this->ResetSelection();
78 
79  // Get a RandomGen instance
81 
82  // Generate (Ev, costheta, phi)
83  double Ev = 0.;
84  double costheta = 0.;
85  double phi = 0;
86  double weight = 0;
87  int nu_pdg = 0;
88 
89  if(fGenWeighted) {
90 
91  //
92  // generate weighted flux
93  //
94 
95  // generate events according to a power law spectrum,
96  // then weight events by flux and inverse power law
97  // (note: cannot use index alpha=1)
98  double alpha = fSpectralIndex;
99 
100  double emin = TMath::Power(fEnergyBins[0],1.0-alpha);
101  double emax = TMath::Power(fEnergyBins[fNumEnergyBins],1.0-alpha);
102  Ev = TMath::Power(emin+(emax-emin)*rnd->RndFlux().Rndm(),1.0/(1.0-alpha));
103  costheta = -1+2*rnd->RndFlux().Rndm();
104  phi = 2.*kPi* rnd->RndFlux().Rndm();
105 
106  unsigned int nnu = fPdgCList->size();
107  unsigned int inu = rnd->RndFlux().Integer(nnu);
108  nu_pdg = (*fPdgCList)[inu];
109 
110  if(Ev < fEnergyBins[0]) {
111  LOG("Flux", pINFO) << "E < Emin";
112  return false;
113  }
114  double flux = this->GetFlux(nu_pdg, Ev, costheta, phi);
115  if(flux<=0) {
116  LOG("Flux", pINFO) << "Flux <= 0";
117  return false;
118  }
119  weight = flux*TMath::Power(Ev,alpha);
120  }
121  else {
122 
123  //
124  // generate nominal flux
125  //
126 
127  Axis_t ax = 0, ay = 0, az = 0;
128  fTotalFluxHisto->GetRandom3(ax, ay, az);
129  Ev = (double)ax;
130  costheta = (double)ay;
131  phi = (double)az;
132  nu_pdg = this->SelectNeutrino(Ev, costheta, phi);
133  weight = 1.0;
134  }
135 
136  // Compute etc trigonometric numbers
137  double sintheta = TMath::Sqrt(1-costheta*costheta);
138  double cosphi = TMath::Cos(phi);
139  double sinphi = TMath::Sin(phi);
140 
141  // Set the neutrino pdg code
142  fgPdgC = nu_pdg;
143 
144  // Set the neutrino weight
145  fWeight = weight;
146 
147  // Compute the neutrino momentum
148  // The `-1' means it is directed towards the detector.
149  double pz = -1.* Ev * costheta;
150  double py = -1.* Ev * sintheta * sinphi;
151  double px = -1.* Ev * sintheta * cosphi;
152 
153  // Default vertex is at the origin
154  double z = 0.0;
155  double y = 0.0;
156  double x = 0.0;
157 
158  // Shift the neutrino position onto the flux generation surface.
159  // The position is computed at the surface of a sphere with R=fRl
160  // at the topocentric horizontal (THZ) coordinate system.
161  if( fRl>0.0 ){
162  z += fRl * costheta;
163  y += fRl * sintheta * sinphi;
164  x += fRl * sintheta * cosphi;
165  }
166 
167  // Apply user-defined rotation from THZ -> user-defined topocentric
168  // coordinate system.
169  if( !fRotTHz2User.IsIdentity() )
170  {
171  TVector3 tx3(x, y, z );
172  TVector3 tp3(px,py,pz);
173 
174  tx3 = fRotTHz2User * tx3;
175  tp3 = fRotTHz2User * tp3;
176 
177  x = tx3.X();
178  y = tx3.Y();
179  z = tx3.Z();
180  px = tp3.X();
181  py = tp3.Y();
182  pz = tp3.Z();
183  }
184 
185  // If the position is left as is, then all generated neutrinos
186  // would point towards the origin.
187  // Displace the position randomly on the surface that is
188  // perpendicular to the selected point P(xo,yo,zo) on the sphere
189  if( fRt>0.0 ){
190  TVector3 vec(x,y,z); // vector towards selected point
191  TVector3 dvec1 = vec.Orthogonal(); // orthogonal vector
192  TVector3 dvec2 = dvec1; // second orthogonal vector
193  dvec2.Rotate(-kPi/2.0,vec); // rotate second vector by 90deg,
194  // now forming a new orthogonal cartesian coordinate system
195  double psi = 2.*kPi* rnd->RndFlux().Rndm(); // rndm angle [0,2pi]
196  double random = rnd->RndFlux().Rndm(); // rndm number [0,1]
197  dvec1.SetMag(TMath::Sqrt(random)*fRt*TMath::Cos(psi));
198  dvec2.SetMag(TMath::Sqrt(random)*fRt*TMath::Sin(psi));
199  x += dvec1.X() + dvec2.X();
200  y += dvec1.Y() + dvec2.Y();
201  z += dvec1.Z() + dvec2.Z();
202  }
203 
204  // Set the neutrino momentum and position 4-vectors with values
205  // calculated at previous steps.
206  fgP4.SetPxPyPzE(px, py, pz, Ev);
207  fgX4.SetXYZT (x, y, z, 0.);
208 
209  // Increment flux neutrino counter used for sample normalization purposes.
210  fNNeutrinos++;
211 
212  // Report and exit
213  LOG("Flux", pINFO)
214  << "Generated neutrino: "
215  << "\n pdg-code: " << fgPdgC
216  << "\n p4: " << utils::print::P4AsShortString(&fgP4)
217  << "\n x4: " << utils::print::X4AsString(&fgX4);
218 
219  return true;
220 }
221 //___________________________________________________________________________
222 long int GAtmoFlux::NFluxNeutrinos(void) const
223 {
224  return fNNeutrinos;
225 }
226 //___________________________________________________________________________
227 void GAtmoFlux::ForceMinEnergy(double emin)
228 {
229  emin = TMath::Max(0., emin);
230  fMinEvCut = emin;
231 }
232 //___________________________________________________________________________
233 void GAtmoFlux::ForceMaxEnergy(double emax)
234 {
235  emax = TMath::Max(0., emax);
236  fMaxEvCut = emax;
237 }
238 //___________________________________________________________________________
239 void GAtmoFlux::Clear(Option_t * opt)
240 {
241 // Dummy clear method needed to conform to GFluxI interface
242 //
243  LOG("Flux", pERROR) << "No clear method implemented for option:"<< opt;
244 }
245 //___________________________________________________________________________
246 void GAtmoFlux::GenerateWeighted(bool gen_weighted)
247 {
248  fGenWeighted = gen_weighted;
249 }
250 //___________________________________________________________________________
251 void GAtmoFlux:: SetSpectralIndex(double index)
252 {
253  if( index != 1.0 ){
254  fSpectralIndex = index;
255  }
256  else {
257  LOG("Flux", pWARN) << "Warning: cannot use a spectral index of unity";
258  }
259 
260  LOG("Flux", pNOTICE) << "Using Spectral Index = " << index;
261 }
262 //___________________________________________________________________________
263 void GAtmoFlux::SetUserCoordSystem(TRotation & rotation)
264 {
265  fRotTHz2User = rotation;
266 }
267 //___________________________________________________________________________
269 {
270  LOG("Flux", pNOTICE) << "Initializing atmospheric flux driver";
271 
272  fMaxEv = -1;
273 
274  fNumPhiBins = 0;
275  fNumCosThetaBins = 0;
276  fNumEnergyBins = 0;
277  fPhiBins = 0;
278  fCosThetaBins = 0;
279  fEnergyBins = 0;
280 
281  fTotalFluxHisto = 0;
282  fTotalFluxHistoIntg = 0;
283 
284  bool allow_dup = false;
285  fPdgCList = new PDGCodeList(allow_dup);
286 
287  // initializing flux TH3D histos [ flux = f(Ev,costheta,phi) ] & files
288  fFluxFile.clear();
289  fFluxHistoMap.clear();
290  fTotalFluxHisto = 0;
291  fTotalFluxHistoIntg = 0;
292 
293  // Default option is to generate unweighted flux neutrinos
294  // (flux = f(E,costheta) will be used as PDFs)
295  // User can enable option to generate weighted neutrinos
296  // (neutrinos will be generated uniformly over costheta,
297  // and using a power law function in neutrino energy.
298  // The input flux = f(E,costheta) will be used for calculating a weight).
299  // Using a weighted flux avoids statistical fluctuations at high energies.
300  fSpectralIndex = 2.0;
301 
302  // weighting switched off by default
303  this->GenerateWeighted(false);
304 
305  // Default: No min/max energy cut
306  this->ForceMinEnergy(0.);
307  this->ForceMaxEnergy(9999999999.);
308 
309  // Default radii
310  fRl = 0.0;
311  fRt = 0.0;
312 
313  // Default detector coord system: Topocentric Horizontal Coordinate system
314  fRotTHz2User.SetToIdentity();
315 
316  // Reset `current' selected flux neutrino
317  this->ResetSelection();
318 
319  // Reset number of neutrinos thrown so far
320  fNNeutrinos = 0;
321 
322  // Done!
323  fInitialized = 1;
324 }
325 //___________________________________________________________________________
326 void GAtmoFlux::ResetSelection(void)
327 {
328 // initializing running neutrino pdg-code, 4-position, 4-momentum
329 
330  fgPdgC = 0;
331  fgP4.SetPxPyPzE (0.,0.,0.,0.);
332  fgX4.SetXYZT (0.,0.,0.,0.);
333  fWeight = 0;
334 }
335 //___________________________________________________________________________
336 void GAtmoFlux::CleanUp(void)
337 {
338  LOG("Flux", pNOTICE) << "Cleaning up...";
339 
340  map<int,TH3D*>::iterator rawiter = fRawFluxHistoMap.begin();
341  for( ; rawiter != fRawFluxHistoMap.end(); ++rawiter) {
342  TH3D * flux_histogram = rawiter->second;
343  if(flux_histogram) {
344  delete flux_histogram;
345  flux_histogram = 0;
346  }
347  }
348  fRawFluxHistoMap.clear();
349 
350  map<int,TH3D*>::iterator iter = fFluxHistoMap.begin();
351  for( ; iter != fFluxHistoMap.end(); ++iter) {
352  TH3D * flux_histogram = iter->second;
353  if(flux_histogram) {
354  delete flux_histogram;
355  flux_histogram = 0;
356  }
357  }
358  fFluxHistoMap.clear();
359 
360  if (fTotalFluxHisto) delete fTotalFluxHisto;
361  if (fPdgCList) delete fPdgCList;
362 
363  if (fPhiBins ) { delete[] fPhiBins ; fPhiBins =NULL; }
364  if (fCosThetaBins) { delete[] fCosThetaBins; fCosThetaBins=NULL; }
365  if (fEnergyBins ) { delete[] fEnergyBins ; fEnergyBins =NULL; }
366 
367 }
368 //___________________________________________________________________________
369 void GAtmoFlux::SetRadii(double Rlongitudinal, double Rtransverse)
370 {
371  LOG ("Flux", pNOTICE) << "Setting R[longitudinal] = " << Rlongitudinal;
372  LOG ("Flux", pNOTICE) << "Setting R[transverse] = " << Rtransverse;
373 
374  fRl = Rlongitudinal;
375  fRt = Rtransverse;
376 }
377 
378 double GAtmoFlux::GetFluxSurfaceArea(void)
379 {
380  return kPi*pow(fRt,2);
381 }
382 
383 double GAtmoFlux::GetLongitudinalRadius(void)
384 {
385  return fRl;
386 }
387 
388 double GAtmoFlux::GetTransverseRadius(void)
389 {
390  return fRt;
391 }
392 
393 //___________________________________________________________________________
394 void GAtmoFlux::AddFluxFile(int nu_pdg, string filename)
395 {
396  // Check file exists
397  std::ifstream f(filename.c_str());
398  if (!f.good()) {
399  LOG("Flux", pFATAL) << "Flux file does not exist "<<filename;
400  exit(-1);
401  }
402  if ( pdg::IsNeutrino(nu_pdg) || pdg::IsAntiNeutrino(nu_pdg) ) {
403  fFluxFlavour.push_back(nu_pdg);
404  fFluxFile.push_back(filename);
405  } else {
406  LOG ("Flux", pWARN)
407  << "Input particle code: " << nu_pdg << " not a neutrino!";
408  }
409 }
410 //___________________________________________________________________________
411 void GAtmoFlux::AddFluxFile(string filename)
412 {
413 // FLUKA and BGLRS provide one file per neutrino species.
414 // HAKKKM provides a single file for all nue,nuebar,numu,numubar.
415 // If no neutrino species is provided, assume that the file contains all 4
416 // but fit it into the franework developed for FLUKA and BGLRS,
417 // i.e. add the file 4 times
418 
419 // Check file exists
420  std::ifstream f(filename.c_str());
421  if (!f.good()) {
422  LOG("Flux", pFATAL) << "Flux file does not exist "<<filename;
423  exit(-1);
424  }
425 
426  fFluxFlavour.push_back(kPdgNuE); fFluxFile.push_back(filename);
427  fFluxFlavour.push_back(kPdgAntiNuE); fFluxFile.push_back(filename);
428  fFluxFlavour.push_back(kPdgNuMu); fFluxFile.push_back(filename);
429  fFluxFlavour.push_back(kPdgAntiNuMu); fFluxFile.push_back(filename);
430 
431 }
432 //___________________________________________________________________________
433 //void GAtmoFlux::SetFluxFile(int nu_pdg, string filename)
434 //{
435 // return AddFluxFile( nu_pdg, filename );
436 //}
437 //___________________________________________________________________________
438 bool GAtmoFlux::LoadFluxData(void)
439 {
440  LOG("Flux", pNOTICE)
441  << "Loading atmospheric neutrino flux simulation data";
442 
443  fFluxHistoMap.clear();
444  fPdgCList->clear();
445 
446  bool loading_status = true;
447 
448  for( unsigned int n=0; n<fFluxFlavour.size(); n++ ){
449  int nu_pdg = fFluxFlavour.at(n);
450  string filename = fFluxFile.at(n);
451  string pname = PDGLibrary::Instance()->Find(nu_pdg)->GetName();
452 
453  LOG("Flux", pNOTICE) << "Loading data for: " << pname;
454 
455  // create histogram
456  TH3D* hist = 0;
457  std::map<int,TH3D*>::iterator myMapEntry = fRawFluxHistoMap.find(nu_pdg);
458  if( myMapEntry == fRawFluxHistoMap.end() ){
459  hist = this->CreateFluxHisto(pname.c_str(), pname.c_str());
460  fRawFluxHistoMap.insert( map<int,TH3D*>::value_type(nu_pdg,hist) );
461  }
462  // now let concrete instances to read the flux-specific data files
463  // and fill the histogram
464  bool loaded = this->FillFluxHisto(nu_pdg, filename);
465 
466  loading_status = loading_status && loaded;
467 
468  if (!loaded) {
469  LOG("Flux", pERROR)
470  << "Error loading atmospheric neutrino flux simulation data from " << filename;
471  break;
472  }
473  }
474 
475  if(loading_status) {
476  map<int,TH3D*>::iterator hist_iter = fRawFluxHistoMap.begin();
477  for ( ; hist_iter != fRawFluxHistoMap.end(); ++hist_iter) {
478  int nu_pdg = hist_iter->first;
479  TH3D* hist = hist_iter->second;
480 
481  TH3D* hnorm = this->CreateNormalisedFluxHisto( hist );
482  fFluxHistoMap.insert( map<int,TH3D*>::value_type(nu_pdg,hnorm) );
483  fPdgCList->push_back(nu_pdg);
484  }
485 
486  LOG("Flux", pNOTICE)
487  << "Atmospheric neutrino flux simulation data loaded!";
488  this->AddAllFluxes();
489  return true;
490  }
491 
492  LOG("Flux", pERROR)
493  << "Error loading atmospheric neutrino flux simulation data";
494  return false;
495 }
496 //___________________________________________________________________________
497 TH3D* GAtmoFlux::CreateNormalisedFluxHisto(TH3D* hist)
498 {
499 // return integrated flux
500 
501  // sanity check
502  if(!hist) return 0;
503 
504  // make new histogram name
505  TString histname = hist->GetName();
506  histname.Append("_IntegratedFlux");
507 
508  // make new histogram
509  TH3D* hist_intg = (TH3D*)(hist->Clone(histname.Data()));
510  hist_intg->Reset();
511 
512  // integrate flux in each bin
513  Double_t dN_dEdS = 0.0;
514  Double_t dS = 0.0;
515  Double_t dE = 0.0;
516  Double_t dN = 0.0;
517 
518  for(Int_t e_bin = 1; e_bin <= hist->GetXaxis()->GetNbins(); e_bin++)
519  {
520  for(Int_t c_bin = 1; c_bin <= hist->GetYaxis()->GetNbins(); c_bin++)
521  {
522  for(Int_t p_bin = 1; p_bin <= hist->GetZaxis()->GetNbins(); p_bin++)
523  {
524  dN_dEdS = hist->GetBinContent(e_bin,c_bin,p_bin);
525 
526  dE = hist->GetXaxis()->GetBinUpEdge(e_bin)
527  - hist->GetXaxis()->GetBinLowEdge(e_bin);
528 
529  dS = ( hist->GetZaxis()->GetBinUpEdge(p_bin)
530  - hist->GetZaxis()->GetBinLowEdge(p_bin) )
531  * ( hist->GetYaxis()->GetBinUpEdge(c_bin)
532  - hist->GetYaxis()->GetBinLowEdge(c_bin) );
533 
534  dN = dN_dEdS*dE*dS;
535 
536  hist_intg->SetBinContent(e_bin,c_bin,p_bin,dN);
537  }
538  }
539  }
540 
541  return hist_intg;
542 }
543 //___________________________________________________________________________
544 void GAtmoFlux::ZeroFluxHisto(TH3D * histo)
545 {
546  LOG("Flux", pDEBUG) << "Forcing flux histogram contents to 0";
547  if(!histo) return;
548  histo->Reset();
549 }
550 //___________________________________________________________________________
551 void GAtmoFlux::AddAllFluxes(void)
552 {
553  LOG("Flux", pNOTICE)
554  << "Computing combined flux & flux normalization factor";
555 
556  if(fTotalFluxHisto) delete fTotalFluxHisto;
557 
558  fTotalFluxHisto = this->CreateFluxHisto("sum", "combined flux" );
559 
560  map<int,TH3D*>::iterator it = fFluxHistoMap.begin();
561  for( ; it != fFluxHistoMap.end(); ++it) {
562  TH3D * flux_histogram = it->second;
563  fTotalFluxHisto->Add(flux_histogram);
564  }
565 
566  fTotalFluxHistoIntg = fTotalFluxHisto->Integral();
567 }
568 //___________________________________________________________________________
569 TH3D * GAtmoFlux::CreateFluxHisto(string name, string title)
570 {
571  LOG("Flux", pNOTICE) << "Instantiating histogram: [" << name << "]";
572  TH3D * hist = new TH3D(
573  name.c_str(), title.c_str(),
574  fNumEnergyBins, fEnergyBins,
575  fNumCosThetaBins, fCosThetaBins,
576  fNumPhiBins, fPhiBins);
577  return hist;
578 }
579 //___________________________________________________________________________
580 int GAtmoFlux::SelectNeutrino(double Ev, double costheta, double phi)
581 {
582 // Select a neutrino species at the input Ev and costheta given their
583 // relatve flux at this bin.
584 // Returns a neutrino PDG code
585 
586  unsigned int n = fPdgCList->size();
587  double * flux = new double[n];
588 
589  unsigned int i=0;
590  map<int,TH3D*>::iterator it = fFluxHistoMap.begin();
591  for( ; it != fFluxHistoMap.end(); ++it) {
592  TH3D * flux_histogram = it->second;
593  int ibin = flux_histogram->FindBin(Ev,costheta,phi);
594  flux[i] = flux_histogram->GetBinContent(ibin);
595  i++;
596  }
597  double flux_sum = 0;
598  for(i=0; i<n; i++) {
599  flux_sum += flux[i];
600  flux[i] = flux_sum;
601  LOG("Flux", pDEBUG)
602  << "Sum{Flux(0->" << i <<")} = " << flux[i];
603  }
604 
605  RandomGen * rnd = RandomGen::Instance();
606  double R = flux_sum * rnd->RndFlux().Rndm();
607 
608  LOG("Flux", pDEBUG) << "R = " << R;
609 
610  for(i=0; i<n; i++) {
611  if( R < flux[i] ) {
612  delete [] flux;
613  int selected_pdg = (*fPdgCList)[i];
614  LOG("Flux", pINFO)
615  << "Selected neutrino PDG code = " << selected_pdg;
616  return selected_pdg;
617  }
618  }
619 
620  // shouldn't be here
621  LOG("Flux", pERROR) << "Could not select a neutrino species!";
622  assert(false);
623 
624  return -1;
625 }
626 
627 //___________________________________________________________________________
628 TH3D* GAtmoFlux::GetFluxHistogram(int flavour)
629 {
630  TH3D* histogram = 0;
631  std::map<int,TH3D*>::iterator it = fRawFluxHistoMap.find(flavour);
632  if(it != fRawFluxHistoMap.end())
633  {
634  histogram = it->second;
635  }
636  return histogram;
637 }
638 
639 /* Returns the total integrated flux in units of 1/(m^2 s). */
640 double GAtmoFlux::GetTotalFlux(void)
641 {
642  double flux = 0.0;
643  map<int,TH3D*>::iterator rawiter;
644 
645  rawiter = fRawFluxHistoMap.begin();
646  for (; rawiter != fRawFluxHistoMap.end(); ++rawiter) {
647  TH3D *h = rawiter->second;
648  if (h) {
649  flux += h->Integral("width");
650  LOG("Flux", pDEBUG) << "Total flux for " << rawiter->first << " equals " << h->Integral("width") << ".";
651  }
652  }
653 
654  return flux;
655 }
656 
657 /* Returns the total integrated flux in units of * 1/(m^2 s) between the
658  * minimum and maximum energy. */
659 double GAtmoFlux::GetTotalFluxInEnergyRange(void)
660 {
661  double flux = 0.0;
662  map<int,TH3D*>::iterator rawiter;
663  int e_min_bin, e_max_bin;
664  double Emin, Emax;
665 
666  Emin = this->MinEnergy();
667  Emax = this->MaxEnergy();
668 
669  if (Emax < Emin) {
670  LOG("Flux", pFATAL) << "Emax = " << Emax << " is less than Emin = " << Emin;
671  exit(-1);
672  }
673 
674  rawiter = fRawFluxHistoMap.begin();
675  for (; rawiter != fRawFluxHistoMap.end(); ++rawiter) {
676  TH3D *h = rawiter->second;
677 
678  if (!h) continue;
679 
680  /* Get the bins containing `emin` and `emax`. */
681  e_min_bin = h->GetXaxis()->FindBin(Emin);
682  e_max_bin = h->GetXaxis()->FindBin(Emax);
683 
684  if (e_min_bin > h->GetXaxis()->GetNbins()) {
685  /* If the minimum bin is past the end, continue. */
686  continue;
687  } else if (e_min_bin == e_max_bin) {
688  /* If they both end up in the same bin, we just take the total bin
689  * contents in that energy bin and multiply by the difference between
690  * the energies. */
691  flux += h->Integral(e_min_bin,e_min_bin,1,h->GetYaxis()->GetNbins(),1,h->GetZaxis()->GetNbins(),"width")*(Emax - Emin)/(h->GetXaxis()->GetBinUpEdge(e_min_bin)-h->GetXaxis()->GetBinLowEdge(e_min_bin));
692  } else {
693  /* First we calculate the integral within that bin from `emin` to the top
694  * edge of the bin. */
695  if (e_min_bin > 0)
696  flux += h->Integral(e_min_bin,e_min_bin,1,h->GetYaxis()->GetNbins(),1,h->GetZaxis()->GetNbins(),"width")*(h->GetXaxis()->GetBinUpEdge(e_min_bin) - Emin)/(h->GetXaxis()->GetBinUpEdge(e_min_bin)-h->GetXaxis()->GetBinLowEdge(e_min_bin));
697  /* Next, we calculate the integral for all the bins between the min and
698  * max bin. */
699  if (e_min_bin < h->GetXaxis()->GetNbins())
700  flux += h->Integral(e_min_bin+1,e_max_bin-1,1,h->GetYaxis()->GetNbins(),1,h->GetZaxis()->GetNbins(),"width");
701  /* Finally, we calculate the integral for the last bin. */
702  if (e_max_bin <= h->GetXaxis()->GetNbins())
703  flux += h->Integral(e_max_bin,e_max_bin,1,h->GetYaxis()->GetNbins(),1,h->GetZaxis()->GetNbins(),"width")*(Emax - h->GetXaxis()->GetBinLowEdge(e_max_bin))/(h->GetXaxis()->GetBinUpEdge(e_max_bin)-h->GetXaxis()->GetBinLowEdge(e_max_bin));
704  }
705  }
706 
707  return flux;
708 }
709 
710 //___________________________________________________________________________
711 double GAtmoFlux::GetFlux(int flavour)
712 {
713  TH3D* flux_hist = this->GetFluxHistogram(flavour);
714  if(!flux_hist) return 0.0;
715 
716  Double_t flux = 0.0;
717  Double_t dN_dEdS = 0.0;
718  Double_t dS = 0.0;
719  Double_t dE = 0.0;
720 
721  for(Int_t e_bin = 1; e_bin <= flux_hist->GetXaxis()->GetNbins(); e_bin++)
722  {
723  for(Int_t c_bin = 1; c_bin <= flux_hist->GetYaxis()->GetNbins(); c_bin++)
724  {
725  for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
726  {
727  dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
728 
729  dE = flux_hist->GetXaxis()->GetBinUpEdge(e_bin)
730  - flux_hist->GetXaxis()->GetBinLowEdge(e_bin);
731 
732  dS = ( flux_hist->GetZaxis()->GetBinUpEdge(p_bin)
733  - flux_hist->GetZaxis()->GetBinLowEdge(p_bin) )
734  * ( flux_hist->GetYaxis()->GetBinUpEdge(c_bin)
735  - flux_hist->GetYaxis()->GetBinLowEdge(c_bin) );
736 
737  flux += dN_dEdS*dE*dS;
738  }
739  }
740  }
741 
742  return flux;
743 }
744 //___________________________________________________________________________
745 double GAtmoFlux::GetFlux(int flavour, double energy)
746 {
747  TH3D* flux_hist = this->GetFluxHistogram(flavour);
748  if(!flux_hist) return 0.0;
749 
750  Double_t flux = 0.0;
751  Double_t dN_dEdS = 0.0;
752  Double_t dS = 0.0;
753 
754  Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);
755 
756  for(Int_t c_bin = 1; c_bin <= flux_hist->GetYaxis()->GetNbins(); c_bin++)
757  {
758  for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
759  {
760  dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
761 
762  dS = ( flux_hist->GetZaxis()->GetBinUpEdge(p_bin)
763  - flux_hist->GetZaxis()->GetBinLowEdge(p_bin) )
764  * ( flux_hist->GetYaxis()->GetBinUpEdge(c_bin)
765  - flux_hist->GetYaxis()->GetBinLowEdge(c_bin) );
766 
767  flux += dN_dEdS*dS;
768  }
769  }
770 
771  return flux;
772 }
773 //___________________________________________________________________________
774 double GAtmoFlux::GetFlux(int flavour, double energy, double costh)
775 {
776  TH3D* flux_hist = this->GetFluxHistogram(flavour);
777  if(!flux_hist) return 0.0;
778 
779  Double_t flux = 0.0;
780  Double_t dN_dEdS = 0.0;
781  Double_t dS = 0.0;
782 
783  Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);
784  Int_t c_bin = flux_hist->GetYaxis()->FindBin(costh);
785 
786  for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
787  {
788  dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
789 
790  dS = ( flux_hist->GetZaxis()->GetBinUpEdge(p_bin)
791  - flux_hist->GetZaxis()->GetBinLowEdge(p_bin) );
792 
793  flux += dN_dEdS*dS;
794  }
795 
796  return flux;
797 }
798 //___________________________________________________________________________
799 double GAtmoFlux::GetFlux(int flavour, double energy, double costh, double phi)
800 {
801  TH3D* flux_hist = this->GetFluxHistogram(flavour);
802  if(!flux_hist) return 0.0;
803 
804  Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);
805  Int_t c_bin = flux_hist->GetYaxis()->FindBin(costh);
806  Int_t p_bin = flux_hist->GetZaxis()->FindBin(phi);
807 
808  Double_t flux = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
809  return flux;
810 }
811 //___________________________________________________________________________
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:45
const int kPdgNuE
Definition: PDGCodes.h:28
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
#define pERROR
Definition: Messenger.h:59
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
const int kPdgAntiNuE
Definition: PDGCodes.h:29
#define pFATAL
Definition: Messenger.h:56
const int kPdgNuMu
Definition: PDGCodes.h:30
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
A list of PDG codes.
Definition: PDGCodeList.h:32
GAtmoFlux * GetFlux(void)
Definition: gAtmoEvGen.cxx:505
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
void Initialize(void)
#define pINFO
Definition: Messenger.h:62
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
#define pWARN
Definition: Messenger.h:60
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
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 pDEBUG
Definition: Messenger.h:63