Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | List of all members
genie::flux::GFLUKAAtmoFlux Class Reference

A flux driver for the FLUKA 3-D Atmospheric Neutrino Flux. More...

#include "/cvmfs/nova.opensciencegrid.org/externals/genie/v3_00_06_p01/Linux64bit+2.6-2.12-e17-debug/GENIE-Generator/src/Tools/Flux/GFLUKAAtmoFlux.h"

Inheritance diagram for genie::flux::GFLUKAAtmoFlux:
genie::flux::GAtmoFlux genie::GFluxI

Public Member Functions

 GFLUKAAtmoFlux ()
 
 ~GFLUKAAtmoFlux ()
 
virtual const PDGCodeListFluxParticles (void)
 declare list of flux neutrinos that can be generated (for init. purposes) More...
 
virtual double MaxEnergy (void)
 declare the max flux neutrino energy that can be generated (for init. purposes) More...
 
virtual bool GenerateNext (void)
 generate the next flux neutrino (return false in err) More...
 
virtual int PdgCode (void)
 returns the flux neutrino pdg code More...
 
virtual double Weight (void)
 returns the flux neutrino weight (if any) More...
 
virtual const TLorentzVector & Momentum (void)
 returns the flux neutrino 4-momentum More...
 
virtual const TLorentzVector & Position (void)
 returns the flux neutrino 4-position (note: expect SI rather than physical units) More...
 
virtual bool End (void)
 true if no more flux nu's can be thrown (eg reaching end of beam sim ntuples) More...
 
virtual long int Index (void)
 returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current entry number) More...
 
virtual void Clear (Option_t *opt)
 reset state variables based on opt More...
 
virtual void GenerateWeighted (bool gen_weighted)
 set whether to generate weighted or unweighted neutrinos More...
 
double Enu (void)
 
double Energy (void)
 
double CosTheta (void)
 
double CosZenith (void)
 
long int NFluxNeutrinos (void) const
 Number of flux nu's generated. Not the same as the number of nu's thrown towards the geometry (if there are cuts). More...
 
void ForceMinEnergy (double emin)
 
void ForceMaxEnergy (double emax)
 
void SetSpectralIndex (double index)
 
void SetRadii (double Rlongitudinal, double Rtransverse)
 
void SetUserCoordSystem (TRotation &rotation)
 Rotation: Topocentric Horizontal -> User-defined Topocentric Coord System. More...
 
void AddFluxFile (int neutrino_pdg, string filename)
 
void AddFluxFile (string filename)
 
bool LoadFluxData (void)
 
TH3D * GetFluxHistogram (int flavour)
 
double GetFlux (int flavour)
 
double GetFlux (int flavour, double energy)
 
double GetFlux (int flavour, double energy, double costh)
 
double GetFlux (int flavour, double energy, double costh, double phi)
 

Protected Member Functions

bool GenerateNext_1try (void)
 
void Initialize (void)
 
void CleanUp (void)
 
void ResetSelection (void)
 
double MinEnergy (void)
 
TH3D * CreateFluxHisto (string name, string title)
 
void ZeroFluxHisto (TH3D *hist)
 
void AddAllFluxes (void)
 
int SelectNeutrino (double Ev, double costheta, double phi)
 
TH3D * CreateNormalisedFluxHisto (TH3D *hist)
 

Protected Attributes

double fMaxEv
 maximum energy (in input flux files) More...
 
PDGCodeListfPdgCList
 input list of neutrino pdg-codes More...
 
int fgPdgC
 current generated nu pdg-code More...
 
TLorentzVector fgP4
 current generated nu 4-momentum More...
 
TLorentzVector fgX4
 current generated nu 4-position More...
 
double fWeight
 current generated nu weight More...
 
long int fNNeutrinos
 number of flux neutrinos thrown so far More...
 
double fMaxEvCut
 user-defined cut: maximum energy More...
 
double fMinEvCut
 user-defined cut: minimum energy More...
 
double fRl
 defining flux neutrino generation surface: longitudinal radius More...
 
double fRt
 defining flux neutrino generation surface: transverse radius More...
 
TRotation fRotTHz2User
 coord. system rotation: THZ -> Topocentric user-defined More...
 
unsigned int fNumPhiBins
 number of phi bins in input flux data files More...
 
unsigned int fNumCosThetaBins
 number of cos(theta) bins in input flux data files More...
 
unsigned int fNumEnergyBins
 number of energy bins in input flux data files More...
 
double * fPhiBins
 phi bins in input flux data files More...
 
double * fCosThetaBins
 cos(theta) bins in input flux data files More...
 
double * fEnergyBins
 energy bins in input flux data files More...
 
bool fGenWeighted
 generate a weighted or unweighted flux? More...
 
double fSpectralIndex
 power law function used for weighted flux More...
 
bool fInitialized
 flag to check that initialization is run More...
 
TH3D * fTotalFluxHisto
 flux = f(Ev,cos8,phi) summed over neutrino species More...
 
double fTotalFluxHistoIntg
 fFluxSum2D integral More...
 
map< int, TH3D * > fFluxHistoMap
 flux = f(Ev,cos8,phi) for each neutrino species More...
 
map< int, TH3D * > fRawFluxHistoMap
 flux = f(Ev,cos8,phi) for each neutrino species More...
 
vector< intfFluxFlavour
 input flux file for each neutrino species More...
 
vector< string > fFluxFile
 input flux file for each neutrino species More...
 

Private Member Functions

void SetBinSizes (void)
 
bool FillFluxHisto (int nu_pdg, string filename)
 

Detailed Description

A flux driver for the FLUKA 3-D Atmospheric Neutrino Flux.

Astrop.Phys.19 (2003) p.269; hep-ph/0207035; hep-ph/9907408 Alfredo.Ferrari Alfre.nosp@m.do.F.nosp@m.errar.nosp@m.i@ce.nosp@m.rn.ch Paola.Sala Paola.nosp@m..Sal.nosp@m.a@cer.nosp@m.n.ch Giuseppe Battistoni Giuse.nosp@m.ppe..nosp@m.Batti.nosp@m.ston.nosp@m.i@mi..nosp@m.infn.nosp@m..it Teresa Montaruli Teres.nosp@m.a.Mo.nosp@m.ntaru.nosp@m.li@b.nosp@m.a.inf.nosp@m.n.it

To be able to use this flux driver you will need to download the flux data from: http://pcbat1.mi.infn.it/~battist/neutrino.html

Please note that this class expects to read flux files formatted as described in the above FLUKA flux page. Each file contains 3 columns:

Author
Costas Andreopoulos <costas.andreopoulos stfc.ac.uk> University of Liverpool & STFC Rutherford Appleton Lab

July 3, 2005

Copyright (c) 2003-2019, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org or see $GENIE/LICENSE

Definition at line 54 of file GFLUKAAtmoFlux.h.

Constructor & Destructor Documentation

GFLUKAAtmoFlux::GFLUKAAtmoFlux ( )

Definition at line 46 of file GFLUKAAtmoFlux.cxx.

References TMVAGlob::Initialize(), LOG, and pNOTICE.

46  :
47 GAtmoFlux()
48 {
49  LOG("Flux", pNOTICE)
50  << "Instantiating the GENIE FLUKA atmospheric neutrino flux driver";
51 
52  this->Initialize();
53  this->SetBinSizes();
54 }
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
#define pNOTICE
Definition: Messenger.h:62
GFLUKAAtmoFlux::~GFLUKAAtmoFlux ( )

Definition at line 56 of file GFLUKAAtmoFlux.cxx.

57 {
58 
59 }

Member Function Documentation

void GAtmoFlux::AddAllFluxes ( void  )
protectedinherited

Definition at line 574 of file GAtmoFlux.cxx.

References it, LOG, and pNOTICE.

Referenced by genie::flux::GAtmoFlux::MinEnergy().

575 {
576  LOG("Flux", pNOTICE)
577  << "Computing combined flux & flux normalization factor";
578 
579  if(fTotalFluxHisto) delete fTotalFluxHisto;
580 
581  fTotalFluxHisto = this->CreateFluxHisto("sum", "combined flux" );
582 
583  map<int,TH3D*>::iterator it = fFluxHistoMap.begin();
584  for( ; it != fFluxHistoMap.end(); ++it) {
585  TH3D * flux_histogram = it->second;
586  fTotalFluxHisto->Add(flux_histogram);
587  }
588 
589  fTotalFluxHistoIntg = fTotalFluxHisto->Integral();
590 }
map< int, TH3D * > fFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:146
set< int >::iterator it
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double fTotalFluxHistoIntg
fFluxSum2D integral
Definition: GAtmoFlux.h:145
TH3D * fTotalFluxHisto
flux = f(Ev,cos8,phi) summed over neutrino species
Definition: GAtmoFlux.h:144
TH3D * CreateFluxHisto(string name, string title)
Definition: GAtmoFlux.cxx:592
#define pNOTICE
Definition: Messenger.h:62
void GAtmoFlux::AddFluxFile ( int  neutrino_pdg,
string  filename 
)
inherited

Definition at line 421 of file GAtmoFlux.cxx.

References exit(), MakeMiniprodValidationCuts::f, shutoffs::filename, genie::pdg::IsAntiNeutrino(), genie::pdg::IsNeutrino(), LOG, pFATAL, and pWARN.

Referenced by genie::flux::GAtmoFlux::CosZenith(), and GetFlux().

422 {
423  // Check file exists
424  std::ifstream f(filename.c_str());
425  if (!f.good()) {
426  LOG("Flux", pFATAL) << "Flux file does not exist "<<filename;
427  exit(-1);
428  }
429  if ( pdg::IsNeutrino(nu_pdg) || pdg::IsAntiNeutrino(nu_pdg) ) {
430  fFluxFlavour.push_back(nu_pdg);
431  fFluxFile.push_back(filename);
432  } else {
433  LOG ("Flux", pWARN)
434  << "Input particle code: " << nu_pdg << " not a neutrino!";
435  }
436 }
vector< int > fFluxFlavour
input flux file for each neutrino species
Definition: GAtmoFlux.h:148
vector< string > fFluxFile
input flux file for each neutrino species
Definition: GAtmoFlux.h:149
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
#define pFATAL
Definition: Messenger.h:57
string filename
Definition: shutoffs.py:106
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
#define pWARN
Definition: Messenger.h:61
exit(0)
void GAtmoFlux::AddFluxFile ( string  filename)
inherited

Definition at line 438 of file GAtmoFlux.cxx.

References exit(), MakeMiniprodValidationCuts::f, shutoffs::filename, genie::kPdgAntiNuE, genie::kPdgAntiNuMu, genie::kPdgNuE, genie::kPdgNuMu, LOG, and pFATAL.

439 {
440 // FLUKA and BGLRS provide one file per neutrino species.
441 // HAKKKM provides a single file for all nue,nuebar,numu,numubar.
442 // If no neutrino species is provided, assume that the file contains all 4
443 // but fit it into the franework developed for FLUKA and BGLRS,
444 // i.e. add the file 4 times
445 
446 // Check file exists
447  std::ifstream f(filename.c_str());
448  if (!f.good()) {
449  LOG("Flux", pFATAL) << "Flux file does not exist "<<filename;
450  exit(-1);
451  }
452 
453  fFluxFlavour.push_back(kPdgNuE); fFluxFile.push_back(filename);
454  fFluxFlavour.push_back(kPdgAntiNuE); fFluxFile.push_back(filename);
455  fFluxFlavour.push_back(kPdgNuMu); fFluxFile.push_back(filename);
456  fFluxFlavour.push_back(kPdgAntiNuMu); fFluxFile.push_back(filename);
457 
458 }
vector< int > fFluxFlavour
input flux file for each neutrino species
Definition: GAtmoFlux.h:148
vector< string > fFluxFile
input flux file for each neutrino species
Definition: GAtmoFlux.h:149
const int kPdgNuE
Definition: PDGCodes.h:28
const int kPdgAntiNuE
Definition: PDGCodes.h:29
#define pFATAL
Definition: Messenger.h:57
const int kPdgNuMu
Definition: PDGCodes.h:30
string filename
Definition: shutoffs.py:106
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
exit(0)
void GAtmoFlux::CleanUp ( void  )
protectedinherited

Definition at line 379 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

Referenced by genie::flux::GAtmoFlux::CosZenith().

380 {
381  LOG("Flux", pNOTICE) << "Cleaning up...";
382 
383  map<int,TH3D*>::iterator rawiter = fRawFluxHistoMap.begin();
384  for( ; rawiter != fRawFluxHistoMap.end(); ++rawiter) {
385  TH3D * flux_histogram = rawiter->second;
386  if(flux_histogram) {
387  delete flux_histogram;
388  flux_histogram = 0;
389  }
390  }
391  fRawFluxHistoMap.clear();
392 
393  map<int,TH3D*>::iterator iter = fFluxHistoMap.begin();
394  for( ; iter != fFluxHistoMap.end(); ++iter) {
395  TH3D * flux_histogram = iter->second;
396  if(flux_histogram) {
397  delete flux_histogram;
398  flux_histogram = 0;
399  }
400  }
401  fFluxHistoMap.clear();
402 
403  if (fTotalFluxHisto) delete fTotalFluxHisto;
404  if (fPdgCList) delete fPdgCList;
405 
406  if (fPhiBins ) { delete[] fPhiBins ; fPhiBins =NULL; }
407  if (fCosThetaBins) { delete[] fCosThetaBins; fCosThetaBins=NULL; }
408  if (fEnergyBins ) { delete[] fEnergyBins ; fEnergyBins =NULL; }
409 
410 }
map< int, TH3D * > fFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:146
double * fPhiBins
phi bins in input flux data files
Definition: GAtmoFlux.h:138
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double * fCosThetaBins
cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:139
double * fEnergyBins
energy bins in input flux data files
Definition: GAtmoFlux.h:140
map< int, TH3D * > fRawFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:147
TH3D * fTotalFluxHisto
flux = f(Ev,cos8,phi) summed over neutrino species
Definition: GAtmoFlux.h:144
PDGCodeList * fPdgCList
input list of neutrino pdg-codes
Definition: GAtmoFlux.h:124
#define pNOTICE
Definition: Messenger.h:62
void GAtmoFlux::Clear ( Option_t *  opt)
virtualinherited

reset state variables based on opt

Implements genie::GFluxI.

Definition at line 282 of file GAtmoFlux.cxx.

References LOG, MECModelEnuComparisons::opt, and pERROR.

Referenced by genie::flux::GAtmoFlux::Index().

283 {
284 // Dummy clear method needed to conform to GFluxI interface
285 //
286  LOG("Flux", pERROR) << "No clear method implemented for option:"<< opt;
287 }
#define pERROR
Definition: Messenger.h:60
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double genie::flux::GAtmoFlux::CosTheta ( void  )
inlineinherited

Definition at line 82 of file GAtmoFlux.h.

References genie::flux::GAtmoFlux::fgP4.

82 { return -fgP4.Pz()/fgP4.Energy(); }
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:126
double genie::flux::GAtmoFlux::CosZenith ( void  )
inlineinherited
TH3D * GAtmoFlux::CreateFluxHisto ( string  name,
string  title 
)
protectedinherited

Definition at line 592 of file GAtmoFlux.cxx.

References analysePickle::hist, LOG, and pNOTICE.

Referenced by genie::flux::GAtmoFlux::MinEnergy().

593 {
594  LOG("Flux", pNOTICE) << "Instantiating histogram: [" << name << "]";
595  TH3D * hist = new TH3D(
596  name.c_str(), title.c_str(),
600  return hist;
601 }
const XML_Char * name
Definition: expat.h:151
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
unsigned int fNumCosThetaBins
number of cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:136
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double * fCosThetaBins
cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:139
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
#define pNOTICE
Definition: Messenger.h:62
TH3D * GAtmoFlux::CreateNormalisedFluxHisto ( TH3D *  hist)
protectedinherited

Definition at line 520 of file GAtmoFlux.cxx.

References dE.

Referenced by genie::flux::GAtmoFlux::MinEnergy().

521 {
522 // return integrated flux
523 
524  // sanity check
525  if(!hist) return 0;
526 
527  // make new histogram name
528  TString histname = hist->GetName();
529  histname.Append("_IntegratedFlux");
530 
531  // make new histogram
532  TH3D* hist_intg = (TH3D*)(hist->Clone(histname.Data()));
533  hist_intg->Reset();
534 
535  // integrate flux in each bin
536  Double_t dN_dEdS = 0.0;
537  Double_t dS = 0.0;
538  Double_t dE = 0.0;
539  Double_t dN = 0.0;
540 
541  for(Int_t e_bin = 1; e_bin <= hist->GetXaxis()->GetNbins(); e_bin++)
542  {
543  for(Int_t c_bin = 1; c_bin <= hist->GetYaxis()->GetNbins(); c_bin++)
544  {
545  for(Int_t p_bin = 1; p_bin <= hist->GetZaxis()->GetNbins(); p_bin++)
546  {
547  dN_dEdS = hist->GetBinContent(e_bin,c_bin,p_bin);
548 
549  dE = hist->GetXaxis()->GetBinUpEdge(e_bin)
550  - hist->GetXaxis()->GetBinLowEdge(e_bin);
551 
552  dS = ( hist->GetZaxis()->GetBinUpEdge(p_bin)
553  - hist->GetZaxis()->GetBinLowEdge(p_bin) )
554  * ( hist->GetYaxis()->GetBinUpEdge(c_bin)
555  - hist->GetYaxis()->GetBinLowEdge(c_bin) );
556 
557  dN = dN_dEdS*dE*dS;
558 
559  hist_intg->SetBinContent(e_bin,c_bin,p_bin,dN);
560  }
561  }
562  }
563 
564  return hist_intg;
565 }
double dE
virtual bool genie::flux::GAtmoFlux::End ( void  )
inlinevirtualinherited

true if no more flux nu's can be thrown (eg reaching end of beam sim ntuples)

Implements genie::GFluxI.

Definition at line 74 of file GAtmoFlux.h.

74 { return false; }
double genie::flux::GAtmoFlux::Energy ( void  )
inlineinherited

Definition at line 81 of file GAtmoFlux.h.

References genie::flux::GAtmoFlux::fgP4.

81 { return fgP4.Energy(); }
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:126
double genie::flux::GAtmoFlux::Enu ( void  )
inlineinherited

Definition at line 80 of file GAtmoFlux.h.

References genie::flux::GAtmoFlux::fgP4.

80 { return fgP4.Energy(); }
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:126
bool GFLUKAAtmoFlux::FillFluxHisto ( int  nu_pdg,
string  filename 
)
privatevirtual

Implements genie::flux::GAtmoFlux.

Definition at line 124 of file GFLUKAAtmoFlux.cxx.

References energy, shutoffs::filename, flux, genie::flux::GAtmoFlux::fRawFluxHistoMap, histo, make_syst_table_plots::ibin, in, kPi, LOG, pERROR, pINFO, pNOTICE, and scale.

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 }
const double kPi
#define pERROR
Definition: Messenger.h:60
string filename
Definition: shutoffs.py:106
Loaders::FluxType flux
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
double energy
Definition: plottest35.C:25
TH2D * histo
#define pINFO
Definition: Messenger.h:63
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
virtual const PDGCodeList& genie::flux::GAtmoFlux::FluxParticles ( void  )
inlinevirtualinherited

declare list of flux neutrinos that can be generated (for init. purposes)

Implements genie::GFluxI.

Definition at line 67 of file GAtmoFlux.h.

References genie::flux::GAtmoFlux::fPdgCList, genie::flux::GAtmoFlux::GenerateNext(), and genie::flux::GAtmoFlux::MaxEnergy().

67 { return *fPdgCList; }
PDGCodeList * fPdgCList
input list of neutrino pdg-codes
Definition: GAtmoFlux.h:124
void GAtmoFlux::ForceMaxEnergy ( double  emax)
inherited

Definition at line 276 of file GAtmoFlux.cxx.

References emax.

Referenced by genie::flux::GAtmoFlux::CosZenith(), and GetFlux().

277 {
278  emax = TMath::Max(0., emax);
279  fMaxEvCut = emax;
280 }
double fMaxEvCut
user-defined cut: maximum energy
Definition: GAtmoFlux.h:130
const double emax
void GAtmoFlux::ForceMinEnergy ( double  emin)
inherited

Definition at line 270 of file GAtmoFlux.cxx.

References emin.

Referenced by genie::flux::GAtmoFlux::CosZenith(), and GetFlux().

271 {
272  emin = TMath::Max(0., emin);
273  fMinEvCut = emin;
274 }
const double emin
double fMinEvCut
user-defined cut: minimum energy
Definition: GAtmoFlux.h:131
bool GAtmoFlux::GenerateNext ( void  )
virtualinherited

generate the next flux neutrino (return false in err)

Implements genie::GFluxI.

Definition at line 92 of file GAtmoFlux.cxx.

References E, Emax, and make_associated_cosmic_defs::p4.

Referenced by genie::flux::GAtmoFlux::FluxParticles().

93 {
94  while(1) {
95  // Attempt to generate next flux neutrino
96  bool nextok = this->GenerateNext_1try();
97  if(!nextok) continue;
98 
99  // Check generated neutrino energy against max energy.
100  // We may have to reject the current neutrino if a user-defined max
101  // energy cut restricts the available range of energies.
102  const TLorentzVector & p4 = this->Momentum();
103  double E = p4.Energy();
104  double Emin = this->MinEnergy();
105  double Emax = this->MaxEnergy();
106  double wght = this->Weight();
107 
108  bool accept = (E<=Emax && E>=Emin && wght>0);
109  if(accept) return true;
110  }
111  return false;
112 }
virtual double Weight(void)
returns the flux neutrino weight (if any)
Definition: GAtmoFlux.h:71
Float_t E
Definition: plot.C:20
virtual double MaxEnergy(void)
declare the max flux neutrino energy that can be generated (for init. purposes)
Definition: GAtmoFlux.cxx:87
double MinEnergy(void)
Definition: GAtmoFlux.h:112
virtual const TLorentzVector & Momentum(void)
returns the flux neutrino 4-momentum
Definition: GAtmoFlux.h:72
double Emax
bool GenerateNext_1try(void)
Definition: GAtmoFlux.cxx:114
bool GAtmoFlux::GenerateNext_1try ( void  )
protectedinherited

Definition at line 114 of file GAtmoFlux.cxx.

References ana::assert(), visualisationForPaperMasterPlot::ax, emax, emin, slidt::fInitialized, flux, GetFlux(), genie::RandomGen::Instance(), kPi, LOG, genie::utils::print::P4AsShortString(), pINFO, generate_hists::rnd, genie::RandomGen::RndFlux(), ana::Sqrt(), ana::weight, submit_syst::x, genie::utils::print::X4AsString(), submit_syst::y, and test::z.

Referenced by genie::flux::GAtmoFlux::CosZenith().

115 {
116  // Must have run intitialization
118 
119  // Reset previously generated neutrino code / 4-p / 4-x
120  this->ResetSelection();
121 
122  // Get a RandomGen instance
124 
125  // Generate (Ev, costheta, phi)
126  double Ev = 0.;
127  double costheta = 0.;
128  double phi = 0;
129  double weight = 0;
130  int nu_pdg = 0;
131 
132  if(fGenWeighted) {
133 
134  //
135  // generate weighted flux
136  //
137 
138  // generate events according to a power law spectrum,
139  // then weight events by flux and inverse power law
140  // (note: cannot use index alpha=1)
141  double alpha = fSpectralIndex;
142 
143  double emin = TMath::Power(fEnergyBins[0],1.0-alpha);
144  double emax = TMath::Power(fEnergyBins[fNumEnergyBins],1.0-alpha);
145  Ev = TMath::Power(emin+(emax-emin)*rnd->RndFlux().Rndm(),1.0/(1.0-alpha));
146  costheta = -1+2*rnd->RndFlux().Rndm();
147  phi = 2.*kPi* rnd->RndFlux().Rndm();
148 
149  unsigned int nnu = fPdgCList->size();
150  unsigned int inu = rnd->RndFlux().Integer(nnu);
151  nu_pdg = (*fPdgCList)[inu];
152 
153  if(Ev < fEnergyBins[0]) {
154  LOG("Flux", pINFO) << "E < Emin";
155  return false;
156  }
157  double flux = this->GetFlux(nu_pdg, Ev, costheta, phi);
158  if(flux<=0) {
159  LOG("Flux", pINFO) << "Flux <= 0";
160  return false;
161  }
162  weight = flux*TMath::Power(Ev,alpha);
163  }
164  else {
165 
166  //
167  // generate nominal flux
168  //
169 
170  Axis_t ax = 0, ay = 0, az = 0;
171  fTotalFluxHisto->GetRandom3(ax, ay, az);
172  Ev = (double)ax;
173  costheta = (double)ay;
174  phi = (double)az;
175  nu_pdg = this->SelectNeutrino(Ev, costheta, phi);
176  weight = 1.0;
177  }
178 
179  // Compute etc trigonometric numbers
180  double sintheta = TMath::Sqrt(1-costheta*costheta);
181  double cosphi = TMath::Cos(phi);
182  double sinphi = TMath::Sin(phi);
183 
184  // Set the neutrino pdg code
185  fgPdgC = nu_pdg;
186 
187  // Set the neutrino weight
188  fWeight = weight;
189 
190  // Compute the neutrino momentum
191  // The `-1' means it is directed towards the detector.
192  double pz = -1.* Ev * costheta;
193  double py = -1.* Ev * sintheta * sinphi;
194  double px = -1.* Ev * sintheta * cosphi;
195 
196  // Default vertex is at the origin
197  double z = 0.0;
198  double y = 0.0;
199  double x = 0.0;
200 
201  // Shift the neutrino position onto the flux generation surface.
202  // The position is computed at the surface of a sphere with R=fRl
203  // at the topocentric horizontal (THZ) coordinate system.
204  if( fRl>0.0 ){
205  z += fRl * costheta;
206  y += fRl * sintheta * sinphi;
207  x += fRl * sintheta * cosphi;
208  }
209 
210  // Apply user-defined rotation from THZ -> user-defined topocentric
211  // coordinate system.
212  if( !fRotTHz2User.IsIdentity() )
213  {
214  TVector3 tx3(x, y, z );
215  TVector3 tp3(px,py,pz);
216 
217  tx3 = fRotTHz2User * tx3;
218  tp3 = fRotTHz2User * tp3;
219 
220  x = tx3.X();
221  y = tx3.Y();
222  z = tx3.Z();
223  px = tp3.X();
224  py = tp3.Y();
225  pz = tp3.Z();
226  }
227 
228  // If the position is left as is, then all generated neutrinos
229  // would point towards the origin.
230  // Displace the position randomly on the surface that is
231  // perpendicular to the selected point P(xo,yo,zo) on the sphere
232  if( fRt>0.0 ){
233  TVector3 vec(x,y,z); // vector towards selected point
234  TVector3 dvec1 = vec.Orthogonal(); // orthogonal vector
235  TVector3 dvec2 = dvec1; // second orthogonal vector
236  dvec2.Rotate(-kPi/2.0,vec); // rotate second vector by 90deg,
237  // now forming a new orthogonal cartesian coordinate system
238  double psi = 2.*kPi* rnd->RndFlux().Rndm(); // rndm angle [0,2pi]
239  double random = rnd->RndFlux().Rndm(); // rndm number [0,1]
240  dvec1.SetMag(TMath::Sqrt(random)*fRt*TMath::Cos(psi));
241  dvec2.SetMag(TMath::Sqrt(random)*fRt*TMath::Sin(psi));
242  x += dvec1.X() + dvec2.X();
243  y += dvec1.Y() + dvec2.Y();
244  z += dvec1.Z() + dvec2.Z();
245  }
246 
247  // Set the neutrino momentum and position 4-vectors with values
248  // calculated at previous steps.
249  fgP4.SetPxPyPzE(px, py, pz, Ev);
250  fgX4.SetXYZT (x, y, z, 0.);
251 
252  // Increment flux neutrino counter used for sample normalization purposes.
253  fNNeutrinos++;
254 
255  // Report and exit
256  LOG("Flux", pINFO)
257  << "Generated neutrino: "
258  << "\n pdg-code: " << fgPdgC
259  << "\n p4: " << utils::print::P4AsShortString(&fgP4)
260  << "\n x4: " << utils::print::X4AsString(&fgX4);
261 
262  return true;
263 }
int SelectNeutrino(double Ev, double costheta, double phi)
Definition: GAtmoFlux.cxx:603
int fgPdgC
current generated nu pdg-code
Definition: GAtmoFlux.h:125
const double kPi
string P4AsShortString(const TLorentzVector *p)
Definition: PrintUtils.cxx:52
unsigned int fNumEnergyBins
number of energy bins in input flux data files
Definition: GAtmoFlux.h:137
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const Var weight
double fSpectralIndex
power law function used for weighted flux
Definition: GAtmoFlux.h:142
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GAtmoFlux.h:129
double fRl
defining flux neutrino generation surface: longitudinal radius
Definition: GAtmoFlux.h:132
bool fInitialized
flag to check that initialization is run
Definition: GAtmoFlux.h:143
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
double GetFlux(int flavour)
Definition: GAtmoFlux.cxx:662
Loaders::FluxType flux
TRotation fRotTHz2User
coord. system rotation: THZ -> Topocentric user-defined
Definition: GAtmoFlux.h:134
const double emin
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
const double emax
bool fGenWeighted
generate a weighted or unweighted flux?
Definition: GAtmoFlux.h:141
#define pINFO
Definition: Messenger.h:63
Eigen::VectorXd vec
double * fEnergyBins
energy bins in input flux data files
Definition: GAtmoFlux.h:140
z
Definition: test.py:28
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:126
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double fRt
defining flux neutrino generation surface: transverse radius
Definition: GAtmoFlux.h:133
assert(nhit_max >=nhit_nbins)
TH3D * fTotalFluxHisto
flux = f(Ev,cos8,phi) summed over neutrino species
Definition: GAtmoFlux.h:144
void ResetSelection(void)
Definition: GAtmoFlux.cxx:369
PDGCodeList * fPdgCList
input list of neutrino pdg-codes
Definition: GAtmoFlux.h:124
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:72
string X4AsString(const TLorentzVector *x)
Definition: PrintUtils.cxx:64
TLorentzVector fgX4
current generated nu 4-position
Definition: GAtmoFlux.h:127
double fWeight
current generated nu weight
Definition: GAtmoFlux.h:128
void GAtmoFlux::GenerateWeighted ( bool  gen_weighted)
virtualinherited

set whether to generate weighted or unweighted neutrinos

Implements genie::GFluxI.

Definition at line 289 of file GAtmoFlux.cxx.

Referenced by GetFlux(), and genie::flux::GAtmoFlux::Index().

290 {
291  fGenWeighted = gen_weighted;
292 }
bool fGenWeighted
generate a weighted or unweighted flux?
Definition: GAtmoFlux.h:141
double GAtmoFlux::GetFlux ( int  flavour)
inherited

Definition at line 662 of file GAtmoFlux.cxx.

References dE, and flux.

Referenced by genie::flux::GAtmoFlux::CosZenith().

663 {
664  TH3D* flux_hist = this->GetFluxHistogram(flavour);
665  if(!flux_hist) return 0.0;
666 
667  Double_t flux = 0.0;
668  Double_t dN_dEdS = 0.0;
669  Double_t dS = 0.0;
670  Double_t dE = 0.0;
671 
672  for(Int_t e_bin = 1; e_bin <= flux_hist->GetXaxis()->GetNbins(); e_bin++)
673  {
674  for(Int_t c_bin = 1; c_bin <= flux_hist->GetYaxis()->GetNbins(); c_bin++)
675  {
676  for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
677  {
678  dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
679 
680  dE = flux_hist->GetXaxis()->GetBinUpEdge(e_bin)
681  - flux_hist->GetXaxis()->GetBinLowEdge(e_bin);
682 
683  dS = ( flux_hist->GetZaxis()->GetBinUpEdge(p_bin)
684  - flux_hist->GetZaxis()->GetBinLowEdge(p_bin) )
685  * ( flux_hist->GetYaxis()->GetBinUpEdge(c_bin)
686  - flux_hist->GetYaxis()->GetBinLowEdge(c_bin) );
687 
688  flux += dN_dEdS*dE*dS;
689  }
690  }
691  }
692 
693  return flux;
694 }
TH3D * GetFluxHistogram(int flavour)
Definition: GAtmoFlux.cxx:651
double dE
Loaders::FluxType flux
double GAtmoFlux::GetFlux ( int  flavour,
double  energy 
)
inherited

Definition at line 696 of file GAtmoFlux.cxx.

References flux.

697 {
698  TH3D* flux_hist = this->GetFluxHistogram(flavour);
699  if(!flux_hist) return 0.0;
700 
701  Double_t flux = 0.0;
702  Double_t dN_dEdS = 0.0;
703  Double_t dS = 0.0;
704 
705  Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);
706 
707  for(Int_t c_bin = 1; c_bin <= flux_hist->GetYaxis()->GetNbins(); c_bin++)
708  {
709  for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
710  {
711  dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
712 
713  dS = ( flux_hist->GetZaxis()->GetBinUpEdge(p_bin)
714  - flux_hist->GetZaxis()->GetBinLowEdge(p_bin) )
715  * ( flux_hist->GetYaxis()->GetBinUpEdge(c_bin)
716  - flux_hist->GetYaxis()->GetBinLowEdge(c_bin) );
717 
718  flux += dN_dEdS*dS;
719  }
720  }
721 
722  return flux;
723 }
TH3D * GetFluxHistogram(int flavour)
Definition: GAtmoFlux.cxx:651
Loaders::FluxType flux
double energy
Definition: plottest35.C:25
double GAtmoFlux::GetFlux ( int  flavour,
double  energy,
double  costh 
)
inherited

Definition at line 725 of file GAtmoFlux.cxx.

References flux.

726 {
727  TH3D* flux_hist = this->GetFluxHistogram(flavour);
728  if(!flux_hist) return 0.0;
729 
730  Double_t flux = 0.0;
731  Double_t dN_dEdS = 0.0;
732  Double_t dS = 0.0;
733 
734  Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);
735  Int_t c_bin = flux_hist->GetYaxis()->FindBin(costh);
736 
737  for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
738  {
739  dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
740 
741  dS = ( flux_hist->GetZaxis()->GetBinUpEdge(p_bin)
742  - flux_hist->GetZaxis()->GetBinLowEdge(p_bin) );
743 
744  flux += dN_dEdS*dS;
745  }
746 
747  return flux;
748 }
TH3D * GetFluxHistogram(int flavour)
Definition: GAtmoFlux.cxx:651
Loaders::FluxType flux
double energy
Definition: plottest35.C:25
double GAtmoFlux::GetFlux ( int  flavour,
double  energy,
double  costh,
double  phi 
)
inherited

Definition at line 750 of file GAtmoFlux.cxx.

References flux.

751 {
752  TH3D* flux_hist = this->GetFluxHistogram(flavour);
753  if(!flux_hist) return 0.0;
754 
755  Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);
756  Int_t c_bin = flux_hist->GetYaxis()->FindBin(costh);
757  Int_t p_bin = flux_hist->GetZaxis()->FindBin(phi);
758 
759  Double_t flux = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
760  return flux;
761 }
TH3D * GetFluxHistogram(int flavour)
Definition: GAtmoFlux.cxx:651
Loaders::FluxType flux
double energy
Definition: plottest35.C:25
TH3D * GAtmoFlux::GetFluxHistogram ( int  flavour)
inherited

Definition at line 651 of file GAtmoFlux.cxx.

References RunSnowGlobes::histogram, and it.

Referenced by genie::flux::GAtmoFlux::CosZenith().

652 {
653  TH3D* histogram = 0;
654  std::map<int,TH3D*>::iterator it = fRawFluxHistoMap.find(flavour);
655  if(it != fRawFluxHistoMap.end())
656  {
657  histogram = it->second;
658  }
659  return histogram;
660 }
set< int >::iterator it
map< int, TH3D * > fRawFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:147
virtual long int genie::flux::GAtmoFlux::Index ( void  )
inlinevirtualinherited

returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current entry number)

Implements genie::GFluxI.

Definition at line 75 of file GAtmoFlux.h.

References genie::flux::GAtmoFlux::Clear(), genie::flux::GAtmoFlux::GenerateWeighted(), and MECModelEnuComparisons::opt.

75 { return -1; }
void GAtmoFlux::Initialize ( void  )
protectedinherited

Definition at line 311 of file GAtmoFlux.cxx.

References slidt::fInitialized, LOG, and pNOTICE.

Referenced by genie::flux::GAtmoFlux::CosZenith().

312 {
313  LOG("Flux", pNOTICE) << "Initializing atmospheric flux driver";
314 
315  fMaxEv = -1;
316 
317  fNumPhiBins = 0;
318  fNumCosThetaBins = 0;
319  fNumEnergyBins = 0;
320  fPhiBins = 0;
321  fCosThetaBins = 0;
322  fEnergyBins = 0;
323 
324  fTotalFluxHisto = 0;
326 
327  bool allow_dup = false;
328  fPdgCList = new PDGCodeList(allow_dup);
329 
330  // initializing flux TH3D histos [ flux = f(Ev,costheta,phi) ] & files
331  fFluxFile.clear();
332  fFluxHistoMap.clear();
333  fTotalFluxHisto = 0;
335 
336  // Default option is to generate unweighted flux neutrinos
337  // (flux = f(E,costheta) will be used as PDFs)
338  // User can enable option to generate weighted neutrinos
339  // (neutrinos will be generated uniformly over costheta,
340  // and using a power law function in neutrino energy.
341  // The input flux = f(E,costheta) will be used for calculating a weight).
342  // Using a weighted flux avoids statistical fluctuations at high energies.
343  fSpectralIndex = 2.0;
344 
345  // weighting switched off by default
346  this->GenerateWeighted(false);
347 
348  // Default: No min/max energy cut
349  this->ForceMinEnergy(0.);
350  this->ForceMaxEnergy(9999999999.);
351 
352  // Default radii
353  fRl = 0.0;
354  fRt = 0.0;
355 
356  // Default detector coord system: Topocentric Horizontal Coordinate system
357  fRotTHz2User.SetToIdentity();
358 
359  // Reset `current' selected flux neutrino
360  this->ResetSelection();
361 
362  // Reset number of neutrinos thrown so far
363  fNNeutrinos = 0;
364 
365  // Done!
366  fInitialized = 1;
367 }
map< int, TH3D * > fFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:146
double * fPhiBins
phi bins in input flux data files
Definition: GAtmoFlux.h:138
vector< string > fFluxFile
input flux file for each neutrino species
Definition: GAtmoFlux.h:149
unsigned int fNumEnergyBins
number of energy bins in input flux data files
Definition: GAtmoFlux.h:137
unsigned int fNumCosThetaBins
number of cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:136
double fSpectralIndex
power law function used for weighted flux
Definition: GAtmoFlux.h:142
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GAtmoFlux.h:129
double fRl
defining flux neutrino generation surface: longitudinal radius
Definition: GAtmoFlux.h:132
bool fInitialized
flag to check that initialization is run
Definition: GAtmoFlux.h:143
virtual void GenerateWeighted(bool gen_weighted)
set whether to generate weighted or unweighted neutrinos
Definition: GAtmoFlux.cxx:289
TRotation fRotTHz2User
coord. system rotation: THZ -> Topocentric user-defined
Definition: GAtmoFlux.h:134
A list of PDG codes.
Definition: PDGCodeList.h:33
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void ForceMaxEnergy(double emax)
Definition: GAtmoFlux.cxx:276
double fTotalFluxHistoIntg
fFluxSum2D integral
Definition: GAtmoFlux.h:145
double * fCosThetaBins
cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:139
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
double fRt
defining flux neutrino generation surface: transverse radius
Definition: GAtmoFlux.h:133
TH3D * fTotalFluxHisto
flux = f(Ev,cos8,phi) summed over neutrino species
Definition: GAtmoFlux.h:144
void ResetSelection(void)
Definition: GAtmoFlux.cxx:369
PDGCodeList * fPdgCList
input list of neutrino pdg-codes
Definition: GAtmoFlux.h:124
#define pNOTICE
Definition: Messenger.h:62
void ForceMinEnergy(double emin)
Definition: GAtmoFlux.cxx:270
double fMaxEv
maximum energy (in input flux files)
Definition: GAtmoFlux.h:123
bool GAtmoFlux::LoadFluxData ( void  )
inherited

Definition at line 465 of file GAtmoFlux.cxx.

References shutoffs::filename, genie::PDGLibrary::Find(), analysePickle::hist, hnorm, genie::PDGLibrary::Instance(), LOG, getGoodRuns4SAM::n, pERROR, eplot::pname, and pNOTICE.

Referenced by genie::flux::GAtmoFlux::CosZenith(), and GetFlux().

466 {
467  LOG("Flux", pNOTICE)
468  << "Loading atmospheric neutrino flux simulation data";
469 
470  fFluxHistoMap.clear();
471  fPdgCList->clear();
472 
473  bool loading_status = true;
474 
475  for( unsigned int n=0; n<fFluxFlavour.size(); n++ ){
476  int nu_pdg = fFluxFlavour.at(n);
477  string filename = fFluxFile.at(n);
478  string pname = PDGLibrary::Instance()->Find(nu_pdg)->GetName();
479 
480  LOG("Flux", pNOTICE) << "Loading data for: " << pname;
481 
482  // create histogram
483  TH3D* hist = 0;
484  std::map<int,TH3D*>::iterator myMapEntry = fRawFluxHistoMap.find(nu_pdg);
485  if( myMapEntry == fRawFluxHistoMap.end() ){
486 // hist = myMapEntry->second;
487 // if(hist==0) {
488  hist = this->CreateFluxHisto(pname.c_str(), pname.c_str());
489  fRawFluxHistoMap.insert( map<int,TH3D*>::value_type(nu_pdg,hist) );
490 // }
491  }
492  // now let concrete instances to read the flux-specific data files
493  // and fill the histogram
494  bool loaded = this->FillFluxHisto(nu_pdg, filename);
495  loading_status = loading_status && loaded;
496  }
497 
498  if(loading_status) {
499  map<int,TH3D*>::iterator hist_iter = fRawFluxHistoMap.begin();
500  for ( ; hist_iter != fRawFluxHistoMap.end(); ++hist_iter) {
501  int nu_pdg = hist_iter->first;
502  TH3D* hist = hist_iter->second;
503 
504  TH3D* hnorm = this->CreateNormalisedFluxHisto( hist );
505  fFluxHistoMap.insert( map<int,TH3D*>::value_type(nu_pdg,hnorm) );
506  fPdgCList->push_back(nu_pdg);
507  }
508 
509  LOG("Flux", pNOTICE)
510  << "Atmospheric neutrino flux simulation data loaded!";
511  this->AddAllFluxes();
512  return true;
513  }
514 
515  LOG("Flux", pERROR)
516  << "Error loading atmospheric neutrino flux simulation data";
517  return false;
518 }
vector< int > fFluxFlavour
input flux file for each neutrino species
Definition: GAtmoFlux.h:148
map< int, TH3D * > fFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:146
vector< string > fFluxFile
input flux file for each neutrino species
Definition: GAtmoFlux.h:149
#define pERROR
Definition: Messenger.h:60
TH3D * CreateNormalisedFluxHisto(TH3D *hist)
Definition: GAtmoFlux.cxx:520
string filename
Definition: shutoffs.py:106
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
TH1F * hnorm
Definition: plots_total.C:4
virtual bool FillFluxHisto(int nu_pdg, string filename)=0
string pname
Definition: eplot.py:33
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
map< int, TH3D * > fRawFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:147
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
TH3D * CreateFluxHisto(string name, string title)
Definition: GAtmoFlux.cxx:592
PDGCodeList * fPdgCList
input list of neutrino pdg-codes
Definition: GAtmoFlux.h:124
#define pNOTICE
Definition: Messenger.h:62
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
double GAtmoFlux::MaxEnergy ( void  )
virtualinherited

declare the max flux neutrino energy that can be generated (for init. purposes)

Implements genie::GFluxI.

Definition at line 87 of file GAtmoFlux.cxx.

Referenced by genie::flux::GAtmoFlux::FluxParticles().

88 {
89  return TMath::Min(fMaxEv, fMaxEvCut);
90 }
double fMaxEvCut
user-defined cut: maximum energy
Definition: GAtmoFlux.h:130
double fMaxEv
maximum energy (in input flux files)
Definition: GAtmoFlux.h:123
double genie::flux::GAtmoFlux::MinEnergy ( void  )
inlineprotectedinherited
virtual const TLorentzVector& genie::flux::GAtmoFlux::Momentum ( void  )
inlinevirtualinherited

returns the flux neutrino 4-momentum

Implements genie::GFluxI.

Definition at line 72 of file GAtmoFlux.h.

References genie::flux::GAtmoFlux::fgP4.

72 { return fgP4; }
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:126
long int GAtmoFlux::NFluxNeutrinos ( void  ) const
inherited

Number of flux nu's generated. Not the same as the number of nu's thrown towards the geometry (if there are cuts).

Definition at line 265 of file GAtmoFlux.cxx.

Referenced by genie::flux::GAtmoFlux::CosZenith().

266 {
267  return fNNeutrinos;
268 }
long int fNNeutrinos
number of flux neutrinos thrown so far
Definition: GAtmoFlux.h:129
virtual int genie::flux::GAtmoFlux::PdgCode ( void  )
inlinevirtualinherited

returns the flux neutrino pdg code

Implements genie::GFluxI.

Definition at line 70 of file GAtmoFlux.h.

References genie::flux::GAtmoFlux::fgPdgC.

70 { return fgPdgC; }
int fgPdgC
current generated nu pdg-code
Definition: GAtmoFlux.h:125
virtual const TLorentzVector& genie::flux::GAtmoFlux::Position ( void  )
inlinevirtualinherited

returns the flux neutrino 4-position (note: expect SI rather than physical units)

Implements genie::GFluxI.

Definition at line 73 of file GAtmoFlux.h.

References genie::flux::GAtmoFlux::fgX4.

73 { return fgX4; }
TLorentzVector fgX4
current generated nu 4-position
Definition: GAtmoFlux.h:127
void GAtmoFlux::ResetSelection ( void  )
protectedinherited

Definition at line 369 of file GAtmoFlux.cxx.

Referenced by genie::flux::GAtmoFlux::CosZenith().

370 {
371 // initializing running neutrino pdg-code, 4-position, 4-momentum
372 
373  fgPdgC = 0;
374  fgP4.SetPxPyPzE (0.,0.,0.,0.);
375  fgX4.SetXYZT (0.,0.,0.,0.);
376  fWeight = 0;
377 }
int fgPdgC
current generated nu pdg-code
Definition: GAtmoFlux.h:125
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:126
TLorentzVector fgX4
current generated nu 4-position
Definition: GAtmoFlux.h:127
double fWeight
current generated nu weight
Definition: GAtmoFlux.h:128
int GAtmoFlux::SelectNeutrino ( double  Ev,
double  costheta,
double  phi 
)
protectedinherited

Definition at line 603 of file GAtmoFlux.cxx.

References ana::assert(), flux, MECModelEnuComparisons::i, make_syst_table_plots::ibin, genie::RandomGen::Instance(), it, LOG, getGoodRuns4SAM::n, pDEBUG, pERROR, pINFO, R, generate_hists::rnd, and genie::RandomGen::RndFlux().

Referenced by genie::flux::GAtmoFlux::MinEnergy().

604 {
605 // Select a neutrino species at the input Ev and costheta given their
606 // relatve flux at this bin.
607 // Returns a neutrino PDG code
608 
609  unsigned int n = fPdgCList->size();
610  double * flux = new double[n];
611 
612  unsigned int i=0;
613  map<int,TH3D*>::iterator it = fFluxHistoMap.begin();
614  for( ; it != fFluxHistoMap.end(); ++it) {
615  TH3D * flux_histogram = it->second;
616  int ibin = flux_histogram->FindBin(Ev,costheta,phi);
617  flux[i] = flux_histogram->GetBinContent(ibin);
618  i++;
619  }
620  double flux_sum = 0;
621  for(i=0; i<n; i++) {
622  flux_sum += flux[i];
623  flux[i] = flux_sum;
624  LOG("Flux", pDEBUG)
625  << "Sum{Flux(0->" << i <<")} = " << flux[i];
626  }
627 
629  double R = flux_sum * rnd->RndFlux().Rndm();
630 
631  LOG("Flux", pDEBUG) << "R = " << R;
632 
633  for(i=0; i<n; i++) {
634  if( R < flux[i] ) {
635  delete [] flux;
636  int selected_pdg = (*fPdgCList)[i];
637  LOG("Flux", pINFO)
638  << "Selected neutrino PDG code = " << selected_pdg;
639  return selected_pdg;
640  }
641  }
642 
643  // shouldn't be here
644  LOG("Flux", pERROR) << "Could not select a neutrino species!";
645  assert(false);
646 
647  return -1;
648 }
map< int, TH3D * > fFluxHistoMap
flux = f(Ev,cos8,phi) for each neutrino species
Definition: GAtmoFlux.h:146
set< int >::iterator it
#define pERROR
Definition: Messenger.h:60
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
Loaders::FluxType flux
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
#define R(x)
#define pINFO
Definition: Messenger.h:63
assert(nhit_max >=nhit_nbins)
PDGCodeList * fPdgCList
input list of neutrino pdg-codes
Definition: GAtmoFlux.h:124
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:72
#define pDEBUG
Definition: Messenger.h:64
void GFLUKAAtmoFlux::SetBinSizes ( void  )
private

Definition at line 61 of file GFLUKAAtmoFlux.cxx.

References genie::flux::GAtmoFlux::fCosThetaBins, genie::flux::GAtmoFlux::fEnergyBins, genie::flux::GAtmoFlux::fMaxEv, genie::flux::GAtmoFlux::fNumCosThetaBins, genie::flux::GAtmoFlux::fNumEnergyBins, genie::flux::GAtmoFlux::fNumPhiBins, genie::flux::GAtmoFlux::fPhiBins, MECModelEnuComparisons::i, genie::flux::kGFlk3DCosThetaMax, genie::flux::kGFlk3DCosThetaMin, genie::flux::kGFlk3DEvMin, genie::flux::kGFlk3DNumCosThetaBins, genie::flux::kGFlk3DNumLogEvBins, genie::flux::kGFlk3DNumLogEvBinsPerDecade, kPi, LOG, and pDEBUG.

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 }
const double kPi
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
unsigned int fNumCosThetaBins
number of cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:136
const double kGFlk3DCosThetaMax
const unsigned int kGFlk3DNumCosThetaBins
const unsigned int kGFlk3DNumLogEvBins
const unsigned int kGFlk3DNumLogEvBinsPerDecade
#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 * fCosThetaBins
cos(theta) bins in input flux data files
Definition: GAtmoFlux.h:139
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
double fMaxEv
maximum energy (in input flux files)
Definition: GAtmoFlux.h:123
#define pDEBUG
Definition: Messenger.h:64
void GAtmoFlux::SetRadii ( double  Rlongitudinal,
double  Rtransverse 
)
inherited

Definition at line 412 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

Referenced by genie::flux::GAtmoFlux::CosZenith(), and GetFlux().

413 {
414  LOG ("Flux", pNOTICE) << "Setting R[longitudinal] = " << Rlongitudinal;
415  LOG ("Flux", pNOTICE) << "Setting R[transverse] = " << Rtransverse;
416 
417  fRl = Rlongitudinal;
418  fRt = Rtransverse;
419 }
double fRl
defining flux neutrino generation surface: longitudinal radius
Definition: GAtmoFlux.h:132
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
double fRt
defining flux neutrino generation surface: transverse radius
Definition: GAtmoFlux.h:133
#define pNOTICE
Definition: Messenger.h:62
void GAtmoFlux::SetSpectralIndex ( double  index)
inherited

Definition at line 294 of file GAtmoFlux.cxx.

References allTimeWatchdog::index, LOG, pNOTICE, and pWARN.

Referenced by genie::flux::GAtmoFlux::CosZenith().

295 {
296  if( index != 1.0 ){
298  }
299  else {
300  LOG("Flux", pWARN) << "Warning: cannot use a spectral index of unity";
301  }
302 
303  LOG("Flux", pNOTICE) << "Using Spectral Index = " << index;
304 }
double fSpectralIndex
power law function used for weighted flux
Definition: GAtmoFlux.h:142
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
#define pWARN
Definition: Messenger.h:61
#define pNOTICE
Definition: Messenger.h:62
void GAtmoFlux::SetUserCoordSystem ( TRotation &  rotation)
inherited

Rotation: Topocentric Horizontal -> User-defined Topocentric Coord System.

Definition at line 306 of file GAtmoFlux.cxx.

Referenced by genie::flux::GAtmoFlux::CosZenith(), and GetFlux().

307 {
308  fRotTHz2User = rotation;
309 }
TRotation fRotTHz2User
coord. system rotation: THZ -> Topocentric user-defined
Definition: GAtmoFlux.h:134
virtual double genie::flux::GAtmoFlux::Weight ( void  )
inlinevirtualinherited

returns the flux neutrino weight (if any)

Implements genie::GFluxI.

Definition at line 71 of file GAtmoFlux.h.

References genie::flux::GAtmoFlux::fWeight.

71 { return fWeight; }
double fWeight
current generated nu weight
Definition: GAtmoFlux.h:128
void GAtmoFlux::ZeroFluxHisto ( TH3D *  hist)
protectedinherited

Definition at line 567 of file GAtmoFlux.cxx.

References LOG, and pDEBUG.

Referenced by genie::flux::GAtmoFlux::MinEnergy().

568 {
569  LOG("Flux", pDEBUG) << "Forcing flux histogram contents to 0";
570  if(!histo) return;
571  histo->Reset();
572 }
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
TH2D * histo
#define pDEBUG
Definition: Messenger.h:64

Member Data Documentation

double* genie::flux::GAtmoFlux::fCosThetaBins
protectedinherited
double* genie::flux::GAtmoFlux::fEnergyBins
protectedinherited
vector<string> genie::flux::GAtmoFlux::fFluxFile
protectedinherited

input flux file for each neutrino species

Definition at line 149 of file GAtmoFlux.h.

vector<int> genie::flux::GAtmoFlux::fFluxFlavour
protectedinherited

input flux file for each neutrino species

Definition at line 148 of file GAtmoFlux.h.

map<int, TH3D*> genie::flux::GAtmoFlux::fFluxHistoMap
protectedinherited

flux = f(Ev,cos8,phi) for each neutrino species

Definition at line 146 of file GAtmoFlux.h.

bool genie::flux::GAtmoFlux::fGenWeighted
protectedinherited

generate a weighted or unweighted flux?

Definition at line 141 of file GAtmoFlux.h.

TLorentzVector genie::flux::GAtmoFlux::fgP4
protectedinherited
int genie::flux::GAtmoFlux::fgPdgC
protectedinherited

current generated nu pdg-code

Definition at line 125 of file GAtmoFlux.h.

Referenced by genie::flux::GAtmoFlux::PdgCode().

TLorentzVector genie::flux::GAtmoFlux::fgX4
protectedinherited

current generated nu 4-position

Definition at line 127 of file GAtmoFlux.h.

Referenced by genie::flux::GAtmoFlux::Position().

bool genie::flux::GAtmoFlux::fInitialized
protectedinherited

flag to check that initialization is run

Definition at line 143 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fMaxEv
protectedinherited

maximum energy (in input flux files)

Definition at line 123 of file GAtmoFlux.h.

Referenced by genie::flux::GHAKKMAtmoFlux::SetBinSizes(), SetBinSizes(), and genie::flux::GBGLRSAtmoFlux::SetBinSizes().

double genie::flux::GAtmoFlux::fMaxEvCut
protectedinherited

user-defined cut: maximum energy

Definition at line 130 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fMinEvCut
protectedinherited

user-defined cut: minimum energy

Definition at line 131 of file GAtmoFlux.h.

Referenced by genie::flux::GAtmoFlux::MinEnergy().

long int genie::flux::GAtmoFlux::fNNeutrinos
protectedinherited

number of flux neutrinos thrown so far

Definition at line 129 of file GAtmoFlux.h.

unsigned int genie::flux::GAtmoFlux::fNumCosThetaBins
protectedinherited

number of cos(theta) bins in input flux data files

Definition at line 136 of file GAtmoFlux.h.

Referenced by genie::flux::GHAKKMAtmoFlux::SetBinSizes(), SetBinSizes(), and genie::flux::GBGLRSAtmoFlux::SetBinSizes().

unsigned int genie::flux::GAtmoFlux::fNumEnergyBins
protectedinherited

number of energy bins in input flux data files

Definition at line 137 of file GAtmoFlux.h.

Referenced by genie::flux::GHAKKMAtmoFlux::SetBinSizes(), SetBinSizes(), and genie::flux::GBGLRSAtmoFlux::SetBinSizes().

unsigned int genie::flux::GAtmoFlux::fNumPhiBins
protectedinherited

number of phi bins in input flux data files

Definition at line 135 of file GAtmoFlux.h.

Referenced by genie::flux::GHAKKMAtmoFlux::SetBinSizes(), SetBinSizes(), and genie::flux::GBGLRSAtmoFlux::SetBinSizes().

PDGCodeList* genie::flux::GAtmoFlux::fPdgCList
protectedinherited

input list of neutrino pdg-codes

Definition at line 124 of file GAtmoFlux.h.

Referenced by genie::flux::GAtmoFlux::FluxParticles().

double* genie::flux::GAtmoFlux::fPhiBins
protectedinherited
map<int, TH3D*> genie::flux::GAtmoFlux::fRawFluxHistoMap
protectedinherited

flux = f(Ev,cos8,phi) for each neutrino species

Definition at line 147 of file GAtmoFlux.h.

Referenced by genie::flux::GHAKKMAtmoFlux::FillFluxHisto(), FillFluxHisto(), and genie::flux::GBGLRSAtmoFlux::FillFluxHisto().

double genie::flux::GAtmoFlux::fRl
protectedinherited

defining flux neutrino generation surface: longitudinal radius

Definition at line 132 of file GAtmoFlux.h.

TRotation genie::flux::GAtmoFlux::fRotTHz2User
protectedinherited

coord. system rotation: THZ -> Topocentric user-defined

Definition at line 134 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fRt
protectedinherited

defining flux neutrino generation surface: transverse radius

Definition at line 133 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fSpectralIndex
protectedinherited

power law function used for weighted flux

Definition at line 142 of file GAtmoFlux.h.

TH3D* genie::flux::GAtmoFlux::fTotalFluxHisto
protectedinherited

flux = f(Ev,cos8,phi) summed over neutrino species

Definition at line 144 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fTotalFluxHistoIntg
protectedinherited

fFluxSum2D integral

Definition at line 145 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fWeight
protectedinherited

current generated nu weight

Definition at line 128 of file GAtmoFlux.h.

Referenced by genie::flux::GAtmoFlux::Weight().


The documentation for this class was generated from the following files: