genie::flux::GJPARCNuFlux Class Reference

A GENIE flux driver encapsulating the JPARC neutrino flux. It reads-in the official JPARC neutrino flux ntuples. More...

#include <GJPARCNuFlux.h>

Inheritance diagram for genie::flux::GJPARCNuFlux:

Public Member Functions
	GJPARCNuFlux ()
	~GJPARCNuFlux ()
const PDGCodeList &	FluxParticles (void)
	declare list of flux neutrinos that can be generated (for init. purposes) More...
double	MaxEnergy (void)
	declare the max flux neutrino energy that can be generated (for init. purposes) More...
bool	GenerateNext (void)
	generate the next flux neutrino (return false in err) More...
int	PdgCode (void)
	returns the flux neutrino pdg code More...
double	Weight (void)
	returns the flux neutrino weight (if any) More...
const TLorentzVector &	Momentum (void)
	returns the flux neutrino 4-momentum More...
const TLorentzVector &	Position (void)
	returns the flux neutrino 4-position (note: expect SI rather than physical units) More...
bool	End (void)
	true if no more flux nu's can be thrown (eg reaching end of beam sim ntuples) More...
long int	Index (void)
	returns corresponding index for current flux neutrino (e.g. for a flux ntuple returns the current entry number) More...
void	Clear (Option_t *opt)
	reset state variables based on opt More...
void	GenerateWeighted (bool gen_weighted=true)
	set whether to generate weighted or unweighted neutrinos More...
bool	LoadBeamSimData (string filename, string det_loc)
	load a jnubeam root flux ntuple More...
void	SetFluxParticles (const PDGCodeList &particles)
	specify list of flux neutrino species More...
void	SetMaxEnergy (double Ev)
	specify maximum flx neutrino energy More...
void	SetFilePOT (double pot)
	flux file norm is in /N POT/det [ND] or /N POT/cm^2 [FD]. Specify N (typically 1E+21) More...
void	SetUpstreamZ (double z0)
	set flux neutrino initial z position (upstream of the detector) More...
void	SetNumOfCycles (int n)
	set how many times to cycle through the ntuple (default: 1 / n=0 means 'infinite') More...
void	DisableOffset (void)
	switch off random offset, must be called before LoadBeamSimData to have any effect More...
void	RandomOffset (void)
	choose a random offset as starting entry in flux ntuple More...
double	POT_1cycle (void)
	flux POT per cycle More...
double	POT_curravg (void)
	current average POT More...
long int	NFluxNeutrinos (void) const
	number of flux neutrinos looped so far More...
double	SumWeight (void) const
	intergated weight for flux neutrinos looped so far More...
const GJPARCNuFluxPassThroughInfo &	PassThroughInfo (void)
	GJPARCNuFluxPassThroughInfo. More...
Public Member Functions inherited from genie::GFluxI
virtual	~GFluxI ()

Private Member Functions
bool	GenerateNext_weighted (void)
void	Initialize (void)
void	SetDefaults (void)
void	CleanUp (void)
void	ResetCurrent (void)
int	DLocName2Id (string name)

Private Attributes
double	fMaxEv
	maximum energy More...
PDGCodeList *	fPdgCList
	list of neutrino pdg-codes More...
int	fgPdgC
	running generated nu pdg-code More...
TLorentzVector	fgP4
	running generated nu 4-momentum More...
TLorentzVector	fgX4
	running generated nu 4-position More...
TFile *	fNuFluxFile
	input flux file More...
TTree *	fNuFluxTree
	input flux ntuple More...
TChain *	fNuFluxChain
	input flux ntuple More...
TTree *	fNuFluxSumTree
	input summary ntuple for flux file. This tree is only present for later flux versions > 10a More...
TChain *	fNuFluxSumChain
	input summary ntuple for flux file. This tree is only present for later flux versions > 10a More...
bool	fNuFluxUsingTree
	are we using a TTree or a TChain to view the input flux file? More...
string	fDetLoc
	input detector location ('sk','nd1','nd2',...) More...
int	fDetLocId
	input detector location id (fDetLoc -> jnubeam idfd) More...
int	fNDetLocIdFound
	per cycle keep track of the number of fDetLocId are found - if this is zero will exit job More...
bool	fIsFDLoc
	input location is a 'far' detector location? More...
bool	fIsNDLoc
	input location is a 'near' detector location? More...
long int	fNEntries
	number of flux ntuple entries More...
long int	fIEntry
	current flux ntuple entry More...
long int	fEntriesThisCycle
	keep track of number of entries used so far for this cycle More...
long int	fOffset
	start looping at entry fOffset More...
double	fNorm
	current flux ntuple normalisation More...
double	fMaxWeight
	max flux neutrino weight in input file for the specified detector location More...
double	fFilePOT
	file POT normalization, typically 1E+21 More...
double	fZ0
	configurable starting z position for each flux neutrino (in detector coord system) More...
int	fNCycles
	how many times to cycle through the flux ntuple More...
int	fICycle
	current cycle More...
double	fSumWeight
	sum of weights for neutrinos thrown so far More...
long int	fNNeutrinos
	number of flux neutrinos thrown so far More...
double	fSumWeightTot1c
	total sum of weights for neutrinos to be thrown / cycle More...
long int	fNNeutrinosTot1c
	total number of flux neutrinos to be thrown / cycle More...
bool	fGenerateWeighted
	generate weighted/deweighted flux neutrinos (default is false) More...
bool	fUseRandomOffset
	whether set random starting point when looping over flux ntuples More...
bool	fLoadedNeutrino
	set to true when GenerateNext_weighted has been called successfully More...
GJPARCNuFluxPassThroughInfo *	fPassThroughInfo

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

Detailed Description

A GENIE flux driver encapsulating the JPARC neutrino flux. It reads-in the official JPARC neutrino flux ntuples.

See http://jnusrv01.kek.jp/internal/t2k/nubeam/flux/ (Note: T2K internal)

Author

Costas Andreopoulos <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

Feb 04, 2008

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

Definition at line 50 of file GJPARCNuFlux.h.

Constructor & Destructor Documentation

genie::flux::GJPARCNuFlux::GJPARCNuFlux

(

)

GJPARCNuFlux::~GJPARCNuFlux

(

)

Definition at line 50 of file GJPARCNuFlux.cxx.

{
  this->CleanUp();
}

Member Function Documentation

void GJPARCNuFlux::CleanUp ( void  )

private

Definition at line 909 of file GJPARCNuFlux.cxx.

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

  if (fPdgCList)        delete fPdgCList;
  if (fPassThroughInfo) delete fPassThroughInfo;

  if (fNuFluxFile) {
        fNuFluxFile->Close();
        delete fNuFluxFile;
  }
}

void GJPARCNuFlux::Clear ( Option_t *  opt )

virtual

reset state variables based on opt

Implements genie::GFluxI.

Definition at line 806 of file GJPARCNuFlux.cxx.

{
// If opt = "CycleHistory" then:
// Reset all counters and state variables from any previous cycles. This
// should be called if, for instance, before event generation this flux
// driver had been used for pre-calculating interaction probabilities. Does
// not reset initial state variables such as flux particle types, POTs etc...
//
  if(std::strcmp(opt, "CycleHistory") == 0){
    // Reset so that generate de-weighted events
    this->GenerateWeighted(false);
    // Reset cycle counters
    fICycle          = 0;
    fSumWeight       = 0;
    fNNeutrinos      = 0;
    fIEntry          = fOffset;
    fEntriesThisCycle = 0;
  }
}

void genie::flux::GJPARCNuFlux::DisableOffset ( void  )

inline

switch off random offset, must be called before LoadBeamSimData to have any effect

Definition at line 83 of file GJPARCNuFlux.h.

int GJPARCNuFlux::DLocName2Id ( string  name )

private

Definition at line 922 of file GJPARCNuFlux.cxx.

{
// detector location: name -> int id
// sk  -> -1
// nd1 -> +1
// nd2 -> +2
// ...

  if(name == "sk"  ) return  -1;

  TString temp;
  for (int i=1; i<51; i++) {
    temp.Form("nd%d",i);
    if(name == temp.Data()) return i;
  }

  return 0;
}

bool genie::flux::GJPARCNuFlux::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 66 of file GJPARCNuFlux.h.

                                              { return  fEntriesThisCycle >= fNEntries
                                                     && fICycle == fNCycles && fNCycles > 0;   }

const PDGCodeList& genie::flux::GJPARCNuFlux::FluxParticles ( void  )

inlinevirtual

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

Implements genie::GFluxI.

Definition at line 59 of file GJPARCNuFlux.h.

{ return *fPdgCList;             }

bool GJPARCNuFlux::GenerateNext ( void  )

virtual

generate the next flux neutrino (return false in err)

Implements genie::GFluxI.

Definition at line 55 of file GJPARCNuFlux.cxx.

{
// Get next (unweighted) flux ntuple entry on the specified detector location
//
  RandomGen * rnd = RandomGen::Instance();
  while(1) {
     // Check for end of flux ntuple
     bool end = this->End();
     if(end) return false;

     // Get next weighted flux ntuple entry
     bool nextok = this->GenerateNext_weighted();
     if(!nextok) continue;

     if(fNCycles==0) {
       LOG("Flux", pNOTICE)
          << "Got flux entry: " << this->Index()
          << " - Cycle: "<< fICycle << "/ infinite";
     } else {
       LOG("Flux", pNOTICE)
          << "Got flux entry: "<< this->Index()
          << " - Cycle: "<< fICycle << "/"<< fNCycles;
     }

     // If de-weighting get fractional weight & decide whether to accept curr flux neutrino
     double f = 1.0;
     if(fGenerateWeighted == false) f = this->Weight();
     LOG("Flux", pNOTICE)
        << "Curr flux neutrino fractional weight = " << f;
     if(f > (1.+controls::kASmallNum)) {
       LOG("Flux", pERROR)
           << "** Fractional weight = " << f << " > 1 !!";
     }
     double r = (f < 1.) ? rnd->RndFlux().Rndm() : 0;
     bool accept = (r<f);
     if(accept) {
       return true;
     }

     LOG("Flux", pNOTICE)
       << "** Rejecting current flux neutrino based on the flux weight only";
  }
  return false;
}

bool GJPARCNuFlux::GenerateNext_weighted ( void  )

private

Definition at line 100 of file GJPARCNuFlux.cxx.

{
// Get next (weighted) flux ntuple entry on the specified detector location
//

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

  // Check whether a jnubeam flux ntuple has been loaded
  if( (!fNuFluxTree && fNuFluxUsingTree) || (!fNuFluxChain && !fNuFluxUsingTree) ) {
     LOG("Flux", pFATAL)
          << "The flux driver has not been properly configured";
     //return false; // don't do this - creates an infinite loop!
     exit(1);
  }

  // Read next flux ntuple entry. Use fEntriesThisCycle to keep track of when
  // in new cycle as fIEntry can now have an offset
  if(fEntriesThisCycle >= fNEntries) {
     // Exit if have not found neutrino at specified location for whole cycle
     if(fNDetLocIdFound == 0){
       LOG("Flux", pFATAL)
         << "The input jnubeam flux ntuple contains no entries for detector id "
         << fDetLocId << ". Terminating job!";
       exit(1);
     }
     fNDetLocIdFound = 0; // reset the counter
     fICycle++;
     fIEntry=fOffset;
     fEntriesThisCycle = 0;
     // Run out of entries @ the current cycle.
     // Check whether more (or infinite) number of cycles is requested
     if(fICycle >= fNCycles && fNCycles != 0){
        LOG("Flux", pWARN)
            << "No more entries in input flux neutrino ntuple";
        return false;
     }
  }

  // In addition to getting info to generate event the following also
  // updates pass-through info (= info on the flux neutrino parent particle
  // that may be stored at an extra branch of the output event tree -alongside
  // with the generated event branch- for use further upstream in the t2k
  // analysis chain -eg for beam reweighting etc-)
  bool found_entry;
  if (fNuFluxUsingTree)
    found_entry = fNuFluxTree->GetEntry(fIEntry) > 0;
  else
    found_entry = fNuFluxChain->GetEntry(fIEntry) > 0;
  assert(found_entry);
  fLoadedNeutrino = true;
  fEntriesThisCycle++;
  fIEntry = (fIEntry+1) % fNEntries;

  if (fNuFluxUsingTree) {
    if(fNuFluxSumTree) fNuFluxSumTree->GetEntry(0); // get entry 0 as only 1 entry in tree
  }
  else {
    // get entry corresponding to current tree number in the chain, as only 1 entry in each tree
    if(fNuFluxSumChain) fNuFluxSumChain->GetEntry(fNuFluxChain->GetTreeNumber());
  }

  // check for negative flux weights
  if(fPassThroughInfo->norm + controls::kASmallNum < 0.0){
    LOG("Flux", pERROR) << "Negative flux weight! Will set weight to 0.0";
    fPassThroughInfo->norm  = 0.0;
  }
  // remember to update fNorm as no longer set to fNuFluxTree branch address
  fNorm = (double) fPassThroughInfo->norm;

  // for 'near detector' flux ntuples make sure that the current entry
  // corresponds to a flux neutrino at the specified detector location
  if(fIsNDLoc           /* nd */  &&
     fDetLocId!=fPassThroughInfo->idfd /* doesn't match specified detector location*/) {
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
        LOG("Flux", pNOTICE)
          << "Current flux neutrino not at specified detector location";
#endif
        return false;
  }


  //
  // Handling neutrinos at specified detector location
  //

  // count the number of times we have neutrinos at specified detector location
  fNDetLocIdFound += 1;


  // update the sum of weights & number of neutrinos
  fSumWeight += this->Weight() * fMaxWeight; // Weight returns fNorm/fMaxWeight
  fNNeutrinos++;

  // Convert the current parent paticle decay mode into a neutrino pdg code
  // See:
  // http://jnusrv01.kek.jp/internal/t2k/nubeam/flux/nemode.h
  //  mode    description
  //  11      numu from pi+
  //  12      numu from K+
  //  13      numu from mu-
  //  21      numu_bar from pi-
  //  22      numu_bar from K-
  //  23      numu_bar from mu+
  //  31      nue from K+ (Ke3)
  //  32      nue from K0L(Ke3)
  //  33      nue from Mu+
  //  41      nue_bar from K- (Ke3)
  //  42      nue_bar from K0L(Ke3)
  //  43      nue_bar from Mu-
  // Since JNuBeam flux version >= 10a the following modes also expected
  //  14      numu from K+(3)
  //  15      numu from K0(3)
  //  24      numu_bar from K-(3)
  //  25      numu_bar from K0(3)
  //  34      nue from pi+
  //  44      nue_bar from pi-
  // In general expect more modes following the rule:
  //  11->19 --> numu
  //  21->29 --> numu_bar
  //  31->39 --> nue
  //  41->49 --> nuebar
  // This is based on example given at:
  //  http://jnusrv01.kek.jp/internal/t2k/nubeam/flux/efill.kumac

  if(fPassThroughInfo->mode >= 11 && fPassThroughInfo->mode <= 19) fgPdgC = kPdgNuMu;
  else if(fPassThroughInfo->mode >= 21 && fPassThroughInfo->mode <= 29) fgPdgC = kPdgAntiNuMu;
  else if(fPassThroughInfo->mode >= 31 && fPassThroughInfo->mode <= 39) fgPdgC = kPdgNuE;
  else if(fPassThroughInfo->mode >= 41 && fPassThroughInfo->mode <= 49) fgPdgC = kPdgAntiNuE;
  else {
    // If here then trying to process a neutrino from an unknown decay mode.
    // Rather than just skipping this flux neutrino the job is aborted to avoid
    // unphysical results.
    LOG("Flux", pFATAL) << "Unexpected decay mode: "<< fPassThroughInfo->mode <<
                           "  --> unable to infer neutrino pdg! Aborting job!";
    exit(1);
  }

  // Check neutrino pdg against declared list of neutrino species declared
  // by the current instance of the JPARC neutrino flux driver.
  // No undeclared neutrino species will be accepted at this point as GENIE
  // has already been configured to handle the specified list.
  // Make sure that the appropriate list of flux neutrino species was set at
  // initialization via GJPARCNuFlux::SetFluxParticles(const PDGCodeList &)

  if( ! fPdgCList->ExistsInPDGCodeList(fgPdgC) ) {
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
     LOG("Flux", pNOTICE)
       << "Current flux neutrino "
       << "not at the list of neutrinos to be considered at this job.";
#endif
     return false;
  }

  // Check current neutrino energy against the maximum flux neutrino energy declared
  // by the current instance of the JPARC neutrino flux driver.
  // No flux neutrino exceeding that maximum energy will be accepted at this point as
  // that maximum energy has already been used for normalizing the interaction probabilities.
  // Make sure that the appropriate maximum flux neutrino energy was set at
  // initialization via GJPARCNuFlux::SetMaxEnergy(double Ev)

  if(fPassThroughInfo->Enu > fMaxEv) {
     LOG("Flux", pWARN)
          << "Flux neutrino energy exceeds declared maximum neutrino energy";
     LOG("Flux", pWARN)
          << "Ev = " << fPassThroughInfo->Enu << "(> Ev{max} = " << fMaxEv << ")";
  }

  // Set the current flux neutrino 4-momentum & 4-position

  double pxnu = fPassThroughInfo->Enu * fPassThroughInfo->nnu[0];
  double pynu = fPassThroughInfo->Enu * fPassThroughInfo->nnu[1];
  double pznu = fPassThroughInfo->Enu * fPassThroughInfo->nnu[2];
  double Enu  = fPassThroughInfo->Enu;
  fgP4.SetPxPyPzE (pxnu, pynu, pznu, Enu);

  if(fIsNDLoc) {
    double cm2m = units::cm / units::m;
    double xnu  = cm2m * fPassThroughInfo->xnu;
    double ynu  = cm2m * fPassThroughInfo->ynu;
    double znu  = 0;

    // projected 4-position (from z=0) back to a configurable plane
    // position (fZ0) upstream of the detector face.
    xnu += (fZ0/fPassThroughInfo->nnu[2])*fPassThroughInfo->nnu[0];
    ynu += (fZ0/fPassThroughInfo->nnu[2])*fPassThroughInfo->nnu[1];
    znu = fZ0;

    fgX4.SetXYZT (xnu,  ynu,  znu,  0.);
  } else {
    fgX4.SetXYZT (0.,0.,0.,0.);

  }

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("Flux", pINFO)
        << "Generated neutrino: "
        << "\n pdg-code: " << fgPdgC
        << "\n p4: " << utils::print::P4AsShortString(&fgP4)
        << "\n x4: " << utils::print::X4AsString(&fgX4);
#endif
  // Update flux pass through info not set as branch addresses of flux ntuples
  if(fNuFluxUsingTree)
    fPassThroughInfo->fluxentry = this->Index();
  else
    fPassThroughInfo->fluxentry = fNuFluxChain->GetTree()->GetReadEntry();

  std::string filename;
  if(fNuFluxUsingTree)
    filename = fNuFluxFile->GetName();
  else
    filename = fNuFluxChain->GetFile()->GetName();
  std::string::size_type start_pos = filename.rfind("/");
  if (start_pos == std::string::npos) start_pos = 0; else ++start_pos;
  std::string basename(filename,start_pos);
  if(fNuFluxUsingTree)
    fPassThroughInfo->fluxfilename = basename + ":" + fNuFluxTree->GetName();
  else
    fPassThroughInfo->fluxfilename = basename + ":" + fNuFluxChain->GetName();
  return true;
}

void GJPARCNuFlux::GenerateWeighted ( bool  gen_weighted = true )

virtual

set whether to generate weighted or unweighted neutrinos

Implements genie::GFluxI.

Definition at line 781 of file GJPARCNuFlux.cxx.

{
// Ignore the flux weights when getting next flux neutrino. Always set to
// true for unweighted event generation but need option to switch off when
// using pre-calculated interaction probabilities for each flux ray to speed
// up event generation.
  fGenerateWeighted = gen_weighted;
}

long int GJPARCNuFlux::Index ( void  )

virtual

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

Implements genie::GFluxI.

Definition at line 363 of file GJPARCNuFlux.cxx.

{
// Return the current flux entry index. If GenerateNext has not yet been
// called then return -1.
//
  if(fLoadedNeutrino){
    // subtract 1 as fIEntry was incremented since call to TTree::GetEntry
    // and deal with special case where fIEntry-1 is last entry in cycle
    return fIEntry == 0 ? fNEntries - 1 : fIEntry-1;
  }
  // return -1 if no neutrino loaded since last call to this->ResetCurrent()
  return -1;
}

void GJPARCNuFlux::Initialize ( void  )

private

Definition at line 826 of file GJPARCNuFlux.cxx.

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

  fMaxEv           = 0;
  fPdgCList        = new PDGCodeList;
  fPassThroughInfo = new GJPARCNuFluxPassThroughInfo;

  fNuFluxFile      = 0;
  fNuFluxTree      = 0;
  fNuFluxChain     = 0;
  fNuFluxSumTree   = 0;
  fNuFluxSumChain  = 0;
  fNuFluxUsingTree = true;
  fDetLoc          = "";
  fDetLocId        = 0;
  fNDetLocIdFound  = 0;
  fIsFDLoc         = false;
  fIsNDLoc         = false;

  fNEntries        = 0;
  fIEntry          = 0;
  fEntriesThisCycle= 0;
  fOffset          = 0;
  fNorm            = 0.;
  fMaxWeight       = 0;
  fFilePOT         = 0;
  fZ0              = 0;
  fNCycles         = 0;
  fICycle          = 0;
  fSumWeight       = 0;
  fNNeutrinos      = 0;
  fSumWeightTot1c  = 0;
  fNNeutrinosTot1c = 0;
  fGenerateWeighted= false;
  fUseRandomOffset = true;
  fLoadedNeutrino  = false;

  this->SetDefaults();
  this->ResetCurrent();
}

bool GJPARCNuFlux::LoadBeamSimData	(	string	filename,
		string	det_loc
	)

load a jnubeam root flux ntuple

Definition at line 377 of file GJPARCNuFlux.cxx.

{
// Loads in a jnubeam beam simulation root file (converted from hbook format)
// into the GJPARCNuFlux driver.
// The detector location can be any of:
//  "sk","nd1" (<-2km),"nd5" (<-nd280),...,"nd10"

  LOG("Flux", pNOTICE)
        << "Loading jnubeam flux tree from ROOT file: " << filename;
  LOG("Flux", pNOTICE)
        << "Detector location: " << detector_location;

  // Check to see if its a single flux file (/dir/root_filename.0.root)
  // or a sequence of files to be tchained (e.g. /dir/root_filename@0@100)
  fNuFluxUsingTree = true;
  if (filename.find('@') != string::npos) { fNuFluxUsingTree = false; }

  vector<string> filenamev = utils::str::Split(filename,"@");
  string fileroot = "";
  int firstfile = -1, lastfile = -1;

  if (!fNuFluxUsingTree) {
    if (filenamev.size() != 3) {
      LOG("Flux", pFATAL)
        << "Flux filename should be specfied as either:\n"
        << "\t For a single input file:  /dir/root_filename.#.root\n"
        << "\t For multiple input files: /dir/root_filename@#@#";
      exit(1);
    }
    fileroot  = filenamev[0];
    firstfile = atoi(filenamev[1].c_str());
    lastfile  = atoi(filenamev[2].c_str());
    LOG("Flux", pNOTICE)
      << "Chaining beam simulation output files with stem: " << fileroot
      << " and run numbers in the range: [" << firstfile << ", " << firstfile << "]";
  }

  if (fNuFluxUsingTree) {
    bool is_accessible = ! (gSystem->AccessPathName( filename.c_str() ));
    if (!is_accessible) {
      LOG("Flux", pFATAL)
        << "The input jnubeam flux file doesn't exist! Initialization failed!";
      exit(1);
    }
  }
  else {
    bool please_exit = false;
    for (int i = firstfile; i < lastfile+1; i++) {
      bool is_accessible = ! (gSystem->AccessPathName( Form("%s.%i.root",fileroot.c_str(),i) ));
      if (!is_accessible) {
        LOG("Flux", pFATAL)
          << "The input jnubeam flux file " << Form("%s.%i.root",fileroot.c_str(),i)
          << "doesn't exist! Initialization failed!";
        please_exit = true;
      }
    }
    if(please_exit)
      exit(1);
  }

  fDetLoc   = detector_location;
  fDetLocId = this->DLocName2Id(fDetLoc);

  if(fDetLocId == 0) {
    LOG("Flux", pERROR)
         << " ** Unknown input detector location: " << fDetLoc;
    return false;
  }

  fIsFDLoc = (fDetLocId==-1);
  fIsNDLoc = (fDetLocId>0);

  if (!fNuFluxUsingTree) {
    string ntuple_name = (fIsNDLoc) ? "h3002" : "h2000";
    fNuFluxChain = new TChain(ntuple_name.c_str());
    int result = fNuFluxChain->Add( Form("%s.%i.root",fileroot.c_str(),firstfile), 0);
    if (result != 1 && fIsNDLoc) {
      LOG("Flux", pINFO)
        << "Could not find tree h3002 in file " << Form("%s.%i.root",fileroot.c_str(),firstfile)
        << " Trying tree h3001";
      delete fNuFluxChain;
      ntuple_name = "h3001";
      fNuFluxChain = new TChain(ntuple_name.c_str());
      result = fNuFluxChain->Add( Form("%s.%i.root",fileroot.c_str(),firstfile), 0);
    }
    if (result != 1) {
      LOG("Flux", pERROR)
        << "** Couldn't get flux tree: " << ntuple_name;
      return false;
    }

    for (int i = firstfile+1; i < lastfile+1; i++) {
      result = fNuFluxChain->Add( Form("%s.%i.root",fileroot.c_str(),i), 0 );
      if (result == 0)
        LOG("Flux", pERROR)
          << "** Couldn't get flux tree " << ntuple_name << " in file " << Form("%s.%i.root",fileroot.c_str(),i);
    fNEntries = fNuFluxChain->GetEntries();
    }
  }

  else {
    fNuFluxFile = new TFile(filename.c_str(), "read");
    if(fNuFluxFile) {
      // nd treename can be h3002 or h3001 depending on fluxfile version
      string ntuple_name = (fIsNDLoc) ? "h3002" : "h2000";
      fNuFluxTree = (TTree*) fNuFluxFile->Get(ntuple_name.c_str());
      if(!fNuFluxTree && fIsNDLoc){
        ntuple_name = "h3001";
        fNuFluxTree = (TTree*) fNuFluxFile->Get(ntuple_name.c_str());
      }
      LOG("Flux", pINFO)
        << "Getting flux tree: " << ntuple_name;
      if(!fNuFluxTree) {
        LOG("Flux", pERROR)
          << "** Couldn't get flux tree: " << ntuple_name;
        return false;
      }
    } else {
      LOG("Flux", pERROR) << "** Couldn't open: " << filename;
      return false;
    }

    fNEntries = fNuFluxTree->GetEntries();
  }

  LOG("Flux", pNOTICE)
    << "Loaded flux tree contains " <<  fNEntries << " entries";

  LOG("Flux", pDEBUG)
    << "Getting tree branches & setting leaf addresses";

  // try to get all the branches that we know about and only set address if
  // they exist
  bool missing_critical = false;
  TBranch * fBr = 0;
  GJPARCNuFluxPassThroughInfo * info = fPassThroughInfo;

  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("norm")) ) fBr->SetAddress(&info->norm);
  else if( (fBr = fNuFluxChain->GetBranch("norm")) ) fNuFluxChain->SetBranchAddress("norm",&info->norm);
  else {
    LOG("Flux", pFATAL) <<"Cannot find flux branch: norm";
    missing_critical = true;
  }
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("Enu")) ) fBr->SetAddress(&info->Enu);
  else if( (fBr = fNuFluxChain->GetBranch("Enu")) ) fNuFluxChain->SetBranchAddress("Enu",&info->Enu);
  else {
    LOG("Flux", pFATAL) <<"Cannot find flux branch: Enu";
    missing_critical = true;
  }
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ppid")) ) fBr->SetAddress(&info->ppid);
  else if( (fBr = fNuFluxChain->GetBranch("ppid")) ) fNuFluxChain->SetBranchAddress("ppid",&info->ppid);
  else {
    LOG("Flux", pFATAL) <<"Cannot find flux branch: ppid";
    missing_critical = true;
  }
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("mode")) ) fBr->SetAddress(&info->mode);
  else if( (fBr = fNuFluxChain->GetBranch("mode")) ) fNuFluxChain->SetBranchAddress("mode",&info->mode);
  else {
    LOG("Flux", pFATAL) <<"Cannot find flux branch: mode";
    missing_critical = true;
  }
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("rnu")) ) fBr->SetAddress(&info->rnu);
  else if( (fBr = fNuFluxChain->GetBranch("rnu")) ) fNuFluxChain->SetBranchAddress("rnu",&info->rnu);
  else if(fIsNDLoc){ // Only required for ND location
    LOG("Flux", pFATAL) <<"Cannot find flux branch: rnu";
    missing_critical = true;
  }
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("xnu")) ) fBr->SetAddress(&info->xnu);
  else if( (fBr = fNuFluxChain->GetBranch("xnu")) ) fNuFluxChain->SetBranchAddress("xnu",&info->xnu);
  else if(fIsNDLoc){ // Only required for ND location
    LOG("Flux", pFATAL) <<"Cannot find flux branch: xnu";
    missing_critical = true;
  }
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ynu")) ) fBr->SetAddress(&info->ynu);
  else if( (fBr = fNuFluxChain->GetBranch("ynu")) ) fNuFluxChain->SetBranchAddress("ynu",&info->ynu);
  else if(fIsNDLoc) { // Only required for ND location
    LOG("Flux", pFATAL) <<"Cannot find flux branch: ynu";
    missing_critical = true;
  }
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("nnu")) ) fBr->SetAddress(info->nnu);
  else if( (fBr = fNuFluxChain->GetBranch("nnu")) ) fNuFluxChain->SetBranchAddress("nnu",info->nnu);
  else if(fIsNDLoc){ // Only required for ND location
    LOG("Flux", pFATAL) <<"Cannot find flux branch: nnu";
    missing_critical = true;
  }
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("idfd")) ) fBr->SetAddress(&info->idfd);
  else if( (fBr = fNuFluxChain->GetBranch("idfd")) ) fNuFluxChain->SetBranchAddress("idfd",&info->idfd);
  else if(fIsNDLoc){ // Only required for ND location
    LOG("Flux", pFATAL) <<"Cannot find flux branch: idfd";
    missing_critical = true;
  }
  // check that have found essential branches
  if(missing_critical){
    LOG("Flux", pFATAL)
     << "Unable to find critical information in the flux ntuple! Initialization failed!";
    exit(1);
  }
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ppi"))     ) fBr->SetAddress(&info->ppi);
  else if( (fBr = fNuFluxChain->GetBranch("ppi"))     ) fNuFluxChain->SetBranchAddress("ppi",&info->ppi);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("xpi"))     ) fBr->SetAddress(info->xpi);
  else if( (fBr = fNuFluxChain->GetBranch("xpi"))     ) fNuFluxChain->SetBranchAddress("xpi",info->xpi);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("npi"))     ) fBr->SetAddress(info->npi);
  else if( (fBr = fNuFluxChain->GetBranch("npi"))     ) fNuFluxChain->SetBranchAddress("npi",info->npi);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ppi0"))    ) fBr->SetAddress(&info->ppi0);
  else if( (fBr = fNuFluxChain->GetBranch("ppi0"))    ) fNuFluxChain->SetBranchAddress("ppi0",&info->ppi0);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("xpi0"))    ) fBr->SetAddress(info->xpi0);
  else if( (fBr = fNuFluxChain->GetBranch("xpi0"))    ) fNuFluxChain->SetBranchAddress("xpi0",info->xpi0);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("npi0"))    ) fBr->SetAddress(info->npi0);
  else if( (fBr = fNuFluxChain->GetBranch("npi0"))    ) fNuFluxChain->SetBranchAddress("npi0",info->npi0);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("nvtx0"))   ) fBr->SetAddress(&info->nvtx0);
  else if( (fBr = fNuFluxChain->GetBranch("nvtx0"))   ) fNuFluxChain->SetBranchAddress("nvtx0",&info->nvtx0);
  // Following branches only present since flux version 10a
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("cospibm")) ) fBr->SetAddress(&info->cospibm);
  else if( (fBr = fNuFluxChain->GetBranch("cospibm")) ) fNuFluxChain->SetBranchAddress("cospibm",&info->cospibm);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("cospi0bm"))) fBr->SetAddress(&info->cospi0bm);
  else if( (fBr = fNuFluxChain->GetBranch("cospi0bm"))) fNuFluxChain->SetBranchAddress("cospi0bm",&info->cospi0bm);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gamom0"))  ) fBr->SetAddress(&info->gamom0);
  else if( (fBr = fNuFluxChain->GetBranch("gamom0"))  ) fNuFluxChain->SetBranchAddress("gamom0",&info->gamom0);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gipart"))  ) fBr->SetAddress(&info->gipart);
  else if( (fBr = fNuFluxChain->GetBranch("gipart"))  ) fNuFluxChain->SetBranchAddress("gipart",&info->gipart);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gvec0"))   ) fBr->SetAddress(info->gvec0);
  else if( (fBr = fNuFluxChain->GetBranch("gvec0"))   ) fNuFluxChain->SetBranchAddress("gvec0",info->gvec0);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpos0"))   ) fBr->SetAddress(info->gpos0);
  else if( (fBr = fNuFluxChain->GetBranch("gpos0"))   ) fNuFluxChain->SetBranchAddress("gpos0",info->gpos0);
  // Following branches only present since flux vesion 10d
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("ng"))      ) fBr->SetAddress(&info->ng);
  else if( (fBr = fNuFluxChain->GetBranch("ng"))      ) fNuFluxChain->SetBranchAddress("ng",&info->ng);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpid"))    ) fBr->SetAddress(info->gpid);
  else if( (fBr = fNuFluxChain->GetBranch("gpid"))    ) fNuFluxChain->SetBranchAddress("gpid",info->gpid);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gmec"))    ) fBr->SetAddress(info->gmec);
  else if( (fBr = fNuFluxChain->GetBranch("gmec"))    ) fNuFluxChain->SetBranchAddress("gmec",info->gmec);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gvx"))     ) fBr->SetAddress(info->gvx);
  else if( (fBr = fNuFluxChain->GetBranch("gvx"))     ) fNuFluxChain->SetBranchAddress("gvx",info->gvx);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gvy"))     ) fBr->SetAddress(info->gvy);
  else if( (fBr = fNuFluxChain->GetBranch("gvy"))     ) fNuFluxChain->SetBranchAddress("gvy",info->gvy);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gvz"))     ) fBr->SetAddress(info->gvz);
  else if( (fBr = fNuFluxChain->GetBranch("gvz"))     ) fNuFluxChain->SetBranchAddress("gvz",info->gvz);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpx"))     ) fBr->SetAddress(info->gpx);
  else if( (fBr = fNuFluxChain->GetBranch("gpx"))     ) fNuFluxChain->SetBranchAddress("gpx",info->gpx);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpy"))     ) fBr->SetAddress(info->gpy);
  else if( (fBr = fNuFluxChain->GetBranch("gpy"))     ) fNuFluxChain->SetBranchAddress("gpy",info->gpy);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gpz"))     ) fBr->SetAddress(info->gpz);
  else if( (fBr = fNuFluxChain->GetBranch("gpz"))     ) fNuFluxChain->SetBranchAddress("gpz",info->gpz);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gmat"))    ) fBr->SetAddress(info->gmat);
  else if( (fBr = fNuFluxChain->GetBranch("gmat"))    ) fNuFluxChain->SetBranchAddress("gmat",info->gmat);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gdistc"))  ) fBr->SetAddress(info->gdistc);
  else if( (fBr = fNuFluxChain->GetBranch("gdistc"))  ) fNuFluxChain->SetBranchAddress("gdistc",info->gdistc);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gdistal")) ) fBr->SetAddress(&info->gdistal);
  else if( (fBr = fNuFluxChain->GetBranch("gdistal")) ) fNuFluxChain->SetBranchAddress("gdistal",&info->gdistal);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gdistti")) ) fBr->SetAddress(&info->gdistti);
  else if( (fBr = fNuFluxChain->GetBranch("gdistti")) ) fNuFluxChain->SetBranchAddress("gdistti",&info->gdistti);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gdistfe")) ) fBr->SetAddress(&info->gdistfe);
  else if( (fBr = fNuFluxChain->GetBranch("gdistfe")) ) fNuFluxChain->SetBranchAddress("gdistfe",&info->gdistfe);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("gcosbm"))  ) fBr->SetAddress(info->gcosbm);
  else if( (fBr = fNuFluxChain->GetBranch("gcosbm"))  ) fNuFluxChain->SetBranchAddress("gcosbm",info->gcosbm);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("Enusk"))   ) fBr->SetAddress(&info->Enusk);
  else if( (fBr = fNuFluxChain->GetBranch("Enusk"))   ) fNuFluxChain->SetBranchAddress("Enusk",&info->Enusk);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("normsk"))  ) fBr->SetAddress(&info->normsk);
  else if( (fBr = fNuFluxChain->GetBranch("normsk"))  ) fNuFluxChain->SetBranchAddress("normsk",&info->normsk);
  if( fNuFluxUsingTree && (fBr = fNuFluxTree->GetBranch("anorm"))   ) fBr->SetAddress(&info->anorm);
  else if( (fBr = fNuFluxChain->GetBranch("anorm"))   ) fNuFluxChain->SetBranchAddress("anorm",&info->anorm);

  // Look for the flux file summary info tree (only expected for > 10a flux versions)
  if( fNuFluxUsingTree && (fNuFluxSumTree = (TTree*) fNuFluxFile->Get("h1000")) ){
    if( (fBr = fNuFluxSumTree->GetBranch("version"))) fBr->SetAddress(&info->version);
    if( (fBr = fNuFluxSumTree->GetBranch("ntrig"))  ) fBr->SetAddress(&info->ntrig);
    if( (fBr = fNuFluxSumTree->GetBranch("tuneid")) ) fBr->SetAddress(&info->tuneid);
    if( (fBr = fNuFluxSumTree->GetBranch("pint"))   ) fBr->SetAddress(&info->pint);
    if( (fBr = fNuFluxSumTree->GetBranch("bpos"))   ) fBr->SetAddress(info->bpos);
    if( (fBr = fNuFluxSumTree->GetBranch("btilt"))  ) fBr->SetAddress(info->btilt);
    if( (fBr = fNuFluxSumTree->GetBranch("brms"))   ) fBr->SetAddress(info->brms);
    if( (fBr = fNuFluxSumTree->GetBranch("emit"))   ) fBr->SetAddress(info->emit);
    if( (fBr = fNuFluxSumTree->GetBranch("alpha"))  ) fBr->SetAddress(info->alpha);
    if( (fBr = fNuFluxSumTree->GetBranch("hcur"))   ) fBr->SetAddress(info->hcur);
    if( (fBr = fNuFluxSumTree->GetBranch("rand"))   ) fBr->SetAddress(&info->rand);
    if( (fBr = fNuFluxSumTree->GetBranch("rseed"))  ) fBr->SetAddress(info->rseed);
  }

  // Look for the flux file summary info tree (only expected for > 10a flux versions)
  if ( !fNuFluxUsingTree ) {
    fNuFluxSumChain = new TChain("h1000");
    int result = fNuFluxSumChain->Add( Form("%s.%i.root",fileroot.c_str(),firstfile), 0 );
    if (result==1) {
      for (int i = firstfile+1; i < lastfile+1; i++) {
        result = fNuFluxSumChain->Add( Form("%s.%i.root",fileroot.c_str(),i), 0 );
      }
      if( (fBr = fNuFluxSumChain->GetBranch("version"))) fNuFluxSumChain->SetBranchAddress("version",&info->version);
      if( (fBr = fNuFluxSumChain->GetBranch("ntrig"))  ) fNuFluxSumChain->SetBranchAddress("ntrig",&info->ntrig);
      if( (fBr = fNuFluxSumChain->GetBranch("tuneid")) ) fNuFluxSumChain->SetBranchAddress("tuneid",&info->tuneid);
      if( (fBr = fNuFluxSumChain->GetBranch("pint"))   ) fNuFluxSumChain->SetBranchAddress("pint",&info->pint);
      if( (fBr = fNuFluxSumChain->GetBranch("bpos"))   ) fNuFluxSumChain->SetBranchAddress("bpos",info->bpos);
      if( (fBr = fNuFluxSumChain->GetBranch("btilt"))  ) fNuFluxSumChain->SetBranchAddress("btilt",info->btilt);
      if( (fBr = fNuFluxSumChain->GetBranch("brms"))   ) fNuFluxSumChain->SetBranchAddress("brms",info->brms);
      if( (fBr = fNuFluxSumChain->GetBranch("emit"))   ) fNuFluxSumChain->SetBranchAddress("emit",info->emit);
      if( (fBr = fNuFluxSumChain->GetBranch("alpha"))  ) fNuFluxSumChain->SetBranchAddress("alpha",info->alpha);
      if( (fBr = fNuFluxSumChain->GetBranch("hcur"))   ) fNuFluxSumChain->SetBranchAddress("hcur",info->hcur);
      if( (fBr = fNuFluxSumChain->GetBranch("rand"))   ) fNuFluxSumChain->SetBranchAddress("rand",&info->rand);
      if( (fBr = fNuFluxSumChain->GetBranch("rseed"))  ) fNuFluxSumChain->SetBranchAddress("rseed",info->rseed);
    }
  }

  // current ntuple cycle # (flux ntuples may be recycled)
  fICycle = 1;

  // sum-up weights & number of neutrinos for the specified location
  // over a complete cycle. Also record the maximum weight as previous
  // method using TTree::GetV1() seg faulted for more than ~1.5E6 entries
  fSumWeightTot1c  = 0;
  fNNeutrinosTot1c = 0;
  fNDetLocIdFound = 0;
  for(int ientry = 0; ientry < fNEntries; ientry++) {
     if (fNuFluxUsingTree)
       fNuFluxTree->GetEntry(ientry);
     else
       fNuFluxChain->GetEntry(ientry);
     // check for negative flux weights
     if(fPassThroughInfo->norm + controls::kASmallNum < 0.0){
       LOG("Flux", pERROR) << "Negative flux weight! Will set weight to 0.0";
       fPassThroughInfo->norm  = 0.0;
     }
     fNorm = (double) fPassThroughInfo->norm;
     // update maximum weight
     fMaxWeight = TMath::Max(fMaxWeight, (double) fPassThroughInfo->norm);
     // compare detector location (see GenerateNext_weighted() for details)
     if(fIsNDLoc && fDetLocId!=fPassThroughInfo->idfd) continue;
     fSumWeightTot1c += fNorm;
     fNNeutrinosTot1c++;
     fNDetLocIdFound++;
  }
  // Exit if have not found neutrino at specified location for whole cycle
  if(fNDetLocIdFound == 0){
    LOG("Flux", pFATAL)
     << "The input jnubeam flux ntuple contains no entries for detector id "
     << fDetLocId << ". Terminating job!";
    exit(1);
  }
  fNDetLocIdFound = 0; // reset the counter

  LOG("Flux", pNOTICE) << "Maximum flux weight = " << fMaxWeight;
  if(fMaxWeight <=0 ) {
      LOG("Flux", pFATAL) << "Non-positive maximum flux weight!";
      exit(1);
  }

  LOG("Flux", pINFO)
    << "Totals / cycle: #neutrinos = " << fNNeutrinosTot1c
    << ", Sum{Weights} = " << fSumWeightTot1c;

  if(fUseRandomOffset){
    this->RandomOffset();  // Random start point when looping over ntuple
  }

  return true;
}

double genie::flux::GJPARCNuFlux::MaxEnergy ( void  )

inlinevirtual

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

Implements genie::GFluxI.

Definition at line 60 of file GJPARCNuFlux.h.

{ return  fMaxEv;                }

const TLorentzVector& genie::flux::GJPARCNuFlux::Momentum ( void  )

inlinevirtual

returns the flux neutrino 4-momentum

Implements genie::GFluxI.

Definition at line 64 of file GJPARCNuFlux.h.

{ return  fgP4;                  }

long int genie::flux::GJPARCNuFlux::NFluxNeutrinos ( void  ) const

inline

number of flux neutrinos looped so far

Definition at line 88 of file GJPARCNuFlux.h.

const GJPARCNuFluxPassThroughInfo& genie::flux::GJPARCNuFlux::PassThroughInfo ( void  )

inline

GJPARCNuFluxPassThroughInfo.

Definition at line 92 of file GJPARCNuFlux.h.

int genie::flux::GJPARCNuFlux::PdgCode ( void  )

inlinevirtual

returns the flux neutrino pdg code

Implements genie::GFluxI.

Definition at line 62 of file GJPARCNuFlux.h.

{ return  fgPdgC;                }

const TLorentzVector& genie::flux::GJPARCNuFlux::Position ( void  )

inlinevirtual

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

Implements genie::GFluxI.

Definition at line 65 of file GJPARCNuFlux.h.

{ return  fgX4;                  }

double GJPARCNuFlux::POT_1cycle

(

void

)

flux POT per cycle

Definition at line 322 of file GJPARCNuFlux.cxx.

{
// Compute number of flux POTs / flux ntuple cycle
//
  if( (!fNuFluxTree && fNuFluxUsingTree) || (!fNuFluxChain && !fNuFluxUsingTree) ) {
     LOG("Flux", pWARN)
          << "The flux driver has not been properly configured";
     return 0;
  }

// double pot = fFilePOT * (fNNeutrinosTot1c/fSumWeightTot1c);
//
// Use the max weight instead, since flux neutrinos get de-weighted
// before thrown to the event generation driver
//
  double pot = fFilePOT / fMaxWeight;
  return pot;
}

double GJPARCNuFlux::POT_curravg

(

void

)

current average POT

Definition at line 341 of file GJPARCNuFlux.cxx.

{
// Compute current number of flux POTs
// On complete cycles, that POT number should be exact.
// Within cycles that is only an average number

  if( (!fNuFluxTree && fNuFluxUsingTree) || (!fNuFluxChain && !fNuFluxUsingTree) ) {
     LOG("Flux", pWARN)
          << "The flux driver has not been properly configured";
     return 0;
  }

// double pot = fNNeutrinos*fFilePOT/fSumWeightTot1c;
//
// See also comment at POT_1cycle()
//
  double cnt   = (double)fNNeutrinos;
  double cnt1c = (double)fNNeutrinosTot1c;
  double pot  = (cnt/cnt1c) * this->POT_1cycle();
  return pot;
}

void GJPARCNuFlux::RandomOffset

(

void

)

choose a random offset as starting entry in flux ntuple

Definition at line 790 of file GJPARCNuFlux.cxx.

{
// Choose a random number between 0-->fNEntries to set as start point for
// looping over flux ntuple. May be necessary when looping over very large
// flux files as always starting fromthe same point may introduce biases
// (inversely proportional to number of cycles). This method resets the
// starting value for fIEntry so must be called before any call to GenerateNext
// is made.
//
  double ran_frac = RandomGen::Instance()->RndFlux().Rndm();
  long int offset = (long int) floor(ran_frac * fNEntries);
  LOG("Flux", pERROR) << "Setting flux driver to start looping over entries "
                      << "with offset of "<< offset;
  fIEntry = fOffset = offset;
}

void GJPARCNuFlux::ResetCurrent ( void  )

private

Definition at line 895 of file GJPARCNuFlux.cxx.

{
// reset running values of neutrino pdg-code, 4-position & 4-momentum
// and the input ntuple leaves

  fLoadedNeutrino = false;

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

  fPassThroughInfo->Reset();
}

void GJPARCNuFlux::SetDefaults ( void  )

private

Definition at line 868 of file GJPARCNuFlux.cxx.

{
// - Set default neutrino species list (nue, nuebar, numu, numubar) and
//   maximum energy (25 GeV).
//   These defaults can be overwritten by user calls (at the driver init) to
//   GJPARCNuFlux::SetMaxEnergy(double Ev) and
//   GJPARCNuFlux::SetFluxParticles(const PDGCodeList & particles)
// - Set the default file normalization to 1E+21 POT
// - Set the default flux neutrino start z position at -5m (z=0 is the
//   detector centre).
// - Set number of cycles to 1

  LOG("Flux", pNOTICE) << "Setting default GJPARCNuFlux driver options";

  PDGCodeList particles;
  particles.push_back(kPdgNuMu);
  particles.push_back(kPdgAntiNuMu);
  particles.push_back(kPdgNuE);
  particles.push_back(kPdgAntiNuE);

  this->SetFluxParticles(particles);
  this->SetMaxEnergy(25./*GeV*/);
  this->SetFilePOT(1E+21);
  this->SetUpstreamZ(-5.0);
  this->SetNumOfCycles(1);
}

void GJPARCNuFlux::SetFilePOT

(

double

pot

)

flux file norm is in /N POT/det [ND] or /N POT/cm^2 [FD]. Specify N (typically 1E+21)

Definition at line 751 of file GJPARCNuFlux.cxx.

{
// The flux normalization is in /N POT/det [ND] or /N POT/cm^2 [FD]
// This method specifies N (typically 1E+21)

  fFilePOT = pot;
}

void GJPARCNuFlux::SetFluxParticles

(

const PDGCodeList &

particles

)

specify list of flux neutrino species

Definition at line 732 of file GJPARCNuFlux.cxx.

{
  if(!fPdgCList) {
     fPdgCList = new PDGCodeList;
  }
  fPdgCList->Copy(particles);

  LOG("Flux", pINFO)
    << "Declared list of neutrino species: " << *fPdgCList;
}

void GJPARCNuFlux::SetMaxEnergy

(

double

Ev

)

specify maximum flx neutrino energy

Definition at line 743 of file GJPARCNuFlux.cxx.

{
  fMaxEv = TMath::Max(0.,Ev);

  LOG("Flux", pINFO)
    << "Declared maximum flux neutrino energy: " << fMaxEv;
}

void GJPARCNuFlux::SetNumOfCycles

(

int

n

)

set how many times to cycle through the ntuple (default: 1 / n=0 means 'infinite')

Definition at line 769 of file GJPARCNuFlux.cxx.

{
// The flux ntuples can be recycled for a number of times to boost generated
// event statistics without requiring enormous beam simulation statistics.
// That option determines how many times the driver is going to cycle through
// the input flux ntuple.
// With n=0 the flux ntuple will be recycled an infinite amount of times so
// that the event generation loop can exit only on a POT or event num check.

  fNCycles = TMath::Max(0, n);
}

void GJPARCNuFlux::SetUpstreamZ

(

double

z0

)

set flux neutrino initial z position (upstream of the detector)

Definition at line 759 of file GJPARCNuFlux.cxx.

{
// The flux neutrino position (x,y) is given at the detector coord system
// at z=0. This method sets the preferred starting z position upstream of
// the upstream detector face. Each flux neutrino will be backtracked from
// z=0 to the input z0.

  fZ0 = z0;
}

double genie::flux::GJPARCNuFlux::SumWeight ( void  ) const

inline

intergated weight for flux neutrinos looped so far

Definition at line 89 of file GJPARCNuFlux.h.

double genie::flux::GJPARCNuFlux::Weight ( void  )

inlinevirtual

returns the flux neutrino weight (if any)

Implements genie::GFluxI.

Definition at line 63 of file GJPARCNuFlux.h.

{ return  fNorm / fMaxWeight;    }

Member Data Documentation

string genie::flux::GJPARCNuFlux::fDetLoc

private

input detector location ('sk','nd1','nd2',...)

Definition at line 120 of file GJPARCNuFlux.h.

int genie::flux::GJPARCNuFlux::fDetLocId

private

input detector location id (fDetLoc -> jnubeam idfd)

Definition at line 121 of file GJPARCNuFlux.h.

long int genie::flux::GJPARCNuFlux::fEntriesThisCycle

private

keep track of number of entries used so far for this cycle

Definition at line 127 of file GJPARCNuFlux.h.

double genie::flux::GJPARCNuFlux::fFilePOT

private

file POT normalization, typically 1E+21

Definition at line 131 of file GJPARCNuFlux.h.

bool genie::flux::GJPARCNuFlux::fGenerateWeighted

private

generate weighted/deweighted flux neutrinos (default is false)

Definition at line 139 of file GJPARCNuFlux.h.

TLorentzVector genie::flux::GJPARCNuFlux::fgP4

private

running generated nu 4-momentum

Definition at line 111 of file GJPARCNuFlux.h.

int genie::flux::GJPARCNuFlux::fgPdgC

private

running generated nu pdg-code

Definition at line 110 of file GJPARCNuFlux.h.

TLorentzVector genie::flux::GJPARCNuFlux::fgX4

private

running generated nu 4-position

Definition at line 112 of file GJPARCNuFlux.h.

int genie::flux::GJPARCNuFlux::fICycle

private

current cycle

Definition at line 134 of file GJPARCNuFlux.h.

long int genie::flux::GJPARCNuFlux::fIEntry

private

current flux ntuple entry

Definition at line 126 of file GJPARCNuFlux.h.

bool genie::flux::GJPARCNuFlux::fIsFDLoc

private

input location is a 'far' detector location?

Definition at line 123 of file GJPARCNuFlux.h.

bool genie::flux::GJPARCNuFlux::fIsNDLoc

private

input location is a 'near' detector location?

Definition at line 124 of file GJPARCNuFlux.h.

bool genie::flux::GJPARCNuFlux::fLoadedNeutrino

private

set to true when GenerateNext_weighted has been called successfully

Definition at line 141 of file GJPARCNuFlux.h.

double genie::flux::GJPARCNuFlux::fMaxEv

private

maximum energy

Definition at line 107 of file GJPARCNuFlux.h.

double genie::flux::GJPARCNuFlux::fMaxWeight

private

max flux neutrino weight in input file for the specified detector location

Definition at line 130 of file GJPARCNuFlux.h.

int genie::flux::GJPARCNuFlux::fNCycles

private

how many times to cycle through the flux ntuple

Definition at line 133 of file GJPARCNuFlux.h.

int genie::flux::GJPARCNuFlux::fNDetLocIdFound

private

per cycle keep track of the number of fDetLocId are found - if this is zero will exit job

Definition at line 122 of file GJPARCNuFlux.h.

long int genie::flux::GJPARCNuFlux::fNEntries

private

number of flux ntuple entries

Definition at line 125 of file GJPARCNuFlux.h.

long int genie::flux::GJPARCNuFlux::fNNeutrinos

private

number of flux neutrinos thrown so far

Definition at line 136 of file GJPARCNuFlux.h.

long int genie::flux::GJPARCNuFlux::fNNeutrinosTot1c

private

total number of flux neutrinos to be thrown / cycle

Definition at line 138 of file GJPARCNuFlux.h.

double genie::flux::GJPARCNuFlux::fNorm

private

current flux ntuple normalisation

Definition at line 129 of file GJPARCNuFlux.h.

TChain* genie::flux::GJPARCNuFlux::fNuFluxChain

private

input flux ntuple

Definition at line 116 of file GJPARCNuFlux.h.

TFile* genie::flux::GJPARCNuFlux::fNuFluxFile

private

input flux file

Definition at line 114 of file GJPARCNuFlux.h.

TChain* genie::flux::GJPARCNuFlux::fNuFluxSumChain

private

input summary ntuple for flux file. This tree is only present for later flux versions > 10a

Definition at line 118 of file GJPARCNuFlux.h.

TTree* genie::flux::GJPARCNuFlux::fNuFluxSumTree

private

input summary ntuple for flux file. This tree is only present for later flux versions > 10a

Definition at line 117 of file GJPARCNuFlux.h.

TTree* genie::flux::GJPARCNuFlux::fNuFluxTree

private

input flux ntuple

Definition at line 115 of file GJPARCNuFlux.h.

bool genie::flux::GJPARCNuFlux::fNuFluxUsingTree

private

are we using a TTree or a TChain to view the input flux file?

Definition at line 119 of file GJPARCNuFlux.h.

long int genie::flux::GJPARCNuFlux::fOffset

private

start looping at entry fOffset

Definition at line 128 of file GJPARCNuFlux.h.

GJPARCNuFluxPassThroughInfo* genie::flux::GJPARCNuFlux::fPassThroughInfo

private

Definition at line 143 of file GJPARCNuFlux.h.

PDGCodeList* genie::flux::GJPARCNuFlux::fPdgCList

private

list of neutrino pdg-codes

Definition at line 108 of file GJPARCNuFlux.h.

double genie::flux::GJPARCNuFlux::fSumWeight

private

sum of weights for neutrinos thrown so far

Definition at line 135 of file GJPARCNuFlux.h.

double genie::flux::GJPARCNuFlux::fSumWeightTot1c

private

total sum of weights for neutrinos to be thrown / cycle

Definition at line 137 of file GJPARCNuFlux.h.

bool genie::flux::GJPARCNuFlux::fUseRandomOffset

private

whether set random starting point when looping over flux ntuples

Definition at line 140 of file GJPARCNuFlux.h.

double genie::flux::GJPARCNuFlux::fZ0

private

configurable starting z position for each flux neutrino (in detector coord system)

Definition at line 132 of file GJPARCNuFlux.h.

---

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

• Generator/src/Tools/Flux/GJPARCNuFlux.h
• Generator/src/Tools/Flux/GJPARCNuFlux.cxx