A flux driver for the Bartol Atmospheric Neutrino Flux. More...

#include <GBGLRSAtmoFlux.h>

Inheritance diagram for genie::flux::GBGLRSAtmoFlux:

Collaboration diagram for genie::flux::GBGLRSAtmoFlux:

Public Member Functions
	GBGLRSAtmoFlux ()

	~GBGLRSAtmoFlux ()

Public Member Functions inherited from genie::flux::GAtmoFlux
virtual	~GAtmoFlux ()

virtual const PDGCodeList &	FluxParticles (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)

double	GetFluxSurfaceArea (void)

double	GetLongitudinalRadius (void)

double	GetTransverseRadius (void)

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	GetTotalFlux (void)

double	GetTotalFluxInEnergyRange (void)

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)

Public Member Functions inherited from genie::GFluxI
virtual	~GFluxI ()

Private Member Functions
void	SetBinSizes (void)

bool	FillFluxHisto (int nu_pdg, string filename)

Additional Inherited Members
Protected Member Functions inherited from genie::flux::GAtmoFlux
	GAtmoFlux ()

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 Member Functions inherited from genie::GFluxI
	GFluxI ()

Protected Attributes inherited from genie::flux::GAtmoFlux
double	fMaxEv
	maximum energy (in input flux files) More...

PDGCodeList *	fPdgCList
	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< int >	fFluxFlavour
	input flux file for each neutrino species More...

vector< string >	fFluxFile
	input flux file for each neutrino species More...

Detailed Description

A flux driver for the Bartol Atmospheric Neutrino Flux.

References:: G. Barr, T.K. Gaisser, P. Lipari, S. Robbins and T. Stanev, astro-ph/0403630

To be able to use this flux driver you will need to download the flux data from: http://www-pnp.physics.ox.ac.uk/~barr/fluxfiles/

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

neutrino energy (GeV) at bin centre
neutrino cos(zenith angle) at bin centre
neutrino flux dN/dlnE (#neutrinos /m^2 /sec /sr)
MC statistical error on the flux (not used here)
Number of unweighted events entering in the bin (not used here) The flux is given in 20 bins of cos(zenith angle) from -1.0 to 1.0 (bin width = 0.1) and 30 equally log-spaced energy bins (10 bins per decade), with Emin = 10.00 GeV.

Note that in the BGLRS input files the flux is defined as dN/dlnE, while in the FLUKA files the flux is defined as dN/dE. We compensate for logarithmic units (dlnE = dE/E) as we read-in the BGLRS files.

Author: Christopher Backhouse c.bac.nosp@m.khou.nosp@m.se1@p.nosp@m.hysi.nosp@m.cs.ox.nosp@m..ac..nosp@m.uk Oxford University

Created:: January 26, 2008

License:: Copyright (c) 2003-2024, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Definition at line 59 of file GBGLRSAtmoFlux.h.

Constructor & Destructor Documentation

GBGLRSAtmoFlux::GBGLRSAtmoFlux ( )

Definition at line 32 of file GBGLRSAtmoFlux.cxx.

References Initialize(), LOG, and pNOTICE.

                                :
 GAtmoFlux()
 {
   LOG("Flux", pNOTICE)
      << "Instantiating the GENIE BGLRS atmospheric neutrino flux driver";
 
   this->Initialize();
   this->SetBinSizes();
 }

GBGLRSAtmoFlux::~GBGLRSAtmoFlux ( )

Definition at line 42 of file GBGLRSAtmoFlux.cxx.

 {
 
 }

Member Function Documentation

bool GBGLRSAtmoFlux::FillFluxHisto	(	int	nu_pdg,
		string	filename
	)

privatevirtual

Implements genie::flux::GAtmoFlux.

Definition at line 119 of file GBGLRSAtmoFlux.cxx.

References genie::flux::GAtmoFlux::fRawFluxHistoMap, genie::constants::kPi, LOG, pERROR, pINFO, and pNOTICE.

 {
   LOG("Flux", pNOTICE)
     << "Loading BGLRS flux for neutrino: " << nu_pdg
     << " from file: " << filename;
 
   TH3D* histo = 0;
   std::map<int,TH3D*>::iterator myMapEntry = fRawFluxHistoMap.find(nu_pdg);
   if( myMapEntry != fRawFluxHistoMap.end() ){
       histo = myMapEntry->second;
   }
   if(!histo) {
      LOG("Flux", pERROR) << "Null flux histogram!";
      return false;
   }
   ifstream flux_stream(filename.c_str(), ios::in);
   if(!flux_stream) {
      LOG("Flux", pERROR) << "Could not open file: " << filename;
      return false;
   }
 
   int ibin;
   double energy, costheta, flux;
   double junkd; // throw away error estimates
 
   double scale = 1.0; // 1.0 [m^2], OR 1.0e-4 [cm^2]
 
   // throw away comment line
   flux_stream.ignore(99999, '\n');
 
   while ( !flux_stream.eof() ) {
     flux = 0.0;
     flux_stream >> energy >> costheta >> flux >> junkd >> junkd;
     if( flux>0.0 ){
       // Compensate for logarithmic units - dlogE=dE/E
       // [Note: should do this explicitly using bin widths]
       flux /= energy;
       LOG("Flux", pINFO)
         << "Flux[Ev = " << energy
         << ", cos8 = " << costheta << "] = " << flux;
       ibin = histo->FindBin( (Axis_t)energy, (Axis_t)costheta, (Axis_t)kPi );
       histo->SetBinContent( ibin, (Stat_t)(scale*flux) );
     }
   }
   return true;
 }

void GBGLRSAtmoFlux::SetBinSizes ( void )

private

Definition at line 47 of file GBGLRSAtmoFlux.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, genie::flux::kBGLRS3DCosThetaMax, genie::flux::kBGLRS3DCosThetaMin, genie::flux::kBGLRS3DEvMin, genie::flux::kBGLRS3DNumCosThetaBins, genie::flux::kBGLRS3DNumLogEvBinsHigh, genie::flux::kBGLRS3DNumLogEvBinsLow, genie::flux::kBGLRS3DNumLogEvBinsPerDecadeHigh, genie::flux::kBGLRS3DNumLogEvBinsPerDecadeLow, genie::constants::kPi, LOG, and pDEBUG.

 {
 // Generate the correct cos(theta) and energy bin sizes.
 //
 // Zenith angle binning: the flux is given in 20 bins of
 // cos(zenith angle) from -1.0 to 1.0 (bin width = 0.1)
 //
 // Neutrino energy binning: the Bartol flux files are
 // provided in two pieces
 //  (1) low energy piece (<10 GeV), solar min or max,
 //      given in 40 log-spaced bins from 0.1 to 10 GeV
 //      (20 bins per decade)
 //  (2) high energy piece (>10 GeV), without solar effects,
 //      given in 30 log-spaced bins from 10 to 1000 GeV
 //      (10 bins per decade)
 
   fPhiBins       = new double [2];
   fCosThetaBins  = new double [kBGLRS3DNumCosThetaBins + 1];
   fEnergyBins    = new double [kBGLRS3DNumLogEvBinsLow + kBGLRS3DNumLogEvBinsHigh + 1];
 
   fPhiBins[0] = 0;
   fPhiBins[1] = 2.*kPi;
 
   double dcostheta =
       (kBGLRS3DCosThetaMax - kBGLRS3DCosThetaMin) /
       (double) kBGLRS3DNumCosThetaBins;
 
   double logEmin = TMath::Log10(kBGLRS3DEvMin);
   double dlogElow = 1.0 / (double) kBGLRS3DNumLogEvBinsPerDecadeLow;
   double dlogEhigh = 1.0 / (double) kBGLRS3DNumLogEvBinsPerDecadeHigh;
 
   double costheta = kBGLRS3DCosThetaMin;
 
   for(unsigned int i=0; i<= kBGLRS3DNumCosThetaBins; i++) {
     if( i==0 ) ; // do nothing
     else costheta += dcostheta;
     fCosThetaBins[i] = costheta;
     if(i != kBGLRS3DNumCosThetaBins) {
       LOG("Flux", pDEBUG)
         << "FLUKA 3d flux: CosTheta bin " << i+1
         << ": lower edge = " << fCosThetaBins[i];
     } else {
       LOG("Flux", pDEBUG)
         << "FLUKA 3d flux: CosTheta bin " << kBGLRS3DNumCosThetaBins
         << ": upper edge = " << fCosThetaBins[kBGLRS3DNumCosThetaBins];
     }
   }
 
   double logE = logEmin;
 
   for(unsigned int i=0; i<=kBGLRS3DNumLogEvBinsLow+kBGLRS3DNumLogEvBinsHigh; i++) {
     if( i==0 ) ; // do nothing
     else if( i<=kBGLRS3DNumLogEvBinsLow ) logE += dlogElow;
     else                                  logE += dlogEhigh;
     fEnergyBins[i] = TMath::Power(10.0, logE);
     if(i != kBGLRS3DNumLogEvBinsLow+kBGLRS3DNumLogEvBinsHigh) {
       LOG("Flux", pDEBUG)
          << "FLUKA 3d flux: Energy bin " << i+1
          << ": lower edge = " << fEnergyBins[i];
     } else {
       LOG("Flux", pDEBUG)
          << "FLUKA 3d flux: Energy bin " << kBGLRS3DNumLogEvBinsLow+kBGLRS3DNumLogEvBinsHigh
          << ": upper edge = " << fEnergyBins[kBGLRS3DNumLogEvBinsLow+kBGLRS3DNumLogEvBinsHigh];
     }
   }
 
   fNumPhiBins      = 1;
   fNumCosThetaBins = kBGLRS3DNumCosThetaBins;
   fNumEnergyBins   = kBGLRS3DNumLogEvBinsLow + kBGLRS3DNumLogEvBinsHigh;
   fMaxEv = fEnergyBins[fNumEnergyBins];
 }

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

Public Member Functions

Private Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation