GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GFLUKAAtmoFlux.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 
20 
23 
24 using std::ifstream;
25 using std::ios;
26 using namespace genie;
27 using namespace genie::flux;
28 using namespace genie::constants;
29 
30 //____________________________________________________________________________
32 GAtmoFlux()
33 {
34  LOG("Flux", pNOTICE)
35  << "Instantiating the GENIE FLUKA atmospheric neutrino flux driver";
36 
37  this->Initialize();
38  this->SetBinSizes();
39 }
40 //___________________________________________________________________________
41 GFLUKAAtmoFlux::~GFLUKAAtmoFlux()
42 {
43 
44 }
45 //___________________________________________________________________________
47 {
48 // Generate the correct cos(theta) and energy bin sizes
49 // The flux is given in 40 bins of cos(zenith angle) from -1.0 to 1.0
50 // (bin width = 0.05) and 61 equally log-spaced energy bins (20 bins
51 // per decade), with Emin = 0.100 GeV.
52 //
53 
54  fPhiBins = new double [2];
55  fCosThetaBins = new double [kGFlk3DNumCosThetaBins + 1];
56  fEnergyBins = new double [kGFlk3DNumLogEvBins + 1];
57 
58  fPhiBins[0] = 0;
59  fPhiBins[1] = 2.*kPi;
60 
61  double dcostheta =
63  (double) kGFlk3DNumCosThetaBins;
64 
65  double logEmax = TMath::Log10(1.);
66  double logEmin = TMath::Log10(kGFlk3DEvMin);
67  double dlogE =
68  (logEmax - logEmin) /
70 
71  for(unsigned int i=0; i<= kGFlk3DNumCosThetaBins; i++) {
72  fCosThetaBins[i] = kGFlk3DCosThetaMin + i * dcostheta;
73  if(i != kGFlk3DNumCosThetaBins) {
74  LOG("Flux", pDEBUG)
75  << "FLUKA 3d flux: CosTheta bin " << i+1
76  << ": lower edge = " << fCosThetaBins[i];
77  } else {
78  LOG("Flux", pDEBUG)
79  << "FLUKA 3d flux: CosTheta bin " << kGFlk3DNumCosThetaBins
80  << ": upper edge = " << fCosThetaBins[kGFlk3DNumCosThetaBins];
81  }
82  }
83 
84  for(unsigned int i=0; i<= kGFlk3DNumLogEvBins; i++) {
85  fEnergyBins[i] = TMath::Power(10., logEmin + i*dlogE);
86  if(i != kGFlk3DNumLogEvBins) {
87  LOG("Flux", pDEBUG)
88  << "FLUKA 3d flux: Energy bin " << i+1
89  << ": lower edge = " << fEnergyBins[i];
90  } else {
91  LOG("Flux", pDEBUG)
92  << "FLUKA 3d flux: Energy bin " << kGFlk3DNumLogEvBins
93  << ": upper edge = " << fEnergyBins[kGFlk3DNumLogEvBins];
94  }
95  }
96 
97  for(unsigned int i=0; i< kGFlk3DNumLogEvBins; i++) {
98  LOG("Flux", pDEBUG)
99  << "FLUKA 3d flux: Energy bin " << i+1
100  << ": bin centre = " << (fEnergyBins[i] + fEnergyBins[i+1])/2.;
101  }
102 
103  fNumPhiBins = 1;
107 }
108 //____________________________________________________________________________
109 bool GFLUKAAtmoFlux::FillFluxHisto(int nu_pdg, string filename)
110 {
111  LOG("Flux", pNOTICE)
112  << "Loading FLUKA atmospheric flux for neutrino: " << nu_pdg
113  << " from file: " << filename;
114 
115  TH3D* histo = 0;
116  std::map<int,TH3D*>::iterator myMapEntry = fRawFluxHistoMap.find(nu_pdg);
117  if( myMapEntry != fRawFluxHistoMap.end() ){
118  histo = myMapEntry->second;
119  }
120  if(!histo) {
121  LOG("Flux", pERROR) << "Null flux histogram!";
122  return false;
123  }
124  ifstream flux_stream(filename.c_str(), ios::in);
125  if(!flux_stream) {
126  LOG("Flux", pERROR) << "Could not open file: " << filename;
127  return false;
128  }
129 
130  int ibin;
131  double energy, costheta, flux;
132  char j1, j2;
133 
134  double scale = 1.0; // 1.0 [m^2], OR 1.0e-4 [cm^2]
135 
136  while ( !flux_stream.eof() ) {
137  flux = 0.0;
138  flux_stream >> energy >> j1 >> costheta >> j2 >> flux;
139  if( flux>0.0 ){
140  LOG("Flux", pINFO)
141  << "Flux[Ev = " << energy
142  << ", cos8 = " << costheta << "] = " << flux;
143  // note: reversing the Fluka sign convention for zenith angle
144  // 1 phi bin
145  ibin = histo->FindBin( (Axis_t)energy, (Axis_t)(-costheta), (Axis_t)kPi );
146  histo->SetBinContent( ibin, (Stat_t)(scale*flux) );
147  }
148  }
149  return true;
150 }
151 //___________________________________________________________________________
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
unsigned int fNumCosThetaBins
number of cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:144
const double kGFlk3DCosThetaMax
const unsigned int kGFlk3DNumCosThetaBins
const unsigned int kGFlk3DNumLogEvBins
const unsigned int kGFlk3DNumLogEvBinsPerDecade
bool FillFluxHisto(int nu_pdg, string filename)
A flux driver for the FLUKA 3-D Atmospheric Neutrino Flux.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const double kGFlk3DCosThetaMin
const double kGFlk3DEvMin
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
#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
#define FLUXDRIVERREG4(_nsa, _nsb, _name, _fqname)
#define pDEBUG
Definition: Messenger.h:63