gFNALExptEvGen.cxx File Reference

#include <cassert>
#include <cstdlib>
#include <csignal>
#include <string>
#include <sstream>
#include <vector>
#include <map>
#include <algorithm>
#include <fstream>
#include <TSystem.h>
#include <TError.h>
#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"

Go to the source code of this file.

Functions
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)
static void	gsSIGTERMhandler (int)
int	main (int argc, char **argv)

Variables
string	kDefOptGeomLUnits = "mm"
string	kDefOptGeomDUnits = "g_cm3"
NtpMCFormat_t	kDefOptNtpFormat = kNFGHEP
string	kDefOptEvFilePrefix = "gntp"
Long_t	gOptRunNu
bool	gOptUsingRootGeom = false
bool	gOptUsingHistFlux = false
string	gOptFluxDriver = ""
map< string, string >	gOptFluxShortNames
PDGCodeList	gOptFluxPdg
map< int, double >	gOptTgtMix
map< int, TH1D * >	gOptFluxHst
string	gOptRootGeom
string	gOptRootGeomTopVol = ""
string	gOptRootGeomMasterVol = ""
double	gOptGeomLUnits = 0
double	gOptGeomDUnits = 0
string	gOptExtMaxPlXml = ""
bool	gOptWriteMaxPlXml = false
string	gOptFluxFile
string	gOptDetectorLocation
int	gOptNev
double	gOptPOT
string	gOptFidCut
int	gOptNScan = 0
double	gOptZmin = -2.0e30
string	gOptEvFilePrefix
int	gOptDebug = 0
long int	gOptRanSeed
string	gOptInpXSecFile
bool	gSigTERM = false

Function Documentation

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

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}

Definition at line 1300 of file gFNALExptEvGen.cxx.

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

Definition at line 1464 of file gFNALExptEvGen.cxx.

{

  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

)

Definition at line 1563 of file gFNALExptEvGen.cxx.

{
  // 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 GetCommandLineArgs	(	int	argc,
		char **	argv
	)

Definition at line 906 of file gFNALExptEvGen.cxx.

{
  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();
}

static void gsSIGTERMhandler ( int  )

static

Definition at line 397 of file gFNALExptEvGen.cxx.

{
  gSigTERM = true;
  std::cerr << "Caught SIGTERM" << std::endl;
}

void LoadExtraOptions

(

void

)

potentially load extra libraries that might extend the list of potential flux drivers, and how to map short names to classes ...

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

Definition at line 853 of file gFNALExptEvGen.cxx.

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

}

int main	(	int	argc,
		char **	argv
	)

Definition at line 404 of file gFNALExptEvGen.cxx.

{
  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 ParseFluxFileConfig

(

string

fopt

)

Definition at line 1725 of file gFNALExptEvGen.cxx.

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

}

void ParseFluxHst

(

string

fopt

)

Definition at line 1630 of file gFNALExptEvGen.cxx.

{
  // 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 PrintSyntax

(

void

)

Definition at line 1274 of file gFNALExptEvGen.cxx.

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

Variable Documentation

int gOptDebug = 0

Definition at line 391 of file gFNALExptEvGen.cxx.

string gOptDetectorLocation

Definition at line 384 of file gFNALExptEvGen.cxx.

string gOptEvFilePrefix

Definition at line 390 of file gFNALExptEvGen.cxx.

string gOptExtMaxPlXml = ""

Definition at line 380 of file gFNALExptEvGen.cxx.

string gOptFidCut

Definition at line 387 of file gFNALExptEvGen.cxx.

string gOptFluxDriver = ""

Definition at line 370 of file gFNALExptEvGen.cxx.

string gOptFluxFile

Definition at line 383 of file gFNALExptEvGen.cxx.

map<int,TH1D*> gOptFluxHst

Definition at line 374 of file gFNALExptEvGen.cxx.

PDGCodeList gOptFluxPdg

Definition at line 372 of file gFNALExptEvGen.cxx.

map<string,string> gOptFluxShortNames

Definition at line 371 of file gFNALExptEvGen.cxx.

double gOptGeomDUnits = 0

Definition at line 379 of file gFNALExptEvGen.cxx.

double gOptGeomLUnits = 0

Definition at line 378 of file gFNALExptEvGen.cxx.

string gOptInpXSecFile

Definition at line 393 of file gFNALExptEvGen.cxx.

int gOptNev

Definition at line 385 of file gFNALExptEvGen.cxx.

int gOptNScan = 0

Definition at line 388 of file gFNALExptEvGen.cxx.

double gOptPOT

Definition at line 386 of file gFNALExptEvGen.cxx.

long int gOptRanSeed

Definition at line 392 of file gFNALExptEvGen.cxx.

string gOptRootGeom

Definition at line 375 of file gFNALExptEvGen.cxx.

string gOptRootGeomMasterVol = ""

Definition at line 377 of file gFNALExptEvGen.cxx.

string gOptRootGeomTopVol = ""

Definition at line 376 of file gFNALExptEvGen.cxx.

Long_t gOptRunNu

Definition at line 367 of file gFNALExptEvGen.cxx.

map<int,double> gOptTgtMix

Definition at line 373 of file gFNALExptEvGen.cxx.

bool gOptUsingHistFlux = false

Definition at line 369 of file gFNALExptEvGen.cxx.

bool gOptUsingRootGeom = false

Definition at line 368 of file gFNALExptEvGen.cxx.

bool gOptWriteMaxPlXml = false

Definition at line 381 of file gFNALExptEvGen.cxx.

double gOptZmin = -2.0e30

Definition at line 389 of file gFNALExptEvGen.cxx.

bool gSigTERM = false

Definition at line 395 of file gFNALExptEvGen.cxx.

string kDefOptEvFilePrefix = "gntp"

Definition at line 363 of file gFNALExptEvGen.cxx.

string kDefOptGeomDUnits = "g_cm3"

Definition at line 361 of file gFNALExptEvGen.cxx.

string kDefOptGeomLUnits = "mm"

Definition at line 360 of file gFNALExptEvGen.cxx.

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 362 of file gFNALExptEvGen.cxx.