gEvGenLArDM.cxx

Go to the documentation of this file.

//______________________________________________________________________________
/*!
 
\program gevgen_lardm
 
\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.
 
         *** Synopsis :
 
           gevgen_lardm [-h]
                       [-r run#]
                        -f flux
                        -g geometry
                        -M dm mass
                       [-c zp_coupling]
                       [-v med_ratio]
                       [-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]
                       [-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_lardm 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'
           -M
              Specifies the DM mass
           -c
              Specifies the Z' coupling constant
              Default: Value in UserPhysicsOptions.xml
           -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[_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[_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]
 
           -n
              Specifies how many events to generate.
 
           -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
              +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
              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_lardm
                       -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 gNuMI 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_lardm
                       -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 <c.andreopoulos \at cern.ch>
 University of Liverpool
 
         Robert Hatcher <rhatcher \at fnal.gov>
         Fermi National Accelerator Laboratory
 
\created August 20, 2008
 
\cpright Copyright (c) 2003-2025, 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/Algorithm/AlgConfigPool.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/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:
//
double          gOptZpCoupling;                // mediator coupling
double          gOptDMMass;                    // DM Mass
double          gOptMedRatio;                  // ratio of mediator to DM mass
Long_t          gOptRunNu;                     // run number
bool            gOptUsingRootGeom = false;     // using root geom or target mix?
map<int,double> gOptTgtMix;                    // target mix  (tgt pdg -> wght frac) / if not using detailed root geom
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 gnumi beam flux ntuple
int             gOptNev;                       // number of events to generate
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);
  PDGLibrary::Instance()->AddDarkMatter(gOptDMMass,gOptMedRatio);
  if (gOptZpCoupling > 0.) {
      Registry * r = AlgConfigPool::Instance()->CommonList("Param", "BoostedDarkMatter");
      r->UnLock();
      r->Set("ZpCoupling", gOptZpCoupling);
      r->Lock();
  }
 
  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
  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
 
  // *************************************************************************
  // * 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_numi", 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("genie::flux::GSimpleNtpFlux");
  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_lardm", pFATAL)
      << "Failed to get any flux driver from GFluxDriverFactory";
    exit(1);
  }
 
  genie::flux::GFluxFileConfigI* flux_file_config =
      dynamic_cast<genie::flux::GFluxFileConfigI*>(flux_driver);
  if ( ! flux_file_config ) {
      LOG("gevgen_lardm", pFATAL)
          << "Failed to get GFluxFileConfigI driver from GFluxDriverFactory";
      exit(1);
  }
 
  //
  // *** Using the detailed ntuple neutrino flux description
  //
  flux_file_config->LoadBeamSimData(gOptFluxFile, "");
  flux_file_config->SetUpstreamZ(gOptZmin);  // was "zmin" from bounding_box
  flux_file_config->SetNumOfCycles(0);
  PDGCodeList dm_pdg;
  dm_pdg.push_back(kPdgDarkMatter);
  flux_file_config->SetFluxParticles(dm_pdg);
 
  // these might come in handy ... avoid repeated dynamic_cast<> calls
  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 ) {
 
    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_lardm", 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_lardm", 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
        << "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_lardm", pERROR)
       << "Inconsistent info back from flux driver "
       << "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_lardm", 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_lardm", 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_lardm", 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;
 
     // 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_lardm", pWARN)
          << "** The flux driver read all the input flux entries: End()==true";
        break;
     }
     if ( ! event ) {
         LOG("gevgen_lardm", pERROR)
             << "Got a null generated neutino event! Retrying ...";
         continue;
     }
     LOG("gevgen_lardm", 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_lardm", pINFO)
    << "The GENIE MC job is done generating events - Cleaning up & exiting...";
 
  // *************************************************************************
  // * 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.
 
  LOG("gevgen_lardm", 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");
 
  ///******* 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_lardm",pNOTICE)
         << statWords << " (" << loadStatus << ") " << extraLibs[ilib];
  }
  // restore the ROOT error message level
  gErrorIgnoreLevel    = prevErrorLevel;
 
}
 
//____________________________________________________________________________
void GetCommandLineArgs(int argc, char ** argv)
{
  LOG("gevgen_lardm", 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_lardm", pDEBUG) << "Reading MC run number";
    gOptRunNu = parser.ArgAsLong('r');
  } else {
    LOG("gevgen_lardm", pDEBUG)
      << "Unspecified run number - Using default";
    gOptRunNu = 0;
  } //-r
 
  //
  // *** geometry
  //
 
  string geom = "";
  string lunits, dunits;
  if( parser.OptionExists('g') ) {
    LOG("gevgen_lardm", 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_lardm", pFATAL)
        << "No geometry option specified - Exiting";
      PrintSyntax();
      exit(1);
  } //-g
 
  // dark matter mass
  if( parser.OptionExists('M') ) {
    LOG("gevgen_lardm", pINFO) << "Reading dark matter mass";
    gOptDMMass = parser.ArgAsDouble('M');
  } else {
    LOG("gevgen_lardm", pFATAL) << "Unspecified dark matter mass - Exiting";
    PrintSyntax();
    exit(1);
  } // -M
 
  // mediator coupling
  if( parser.OptionExists('c') ) {
    LOG("gevgen_lardm", pINFO) << "Reading mediator coupling";
    gOptZpCoupling = parser.ArgAsDouble('c');
  } else {
    LOG("gevgen_lardm", pINFO) << "Unspecified mediator coupling - Using value from config file";
    gOptZpCoupling = -1.;
  } // -c
 
  // mediator mass ratio
  if( parser.OptionExists('v') ) {
    LOG("gevgen_lardm", pINFO) << "Reading mediator mass ratio";
    gOptMedRatio = parser.ArgAsDouble('v');
  } else {
    LOG("gevgen_lardm", pINFO) << "Unspecified mediator mass ratio - Using default";
    gOptMedRatio = 0.5;
  } // -v
 
  if(gOptUsingRootGeom) {
     // using a ROOT geometry - get requested geometry units
 
     // legth units:
     if( parser.OptionExists('L') ) {
        LOG("gevgen_lardm", pDEBUG)
           << "Checking for input geometry length units";
        lunits = parser.ArgAsString('L');
     } else {
        LOG("gevgen_lardm", pDEBUG) << "Using default geometry length units";
        lunits = kDefOptGeomLUnits;
     } // -L
     // density units:
     if( parser.OptionExists('D') ) {
        LOG("gevgen_lardm", pDEBUG)
           << "Checking for input geometry density units";
        dunits = parser.ArgAsString('D');
     } else {
        LOG("gevgen_lardm", 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_lardm", pDEBUG) << "Checking for input volume name";
        gOptRootGeomTopVol = parser.ArgAsString('t');
     } else {
        LOG("gevgen_lardm", 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_lardm", 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_lardm", pINFO)
            << "Will write maximum path lengths XML file: " << gOptExtMaxPlXml;
        }
     } else {
        LOG("gevgen_lardm", pDEBUG)
               << "Will compute the maximum path lengths at job init";
        gOptExtMaxPlXml = "";
     } // -m
 
     // fidcut:
     if( parser.OptionExists('F') ) {
       LOG("gevgen_lardm", pDEBUG) << "Using Fiducial cut?";
       gOptFidCut = parser.ArgAsString('F');
     } else {
       LOG("gevgen_lardm", pDEBUG) << "No fiducial volume cut";
       gOptFidCut = "";
     } //-F
 
     // how to scan the geometry (if relevant)
     if( parser.OptionExists('S') ) {
         LOG("gevgen_lardm", pDEBUG)  << "Reading requested geom scan count";
         gOptNScan = parser.ArgAsInt('S');
     } else {
         LOG("gevgen_lardm", pDEBUG) << "No geom scan count was requested";
         gOptNScan = 0;
     } //-S
 
     // z for flux rays to start
     if( parser.OptionExists('z') ) {
         LOG("gevgen_lardm", pDEBUG)  << "Reading requested zmin";
         gOptZmin = parser.ArgAsDouble('z');
     } else {
         LOG("gevgen_lardm", pDEBUG) << "No zmin was requested";
         gOptZmin = -2.0e30; // < -1.0e30 ==> leave it on flux window
     } //-z
 
     // debug flags
     if ( parser.OptionExists('d') ) {
         LOG("gevgen_lardm", pDEBUG) << "Reading debug flag value";
         gOptDebug = parser.ArgAsInt('d');
       } else {
         LOG("gevgen_lardm", pDEBUG) << "Unspecified debug flags - Using default";
         gOptDebug = 0;
       } //-d
 
  } // 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_lardm", 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_lardm", 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_lardm", pDEBUG) << "Getting input flux";
    gOptFluxFile = parser.ArgAsString('f');
  } else {
    LOG("gevgen_lardm", pFATAL) << "No flux info was specified - Exiting";
    PrintSyntax();
    exit(1);
  }
 
  // number of events to generate
  if( parser.OptionExists('n') ) {
    LOG("gevgen_lardm", pDEBUG)
        << "Reading limit on number of events to generate";
    gOptNev = parser.ArgAsInt('n');
  } else {
    LOG("gevgen_lardm", pDEBUG)
       << "Will keep on generating events till the flux driver stops";
    gOptNev = -1;
  } //-n
 
  // event file prefix
  if( parser.OptionExists('o') ) {
    LOG("gevgen_lardm", pDEBUG) << "Reading the event filename prefix";
    gOptEvFilePrefix = parser.ArgAsString('o');
  } else {
    LOG("gevgen_lardm", pDEBUG)
      << "Will set the default event filename prefix";
    gOptEvFilePrefix = kDefOptEvFilePrefix;
  } //-o
 
 
  // random number seed
  if( parser.OptionExists("seed") ) {
    LOG("gevgen_lardm", pINFO) << "Reading random number seed";
    gOptRanSeed = parser.ArgAsLong("seed");
  } else {
    LOG("gevgen_lardm", pINFO) << "Unspecified random number seed - Using default";
    gOptRanSeed = -1;
  }
 
  // input cross-section file
  if( parser.OptionExists("cross-sections") ) {
    LOG("gevgen_lardm", pINFO) << "Reading cross-section file";
    gOptInpXSecFile = parser.ArgAsString("cross-sections");
  } else {
    LOG("gevgen_lardm", pINFO) << "Unspecified cross-section file";
    gOptInpXSecFile = "";
  }
 
  //
  // >>> 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;
  fluxinfo << "file = "        << gOptFluxFile;
 
  ostringstream exposure;
  if(gOptNev > 0)
      exposure << "Number of events = " << gOptNev;
 
  LOG("gevgen_lardm", pNOTICE)
      << "\n\n"
      << utils::print::PrintFramedMesg("NuMI expt. event generation job configuration");
 
  LOG("gevgen_lardm", 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_lardm", pNOTICE) << *RunOpt::Instance();
}
//____________________________________________________________________________
void PrintSyntax(void)
{
  LOG("gevgen_lardm", pFATAL)
   << "\n **Syntax**"
   << "\n gevgen_lardm [-h] [-r run#]"
   << "\n            -f flux -g geometry -M dm_mass"
   << "\n            [-c zp_coupling] [-v med_ratio]"
   << "\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] "
   << "\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:
  ///      zcly:x0,y0,radius,zmin,zmax           - cylinder along z at (x0,y0) capped at z's
  ///      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) zcly:(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_lardm", pWARN)
      << "Can not create GeomVolSelectorFiduction,"
      << " geometry driver is not ROOTGeomAnalyzer";
    return;
  }
 
  LOG("gevgen_lardm", 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_lardm", pWARN)
      << "Can not create GeomVolSelectorFiduction,"
      << " no \":\" separating type from values.  nsplit=" << strtok.size();
    for ( unsigned int i=0; i < strtok.size(); ++i )
      LOG("gevgen_lardm",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_lardm", 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("box")    != string::npos ) {
    // box (xmin,ymin,zmin) (xmax,ymax,zmax)
    if ( nvals < 6 )
      LOG("gevgen_lardm", 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_lardm", pFATAL) << "MakeZPolygon needs 7 values, not " << nvals
                                << " fidcut=\"" << fidcut << "\"";
    int nfaces = (int)vals[0];
    if ( nfaces < 3 )
      LOG("gevgen_lardm", 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_lardm", pFATAL) << "MakeZSphere needs 4 values, not " << nvals
                                << " fidcut=\"" << fidcut << "\"";
    fidsel->MakeSphere(vals[0],vals[1],vals[2],vals[3]);
 
  } else {
    LOG("gevgen_lardm", pFATAL)
      << "Can not create GeomVolSelectorFiduction for shape \"" << stype << "\"";
  }
 
  if ( master  ) {
    fidsel->ConvertShapeMaster2Top(rgeom);
    LOG("gevgen_lardm", pNOTICE) << "Convert fiducial volume from master to topvol coords";
  }
  if ( reverse ) {
    fidsel->SetReverseFiducial(true);
    LOG("gevgen_lardm", 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_numi", pWARN)
      << "Can not create GeomVolSelectorRockBox,"
      << " geometry driver is not ROOTGeomAnalyzer";
    return;
  }
 
  LOG("gevgen_numi", 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_numi", pWARN)
      << "Can not create GeomVolSelectorRockBox,"
      << " no \":\" separating type from values.  nsplit=" << strtok.size();
    for ( unsigned int i=0; i < strtok.size(); ++i )
      LOG("gevgen_numi", 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_numi", 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_numi", 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);
 
}
//____________________________________________________________________________