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)
virtual double	MaxEnergy (void)
	declare the max flux neutrino energy that can be generated (for init. purposes)
virtual bool	GenerateNext (void)
	generate the next flux neutrino (return false in err)
virtual int	PdgCode (void)
	returns the flux neutrino pdg code
virtual double	Weight (void)
	returns the flux neutrino weight (if any)
virtual const TLorentzVector &	Momentum (void)
	returns the flux neutrino 4-momentum
virtual const TLorentzVector &	Position (void)
	returns the flux neutrino 4-position (note: expect SI rather than physical units)
virtual bool	End (void)
	true if no more flux nu's can be thrown (eg reaching end of beam sim ntuples)
virtual long int	Index (void)
	returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current entry number)
virtual void	Clear (Option_t *opt)
	reset state variables based on opt
virtual void	GenerateWeighted (bool gen_weighted)
	set whether to generate weighted or unweighted neutrinos
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).
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.
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)
PDGCodeList *	fPdgCList
	input list of neutrino pdg-codes
int	fgPdgC
	current generated nu pdg-code
TLorentzVector	fgP4
	current generated nu 4-momentum
TLorentzVector	fgX4
	current generated nu 4-position
double	fWeight
	current generated nu weight
long int	fNNeutrinos
	number of flux neutrinos thrown so far
double	fMaxEvCut
	user-defined cut: maximum energy
double	fMinEvCut
	user-defined cut: minimum energy
double	fRl
	defining flux neutrino generation surface: longitudinal radius
double	fRt
	defining flux neutrino generation surface: transverse radius
TRotation	fRotTHz2User
	coord. system rotation: THZ -> Topocentric user-defined
unsigned int	fNumPhiBins
	number of phi bins in input flux data files
unsigned int	fNumCosThetaBins
	number of cos(theta) bins in input flux data files
unsigned int	fNumEnergyBins
	number of energy bins in input flux data files
double *	fPhiBins
	phi bins in input flux data files
double *	fCosThetaBins
	cos(theta) bins in input flux data files
double *	fEnergyBins
	energy bins in input flux data files
bool	fGenWeighted
	generate a weighted or unweighted flux?
double	fSpectralIndex
	power law function used for weighted flux
bool	fInitialized
	flag to check that initialization is run
TH3D *	fTotalFluxHisto
	flux = f(Ev,cos8,phi) summed over neutrino species
double	fTotalFluxHistoIntg
	fFluxSum2D integral
map< int, TH3D * >	fFluxHistoMap
	flux = f(Ev,cos8,phi) for each neutrino species
map< int, TH3D * >	fRawFluxHistoMap
	flux = f(Ev,cos8,phi) for each neutrino species
vector< int >	fFluxFlavour
	input flux file for each neutrino species
vector< string >	fFluxFile
	input flux file for each neutrino species

Detailed Description

A driver for the HAKKM 3-D atmospheric neutrino flux (commonly known as the ‘Honda flux’)

References:\n 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 \at cern.ch> University of Liverpool

Created:\n July 9, 2015

License:\n Copyright (c) 2003-2025, 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::GHAKKMAtmoFlux ( )

Definition at line 32 of file GHAKKMAtmoFlux.cxx.

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

References genie::flux::GAtmoFlux::GAtmoFlux(), genie::flux::GAtmoFlux::Initialize(), LOG, pNOTICE, and SetBinSizes().

◆ ~GHAKKMAtmoFlux()

GHAKKMAtmoFlux::~GHAKKMAtmoFlux ( )

Definition at line 42 of file GHAKKMAtmoFlux.cxx.

{
 
}

Member Function Documentation

◆ FillFluxHisto()

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

privatevirtual

Implements genie::flux::GAtmoFlux.

Definition at line 135 of file GHAKKMAtmoFlux.cxx.

{
  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;
}

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.

◆ SetBinSizes()

void GHAKKMAtmoFlux::SetBinSizes ( void )

private

Definition at line 47 of file GHAKKMAtmoFlux.cxx.

{
  //
  // 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;
 
}

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.

Referenced by GHAKKMAtmoFlux().

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

Public Member Functions

Private Member Functions