GHAKKMAtmoFlux.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 */
10 //____________________________________________________________________________
11 
12 #include <cassert>
13 #include <fstream>
14 
15 #include <TH3D.h>
16 #include <TMath.h>
17 
22 
25 
26 using std::ifstream;
27 using std::getline;
28 using std::ios;
29 using namespace genie;
30 using namespace genie::flux;
31 
32 //____________________________________________________________________________
34 GAtmoFlux()
35 {
36  LOG("Flux", pNOTICE)
37  << "Instantiating the GENIE HAKKM atmospheric neutrino flux driver";
38 
39  this->Initialize();
40  this->SetBinSizes();
41 }
42 //___________________________________________________________________________
43 GHAKKMAtmoFlux::~GHAKKMAtmoFlux()
44 {
45 
46 }
47 //___________________________________________________________________________
49 {
50  //
51  // cos(theta)
52  //
53 
54  fCosThetaBins = new double [kGHnd3DNumCosThetaBins + 1];
56 
57  double dcostheta =
59  (double) kGHnd3DNumCosThetaBins;
60 
61  for(unsigned int i=0; i<= kGHnd3DNumCosThetaBins; i++) {
62  fCosThetaBins[i] = kGHnd3DCosThetaMin + i * dcostheta;
63  if(i != kGHnd3DNumCosThetaBins) {
64  LOG("Flux", pDEBUG)
65  << "HAKKM 3D flux: CosTheta bin " << i+1
66  << ": lower edge = " << fCosThetaBins[i];
67  } else {
68  LOG("Flux", pDEBUG)
69  << "HAKKM 3D flux: CosTheta bin " << kGHnd3DNumCosThetaBins
70  << ": upper edge = " << fCosThetaBins[kGHnd3DNumCosThetaBins];
71  }
72  }
73 
74  //
75  // phi
76  //
77 
78  fPhiBins = new double [kGHnd3DNumPhiBins + 1];
80 
81  double d2r = constants::kPi/180.;
82 
83  double dphi =
84  d2r * (kGHnd3DPhiMax - kGHnd3DPhiMin) /
85  (double) kGHnd3DNumPhiBins;
86 
87  for(unsigned int i=0; i<= kGHnd3DNumPhiBins; i++) {
88  fPhiBins[i] = kGHnd3DPhiMin + i * dphi;
89  if(i != kGHnd3DNumPhiBins) {
90  LOG("Flux", pDEBUG)
91  << "HAKKM 3D flux: Phi bin " << i+1
92  << ": lower edge = " << fPhiBins[i];
93  } else {
94  LOG("Flux", pDEBUG)
95  << "HAKKM 3D flux: Phi bin " << kGHnd3DNumPhiBins
96  << ": upper edge = " << fPhiBins[kGHnd3DNumPhiBins];
97  }
98  }
99 
100  //
101  // log(E)
102  //
103 
104  // For each costheta,phi pair there are N logarithmically spaced
105  // neutrino energy values (starting at 0.1 GeV with 20 values per decade
106  // up to 10000 GeV) each with corresponding flux values.
107  // To construct a flux histogram, use N+1 bins from 0 up to maximum
108  // value. Assume that the flux value given for E=E0 is the flux value
109  // at the bin that has E0 as its upper edge.
110  //
111  fEnergyBins = new double [kGHnd3DNumLogEvBins + 1]; // bin edges
113 
114  double logEmax = TMath::Log10(1.);
115  double logEmin = TMath::Log10(0.1);
116  double dlogE =
117  (logEmax - logEmin) /
119 
120  fEnergyBins[0] = 0;
121  for(unsigned int i=0; i < fNumEnergyBins; i++) {
122  fEnergyBins[i+1] = TMath::Power(10., logEmin + i*dlogE);
123  if(i < kGHnd3DNumLogEvBins) {
124  LOG("Flux", pDEBUG)
125  << "HAKKM 3D flux: Energy bin " << i+1
126  << ": upper edge = " << fEnergyBins[i+1] << " GeV";
127  }
128  }
129 
131  LOG("Flux", pDEBUG)
132  << "HAKKM 3D flux: Maximum energy = " << fMaxEv;
133 
134 }
135 //____________________________________________________________________________
136 bool GHAKKMAtmoFlux::FillFluxHisto(int nu_pdg, string filename)
137 {
138  LOG("Flux", pNOTICE)
139  << "Loading HAKKM flux for neutrino: " << nu_pdg
140  << " from file: " << filename;
141 
142  TH3D* histo = 0;
143  std::map<int,TH3D*>::iterator myMapEntry = fRawFluxHistoMap.find(nu_pdg);
144  if( myMapEntry != fRawFluxHistoMap.end() ){
145  histo = myMapEntry->second;
146  }
147  if(!histo) {
148  LOG("Flux", pERROR) << "Null flux histogram!";
149  return false;
150  }
151 
152  ifstream flux_stream(filename.c_str(), ios::in);
153  if(!flux_stream) {
154  LOG("Flux", pERROR) << "Could not open file: " << filename;
155  return false;
156  }
157 
158  double r2d = 180./constants::kPi;
159 
160  int icostheta = kGHnd3DNumCosThetaBins; // starting from last bin [0.9 - 1.0]
161  int iphi = 0;
162 
163  while (!flux_stream.eof()) {
164 
165  string comment = "";
166  getline(flux_stream,comment);
167  LOG("Flux", pDEBUG) << "Comment line from HAKKM input file: " << comment;
168  getline(flux_stream,comment);
169  LOG("Flux", pDEBUG) << "Comment line from HAKKM input file: " << comment;
170 
171  iphi++;
172  if(iphi == kGHnd3DNumPhiBins+1) {
173  icostheta--;
174  iphi=1;
175  }
176 
177  LOG("Flux", pDEBUG)
178  << "icostheta = " << icostheta << ", iphi = " << iphi << " / "
179  << "costheta bins = " << kGHnd3DNumCosThetaBins << ", phi bins = " << kGHnd3DNumPhiBins;
180  LOG("Flux", pDEBUG)
181  << "The following set of (energy,flux) values corresponds to "
182  << "costheta = [" << fCosThetaBins[icostheta-1] << ", " << fCosThetaBins[icostheta] << "]"
183  << ", phi = [" << fPhiBins[iphi-1] << ", " << fPhiBins[iphi] << "] rad"
184  << " or [" << r2d * fPhiBins[iphi-1] << ", " << r2d * fPhiBins[iphi] << "] deg) ";
185 
186  for(unsigned int i=0; i < kGHnd3DNumLogEvBins; i++) {
187 
188  int ienergy = i+1;
189 
190  double energy = 0;
191  double fnumu = 0;
192  double fnumubar = 0;
193  double fnue = 0;
194  double fnuebar = 0;
195 
196  flux_stream >> energy >> fnumu >> fnumubar >> fnue >> fnuebar;
197 
198  // fitting this easily into what is done for FLUKA, BGLRS where a
199  // different file is specified for each neurtino species means that
200  // the input file for HAKKM has to be read 4 times (at most).
201  // However, this maintains the ability to switch off individual
202  // components at source and generate interactions for some species only
203 
204  double flux = 0.;
205  if(nu_pdg == kPdgNuMu ) flux = fnumu;
206  if(nu_pdg == kPdgAntiNuMu) flux = fnumubar;
207  if(nu_pdg == kPdgNuE ) flux = fnue;
208  if(nu_pdg == kPdgAntiNuE ) flux = fnuebar;
209  LOG("Flux", pDEBUG)
210  << "Flux (nu_pdg = " << nu_pdg
211  << "; Ev = " << energy << " GeV / bin used = ["
212  << fEnergyBins[ienergy-1] << ", " << fEnergyBins[ienergy] << "] GeV"
213  << ") = " << flux << " (m^2 sec sr GeV)^-1";
214  if(flux > 0.) {
215  histo->SetBinContent(ienergy,icostheta,iphi,flux);
216  }
217  }
218  getline(flux_stream,comment);
219  LOG("Flux", pDEBUG) << comment;
220  }
221  return true;
222 }
223 //___________________________________________________________________________
A class for generating concrete GFluxI derived classes based on the factory pattern. This code supplies a CPP macro which allows the classes to self-register and thus no modification of this class is needed in order to expand the list of classes it knows about.
double * fPhiBins
phi bins in input flux data files
Definition: GAtmoFlux.h:138
const int kPdgNuE
Definition: PDGCodes.h:28
unsigned int fNumEnergyBins
number of energy bins in input flux data files
Definition: GAtmoFlux.h:137
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
const unsigned int kGHnd3DNumPhiBins
unsigned int fNumCosThetaBins
number of cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:136
const int kPdgAntiNuE
Definition: PDGCodes.h:29
const int kPdgNuMu
Definition: PDGCodes.h:30
bool FillFluxHisto(int nu_pdg, string filename)
string filename
Definition: shutoffs.py:106
const double kGHnd3DCosThetaMin
void Initialize(Bool_t useTMVAStyle=kTRUE)
Definition: tmvaglob.C:184
Loaders::FluxType flux
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:97
const double kGHnd3DCosThetaMax
double energy
Definition: plottest35.C:25
double * fCosThetaBins
cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:139
TH2D * histo
unsigned int fNumPhiBins
number of phi bins in input flux data files
Definition: GAtmoFlux.h:135
double * fEnergyBins
energy bins in input flux data files
Definition: GAtmoFlux.h:140
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
const double kGHnd3DPhiMax
ifstream in
Definition: comparison.C:7
map< int, TH3D * > fRawFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:147
A driver for the HAKKM 3-D atmospheric neutrino flux (commonly known as the `Honda flux&#39;) ...
const unsigned int kGHnd3DNumCosThetaBins
#define pNOTICE
Definition: Messenger.h:62
A base class for the FLUKA, BGLRS and ATMNC atmo. nu. flux drivers. The driver depends on data files ...
Definition: GAtmoFlux.h:61
double fMaxEv
maximum energy (in input flux files)
Definition: GAtmoFlux.h:123
static const double kPi
Definition: Constants.h:38
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
#define FLUXDRIVERREG4(_nsa, _nsb, _name, _fqname)
const double kGHnd3DPhiMin
#define pDEBUG
Definition: Messenger.h:64
const unsigned int kGHnd3DNumLogEvBinsPerDecade