gT2KEvGen.cxx File Reference

#include <cassert>
#include <cstdlib>
#include <string>
#include <sstream>
#include <vector>
#include <map>
#include <TSystem.h>
#include <TTree.h>
#include <TFile.h>
#include <TH1D.h>
#include <TMath.h>
#include <TGeoVolume.h>
#include <TGeoShape.h>
#include <TList.h>
#include <TObject.h>
#include "Framework/Conventions/Units.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/EventGen/GFluxI.h"
#include "Framework/EventGen/GMCJDriver.h"
#include "Framework/EventGen/GMCJMonitor.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Ntuple/NtpWriter.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGCodeList.h"
#include "Framework/Ntuple/NtpMCFormat.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/Utils/XSecSplineList.h"
#include "Framework/Utils/StringUtils.h"
#include "Framework/Utils/UnitUtils.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"
#include "Framework/Utils/CmdLnArgParser.h"
#include "Framework/Utils/T2KEvGenMetaData.h"
#include "Framework/Utils/SystemUtils.h"
#include "Framework/Utils/PrintUtils.h"

Go to the source code of this file.

Functions
void	GetCommandLineArgs (int argc, char **argv)
void	PrintSyntax (void)
int	main (int argc, char **argv)

Variables
string	kDefOptGeomLUnits = "mm"
string	kDefOptGeomDUnits = "g_cm3"
NtpMCFormat_t	kDefOptNtpFormat = kNFGHEP
double	kDefOptFluxNorm = 1E+21
string	kDefOptEvFilePrefix = "gntp"
Long_t	gOptRunNu
bool	gOptUsingRootGeom = false
bool	gOptUsingHistFlux = false
map< int, double >	gOptTgtMix
map< int, TH1D * >	gOptFluxHst
string	gOptRootGeom
string	gOptRootGeomTopVol = ""
double	gOptGeomLUnits = 0
double	gOptGeomDUnits = 0
string	gOptExtMaxPlXml
string	gOptFluxFile
string	gOptDetectorLocation
double	gOptFluxNorm
PDGCodeList	gOptFluxNtpNuList (false)
int	gOptFluxNCycles
int	gOptNev
double	gOptPOT
bool	gOptExitAtEndOfFullFluxCycles
string	gOptEvFilePrefix
bool	gOptUseFluxProbs = false
bool	gOptSaveFluxProbsFile = false
string	gOptFluxProbFileName
string	gOptSaveFluxProbsFileName
bool	gOptRandomFluxOffset = false
long int	gOptRanSeed
string	gOptInpXSecFile

Function Documentation

void GetCommandLineArgs	(	int	argc,
		char **	argv
	)

Definition at line 936 of file gT2KEvGen.cxx.

{
  LOG("gevgen_t2k", pINFO) << "Parsing command line arguments";

  // Common run options. Set defaults and read.
  RunOpt::Instance()->EnableBareXSecPreCalc(true);
  RunOpt::Instance()->ReadFromCommandLine(argc,argv);

  // Parse run options for this app

  CmdLnArgParser parser(argc,argv);

  // help?
  bool help = parser.OptionExists('h');
  if(help) {
      PrintSyntax();
      exit(0);
  }

  // run number:
  if( parser.OptionExists('r') ) {
    LOG("gevgen_t2k", pDEBUG) << "Reading MC run number";
    gOptRunNu = parser.ArgAsLong('r');
  } else {
    LOG("gevgen_t2k", pDEBUG) << "Unspecified run number - Using default";
    gOptRunNu = 0;
  } //-r

  //
  // *** geometry
  //

  string geom = "";
  string lunits, dunits;
  if( parser.OptionExists('g') ) {
    LOG("gevgen_t2k", pDEBUG) << "Getting input geometry";
    geom = parser.ArgAsString('g');

    // is it a ROOT file that contains a ROOT geometry?
    bool accessible_geom_file =
            ! (gSystem->AccessPathName(geom.c_str()));
    if (accessible_geom_file) {
      gOptRootGeom      = geom;
      gOptUsingRootGeom = true;
    }
  } else {
      LOG("gevgen_t2k", pFATAL)
        << "No geometry option specified - Exiting";
      PrintSyntax();
      exit(1);
  } //-g

  if(gOptUsingRootGeom) {
     // using a ROOT geometry - get requested geometry units

     // legth units:
     if( parser.OptionExists('L') ) {
        LOG("gevgen_t2k", pDEBUG)
           << "Checking for input geometry length units";
        lunits = parser.ArgAsString('L');
     } else {
        LOG("gevgen_t2k", pDEBUG) << "Using default geometry length units";
        lunits = kDefOptGeomLUnits;
     } // -L
     // density units:
     if( parser.OptionExists('D') ) {
        LOG("gevgen_t2k", pDEBUG)
           << "Checking for input geometry density units";
        dunits = parser.ArgAsString('D');
     } else {
        LOG("gevgen_t2k", pDEBUG) << "Using default geometry density units";
        dunits = kDefOptGeomDUnits;
     } // -D
     gOptGeomLUnits = genie::utils::units::UnitFromString(lunits);
     gOptGeomDUnits = genie::utils::units::UnitFromString(dunits);

     // check whether an event generation volume name has been
     // specified -- default is the 'top volume'
     if( parser.OptionExists('t') ) {
        LOG("gevgen_t2k", pDEBUG) << "Checking for input volume name";
        gOptRootGeomTopVol = parser.ArgAsString('t');
     } else {
        LOG("gevgen_t2k", pDEBUG) << "Using the <master volume>";
     } // -t

     // check whether an XML file with the maximum (density weighted)
     // path lengths for each detector material is specified -
     // otherwise will compute the max path lengths at job init
     if( parser.OptionExists('m') ) {
        LOG("gevgen_t2k", pDEBUG)
              << "Checking for maximum path lengths XML file";
        gOptExtMaxPlXml = parser.ArgAsString('m');
     } else {
        LOG("gevgen_t2k", pDEBUG)
           << "Will compute the maximum path lengths at job init";
        gOptExtMaxPlXml = "";
     } // -m
  } // using root geom?

  else {
    // User has specified a target mix.
    // Decode the list of target pdf codes & their corresponding weight fraction
    // (specified as 'pdg_code_1[fraction_1],pdg_code_2[fraction_2],...')
    // See documentation on top section of this file.
    //
    gOptTgtMix.clear();
    vector<string> tgtmix = utils::str::Split(geom,",");
    if(tgtmix.size()==1) {
         int    pdg = atoi(tgtmix[0].c_str());
         double wgt = 1.0;
         gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
    } else {
      vector<string>::const_iterator tgtmix_iter = tgtmix.begin();
      for( ; tgtmix_iter != tgtmix.end(); ++tgtmix_iter) {
         string tgt_with_wgt = *tgtmix_iter;
         string::size_type open_bracket  = tgt_with_wgt.find("[");
         string::size_type close_bracket = tgt_with_wgt.find("]");
         if (open_bracket ==string::npos ||
             close_bracket==string::npos)
         {
             LOG("gevgen_t2k", pFATAL)
                << "You made an error in specifying the target mix";
             PrintSyntax();
             exit(1);
         }
         string::size_type ibeg = 0;
         string::size_type iend = open_bracket;
         string::size_type jbeg = open_bracket+1;
         string::size_type jend = close_bracket;
         int    pdg = atoi(tgt_with_wgt.substr(ibeg,iend-ibeg).c_str());
         double wgt = atof(tgt_with_wgt.substr(jbeg,jend-jbeg).c_str());
         LOG("gevgen_t2k", pDEBUG)
            << "Adding to target mix: pdg = " << pdg << ", wgt = " << wgt;
         gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));

      }// tgtmix_iter
    } // >1 materials in mix
  } // using tgt mix?

  //
  // *** flux
  //

  if( parser.OptionExists('f') ) {
    LOG("gevgen_t2k", pDEBUG) << "Getting input flux";
    string flux = parser.ArgAsString('f');
    gOptUsingHistFlux = (flux.find("[") != string::npos);

    if(!gOptUsingHistFlux) {
        // Using JNUBEAM flux ntuples
        // Extract JNUBEAM flux (root) file name & detector location.
        // Also extract the list of JNUBEAM neutrinos to consider (if none
        // is specified here then I will consider all flavours)
        //
        vector<string> fluxv = utils::str::Split(flux,",");
        if(fluxv.size()<2) {
           LOG("gevgen_t2k", pFATAL)
             << "You need to specify both a flux ntuple ROOT file "
             << " _AND_ a detector location";
           PrintSyntax();
           exit(1);
        }
        gOptFluxFile         = fluxv[0];
        gOptDetectorLocation = fluxv[1];

        // Extract the list of neutrinos to consider (if any).
        //
        for(unsigned int inu = 2; inu < fluxv.size(); inu++)
        {
           gOptFluxNtpNuList.push_back( atoi(fluxv[inu].c_str()) );
        }

    } else {
        // Using flux from histograms
        // Extract the root file name & the list of histogram names & neutrino
        // species (specified as 'filename,histo1[species1],histo2[species2],...')
        // See documentation on top section of this file.
        //
        vector<string> fluxv = utils::str::Split(flux,",");
        if(fluxv.size()<2) {
           LOG("gevgen_t2k", pFATAL)
             << "You need to specify both a flux ntuple ROOT file "
             << " _AND_ a detector location";
           PrintSyntax();
           exit(1);
        }
        gOptFluxFile = fluxv[0];
        bool accessible_flux_file =
               !(gSystem->AccessPathName(gOptFluxFile.c_str()));
        if (!accessible_flux_file) {
            LOG("gevgen_t2k", pFATAL)
              << "Can not access flux file: " << gOptFluxFile;
            PrintSyntax();
            exit(1);
        }
        // Extract energy spectra for all specified neutrino species
        TFile flux_file(gOptFluxFile.c_str(), "read");
        for(unsigned int inu=1; inu<fluxv.size(); inu++) {
            string nutype_and_histo = fluxv[inu];
            string::size_type open_bracket  = nutype_and_histo.find("[");
            string::size_type close_bracket = nutype_and_histo.find("]");
            if (open_bracket ==string::npos ||
                close_bracket==string::npos)
            {
                LOG("gevgen_t2k", pFATAL)
                   << "You made an error in specifying the flux histograms";
                PrintSyntax();
                exit(1);
            }
            string::size_type ibeg = 0;
            string::size_type iend = open_bracket;
            string::size_type jbeg = open_bracket+1;
            string::size_type jend = close_bracket;
            string nutype = nutype_and_histo.substr(ibeg,iend-ibeg);
            string histo  = nutype_and_histo.substr(jbeg,jend-jbeg);
            // access specified histogram from the input root file
            TH1D * ihst = (TH1D*) flux_file.Get(histo.c_str());
            if(!ihst) {
                LOG("gevgen_t2k", pFATAL)
                  << "Can not find histogram: " << histo
                  << " in flux file: " << gOptFluxFile;
                PrintSyntax();
                exit(1);
            }
            // create a local copy of the input histogram
            TString origname = ihst->GetName();
            TString tmpname; tmpname.Form("%s_", origname.Data());
            // Copy in the flux histogram from the root file
            // use Clone rather than assuming fix bin widths and rebooking
            TH1D* spectrum = (TH1D*)ihst->Clone();
            spectrum->SetNameTitle(tmpname.Data(),ihst->GetName());
            spectrum->SetDirectory(0);

            // get rid of original
            delete ihst;
            // rename copy
            spectrum->SetName(origname.Data());

            // convert neutrino name -> pdg code
            int pdg = atoi(nutype.c_str());
            if(!pdg::IsNeutrino(pdg) && !pdg::IsAntiNeutrino(pdg)) {
                LOG("gevgen_t2k", pFATAL)
                    << "Unknown neutrino type: " << nutype;
                PrintSyntax();
                exit(1);
            }
            // store flux neutrino code / energy spectrum
            LOG("gevgen_t2k", pDEBUG)
              << "Adding energy spectrum for flux neutrino: pdg = " << pdg;
            gOptFluxHst.insert(map<int, TH1D*>::value_type(pdg, spectrum));
        }//inu
        if(gOptFluxHst.size()<1) {
           LOG("gevgen_t2k", pFATAL)
               << "You have not specified any flux histogram!";
           PrintSyntax();
           exit(1);
        }
        flux_file.Close();
    } // flux from histograms or from JNUBEAM ntuples?

  } else {
      LOG("gevgen_t2k", pFATAL) << "No flux info was specified - Exiting";
      PrintSyntax();
      exit(1);
  }

  // Use a random offset when looping over flux entries
  if(parser.OptionExists('R')) gOptRandomFluxOffset = true;

  //
  // *** pre-calculated flux interaction probabilities
  //

  // using pre-calculated flux interaction probabilities
  if( parser.OptionExists('P') ){
    gOptFluxProbFileName = parser.ArgAsString('P');
    if(gOptFluxProbFileName.length() > 0){
      gOptUseFluxProbs = true;
      bool accessible =
              !(gSystem->AccessPathName(gOptFluxProbFileName.c_str()));
      if(!accessible){
        LOG("gevgen_t2k", pFATAL)
          << "Can not access pre-calculated flux probabilities file: " << gOptFluxProbFileName;
        PrintSyntax();
        exit(1);
      }
    }
    else {
      LOG("gevgen_t2k", pFATAL)
        << "No flux interaction probabilites were specified - exiting";
      PrintSyntax();
      exit(1);
    }
  }

  // pre-generating interaction probs and saving to output file
  if( parser.OptionExists('S') ){
    gOptSaveFluxProbsFile = true;
    gOptSaveFluxProbsFileName = parser.ArgAsString('S');
  }

  // cannot save and run at the same time
  if(gOptUseFluxProbs && gOptSaveFluxProbsFile){
    LOG("gevgen_t2k", pFATAL)
     << "Cannot specify both the -P and -S options at the same time!";
    exit(1);
  }

  // only makes sense to be setting these options for a realistic flux
  if(gOptUsingHistFlux && (gOptUseFluxProbs || gOptSaveFluxProbsFile)){
    LOG("gevgen_t2k", pFATAL)
     << "Using pre-calculated flux interaction probabilities only makes "
     << "sense when using JNUBEAM flux option!";
    exit(1);
  }

  // flux file POT normalization
  // only relevant when using the JNUBEAM flux ntuples
  if( parser.OptionExists('p') ) {
    LOG("gevgen_t2k", pDEBUG)  << "Reading flux file normalization";
    gOptFluxNorm = parser.ArgAsDouble('p');
  } else {
    LOG("gevgen_t2k", pDEBUG)
        << "Setting standard normalization for JNUBEAM flux ntuples";
    gOptFluxNorm = kDefOptFluxNorm;
  } //-p

  // number of times to cycle through the JNUBEAM flux ntuple contents
  if( parser.OptionExists('c') ) {
    LOG("gevgen_t2k", pDEBUG) << "Reading number of flux ntuple cycles";
    gOptFluxNCycles = parser.ArgAsInt('c');
  } else {
    LOG("gevgen_t2k", pDEBUG)
        << "Setting standard number of cycles for JNUBEAM flux ntuples";
    gOptFluxNCycles = -1;
  } //-c

  // limit on max number of events that can be generated
  if( parser.OptionExists('n') ) {
    LOG("gevgen_t2k", pDEBUG)
        << "Reading limit on number of events to generate";
    gOptNev = parser.ArgAsInt('n');
  } else {
    LOG("gevgen_t2k", pDEBUG)
      << "Will keep on generating events till the flux driver stops";
    gOptNev = -1;
  } //-n

  // exposure (in POT)
  bool uppc_e = parser.OptionExists('E');
  char pot_args = 'e';
  bool pot_exit = true;
  if(uppc_e) {
    pot_args = 'E';
    pot_exit = false;
  }
  gOptExitAtEndOfFullFluxCycles = pot_exit;
  if( parser.OptionExists(pot_args) ) {
    LOG("gevgen_t2k", pDEBUG)  << "Reading requested exposure in POT";
    gOptPOT = parser.ArgAsDouble(pot_args);
  } else {
    LOG("gevgen_t2k", pDEBUG) << "No POT exposure was requested";
    gOptPOT = -1;
  } //-e, -E

  // event file prefix
  if( parser.OptionExists('o') ) {
    LOG("gevgen_t2k", pDEBUG) << "Reading the event filename prefix";
    gOptEvFilePrefix = parser.ArgAsString('o');
  } else {
    LOG("gevgen_t2k", pDEBUG)
      << "Will set the default event filename prefix";
    gOptEvFilePrefix = kDefOptEvFilePrefix;
  } //-o


  // random number seed
  if( parser.OptionExists("seed") ) {
    LOG("gevgen_t2k", pINFO) << "Reading random number seed";
    gOptRanSeed = parser.ArgAsLong("seed");
  } else {
    LOG("gevgen_t2k", pINFO) << "Unspecified random number seed - Using default";
    gOptRanSeed = -1;
  }

  // input cross-section file
  if( parser.OptionExists("cross-sections") ) {
    LOG("gevgen_t2k", pINFO) << "Reading cross-section file";
    gOptInpXSecFile = parser.ArgAsString("cross-sections");
  } else {
    LOG("gevgen_t2k", pINFO) << "Unspecified cross-section file";
    gOptInpXSecFile = "";
  }

  //
  // >>> perform 'sanity' checks on command line arguments
  //

  // If we use a JNUBEAM flux ntuple, the 'exposure' may be set
  // either as:
  // - a number of POTs (whichever number of flux ntuple cycles that corresponds to)
  // - a number of generated events (whichever number of POTs that corresponds to)
  // - a number of flux ntuple cycles (whichever number of POTs that corresponds to)
  // Only one of those options can be set.
  if(!gOptUsingHistFlux) {
    int nset=0;
    if(gOptPOT         > 0) nset++;
    if(gOptFluxNCycles > 0) nset++;
    if(gOptNev         > 0) nset++;
    if(nset==0) {
       LOG("gevgen_t2k", pWARN)
        << "** To use a JNUBEAM flux ntuple you need to specify an exposure, "
        << "either via the -c, -e or -n options";
       LOG("gevgen_t2k", pWARN)
        << "** gevgen_t2k automatically sets the exposure via '-c 1'";
       gOptFluxNCycles = 1;
    }
    if(nset>1) {
       LOG("gevgen_t2k", pFATAL)
         << "You can not specify more than one of the -c, -e or -n options";
       PrintSyntax();
       exit(1);
    }
  }
  // If we use a flux histograms (not JNUBEAM flux ntuples) then -currently- the
  // only way to control exposure is via a number of events
  if(gOptUsingHistFlux) {
     if(gOptNev < 0) {
       LOG("gevgen_t2k", pFATAL)
         << "If you're using flux from histograms you need to specify the -n option";
       PrintSyntax();
       exit(1);
     }
  }
  // If we don't use a detailed ROOT detector geometry (but just a target mix) then
  // don't accept POT as a way to control job statistics (not enough info is passed
  // in the target mix to compute POT & the calculation can be easily done offline)
  if(!gOptUsingRootGeom) {
    if(gOptPOT > 0) {
       LOG("gevgen_t2k", pFATAL)
         << "You may not use the -e, -E options "
         << "without a detailed detector geometry description input";
       exit(1);
    }
  }

  //
  // >>> print the command line options
  //
  PDGLibrary * pdglib = PDGLibrary::Instance();

  ostringstream gminfo;
  if (gOptUsingRootGeom) {
    gminfo << "Using ROOT geometry - file: " << gOptRootGeom
           << ", top volume: "
           << ((gOptRootGeomTopVol.size()==0) ? "<master volume>" : gOptRootGeomTopVol)
           << ", max{PL} file: "
           << ((gOptExtMaxPlXml.size()==0) ? "<none>" : gOptExtMaxPlXml)
           << ", length  units: " << lunits
           << ", density units: " << dunits;
  } else {
    gminfo << "Using target mix - ";
    map<int,double>::const_iterator iter;
    for(iter = gOptTgtMix.begin(); iter != gOptTgtMix.end(); ++iter) {
          int    pdg_code = iter->first;
          double wgt      = iter->second;
          TParticlePDG * p = pdglib->Find(pdg_code);
          if(p) {
            string name = p->GetName();
            gminfo << "(" << name << ") -> " << 100*wgt << "% / ";
          }//p?
    }
  }

  ostringstream fluxinfo;
  if(gOptUsingHistFlux) {
    fluxinfo << "Using flux histograms - ";
    map<int,TH1D*>::const_iterator iter;
    for(iter = gOptFluxHst.begin(); iter != gOptFluxHst.end(); ++iter) {
          int    pdg_code = iter->first;
          TH1D * spectrum = iter->second;
          TParticlePDG * p = pdglib->Find(pdg_code);
          if(p) {
            string name = p->GetName();
            fluxinfo << "(" << name << ") -> " << spectrum->GetName() << " / ";
          }//p?
    }
  } else {
    fluxinfo << "Using JNUBEAM flux ntuple - "
             << "file: "        << gOptFluxFile
             << ", location: "  << gOptDetectorLocation
             << ", pot norm: "  << gOptFluxNorm;

     if( gOptFluxNtpNuList.size() > 0 ) {
       fluxinfo << ", ** this job is considering only: ";
       for(unsigned int inu = 0; inu < gOptFluxNtpNuList.size(); inu++) {
         fluxinfo << gOptFluxNtpNuList[inu];
         if(inu < gOptFluxNtpNuList.size()-1) fluxinfo << ",";
       }
     }

  }

  ostringstream exposure;
  if(gOptPOT > 0)
      exposure << "Number of POTs = " << gOptPOT;
  if(gOptFluxNCycles > 0)
      exposure << "Number of flux cycles = " << gOptFluxNCycles;
  if(gOptNev > 0)
      exposure << "Number of events = " << gOptNev;

  LOG("gevgen_t2k", pNOTICE)
     << "\n\n"
     << utils::print::PrintFramedMesg("T2K event generation job configuration");

  LOG("gevgen_t2k", pNOTICE)
     << "\n -  Run number: " << gOptRunNu
     << "\n -  Random number seed: " << gOptRanSeed
     << "\n - Using cross-section file: " << gOptInpXSecFile
     << "\n - Flux     @ " << fluxinfo.str()
     << "\n - Geometry @ " << gminfo.str()
     << "\n - Exposure @ " << exposure.str();

  LOG("gevgen_t2k", pNOTICE) << *RunOpt::Instance();
}

int main	(	int	argc,
		char **	argv
	)

Definition at line 512 of file gT2KEvGen.cxx.

{
  // Parse command line arguments
  GetCommandLineArgs(argc,argv);


  if ( ! RunOpt::Instance()->Tune() ) {
    LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
    exit(-1);
  }
  RunOpt::Instance()->BuildTune();

  // Initialization of random number generators, cross-section table,
  // messenger thresholds, cache file
  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
  utils::app_init::CacheFile(RunOpt::Instance()->CacheFile());
  utils::app_init::RandGen(gOptRanSeed);
  utils::app_init::XSecTable(gOptInpXSecFile, true);

  // Set GHEP print level
  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());

  // *************************************************************************
  // * Create / configure the geometry driver
  // *************************************************************************
  GeomAnalyzerI * geom_driver = 0;
  double zmin=0, zmax=0;

  if(gOptUsingRootGeom) {
    //
    // *** Using a realistic root-based detector geometry description
    //

    // creating & configuring a root geometry driver
    geometry::ROOTGeomAnalyzer * rgeom =
            new geometry::ROOTGeomAnalyzer(gOptRootGeom);
    rgeom -> SetLengthUnits  (gOptGeomLUnits);
    rgeom -> SetDensityUnits (gOptGeomDUnits);
    // getting the bounding box dimensions, before setting topvolume,
    // along z so as to set the appropriate upstream generation surface
    // for the JPARC flux driver
    TGeoVolume * topvol = rgeom->GetGeometry()->GetTopVolume();
    TGeoShape * bounding_box = topvol->GetShape();
    bounding_box->GetAxisRange(3, zmin, zmax);
    zmin *= rgeom->LengthUnits();
    zmax *= rgeom->LengthUnits();
    // now update to the requested topvolume for use in recursive exhaust method
    rgeom -> SetTopVolName   (gOptRootGeomTopVol);
    topvol = rgeom->GetGeometry()->GetTopVolume();
    if(!topvol) {
      LOG("gevgen_t2k", pFATAL) << "Null top ROOT geometry volume!";
      exit(1);
    }
    // switch on/off volumes as requested
    if ( (gOptRootGeomTopVol[0] == '+') || (gOptRootGeomTopVol[0] == '-') ) {
      bool exhaust = (*gOptRootGeomTopVol.c_str() == '+');
      utils::geometry::RecursiveExhaust(topvol, gOptRootGeomTopVol, exhaust);
    }

    // casting to the GENIE geometry driver interface
    geom_driver = dynamic_cast<GeomAnalyzerI *> (rgeom);
  }
  else {
    //
    // *** Using a 'point' geometry with the specified target mix
    // *** ( = a list of targets with their corresponding weight fraction)
    //

    // creating & configuring a point geometry driver
    geometry::PointGeomAnalyzer * pgeom =
              new geometry::PointGeomAnalyzer(gOptTgtMix);
    // casting to the GENIE geometry driver interface
    geom_driver = dynamic_cast<GeomAnalyzerI *> (pgeom);
  }

  // *************************************************************************
  // * Create / configure the flux driver
  // *************************************************************************
  GFluxI * flux_driver = 0;

  flux::GJPARCNuFlux *   jparc_flux_driver = 0;
  flux::GCylindTH1Flux * hst_flux_driver   = 0;

  if(!gOptUsingHistFlux) {
    //
    // *** Using the detailed JPARC neutrino flux desription by feeding-in
    // *** the JNUBEAM flux simulation ntuples
    //

    // creating JPARC neutrino flux driver
    jparc_flux_driver = new flux::GJPARCNuFlux;
    // before loading the beam sim data set whether to use a random offset when looping
    if(gOptRandomFluxOffset == false) jparc_flux_driver->DisableOffset();
    // specify input JNUBEAM file & detector location
    bool beam_sim_data_success = jparc_flux_driver->LoadBeamSimData(gOptFluxFile, gOptDetectorLocation);
    if(!beam_sim_data_success) {
      LOG("gevgen_t2k", pFATAL)
        << "The flux driver has not been properly configured. Exiting";
      exit(1);
    }
    // specify JNUBEAM normalization
    jparc_flux_driver->SetFilePOT(gOptFluxNorm);
    // specify upstream generation surface
    jparc_flux_driver->SetUpstreamZ(zmin);
    // specify which flavours to consider -
    // if no neutrino list was specified then the MC job will consider all flavours
    if( gOptFluxNtpNuList.size() > 0 ) {
       jparc_flux_driver->SetFluxParticles(gOptFluxNtpNuList);
    }

    // casting to the GENIE flux driver interface
    flux_driver = dynamic_cast<GFluxI *> (jparc_flux_driver);
  }
  else {
    //
    // *** Using fluxes from histograms (for all specified neutrino species)
    //

    // creating & configuring a generic GCylindTH1Flux flux driver
    TVector3 bdir (0,0,1); // dir along +z
    TVector3 bspot(0,0,0);
    hst_flux_driver = new flux::GCylindTH1Flux;
    hst_flux_driver->SetNuDirection      (bdir);
    hst_flux_driver->SetBeamSpot         (bspot);
    hst_flux_driver->SetTransverseRadius (-1);
    map<int,TH1D*>::iterator it = gOptFluxHst.begin();
    for( ; it != gOptFluxHst.end(); ++it) {
        int    pdg_code = it->first;
        TH1D * spectrum = new TH1D(*(it->second));
        hst_flux_driver->AddEnergySpectrum(pdg_code, spectrum);
    }
    // casting to the GENIE flux driver interface
    flux_driver = dynamic_cast<GFluxI *> (hst_flux_driver);
  }

  // *************************************************************************
  // * Create/configure the event generation driver
  // *************************************************************************
  GMCJDriver * mcj_driver = new GMCJDriver;
  mcj_driver->SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
  mcj_driver->UseFluxDriver(flux_driver);
  mcj_driver->UseGeomAnalyzer(geom_driver);
  mcj_driver->UseMaxPathLengths(gOptExtMaxPlXml);
  // do not calculate probability scales if using pre-generated flux probs
  bool calc_prob_scales = (gOptSaveFluxProbsFile || gOptUseFluxProbs) ? false : true;
  mcj_driver->Configure(calc_prob_scales);
  mcj_driver->UseSplines();
  mcj_driver->ForceSingleProbScale();

  // *************************************************************************
  // * If specified use pre-calculated flux interaction probabilities instead
  // * of preselecting based on max path lengths
  // *************************************************************************

  if(gOptUseFluxProbs || gOptSaveFluxProbsFile){

    bool success = false;

    // set flux probs output file name
    if(gOptSaveFluxProbsFile){
      // default output name is ${FLUFILENAME}.${TOPVOL}.flxprobs.root
      string basename = gOptFluxFile.substr(gOptFluxFile.rfind("/")+1);
      string name = basename.substr(0, basename.rfind("."));
      if(gOptRootGeomTopVol.length()>0)
        name += "."+gOptRootGeomTopVol+".flxprobs.root";
      else
        name += ".master.flxprobs.root";
      // if specified override with cmd line option
      if(gOptSaveFluxProbsFileName.size()>0) name = gOptSaveFluxProbsFileName;
      // Tell the driver save pre-generated probabilities to an output file
      mcj_driver->SaveFluxProbabilities(name);
    }

    // Either load pre-generated flux probabilities
    if(gOptFluxProbFileName.size() > 0){
      success = mcj_driver->LoadFluxProbabilities(gOptFluxProbFileName);
    }
    // Or pre-calculate them
    else success = mcj_driver->PreCalcFluxProbabilities();

    if(success){
      LOG("gevgen_t2k", pNOTICE)
       << "Successfully calculated/loaded flux interaction probabilities!";
      // Print out a list of expected number of events per POT and per cycle
      // based on the pre-generated flux interaction probabilities
      map<int, double> sum_probs_map = mcj_driver->SumFluxIntProbs();
      map<int, double>::const_iterator sum_probs_it = sum_probs_map.begin();
      double ntot_per_pot = 0.0;
      double ntot_per_cycle = 0.0;
      double pscale = mcj_driver->GlobProbScale();
      double pot_1cycle = jparc_flux_driver->POT_1cycle();
      LOG("T2KProdInfo", pNOTICE) <<
          "Expected event rates based on flux interaction probabilities:";
      for(; sum_probs_it != sum_probs_map.end(); sum_probs_it++){
        double sum_probs = sum_probs_it->second;
        double nevts_per_cycle = sum_probs / pscale;  // take into account rescale
        double nevts_per_pot = sum_probs/pot_1cycle;  // == (sum_probs*pscale)/(pot_1cycle*pscale)
        ntot_per_pot += nevts_per_pot;
        ntot_per_cycle += nevts_per_cycle;
        LOG("T2KProdInfo", pNOTICE) <<
            "     PDG "<< sum_probs_it->first << ": " << nevts_per_cycle <<
            " Events/Cycle, "<< nevts_per_pot << " Events/POT";
      }
      LOG("T2KProdInfo", pNOTICE) << "      -----------";
      LOG("T2KProdInfo", pNOTICE) <<
          "     All neutrino species: " << ntot_per_cycle <<
          " Events/Cycle, "<< ntot_per_pot << " Events/POT";
      LOG("T2KProdInfo", pNOTICE) <<
          "N.B. This assumes input flux file corresponds to "<< pot_1cycle <<
          "POT, ensure this is correct if using these numbers!";
    }
    else {
      LOG("gevgen_t2k", pFATAL)
       << "Failed to calculated/loaded flux interaction probabilities!";
      return 1;
    }

    // Exit now if just pre-generating interaction probabilities
    if(gOptSaveFluxProbsFile){
      LOG("gevgen_t2k", pNOTICE)
       << "Will not generate events - just pre-calculating flux interaction"
       << "probabilities";
      return 0;
    }
  } // Pre-calculated flux interaction probabilities

  // *************************************************************************
  // * Work out number of cycles for current exposure settings
  // *************************************************************************

  if(!gOptUsingHistFlux) {
    // If a number of POT was requested, then work out how many flux ntuple
    // cycles are required for accumulating those statistics
    if(gOptPOT>0) {
      double fpot_1c = jparc_flux_driver->POT_1cycle(); // flux POT / cycle
      double psc     = mcj_driver->GlobProbScale();     // interaction prob. scale
      double pot_1c  = fpot_1c / psc;                   // actual POT / cycle
      int    ncycles = (int) TMath::Max(1., TMath::Ceil(gOptPOT/pot_1c));

      LOG("gevgen_t2k", pNOTICE)
         << " *** POT/cycle:  " << pot_1c;
      LOG("gevgen_t2k", pNOTICE)
         << " *** Requested POT will be accumulated in: "
         << ncycles << " flux ntuple cycles";

      jparc_flux_driver->SetNumOfCycles(ncycles);
    }
    // If a number of events was requested, then set the number of flux
    // ntuple cycles to 'infinite'
    else if(gOptNev>0) {
       jparc_flux_driver->SetNumOfCycles(0);
    }
    // Just set the number of cycles to the requested value
    else {
       jparc_flux_driver->SetNumOfCycles(gOptFluxNCycles);
    }
  }

  // *************************************************************************
  // * Prepare for writing the output event tree & status file
  // *************************************************************************

  // Initialize an Ntuple Writer to save GHEP records into a TTree
  NtpWriter ntpw(kDefOptNtpFormat, gOptRunNu, gOptRanSeed);
  ntpw.CustomizeFilenamePrefix(gOptEvFilePrefix);
  ntpw.Initialize();

  // Add a custom-branch at the standard GENIE event tree so that
  // info on the flux neutrino parent particle can be passed-through
  flux::GJPARCNuFluxPassThroughInfo * flux_info = 0;
  if(!gOptUsingHistFlux) {
    TBranch * flux = ntpw.EventTree()->Branch("flux",
       "genie::flux::GJPARCNuFluxPassThroughInfo", &flux_info, 32000, 1);
    assert(flux);
    flux->SetAutoDelete(kFALSE);
  }

  // Create a MC job monitor for a periodically updated status file
  GMCJMonitor mcjmonitor(gOptRunNu);
  mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());

  // *************************************************************************
  // * Event generation loop
  // *************************************************************************

  int ievent = 0;
  while (1)
  {
     LOG("gevgen_t2k", pNOTICE)
          << " *** Generating event............ " << ievent;

     // In case the required statistics was expressed as 'number of events'
     // then quit if that number has been generated
     if(ievent == gOptNev) break;

     // In case the required statistics was expressed as 'number of POT' and
     // the user does not want to wait till the end of the flux cycle to exit
     // the event loop, then quit if the requested POT has been generated.
     // In this case the computed POT may not be as accurate as if the program
     // was waiting for the current flux cycle to be completed.
     if(!gOptExitAtEndOfFullFluxCycles && gOptPOT>0) {
        double fpot = jparc_flux_driver->POT_curravg(); // current POT in flux file
        double psc  = mcj_driver->GlobProbScale();      // interaction prob. scale
        double pot  = fpot / psc;                       // POT for generated sample
        if(pot >= gOptPOT) break;
     }

     // Generate a single event using neutrinos coming from the specified flux
     // and hitting the specified geometry or target mix
     EventRecord * event = mcj_driver->GenerateEvent();

     // Check whether a null event was returned due to the flux driver reaching
     // the end of the input flux ntuple - exit the event generation loop
     if(!event && jparc_flux_driver->End()) {
        LOG("gevgen_t2k", pNOTICE)
          << "** The JPARC flux driver read all the input flux ntuple entries";
        break;
     }
     if(!event) {
         LOG("gevgen_t2k", pERROR)
             << "Got a null generated neutino event! Retrying ...";
         continue;
     }
     LOG("gevgen_t2k", pINFO)
         << "Generated event: " << *event;

     // A valid event was generated: extract flux info (parent decay/prod
     // position/kinematics) for that simulated event so that it can be
     // passed-through.
     // Can only do so if I am generating events using the JNUBEAM flux
     // ntuples, not simple histograms
     if(!gOptUsingHistFlux) {
        flux_info = new flux::GJPARCNuFluxPassThroughInfo(
            jparc_flux_driver->PassThroughInfo());
        LOG("gevgen_t2k", pINFO)
          << "Pass-through flux info associated with generated event: "
          << *flux_info;
     }

     // Add event at the output ntuple, refresh the mc job monitor & clean-up
     ntpw.AddEventRecord(ievent, event);
     mcjmonitor.Update(ievent,event);
     delete event;
     if(flux_info) delete flux_info;
     ievent++;
  } //1

  LOG("gevgen_t2k", pNOTICE)
    << "The GENIE MC job is done generaing events - Cleaning up & exiting...";

  // *************************************************************************
  // * Print job statistics &
  // * calculate normalization factor for the generated sample
  // *************************************************************************
  if(!gOptUsingHistFlux && gOptUsingRootGeom)
  {
    // POT normalization will only be calculated if event generation was based
    // on beam simulation  outputs (not just histograms) & a detailed detector
    // geometry description.
    double fpot = jparc_flux_driver->POT_curravg(); // current POT in flux file
    double psc  = mcj_driver->GlobProbScale();      // interaction prob. scale
    double pot  = fpot / psc;                       // POT for generated sample
    // Get nunber of flux neutrinos read-in by flux friver, number of flux
    // neutrinos actually thrown to the event generation driver and number
    // of neutrino interactions actually generated
    long int nflx_evg = mcj_driver        -> NFluxNeutrinos();
    long int nflx     = jparc_flux_driver -> NFluxNeutrinos();
    long int nev      = ievent;

    LOG("gevgen_t2k", pNOTICE)
        << "\n >> Actual JNUBEAM flux file normalization:  " << fpot
            << " POT * " << ((gOptDetectorLocation == "sk") ? "cm^2" : "det")
        << "\n >> Interaction probability scaling factor:  " << psc
        << "\n >> N of flux v read-in by flux driver:      " << nflx
        << "\n >> N of flux v thrown to event gen driver:  " << nflx_evg
        << "\n >> N of generated v interactions:           " << nev
        << "\n ** Normalization for generated sample:      " << pot
            << " POT * " << ((gOptDetectorLocation == "sk") ? "cm^2" : "det");

    ntpw.EventTree()->SetWeight(pot); // POT
  }

  // *************************************************************************
  // * MC job meta-data
  // *************************************************************************

  genie::utils::T2KEvGenMetaData * metadata = new genie::utils::T2KEvGenMetaData;

  metadata -> jnubeam_file       = gOptFluxFile;
  metadata -> detector_location  = gOptDetectorLocation;
  metadata -> geom_file          = gOptRootGeom;
  metadata -> geom_top_volume    = gOptRootGeomTopVol;
  metadata -> geom_length_units  = gOptGeomLUnits;
  metadata -> geom_density_units = gOptGeomDUnits;
  metadata -> using_root_geom    = gOptUsingRootGeom;
  metadata -> target_mix         = gOptTgtMix;
  metadata -> using_hist_flux    = gOptUsingHistFlux;
  metadata -> flux_hists         = gOptFluxHst;

  ntpw.EventTree()->GetUserInfo()->Add(metadata);

  // *************************************************************************
  // * Save & clean-up
  // *************************************************************************

  // Save the generated event tree & close the output file
  ntpw.Save();

  // Clean-up
  delete geom_driver;
  delete flux_driver;
  delete mcj_driver;
  map<int,TH1D*>::iterator it = gOptFluxHst.begin();
  for( ; it != gOptFluxHst.end(); ++it) {
    TH1D * spectrum = it->second;
    if(spectrum) delete spectrum;
  }
  gOptFluxHst.clear();

  LOG("gevgen_t2k", pNOTICE) << "Done!";

  return 0;
}

void PrintSyntax

(

void

)

Definition at line 1462 of file gT2KEvGen.cxx.

{
  LOG("gevgen_t2k", pFATAL)
   << "\n **Syntax**"
   << "\n gevgen_t2k [-h] "
   << "\n           [-r run#]"
   << "\n            -f flux"
   << "\n            -g geometry"
   << "\n           [-p pot_normalization_of_flux_file]"
   << "\n           [-t top_volume_name_at_geom]"
   << "\n           [-P pre_gen_prob_file]"
   << "\n           [-S] [output_name]"
   << "\n           [-m max_path_lengths_xml_file]"
   << "\n           [-L length_units_at_geom]"
   << "\n           [-D density_units_at_geom]"
   << "\n           [-n n_of_events]"
   << "\n           [-c flux_cycles]"
   << "\n           [-e, -E exposure_in_POTs]"
   << "\n           [-o output_event_file_prefix]"
   << "\n           [-R]"
   << "\n           [--seed random_number_seed]"
   << "\n            --cross-sections xml_file"
   << "\n           [--event-generator-list list_name]"
   << "\n           [--message-thresholds xml_file]"
   << "\n           [--unphysical-event-mask mask]"
   << "\n           [--event-record-print-level level]"
   << "\n           [--mc-job-status-refresh-rate  rate]"
   << "\n           [--cache-file root_file]"
   << "\n"
   << " Please also read the detailed documentation at http://www.genie-mc.org"
   << " or look at the source code: $GENIE/src/Apps/gT2KEvGen.cxx"
   << "\n";
}

Variable Documentation

string gOptDetectorLocation

Definition at line 495 of file gT2KEvGen.cxx.

string gOptEvFilePrefix

Definition at line 502 of file gT2KEvGen.cxx.

bool gOptExitAtEndOfFullFluxCycles

Definition at line 501 of file gT2KEvGen.cxx.

string gOptExtMaxPlXml

Definition at line 493 of file gT2KEvGen.cxx.

string gOptFluxFile

Definition at line 494 of file gT2KEvGen.cxx.

map<int,TH1D*> gOptFluxHst

Definition at line 488 of file gT2KEvGen.cxx.

int gOptFluxNCycles

Definition at line 498 of file gT2KEvGen.cxx.

double gOptFluxNorm

Definition at line 496 of file gT2KEvGen.cxx.

PDGCodeList gOptFluxNtpNuList(false)

string gOptFluxProbFileName

Definition at line 505 of file gT2KEvGen.cxx.

double gOptGeomDUnits = 0

Definition at line 492 of file gT2KEvGen.cxx.

double gOptGeomLUnits = 0

Definition at line 491 of file gT2KEvGen.cxx.

string gOptInpXSecFile

Definition at line 509 of file gT2KEvGen.cxx.

int gOptNev

Definition at line 499 of file gT2KEvGen.cxx.

double gOptPOT

Definition at line 500 of file gT2KEvGen.cxx.

bool gOptRandomFluxOffset = false

Definition at line 507 of file gT2KEvGen.cxx.

long int gOptRanSeed

Definition at line 508 of file gT2KEvGen.cxx.

string gOptRootGeom

Definition at line 489 of file gT2KEvGen.cxx.

string gOptRootGeomTopVol = ""

Definition at line 490 of file gT2KEvGen.cxx.

Long_t gOptRunNu

Definition at line 484 of file gT2KEvGen.cxx.

bool gOptSaveFluxProbsFile = false

Definition at line 504 of file gT2KEvGen.cxx.

string gOptSaveFluxProbsFileName

Definition at line 506 of file gT2KEvGen.cxx.

map<int,double> gOptTgtMix

Definition at line 487 of file gT2KEvGen.cxx.

bool gOptUseFluxProbs = false

Definition at line 503 of file gT2KEvGen.cxx.

bool gOptUsingHistFlux = false

Definition at line 486 of file gT2KEvGen.cxx.

bool gOptUsingRootGeom = false

Definition at line 485 of file gT2KEvGen.cxx.

string kDefOptEvFilePrefix = "gntp"

Definition at line 480 of file gT2KEvGen.cxx.

double kDefOptFluxNorm = 1E+21

Definition at line 479 of file gT2KEvGen.cxx.

string kDefOptGeomDUnits = "g_cm3"

Definition at line 477 of file gT2KEvGen.cxx.

string kDefOptGeomLUnits = "mm"

Definition at line 476 of file gT2KEvGen.cxx.

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 478 of file gT2KEvGen.cxx.