GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GBGLRSAtmoFlux.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  Christopher Backhouse <c.backhouse1@physics.ox.ac.uk>
7  Oxford University
8 */
9 //____________________________________________________________________________
10 
11 #include <fstream>
12 #include <cassert>
13 
14 #include <TH3D.h>
15 #include <TMath.h>
16 
20 
23 
24 using std::ifstream;
25 using std::ios;
26 
27 using namespace genie;
28 using namespace genie::flux;
29 using namespace genie::constants;
30 
31 //____________________________________________________________________________
33 GAtmoFlux()
34 {
35  LOG("Flux", pNOTICE)
36  << "Instantiating the GENIE BGLRS atmospheric neutrino flux driver";
37 
38  this->Initialize();
39  this->SetBinSizes();
40 }
41 //___________________________________________________________________________
42 GBGLRSAtmoFlux::~GBGLRSAtmoFlux()
43 {
44 
45 }
46 //___________________________________________________________________________
48 {
49 // Generate the correct cos(theta) and energy bin sizes.
50 //
51 // Zenith angle binning: the flux is given in 20 bins of
52 // cos(zenith angle) from -1.0 to 1.0 (bin width = 0.1)
53 //
54 // Neutrino energy binning: the Bartol flux files are
55 // provided in two pieces
56 // (1) low energy piece (<10 GeV), solar min or max,
57 // given in 40 log-spaced bins from 0.1 to 10 GeV
58 // (20 bins per decade)
59 // (2) high energy piece (>10 GeV), without solar effects,
60 // given in 30 log-spaced bins from 10 to 1000 GeV
61 // (10 bins per decade)
62 
63  fPhiBins = new double [2];
64  fCosThetaBins = new double [kBGLRS3DNumCosThetaBins + 1];
66 
67  fPhiBins[0] = 0;
68  fPhiBins[1] = 2.*kPi;
69 
70  double dcostheta =
72  (double) kBGLRS3DNumCosThetaBins;
73 
74  double logEmin = TMath::Log10(kBGLRS3DEvMin);
75  double dlogElow = 1.0 / (double) kBGLRS3DNumLogEvBinsPerDecadeLow;
76  double dlogEhigh = 1.0 / (double) kBGLRS3DNumLogEvBinsPerDecadeHigh;
77 
78  double costheta = kBGLRS3DCosThetaMin;
79 
80  for(unsigned int i=0; i<= kBGLRS3DNumCosThetaBins; i++) {
81  if( i==0 ) ; // do nothing
82  else costheta += dcostheta;
83  fCosThetaBins[i] = costheta;
84  if(i != kBGLRS3DNumCosThetaBins) {
85  LOG("Flux", pDEBUG)
86  << "FLUKA 3d flux: CosTheta bin " << i+1
87  << ": lower edge = " << fCosThetaBins[i];
88  } else {
89  LOG("Flux", pDEBUG)
90  << "FLUKA 3d flux: CosTheta bin " << kBGLRS3DNumCosThetaBins
91  << ": upper edge = " << fCosThetaBins[kBGLRS3DNumCosThetaBins];
92  }
93  }
94 
95  double logE = logEmin;
96 
97  for(unsigned int i=0; i<=kBGLRS3DNumLogEvBinsLow+kBGLRS3DNumLogEvBinsHigh; i++) {
98  if( i==0 ) ; // do nothing
99  else if( i<=kBGLRS3DNumLogEvBinsLow ) logE += dlogElow;
100  else logE += dlogEhigh;
101  fEnergyBins[i] = TMath::Power(10.0, logE);
102  if(i != kBGLRS3DNumLogEvBinsLow+kBGLRS3DNumLogEvBinsHigh) {
103  LOG("Flux", pDEBUG)
104  << "FLUKA 3d flux: Energy bin " << i+1
105  << ": lower edge = " << fEnergyBins[i];
106  } else {
107  LOG("Flux", pDEBUG)
108  << "FLUKA 3d flux: Energy bin " << kBGLRS3DNumLogEvBinsLow+kBGLRS3DNumLogEvBinsHigh
110  }
111  }
112 
113  fNumPhiBins = 1;
117 }
118 //___________________________________________________________________________
119 bool GBGLRSAtmoFlux::FillFluxHisto(int nu_pdg, string filename)
120 {
121  LOG("Flux", pNOTICE)
122  << "Loading BGLRS flux for neutrino: " << nu_pdg
123  << " from file: " << filename;
124 
125  TH3D* histo = 0;
126  std::map<int,TH3D*>::iterator myMapEntry = fRawFluxHistoMap.find(nu_pdg);
127  if( myMapEntry != fRawFluxHistoMap.end() ){
128  histo = myMapEntry->second;
129  }
130  if(!histo) {
131  LOG("Flux", pERROR) << "Null flux histogram!";
132  return false;
133  }
134  ifstream flux_stream(filename.c_str(), ios::in);
135  if(!flux_stream) {
136  LOG("Flux", pERROR) << "Could not open file: " << filename;
137  return false;
138  }
139 
140  int ibin;
141  double energy, costheta, flux;
142  double junkd; // throw away error estimates
143 
144  double scale = 1.0; // 1.0 [m^2], OR 1.0e-4 [cm^2]
145 
146  // throw away comment line
147  flux_stream.ignore(99999, '\n');
148 
149  while ( !flux_stream.eof() ) {
150  flux = 0.0;
151  flux_stream >> energy >> costheta >> flux >> junkd >> junkd;
152  if( flux>0.0 ){
153  // Compensate for logarithmic units - dlogE=dE/E
154  // [Note: should do this explicitly using bin widths]
155  flux /= energy;
156  LOG("Flux", pINFO)
157  << "Flux[Ev = " << energy
158  << ", cos8 = " << costheta << "] = " << flux;
159  ibin = histo->FindBin( (Axis_t)energy, (Axis_t)costheta, (Axis_t)kPi );
160  histo->SetBinContent( ibin, (Stat_t)(scale*flux) );
161  }
162  }
163  return true;
164 }
165 //___________________________________________________________________________
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
#define pERROR
Definition: Messenger.h:59
bool FillFluxHisto(int nu_pdg, string filename)
const double kBGLRS3DCosThetaMin
unsigned int fNumCosThetaBins
number of cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:144
const unsigned int kBGLRS3DNumLogEvBinsHigh
const unsigned int kBGLRS3DNumLogEvBinsLow
const unsigned int kBGLRS3DNumLogEvBinsPerDecadeLow
const double kBGLRS3DEvMin
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const double kBGLRS3DCosThetaMax
void Initialize(void)
double * fCosThetaBins
cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:147
#define pINFO
Definition: Messenger.h:62
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
map< int, TH3D * > fRawFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:155
const unsigned int kBGLRS3DNumLogEvBinsPerDecadeHigh
#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
const unsigned int kBGLRS3DNumCosThetaBins
A flux driver for the Bartol Atmospheric Neutrino Flux.
#define FLUXDRIVERREG4(_nsa, _nsb, _name, _fqname)
#define pDEBUG
Definition: Messenger.h:63