gFNALExptEvGen.cxx

Go to the documentation of this file.

//______________________________________________________________________________
/*!

\program gevgen_fnal

\brief   A GENIE event generation driver `customized' for the FNAL neutrino
         experiments using flux drivers for file types used by those expts.
         This program was adapted from gevgen_t2k used at T2K.

         This driver can use either the FNAL-supported neutrino flux ntuples
         (generated by gNuMI, g4numi, g4numi_flugg), or "dk2nu" format,
         or plain flux histograms for all neutrino species you are considering.
         You can specify either a ROOT-based detailed detector geometry
         description or a simple target mix. See below for further details.

         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 casesrequired 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_fnal [-h]
                       [-r run#]
                        -f flux
                        -g geometry
                       [-t top_volume_name_at_geom || -t +Vol1-Vol2...]
                       [-m max_path_lengths_xml_file]
                       [-L length_units_at_geom]
                       [-D density_units_at_geom]
                       [-n n_of_events]
                       [-e exposure_in_POTs]
                       [-o output_event_file_prefix]
                       [-F fid_cut_string]
                       [-S nrays]
                       [-z zmin]
                       [-d debug flags]
                       [--seed random_number_seed]
                        --cross-sections xml_file
                       [--event-generator-list list_name]
                       [--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_fnal 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
                  [Examples]
                  - To use the master volume from the ROOT geometry stored
                    in the /some/path/nova-geom.root file, type:
                    '-g /some/path/nova-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.
                  [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).
           -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]
           -D
              Input geometry density units, eg "g_cm3", "clhep_def_density_unit",...
              [default: g_cm3]
           -f
              Input 'neutrino flux'.
              This option can be used to specify any of:
              1 > A g[4][numi|lbne][_flugg] beam simulation output file
                  and the detector location
                  The general sytax is:
                      -f /full/path/flux_file.root,detector,flavor1,flavor2...
                  [Notes]
                  - For more information on the flux ntuples, see the gNuMI doc.
                  - The original hbook ntuples need to be converted to a ROOT
                    format using the h2root ROOT utility.
                  - If flavors aren't specified then use default (12,-12,14,-14)
                  - See GNuMIFlux.xml or Dk2Nu.xml for all supported detector
                    locations
                  - The g[4][NuMI|LBNE][_flugg] flux ntuples are read via GENIE's
                    GNuMIFlux driver, and dk2nu file via an external
                    product w/ the GDk2NuFlux driver (if it can be found).
                    This customized GENIE event generation driver passes-through
                    the complete gNuMI input flux information (eg parent decay
                    kinematics / position etc) for each neutrino event it
                    generates (as an additional 'flux' branch is added on the
                    output event tree).
                  [Examples]
                  - To use the gNuMI flux ntuple flux.root at MINOS near
                    detector location, type:
                     '-f /path/flux.root,MINOS-NearDet'
              1a> Similar to 1 above, but filename contains "dk2nu", then use
                  the GDk2NuFlux driver
              1b> Similar to 1 above, but filename contains "gsimple", then
                  use GSimpleNtpFlux driver
              2 > 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 using a driver that
                    supplies a detailed simulated beam flux.
                  - When using flux from histograms there is no branch with
                    neutrino parent information added in the output event
                    tree as your flux input contains no such information.
                  - 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]

           -e
              Specifies how many POTs to generate.
           -n
              Specifies how many events to generate.

             -------
             [Note on exposure / statistics]
              Both -e and -n options can be used to set the exposure.
              - If the input flux is a non-histogram driver then any of these
                options can be used (one at a time).
              - If the input flux is described with histograms then only the -n
                option is available.
             -------

           -F
              Apply a fiducial cut (for now hard coded ... generalize)
              Only used with ROOTGeomAnalyzer
              if string starts with "-" then reverses sense (ie. anti-fiducial)
           -S
              Number of rays to use to scan geometry for max path length
              Only used with ROOTGeomAnalyzer & { GNuMIFlux, GSimpleNtpFlux, GDk2NuFlux }
              +N  Use flux to scan geometry for max path length
              -N  Use N rays x N points on each face of a box
           -z
              Z from which to start flux ray in user world coordinates
              Only use with ROOTGeomAnalyzer & { GNuMIFlux, GSimpleNtpFlux, GDk2NuFlux }
              If left unset then flux originates on the flux window
              [No longer attempts to determine z from geometry, generally
              got this wrong]
           -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'
           --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_fnal
                       -r 1001
                       -f /data/mc_inputs/flux/flux_00001.root,MINOS-NearDet,12,-12
                       -g /data/mc_inputs/geom/minos.root
                       -L mm -D g_cm3
                       -e 5E+17
                       --cross-sections /data/xsec.xml

             Generate events (run number 1001) using the flux ntuple in
             /data/mc_inputs/flux/v1/flux_00001.root
             The job will load the MINOS near detector detector geometry
             description from /data/mc_inputs/geom/minos.root and interpret it
             using 'mm' as the length unit and 'g/cm^3' as the density unit.
             The job will stop as it accumulates a sample corresponding to
             5E+17 POT.
             Pre-computed cross-section data are loaded from /data/xsec.xml

         (2) shell% gevgen_fnal
                       -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/xsec.xml

         Please read the GENIE user manual for more information.

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

         Robert Hatcher <rhatcher \at fnal.gov>
         Fermi National Accelerator Laboratory

\created August 20, 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 <csignal>

#include <string>
#include <sstream>
#include <vector>
#include <map>
#include <algorithm>  // for transform()
#include <fstream>

#include <TSystem.h>
#include <TError.h>  // for gErrorIgnoreLevel
#include <TTree.h>
#include <TFile.h>
#include <TH1D.h>
#include <TMath.h>
#include <TGeoVolume.h>
#include <TGeoShape.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/Numerical/RandomGen.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/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/PrintUtils.h"
#include "Framework/Utils/SystemUtils.h"

#ifdef __GENIE_FLUX_DRIVERS_ENABLED__
#include "Tools/Flux/GFluxDriverFactory.h"
#include "Tools/Flux/GCylindTH1Flux.h"
#include "Tools/Flux/GFluxFileConfigI.h"
#include "Tools/Flux/GFluxExposureI.h"

//#include "Tools/Flux/GNuMIFlux.h"
//#include "Tools/Flux/GSimpleNtpFlux.h"
//  #ifdef __DK2NU_FLUX_DRIVER_AVAILABLE__
//    #include "dk2nu/tree/dk2nu.h"
//    #include "dk2nu/tree/dkmeta.h"
//    #include "dk2nu/tree/NuChoice.h"
//    #include "dk2nu/genie/GDk2NuFlux.h"
//  #endif

#endif

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

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

using namespace genie;

// Forward function declarations
//
void LoadExtraOptions   (void);
void GetCommandLineArgs (int argc, char ** argv);
void PrintSyntax        (void);
void CreateFidSelection (string fidcut, GeomAnalyzerI* geom_driver);
void CreateRockBoxSelection (string fidcut, GeomAnalyzerI* geom_driver);
void DetermineFluxDriver(string fopt);
void ParseFluxHst       (string fopt);
void ParseFluxFileConfig(string fopt);

// 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
string          kDefOptEvFilePrefix = "gntp";

// User-specified options:
//
Long_t          gOptRunNu;                     // run number
bool            gOptUsingRootGeom = false;     // using root geom or target mix?
bool            gOptUsingHistFlux = false;     // using beam flux ntuples or flux from histograms?
string          gOptFluxDriver = "";           // flux driver class to use
map<string,string> gOptFluxShortNames;
PDGCodeList     gOptFluxPdg;                   // list of neutrino flavors to accept
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
string          gOptRootGeomMasterVol = "";    // (highest level of geometry)
double          gOptGeomLUnits = 0;            // input geometry length units
double          gOptGeomDUnits = 0;            // input geometry density units
string          gOptExtMaxPlXml = "";          // max path lengths XML file for input geometry
bool            gOptWriteMaxPlXml = false;     // rather than read file, write the file
                                               //   triggered by leading '+' on given ofilename
string          gOptFluxFile;                  // ROOT file with beam flux ntuple
string          gOptDetectorLocation;          // detector location (see GNuMIFlux.xml for supported locations))
int             gOptNev;                       // number of events to generate
double          gOptPOT;                       // exposure (in POT)
string          gOptFidCut;                    // fiducial cut selection
int             gOptNScan = 0;                 // # of geometry scan rays
double          gOptZmin = -2.0e30;            // starting z position [ if abs() < 1e30 ]
string          gOptEvFilePrefix;              // event file prefix
int             gOptDebug = 0;                 // debug flags
long int        gOptRanSeed;                   // random number seed
string          gOptInpXSecFile;               // cross-section splines

bool            gSigTERM = false;              // was TERM signal sent?

static void gsSIGTERMhandler(int /* s */)
{
  gSigTERM = true;
  std::cerr << "Caught SIGTERM" << std::endl;
}

//____________________________________________________________________________
int main(int argc, char ** argv)
{
  LoadExtraOptions();
  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, false);

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

  // *************************************************************************
  // * Create / configure the geometry driver
  // *************************************************************************
  GeomAnalyzerI * geom_driver = 0;

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

    // creating & configuring a root geometry driver
    geometry::ROOTGeomAnalyzer * rgeom =
            new geometry::ROOTGeomAnalyzer(gOptRootGeom);
    TGeoVolume * topvol = rgeom->GetGeometry()->GetTopVolume();
    if ( ! topvol ) {
      LOG("gevgen_fnal", pFATAL) << "Null top ROOT geometry volume!";
      exit(1);
    }
    // retrieve the master volume name
    gOptRootGeomMasterVol = topvol->GetName();

    rgeom -> SetLengthUnits  (gOptGeomLUnits);
    rgeom -> SetDensityUnits (gOptGeomDUnits);
    rgeom -> SetTopVolName   (gOptRootGeomTopVol); // set user defined "topvol"

    // getting the bounding box dimensions along z so as to set the
    // appropriate upstream generation surface for the NuMI flux driver

    // RWH 2010-07-16:  do not try to automatically get zmin from geometry, rather
    // by default let the flux start from the window.  If the user wants to
    // override this then they need to explicitly set a "zmin".   Trying to use
    // the geometry is fraught with problems in local vs. global coordinates and
    // units where it can appear to work in some cases but it actually isn't really
    // universally correct.
    //was// TGeoShape * bounding_box = topvol->GetShape();
    //was// bounding_box->GetAxisRange(3, zmin, zmax);
    //was// zmin *= rgeom->LengthUnits();
    //was// zmax *= rgeom->LengthUnits();

    // 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);

    // user specifid a fiducial volume cut ... parse that out
    if ( gOptFidCut.find("rock") != std::string::npos )
                                 CreateRockBoxSelection(gOptFidCut,rgeom);
    else if ( gOptFidCut != "" ) CreateFidSelection(gOptFidCut,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 =
    genie::flux::GFluxDriverFactory::Instance().GetFluxDriver(gOptFluxDriver);
  if ( ! flux_driver ) {
    // couldn't get the requested flux driver ?
    std::ostringstream s;
    s << "Known FluxDrivers:\n";
    const std::vector<std::string>& known =
      genie::flux::GFluxDriverFactory::Instance().AvailableFluxDrivers();
    std::vector<std::string>::const_iterator itr = known.begin();
    for ( ; itr != known.end(); ++itr ) s << "  " << (*itr) << "\n";
    LOG("gevgen_fnal", pFATAL)
      << "Failed to get any flux driver from GFluxDriverFactory\n"
      << "when using \"" << gOptFluxDriver << "\"\n" << s.str();
    exit(1);
  }

  if ( ! gOptUsingHistFlux ) {
    genie::flux::GFluxFileConfigI* flux_file_config =
      dynamic_cast<genie::flux::GFluxFileConfigI*>(flux_driver);
    if ( ! flux_file_config ) {
      LOG("gevgen_fnal", pFATAL)
        << "Failed to get GFluxFileConfigI driver from GFluxDriverFactory\n"
        << "when using \"" << gOptFluxDriver << "\"";
      exit(1);
    }

    //
    // *** Using the detailed ntuple neutrino flux description
    //
    flux_file_config->LoadBeamSimData(gOptFluxFile, gOptDetectorLocation);
    flux_file_config->SetUpstreamZ(gOptZmin);  // was "zmin" from bounding_box
    flux_file_config->SetNumOfCycles(0);

    if ( gOptFluxPdg.size() > 0 ) {
      // user specified list of neutrino PDGs
      flux_file_config->SetFluxParticles(gOptFluxPdg);
      std::ostringstream s;
      PDGCodeList::const_iterator itr = gOptFluxPdg.begin();
      for ( ; itr != gOptFluxPdg.end(); ++itr) s << (*itr) << " ";
      LOG("gevgen_fnal", pNOTICE)
        << "Limiting to nu PDGs: " << s.str();
    }
  }
  else {
    //
    // *** Using fluxes from histograms (for all specified neutrino species)
    //
    genie::flux::GCylindTH1Flux * hist_flux_driver =
      dynamic_cast<genie::flux::GCylindTH1Flux*>(flux_driver);
    if ( ! hist_flux_driver ) {
      LOG("gevgen_fnal", pFATAL)
        << "Failed to get GCylinderTH1Flux driver from GFluxDriverFactory\n"
        << "when using " << gOptFluxDriver;
      exit(1);
    }

    // creating & configuring a generic GCylindTH1Flux flux driver
    TVector3 bdir (0,0,1); // dir along +z
    TVector3 bspot(0,0,0);

    hist_flux_driver->SetNuDirection      (bdir);
    hist_flux_driver->SetBeamSpot         (bspot);
    hist_flux_driver->SetTransverseRadius (-1);
    map<int,TH1D*>::iterator it = gOptFluxHst.begin();
    for( ; it != gOptFluxHst.end(); ++it) {
        int    pdg_code = it->first;
        TH1D * spectrum = it->second;
        hist_flux_driver->AddEnergySpectrum(pdg_code, spectrum);
        // once the histogram has been added to the GCylindTH1Flux driver
        // it is owned by the driver and it is up to the the driver
        // to clean up (i.e. delete it).
        // remove it from this map to avoid double deletion.
        it->second = 0;
    }
  }

  // these might come in handy ... avoid repeated dynamic_cast<> calls
  genie::flux::GFluxExposureI* fluxExposureI =
    dynamic_cast<genie::flux::GFluxExposureI*>(flux_driver);
  genie::flux::GFluxFileConfigI* fluxFileConfigI =
    dynamic_cast<genie::flux::GFluxFileConfigI*>(flux_driver);


  // *************************************************************************
  // * Handle chicken/egg problem: geom analyzer vs. flux.
  // * Need both at this point change geom scan defaults.
  // *************************************************************************
  if ( gOptUsingRootGeom && ! gOptUsingHistFlux ) {

    geometry::ROOTGeomAnalyzer * rgeom =
      dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
    if ( ! rgeom ) assert(0);

    rgeom -> SetDebugFlags(gOptDebug);

    // even if user doesn't specify gOptNScan configure to scan using flux
    if ( gOptNScan >= 0 ) {
      LOG("gevgen_fnal", pNOTICE)
        << "Using ROOTGeomAnalyzer: geom scan using flux: nparticles=" << gOptNScan;
      rgeom->SetScannerFlux(flux_driver);
      if ( gOptNScan > 0 ) rgeom->SetScannerNParticles(gOptNScan);
    } else {
      int nabs = TMath::Abs(gOptNScan);
      LOG("gevgen_fnal", pNOTICE)
        << "Using ROOTGeomAnalyzer: geom scan using box: npoints=nrays=" << nabs;
      rgeom->SetScannerNPoints(nabs);
      rgeom->SetScannerNRays(nabs);
    }
  }

  // *************************************************************************
  // * 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);
  if ( ( gOptExtMaxPlXml != "" ) && ! gOptWriteMaxPlXml ) {
    mcj_driver->UseMaxPathLengths(gOptExtMaxPlXml);
  }
  mcj_driver->Configure();
  mcj_driver->UseSplines();
  mcj_driver->ForceSingleProbScale();

  if ( ( gOptExtMaxPlXml != "" ) && gOptWriteMaxPlXml ) {
    geometry::ROOTGeomAnalyzer * rgeom =
      dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
    if ( rgeom ) {
      const genie::PathLengthList& maxpath = rgeom->GetMaxPathLengths();
      std::string maxplfile = gOptExtMaxPlXml;
      maxpath.SaveAsXml(maxplfile);
      // append extra info to file
      std::ofstream mpfile(maxplfile.c_str(), std::ios_base::app);
      mpfile
        << std::endl
        << "<!-- this file is only relevant for a setup compatible with:"
        << std::endl
        << "geom: " << gOptRootGeom << " top: \"" << gOptRootGeomTopVol << "\""
        << std::endl
        << "flux: " << gOptFluxFile
        << std::endl
        << "location: " << gOptDetectorLocation
        << std::endl
        << "fidcut: " << gOptFidCut
        << std::endl
        << "nscan: " << gOptNScan << " (0>= use flux, <0 use box |nscan| points/rays)"
        << std::endl
        << "zmin: " << gOptZmin << " (if |zmin| > 1e30, leave ray on flux window)"
        << std::endl
        << "-->"
        << std::endl;
      mpfile.close();
    }
  }

  // *************************************************************************
  // * 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();


  std::vector<TBranch*>    extraBranches;
  std::vector<std::string> branchNames;
  std::vector<std::string> branchClassNames;
  std::vector<void**>      branchObjPointers;

  // Add custom branch(s) to the standard GENIE event tree so that
  // info on the flux neutrino parent particle can be passed-through
  if ( fluxFileConfigI ) {
    fluxFileConfigI->GetBranchInfo(branchNames,branchClassNames,
                                    branchObjPointers);
    size_t nn = branchNames.size();
    size_t nc = branchClassNames.size();
    size_t np = branchObjPointers.size();
    if ( nn != nc || nc != np ) {
     LOG("gevgen_fnal", pERROR)
       << "Inconsistent info back from \"" << gOptFluxDriver << "\" "
       << "for branch info:  " << nn << " " << nc << " " << np;
    } else {
      for (size_t ii = 0; ii < nn; ++ii) {
        const char* bname = branchNames[ii].c_str();
        const char* cname = branchClassNames[ii].c_str();
        void**&     optr  = branchObjPointers[ii];  // note critical '&' !
        if ( ! optr || ! *optr ) continue;  // no pointer supplied, skip it
        int split = 99;  // 1
        LOG("gevgen_fnal", pNOTICE)
          << "Adding extra branch \"" << bname << "\" of type \""
          << cname << "\" (" << optr << ") to output tree";
        TBranch* bptr = ntpw.EventTree()->Branch(bname,cname,optr,32000,split);
        extraBranches.push_back(bptr);

        if ( bptr ) {
          // don't delete this !!! we're sharing
          bptr->SetAutoDelete(false);
        } else {
          LOG("gevgen_fnal", pERROR)
          << "FAILED to add extra branch \"" << bname << "\" of type \""
          << cname << "\" to output tree";
        }
      } // loop over additions
    } // same # of entries
  } // of genie::flux::GFluxFileConfigI type

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

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

  // define handler to allow signal to end job gracefully
  signal(SIGTERM,gsSIGTERMhandler);

  int ievent = 0;
  while ( ! gSigTERM )
  {
     LOG("gevgen_fnal", pINFO)
          << " *** 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'
     // then exit the event loop if the requested POT has been generated.
     if ( gOptPOT > 0 && fluxExposureI ) {
        double fpot = fluxExposureI->GetTotalExposure(); // current POTs used
        double psc  = mcj_driver->GlobProbScale();  // interaction prob. scale
        double pot  = fpot / psc;                   // POTs 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 && flux_driver->End() ) {
        LOG("gevgen_fnal", pWARN)
          << "** The flux driver read all the input flux entries: End()==true";
        break;
     }
     if ( ! event ) {
         LOG("gevgen_fnal", pERROR)
             << "Got a null generated neutino event! Retrying ...";
         continue;
     }
     if ( print_level >= 0 ) {
       LOG("gevgen_fnal", pINFO)
         << "Generated event: " << *event;
     }
     // A valid event was generated: flux info (parent decay/prod
     // position/kinematics) for that simulated event should already
     // be connected to the right output tree branch

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

  } //1

  // Copy metadata tree, if available
  if ( fluxFileConfigI ) {
    TTree* t1 = fluxFileConfigI->GetMetaDataTree();
    if ( t1 ) {
      size_t nmeta = t1->GetEntries();
      TTree* t2 = (TTree*)t1->Clone(0);
      for (size_t i = 0; i < nmeta; ++i) {
        t1->GetEntry(i);
        t2->Fill();
      }
      t2->Write();
    }
  }

  LOG("gevgen_fnal", pINFO)
    << "The GENIE MC job is done generating 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 ntuples (not just histograms) & a detailed detector
    // geometry description.
    // Get nunber of flux neutrinos read-in by flux driver, number of flux
    // neutrinos actually thrown to the event generation driver and number
    // of neutrino interactions actually generated
    long int nflx     = 0;
    long int nflx_evg = mcj_driver-> NFluxNeutrinos();
    double   fpot     = 0;
    const char* exposureUnits = "(unknown units)";
    if ( fluxExposureI ) {
      fpot = fluxExposureI->GetTotalExposure(); // POTs used so far
      nflx = fluxExposureI->NFluxNeutrinos();
      //genie::flux::Exposure_t etype = fluxExposureI->GetExposureType();
      //exposureUnits = genie::flux::GFluxExposureI::AsString(etype);
      exposureUnits = fluxExposureI->GetExposureUnits();
    }
    if ( fluxFileConfigI ) {
      fluxFileConfigI->PrintConfig();
    }
    double psc   = mcj_driver->GlobProbScale();      // interaction prob. scale
    if ( psc <= 0.0 ) {
       LOG("gevgen_fnal", pFATAL) << "MCJobDriver GlobalProbScale was " << psc;
    }
    double pot   = fpot / psc;                       // POT for generated sample
    long int nev = ievent;

    LOG("gevgen_fnal", pNOTICE)
        << "\n >> Interaction probability scaling factor:  " << psc
        << "\n >> using: " << gOptFluxDriver
        << "\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
        << " " << exposureUnits << " * detector";

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

  }

  // *************************************************************************
  // * 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;
  // this list should only be histograms that have (for some reason)
  // not been handed over to the GCylindTH1Flux 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_fnal", pNOTICE) << "Done!";

  return 0;
}

//____________________________________________________________________________
void LoadExtraOptions(void)
{
  /// potentially load extra libraries that might extend the list of
  /// potential flux drivers, and how to map short names to classes ...
  // we really should at this point read some external file to get
  // an expandable list of libraries ... but for now just hard code it

  vector<string> extraLibs;

  ///***** this part should come from reading an external file
  /// placeholder file $GENIE/config/FluxDriverExpansion.xml

  extraLibs.push_back("libdk2nuTree");
  extraLibs.push_back("libdk2nuGenie");

  // what one might expect to find at the beginning of -f <arg>

  gOptFluxShortNames["histo"]   = "genie::flux::GCylindTH1Flux";
  gOptFluxShortNames["hist"]    = "genie::flux::GCylindTH1Flux";

  gOptFluxShortNames["simple"]  = "genie::flux::GSimpleNtpFlux";
  gOptFluxShortNames["numi"]    = "genie::flux::GNuMIFlux";
  gOptFluxShortNames["dk2nu"]   = "genie::flux::GDk2NuFlux";

  ///******* done with fake "read"

  // see if there are any libraries to load
  // just let ROOT do it ... check error code on return
  // but tweak ROOT's ERROR message output so we don't see huge messages
  // for failures
  //  gErrorIgnoreLevel to kError (from TError.h)

  Int_t prevErrorLevel = gErrorIgnoreLevel;
  gErrorIgnoreLevel    = kFatal;
  for (size_t ilib=0; ilib < extraLibs.size(); ++ilib) {
    // library names should be like libXYZZY  without extension [e.g. .so]
    // but with the leading "lib"
    int loadStatus = gSystem->Load(extraLibs[ilib].c_str());
    const char* statWords = "failed to load";
    if      ( loadStatus ==  0 ) { statWords = "successfully loaded"; }
    else if ( loadStatus ==  1 ) { statWords = "already had loaded"; }
    else if ( loadStatus == -1 ) { statWords = "couldn't find library"; }
    else if ( loadStatus == -2 ) { statWords = "mismatched version"; }

    LOG("gevgen_fnal",pNOTICE)
         << statWords << " (" << loadStatus << ") " << extraLibs[ilib];
  }
  // restore the ROOT error message level
  gErrorIgnoreLevel    = prevErrorLevel;

}

//____________________________________________________________________________
void GetCommandLineArgs(int argc, char ** argv)
{
  LOG("gevgen_fnal", 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_fnal", pDEBUG) << "Reading MC run number";
    gOptRunNu = parser.ArgAsLong('r');
  } else {
    LOG("gevgen_fnal", pDEBUG)
      << "Unspecified run number - Using default";
    gOptRunNu = 0;
  } //-r

  //
  // *** geometry
  //

  string geom = "";
  string lunits, dunits;
  if( parser.OptionExists('g') ) {
    LOG("gevgen_fnal", 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_fnal", 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_fnal", pDEBUG)
           << "Checking for input geometry length units";
        lunits = parser.ArgAsString('L');
     } else {
        LOG("gevgen_fnal", pDEBUG) << "Using default geometry length units";
        lunits = kDefOptGeomLUnits;
     } // -L
     // density units:
     if( parser.OptionExists('D') ) {
        LOG("gevgen_fnal", pDEBUG)
           << "Checking for input geometry density units";
        dunits = parser.ArgAsString('D');
     } else {
        LOG("gevgen_fnal", 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_fnal", pDEBUG) << "Checking for input volume name";
        gOptRootGeomTopVol = parser.ArgAsString('t');
     } else {
        LOG("gevgen_fnal", 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 passed name starts with '+', then compute max at job init, but write out the result
     if ( parser.OptionExists('m') ) {
        LOG("gevgen_fnal", pDEBUG)
              << "Checking for maximum path lengths XML file";
        gOptExtMaxPlXml   = parser.ArgAsString('m');
        gOptWriteMaxPlXml = false;
        if ( gOptExtMaxPlXml[0] == '+' ) {
          gOptExtMaxPlXml   = gOptExtMaxPlXml.substr(1,std::string::npos);
          gOptWriteMaxPlXml = true;
          LOG("gevgen_fnal", pINFO)
            << "Will write maximum path lengths XML file: " << gOptExtMaxPlXml;
        }
     } else {
        LOG("gevgen_fnal", pDEBUG)
               << "Will compute the maximum path lengths at job init";
        gOptExtMaxPlXml = "";
     } // -m

     // fidcut:
     if( parser.OptionExists('F') ) {
       LOG("gevgen_fnal", pDEBUG) << "Using Fiducial cut?";
       gOptFidCut = parser.ArgAsString('F');
     } else {
       LOG("gevgen_fnal", pDEBUG) << "No fiducial volume cut";
       gOptFidCut = "";
     } //-F

     if(!gOptUsingHistFlux) {
       // how to scan the geometry (if relevant)
       if( parser.OptionExists('S') ) {
         LOG("gevgen_fnal", pDEBUG)  << "Reading requested geom scan count";
         gOptNScan = parser.ArgAsInt('S');
       } else {
         LOG("gevgen_fnal", pDEBUG) << "No geom scan count was requested";
         gOptNScan = 0;
       } //-S

       // z for flux rays to start
       if( parser.OptionExists('z') ) {
         LOG("gevgen_fnal", pDEBUG)  << "Reading requested zmin";
         gOptZmin = parser.ArgAsDouble('z');
       } else {
         LOG("gevgen_fnal", pDEBUG) << "No zmin was requested";
         gOptZmin = -2.0e30; // < -1.0e30 ==> leave it on flux window
       } //-z

       // debug flags
       if ( parser.OptionExists('d') ) {
         LOG("gevgen_fnal", pDEBUG) << "Reading debug flag value";
         gOptDebug = parser.ArgAsInt('d');
       } else {
         LOG("gevgen_fnal", pDEBUG) << "Unspecified debug flags - Using default";
         gOptDebug = 0;
       } //-d

     } // root geom && gnumi flux

  } // 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_fnal", 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_fnal", 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_fnal", pDEBUG) << "Getting input flux";
    DetermineFluxDriver(parser.ArgAsString('f'));
  } else {
    LOG("gevgen_fnal", pFATAL) << "No flux info was specified - Exiting";
    PrintSyntax();
    exit(1);
  }

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

  // statistics to generate in terms of POT
  if( parser.OptionExists('e') ) {
    LOG("gevgen_fnal", pDEBUG)  << "Reading requested exposure in POT";
    gOptPOT = parser.ArgAsDouble('e');
  } else {
    LOG("gevgen_fnal", pDEBUG) << "No POT exposure was requested";
    gOptPOT = -1;
  } //-e

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


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

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


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

  // Tthe 'exposure' may be set either as:
  // - a number of POTs
  // - a number of generated events
  // Only one of those options can be set.
  if(!gOptUsingHistFlux) {
    int nset=0;
    if(gOptPOT > 0) nset++;
    if(gOptNev > 0) nset++;
    if(nset==0) {
       LOG("gevgen_fnal", pFATAL)
        << "** To use a gNuMI flux ntuple you need to specify an exposure, "
        << "either via the -e or -n options";
       PrintSyntax();
       exit(1);
    }
    if(nset>1) {
       LOG("gevgen_fnal", pFATAL)
         << "You can not specify more than one of the -e or -n options";
       PrintSyntax();
       exit(1);
    }
  }
  // If we use a flux histograms (not flux ntuples) then -currently- the
  // only way to control exposure is via a number of events
  if(gOptUsingHistFlux) {
     if(gOptNev < 0) {
       LOG("gevgen_fnal", 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_fnal", pFATAL)
         << "You may not use the -e option without a detector geometry description";
       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 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 " << gOptFluxDriver << " flux driver- "
             << "file = "        << gOptFluxFile
             << ", location = "  << gOptDetectorLocation;
  }

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


  LOG("gevgen_fnal", pNOTICE)
     << "\n\n"
     << utils::print::PrintFramedMesg("FNAL expt. event generation job configuration");

  LOG("gevgen_fnal", 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_fnal", pNOTICE) << *RunOpt::Instance();
}
//____________________________________________________________________________
void PrintSyntax(void)
{
  LOG("gevgen_fnal", pFATAL)
   << "\n **Syntax**"
   << "\n gevgen_fnal [-h] [-r run#]"
   << "\n            -f flux -g geometry"
   << "\n            [-t top_volume_name_at_geom] [-m max_path_lengths_xml_file]"
   << "\n            [-L length_units_at_geom] [-D density_units_at_geom]"
   << "\n            [-n n_of_events] [-e exposure_in_POTs]"
   << "\n            [-o output_event_file_prefix]"
   << "\n            [-F fid_cut_string] [-S nrays_scan]"
   << "\n            [-z zmin_start]"
   << "\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 "
   << "$GENIE/src/Apps/gFNALExptEvGen.cxx"
   << "\n";
}
//____________________________________________________________________________
void CreateFidSelection (string fidcut, GeomAnalyzerI* geom_driver)
{
  ///
  /// User defined fiducial volume cut
  ///      [0][M]<SHAPE>:val1,val2,...
  ///   "0" means reverse the cut (i.e. exclude the volume)
  ///   "M" means the coordinates are given in the ROOT geometry
  ///       "master" system and need to be transformed to "top vol" system
  ///   <SHAPE> can be any of "zcyl" "box" "zpoly" "sphere"
  ///       [each takes different # of args]
  ///   This must be followed by a ":" and a list of values separated by punctuation
  ///       (allowed separators: commas , parentheses () braces {} or brackets [] )
  ///   Value mapping:
  ///      zcyl:x0,y0,radius,zmin,zmax           - cylinder along z at (x0,y0) capped at z's
  ///      xcyl:y0,z0,radius,xmin,ymax           - cylinder along x
  ///      ycyl:x0,z0,radius,ymin,ymax           - cylinder along y
  ///      gcyl:{x0,y0,z0}{dx,dy,dz},radius,{plane1}{plane2} -- generic cylinder w/ arbitrary orientation and caps
  ///                                              {planeX} = 4 values to define plane and orientation
  ///      box:xmin,ymin,zmin,xmax,ymax,zmax     - box w/ upper & lower extremes
  ///      zpoly:nfaces,x0,y0,r_in,phi,zmin,zmax - nfaces sided polygon in x-y plane
  //       sphere:x0,y0,z0,radius                - sphere of fixed radius at (x0,y0,z0)
  ///   Examples:
  ///      1) 0mbox:0,0,0.25,1,1,8.75
  ///         exclude (i.e. reverse) a box in master coordinates w/ corners (0,0,0.25) (1,1,8.75)
  ///      2) mzpoly:6,(2,-1),1.75,0,{0.25,8.75}
  ///         six sided polygon in x-y plane, centered at x,y=(2,-1) w/ inscribed radius 1.75
  ///         no rotation (so first face is in y-z plane +r from center, i.e. hex sits on point)
  ///         limited to the z range of {0.25,8.75} in the master ROOT geom coordinates
  ///      3) zcyl:(3,4),5.5,-2,10
  ///         a cylinder oriented parallel to the z axis in the "top vol" coordinates
  ///         at x,y=(3,4) with radius 5.5 and z range of {-2,10}
  ///
  geometry::ROOTGeomAnalyzer * rgeom =
    dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
  if ( ! rgeom ) {
    LOG("gevgen_fnal", pWARN)
      << "Can not create GeomVolSelectorFiduction,"
      << " geometry driver is not ROOTGeomAnalyzer";
    return;
  }

  LOG("gevgen_fnal", pNOTICE) << "-F " << fidcut;

  genie::geometry::GeomVolSelectorFiducial* fidsel =
    new genie::geometry::GeomVolSelectorFiducial();

  fidsel->SetRemoveEntries(true);  // drop segments that won't be considered

  // convert string to lowercase
  std::transform(fidcut.begin(),fidcut.end(),fidcut.begin(),::tolower);

  vector<string> strtok = genie::utils::str::Split(fidcut,":");
  if ( strtok.size() != 2 ) {
    LOG("gevgen_fnal", pWARN)
      << "Can not create GeomVolSelectorFiduction,"
      << " no \":\" separating type from values.  nsplit=" << strtok.size();
    for ( unsigned int i=0; i < strtok.size(); ++i )
      LOG("gevgen_fnal",pNOTICE)
        << "strtok[" << i << "] = \"" << strtok[i] << "\"";
    return;
  }

  // parse out optional "x" and "m"
  string stype = strtok[0];
  bool reverse = ( stype.find("0") != string::npos );
  bool master  = ( stype.find("m") != string::npos );  // action after values are set

  // parse out values
  vector<double> vals;
  vector<string> valstrs = genie::utils::str::Split(strtok[1]," ,;(){}[]");
  vector<string>::const_iterator iter = valstrs.begin();
  for ( ; iter != valstrs.end(); ++iter ) {
    const string& valstr1 = *iter;
    if ( valstr1 != "" ) vals.push_back(atof(valstr1.c_str()));
  }
  size_t nvals = vals.size();

  std::cout << "ivals = [";
  for (unsigned int i=0; i < nvals; ++i) {
    if (i>0) cout << ",";
    std::cout << vals[i];
  }
  std::cout << "]" << std::endl;

  // std::vector elements are required to be adjacent so we can treat address as ptr

  if        ( stype.find("zcyl")   != string::npos ) {
    // cylinder along z direction at (x0,y0) radius zmin zmax
    if ( nvals < 5 )
      LOG("gevgen_fnal", pFATAL) << "MakeZCylinder needs 5 values, not " << nvals
                                << " fidcut=\"" << fidcut << "\"";
    fidsel->MakeZCylinder(vals[0],vals[1],vals[2],vals[3],vals[4]);

  } else if ( stype.find("xcyl")   != string::npos ) {
    // cylinder along x direction at (y0,z0) radius xmin xmax
    if ( nvals < 5 )
      LOG("gevgen_fnal", pFATAL) << "MakeXCylinder needs 5 values, not " << nvals
                                << " fidcut=\"" << fidcut << "\"";
    fidsel->MakeXCylinder(vals[0],vals[1],vals[2],vals[3],vals[4]);

  } else if ( stype.find("ycyl")   != string::npos ) {
    // cylinder along y direction at (x0,z0) radius ymin ymax
    if ( nvals < 5 )
      LOG("gevgen_fnal", pFATAL) << "MakeYCylinder needs 5 values, not " << nvals
                                << " fidcut=\"" << fidcut << "\"";
    fidsel->MakeYCylinder(vals[0],vals[1],vals[2],vals[3],vals[4]);

  } else if ( stype.find("gcyl")   != string::npos ) {
    // cylinder along arbitrary direction at (x0,y0,z0) radius {plane1} {plane2}
    if ( nvals < 14 )
      LOG("gevgen_fnal", pFATAL) << "MakeYCylinder needs 14 values, not " << nvals
                                << " fidcut=\"" << fidcut << "\"";
    Double_t base[3] = { vals[0], vals[1], vals[2] };
    Double_t axis[3] = { vals[3], vals[4], vals[5] };
    Double_t radius = vals[6];
    Double_t cap1[4] = { vals[ 7], vals[ 8], vals[ 9], vals[10] };
    Double_t cap2[4] = { vals[11], vals[12], vals[13], vals[14] };

    fidsel->MakeCylinder(base,axis,radius,cap1,cap2);

  } else if ( stype.find("box")    != string::npos ) {
    // box (xmin,ymin,zmin) (xmax,ymax,zmax)
    if ( nvals < 6 )
      LOG("gevgen_fnal", pFATAL) << "MakeBox needs 6 values, not " << nvals
                                << " fidcut=\"" << fidcut << "\"";
    double xyzmin[3] = { vals[0], vals[1], vals[2] };
    double xyzmax[3] = { vals[3], vals[4], vals[5] };
    fidsel->MakeBox(xyzmin,xyzmax);

  } else if ( stype.find("zpoly")  != string::npos ) {
    // polygon along z direction nfaces at (x0,y0) radius phi zmin zmax
    if ( nvals < 7 )
      LOG("gevgen_fnal", pFATAL) << "MakeZPolygon needs 7 values, not " << nvals
                                << " fidcut=\"" << fidcut << "\"";
    int nfaces = (int)vals[0];
    if ( nfaces < 3 )
      LOG("gevgen_fnal", pFATAL) << "MakeZPolygon needs nfaces>=3, not " << nfaces
                                << " fidcut=\"" << fidcut << "\"";
    fidsel->MakeZPolygon(nfaces,vals[1],vals[2],vals[3],vals[4],vals[5],vals[6]);

  } else if ( stype.find("sphere") != string::npos ) {
    // sphere at (x0,y0,z0) radius
    if ( nvals < 4 )
      LOG("gevgen_fnal", pFATAL) << "MakeZSphere needs 4 values, not " << nvals
                                << " fidcut=\"" << fidcut << "\"";
    fidsel->MakeSphere(vals[0],vals[1],vals[2],vals[3]);

  } else {
    LOG("gevgen_fnal", pFATAL)
      << "Can not create GeomVolSelectorFiduction for shape \"" << stype << "\"";
  }

  if ( master  ) {
    fidsel->ConvertShapeMaster2Top(rgeom);
    LOG("gevgen_fnal", pNOTICE) << "Convert fiducial volume from master to topvol coords";
  }
  if ( reverse ) {
    fidsel->SetReverseFiducial(true);
    LOG("gevgen_fnal", pNOTICE) << "Reverse sense of fiducial volume cut";
  }
  rgeom->AdoptGeomVolSelector(fidsel);

}
//____________________________________________________________________________
void CreateRockBoxSelection (string fidcut, GeomAnalyzerI* geom_driver)
{

  if( fidcut.find_first_not_of(" \t\n") != 0) // trim any leading whitespace
    fidcut.erase( 0, fidcut.find_first_not_of(" \t\n")  );

  // convert string to lowercase
  std::transform(fidcut.begin(),fidcut.end(),fidcut.begin(),::tolower);

  genie::geometry::ROOTGeomAnalyzer * rgeom =
    dynamic_cast<genie::geometry::ROOTGeomAnalyzer *>(geom_driver);
  if ( ! rgeom ) {
    LOG("gevgen_fnal", pWARN)
      << "Can not create GeomVolSelectorRockBox,"
      << " geometry driver is not ROOTGeomAnalyzer";
    return;
  }

  LOG("gevgen_fnal", pWARN) << "fiducial (rock) cut: " << fidcut;

  // for now, only fiducial no "rock box"
  genie::geometry::GeomVolSelectorRockBox* rocksel =
    new genie::geometry::GeomVolSelectorRockBox();

  vector<string> strtok = genie::utils::str::Split(fidcut,":");
  if ( strtok.size() != 2 ) {
    LOG("gevgen_fnal", pWARN)
      << "Can not create GeomVolSelectorRockBox,"
      << " no \":\" separating type from values.  nsplit=" << strtok.size();
    for ( unsigned int i=0; i < strtok.size(); ++i )
      LOG("gevgen_fnal", pWARN)
        << "strtok[" << i << "] = \"" << strtok[i] << "\"";
    return;
  }

  string stype = strtok[0];

  // parse out values
  vector<double> vals;
  vector<string> valstrs = genie::utils::str::Split(strtok[1]," ,;(){}[]\t\n\r");
  vector<string>::const_iterator iter = valstrs.begin();
  for ( ; iter != valstrs.end(); ++iter ) {
    const string& valstr1 = *iter;
    if ( valstr1 != "" ) {
      double aval = atof(valstr1.c_str());
      LOG("gevgen_fnal", pWARN) << "rock value [" << vals.size() << "] "
                                << aval;
      vals.push_back(aval);
    }
  }
  size_t nvals = vals.size();

  rocksel->SetRemoveEntries(true);  // drop segments that won't be considered

  // assume coordinates are in the *master* (not "top volume") system
  // need to set fTopVolume to fWorldVolume
  //fTopVolume = fWorldVolume;
  //rgeom->SetTopVolName(fTopVolume.c_str());
  gOptRootGeomTopVol = gOptRootGeomMasterVol;
  rgeom->SetTopVolName(gOptRootGeomMasterVol);

  if ( nvals < 6 ) {
    LOG("gevgen_fnal", pFATAL)  << "rockbox needs at "
                                << "least 6 values, found "
                                << nvals << "in \""
                                << strtok[1] << "\"";
    exit(1);

  }
  double xyzmin[3] = { vals[0], vals[1], vals[2] };
  double xyzmax[3] = { vals[3], vals[4], vals[5] };

  bool   rockonly  = true;
  double wallmin   = 800.;   // geometry in cm, ( 8 meter buffer)
  double dedx      = 2.5 * 1.7e-3; // GeV/cm, rho=2.5, 1.7e-3 ~ rock like loss
  double fudge     = 1.05;

  if ( nvals >=  7 ) rockonly = vals[6];
  if ( nvals >=  8 ) wallmin  = vals[7];
  if ( nvals >=  9 ) dedx     = vals[8];
  if ( nvals >= 10 ) fudge    = vals[9];

  rocksel->SetRockBoxMinimal(xyzmin,xyzmax);
  rocksel->SetMinimumWall(wallmin);
  rocksel->SetDeDx(dedx/fudge);

  if ( nvals >= 11 ) rocksel->SetExpandFromInclusion((int)vals[10]);

  // if not rock-only then make a tiny exclusion bubble
  // call to MakeBox shouldn't be necessary
  //  should be done by SetRockBoxMinimal but for some GENIE versions isn't
  if ( ! rockonly ) rocksel->MakeSphere(0,0,0,1.0e-10);
  else              rocksel->MakeBox(xyzmin,xyzmax);

  rgeom->AdoptGeomVolSelector(rocksel);

}

//____________________________________________________________________________
void DetermineFluxDriver(string fopt)
{
  // based on the -f option string determine which flux driver to use
  // this may take some guessing

  // first look for strings that look like "<proto>:..."
  // or ....<proto>.root,....
  // where "<proto>" is a key the gOptFluxShortNames map

  map<string,string>::const_iterator mitr     = gOptFluxShortNames.begin();
  map<string,string>::const_iterator mitr_end = gOptFluxShortNames.end();
  for ( ; mitr != mitr_end; ++mitr ) {
    string proto  = mitr->first + string(":");
    string gproto = string("g") + proto;
    string protor = proto + ".root,";
    string full  = mitr->second;
    if ( fopt.find(proto) == 0 ) {
      fopt.erase(0,proto.size());
      gOptFluxDriver = full;
      break;
    } else if ( fopt.find(gproto) == 0 ) {
      fopt.erase(0,gproto.size());
      gOptFluxDriver = full;
      break;
    } else if ( fopt.find(protor) != std::string::npos ) {
      gOptFluxDriver = full;
      break;
    }
  }
  // tested all cases where user might have specified explicitly
  // or been part of an extended file extension
  // this is where it gets messy
  if ( gOptFluxDriver == "" ) {

    // not specified ? guess from file name itself
    if ( fopt.find("gsimple") != std::string::npos ) {
      // put dk2nu after gsimple in case simple files are derived from dk2nu
      // then both are in name we should choose gsimple
      gOptFluxDriver = "genie::flux::GSimpleNtpFlux";
    } else if ( fopt.find("dk2nu") != std::string::npos ) {
      gOptFluxDriver = "genie::flux::GDk2NuFlux";
    } else {
      // does it look like the histogram format
      const char* hstrings[] = { ",12[", ",+12[", ",-12[",
                                 ",14[", ",+14[", ",-14[",
                                 ",16[", ",+16[", ",-16[" };
      size_t nh = sizeof(hstrings)/sizeof(const char*);
      for (size_t ih=0; ih<nh; ++ih) {
        if ( fopt.find(hstrings[ih]) != std::string::npos ) {
          // hey!
          gOptFluxDriver = "genie::flux::GCylindTH1Flux";
          break;
        }
      } // loop over possible histogram specifiers
    }

    // fall through default ... hope it works
    if ( gOptFluxDriver == "" ) {
      gOptFluxDriver = "genie::flux::GNuMIFlux";
    }
  }

  gOptUsingHistFlux = ( gOptFluxDriver == "genie::flux::GCylindTH1Flux" );
  if ( gOptUsingHistFlux ) ParseFluxHst(fopt);
  else                     ParseFluxFileConfig(fopt);
}
//____________________________________________________________________________
void ParseFluxHst(string flux)
{
  // 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_fnal", 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_fnal", 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_fnal", 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_fnal", pFATAL)
        << "Can not find histogram: " << histo
        << " in flux file: " << gOptFluxFile;
      PrintSyntax();
      exit(1);
    }

    // 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("spectrum","neutrino_flux");
    spectrum->SetDirectory(0);
    //// and remove bins outside the emin,emax range (not for gFNALExptEvGen)
    //for(int ibin = 1; ibin <= hst->GetNbinsX(); ibin++) {
    //  if(hst->GetBinLowEdge(ibin) + hst->GetBinWidth(ibin) > emax ||
    //     hst->GetBinLowEdge(ibin) < emin) {
    //    spectrum->SetBinContent(ibin, 0);
    //  }
    //}

    // get rid of original
    delete ihst;

    // convert neutrino name -> pdg code
    int pdg = atoi(nutype.c_str());
    if(!pdg::IsNeutrino(pdg) && !pdg::IsAntiNeutrino(pdg)) {
      LOG("gevgen_fnal", pFATAL)
        << "Unknown neutrino type: " << nutype;
      PrintSyntax();
      exit(1);
    }
    // store flux neutrino code / energy spectrum
    LOG("gevgen_fnal", 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_fnal", pFATAL)
      << "You have not specified any flux histogram!";
    PrintSyntax();
    exit(1);
  }

  flux_file.Close();
}

//____________________________________________________________________________
void ParseFluxFileConfig(string flux)
{
  // Using gnumi/gsimple/dk2nu beam flux ntuples
  // Extract beam flux (root) file name & detector location
  //
  vector<string> fluxv = utils::str::Split(flux,",");
  if(fluxv.size()<2) {
    LOG("gevgen_fnal", pFATAL)
      << "You need to specify both a flux ntuple ROOT file "
      << " _AND_ a detector location";
    PrintSyntax();
    exit(1);
  }
  gOptFluxFile         = fluxv[0];
  gOptDetectorLocation = fluxv[1];

  for ( size_t j = 2; j < fluxv.size(); ++j ) {
    int ipdg = atoi(fluxv[j].c_str());
    gOptFluxPdg.push_back(ipdg);
  }

}

//____________________________________________________________________________