Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
genie::flux::GAtmoFlux Class Referenceabstract

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 "/cvmfs/nova.opensciencegrid.org/externals/genie/v3_00_06_p01/Linux64bit+2.6-2.12-e17-debug/GENIE-Generator/src/Tools/Flux/GAtmoFlux.h"

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

Public Member Functions

virtual ~GAtmoFlux ()
 
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

 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 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...
 

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 <costas.andreopoulos stfc.ac.uk> University of Liverpool & STFC Rutherford Appleton Lab

January 26, 2008

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

Constructor & Destructor Documentation

GAtmoFlux::~GAtmoFlux ( )
virtual

Definition at line 82 of file GAtmoFlux.cxx.

83 {
84  this->CleanUp();
85 }
GAtmoFlux::GAtmoFlux ( )
protected

Definition at line 77 of file GAtmoFlux.cxx.

References slidt::fInitialized.

Referenced by CosZenith().

78 {
79  fInitialized = 0;
80 }
bool fInitialized
flag to check that initialization is run
Definition: GAtmoFlux.h:143

Member Function Documentation

void GAtmoFlux::AddAllFluxes ( void  )
protected

Definition at line 574 of file GAtmoFlux.cxx.

References it, LOG, and pNOTICE.

Referenced by 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 
)

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 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)

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  )
protected

Definition at line 379 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

Referenced by 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)
virtual

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 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  )
inline

Definition at line 82 of file GAtmoFlux.h.

References fgP4.

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

Definition at line 592 of file GAtmoFlux.cxx.

References analysePickle::hist, LOG, and pNOTICE.

Referenced by 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)
protected

Definition at line 520 of file GAtmoFlux.cxx.

References dE.

Referenced by 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  )
inlinevirtual

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  )
inline

Definition at line 81 of file GAtmoFlux.h.

References fgP4.

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

Definition at line 80 of file GAtmoFlux.h.

References fgP4.

80 { return fgP4.Energy(); }
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:126
virtual bool genie::flux::GAtmoFlux::FillFluxHisto ( int  nu_pdg,
string  filename 
)
protectedpure virtual
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 67 of file GAtmoFlux.h.

References fPdgCList, GenerateNext(), and MaxEnergy().

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

Definition at line 276 of file GAtmoFlux.cxx.

References emax.

Referenced by 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)

Definition at line 270 of file GAtmoFlux.cxx.

References emin.

Referenced by 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  )
virtual

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 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  )
protected

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 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)
virtual

set whether to generate weighted or unweighted neutrinos

Implements genie::GFluxI.

Definition at line 289 of file GAtmoFlux.cxx.

Referenced by GetFlux(), and Index().

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

Definition at line 662 of file GAtmoFlux.cxx.

References dE, and flux.

Referenced by 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 
)

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 
)

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 
)

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)

Definition at line 651 of file GAtmoFlux.cxx.

References RunSnowGlobes::histogram, and it.

Referenced by 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  )
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 75 of file GAtmoFlux.h.

References Clear(), GenerateWeighted(), and MECModelEnuComparisons::opt.

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

Definition at line 311 of file GAtmoFlux.cxx.

References slidt::fInitialized, LOG, and pNOTICE.

Referenced by 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  )

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 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  )
virtual

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 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  )
inlineprotected
virtual const TLorentzVector& genie::flux::GAtmoFlux::Momentum ( void  )
inlinevirtual

returns the flux neutrino 4-momentum

Implements genie::GFluxI.

Definition at line 72 of file GAtmoFlux.h.

References fgP4.

72 { return fgP4; }
TLorentzVector fgP4
current generated nu 4-momentum
Definition: GAtmoFlux.h:126
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 265 of file GAtmoFlux.cxx.

Referenced by 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  )
inlinevirtual

returns the flux neutrino pdg code

Implements genie::GFluxI.

Definition at line 70 of file GAtmoFlux.h.

References fgPdgC.

70 { return fgPdgC; }
int fgPdgC
current generated nu pdg-code
Definition: GAtmoFlux.h:125
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 73 of file GAtmoFlux.h.

References fgX4.

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

Definition at line 369 of file GAtmoFlux.cxx.

Referenced by 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 
)
protected

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 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 GAtmoFlux::SetRadii ( double  Rlongitudinal,
double  Rtransverse 
)

Definition at line 412 of file GAtmoFlux.cxx.

References LOG, and pNOTICE.

Referenced by 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)

Definition at line 294 of file GAtmoFlux.cxx.

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

Referenced by 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)

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

Definition at line 306 of file GAtmoFlux.cxx.

Referenced by 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  )
inlinevirtual

returns the flux neutrino weight (if any)

Implements genie::GFluxI.

Definition at line 71 of file GAtmoFlux.h.

References fWeight.

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

Definition at line 567 of file GAtmoFlux.cxx.

References LOG, and pDEBUG.

Referenced by 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
protected
double* genie::flux::GAtmoFlux::fEnergyBins
protected
vector<string> genie::flux::GAtmoFlux::fFluxFile
protected

input flux file for each neutrino species

Definition at line 149 of file GAtmoFlux.h.

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

input flux file for each neutrino species

Definition at line 148 of file GAtmoFlux.h.

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

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

Definition at line 146 of file GAtmoFlux.h.

bool genie::flux::GAtmoFlux::fGenWeighted
protected

generate a weighted or unweighted flux?

Definition at line 141 of file GAtmoFlux.h.

TLorentzVector genie::flux::GAtmoFlux::fgP4
protected

current generated nu 4-momentum

Definition at line 126 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 125 of file GAtmoFlux.h.

Referenced by PdgCode().

TLorentzVector genie::flux::GAtmoFlux::fgX4
protected

current generated nu 4-position

Definition at line 127 of file GAtmoFlux.h.

Referenced by Position().

bool genie::flux::GAtmoFlux::fInitialized
protected

flag to check that initialization is run

Definition at line 143 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fMaxEv
protected
double genie::flux::GAtmoFlux::fMaxEvCut
protected

user-defined cut: maximum energy

Definition at line 130 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fMinEvCut
protected

user-defined cut: minimum energy

Definition at line 131 of file GAtmoFlux.h.

Referenced by MinEnergy().

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

number of flux neutrinos thrown so far

Definition at line 129 of file GAtmoFlux.h.

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

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

Definition at line 136 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 137 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 135 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 124 of file GAtmoFlux.h.

Referenced by FluxParticles().

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

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

Definition at line 147 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 132 of file GAtmoFlux.h.

TRotation genie::flux::GAtmoFlux::fRotTHz2User
protected

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

Definition at line 134 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fRt
protected

defining flux neutrino generation surface: transverse radius

Definition at line 133 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fSpectralIndex
protected

power law function used for weighted flux

Definition at line 142 of file GAtmoFlux.h.

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

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

Definition at line 144 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fTotalFluxHistoIntg
protected

fFluxSum2D integral

Definition at line 145 of file GAtmoFlux.h.

double genie::flux::GAtmoFlux::fWeight
protected

current generated nu weight

Definition at line 128 of file GAtmoFlux.h.

Referenced by Weight().


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