genie::flux::GAtmoFlux Class Reference

A base class for the FLUKA, BGLRS and ATMNC atmo. nu. flux drivers. The driver depends on data files provided by the atmospheric neutrino flux simulation authors in order to determine the angular and energy dependence for each neutrino species. The position of each flux neutrino [going towards a detector centered at (0,0,0)] is generated uniformly on a plane that is perpendicular to a sphere of radius Rl at the point that is determined by the generated neutrino direction (theta,phi). The size of the area of that plane, where flux neutrinos are generated, is determined by the transverse radius Rt. You can tweak Rl, Rt to match the size of your detector. Initially, neutrino coordinates are generated in a default detector coordinate system (Topocentric Horizontal Coordinate -THZ-): +z: Points towards the local zenith. +x: On same plane as local meridian, pointing south. +y: As needed to make a right-handed coordinate system. origin: detector centre Alternative user-defined topocentric systems can be defined by specifying the appropriate rotation from THZ. The driver allows minimum and maximum energy cuts. Also it provides the options to generate wither unweighted or weighted flux neutrinos (the latter giving smoother distributions at the tails). More...

#include <GAtmoFlux.h>

Inheritance diagram for genie::flux::GAtmoFlux:

Collaboration diagram for genie::flux::GAtmoFlux:

Public Member Functions
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 ()

Protected Member Functions
	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)
virtual bool	FillFluxHisto (int nu_pdg, string filename)=0
Protected Member Functions inherited from genie::GFluxI
	GFluxI ()

Protected Attributes
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 base class for the FLUKA, BGLRS and ATMNC atmo. nu. flux drivers. The driver depends on data files provided by the atmospheric neutrino flux simulation authors in order to determine the angular and energy dependence for each neutrino species. The position of each flux neutrino [going towards a detector centered at (0,0,0)] is generated uniformly on a plane that is perpendicular to a sphere of radius Rl at the point that is determined by the generated neutrino direction (theta,phi). The size of the area of that plane, where flux neutrinos are generated, is determined by the transverse radius Rt. You can tweak Rl, Rt to match the size of your detector. Initially, neutrino coordinates are generated in a default detector coordinate system (Topocentric Horizontal Coordinate -THZ-): +z: Points towards the local zenith. +x: On same plane as local meridian, pointing south. +y: As needed to make a right-handed coordinate system. origin: detector centre Alternative user-defined topocentric systems can be defined by specifying the appropriate rotation from THZ. The driver allows minimum and maximum energy cuts. Also it provides the options to generate wither unweighted or weighted flux neutrinos (the latter giving smoother distributions at the tails).

Author

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

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 60 of file GAtmoFlux.h.

Constructor & Destructor Documentation

GAtmoFlux::~GAtmoFlux ( )

virtual

Definition at line 39 of file GAtmoFlux.cxx.

{
  this->CleanUp();
}

GAtmoFlux::GAtmoFlux ( )

protected

Definition at line 34 of file GAtmoFlux.cxx.

{
  fInitialized = 0;
}

Member Function Documentation

void GAtmoFlux::AddAllFluxes ( void  )

protected

Definition at line 551 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

{
  LOG("Flux", pNOTICE)
       << "Computing combined flux & flux normalization factor";

  if(fTotalFluxHisto) delete fTotalFluxHisto;

  fTotalFluxHisto = this->CreateFluxHisto("sum", "combined flux" );

  map<int,TH3D*>::iterator it = fFluxHistoMap.begin();
  for( ; it != fFluxHistoMap.end(); ++it) {
    TH3D * flux_histogram = it->second;
    fTotalFluxHisto->Add(flux_histogram);
  }

  fTotalFluxHistoIntg = fTotalFluxHisto->Integral();
}

void GAtmoFlux::AddFluxFile	(	int	neutrino_pdg,
		string	filename
	)

Definition at line 394 of file GAtmoFlux.cxx.

References genie::pdg::IsAntiNeutrino(), genie::pdg::IsNeutrino(), LOG, pFATAL, and pWARN.

Referenced by GetFlux(), testGetTotalFlux(), and testGetTotalFluxInEnergyRange().

{
  // Check file exists
  std::ifstream f(filename.c_str());
  if (!f.good()) {
    LOG("Flux", pFATAL) << "Flux file does not exist "<<filename;
    exit(-1);
  }
  if ( pdg::IsNeutrino(nu_pdg) || pdg::IsAntiNeutrino(nu_pdg) ) {
    fFluxFlavour.push_back(nu_pdg);
    fFluxFile.push_back(filename);
  } else {
    LOG ("Flux", pWARN)
      << "Input particle code: " << nu_pdg << " not a neutrino!";
  }
}

void GAtmoFlux::AddFluxFile

(

string

filename

)

Definition at line 411 of file GAtmoFlux.cxx.

References genie::kPdgAntiNuE, genie::kPdgAntiNuMu, genie::kPdgNuE, genie::kPdgNuMu, LOG, and pFATAL.

{
// FLUKA and BGLRS provide one file per neutrino species.
// HAKKKM provides a single file for all nue,nuebar,numu,numubar.
// If no neutrino species is provided, assume that the file contains all 4
// but fit it into the franework developed for FLUKA and BGLRS,
// i.e. add the file 4 times

// Check file exists
  std::ifstream f(filename.c_str());
  if (!f.good()) {
    LOG("Flux", pFATAL) << "Flux file does not exist "<<filename;
    exit(-1);
  }

  fFluxFlavour.push_back(kPdgNuE);      fFluxFile.push_back(filename);
  fFluxFlavour.push_back(kPdgAntiNuE);  fFluxFile.push_back(filename);
  fFluxFlavour.push_back(kPdgNuMu);     fFluxFile.push_back(filename);
  fFluxFlavour.push_back(kPdgAntiNuMu); fFluxFile.push_back(filename);

}

void GAtmoFlux::CleanUp ( void  )

protected

Definition at line 336 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

{
  LOG("Flux", pNOTICE) << "Cleaning up...";

  map<int,TH3D*>::iterator rawiter = fRawFluxHistoMap.begin();
  for( ; rawiter != fRawFluxHistoMap.end(); ++rawiter) {
    TH3D * flux_histogram = rawiter->second;
    if(flux_histogram) {
       delete flux_histogram;
       flux_histogram = 0;
    }
  }
  fRawFluxHistoMap.clear();

  map<int,TH3D*>::iterator iter = fFluxHistoMap.begin();
  for( ; iter != fFluxHistoMap.end(); ++iter) {
    TH3D * flux_histogram = iter->second;
    if(flux_histogram) {
       delete flux_histogram;
       flux_histogram = 0;
    }
  }
  fFluxHistoMap.clear();

  if (fTotalFluxHisto) delete fTotalFluxHisto;
  if (fPdgCList) delete fPdgCList;

  if (fPhiBins     ) { delete[] fPhiBins     ; fPhiBins     =NULL; }
  if (fCosThetaBins) { delete[] fCosThetaBins; fCosThetaBins=NULL; }
  if (fEnergyBins  ) { delete[] fEnergyBins  ; fEnergyBins  =NULL; }

}

void GAtmoFlux::Clear ( Option_t *  opt )

virtual

reset state variables based on opt

Implements genie::GFluxI.

Definition at line 239 of file GAtmoFlux.cxx.

References LOG, and pERROR.

{
// Dummy clear method needed to conform to GFluxI interface
//
  LOG("Flux", pERROR) << "No clear method implemented for option:"<< opt;
}

double genie::flux::GAtmoFlux::CosTheta ( void  )

inline

Definition at line 81 of file GAtmoFlux.h.

References fgP4.

{ return -fgP4.Pz()/fgP4.Energy(); }

double genie::flux::GAtmoFlux::CosZenith ( void  )

inline

Definition at line 82 of file GAtmoFlux.h.

References fgP4.

{ return -fgP4.Pz()/fgP4.Energy(); }

TH3D * GAtmoFlux::CreateFluxHisto ( string  name, string  title  )

protected

Definition at line 569 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

{
  LOG("Flux", pNOTICE) << "Instantiating histogram: [" << name << "]";
  TH3D * hist = new TH3D(
     name.c_str(), title.c_str(),
     fNumEnergyBins,   fEnergyBins,
     fNumCosThetaBins, fCosThetaBins,
     fNumPhiBins,      fPhiBins);
  return hist;
}

TH3D * GAtmoFlux::CreateNormalisedFluxHisto ( TH3D *  hist )

protected

Definition at line 497 of file GAtmoFlux.cxx.

{
// return integrated flux

  // sanity check
  if(!hist) return 0;

  // make new histogram name
  TString histname = hist->GetName();
  histname.Append("_IntegratedFlux");

  // make new histogram
  TH3D* hist_intg = (TH3D*)(hist->Clone(histname.Data()));
  hist_intg->Reset();

  // integrate flux in each bin
  Double_t dN_dEdS = 0.0;
  Double_t dS = 0.0;
  Double_t dE = 0.0;
  Double_t dN = 0.0;

  for(Int_t e_bin = 1; e_bin <= hist->GetXaxis()->GetNbins(); e_bin++)
  {
    for(Int_t c_bin = 1; c_bin <= hist->GetYaxis()->GetNbins(); c_bin++)
    {
      for(Int_t p_bin = 1; p_bin <= hist->GetZaxis()->GetNbins(); p_bin++)
      {
         dN_dEdS = hist->GetBinContent(e_bin,c_bin,p_bin);

         dE = hist->GetXaxis()->GetBinUpEdge(e_bin)
            - hist->GetXaxis()->GetBinLowEdge(e_bin);

         dS = ( hist->GetZaxis()->GetBinUpEdge(p_bin)
              - hist->GetZaxis()->GetBinLowEdge(p_bin)  )
            * ( hist->GetYaxis()->GetBinUpEdge(c_bin)
              - hist->GetYaxis()->GetBinLowEdge(c_bin) );

         dN = dN_dEdS*dE*dS;

         hist_intg->SetBinContent(e_bin,c_bin,p_bin,dN);
      }
    }
  }

  return hist_intg;
}

virtual bool genie::flux::GAtmoFlux::End ( void  )

inlinevirtual

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

Implements genie::GFluxI.

Definition at line 73 of file GAtmoFlux.h.

{ return false;      }

double genie::flux::GAtmoFlux::Energy ( void  )

inline

Definition at line 80 of file GAtmoFlux.h.

References fgP4.

{ return fgP4.Energy(); }

double genie::flux::GAtmoFlux::Enu ( void  )

inline

Definition at line 79 of file GAtmoFlux.h.

References fgP4.

{ return fgP4.Energy(); }

virtual bool genie::flux::GAtmoFlux::FillFluxHisto ( int  nu_pdg, string  filename  )

protectedpure virtual

Implemented in genie::flux::GBGLRSAtmoFlux, genie::flux::GFLUKAAtmoFlux, and genie::flux::GHAKKMAtmoFlux.

virtual const PDGCodeList& genie::flux::GAtmoFlux::FluxParticles ( void  )

inlinevirtual

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

Implements genie::GFluxI.

Definition at line 66 of file GAtmoFlux.h.

References fPdgCList.

{ return *fPdgCList; }

void GAtmoFlux::ForceMaxEnergy

(

double

emax

)

Definition at line 233 of file GAtmoFlux.cxx.

Referenced by GetFlux(), testGetTotalFlux(), and testGetTotalFluxInEnergyRange().

{
  emax = TMath::Max(0., emax);
  fMaxEvCut = emax;
}

void GAtmoFlux::ForceMinEnergy

(

double

emin

)

Definition at line 227 of file GAtmoFlux.cxx.

Referenced by GetFlux(), testGetTotalFlux(), and testGetTotalFluxInEnergyRange().

{
  emin = TMath::Max(0., emin);
  fMinEvCut = emin;
}

bool GAtmoFlux::GenerateNext ( void  )

virtual

generate the next flux neutrino (return false in err)

Implements genie::GFluxI.

Definition at line 49 of file GAtmoFlux.cxx.

{
  while(1) {
     // Attempt to generate next flux neutrino
     bool nextok = this->GenerateNext_1try();
     if(!nextok) continue;

     // Check generated neutrino energy against max energy.
     // We may have to reject the current neutrino if a user-defined max
     // energy cut restricts the available range of energies.
     const TLorentzVector & p4 = this->Momentum();
     double E    = p4.Energy();
     double Emin = this->MinEnergy();
     double Emax = this->MaxEnergy();
     double wght = this->Weight();

     bool accept = (E<=Emax && E>=Emin && wght>0);
     if(accept) return true;
  }
  return false;
}

bool GAtmoFlux::GenerateNext_1try ( void  )

protected

Definition at line 71 of file GAtmoFlux.cxx.

References GetFlux(), genie::RandomGen::Instance(), genie::constants::kPi, LOG, genie::utils::print::P4AsShortString(), pINFO, genie::RandomGen::RndFlux(), and genie::utils::print::X4AsString().

{
  // Must have run intitialization
  assert(fInitialized);

  // Reset previously generated neutrino code / 4-p / 4-x
  this->ResetSelection();

  // Get a RandomGen instance
  RandomGen * rnd = RandomGen::Instance();

  // Generate (Ev, costheta, phi)
  double Ev       = 0.;
  double costheta = 0.;
  double phi      = 0;
  double weight   = 0;
  int    nu_pdg   = 0;

  if(fGenWeighted) {

     //
     // generate weighted flux
     //

     // generate events according to a power law spectrum,
     // then weight events by flux and inverse power law
     // (note: cannot use index alpha=1)
     double alpha = fSpectralIndex;

     double emin = TMath::Power(fEnergyBins[0],1.0-alpha);
     double emax = TMath::Power(fEnergyBins[fNumEnergyBins],1.0-alpha);
     Ev          = TMath::Power(emin+(emax-emin)*rnd->RndFlux().Rndm(),1.0/(1.0-alpha));
     costheta    = -1+2*rnd->RndFlux().Rndm();
     phi         = 2.*kPi* rnd->RndFlux().Rndm();

     unsigned int nnu = fPdgCList->size();
     unsigned int inu = rnd->RndFlux().Integer(nnu);
     nu_pdg   = (*fPdgCList)[inu];

     if(Ev < fEnergyBins[0]) {
        LOG("Flux", pINFO) << "E < Emin";
        return false;
     }
     double flux = this->GetFlux(nu_pdg, Ev, costheta, phi);
     if(flux<=0) {
        LOG("Flux", pINFO) << "Flux <= 0";
        return false;
     }
     weight = flux*TMath::Power(Ev,alpha);
  }
  else {

     //
     // generate nominal flux
     //

     Axis_t ax = 0, ay = 0, az = 0;
     fTotalFluxHisto->GetRandom3(ax, ay, az);
     Ev       = (double)ax;
     costheta = (double)ay;
     phi      = (double)az;
     nu_pdg   = this->SelectNeutrino(Ev, costheta, phi);
     weight   = 1.0;
  }

  // Compute etc trigonometric numbers
  double sintheta  = TMath::Sqrt(1-costheta*costheta);
  double cosphi    = TMath::Cos(phi);
  double sinphi    = TMath::Sin(phi);

  // Set the neutrino pdg code
  fgPdgC = nu_pdg;

  // Set the neutrino weight
  fWeight = weight;

  // Compute the neutrino momentum
  // The `-1' means it is directed towards the detector.
  double pz = -1.* Ev * costheta;
  double py = -1.* Ev * sintheta * sinphi;
  double px = -1.* Ev * sintheta * cosphi;

  // Default vertex is at the origin
  double z = 0.0;
  double y = 0.0;
  double x = 0.0;

  // Shift the neutrino position onto the flux generation surface.
  // The position is computed at the surface of a sphere with R=fRl
  // at the topocentric horizontal (THZ) coordinate system.
  if( fRl>0.0 ){
    z += fRl * costheta;
    y += fRl * sintheta * sinphi;
    x += fRl * sintheta * cosphi;
  }

  // Apply user-defined rotation from THZ -> user-defined topocentric
  // coordinate system.
  if( !fRotTHz2User.IsIdentity() )
  {
    TVector3 tx3(x, y, z );
    TVector3 tp3(px,py,pz);

    tx3 = fRotTHz2User * tx3;
    tp3 = fRotTHz2User * tp3;

    x  = tx3.X();
    y  = tx3.Y();
    z  = tx3.Z();
    px = tp3.X();
    py = tp3.Y();
    pz = tp3.Z();
  }

  // If the position is left as is, then all generated neutrinos
  // would point towards the origin.
  // Displace the position randomly on the surface that is
  // perpendicular to the selected point P(xo,yo,zo) on the sphere
  if( fRt>0.0 ){
    TVector3 vec(x,y,z);               // vector towards selected point
    TVector3 dvec1 = vec.Orthogonal(); // orthogonal vector
    TVector3 dvec2 = dvec1;            // second orthogonal vector
    dvec2.Rotate(-kPi/2.0,vec);        // rotate second vector by 90deg,
                                       // now forming a new orthogonal cartesian coordinate system
    double psi = 2.*kPi* rnd->RndFlux().Rndm(); // rndm angle [0,2pi]
    double random = rnd->RndFlux().Rndm();      // rndm number  [0,1]
    dvec1.SetMag(TMath::Sqrt(random)*fRt*TMath::Cos(psi));
    dvec2.SetMag(TMath::Sqrt(random)*fRt*TMath::Sin(psi));
    x += dvec1.X() + dvec2.X();
    y += dvec1.Y() + dvec2.Y();
    z += dvec1.Z() + dvec2.Z();
  }

  // Set the neutrino momentum and position 4-vectors with values
  // calculated at previous steps.
  fgP4.SetPxPyPzE(px, py, pz, Ev);
  fgX4.SetXYZT   (x,  y,  z,  0.);

  // Increment flux neutrino counter used for sample normalization purposes.
  fNNeutrinos++;

  // Report and exit
  LOG("Flux", pINFO)
       << "Generated neutrino: "
       << "\n pdg-code: " << fgPdgC
       << "\n p4: " << utils::print::P4AsShortString(&fgP4)
       << "\n x4: " << utils::print::X4AsString(&fgX4);

  return true;
}

void GAtmoFlux::GenerateWeighted ( bool  gen_weighted )

virtual

set whether to generate weighted or unweighted neutrinos

Implements genie::GFluxI.

Definition at line 246 of file GAtmoFlux.cxx.

{
  fGenWeighted = gen_weighted;
}

double GAtmoFlux::GetFlux

(

int

flavour

)

Definition at line 711 of file GAtmoFlux.cxx.

Referenced by testGetTotalFlux(), and testGetTotalFluxInEnergyRange().

{
  TH3D* flux_hist = this->GetFluxHistogram(flavour);
  if(!flux_hist) return 0.0;

  Double_t flux    = 0.0;
  Double_t dN_dEdS = 0.0;
  Double_t dS      = 0.0;
  Double_t dE      = 0.0;

  for(Int_t e_bin = 1; e_bin <= flux_hist->GetXaxis()->GetNbins(); e_bin++)
  {
    for(Int_t c_bin = 1; c_bin <= flux_hist->GetYaxis()->GetNbins(); c_bin++)
    {
      for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
      {
         dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);

         dE = flux_hist->GetXaxis()->GetBinUpEdge(e_bin)
            - flux_hist->GetXaxis()->GetBinLowEdge(e_bin);

         dS = ( flux_hist->GetZaxis()->GetBinUpEdge(p_bin)
              - flux_hist->GetZaxis()->GetBinLowEdge(p_bin) )
            * ( flux_hist->GetYaxis()->GetBinUpEdge(c_bin)
              - flux_hist->GetYaxis()->GetBinLowEdge(c_bin) );

         flux += dN_dEdS*dE*dS;
      }
    }
  }

  return flux;
}

double GAtmoFlux::GetFlux	(	int	flavour,
		double	energy
	)

Definition at line 745 of file GAtmoFlux.cxx.

{
  TH3D* flux_hist = this->GetFluxHistogram(flavour);
  if(!flux_hist) return 0.0;

  Double_t flux    = 0.0;
  Double_t dN_dEdS = 0.0;
  Double_t dS      = 0.0;

  Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);

  for(Int_t c_bin = 1; c_bin <= flux_hist->GetYaxis()->GetNbins(); c_bin++)
  {
    for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
    {
       dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);

       dS = ( flux_hist->GetZaxis()->GetBinUpEdge(p_bin)
              - flux_hist->GetZaxis()->GetBinLowEdge(p_bin) )
          * ( flux_hist->GetYaxis()->GetBinUpEdge(c_bin)
              - flux_hist->GetYaxis()->GetBinLowEdge(c_bin) );

       flux += dN_dEdS*dS;
    }
  }

  return flux;
}

double GAtmoFlux::GetFlux	(	int	flavour,
		double	energy,
		double	costh
	)

Definition at line 774 of file GAtmoFlux.cxx.

{
  TH3D* flux_hist = this->GetFluxHistogram(flavour);
  if(!flux_hist) return 0.0;

  Double_t flux    = 0.0;
  Double_t dN_dEdS = 0.0;
  Double_t dS      = 0.0;

  Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);
  Int_t c_bin = flux_hist->GetYaxis()->FindBin(costh);

  for( Int_t p_bin = 1; p_bin <= flux_hist->GetZaxis()->GetNbins(); p_bin++)
  {
     dN_dEdS = flux_hist->GetBinContent(e_bin,c_bin,p_bin);

     dS = ( flux_hist->GetZaxis()->GetBinUpEdge(p_bin)
          - flux_hist->GetZaxis()->GetBinLowEdge(p_bin) );

     flux += dN_dEdS*dS;
  }

  return flux;
}

double GAtmoFlux::GetFlux	(	int	flavour,
		double	energy,
		double	costh,
		double	phi
	)

Definition at line 799 of file GAtmoFlux.cxx.

{
  TH3D* flux_hist = this->GetFluxHistogram(flavour);
  if(!flux_hist) return 0.0;

  Int_t e_bin = flux_hist->GetXaxis()->FindBin(energy);
  Int_t c_bin = flux_hist->GetYaxis()->FindBin(costh);
  Int_t p_bin = flux_hist->GetZaxis()->FindBin(phi);

  Double_t flux = flux_hist->GetBinContent(e_bin,c_bin,p_bin);
  return flux;
}

TH3D * GAtmoFlux::GetFluxHistogram

(

int

flavour

)

Definition at line 628 of file GAtmoFlux.cxx.

{
  TH3D* histogram = 0;
  std::map<int,TH3D*>::iterator it = fRawFluxHistoMap.find(flavour);
  if(it != fRawFluxHistoMap.end())
  {
    histogram = it->second;
  }
  return histogram;
}

double GAtmoFlux::GetFluxSurfaceArea

(

void

)

Definition at line 378 of file GAtmoFlux.cxx.

References genie::constants::kPi.

Referenced by main().

{
  return kPi*pow(fRt,2);
}

double GAtmoFlux::GetLongitudinalRadius

(

void

)

Definition at line 383 of file GAtmoFlux.cxx.

{
  return fRl;
}

double GAtmoFlux::GetTotalFlux

(

void

)

Definition at line 640 of file GAtmoFlux.cxx.

References LOG, and pDEBUG.

Referenced by testGetTotalFlux(), and testGetTotalFluxInEnergyRange().

{
  double flux = 0.0;
  map<int,TH3D*>::iterator rawiter;

  rawiter = fRawFluxHistoMap.begin();
  for (; rawiter != fRawFluxHistoMap.end(); ++rawiter) {
    TH3D *h = rawiter->second;
    if (h) {
      flux += h->Integral("width");
      LOG("Flux", pDEBUG) << "Total flux for " << rawiter->first << " equals " << h->Integral("width") << ".";
    }
  }

  return flux;
}

double GAtmoFlux::GetTotalFluxInEnergyRange

(

void

)

Definition at line 659 of file GAtmoFlux.cxx.

References LOG, and pFATAL.

Referenced by main(), and testGetTotalFluxInEnergyRange().

{
  double flux = 0.0;
  map<int,TH3D*>::iterator rawiter;
  int e_min_bin, e_max_bin;
  double Emin, Emax;

  Emin = this->MinEnergy();
  Emax = this->MaxEnergy();

  if (Emax < Emin) {
      LOG("Flux", pFATAL) << "Emax = " << Emax << " is less than Emin = " << Emin;
      exit(-1);
  }

  rawiter = fRawFluxHistoMap.begin();
  for (; rawiter != fRawFluxHistoMap.end(); ++rawiter) {
    TH3D *h = rawiter->second;

    if (!h) continue;

    /* Get the bins containing `emin` and `emax`. */
    e_min_bin = h->GetXaxis()->FindBin(Emin);
    e_max_bin = h->GetXaxis()->FindBin(Emax);

    if (e_min_bin > h->GetXaxis()->GetNbins()) {
      /* If the minimum bin is past the end, continue. */
      continue;
    } else if (e_min_bin == e_max_bin) {
      /* If they both end up in the same bin, we just take the total bin
       * contents in that energy bin and multiply by the difference between
       * the energies. */
      flux += h->Integral(e_min_bin,e_min_bin,1,h->GetYaxis()->GetNbins(),1,h->GetZaxis()->GetNbins(),"width")*(Emax - Emin)/(h->GetXaxis()->GetBinUpEdge(e_min_bin)-h->GetXaxis()->GetBinLowEdge(e_min_bin));
    } else {
      /* First we calculate the integral within that bin from `emin` to the top
       * edge of the bin. */
      if (e_min_bin > 0)
          flux += h->Integral(e_min_bin,e_min_bin,1,h->GetYaxis()->GetNbins(),1,h->GetZaxis()->GetNbins(),"width")*(h->GetXaxis()->GetBinUpEdge(e_min_bin) - Emin)/(h->GetXaxis()->GetBinUpEdge(e_min_bin)-h->GetXaxis()->GetBinLowEdge(e_min_bin));
      /* Next, we calculate the integral for all the bins between the min and
       * max bin. */
      if (e_min_bin < h->GetXaxis()->GetNbins())
          flux += h->Integral(e_min_bin+1,e_max_bin-1,1,h->GetYaxis()->GetNbins(),1,h->GetZaxis()->GetNbins(),"width");
      /* Finally, we calculate the integral for the last bin. */
      if (e_max_bin <= h->GetXaxis()->GetNbins())
          flux += h->Integral(e_max_bin,e_max_bin,1,h->GetYaxis()->GetNbins(),1,h->GetZaxis()->GetNbins(),"width")*(Emax - h->GetXaxis()->GetBinLowEdge(e_max_bin))/(h->GetXaxis()->GetBinUpEdge(e_max_bin)-h->GetXaxis()->GetBinLowEdge(e_max_bin));
    }
  }

  return flux;
}

double GAtmoFlux::GetTransverseRadius

(

void

)

Definition at line 388 of file GAtmoFlux.cxx.

{
  return fRt;
}

virtual long int genie::flux::GAtmoFlux::Index ( void  )

inlinevirtual

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

Implements genie::GFluxI.

Definition at line 74 of file GAtmoFlux.h.

{ return -1;         }

void GAtmoFlux::Initialize ( void  )

protected

Definition at line 268 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

{
  LOG("Flux", pNOTICE) << "Initializing atmospheric flux driver";

  fMaxEv = -1;

  fNumPhiBins = 0;
  fNumCosThetaBins = 0;
  fNumEnergyBins = 0;
  fPhiBins = 0;
  fCosThetaBins = 0;
  fEnergyBins = 0;

  fTotalFluxHisto = 0;
  fTotalFluxHistoIntg = 0;

  bool allow_dup = false;
  fPdgCList = new PDGCodeList(allow_dup);

  // initializing flux TH3D histos [ flux = f(Ev,costheta,phi) ] & files
  fFluxFile.clear();
  fFluxHistoMap.clear();
  fTotalFluxHisto = 0;
  fTotalFluxHistoIntg = 0;

  // Default option is to generate unweighted flux neutrinos
  // (flux = f(E,costheta) will be used as PDFs)
  // User can enable option to generate weighted neutrinos
  // (neutrinos will be generated uniformly over costheta,
  // and using a power law function in neutrino energy.
  // The input flux = f(E,costheta) will be used for calculating a weight).
  // Using a weighted flux avoids statistical fluctuations at high energies.
  fSpectralIndex = 2.0;

  // weighting switched off by default
  this->GenerateWeighted(false);

  // Default: No min/max energy cut
  this->ForceMinEnergy(0.);
  this->ForceMaxEnergy(9999999999.);

  // Default radii
  fRl = 0.0;
  fRt = 0.0;

  // Default detector coord system: Topocentric Horizontal Coordinate system
  fRotTHz2User.SetToIdentity();

  // Reset `current' selected flux neutrino
  this->ResetSelection();

  // Reset number of neutrinos thrown so far
  fNNeutrinos = 0;

  // Done!
  fInitialized = 1;
}

bool GAtmoFlux::LoadFluxData

(

void

)

Definition at line 438 of file GAtmoFlux.cxx.

References genie::PDGLibrary::Find(), genie::PDGLibrary::Instance(), LOG, pERROR, and pNOTICE.

Referenced by GetFlux(), testGetTotalFlux(), and testGetTotalFluxInEnergyRange().

{
  LOG("Flux", pNOTICE)
        << "Loading atmospheric neutrino flux simulation data";

  fFluxHistoMap.clear();
  fPdgCList->clear();

  bool loading_status = true;

  for( unsigned int n=0; n<fFluxFlavour.size(); n++ ){
    int nu_pdg      = fFluxFlavour.at(n);
    string filename = fFluxFile.at(n);
    string pname = PDGLibrary::Instance()->Find(nu_pdg)->GetName();

    LOG("Flux", pNOTICE) << "Loading data for: " << pname;

    // create histogram
    TH3D* hist = 0;
    std::map<int,TH3D*>::iterator myMapEntry = fRawFluxHistoMap.find(nu_pdg);
    if( myMapEntry == fRawFluxHistoMap.end() ){
        hist = this->CreateFluxHisto(pname.c_str(), pname.c_str());
        fRawFluxHistoMap.insert( map<int,TH3D*>::value_type(nu_pdg,hist) );
    }
    // now let concrete instances to read the flux-specific data files
    // and fill the histogram
    bool loaded = this->FillFluxHisto(nu_pdg, filename);

    loading_status = loading_status && loaded;

    if (!loaded) {
        LOG("Flux", pERROR)
          << "Error loading atmospheric neutrino flux simulation data from " << filename;
        break;
    }
  }

  if(loading_status) {
    map<int,TH3D*>::iterator hist_iter = fRawFluxHistoMap.begin();
    for ( ; hist_iter != fRawFluxHistoMap.end(); ++hist_iter) {
      int   nu_pdg = hist_iter->first;
      TH3D* hist   = hist_iter->second;

      TH3D* hnorm = this->CreateNormalisedFluxHisto( hist );
      fFluxHistoMap.insert( map<int,TH3D*>::value_type(nu_pdg,hnorm) );
      fPdgCList->push_back(nu_pdg);
    }

    LOG("Flux", pNOTICE)
          << "Atmospheric neutrino flux simulation data loaded!";
    this->AddAllFluxes();
    return true;
  }

  LOG("Flux", pERROR)
    << "Error loading atmospheric neutrino flux simulation data";
  return false;
}

double GAtmoFlux::MaxEnergy ( void  )

virtual

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

Implements genie::GFluxI.

Definition at line 44 of file GAtmoFlux.cxx.

{
  return TMath::Min(fMaxEv, fMaxEvCut);
}

double genie::flux::GAtmoFlux::MinEnergy ( void  )

inlineprotected

Definition at line 120 of file GAtmoFlux.h.

References fMinEvCut.

{ return fMinEvCut; }

virtual const TLorentzVector& genie::flux::GAtmoFlux::Momentum ( void  )

inlinevirtual

returns the flux neutrino 4-momentum

Implements genie::GFluxI.

Definition at line 71 of file GAtmoFlux.h.

References fgP4.

{ return fgP4;       }

long int GAtmoFlux::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).

Definition at line 222 of file GAtmoFlux.cxx.

{
  return fNNeutrinos;
}

virtual int genie::flux::GAtmoFlux::PdgCode ( void  )

inlinevirtual

returns the flux neutrino pdg code

Implements genie::GFluxI.

Definition at line 69 of file GAtmoFlux.h.

References fgPdgC.

{ return fgPdgC;     }

virtual const TLorentzVector& genie::flux::GAtmoFlux::Position ( void  )

inlinevirtual

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

Implements genie::GFluxI.

Definition at line 72 of file GAtmoFlux.h.

References fgX4.

{ return fgX4;       }

void GAtmoFlux::ResetSelection ( void  )

protected

Definition at line 326 of file GAtmoFlux.cxx.

{
// initializing running neutrino pdg-code, 4-position, 4-momentum

  fgPdgC = 0;
  fgP4.SetPxPyPzE (0.,0.,0.,0.);
  fgX4.SetXYZT    (0.,0.,0.,0.);
  fWeight = 0;
}

int GAtmoFlux::SelectNeutrino ( double  Ev, double  costheta, double  phi  )

protected

Definition at line 580 of file GAtmoFlux.cxx.

References genie::RandomGen::Instance(), LOG, pDEBUG, pERROR, pINFO, and genie::RandomGen::RndFlux().

{
// Select a neutrino species at the input Ev and costheta given their
// relatve flux at this bin.
// Returns a neutrino PDG code

  unsigned int n = fPdgCList->size();
  double * flux = new double[n];

  unsigned int i=0;
  map<int,TH3D*>::iterator it = fFluxHistoMap.begin();
  for( ; it != fFluxHistoMap.end(); ++it) {
     TH3D * flux_histogram = it->second;
     int ibin = flux_histogram->FindBin(Ev,costheta,phi);
     flux[i]  = flux_histogram->GetBinContent(ibin);
     i++;
  }
  double flux_sum = 0;
  for(i=0; i<n; i++) {
     flux_sum  += flux[i];
     flux[i]    = flux_sum;
     LOG("Flux", pDEBUG)
       << "Sum{Flux(0->" << i <<")} = " << flux[i];
  }

  RandomGen * rnd = RandomGen::Instance();
  double R = flux_sum * rnd->RndFlux().Rndm();

  LOG("Flux", pDEBUG) << "R = " << R;

  for(i=0; i<n; i++) {
     if( R < flux[i] ) {
        delete [] flux;
        int selected_pdg = (*fPdgCList)[i];
        LOG("Flux", pINFO)
          << "Selected neutrino PDG code = " << selected_pdg;
        return selected_pdg;
     }
  }

  // shouldn't be here
  LOG("Flux", pERROR) << "Could not select a neutrino species!";
  assert(false);

  return -1;
}

void GAtmoFlux::SetRadii	(	double	Rlongitudinal,
		double	Rtransverse
	)

Definition at line 369 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

Referenced by GetFlux().

{
  LOG ("Flux", pNOTICE) << "Setting R[longitudinal] = " << Rlongitudinal;
  LOG ("Flux", pNOTICE) << "Setting R[transverse]   = " << Rtransverse;

  fRl = Rlongitudinal;
  fRt = Rtransverse;
}

void GAtmoFlux::SetSpectralIndex

(

double

index

)

Definition at line 251 of file GAtmoFlux.cxx.

References LOG, pNOTICE, and pWARN.

{
  if( index != 1.0 ){
    fSpectralIndex = index;
  }
  else {
    LOG("Flux", pWARN) << "Warning: cannot use a spectral index of unity";
  }

  LOG("Flux", pNOTICE) << "Using Spectral Index = " << index;
}

void GAtmoFlux::SetUserCoordSystem

(

TRotation &

rotation

)

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

Definition at line 263 of file GAtmoFlux.cxx.

Referenced by GetFlux().

{
  fRotTHz2User = rotation;
}

virtual double genie::flux::GAtmoFlux::Weight ( void  )

inlinevirtual

returns the flux neutrino weight (if any)

Implements genie::GFluxI.

Definition at line 70 of file GAtmoFlux.h.

References fWeight.

{ return fWeight;    }

void GAtmoFlux::ZeroFluxHisto ( TH3D *  hist )

protected

Definition at line 544 of file GAtmoFlux.cxx.

References LOG, and pDEBUG.

{
  LOG("Flux", pDEBUG) << "Forcing flux histogram contents to 0";
  if(!histo) return;
  histo->Reset();
}

Member Data Documentation

double* genie::flux::GAtmoFlux::fCosThetaBins

protected

cos(theta) bins in input flux data files

Definition at line 147 of file GAtmoFlux.h.

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

double* genie::flux::GAtmoFlux::fEnergyBins

protected

energy bins in input flux data files

Definition at line 148 of file GAtmoFlux.h.

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

vector<string> genie::flux::GAtmoFlux::fFluxFile

protected

input flux file for each neutrino species

Definition at line 157 of file GAtmoFlux.h.

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

protected

input flux file for each neutrino species

Definition at line 156 of file GAtmoFlux.h.

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

protected

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

Definition at line 154 of file GAtmoFlux.h.

bool genie::flux::GAtmoFlux::fGenWeighted

protected

generate a weighted or unweighted flux?

Definition at line 149 of file GAtmoFlux.h.

TLorentzVector genie::flux::GAtmoFlux::fgP4

protected

current generated nu 4-momentum

Definition at line 134 of file GAtmoFlux.h.

Referenced by CosTheta(), CosZenith(), Energy(), Enu(), and Momentum().

int genie::flux::GAtmoFlux::fgPdgC

protected

current generated nu pdg-code

Definition at line 133 of file GAtmoFlux.h.

Referenced by PdgCode().

TLorentzVector genie::flux::GAtmoFlux::fgX4

protected

current generated nu 4-position

Definition at line 135 of file GAtmoFlux.h.

Referenced by Position().

bool genie::flux::GAtmoFlux::fInitialized

protected

flag to check that initialization is run

Definition at line 151 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fMaxEv

protected

maximum energy (in input flux files)

Definition at line 131 of file GAtmoFlux.h.

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

double genie::flux::GAtmoFlux::fMaxEvCut

protected

user-defined cut: maximum energy

Definition at line 138 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fMinEvCut

protected

user-defined cut: minimum energy

Definition at line 139 of file GAtmoFlux.h.

Referenced by MinEnergy().

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

protected

number of flux neutrinos thrown so far

Definition at line 137 of file GAtmoFlux.h.

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

protected

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

Definition at line 144 of file GAtmoFlux.h.

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

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

protected

number of energy bins in input flux data files

Definition at line 145 of file GAtmoFlux.h.

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

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

protected

number of phi bins in input flux data files

Definition at line 143 of file GAtmoFlux.h.

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

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

protected

input list of neutrino pdg-codes

Definition at line 132 of file GAtmoFlux.h.

Referenced by FluxParticles().

double* genie::flux::GAtmoFlux::fPhiBins

protected

phi bins in input flux data files

Definition at line 146 of file GAtmoFlux.h.

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

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

protected

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

Definition at line 155 of file GAtmoFlux.h.

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

double genie::flux::GAtmoFlux::fRl

protected

defining flux neutrino generation surface: longitudinal radius

Definition at line 140 of file GAtmoFlux.h.

TRotation genie::flux::GAtmoFlux::fRotTHz2User

protected

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

Definition at line 142 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fRt

protected

defining flux neutrino generation surface: transverse radius

Definition at line 141 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fSpectralIndex

protected

power law function used for weighted flux

Definition at line 150 of file GAtmoFlux.h.

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

protected

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

Definition at line 152 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fTotalFluxHistoIntg

protected

fFluxSum2D integral

Definition at line 153 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fWeight

protected

current generated nu weight

Definition at line 136 of file GAtmoFlux.h.

Referenced by Weight().

---

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

• GAtmoFlux.h
• GAtmoFlux.cxx