GFLUKAAtmoFlux.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  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Feb 05, 2008 - CA
14  In vrs 2.3.1, most of the driver implementation code was factored out
15  to the new GAtmoFlux base class so as to be shared by the functionally
16  identical GBGLRSAtmoFlux driver
17  @ Feb 05, 2008 - Chris Backhouse (Oxford)
18  Fixed a bug in bin definitions (the TH2 constructor takes an array of
19  bin lower edges, but the bin centres were being passed in instead).
20  @ Feb 23, 2010 - CA
21  Build bin arrays at ctor. Re-structuring and clean-up.
22 
23 */
24 //____________________________________________________________________________
25 
26 #include <cassert>
27 #include <fstream>
28 
29 #include <TH3D.h>
30 #include <TMath.h>
31 
35 
38 
39 using std::ifstream;
40 using std::ios;
41 using namespace genie;
42 using namespace genie::flux;
43 using namespace genie::constants;
44 
45 //____________________________________________________________________________
47 GAtmoFlux()
48 {
49  LOG("Flux", pNOTICE)
50  << "Instantiating the GENIE FLUKA atmospheric neutrino flux driver";
51 
52  this->Initialize();
53  this->SetBinSizes();
54 }
55 //___________________________________________________________________________
56 GFLUKAAtmoFlux::~GFLUKAAtmoFlux()
57 {
58 
59 }
60 //___________________________________________________________________________
62 {
63 // Generate the correct cos(theta) and energy bin sizes
64 // The flux is given in 40 bins of cos(zenith angle) from -1.0 to 1.0
65 // (bin width = 0.05) and 61 equally log-spaced energy bins (20 bins
66 // per decade), with Emin = 0.100 GeV.
67 //
68 
69  fPhiBins = new double [2];
70  fCosThetaBins = new double [kGFlk3DNumCosThetaBins + 1];
71  fEnergyBins = new double [kGFlk3DNumLogEvBins + 1];
72 
73  fPhiBins[0] = 0;
74  fPhiBins[1] = 2.*kPi;
75 
76  double dcostheta =
78  (double) kGFlk3DNumCosThetaBins;
79 
80  double logEmax = TMath::Log10(1.);
81  double logEmin = TMath::Log10(kGFlk3DEvMin);
82  double dlogE =
83  (logEmax - logEmin) /
85 
86  for(unsigned int i=0; i<= kGFlk3DNumCosThetaBins; i++) {
87  fCosThetaBins[i] = kGFlk3DCosThetaMin + i * dcostheta;
88  if(i != kGFlk3DNumCosThetaBins) {
89  LOG("Flux", pDEBUG)
90  << "FLUKA 3d flux: CosTheta bin " << i+1
91  << ": lower edge = " << fCosThetaBins[i];
92  } else {
93  LOG("Flux", pDEBUG)
94  << "FLUKA 3d flux: CosTheta bin " << kGFlk3DNumCosThetaBins
95  << ": upper edge = " << fCosThetaBins[kGFlk3DNumCosThetaBins];
96  }
97  }
98 
99  for(unsigned int i=0; i<= kGFlk3DNumLogEvBins; i++) {
100  fEnergyBins[i] = TMath::Power(10., logEmin + i*dlogE);
101  if(i != kGFlk3DNumLogEvBins) {
102  LOG("Flux", pDEBUG)
103  << "FLUKA 3d flux: Energy bin " << i+1
104  << ": lower edge = " << fEnergyBins[i];
105  } else {
106  LOG("Flux", pDEBUG)
107  << "FLUKA 3d flux: Energy bin " << kGFlk3DNumLogEvBins
108  << ": upper edge = " << fEnergyBins[kGFlk3DNumLogEvBins];
109  }
110  }
111 
112  for(unsigned int i=0; i< kGFlk3DNumLogEvBins; i++) {
113  LOG("Flux", pDEBUG)
114  << "FLUKA 3d flux: Energy bin " << i+1
115  << ": bin centre = " << (fEnergyBins[i] + fEnergyBins[i+1])/2.;
116  }
117 
118  fNumPhiBins = 1;
122 }
123 //____________________________________________________________________________
124 bool GFLUKAAtmoFlux::FillFluxHisto(int nu_pdg, string filename)
125 {
126  LOG("Flux", pNOTICE)
127  << "Loading FLUKA atmospheric flux for neutrino: " << nu_pdg
128  << " from file: " << filename;
129 
130  TH3D* histo = 0;
131  std::map<int,TH3D*>::iterator myMapEntry = fRawFluxHistoMap.find(nu_pdg);
132  if( myMapEntry != fRawFluxHistoMap.end() ){
133  histo = myMapEntry->second;
134  }
135  if(!histo) {
136  LOG("Flux", pERROR) << "Null flux histogram!";
137  return false;
138  }
139  ifstream flux_stream(filename.c_str(), ios::in);
140  if(!flux_stream) {
141  LOG("Flux", pERROR) << "Could not open file: " << filename;
142  return false;
143  }
144 
145  int ibin;
146  double energy, costheta, flux;
147  char j1, j2;
148 
149  double scale = 1.0; // 1.0 [m^2], OR 1.0e-4 [cm^2]
150 
151  while ( !flux_stream.eof() ) {
152  flux = 0.0;
153  flux_stream >> energy >> j1 >> costheta >> j2 >> flux;
154  if( flux>0.0 ){
155  LOG("Flux", pINFO)
156  << "Flux[Ev = " << energy
157  << ", cos8 = " << costheta << "] = " << flux;
158  // note: reversing the Fluka sign convention for zenith angle
159  // 1 phi bin
160  ibin = histo->FindBin( (Axis_t)energy, (Axis_t)(-costheta), (Axis_t)kPi );
161  histo->SetBinContent( ibin, (Stat_t)(scale*flux) );
162  }
163  }
164  return true;
165 }
166 //___________________________________________________________________________
const double kPi
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
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
unsigned int fNumCosThetaBins
number of cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:136
const double kGFlk3DCosThetaMax
string filename
Definition: shutoffs.py:106
void Initialize(Bool_t useTMVAStyle=kTRUE)
Definition: tmvaglob.C:184
Loaders::FluxType flux
const unsigned int kGFlk3DNumCosThetaBins
const unsigned int kGFlk3DNumLogEvBins
Double_t scale
Definition: plot.C:25
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:97
const double kGFlk3DCosThetaMin
const double kGFlk3DEvMin
double energy
Definition: plottest35.C:25
double * fCosThetaBins
cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:139
TH2D * histo
#define pINFO
Definition: Messenger.h:63
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
ifstream in
Definition: comparison.C:7
map< int, TH3D * > fRawFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:147
#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
#define FLUXDRIVERREG4(_nsa, _nsb, _name, _fqname)
#define pDEBUG
Definition: Messenger.h:64