A driver for the HAKKM 3-D atmospheric neutrino flux (commonly known as the `Honda flux') More...

#include <GHAKKMAtmoFlux.h>

Inheritance diagram for genie::flux::GHAKKMAtmoFlux:

Collaboration diagram for genie::flux::GHAKKMAtmoFlux:

Public Member Functions
	GHAKKMAtmoFlux ()

	~GHAKKMAtmoFlux ()

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 driver for the HAKKM 3-D atmospheric neutrino flux (commonly known as the `Honda flux')

References:: M. Honda, M. Sajjad Athar, T. Kajita, K. Kasahara, and S. Midorikawa Phys. Rev. D 92 (2015) 023004

The flux files necessary for running this flux driver can be obtained from:http://www.icrr.u-tokyo.ac.jp/~mhonda/

Author: Costas Andreopoulos <c.andreopoulos cern.ch> University of Liverpool

Created:: July 9, 2015

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

Definition at line 45 of file GHAKKMAtmoFlux.h.

Constructor & Destructor Documentation

GHAKKMAtmoFlux::GHAKKMAtmoFlux ( )

Definition at line 32 of file GHAKKMAtmoFlux.cxx.

References Initialize(), LOG, and pNOTICE.

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

GHAKKMAtmoFlux::~GHAKKMAtmoFlux ( )

Definition at line 42 of file GHAKKMAtmoFlux.cxx.

 {
 
 }

Member Function Documentation

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

privatevirtual

Implements genie::flux::GAtmoFlux.

Definition at line 135 of file GHAKKMAtmoFlux.cxx.

References genie::flux::GAtmoFlux::fCosThetaBins, genie::flux::GAtmoFlux::fEnergyBins, genie::flux::GAtmoFlux::fPhiBins, genie::flux::GAtmoFlux::fRawFluxHistoMap, genie::flux::kGHnd3DNumCosThetaBins, genie::flux::kGHnd3DNumLogEvBins, genie::flux::kGHnd3DNumPhiBins, genie::kPdgAntiNuE, genie::kPdgAntiNuMu, genie::kPdgNuE, genie::kPdgNuMu, genie::constants::kPi, LOG, pDEBUG, pERROR, and pNOTICE.

 {
   LOG("Flux", pNOTICE)
     << "Loading HAKKM 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;
   }
 
   double r2d = 180./constants::kPi;
 
   int icostheta = kGHnd3DNumCosThetaBins; // starting from last bin [0.9 - 1.0]
   int iphi      = 0;
 
   while (!flux_stream.eof()) {
 
     string comment = "";
     getline(flux_stream,comment);
     LOG("Flux", pDEBUG) << "Comment line from HAKKM input file: " << comment;
     getline(flux_stream,comment);
     LOG("Flux", pDEBUG) << "Comment line from HAKKM input file: " << comment;
 
     iphi++;
     if(iphi == kGHnd3DNumPhiBins+1) {
       icostheta--;
       iphi=1;
     }
 
     LOG("Flux", pDEBUG)
        << "icostheta = " << icostheta << ", iphi = " << iphi << " / "
        << "costheta bins = " << kGHnd3DNumCosThetaBins << ", phi bins = " << kGHnd3DNumPhiBins;
     LOG("Flux", pDEBUG)
        << "The following set of (energy,flux) values corresponds to "
        << "costheta = [" << fCosThetaBins[icostheta-1] << ", " << fCosThetaBins[icostheta] << "]"
        << ", phi = [" << fPhiBins[iphi-1] << ", " << fPhiBins[iphi] << "] rad"
        << " or [" << r2d * fPhiBins[iphi-1] << ", " << r2d * fPhiBins[iphi] << "] deg) ";
 
     for(unsigned int i=0; i < kGHnd3DNumLogEvBins; i++) {
 
       int ienergy = i+1;
 
       double energy   = 0;
       double fnumu    = 0;
       double fnumubar = 0;
       double fnue     = 0;
       double fnuebar  = 0;
 
       flux_stream >> energy >> fnumu >> fnumubar >> fnue >> fnuebar;
 
       // fitting this easily into what is done for FLUKA, BGLRS where a
       // different file is specified for each neurtino species means that
       // the input file for HAKKM has to be read 4 times (at most).
       // However, this maintains the ability to switch off individual
       // components at source and generate interactions for some species only
 
       double flux = 0.;
       if(nu_pdg == kPdgNuMu    ) flux = fnumu;
       if(nu_pdg == kPdgAntiNuMu) flux = fnumubar;
       if(nu_pdg == kPdgNuE     ) flux = fnue;
       if(nu_pdg == kPdgAntiNuE ) flux = fnuebar;
       LOG("Flux", pDEBUG)
          << "Flux (nu_pdg = " << nu_pdg
          << "; Ev = " << energy << " GeV / bin used = ["
          << fEnergyBins[ienergy-1] << ", " << fEnergyBins[ienergy] << "] GeV"
          << ") = " << flux << " (m^2 sec sr GeV)^-1";
       if(flux > 0.) {
         histo->SetBinContent(ienergy,icostheta,iphi,flux);
       }
     }
     getline(flux_stream,comment);
     LOG("Flux", pDEBUG) << comment;
   }
   return true;
 }

void GHAKKMAtmoFlux::SetBinSizes ( void )

private

Definition at line 47 of file GHAKKMAtmoFlux.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::kGHnd3DCosThetaMax, genie::flux::kGHnd3DCosThetaMin, genie::flux::kGHnd3DNumCosThetaBins, genie::flux::kGHnd3DNumLogEvBins, genie::flux::kGHnd3DNumLogEvBinsPerDecade, genie::flux::kGHnd3DNumPhiBins, genie::flux::kGHnd3DPhiMax, genie::flux::kGHnd3DPhiMin, genie::constants::kPi, LOG, and pDEBUG.

 {
   //
   // cos(theta)
   //
 
   fCosThetaBins    = new double [kGHnd3DNumCosThetaBins + 1];
   fNumCosThetaBins = kGHnd3DNumCosThetaBins;
 
   double dcostheta =
       (kGHnd3DCosThetaMax - kGHnd3DCosThetaMin) /
       (double) kGHnd3DNumCosThetaBins;
 
   for(unsigned int i=0; i<= kGHnd3DNumCosThetaBins; i++) {
      fCosThetaBins[i] = kGHnd3DCosThetaMin + i * dcostheta;
      if(i != kGHnd3DNumCosThetaBins) {
        LOG("Flux", pDEBUG)
          << "HAKKM 3D flux: CosTheta bin " << i+1
          << ": lower edge = " << fCosThetaBins[i];
      } else {
        LOG("Flux", pDEBUG)
          << "HAKKM 3D flux: CosTheta bin " << kGHnd3DNumCosThetaBins
          << ": upper edge = " << fCosThetaBins[kGHnd3DNumCosThetaBins];
      }
   }
 
   //
   // phi
   //
 
   fPhiBins    = new double [kGHnd3DNumPhiBins + 1];
   fNumPhiBins = kGHnd3DNumPhiBins;
 
   double d2r = constants::kPi/180.;
 
   double dphi =
       d2r * (kGHnd3DPhiMax - kGHnd3DPhiMin) /
       (double) kGHnd3DNumPhiBins;
 
   for(unsigned int i=0; i<= kGHnd3DNumPhiBins; i++) {
      fPhiBins[i] = kGHnd3DPhiMin + i * dphi;
      if(i != kGHnd3DNumPhiBins) {
        LOG("Flux", pDEBUG)
          << "HAKKM 3D flux: Phi bin " << i+1
          << ": lower edge = " << fPhiBins[i];
      } else {
        LOG("Flux", pDEBUG)
          << "HAKKM 3D flux: Phi bin " << kGHnd3DNumPhiBins
          << ": upper edge = " << fPhiBins[kGHnd3DNumPhiBins];
      }
   }
 
   //
   // log(E)
   //
 
   // For each costheta,phi pair there are N logarithmically spaced
   // neutrino energy values (starting at 0.1 GeV with 20 values per decade
   // up to 10000 GeV) each with corresponding flux values.
   // To construct a flux histogram, use N+1 bins from 0 up to maximum
   // value. Assume that the flux value given for E=E0 is the flux value
   // at the bin that has E0 as its upper edge.
   //
   fEnergyBins    = new double [kGHnd3DNumLogEvBins + 1]; // bin edges
   fNumEnergyBins = kGHnd3DNumLogEvBins;
 
   double logEmax = TMath::Log10(1.);
   double logEmin = TMath::Log10(0.1);
   double dlogE =
       (logEmax - logEmin) /
       (double) kGHnd3DNumLogEvBinsPerDecade;
 
   fEnergyBins[0] = 0;
   for(unsigned int i=0; i < fNumEnergyBins; i++) {
      fEnergyBins[i+1] = TMath::Power(10., logEmin + i*dlogE);
      if(i < kGHnd3DNumLogEvBins) {
        LOG("Flux", pDEBUG)
           << "HAKKM 3D flux: Energy bin " << i+1
           << ": upper edge = " << fEnergyBins[i+1] << " GeV";
      }
   }
 
   fMaxEv = fEnergyBins[fNumEnergyBins];
   LOG("Flux", pDEBUG)
     << "HAKKM 3D flux: Maximum energy = " << fMaxEv;
 
 }

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