GBGLRSAtmoFlux.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: Christopher Backhouse <c.backhouse1@physics.ox.ac.uk>
8  Oxford University
9 
10  For the class documentation see the corresponding header file.
11  @ Feb 05, 2008 - CB
12  This concrete flux driver was added in 2.3.1 by C.Backhouse (Oxford U.)
13  @ Feb 23, 2010 - CA
14  Build bin arrays at ctor. Re-structuring and clean-up.
15  @ Feb 23, 2012 - AB
16  Combine the flux calculations at low and high energies.
17 
18 */
19 //____________________________________________________________________________
20 
21 #include <fstream>
22 #include <cassert>
23 
24 #include <TH3D.h>
25 #include <TMath.h>
26 
30 
33 
34 using std::ifstream;
35 using std::ios;
36 
37 using namespace genie;
38 using namespace genie::flux;
39 using namespace genie::constants;
40 
41 //____________________________________________________________________________
43 GAtmoFlux()
44 {
45  LOG("Flux", pNOTICE)
46  << "Instantiating the GENIE BGLRS atmospheric neutrino flux driver";
47 
48  this->Initialize();
49  this->SetBinSizes();
50 }
51 //___________________________________________________________________________
52 GBGLRSAtmoFlux::~GBGLRSAtmoFlux()
53 {
54 
55 }
56 //___________________________________________________________________________
58 {
59 // Generate the correct cos(theta) and energy bin sizes.
60 //
61 // Zenith angle binning: the flux is given in 20 bins of
62 // cos(zenith angle) from -1.0 to 1.0 (bin width = 0.1)
63 //
64 // Neutrino energy binning: the Bartol flux files are
65 // provided in two pieces
66 // (1) low energy piece (<10 GeV), solar min or max,
67 // given in 40 log-spaced bins from 0.1 to 10 GeV
68 // (20 bins per decade)
69 // (2) high energy piece (>10 GeV), without solar effects,
70 // given in 30 log-spaced bins from 10 to 1000 GeV
71 // (10 bins per decade)
72 
73  fPhiBins = new double [2];
74  fCosThetaBins = new double [kBGLRS3DNumCosThetaBins + 1];
76 
77  fPhiBins[0] = 0;
78  fPhiBins[1] = 2.*kPi;
79 
80  double dcostheta =
82  (double) kBGLRS3DNumCosThetaBins;
83 
84  double logEmin = TMath::Log10(kBGLRS3DEvMin);
85  double dlogElow = 1.0 / (double) kBGLRS3DNumLogEvBinsPerDecadeLow;
86  double dlogEhigh = 1.0 / (double) kBGLRS3DNumLogEvBinsPerDecadeHigh;
87 
88  double costheta = kBGLRS3DCosThetaMin;
89 
90  for(unsigned int i=0; i<= kBGLRS3DNumCosThetaBins; i++) {
91  if( i==0 ) ; // do nothing
92  else costheta += dcostheta;
93  fCosThetaBins[i] = costheta;
94  if(i != kBGLRS3DNumCosThetaBins) {
95  LOG("Flux", pDEBUG)
96  << "FLUKA 3d flux: CosTheta bin " << i+1
97  << ": lower edge = " << fCosThetaBins[i];
98  } else {
99  LOG("Flux", pDEBUG)
100  << "FLUKA 3d flux: CosTheta bin " << kBGLRS3DNumCosThetaBins
101  << ": upper edge = " << fCosThetaBins[kBGLRS3DNumCosThetaBins];
102  }
103  }
104 
105  double logE = logEmin;
106 
107  for(unsigned int i=0; i<=kBGLRS3DNumLogEvBinsLow+kBGLRS3DNumLogEvBinsHigh; i++) {
108  if( i==0 ) ; // do nothing
109  else if( i<=kBGLRS3DNumLogEvBinsLow ) logE += dlogElow;
110  else logE += dlogEhigh;
111  fEnergyBins[i] = TMath::Power(10.0, logE);
112  if(i != kBGLRS3DNumLogEvBinsLow+kBGLRS3DNumLogEvBinsHigh) {
113  LOG("Flux", pDEBUG)
114  << "FLUKA 3d flux: Energy bin " << i+1
115  << ": lower edge = " << fEnergyBins[i];
116  } else {
117  LOG("Flux", pDEBUG)
118  << "FLUKA 3d flux: Energy bin " << kBGLRS3DNumLogEvBinsLow+kBGLRS3DNumLogEvBinsHigh
120  }
121  }
122 
123  fNumPhiBins = 1;
127 }
128 //___________________________________________________________________________
129 bool GBGLRSAtmoFlux::FillFluxHisto(int nu_pdg, string filename)
130 {
131  LOG("Flux", pNOTICE)
132  << "Loading BGLRS flux for neutrino: " << nu_pdg
133  << " from file: " << filename;
134 
135  TH3D* histo = 0;
136  std::map<int,TH3D*>::iterator myMapEntry = fRawFluxHistoMap.find(nu_pdg);
137  if( myMapEntry != fRawFluxHistoMap.end() ){
138  histo = myMapEntry->second;
139  }
140  if(!histo) {
141  LOG("Flux", pERROR) << "Null flux histogram!";
142  return false;
143  }
144  ifstream flux_stream(filename.c_str(), ios::in);
145  if(!flux_stream) {
146  LOG("Flux", pERROR) << "Could not open file: " << filename;
147  return false;
148  }
149 
150  int ibin;
151  double energy, costheta, flux;
152  double junkd; // throw away error estimates
153 
154  double scale = 1.0; // 1.0 [m^2], OR 1.0e-4 [cm^2]
155 
156  // throw away comment line
157  flux_stream.ignore(99999, '\n');
158 
159  while ( !flux_stream.eof() ) {
160  flux = 0.0;
161  flux_stream >> energy >> costheta >> flux >> junkd >> junkd;
162  if( flux>0.0 ){
163  // Compensate for logarithmic units - dlogE=dE/E
164  // [Note: should do this explicitly using bin widths]
165  flux /= energy;
166  LOG("Flux", pINFO)
167  << "Flux[Ev = " << energy
168  << ", cos8 = " << costheta << "] = " << flux;
169  ibin = histo->FindBin( (Axis_t)energy, (Axis_t)costheta, (Axis_t)kPi );
170  histo->SetBinContent( ibin, (Stat_t)(scale*flux) );
171  }
172  }
173  return true;
174 }
175 //___________________________________________________________________________
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
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:136
const unsigned int kBGLRS3DNumLogEvBinsHigh
string filename
Definition: shutoffs.py:106
const unsigned int kBGLRS3DNumLogEvBinsLow
void Initialize(Bool_t useTMVAStyle=kTRUE)
Definition: tmvaglob.C:184
Loaders::FluxType flux
const unsigned int kBGLRS3DNumLogEvBinsPerDecadeLow
const double kBGLRS3DEvMin
Double_t scale
Definition: plot.C:25
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const double kBGLRS3DCosThetaMax
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
const unsigned int kBGLRS3DNumLogEvBinsPerDecadeHigh
#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
const unsigned int kBGLRS3DNumCosThetaBins
A flux driver for the Bartol Atmospheric Neutrino Flux.
#define FLUXDRIVERREG4(_nsa, _nsb, _name, _fqname)
#define pDEBUG
Definition: Messenger.h:64