gT2KEvGen.cxx

Go to the documentation of this file.

//________________________________________________________________________________________
/*!

\program gevgen_t2k

\brief   A GENIE event generation driver 'customized' for T2K.

         This driver can handle the JPARC neutrino flux files generated by JNUBEAM and the
         realistic detector geometries / target mix of T2K detectors.

         T2K users should note that the generic GENIE event generation driver (gevgen) may
         still be a more appropriate tool to use for the simpler event generation cases
         required for many 4-vector level / systematic studies.
         Please see the GENIE documentation (http://www.genie-mc.org) and contact me
         <constantinos.andreopoulos \at cern.ch> if in doubt.

         *** Synopsis :

           gevgen_t2k [-h]
                      [-r run#]
                       -f flux
                       -g geometry
                      [-p pot_normalization_of_flux_file]
                      [-t top_volume_name_at_geom || -t +Vol1-Vol2...]
                      [-P pre_gen_prob_file_name]
                      [-S] [output_name]
                      [-m max_path_lengths_xml_file]
                      [-L length_units_at_geom]
                      [-D density_units_at_geom]
                      [-n n_of_events]
                      [-c flux_cycles]
                      [-e, -E exposure_in_POTs]
                      [-o output_event_file_prefix]
                      [-R]
                      [--seed random_number_seed]
                       --cross-sections xml_file
                      [--tune genie_tune]
                      [--message-thresholds xml_file]
                      [--unphysical-event-mask mask]
                      [--event-record-print-level level]
                      [--mc-job-status-refresh-rate  rate]
                      [--cache-file root_file]

         *** Options :

           [] Denotes an optional argument

           -h
              Prints out the gevgen_t2k syntax and exits.
           -r
              Specifies the MC run number [default: 1000].
           -g
              Input 'geometry'.
              This option can be used to specify any of:
              1 > A ROOT file containing a ROOT/GEANT geometry description
                  [Note]
                  - This is the standard option for generating events in the
                    nd280, 2km and INGRID detectors.
                  [Examples]
                  - To use the master volume from the ROOT geometry stored
                    in the nd280-geom.root file, type:
                    '-g /some/path/nd280-geom.root'
              2 > A mix of target materials, each with its corresponding weight,
                  typed as a comma-separated list of nuclear PDG codes (in the
                  std PDG2006 convention: 10LZZZAAAI) with the weight fractions
                  in brackets, eg code1[fraction1],code2[fraction2],...
                  If that option is used (no detailed input geometry description)
                  then the interaction vertices are distributed in the detector
                  by the detector MC.
                  [Note]
                  - This is the standard option for generating events in the
                    SuperK detector.
                  [Examples]
                  - To use a target mix of 95% O16 and 5% H type:
                    '-g 1000080160[0.95],1000010010[0.05]'
                  - To use a target which is 100% C12, type:
                    '-g 1000060120'
           -t
              Input 'top volume' for event generation -
              can be used to force event generation in given sub-detector.
              [default: the 'master volume' of the input geometry]
              You can also use the -t option to switch generation on/off at
              multiple volumes as, for example, in:
              `-t +Vol1-Vol2+Vol3-Vol4',
              `-t "+Vol1 -Vol2 +Vol3 -Vol4"',
              `-t -Vol2-Vol4+Vol1+Vol3',
              `-t "-Vol2 -Vol4 +Vol1 +Vol3"'m
              where:
              "+Vol1" and "+Vol3" tells GENIE to `switch on'  Vol1 and Vol3, while
              "-Vol2" and "-Vol4" tells GENIE to `switch off' Vol2 and Vol4.
              If the very first character is a '+', GENIE will neglect all volumes
              except the ones explicitly turned on. Vice versa, if the very first
              character is a `-', GENIE will keep all volumes except the ones
              explicitly turned off (feature contributed by J.Holeczek).
           -P
              Use exact interaction probabilities (stored in the file specified
              via this option) corresponding to the flux file input in this job.
              For complex geometries this will dramatically speed up event generation!
              If this option is chosen then no max path length file needs to be provided.
              Instead of calculating the interaction probabilities on the fly
              per job you can pre-generate them using a dedicated job (see the
              -S option) and tell the event generator to use them via the -P option.
              This is advised when processing flux files with more than ~100k entries
              as the time to pre-calculate the interaction probabilities becomes
              comparable to the event generation time. For smaller flux files there is
              less book-keeping if you just calculate them per job and on the fly.
           -S [output_name]
              Pre-generate flux interaction probabilities and save to root
              output file for use with future event generation jobs. With this
              option no events are generated, it is just a preparatory step
              before actual event generation. When using this option care
              should be taken to run with same arguments used for the actual
              event generation job (same ROOT geomtery, flux file, probe
              particle  types, etc...).
              The output root file is specified as the input for event
              generation using the [pre_gen_prob_file] optional value of the
              -P option. The default output interaction probabilities file
              name is constructed as: [FLUXFILENAME].[TOPVOL].flxprobs.root.
              Specifying [output_name] will override this.
              Introducing multiple functionality to the executable is not
              desirable but is less error prone than duplicating a lot of the
              functionality in a separate application.
           -m
              An XML file (generated by gmxpl) with the max (density weighted)
              path-lengths for each target material in the input ROOT geometry.
              If no file is input, then the geometry will be scanned at MC job
              initialization to determine those max path lengths.
              Supplying this file can speed-up the MC job initialization.
           -L
              Input geometry length units, eg 'm', 'cm', 'mm', ...
              [default: 'mm']
              Note that typically:
               - nd280m uses: 'mm'
               - ...
           -D
              Input geometry density units, eg 'g_cm3', 'clhep_def_density_unit',...
              [default: 'g_cm3']
              Note that typically:
               - nd280m uses: 'clhep_def_density_unit'
               - ...
           -f
              Input 'neutrino flux'.
              This option can be used to specify any of:
              1 > A JNUBEAM beam simulation output file and the detector location.
                  The general sytax is:
                      -f /full/path/flux_file.root,detector_location(,list_of_neutrino_codes)
                  [Notes]
                  - For more information on the flux ntuples, see (T2K internal):
                    http://jnusrv01.kek.jp/internal/t2k/nubeam/flux/
                  - The original HBOOK JNUBAM ntuples need to be converted to a ROOT
                    format using the h2root ROOT utility.
                  - The detector location can be any of 'sk' or the near detector
                    positions 'nd1',...,'nd6','nd13' simulated in JNUBEAM.
                    See the above JNUBEAM web page for more info.
                  - The neutrino codes are the PDG ones.
                  - The (comma separated) list of neutrino codes is optional.
                    It may be used for considering only certain neutrino flavours
                    (eg. nu_e only).
                    If no neutrino list is specified then GENIE will consider all
                    possible flavours (nu_e, nu_e_bar, nu_mu, nu_mu_bar).
                    See examples below.
                  - The JNUBEAM flux ntuples are read via GENIE's GJPARCNuFlux
                    driver. This customized GENIE event generation driver
                    passes-through the complete JPARC input flux information
                    (eg parent decay kinematics / position etc) for each neutrino
                    event it generates (an additional 'flux' branch is added at
                    the output event tree).
                  [Examples]
                  - To use the SuperK flux ntuple from the flux.root file,
                    type:
                       '-f /path/flux.root,sk'
                  - To do the same as above, but considering only nu_mu and
                    nu_mu_bar, type:
                       '-f /path/flux.root,sk,14,-14'
                  - To use the 2km flux ntuple [near detector position 'nd1'
                    in the JNUBEAM flux simulation] from the flux.root file,
                    type:
                       '-f /path/flux.root,nd1'
                  - To use the nd280 flux ntuple [near detector position 'nd5'
                    in the JNUBEAM flux simulation] from the flux.root file,
                    type:
                       '-f /path/flux.root,nd5'
                  - To do the same as above, but considering only nu_e
                    type:
                       '-f /path/flux.root,nd5,12'
              2 > A list of JNUBEAM beam simulation output files and the detector location.
                  The general sytax is:
                      -f /full/path/flux_file_prefix@first_file_number@last_file_number,detector_location(,list_of_neutrino_codes)
                  [Notes]
                  - The ".root" is assumed.
                  - All the files in the series between flux_file_prefix.first_file_number.root
                    and flux_file_prefix.last_file_number.root should exist.
                  - It is important that you set the -p option correctly. See note below.
                  - Also see notes from option 1.
                  [Examples]
                  - To use the SuperK flux ntuples from the files: flux.0.root, flux.1.root, flux.2.root, flux.3.root
                    type:
                       '-f /path/flux.@0@3,sk'
                  - To use the nd280 flux ntuple [near detector position 'nd5'
                    in the JNUBEAM flux simulation] from the files in the series
                    flux.0.root file to flux.100.root, considering only nu_e, type:
                       '-f /path/flux.@0@100,nd5,12'
              3 > A set of histograms stored in a ROOT file.
                  The general syntax is:
                      -f /path/histogram_file.root,neutrino_code[histo_name],...
                  [Notes]
                  - The neutrino codes are the PDG ones.
                  - The 'neutrino_code[histogram_name]' part of the option can be
                    repeated multiple times (separated by commas), once for each
                    flux neutrino species you want to consider, eg
                    '-f somefile.root,12[nuehst],-12[nuebarhst],14[numuhst]'
                  - When using flux from histograms then there is no point in using
                    a 'detailed detector geometry description' as your flux input
                    contains no directional information for those flux neutrinos.
                    The neutrino direction is conventionally set to be +z {x=0,y=0}.
                    So, when using this option you must be using a simple 'target mix'
                    See the -g option for possible geometry settings.
                    If you want to use the detailed detector geometry description
                    then you should be feeding this driver with the JNUBEAM flux
                    simulation outputs.
                  - When using flux from histograms there is no branch with neutrino
                    parent information (code, decay mode, 4-momentum at prod & decay)
                    added in the output event tree as your flux input contains no
                    such information. If you want to be getting the neutrino parent
                    information written out as well then you should be feeding this
                    driver with the JNUBEAM flux simulation outputs.
                  - Note that the relative normalization of the flux histograms is
                    taken into account and is reflected in the relative frequency
                    of flux neutrinos thrown by the flux driver
                  [Examples]
                  - To use the histogram 'h100' (representing the nu_mu flux) and
                    the histogram 'h300' (representing the nu_e flux) and the
                    histogram 'h301' (representing the nu_e_bar flux) from the
                    flux.root file in /path/
                    type:
                      '-f /path/flux.root,14[h100],12[h300],-12[h301]
           -p
              POT normalization of the input flux file.
              [default: The 'standard' JNUBEAM flux ntuple normalization of
               1E+21 POT/detector for the near detectors and
               1E+21 POT/cm2 for the far detector]
              That will be used to interpret the flux weights & calculate the actual
              POT normalization for the generated event sample.
              [Note]
              - If you are using the multiple JNUBEAM flux file entry method it is
                very important that you set this. It should be set to the total POT
                of all input flux files.
              [Examples]
              - If you have 10 standard JNUBEAM files use '-p 10E+21'
           -c
              Specifies how many times to cycle a JNUBEAM flux ntuple.
              Due to the large rejection factor when generating unweighted events
              in the full energy range (approximately ~500 flux neutrinos will be
              rejected before getting an interaction in nd280), an option is
              provided to recycle the flux ntuples for a number of times.
              That option can be used to boost the generated statistics without
              requiring enormous flux files.
              See also 'Note on exposure / statistics' below.
           -e
              Specifies how many POTs to generate.
              If that option is set, gevgen_t2k will work out how many times it has
              to cycle through the input flux ntuple in order to accumulate the
              requested statistics. The program will stop at the earliest complete
              flux ntuple cycle after accumulating the required statistics, so the
              actual statistics will 'slightly' overshoot that number.
              See also 'Note on exposure / statistics' below.
           -E
              Specifies how many POTs to generate.
              That option is similar to -e but the program will stop immediatelly
              after the requested POT has been accumulated. That reduces the
              generated POT overshoot of the requested POT, but the POT calculation
              may not be as exact as with -e.
              See also 'Note on exposure / statistics' below.
           -n
             Specifies how many events to generate.

             --------------------------
             [Note on setting the exposure / statistics]
              All -c, -e (-E) and -n options can be used to set the exposure.
              - If the input flux is a JNUBEAM ntuple then any of these options can
                be used (only one at a time).
                If no option is set, then the program will automatically set '-c 1'
              - If the input flux is described with histograms then only the -n
                option is available.
             --------------------------

           -o
              Sets the prefix of the output event file.
              The output filename is built as:
              [prefix].[run_number].[event_tree_format].[file_format]
              The default output filename is:
              gntp.[run_number].ghep.root
              This cmd line arguments lets you override 'gntp'
           -R
              Tell the flux driver to start looping over the flux ntuples with a
              random offset. May be necessary to avoid biases introduced by always
              starting at the same point when using very large flux input files.
           --seed
              Random number seed.
           --cross-sections
              Name (incl. full path) of an XML file with pre-computed
              cross-section values used for constructing splines.
           --tune
              Specifies a GENIE comprehensive neutrino interaction model tune.
              [default: "Default"].
           --message-thresholds
              Allows users to customize the message stream thresholds.
              The thresholds are specified using an XML file.
              See $GENIE/config/Messenger.xml for the XML schema.
           --unphysical-event-mask
              Allows users to specify a 16-bit mask to allow certain types of
              unphysical events to be written in the output file.
              [default: all unphysical events are rejected]
           --event-record-print-level
              Allows users to set the level of information shown when the event
              record is printed in the screen. See GHepRecord::Print().
           --mc-job-status-refresh-rate
              Allows users to customize the refresh rate of the status file.
           --cache-file
              Allows users to specify a cache file so that the cache can be
              re-used in subsequent MC jobs.

         *** Examples:

         (1) shell% gevgen_t2k
                       -r 1001
                       -f /data/t2k/flux/07a/jnubeam001.root,nd5
                       -g /data/t2k/geom/nd280.root
                       -L mm -D clhep_def_density_unit
                       -e 5E+17
                       --cross-sections /data/t2k/xsec/xsec.xml

             Generate events (run number 1001) using the JNUBEAM flux ntuple in
             /data/t2k/flux/07a/jnubeam001.root & picking up the flux entries for
             the detector location nd5 (:nd280m). The job will load the nd280
             detector geometry description from /data/t2k/geom/nd280.root and
             use it thinking that the geometry length unit is 'mm' and the density
             unit is 'clhep_def_density_unit' (g_cm3 / 0.62415185185E+19)
             The job will stop on the first complete flux ntuple cycle after
             generating 5E+17 POT. Pre-computed cross-sections for all relevant
             initial states are loaded from /data/t2k/xsec/xsec.xml.

         (2) shell% gevgen_t2k
                       -r 1001
                       -f /data/t2k/flux/07a/jnubeam001.root,nd5
                       -g /data/t2k/geom/nd280.root
                       -L mm -D clhep_def_density_unit
                       -c 100
                       --cross-sections /data/t2k/xsec/xsec.xml

             As before, but now the job will stop after 100 flux ntuple cycles -
             whatever POT & number of events that may correspond to.

         (3) shell% gevgen_t2k
                       -r 1001
                       -f /data/t2k/flux/07a/jnubeam001.root,nd5
                       -g /data/t2k/geom/nd280.root
                       -L mm -D clhep_def_density_unit
                       -n 100000
                       --cross-sections /data/t2k/xsec/xsec.xml

             As before, but now the job will stop after generating 100000 events -
             whatever POT & number of flux ntuple cycles that may correspond to.

         (4) shell% gevgen_t2k
                       -r 1001
                       -f /data/t2k/flux/07a/jnubeam001.root,nd5,-12,12
                       -g /data/t2k/geom/nd280.root
                       -L mm -D clhep_def_density_unit
                       -n 100000
                       --cross-sections /data/t2k/xsec/xsec.xml

             As before, but now the job will consider flux nu_e and nu_e_bar only!

         (5) shell% gevgen_t2k
                      -r 1001
                      -f /data/t2k/flux/07a/jnubeam001.root,sk
                      -g 1000080160[0.95],1000010010[0.05]
                      -n 50000
                       --cross-sections /data/t2k/xsec/xsec.xml

             Generate events (run number 1001) using the JNUBEAM flux ntuple in
             /data/t2k/flux/07a/jnubeam001.root & picking up the flux entries for
             the SuperK detector location. This time, the job will not use any
             detailed detector geometry description but just (95% O16 + 5% H)
             target mix. The job will stop after generating 50000 events.

         (6) shell% gevgen_t2k
                       -r 1001
                       -f /data/t2k/flux/hst/flux.root,12[h100],-12[h101],14[h200]
                       -g 1000080160[0.95],1000010010[0.05]
                       -n 50000
                       --cross-sections /data/t2k/xsec/xsec.xml

             As before, but in this case the flux description is not based on a JNUBEAM
             ntuple but a set of histograms at the /data/t2k/flux/hst/flux.root file:
             The histogram named 'h100' will be used for the nu_e flux, 'h101' will
             will be used for the nu_e_bar flux, and 'h200' for the nu_mu flux.


         Please read the GENIE User Manual for more information.

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

\created February 05, 2008

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

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

#ifdef __GENIE_FLUX_DRIVERS_ENABLED__
#include "Tools/Flux/GJPARCNuFlux.h"
#include "Tools/Flux/GCylindTH1Flux.h"
#endif

#ifdef __GENIE_GEOM_DRIVERS_ENABLED__
#include "Tools/Geometry/GeoUtils.h"
#include "Tools/Geometry/ROOTGeomAnalyzer.h"
#include "Tools/Geometry/PointGeomAnalyzer.h"
#endif

using std::string;
using std::vector;
using std::map;
using std::ostringstream;

using namespace genie;

void GetCommandLineArgs (int argc, char ** argv);
void PrintSyntax        (void);

// Default options (override them using the command line arguments):
//
string          kDefOptGeomLUnits   = "mm";    // default geometry length units
string          kDefOptGeomDUnits   = "g_cm3"; // default geometry density units
NtpMCFormat_t   kDefOptNtpFormat    = kNFGHEP; // default event tree format
double          kDefOptFluxNorm     = 1E+21;   // std JNUBEAM flux ntuple norm. (POT*detector [NDs] or POT*cm^2 [SK])
string          kDefOptEvFilePrefix = "gntp";

// User-specified options:
//
Long_t          gOptRunNu;                     // run number
bool            gOptUsingRootGeom = false;     // using root geom or target mix?
bool            gOptUsingHistFlux = false;     // using JNUBEAM flux ntuples or flux from histograms?
map<int,double> gOptTgtMix;                    // target mix  (tgt pdg -> wght frac) / if not using detailed root geom
map<int,TH1D*>  gOptFluxHst;                   // flux histos (nu pdg  -> spectrum)  / if not using beam sim ntuples
string          gOptRootGeom;                  // input ROOT file with realistic detector geometry
string          gOptRootGeomTopVol = "";       // input geometry top event generation volume
double          gOptGeomLUnits = 0;            // input geometry length units
double          gOptGeomDUnits = 0;            // input geometry density units
string          gOptExtMaxPlXml;               // max path lengths XML file for input geometry
string          gOptFluxFile;                  // ROOT file with JNUBEAM flux ntuple
string          gOptDetectorLocation;          // detector location ('sk','nd1','nd2',...)
double          gOptFluxNorm;                  // JNUBEAM flux ntuple normalization
PDGCodeList     gOptFluxNtpNuList(false);      // JNUBEAM flux ntuple neutrinos to consider (can be used for considering certain flavours only)
int             gOptFluxNCycles;               // number of flux ntuple cycles
int             gOptNev;                       // number of events to generate
double          gOptPOT;                       // exposure (in POT)
bool            gOptExitAtEndOfFullFluxCycles; // once POT >= requested_POT, stop at once or at the end of the flux cycle?
string          gOptEvFilePrefix;              // event file prefix
bool            gOptUseFluxProbs = false;      // use pre-calculated flux interaction probs instead of estimating them using the max paths
bool            gOptSaveFluxProbsFile = false; // special mode: no events generated, calculate and save flux interaction probs to root file
string          gOptFluxProbFileName;          // filename for file containg flux probs
string          gOptSaveFluxProbsFileName;     // output filename for pre-generated flux probabilities
bool            gOptRandomFluxOffset = false;  // start looping over flux file from random start entry
long int        gOptRanSeed;                   // random number seed
string          gOptInpXSecFile;               // cross-section splines

//____________________________________________________________________________
int main(int argc, char ** argv)
{
  // 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 GetCommandLineArgs(int argc, char ** argv)
{
  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();
}
//____________________________________________________________________________
void PrintSyntax(void)
{
  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";
}
//____________________________________________________________________________