GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GHAKKMAtmoFlux.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 <fstream>
13 
14 #include <TH3D.h>
15 #include <TMath.h>
16 
21 
24 
25 using std::ifstream;
26 using std::getline;
27 using std::ios;
28 using namespace genie;
29 using namespace genie::flux;
30 
31 //____________________________________________________________________________
33 GAtmoFlux()
34 {
35  LOG("Flux", pNOTICE)
36  << "Instantiating the GENIE HAKKM atmospheric neutrino flux driver";
37 
38  this->Initialize();
39  this->SetBinSizes();
40 }
41 //___________________________________________________________________________
42 GHAKKMAtmoFlux::~GHAKKMAtmoFlux()
43 {
44 
45 }
46 //___________________________________________________________________________
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 }
134 //____________________________________________________________________________
135 bool GHAKKMAtmoFlux::FillFluxHisto(int nu_pdg, string filename)
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 }
222 //___________________________________________________________________________
double * fPhiBins
phi bins in input flux data files
Definition: GAtmoFlux.h:146
const int kPdgNuE
Definition: PDGCodes.h:28
unsigned int fNumEnergyBins
number of energy bins in input flux data files
Definition: GAtmoFlux.h:145
#define pERROR
Definition: Messenger.h:59
const unsigned int kGHnd3DNumPhiBins
unsigned int fNumCosThetaBins
number of cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:144
const int kPdgAntiNuE
Definition: PDGCodes.h:29
const int kPdgNuMu
Definition: PDGCodes.h:30
bool FillFluxHisto(int nu_pdg, string filename)
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
void Initialize(void)
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 int kPdgAntiNuMu
Definition: PDGCodes.h:31
const double kGHnd3DPhiMax
map< int, TH3D * > fRawFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:155
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:61
A base class for the FLUKA, BGLRS and ATMNC atmo. nu. flux drivers. The driver depends on data files ...
Definition: GAtmoFlux.h:60
double fMaxEv
maximum energy (in input flux files)
Definition: GAtmoFlux.h:131
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:63
const unsigned int kGHnd3DNumLogEvBinsPerDecade