GJPARCNuFlux.cxx

Go to the documentation of this file.

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

 Costas Andreopoulos <c.andreopoulos \at cern.ch>
 University of Liverpool 
*/
//____________________________________________________________________________

#include <cstdlib>
#include <iostream>
#include <cassert>

#include <TFile.h>
#include <TTree.h>
#include <TChain.h>
#include <TString.h>
#include <TSystem.h>

#include "Framework/Conventions/Units.h"
#include "Framework/Conventions/GBuild.h"
#include "Framework/Conventions/Controls.h"
#include "Tools/Flux/GJPARCNuFlux.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGCodeList.h"
#include "Framework/Numerical/MathUtils.h"
#include "Framework/Utils/PrintUtils.h"
#include "Framework/Utils/StringUtils.h"

#include "Tools/Flux/GFluxDriverFactory.h"
FLUXDRIVERREG4(genie,flux,GJPARCNuFlux,genie::flux::GJPARCNuFlux)

using std::endl;
using std::cout;
using std::endl;
using namespace genie;
using namespace genie::flux;

ClassImp(GJPARCNuFluxPassThroughInfo)

//____________________________________________________________________________
GJPARCNuFlux::GJPARCNuFlux()
{
  this->Initialize();
}
//___________________________________________________________________________
GJPARCNuFlux::~GJPARCNuFlux()
{
  this->CleanUp();
}
//___________________________________________________________________________
bool GJPARCNuFlux::GenerateNext(void)
{
// 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)
{
// 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;
}
//___________________________________________________________________________
double GJPARCNuFlux::POT_1cycle(void)
{
// 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)
{
// 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;
}
//___________________________________________________________________________
long int GJPARCNuFlux::Index(void)
{
// 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;
}
//___________________________________________________________________________
bool GJPARCNuFlux::LoadBeamSimData(string filename, string detector_location)
{
// 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;
}
//___________________________________________________________________________
void GJPARCNuFlux::SetFluxParticles(const PDGCodeList & particles)
{
  if(!fPdgCList) {
     fPdgCList = new PDGCodeList;
  }
  fPdgCList->Copy(particles);

  LOG("Flux", pINFO)
    << "Declared list of neutrino species: " << *fPdgCList;
}
//___________________________________________________________________________
void GJPARCNuFlux::SetMaxEnergy(double Ev)
{
  fMaxEv = TMath::Max(0.,Ev);

  LOG("Flux", pINFO)
    << "Declared maximum flux neutrino energy: " << fMaxEv;
}
//___________________________________________________________________________
void GJPARCNuFlux::SetFilePOT(double pot)
{
// 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::SetUpstreamZ(double z0)
{
// 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;
}
//___________________________________________________________________________
void GJPARCNuFlux::SetNumOfCycles(int n)
{
// 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::GenerateWeighted(bool gen_weighted)
{
// 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;
}
//___________________________________________________________________________
void GJPARCNuFlux::RandomOffset()
{
// 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::Clear(Option_t * opt)
{
// 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 GJPARCNuFlux::Initialize(void)
{
  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();
}
//___________________________________________________________________________
void GJPARCNuFlux::SetDefaults(void)
{
// - 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::ResetCurrent(void)
{
// 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::CleanUp(void)
{
  LOG("Flux", pNOTICE) << "Cleaning up...";

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

  if (fNuFluxFile) {
        fNuFluxFile->Close();
        delete fNuFluxFile;
  }
}
//___________________________________________________________________________
int GJPARCNuFlux::DLocName2Id(string name)
{
// 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;
}
//___________________________________________________________________________
GJPARCNuFluxPassThroughInfo::GJPARCNuFluxPassThroughInfo() :
TObject()
{
  this->Reset();
}
//___________________________________________________________________________
GJPARCNuFluxPassThroughInfo::GJPARCNuFluxPassThroughInfo(
                        const GJPARCNuFluxPassThroughInfo & info) :
TObject()
{
  fluxentry  = info.fluxentry;
  fluxfilename  = info.fluxfilename;
  Enu        = info.Enu;
  ppid       = info.ppid;
  mode       = info.mode;
  ppi        = info.ppi;
  ppi0       = info.ppi0;
  nvtx0      = info.nvtx0;
  cospibm    = info.cospibm;
  cospi0bm   = info.cospi0bm;
  idfd       = info.idfd;
  gamom0     = info.gamom0;
  gipart     = info.gipart;
  xnu        = info.xnu;
  ynu        = info.ynu;
  rnu        = info.rnu;
  for(int i = 0; i<3; i++){
    nnu[i] = info.nnu[i];
    xpi[i] = info.xpi[i];
    xpi0[i] = info.xpi0[i];
    gpos0[i] = info.gpos0[i];
    npi[i]  = info.npi[i];
    npi0[i] = info.npi0[i];
    gvec0[i] = info.gvec0[i];
  }
  ng = info.ng;
  for(int ip = 0; ip<fNgmax; ip++){
    gpid[ip] = info.gpid[ip];
    gmec[ip] = info.gmec[ip];
    gcosbm[ip] = info.gcosbm[ip];
    gvx[ip] = info.gvx[ip];
    gvy[ip] = info.gvy[ip];
    gvz[ip] = info.gvz[ip];
    gpx[ip] = info.gpx[ip];
    gpy[ip] = info.gpy[ip];
    gpz[ip] = info.gpz[ip];
    gmat[ip] = info.gmat[ip];
    gdistc[ip] = info.gdistc[ip];
    gdistal[ip] = info.gdistal[ip];
    gdistti[ip] = info.gdistti[ip];
    gdistfe[ip] = info.gdistfe[ip];
  }
  norm = info.norm;
  Enusk = info.Enusk;
  normsk = info.normsk;
  anorm = info.anorm;
  version = info.version;
  ntrig = info.ntrig;
  tuneid = info.tuneid;
  pint = info.pint;
  rand = info.rand;
  for(int i = 0; i < 2; i++){
    bpos[i] = info.bpos[i];
    btilt[i] = info.btilt[i];
    brms[i] = info.brms[i];
    emit[i] = info.emit[i];
    alpha[i] = info.alpha[i];
    rseed[i] = info.rseed[i];
  }
  for(int i = 0; i < 3; i++) hcur[i] = info.hcur[i];
}
//___________________________________________________________________________
void GJPARCNuFluxPassThroughInfo::Reset(){
  fluxentry= -1;
  fluxfilename = "Not-set";
  Enu      = 0;
  ppid     = 0;
  mode     = 0;
  ppi      = 0.;
  norm     = 0;
  cospibm  = 0.;
  nvtx0    = 0;
  ppi0     = 0.;
  idfd     = 0;
  rnu      = -999999.;
  xnu      = -999999.;
  ynu      = -999999.;
  cospi0bm = -999999.;
  gipart   = -1;
  gamom0   = 0.;
  for(int i=0; i<3; i++){
    nnu[i] = 0.;
    xpi[i] = -999999.;
    npi[i] = 0.;
    xpi0[i] = -999999.;
    npi0[i] = 0.;
    gpos0[i] = -999999.;
    gvec0[i] = 0.;
  }
  ng = -1;
  for(int ip = 0; ip<fNgmax; ip++){
    gpid[ip] = -999999;
    gmec[ip] = -999999;
    gcosbm[ip] = -999999.;
    gvx[ip] = -999999.;
    gvy[ip] = -999999.;
    gvz[ip] = -999999.;
    gpx[ip] = -999999.;
    gpy[ip] = -999999.;
    gpz[ip] = -999999.;
    gmat[ip] = -999999;
    gdistc[ip] = -999999.;
    gdistal[ip] = -999999.;
    gdistti[ip] = -999999.;
    gdistfe[ip] = -999999.;
  }
  Enusk = -999999.;
  normsk = -999999.;
  anorm = -999999.;
  version = -999999.;
  ntrig = -999999;
  tuneid = -999999;
  pint = -999999;
  rand = -999999;
  for(int i = 0; i < 2; i++){
    bpos[i] = -999999.;
    btilt[i] = -999999.;
    brms[i] = -999999.;
    emit[i] = -999999.;
    alpha[i] = -999999.;
    rseed[i] = -999999;
  }
  for(int i = 0; i < 3; i++) hcur[i] = -999999.;
}
//___________________________________________________________________________
namespace genie {
namespace flux  {
  ostream & operator << (
    ostream & stream, const genie::flux::GJPARCNuFluxPassThroughInfo & info)
    {
      stream << "\n idfd       = " << info.idfd
             << "\n norm       = " << info.norm
             << "\n flux entry = " << info.fluxentry
             << "\n flux file  = " << info.fluxfilename
             << "\n Enu        = " << info.Enu
             << "\n geant code = " << info.ppid
             << "\n (pdg code) = " << pdg::GeantToPdg(info.ppid)
             << "\n decay mode = " << info.mode
             << "\n nvtx0      = " << info.nvtx0
             << "\n |momentum| @ decay       = " << info.ppi
             << "\n position_vector @ decay  = ("
                 << info.xpi[0] << ", "
                 << info.xpi[1] << ", "
                 << info.xpi[2]  << ")"
             << "\n direction_vector @ decay = ("
                 << info.npi[0] << ", "
                 << info.npi[1] << ", "
                 << info.npi[2]  << ")"
             << "\n |momentum| @ prod.       = " << info.ppi0
             << "\n position_vector @ prod.  = ("
                 << info.xpi0[0] << ", "
                 << info.xpi0[1] << ", "
                 << info.xpi0[2]  << ")"
             << "\n direction_vector @ prod. = ("
                 << info.npi0[0] << ", "
                 << info.npi0[1] << ", "
                 << info.npi0[2]  << ")"
             << "\n cospibm = " << info.cospibm
             << "\n cospi0bm = " << info.cospi0bm
             << "\n Plus additional info if flux version is later than 07a"
             << endl;

    return stream;
  }
}//flux
}//genie
//___________________________________________________________________________