GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
genie::flux::GAtmoFlux Class Referenceabstract

A base class for the FLUKA, BGLRS and ATMNC atmo. nu. flux drivers. The driver depends on data files provided by the atmospheric neutrino flux simulation authors in order to determine the angular and energy dependence for each neutrino species. The position of each flux neutrino [going towards a detector centered at (0,0,0)] is generated uniformly on a plane that is perpendicular to a sphere of radius Rl at the point that is determined by the generated neutrino direction (theta,phi). The size of the area of that plane, where flux neutrinos are generated, is determined by the transverse radius Rt. You can tweak Rl, Rt to match the size of your detector. Initially, neutrino coordinates are generated in a default detector coordinate system (Topocentric Horizontal Coordinate -THZ-): +z: Points towards the local zenith. +x: On same plane as local meridian, pointing south. +y: As needed to make a right-handed coordinate system. origin: detector centre Alternative user-defined topocentric systems can be defined by specifying the appropriate rotation from THZ. The driver allows minimum and maximum energy cuts. Also it provides the options to generate wither unweighted or weighted flux neutrinos (the latter giving smoother distributions at the tails). More...

#include <GAtmoFlux.h>

Inheritance diagram for genie::flux::GAtmoFlux:
Inheritance graph
[legend]
Collaboration diagram for genie::flux::GAtmoFlux:
Collaboration graph
[legend]

Public Member Functions

virtual ~GAtmoFlux ()
 
virtual const PDGCodeListFluxParticles (void)
 declare list of flux neutrinos that can be generated (for init. purposes) More...
 
virtual double MaxEnergy (void)
 declare the max flux neutrino energy that can be generated (for init. purposes) More...
 
virtual bool GenerateNext (void)
 generate the next flux neutrino (return false in err) More...
 
virtual int PdgCode (void)
 returns the flux neutrino pdg code More...
 
virtual double Weight (void)
 returns the flux neutrino weight (if any) More...
 
virtual const TLorentzVector & Momentum (void)
 returns the flux neutrino 4-momentum More...
 
virtual const TLorentzVector & Position (void)
 returns the flux neutrino 4-position (note: expect SI rather than physical units) More...
 
virtual bool End (void)
 true if no more flux nu's can be thrown (eg reaching end of beam sim ntuples) More...
 
virtual long int Index (void)
 returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current entry number) More...
 
virtual void Clear (Option_t *opt)
 reset state variables based on opt More...
 
virtual void GenerateWeighted (bool gen_weighted)
 set whether to generate weighted or unweighted neutrinos More...
 
double Enu (void)
 
double Energy (void)
 
double CosTheta (void)
 
double CosZenith (void)
 
long int NFluxNeutrinos (void) const
 Number of flux nu's generated. Not the same as the number of nu's thrown towards the geometry (if there are cuts). More...
 
void ForceMinEnergy (double emin)
 
void ForceMaxEnergy (double emax)
 
void SetSpectralIndex (double index)
 
void SetRadii (double Rlongitudinal, double Rtransverse)
 
double GetFluxSurfaceArea (void)
 
double GetLongitudinalRadius (void)
 
double GetTransverseRadius (void)
 
void SetUserCoordSystem (TRotation &rotation)
 Rotation: Topocentric Horizontal -> User-defined Topocentric Coord System. More...
 
void AddFluxFile (int neutrino_pdg, string filename)
 
void AddFluxFile (string filename)
 
bool LoadFluxData (void)
 
TH3D * GetFluxHistogram (int flavour)
 
double GetTotalFlux (void)
 
double GetTotalFluxInEnergyRange (void)
 
double GetFlux (int flavour)
 
double GetFlux (int flavour, double energy)
 
double GetFlux (int flavour, double energy, double costh)
 
double GetFlux (int flavour, double energy, double costh, double phi)
 
- Public Member Functions inherited from genie::GFluxI
virtual ~GFluxI ()
 

Protected Member Functions

 GAtmoFlux ()
 
bool GenerateNext_1try (void)
 
void Initialize (void)
 
void CleanUp (void)
 
void ResetSelection (void)
 
double MinEnergy (void)
 
TH3D * CreateFluxHisto (string name, string title)
 
void ZeroFluxHisto (TH3D *hist)
 
void AddAllFluxes (void)
 
int SelectNeutrino (double Ev, double costheta, double phi)
 
TH3D * CreateNormalisedFluxHisto (TH3D *hist)
 
virtual bool FillFluxHisto (int nu_pdg, string filename)=0
 
- Protected Member Functions inherited from genie::GFluxI
 GFluxI ()
 

Protected Attributes

double fMaxEv
 maximum energy (in input flux files) More...
 
PDGCodeListfPdgCList
 input list of neutrino pdg-codes More...
 
int fgPdgC
 current generated nu pdg-code More...
 
TLorentzVector fgP4
 current generated nu 4-momentum More...
 
TLorentzVector fgX4
 current generated nu 4-position More...
 
double fWeight
 current generated nu weight More...
 
long int fNNeutrinos
 number of flux neutrinos thrown so far More...
 
double fMaxEvCut
 user-defined cut: maximum energy More...
 
double fMinEvCut
 user-defined cut: minimum energy More...
 
double fRl
 defining flux neutrino generation surface: longitudinal radius More...
 
double fRt
 defining flux neutrino generation surface: transverse radius More...
 
TRotation fRotTHz2User
 coord. system rotation: THZ -> Topocentric user-defined More...
 
unsigned int fNumPhiBins
 number of phi bins in input flux data files More...
 
unsigned int fNumCosThetaBins
 number of cos(theta) bins in input flux data files More...
 
unsigned int fNumEnergyBins
 number of energy bins in input flux data files More...
 
double * fPhiBins
 phi bins in input flux data files More...
 
double * fCosThetaBins
 cos(theta) bins in input flux data files More...
 
double * fEnergyBins
 energy bins in input flux data files More...
 
bool fGenWeighted
 generate a weighted or unweighted flux? More...
 
double fSpectralIndex
 power law function used for weighted flux More...
 
bool fInitialized
 flag to check that initialization is run More...
 
TH3D * fTotalFluxHisto
 flux = f(Ev,cos8,phi) summed over neutrino species More...
 
double fTotalFluxHistoIntg
 fFluxSum2D integral More...
 
map< int, TH3D * > fFluxHistoMap
 flux = f(Ev,cos8,phi) for each neutrino species More...
 
map< int, TH3D * > fRawFluxHistoMap
 flux = f(Ev,cos8,phi) for each neutrino species More...
 
vector< int > fFluxFlavour
 input flux file for each neutrino species More...
 
vector< string > fFluxFile
 input flux file for each neutrino species More...
 

Detailed Description

A base class for the FLUKA, BGLRS and ATMNC atmo. nu. flux drivers. The driver depends on data files provided by the atmospheric neutrino flux simulation authors in order to determine the angular and energy dependence for each neutrino species. The position of each flux neutrino [going towards a detector centered at (0,0,0)] is generated uniformly on a plane that is perpendicular to a sphere of radius Rl at the point that is determined by the generated neutrino direction (theta,phi). The size of the area of that plane, where flux neutrinos are generated, is determined by the transverse radius Rt. You can tweak Rl, Rt to match the size of your detector. Initially, neutrino coordinates are generated in a default detector coordinate system (Topocentric Horizontal Coordinate -THZ-): +z: Points towards the local zenith. +x: On same plane as local meridian, pointing south. +y: As needed to make a right-handed coordinate system. origin: detector centre Alternative user-defined topocentric systems can be defined by specifying the appropriate rotation from THZ. The driver allows minimum and maximum energy cuts. Also it provides the options to generate wither unweighted or weighted flux neutrinos (the latter giving smoother distributions at the tails).

Author
Costas Andreopoulos <c.andreopoulos cern.ch> University of Liverpool
Created:
January 26, 2008
License:
Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 60 of file GAtmoFlux.h.

Constructor & Destructor Documentation

GAtmoFlux::~GAtmoFlux ( )
virtual

Definition at line 39 of file GAtmoFlux.cxx.

40 {
41  this->CleanUp();
42 }
GAtmoFlux::GAtmoFlux ( )
protected

Definition at line 34 of file GAtmoFlux.cxx.

35 {
36  fInitialized = 0;
37 }
bool fInitialized
flag to check that initialization is run
Definition: GAtmoFlux.h:151

Member Function Documentation

void GAtmoFlux::AddAllFluxes ( void  )
protected

Definition at line 551 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

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 }
map< int, TH3D * > fFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:154
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double fTotalFluxHistoIntg
fFluxSum2D integral
Definition: GAtmoFlux.h:153
TH3D * fTotalFluxHisto
flux = f(Ev,cos8,phi) summed over neutrino species
Definition: GAtmoFlux.h:152
TH3D * CreateFluxHisto(string name, string title)
Definition: GAtmoFlux.cxx:569
#define pNOTICE
Definition: Messenger.h:61
void GAtmoFlux::AddFluxFile ( int  neutrino_pdg,
string  filename 
)

Definition at line 394 of file GAtmoFlux.cxx.

References genie::pdg::IsAntiNeutrino(), genie::pdg::IsNeutrino(), LOG, pFATAL, and pWARN.

Referenced by GetFlux(), testGetTotalFlux(), and testGetTotalFluxInEnergyRange().

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 }
vector< int > fFluxFlavour
input flux file for each neutrino species
Definition: GAtmoFlux.h:156
vector< string > fFluxFile
input flux file for each neutrino species
Definition: GAtmoFlux.h:157
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
#define pFATAL
Definition: Messenger.h:56
#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
#define pWARN
Definition: Messenger.h:60
void GAtmoFlux::AddFluxFile ( string  filename)

Definition at line 411 of file GAtmoFlux.cxx.

References genie::kPdgAntiNuE, genie::kPdgAntiNuMu, genie::kPdgNuE, genie::kPdgNuMu, LOG, and pFATAL.

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 }
vector< int > fFluxFlavour
input flux file for each neutrino species
Definition: GAtmoFlux.h:156
vector< string > fFluxFile
input flux file for each neutrino species
Definition: GAtmoFlux.h:157
const int kPdgNuE
Definition: PDGCodes.h:28
const int kPdgAntiNuE
Definition: PDGCodes.h:29
#define pFATAL
Definition: Messenger.h:56
const int kPdgNuMu
Definition: PDGCodes.h:30
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
void GAtmoFlux::CleanUp ( void  )
protected

Definition at line 336 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

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 }
map< int, TH3D * > fFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:154
double * fPhiBins
phi bins in input flux data files
Definition: GAtmoFlux.h:146
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double * fCosThetaBins
cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:147
double * fEnergyBins
energy bins in input flux data files
Definition: GAtmoFlux.h:148
map< int, TH3D * > fRawFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:155
TH3D * fTotalFluxHisto
flux = f(Ev,cos8,phi) summed over neutrino species
Definition: GAtmoFlux.h:152
PDGCodeList * fPdgCList
input list of neutrino pdg-codes
Definition: GAtmoFlux.h:132
#define pNOTICE
Definition: Messenger.h:61
void GAtmoFlux::Clear ( Option_t *  opt)
virtual

reset state variables based on opt

Implements genie::GFluxI.

Definition at line 239 of file GAtmoFlux.cxx.

References LOG, and pERROR.

240 {
241 // Dummy clear method needed to conform to GFluxI interface
242 //
243  LOG("Flux", pERROR) << "No clear method implemented for option:"<< opt;
244 }
#define pERROR
Definition: Messenger.h:59
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double genie::flux::GAtmoFlux::CosTheta ( void  )
inline

Definition at line 81 of file GAtmoFlux.h.

References fgP4.

81 { return -fgP4.Pz()/fgP4.Energy(); }
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:134
double genie::flux::GAtmoFlux::CosZenith ( void  )
inline

Definition at line 82 of file GAtmoFlux.h.

References fgP4.

82 { return -fgP4.Pz()/fgP4.Energy(); }
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:134
TH3D * GAtmoFlux::CreateFluxHisto ( string  name,
string  title 
)
protected

Definition at line 569 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

570 {
571  LOG("Flux", pNOTICE) << "Instantiating histogram: [" << name << "]";
572  TH3D * hist = new TH3D(
573  name.c_str(), title.c_str(),
577  return hist;
578 }
double * fPhiBins
phi bins in input flux data files
Definition: GAtmoFlux.h:146
unsigned int fNumEnergyBins
number of energy bins in input flux data files
Definition: GAtmoFlux.h:145
unsigned int fNumCosThetaBins
number of cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:144
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double * fCosThetaBins
cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:147
unsigned int fNumPhiBins
number of phi bins in input flux data files
Definition: GAtmoFlux.h:143
double * fEnergyBins
energy bins in input flux data files
Definition: GAtmoFlux.h:148
const char * name
#define pNOTICE
Definition: Messenger.h:61
TH3D * GAtmoFlux::CreateNormalisedFluxHisto ( TH3D *  hist)
protected

Definition at line 497 of file GAtmoFlux.cxx.

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 }
virtual bool genie::flux::GAtmoFlux::End ( void  )
inlinevirtual

true if no more flux nu's can be thrown (eg reaching end of beam sim ntuples)

Implements genie::GFluxI.

Definition at line 73 of file GAtmoFlux.h.

73 { return false; }
double genie::flux::GAtmoFlux::Energy ( void  )
inline

Definition at line 80 of file GAtmoFlux.h.

References fgP4.

80 { return fgP4.Energy(); }
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:134
double genie::flux::GAtmoFlux::Enu ( void  )
inline

Definition at line 79 of file GAtmoFlux.h.

References fgP4.

79 { return fgP4.Energy(); }
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:134
virtual bool genie::flux::GAtmoFlux::FillFluxHisto ( int  nu_pdg,
string  filename 
)
protectedpure virtual
virtual const PDGCodeList& genie::flux::GAtmoFlux::FluxParticles ( void  )
inlinevirtual

declare list of flux neutrinos that can be generated (for init. purposes)

Implements genie::GFluxI.

Definition at line 66 of file GAtmoFlux.h.

References fPdgCList.

66 { return *fPdgCList; }
PDGCodeList * fPdgCList
input list of neutrino pdg-codes
Definition: GAtmoFlux.h:132
void GAtmoFlux::ForceMaxEnergy ( double  emax)

Definition at line 233 of file GAtmoFlux.cxx.

Referenced by GetFlux(), testGetTotalFlux(), and testGetTotalFluxInEnergyRange().

234 {
235  emax = TMath::Max(0., emax);
236  fMaxEvCut = emax;
237 }
double fMaxEvCut
user-defined cut: maximum energy
Definition: GAtmoFlux.h:138
void GAtmoFlux::ForceMinEnergy ( double  emin)

Definition at line 227 of file GAtmoFlux.cxx.

Referenced by GetFlux(), testGetTotalFlux(), and testGetTotalFluxInEnergyRange().

228 {
229  emin = TMath::Max(0., emin);
230  fMinEvCut = emin;
231 }
double fMinEvCut
user-defined cut: minimum energy
Definition: GAtmoFlux.h:139
bool GAtmoFlux::GenerateNext ( void  )
virtual

generate the next flux neutrino (return false in err)

Implements genie::GFluxI.

Definition at line 49 of file GAtmoFlux.cxx.

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 }
virtual double Weight(void)
returns the flux neutrino weight (if any)
Definition: GAtmoFlux.h:70
virtual double MaxEnergy(void)
declare the max flux neutrino energy that can be generated (for init. purposes)
Definition: GAtmoFlux.cxx:44
double MinEnergy(void)
Definition: GAtmoFlux.h:120
virtual const TLorentzVector & Momentum(void)
returns the flux neutrino 4-momentum
Definition: GAtmoFlux.h:71
bool GenerateNext_1try(void)
Definition: GAtmoFlux.cxx:71
bool GAtmoFlux::GenerateNext_1try ( void  )
protected

Definition at line 71 of file GAtmoFlux.cxx.

References GetFlux(), genie::RandomGen::Instance(), genie::constants::kPi, LOG, genie::utils::print::P4AsShortString(), pINFO, genie::RandomGen::RndFlux(), and genie::utils::print::X4AsString().

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 }
int SelectNeutrino(double Ev, double costheta, double phi)
Definition: GAtmoFlux.cxx:580
int fgPdgC
current generated nu pdg-code
Definition: GAtmoFlux.h:133
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:45
unsigned int fNumEnergyBins
number of energy bins in input flux data files
Definition: GAtmoFlux.h:145
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
double fSpectralIndex
power law function used for weighted flux
Definition: GAtmoFlux.h:150
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GAtmoFlux.h:137
double fRl
defining flux neutrino generation surface: longitudinal radius
Definition: GAtmoFlux.h:140
bool fInitialized
flag to check that initialization is run
Definition: GAtmoFlux.h:151
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
double GetFlux(int flavour)
Definition: GAtmoFlux.cxx:711
TRotation fRotTHz2User
coord. system rotation: THZ -&gt; Topocentric user-defined
Definition: GAtmoFlux.h:142
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool fGenWeighted
generate a weighted or unweighted flux?
Definition: GAtmoFlux.h:149
#define pINFO
Definition: Messenger.h:62
double * fEnergyBins
energy bins in input flux data files
Definition: GAtmoFlux.h:148
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:134
double fRt
defining flux neutrino generation surface: transverse radius
Definition: GAtmoFlux.h:141
TH3D * fTotalFluxHisto
flux = f(Ev,cos8,phi) summed over neutrino species
Definition: GAtmoFlux.h:152
void ResetSelection(void)
Definition: GAtmoFlux.cxx:326
PDGCodeList * fPdgCList
input list of neutrino pdg-codes
Definition: GAtmoFlux.h:132
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:71
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:57
TLorentzVector fgX4
current generated nu 4-position
Definition: GAtmoFlux.h:135
double fWeight
current generated nu weight
Definition: GAtmoFlux.h:136
void GAtmoFlux::GenerateWeighted ( bool  gen_weighted)
virtual

set whether to generate weighted or unweighted neutrinos

Implements genie::GFluxI.

Definition at line 246 of file GAtmoFlux.cxx.

247 {
248  fGenWeighted = gen_weighted;
249 }
bool fGenWeighted
generate a weighted or unweighted flux?
Definition: GAtmoFlux.h:149
double GAtmoFlux::GetFlux ( int  flavour)

Definition at line 711 of file GAtmoFlux.cxx.

Referenced by testGetTotalFlux(), and testGetTotalFluxInEnergyRange().

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 }
TH3D * GetFluxHistogram(int flavour)
Definition: GAtmoFlux.cxx:628
double GAtmoFlux::GetFlux ( int  flavour,
double  energy 
)

Definition at line 745 of file GAtmoFlux.cxx.

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 }
TH3D * GetFluxHistogram(int flavour)
Definition: GAtmoFlux.cxx:628
double GAtmoFlux::GetFlux ( int  flavour,
double  energy,
double  costh 
)

Definition at line 774 of file GAtmoFlux.cxx.

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 }
TH3D * GetFluxHistogram(int flavour)
Definition: GAtmoFlux.cxx:628
double GAtmoFlux::GetFlux ( int  flavour,
double  energy,
double  costh,
double  phi 
)

Definition at line 799 of file GAtmoFlux.cxx.

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 }
TH3D * GetFluxHistogram(int flavour)
Definition: GAtmoFlux.cxx:628
TH3D * GAtmoFlux::GetFluxHistogram ( int  flavour)

Definition at line 628 of file GAtmoFlux.cxx.

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 }
map< int, TH3D * > fRawFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:155
double GAtmoFlux::GetFluxSurfaceArea ( void  )

Definition at line 378 of file GAtmoFlux.cxx.

References genie::constants::kPi.

Referenced by main().

379 {
380  return kPi*pow(fRt,2);
381 }
double fRt
defining flux neutrino generation surface: transverse radius
Definition: GAtmoFlux.h:141
double GAtmoFlux::GetLongitudinalRadius ( void  )

Definition at line 383 of file GAtmoFlux.cxx.

384 {
385  return fRl;
386 }
double fRl
defining flux neutrino generation surface: longitudinal radius
Definition: GAtmoFlux.h:140
double GAtmoFlux::GetTotalFlux ( void  )

Definition at line 640 of file GAtmoFlux.cxx.

References LOG, and pDEBUG.

Referenced by testGetTotalFlux(), and testGetTotalFluxInEnergyRange().

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 }
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
map< int, TH3D * > fRawFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:155
#define pDEBUG
Definition: Messenger.h:63
double GAtmoFlux::GetTotalFluxInEnergyRange ( void  )

Definition at line 659 of file GAtmoFlux.cxx.

References LOG, and pFATAL.

Referenced by main(), and testGetTotalFluxInEnergyRange().

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 }
#define pFATAL
Definition: Messenger.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
virtual double MaxEnergy(void)
declare the max flux neutrino energy that can be generated (for init. purposes)
Definition: GAtmoFlux.cxx:44
double MinEnergy(void)
Definition: GAtmoFlux.h:120
map< int, TH3D * > fRawFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:155
double GAtmoFlux::GetTransverseRadius ( void  )

Definition at line 388 of file GAtmoFlux.cxx.

389 {
390  return fRt;
391 }
double fRt
defining flux neutrino generation surface: transverse radius
Definition: GAtmoFlux.h:141
virtual long int genie::flux::GAtmoFlux::Index ( void  )
inlinevirtual

returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current entry number)

Implements genie::GFluxI.

Definition at line 74 of file GAtmoFlux.h.

74 { return -1; }
void GAtmoFlux::Initialize ( void  )
protected

Definition at line 268 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

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;
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;
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 }
map< int, TH3D * > fFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:154
double * fPhiBins
phi bins in input flux data files
Definition: GAtmoFlux.h:146
vector< string > fFluxFile
input flux file for each neutrino species
Definition: GAtmoFlux.h:157
unsigned int fNumEnergyBins
number of energy bins in input flux data files
Definition: GAtmoFlux.h:145
unsigned int fNumCosThetaBins
number of cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:144
double fSpectralIndex
power law function used for weighted flux
Definition: GAtmoFlux.h:150
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GAtmoFlux.h:137
double fRl
defining flux neutrino generation surface: longitudinal radius
Definition: GAtmoFlux.h:140
bool fInitialized
flag to check that initialization is run
Definition: GAtmoFlux.h:151
virtual void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
Definition: GAtmoFlux.cxx:246
TRotation fRotTHz2User
coord. system rotation: THZ -&gt; Topocentric user-defined
Definition: GAtmoFlux.h:142
A list of PDG codes.
Definition: PDGCodeList.h:32
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void ForceMaxEnergy(double emax)
Definition: GAtmoFlux.cxx:233
double fTotalFluxHistoIntg
fFluxSum2D integral
Definition: GAtmoFlux.h:153
double * fCosThetaBins
cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:147
unsigned int fNumPhiBins
number of phi bins in input flux data files
Definition: GAtmoFlux.h:143
double * fEnergyBins
energy bins in input flux data files
Definition: GAtmoFlux.h:148
double fRt
defining flux neutrino generation surface: transverse radius
Definition: GAtmoFlux.h:141
TH3D * fTotalFluxHisto
flux = f(Ev,cos8,phi) summed over neutrino species
Definition: GAtmoFlux.h:152
void ResetSelection(void)
Definition: GAtmoFlux.cxx:326
PDGCodeList * fPdgCList
input list of neutrino pdg-codes
Definition: GAtmoFlux.h:132
#define pNOTICE
Definition: Messenger.h:61
void ForceMinEnergy(double emin)
Definition: GAtmoFlux.cxx:227
double fMaxEv
maximum energy (in input flux files)
Definition: GAtmoFlux.h:131
bool GAtmoFlux::LoadFluxData ( void  )

Definition at line 438 of file GAtmoFlux.cxx.

References genie::PDGLibrary::Find(), genie::PDGLibrary::Instance(), LOG, pERROR, and pNOTICE.

Referenced by GetFlux(), testGetTotalFlux(), and testGetTotalFluxInEnergyRange().

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 }
vector< int > fFluxFlavour
input flux file for each neutrino species
Definition: GAtmoFlux.h:156
map< int, TH3D * > fFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:154
vector< string > fFluxFile
input flux file for each neutrino species
Definition: GAtmoFlux.h:157
#define pERROR
Definition: Messenger.h:59
TH3D * CreateNormalisedFluxHisto(TH3D *hist)
Definition: GAtmoFlux.cxx:497
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
virtual bool FillFluxHisto(int nu_pdg, string filename)=0
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
map< int, TH3D * > fRawFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:155
TH3D * CreateFluxHisto(string name, string title)
Definition: GAtmoFlux.cxx:569
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
PDGCodeList * fPdgCList
input list of neutrino pdg-codes
Definition: GAtmoFlux.h:132
#define pNOTICE
Definition: Messenger.h:61
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
double GAtmoFlux::MaxEnergy ( void  )
virtual

declare the max flux neutrino energy that can be generated (for init. purposes)

Implements genie::GFluxI.

Definition at line 44 of file GAtmoFlux.cxx.

45 {
46  return TMath::Min(fMaxEv, fMaxEvCut);
47 }
double fMaxEvCut
user-defined cut: maximum energy
Definition: GAtmoFlux.h:138
double fMaxEv
maximum energy (in input flux files)
Definition: GAtmoFlux.h:131
double genie::flux::GAtmoFlux::MinEnergy ( void  )
inlineprotected

Definition at line 120 of file GAtmoFlux.h.

References fMinEvCut.

120 { return fMinEvCut; }
double fMinEvCut
user-defined cut: minimum energy
Definition: GAtmoFlux.h:139
virtual const TLorentzVector& genie::flux::GAtmoFlux::Momentum ( void  )
inlinevirtual

returns the flux neutrino 4-momentum

Implements genie::GFluxI.

Definition at line 71 of file GAtmoFlux.h.

References fgP4.

71 { return fgP4; }
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:134
long int GAtmoFlux::NFluxNeutrinos ( void  ) const

Number of flux nu's generated. Not the same as the number of nu's thrown towards the geometry (if there are cuts).

Definition at line 222 of file GAtmoFlux.cxx.

223 {
224  return fNNeutrinos;
225 }
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GAtmoFlux.h:137
virtual int genie::flux::GAtmoFlux::PdgCode ( void  )
inlinevirtual

returns the flux neutrino pdg code

Implements genie::GFluxI.

Definition at line 69 of file GAtmoFlux.h.

References fgPdgC.

69 { return fgPdgC; }
int fgPdgC
current generated nu pdg-code
Definition: GAtmoFlux.h:133
virtual const TLorentzVector& genie::flux::GAtmoFlux::Position ( void  )
inlinevirtual

returns the flux neutrino 4-position (note: expect SI rather than physical units)

Implements genie::GFluxI.

Definition at line 72 of file GAtmoFlux.h.

References fgX4.

72 { return fgX4; }
TLorentzVector fgX4
current generated nu 4-position
Definition: GAtmoFlux.h:135
void GAtmoFlux::ResetSelection ( void  )
protected

Definition at line 326 of file GAtmoFlux.cxx.

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 }
int fgPdgC
current generated nu pdg-code
Definition: GAtmoFlux.h:133
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:134
TLorentzVector fgX4
current generated nu 4-position
Definition: GAtmoFlux.h:135
double fWeight
current generated nu weight
Definition: GAtmoFlux.h:136
int GAtmoFlux::SelectNeutrino ( double  Ev,
double  costheta,
double  phi 
)
protected

Definition at line 580 of file GAtmoFlux.cxx.

References genie::RandomGen::Instance(), LOG, pDEBUG, pERROR, pINFO, and genie::RandomGen::RndFlux().

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 }
map< int, TH3D * > fFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:154
#define pERROR
Definition: Messenger.h:59
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pINFO
Definition: Messenger.h:62
PDGCodeList * fPdgCList
input list of neutrino pdg-codes
Definition: GAtmoFlux.h:132
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:71
#define pDEBUG
Definition: Messenger.h:63
void GAtmoFlux::SetRadii ( double  Rlongitudinal,
double  Rtransverse 
)

Definition at line 369 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

Referenced by GetFlux().

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 }
double fRl
defining flux neutrino generation surface: longitudinal radius
Definition: GAtmoFlux.h:140
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double fRt
defining flux neutrino generation surface: transverse radius
Definition: GAtmoFlux.h:141
#define pNOTICE
Definition: Messenger.h:61
void GAtmoFlux::SetSpectralIndex ( double  index)

Definition at line 251 of file GAtmoFlux.cxx.

References LOG, pNOTICE, and pWARN.

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 }
double fSpectralIndex
power law function used for weighted flux
Definition: GAtmoFlux.h:150
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pWARN
Definition: Messenger.h:60
#define pNOTICE
Definition: Messenger.h:61
void GAtmoFlux::SetUserCoordSystem ( TRotation &  rotation)

Rotation: Topocentric Horizontal -> User-defined Topocentric Coord System.

Definition at line 263 of file GAtmoFlux.cxx.

Referenced by GetFlux().

264 {
265  fRotTHz2User = rotation;
266 }
TRotation fRotTHz2User
coord. system rotation: THZ -&gt; Topocentric user-defined
Definition: GAtmoFlux.h:142
virtual double genie::flux::GAtmoFlux::Weight ( void  )
inlinevirtual

returns the flux neutrino weight (if any)

Implements genie::GFluxI.

Definition at line 70 of file GAtmoFlux.h.

References fWeight.

70 { return fWeight; }
double fWeight
current generated nu weight
Definition: GAtmoFlux.h:136
void GAtmoFlux::ZeroFluxHisto ( TH3D *  hist)
protected

Definition at line 544 of file GAtmoFlux.cxx.

References LOG, and pDEBUG.

545 {
546  LOG("Flux", pDEBUG) << "Forcing flux histogram contents to 0";
547  if(!histo) return;
548  histo->Reset();
549 }
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pDEBUG
Definition: Messenger.h:63

Member Data Documentation

double* genie::flux::GAtmoFlux::fCosThetaBins
protected
double* genie::flux::GAtmoFlux::fEnergyBins
protected
vector<string> genie::flux::GAtmoFlux::fFluxFile
protected

input flux file for each neutrino species

Definition at line 157 of file GAtmoFlux.h.

vector<int> genie::flux::GAtmoFlux::fFluxFlavour
protected

input flux file for each neutrino species

Definition at line 156 of file GAtmoFlux.h.

map<int, TH3D*> genie::flux::GAtmoFlux::fFluxHistoMap
protected

flux = f(Ev,cos8,phi) for each neutrino species

Definition at line 154 of file GAtmoFlux.h.

bool genie::flux::GAtmoFlux::fGenWeighted
protected

generate a weighted or unweighted flux?

Definition at line 149 of file GAtmoFlux.h.

TLorentzVector genie::flux::GAtmoFlux::fgP4
protected

current generated nu 4-momentum

Definition at line 134 of file GAtmoFlux.h.

Referenced by CosTheta(), CosZenith(), Energy(), Enu(), and Momentum().

int genie::flux::GAtmoFlux::fgPdgC
protected

current generated nu pdg-code

Definition at line 133 of file GAtmoFlux.h.

Referenced by PdgCode().

TLorentzVector genie::flux::GAtmoFlux::fgX4
protected

current generated nu 4-position

Definition at line 135 of file GAtmoFlux.h.

Referenced by Position().

bool genie::flux::GAtmoFlux::fInitialized
protected

flag to check that initialization is run

Definition at line 151 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fMaxEv
protected
double genie::flux::GAtmoFlux::fMaxEvCut
protected

user-defined cut: maximum energy

Definition at line 138 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fMinEvCut
protected

user-defined cut: minimum energy

Definition at line 139 of file GAtmoFlux.h.

Referenced by MinEnergy().

long int genie::flux::GAtmoFlux::fNNeutrinos
protected

number of flux neutrinos thrown so far

Definition at line 137 of file GAtmoFlux.h.

unsigned int genie::flux::GAtmoFlux::fNumCosThetaBins
protected

number of cos(theta) bins in input flux data files

Definition at line 144 of file GAtmoFlux.h.

Referenced by genie::flux::GHAKKMAtmoFlux::SetBinSizes(), genie::flux::GFLUKAAtmoFlux::SetBinSizes(), and genie::flux::GBGLRSAtmoFlux::SetBinSizes().

unsigned int genie::flux::GAtmoFlux::fNumEnergyBins
protected

number of energy bins in input flux data files

Definition at line 145 of file GAtmoFlux.h.

Referenced by genie::flux::GHAKKMAtmoFlux::SetBinSizes(), genie::flux::GFLUKAAtmoFlux::SetBinSizes(), and genie::flux::GBGLRSAtmoFlux::SetBinSizes().

unsigned int genie::flux::GAtmoFlux::fNumPhiBins
protected

number of phi bins in input flux data files

Definition at line 143 of file GAtmoFlux.h.

Referenced by genie::flux::GHAKKMAtmoFlux::SetBinSizes(), genie::flux::GFLUKAAtmoFlux::SetBinSizes(), and genie::flux::GBGLRSAtmoFlux::SetBinSizes().

PDGCodeList* genie::flux::GAtmoFlux::fPdgCList
protected

input list of neutrino pdg-codes

Definition at line 132 of file GAtmoFlux.h.

Referenced by FluxParticles().

double* genie::flux::GAtmoFlux::fPhiBins
protected
map<int, TH3D*> genie::flux::GAtmoFlux::fRawFluxHistoMap
protected

flux = f(Ev,cos8,phi) for each neutrino species

Definition at line 155 of file GAtmoFlux.h.

Referenced by genie::flux::GHAKKMAtmoFlux::FillFluxHisto(), genie::flux::GFLUKAAtmoFlux::FillFluxHisto(), and genie::flux::GBGLRSAtmoFlux::FillFluxHisto().

double genie::flux::GAtmoFlux::fRl
protected

defining flux neutrino generation surface: longitudinal radius

Definition at line 140 of file GAtmoFlux.h.

TRotation genie::flux::GAtmoFlux::fRotTHz2User
protected

coord. system rotation: THZ -> Topocentric user-defined

Definition at line 142 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fRt
protected

defining flux neutrino generation surface: transverse radius

Definition at line 141 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fSpectralIndex
protected

power law function used for weighted flux

Definition at line 150 of file GAtmoFlux.h.

TH3D* genie::flux::GAtmoFlux::fTotalFluxHisto
protected

flux = f(Ev,cos8,phi) summed over neutrino species

Definition at line 152 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fTotalFluxHistoIntg
protected

fFluxSum2D integral

Definition at line 153 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fWeight
protected

current generated nu weight

Definition at line 136 of file GAtmoFlux.h.

Referenced by Weight().


The documentation for this class was generated from the following files: