GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Private Member Functions | List of all members
genie::flux::GHAKKMAtmoFlux Class Reference

A driver for the HAKKM 3-D atmospheric neutrino flux (commonly known as the `Honda flux') More...

#include <GHAKKMAtmoFlux.h>

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

Public Member Functions

 GHAKKMAtmoFlux ()
 
 ~GHAKKMAtmoFlux ()
 
- Public Member Functions inherited from genie::flux::GAtmoFlux
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 ()
 

Private Member Functions

void SetBinSizes (void)
 
bool FillFluxHisto (int nu_pdg, string filename)
 

Additional Inherited Members

- Protected Member Functions inherited from genie::flux::GAtmoFlux
 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)
 
- Protected Member Functions inherited from genie::GFluxI
 GFluxI ()
 
- Protected Attributes inherited from genie::flux::GAtmoFlux
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 driver for the HAKKM 3-D atmospheric neutrino flux (commonly known as the `Honda flux')

References:
M. Honda, M. Sajjad Athar, T. Kajita, K. Kasahara, and S. Midorikawa Phys. Rev. D 92 (2015) 023004

The flux files necessary for running this flux driver can be obtained from:​http://www.icrr.u-tokyo.ac.jp/~mhonda/

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

Definition at line 45 of file GHAKKMAtmoFlux.h.

Constructor & Destructor Documentation

GHAKKMAtmoFlux::GHAKKMAtmoFlux ( )

Definition at line 32 of file GHAKKMAtmoFlux.cxx.

References Initialize(), LOG, and pNOTICE.

32  :
33 GAtmoFlux()
34 {
35  LOG("Flux", pNOTICE)
36  << "Instantiating the GENIE HAKKM atmospheric neutrino flux driver";
37 
38  this->Initialize();
39  this->SetBinSizes();
40 }
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pNOTICE
Definition: Messenger.h:61
GHAKKMAtmoFlux::~GHAKKMAtmoFlux ( )

Definition at line 42 of file GHAKKMAtmoFlux.cxx.

43 {
44 
45 }

Member Function Documentation

bool GHAKKMAtmoFlux::FillFluxHisto ( int  nu_pdg,
string  filename 
)
privatevirtual

Implements genie::flux::GAtmoFlux.

Definition at line 135 of file GHAKKMAtmoFlux.cxx.

References genie::flux::GAtmoFlux::fCosThetaBins, genie::flux::GAtmoFlux::fEnergyBins, genie::flux::GAtmoFlux::fPhiBins, genie::flux::GAtmoFlux::fRawFluxHistoMap, genie::flux::kGHnd3DNumCosThetaBins, genie::flux::kGHnd3DNumLogEvBins, genie::flux::kGHnd3DNumPhiBins, genie::kPdgAntiNuE, genie::kPdgAntiNuMu, genie::kPdgNuE, genie::kPdgNuMu, genie::constants::kPi, LOG, pDEBUG, pERROR, and pNOTICE.

136 {
137  LOG("Flux", pNOTICE)
138  << "Loading HAKKM flux for neutrino: " << nu_pdg
139  << " from file: " << filename;
140 
141  TH3D* histo = 0;
142  std::map<int,TH3D*>::iterator myMapEntry = fRawFluxHistoMap.find(nu_pdg);
143  if( myMapEntry != fRawFluxHistoMap.end() ){
144  histo = myMapEntry->second;
145  }
146  if(!histo) {
147  LOG("Flux", pERROR) << "Null flux histogram!";
148  return false;
149  }
150 
151  ifstream flux_stream(filename.c_str(), ios::in);
152  if(!flux_stream) {
153  LOG("Flux", pERROR) << "Could not open file: " << filename;
154  return false;
155  }
156 
157  double r2d = 180./constants::kPi;
158 
159  int icostheta = kGHnd3DNumCosThetaBins; // starting from last bin [0.9 - 1.0]
160  int iphi = 0;
161 
162  while (!flux_stream.eof()) {
163 
164  string comment = "";
165  getline(flux_stream,comment);
166  LOG("Flux", pDEBUG) << "Comment line from HAKKM input file: " << comment;
167  getline(flux_stream,comment);
168  LOG("Flux", pDEBUG) << "Comment line from HAKKM input file: " << comment;
169 
170  iphi++;
171  if(iphi == kGHnd3DNumPhiBins+1) {
172  icostheta--;
173  iphi=1;
174  }
175 
176  LOG("Flux", pDEBUG)
177  << "icostheta = " << icostheta << ", iphi = " << iphi << " / "
178  << "costheta bins = " << kGHnd3DNumCosThetaBins << ", phi bins = " << kGHnd3DNumPhiBins;
179  LOG("Flux", pDEBUG)
180  << "The following set of (energy,flux) values corresponds to "
181  << "costheta = [" << fCosThetaBins[icostheta-1] << ", " << fCosThetaBins[icostheta] << "]"
182  << ", phi = [" << fPhiBins[iphi-1] << ", " << fPhiBins[iphi] << "] rad"
183  << " or [" << r2d * fPhiBins[iphi-1] << ", " << r2d * fPhiBins[iphi] << "] deg) ";
184 
185  for(unsigned int i=0; i < kGHnd3DNumLogEvBins; i++) {
186 
187  int ienergy = i+1;
188 
189  double energy = 0;
190  double fnumu = 0;
191  double fnumubar = 0;
192  double fnue = 0;
193  double fnuebar = 0;
194 
195  flux_stream >> energy >> fnumu >> fnumubar >> fnue >> fnuebar;
196 
197  // fitting this easily into what is done for FLUKA, BGLRS where a
198  // different file is specified for each neurtino species means that
199  // the input file for HAKKM has to be read 4 times (at most).
200  // However, this maintains the ability to switch off individual
201  // components at source and generate interactions for some species only
202 
203  double flux = 0.;
204  if(nu_pdg == kPdgNuMu ) flux = fnumu;
205  if(nu_pdg == kPdgAntiNuMu) flux = fnumubar;
206  if(nu_pdg == kPdgNuE ) flux = fnue;
207  if(nu_pdg == kPdgAntiNuE ) flux = fnuebar;
208  LOG("Flux", pDEBUG)
209  << "Flux (nu_pdg = " << nu_pdg
210  << "; Ev = " << energy << " GeV / bin used = ["
211  << fEnergyBins[ienergy-1] << ", " << fEnergyBins[ienergy] << "] GeV"
212  << ") = " << flux << " (m^2 sec sr GeV)^-1";
213  if(flux > 0.) {
214  histo->SetBinContent(ienergy,icostheta,iphi,flux);
215  }
216  }
217  getline(flux_stream,comment);
218  LOG("Flux", pDEBUG) << comment;
219  }
220  return true;
221 }
double * fPhiBins
phi bins in input flux data files
Definition: GAtmoFlux.h:146
const int kPdgNuE
Definition: PDGCodes.h:28
#define pERROR
Definition: Messenger.h:59
const unsigned int kGHnd3DNumPhiBins
const int kPdgAntiNuE
Definition: PDGCodes.h:29
const int kPdgNuMu
Definition: PDGCodes.h:30
const unsigned int kGHnd3DNumLogEvBins
#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
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
map< int, TH3D * > fRawFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:155
const unsigned int kGHnd3DNumCosThetaBins
#define pNOTICE
Definition: Messenger.h:61
#define pDEBUG
Definition: Messenger.h:63
void GHAKKMAtmoFlux::SetBinSizes ( void  )
private

Definition at line 47 of file GHAKKMAtmoFlux.cxx.

References genie::flux::GAtmoFlux::fCosThetaBins, genie::flux::GAtmoFlux::fEnergyBins, genie::flux::GAtmoFlux::fMaxEv, genie::flux::GAtmoFlux::fNumCosThetaBins, genie::flux::GAtmoFlux::fNumEnergyBins, genie::flux::GAtmoFlux::fNumPhiBins, genie::flux::GAtmoFlux::fPhiBins, genie::flux::kGHnd3DCosThetaMax, genie::flux::kGHnd3DCosThetaMin, genie::flux::kGHnd3DNumCosThetaBins, genie::flux::kGHnd3DNumLogEvBins, genie::flux::kGHnd3DNumLogEvBinsPerDecade, genie::flux::kGHnd3DNumPhiBins, genie::flux::kGHnd3DPhiMax, genie::flux::kGHnd3DPhiMin, genie::constants::kPi, LOG, and pDEBUG.

48 {
49  //
50  // cos(theta)
51  //
52 
53  fCosThetaBins = new double [kGHnd3DNumCosThetaBins + 1];
55 
56  double dcostheta =
58  (double) kGHnd3DNumCosThetaBins;
59 
60  for(unsigned int i=0; i<= kGHnd3DNumCosThetaBins; i++) {
61  fCosThetaBins[i] = kGHnd3DCosThetaMin + i * dcostheta;
62  if(i != kGHnd3DNumCosThetaBins) {
63  LOG("Flux", pDEBUG)
64  << "HAKKM 3D flux: CosTheta bin " << i+1
65  << ": lower edge = " << fCosThetaBins[i];
66  } else {
67  LOG("Flux", pDEBUG)
68  << "HAKKM 3D flux: CosTheta bin " << kGHnd3DNumCosThetaBins
69  << ": upper edge = " << fCosThetaBins[kGHnd3DNumCosThetaBins];
70  }
71  }
72 
73  //
74  // phi
75  //
76 
77  fPhiBins = new double [kGHnd3DNumPhiBins + 1];
79 
80  double d2r = constants::kPi/180.;
81 
82  double dphi =
83  d2r * (kGHnd3DPhiMax - kGHnd3DPhiMin) /
84  (double) kGHnd3DNumPhiBins;
85 
86  for(unsigned int i=0; i<= kGHnd3DNumPhiBins; i++) {
87  fPhiBins[i] = kGHnd3DPhiMin + i * dphi;
88  if(i != kGHnd3DNumPhiBins) {
89  LOG("Flux", pDEBUG)
90  << "HAKKM 3D flux: Phi bin " << i+1
91  << ": lower edge = " << fPhiBins[i];
92  } else {
93  LOG("Flux", pDEBUG)
94  << "HAKKM 3D flux: Phi bin " << kGHnd3DNumPhiBins
95  << ": upper edge = " << fPhiBins[kGHnd3DNumPhiBins];
96  }
97  }
98 
99  //
100  // log(E)
101  //
102 
103  // For each costheta,phi pair there are N logarithmically spaced
104  // neutrino energy values (starting at 0.1 GeV with 20 values per decade
105  // up to 10000 GeV) each with corresponding flux values.
106  // To construct a flux histogram, use N+1 bins from 0 up to maximum
107  // value. Assume that the flux value given for E=E0 is the flux value
108  // at the bin that has E0 as its upper edge.
109  //
110  fEnergyBins = new double [kGHnd3DNumLogEvBins + 1]; // bin edges
112 
113  double logEmax = TMath::Log10(1.);
114  double logEmin = TMath::Log10(0.1);
115  double dlogE =
116  (logEmax - logEmin) /
118 
119  fEnergyBins[0] = 0;
120  for(unsigned int i=0; i < fNumEnergyBins; i++) {
121  fEnergyBins[i+1] = TMath::Power(10., logEmin + i*dlogE);
122  if(i < kGHnd3DNumLogEvBins) {
123  LOG("Flux", pDEBUG)
124  << "HAKKM 3D flux: Energy bin " << i+1
125  << ": upper edge = " << fEnergyBins[i+1] << " GeV";
126  }
127  }
128 
130  LOG("Flux", pDEBUG)
131  << "HAKKM 3D flux: Maximum energy = " << fMaxEv;
132 
133 }
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
const unsigned int kGHnd3DNumPhiBins
unsigned int fNumCosThetaBins
number of cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:144
const double kGHnd3DCosThetaMin
const unsigned int kGHnd3DNumLogEvBins
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const double kGHnd3DCosThetaMax
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 double kGHnd3DPhiMax
const unsigned int kGHnd3DNumCosThetaBins
double fMaxEv
maximum energy (in input flux files)
Definition: GAtmoFlux.h:131
const double kGHnd3DPhiMin
#define pDEBUG
Definition: Messenger.h:63
const unsigned int kGHnd3DNumLogEvBinsPerDecade

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