#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"

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.

References genie::CmdLnArgParser::ArgAsDouble(), genie::CmdLnArgParser::ArgAsInt(), genie::CmdLnArgParser::ArgAsLong(), genie::CmdLnArgParser::ArgAsString(), genie::RunOpt::EnableBareXSecPreCalc(), exit(), genie::PDGLibrary::Find(), flux, geom(), gOptDetectorLocation, gOptEvFilePrefix, gOptExitAtEndOfFullFluxCycles, gOptExtMaxPlXml, gOptFluxFile, gOptFluxHst, gOptFluxNCycles, gOptFluxNorm, gOptFluxNtpNuList, gOptFluxProbFileName, gOptGeomDUnits, gOptGeomLUnits, gOptInpXSecFile, gOptNev, gOptPOT, gOptRandomFluxOffset, gOptRanSeed, gOptRootGeom, gOptRootGeomTopVol, gOptRunNu, gOptSaveFluxProbsFile, gOptSaveFluxProbsFileName, gOptTgtMix, gOptUseFluxProbs, gOptUsingHistFlux, gOptUsingRootGeom, plot_validation_datamc::help, histo, make_syst_table_plots::ibin, genie::PDGLibrary::Instance(), genie::RunOpt::Instance(), genie::pdg::IsAntiNeutrino(), genie::pdg::IsNeutrino(), kDefOptEvFilePrefix, kDefOptFluxNorm, kDefOptGeomDUnits, kDefOptGeomLUnits, LOG, nutype(), genie::CmdLnArgParser::OptionExists(), plot_validation_datamc::parser, pDEBUG, make_root_from_grid_output::pdg, pFATAL, pINFO, pNOTICE, genie::utils::print::PrintFramedMesg(), PrintSyntax(), genie::PDGCodeList::push_back(), pWARN, genie::RunOpt::ReadFromCommandLine(), genie::utils::str::Split(), genie::utils::units::UnitFromString(), and wgt.

Referenced by main().

 {
   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());
             TH1D * spectrum = new TH1D(
                  tmpname.Data(), ihst->GetName(), ihst->GetNbinsX(),
                  ihst->GetXaxis()->GetXmin(), ihst->GetXaxis()->GetXmax());
             spectrum->SetDirectory(0);
             for(int ibin = 1; ibin <= ihst->GetNbinsX(); ibin++) {
                spectrum->SetBinContent(ibin, ihst->GetBinContent(ibin));
             }
             // 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);
   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 1463 of file gT2KEvGen.cxx.

References LOG, and pFATAL.

Referenced by GetCommandLineArgs().

 {
   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";
 }