gEvGenLArDM.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/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"

Include dependency graph for gEvGenLArDM.cxx:

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"
double	gOptZpCoupling
double	gOptDMMass
double	gOptMedRatio
Long_t	gOptRunNu
bool	gOptUsingRootGeom = false
map< int, double >	gOptTgtMix
string	gOptRootGeom
string	gOptRootGeomTopVol = ""
string	gOptRootGeomMasterVol = ""
double	gOptGeomLUnits = 0
double	gOptGeomDUnits = 0
string	gOptExtMaxPlXml = ""
bool	gOptWriteMaxPlXml = false
string	gOptFluxFile
int	gOptNev
string	gOptFidCut
int	gOptNScan = 0
double	gOptZmin = -2.0e30
string	gOptEvFilePrefix
int	gOptDebug = 0
long int	gOptRanSeed
string	gOptInpXSecFile
bool	gSigTERM = false

Function Documentation

CreateFidSelection()

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

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}

Definition at line 1117 of file gEvGenLArDM.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:
  ///      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);
 
}

References genie::geometry::ROOTGeomAnalyzer::AdoptGeomVolSelector(), genie::geometry::GeomVolSelectorFiducial::ConvertShapeMaster2Top(), LOG, genie::geometry::GeomVolSelectorFiducial::MakeBox(), genie::geometry::GeomVolSelectorFiducial::MakeSphere(), genie::geometry::GeomVolSelectorFiducial::MakeZCylinder(), genie::geometry::GeomVolSelectorFiducial::MakeZPolygon(), pFATAL, pNOTICE, pWARN, genie::geometry::GeomVolSelectorI::SetRemoveEntries(), genie::geometry::GeomVolSelectorFiducial::SetReverseFiducial(), and genie::utils::str::Split().

Referenced by main().

CreateRockBoxSelection()

void CreateRockBoxSelection	(	string	fidcut,
		GeomAnalyzerI *	geom_driver )

Definition at line 1250 of file gEvGenLArDM.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_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);
 
}

References genie::geometry::ROOTGeomAnalyzer::AdoptGeomVolSelector(), gOptRootGeomMasterVol, gOptRootGeomTopVol, LOG, genie::geometry::GeomVolSelectorFiducial::MakeBox(), genie::geometry::GeomVolSelectorFiducial::MakeSphere(), pFATAL, pWARN, genie::geometry::GeomVolSelectorRockBox::SetDeDx(), genie::geometry::GeomVolSelectorRockBox::SetExpandFromInclusion(), genie::geometry::GeomVolSelectorRockBox::SetMinimumWall(), genie::geometry::GeomVolSelectorI::SetRemoveEntries(), genie::geometry::GeomVolSelectorRockBox::SetRockBoxMinimal(), genie::geometry::ROOTGeomAnalyzer::SetTopVolName(), and genie::utils::str::Split().

Referenced by main().

DetermineFluxDriver()

void DetermineFluxDriver

(

string

fopt

)

Definition at line 1574 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);
}

References gOptFluxDriver, gOptFluxShortNames, gOptUsingHistFlux, ParseFluxFileConfig(), and ParseFluxHst().

Referenced by GetCommandLineArgs().

GetCommandLineArgs()

void GetCommandLineArgs	(	int	argc,
		char **	argv )

Definition at line 774 of file gEvGenLArDM.cxx.

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

References genie::CmdLnArgParser::ArgAsDouble(), genie::CmdLnArgParser::ArgAsInt(), genie::CmdLnArgParser::ArgAsLong(), genie::CmdLnArgParser::ArgAsString(), genie::RunOpt::EnableBareXSecPreCalc(), genie::PDGLibrary::Find(), geom, gOptDebug, gOptDMMass, gOptEvFilePrefix, gOptExtMaxPlXml, gOptFidCut, gOptFluxFile, gOptGeomDUnits, gOptGeomLUnits, gOptInpXSecFile, gOptMedRatio, gOptNev, gOptNScan, gOptRanSeed, gOptRootGeom, gOptRootGeomTopVol, gOptRunNu, gOptTgtMix, gOptUsingRootGeom, gOptWriteMaxPlXml, gOptZmin, gOptZpCoupling, genie::PDGLibrary::Instance(), genie::RunOpt::Instance(), kDefOptEvFilePrefix, kDefOptGeomDUnits, kDefOptGeomLUnits, LOG, lunits, genie::CmdLnArgParser::OptionExists(), pDEBUG, pFATAL, pINFO, pNOTICE, genie::utils::print::PrintFramedMesg(), PrintSyntax(), genie::RunOpt::ReadFromCommandLine(), genie::utils::str::Split(), and genie::utils::units::UnitFromString().

Referenced by main().

gsSIGTERMhandler()

void gsSIGTERMhandler ( int )

static

Definition at line 376 of file gEvGenLArDM.cxx.

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

References gSigTERM.

Referenced by main().

LoadExtraOptions()

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 730 of file gEvGenLArDM.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");
 
  ///******* 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;
 
}

References LOG, and pNOTICE.

Referenced by main().

main()

int main	(	int	argc,
		char **	argv )

Definition at line 383 of file gEvGenLArDM.cxx.

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

References genie::PDGLibrary::AddDarkMatter(), genie::NtpWriter::AddEventRecord(), genie::flux::GFluxDriverFactory::AvailableFluxDrivers(), genie::RunOpt::BuildTune(), genie::utils::app_init::CacheFile(), genie::AlgConfigPool::CommonList(), genie::GMCJDriver::Configure(), CreateFidSelection(), CreateRockBoxSelection(), genie::NtpWriter::CustomizeFilenamePrefix(), genie::GFluxI::End(), genie::NtpWriter::EventTree(), genie::GMCJDriver::ForceSingleProbScale(), genie::GMCJDriver::GenerateEvent(), genie::flux::GFluxFileConfigI::GetBranchInfo(), GetCommandLineArgs(), genie::flux::GFluxDriverFactory::GetFluxDriver(), genie::geometry::ROOTGeomAnalyzer::GetGeometry(), genie::geometry::ROOTGeomAnalyzer::GetMaxPathLengths(), genie::flux::GFluxFileConfigI::GetMetaDataTree(), gOptDebug, gOptDMMass, gOptEvFilePrefix, gOptExtMaxPlXml, gOptFidCut, gOptFluxFile, gOptGeomDUnits, gOptGeomLUnits, gOptInpXSecFile, gOptMedRatio, gOptNev, gOptNScan, gOptRanSeed, gOptRootGeom, gOptRootGeomMasterVol, gOptRootGeomTopVol, gOptRunNu, gOptTgtMix, gOptUsingRootGeom, gOptWriteMaxPlXml, gOptZmin, gOptZpCoupling, gSigTERM, gsSIGTERMhandler(), genie::NtpWriter::Initialize(), genie::AlgConfigPool::Instance(), genie::flux::GFluxDriverFactory::Instance(), genie::PDGLibrary::Instance(), genie::RunOpt::Instance(), kDefOptNtpFormat, genie::kPdgDarkMatter, genie::flux::GFluxFileConfigI::LoadBeamSimData(), LoadExtraOptions(), genie::Registry::Lock(), LOG, genie::utils::app_init::MesgThresholds(), pERROR, pFATAL, pINFO, pNOTICE, genie::PDGCodeList::push_back(), pWARN, genie::utils::app_init::RandGen(), genie::utils::geometry::RecursiveExhaust(), genie::NtpWriter::Save(), genie::PathLengthList::SaveAsXml(), genie::Registry::Set(), genie::GMCJDriver::SetEventGeneratorList(), genie::flux::GFluxFileConfigI::SetFluxParticles(), genie::flux::GFluxFileConfigI::SetNumOfCycles(), genie::GHepRecord::SetPrintLevel(), genie::GMCJMonitor::SetRefreshRate(), genie::geometry::ROOTGeomAnalyzer::SetScannerFlux(), genie::geometry::ROOTGeomAnalyzer::SetScannerNParticles(), genie::geometry::ROOTGeomAnalyzer::SetScannerNPoints(), genie::geometry::ROOTGeomAnalyzer::SetScannerNRays(), genie::flux::GFluxFileConfigI::SetUpstreamZ(), genie::Registry::UnLock(), genie::GMCJMonitor::Update(), genie::GMCJDriver::UseFluxDriver(), genie::GMCJDriver::UseGeomAnalyzer(), genie::GMCJDriver::UseMaxPathLengths(), genie::GMCJDriver::UseSplines(), and genie::utils::app_init::XSecTable().

ParseFluxFileConfig()

void ParseFluxFileConfig

(

string

fopt

)

Definition at line 1813 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);
  }
 
}

References gOptDetectorLocation, gOptFluxFile, gOptFluxPdg, LOG, pFATAL, PrintSyntax(), and genie::utils::str::Split().

Referenced by DetermineFluxDriver().

ParseFluxHst()

void ParseFluxHst

(

string

fopt

)

Definition at line 1641 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],...')
  // for variable width histograms, default to multiply in the width
  //   histo1[species1]WIDTH = multiply in the width
  //   histo1[species1]NOWIDTH = don't multiply in the width
  // possibly with configuration ";r=1.1,dir=(0.1,0.2,0.3),spot=(-1,-2,-3)"
  // appended
  // See documentation on top section of this file.
 
  vector<string> fluxtop = utils::str::Split(flux,";");
  string beamsettings;
  string histosettings;
  if ( fluxtop.size() == 1 ) {
    histosettings = fluxtop[0];
    LOG("gevgen_fnal", pNOTICE)
      << "ParseFluxHst: no settings for radius, direction, spot position"
      << " using defaults, r=" << gOptBeamRadius;
  } else if ( fluxtop.size() > 2 ) {
    LOG("gevgen_fnal", pFATAL)
      << "ParseFluxHst: too many ';' separated fields";
    PrintSyntax();
    exit(1);
  } else {
    // decide which is which depending on which has "=" in it
    string::size_type eqpos0 = fluxtop[0].find("=");
    string::size_type eqpos1 = fluxtop[1].find("=");
    if        (eqpos0 == string::npos && eqpos1 != string::npos ) {
      histosettings = fluxtop[0];
      beamsettings  = fluxtop[1];
    } else if (eqpos0 != string::npos && eqpos1 == string::npos ) {
      beamsettings  = fluxtop[0];
      histosettings = fluxtop[1];
    } else {
      LOG("gevgen_fnal", pFATAL)
        << "ParseFluxHst: too many / not enough fields with '=' "
        << "\n" << fluxtop[0] << "\n" << fluxtop[1];
      PrintSyntax();
      exit(1);
    }
    // now parse the string we have
    double r=-1, dx=0, dy=0, dz=1, x=0, y=0, z=0;
    int nscan = sscanf(beamsettings.c_str(),
                       "r=%lf,dir=(%lf,%lf,%lf),spot=(%lf,%lf,%lf)",
                       &r,&dx,&dy,&dz,&x,&y,&z);
    cout << "nscan = " << nscan << endl;
    gOptBeamRadius = r;
    gOptBeamDir    = TVector3(dx,dy,dz);
    gOptBeamSpot   = TVector3(x,y,z);
 
  }
  LOG("gevgen_fnal", pNOTICE)
       << "ParseFluxHst: "
       << " using r=" << gOptBeamRadius
       << ",dir=("
       << gOptBeamDir.Px() << ","
       << gOptBeamDir.Py() << ","
       << gOptBeamDir.Pz() << ")"
       << ",spot=("
       << gOptBeamSpot.X() << ","
       << gOptBeamSpot.Y() << ","
       << gOptBeamSpot.Z() << ")";
 
 
  vector<string> fluxv = utils::str::Split(histosettings,",");
  if(fluxv.size()<2) {
    LOG("gevgen_fnal", pFATAL)
      << "You need to specify both a flux histogram ROOT file "
      << " _AND_ at least one histogram[pdg] mapping";
    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);
    string extra  = nutype_and_histo.substr(jend+1,string::npos);
    std::transform(extra.begin(),extra.end(),extra.begin(),::toupper);
    LOG("gevgen_fnal", pNOTICE) // pDEBUG
      << " =======> nutype " << nutype << " histo " << histo << " extra " << extra;
    // 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);
    // get rid of original
    delete ihst;
 
    bool force_binwidth = false;
#if ROOT_VERSION_CODE <= ROOT_VERSION(9,99,99)
    // GetRandom() sampling on variable bin width histograms does not
    // correctly account for bin widths for all versions of ROOT prior
    // to (currently forever).  At some point this might change and
    // the necessity of this code snippet will go away
    TAxis* xaxis = spectrum->GetXaxis();
    if (xaxis->IsVariableBinSize()) force_binwidth = true;
#endif
    if ( extra == "WIDTH"   ) force_binwidth = true;
    if ( extra == "NOWIDTH" ) force_binwidth = false;
    if ( force_binwidth ) {
      LOG("gevgen_fnal", pNOTICE)
        << "multiplying by bin width for histogram " << histo
        << " as " << spectrum->GetName() << " for nutype " << nutype
        << " from " << gOptFluxFile;
      for(int ibin = 1; ibin <= spectrum->GetNbinsX(); ++ibin) {
        double data = spectrum->GetBinContent(ibin);
        double width = spectrum->GetBinWidth(ibin);
        spectrum->SetBinContent(ibin,data*width);
      }
    }
 
    // 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();
}

References gOptBeamDir, gOptBeamRadius, gOptBeamSpot, gOptFluxFile, gOptFluxHst, genie::pdg::IsAntiNeutrino(), genie::pdg::IsNeutrino(), LOG, pDEBUG, pFATAL, pNOTICE, PrintSyntax(), and genie::utils::str::Split().

Referenced by DetermineFluxDriver().

PrintSyntax()

void PrintSyntax

(

void

)

Definition at line 1090 of file gEvGenLArDM.cxx.

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

References LOG, and pFATAL.

Referenced by GetCommandLineArgs().

Variable Documentation

gOptDebug

int gOptDebug = 0

Definition at line 370 of file gEvGenLArDM.cxx.

Referenced by GetCommandLineArgs(), and main().

gOptDMMass

double gOptDMMass

Definition at line 351 of file gEvGenLArDM.cxx.

gOptEvFilePrefix

string gOptEvFilePrefix

Definition at line 369 of file gEvGenLArDM.cxx.

gOptExtMaxPlXml

string gOptExtMaxPlXml = ""

Definition at line 361 of file gEvGenLArDM.cxx.

gOptFidCut

string gOptFidCut

Definition at line 366 of file gEvGenLArDM.cxx.

Referenced by GetCommandLineArgs(), and main().

gOptFluxFile

string gOptFluxFile

Definition at line 364 of file gEvGenLArDM.cxx.

Referenced by GetCommandLineArgs(), main(), ParseFluxFileConfig(), and ParseFluxHst().

gOptGeomDUnits

double gOptGeomDUnits = 0

Definition at line 360 of file gEvGenLArDM.cxx.

gOptGeomLUnits

double gOptGeomLUnits = 0

Definition at line 359 of file gEvGenLArDM.cxx.

gOptInpXSecFile

string gOptInpXSecFile

Definition at line 372 of file gEvGenLArDM.cxx.

gOptMedRatio

double gOptMedRatio

Definition at line 352 of file gEvGenLArDM.cxx.

gOptNev

int gOptNev

Definition at line 365 of file gEvGenLArDM.cxx.

gOptNScan

int gOptNScan = 0

Definition at line 367 of file gEvGenLArDM.cxx.

Referenced by GetCommandLineArgs(), and main().

gOptRanSeed

long int gOptRanSeed

Definition at line 371 of file gEvGenLArDM.cxx.

gOptRootGeom

string gOptRootGeom

Definition at line 356 of file gEvGenLArDM.cxx.

gOptRootGeomMasterVol

string gOptRootGeomMasterVol = ""

Definition at line 358 of file gEvGenLArDM.cxx.

Referenced by CreateRockBoxSelection(), and main().

gOptRootGeomTopVol

string gOptRootGeomTopVol = ""

Definition at line 357 of file gEvGenLArDM.cxx.

gOptRunNu

Long_t gOptRunNu

Definition at line 353 of file gEvGenLArDM.cxx.

gOptTgtMix

map<int,double> gOptTgtMix

Definition at line 355 of file gEvGenLArDM.cxx.

gOptUsingRootGeom

bool gOptUsingRootGeom = false

Definition at line 354 of file gEvGenLArDM.cxx.

gOptWriteMaxPlXml

bool gOptWriteMaxPlXml = false

Definition at line 362 of file gEvGenLArDM.cxx.

Referenced by GetCommandLineArgs(), and main().

gOptZmin

double gOptZmin = -2.0e30

Definition at line 368 of file gEvGenLArDM.cxx.

Referenced by GetCommandLineArgs(), and main().

gOptZpCoupling

double gOptZpCoupling

Definition at line 350 of file gEvGenLArDM.cxx.

gSigTERM

bool gSigTERM = false

Definition at line 374 of file gEvGenLArDM.cxx.

Referenced by gsSIGTERMhandler(), gsSIGTERMhandler(), and main().

kDefOptEvFilePrefix

string kDefOptEvFilePrefix = "gntp"

Definition at line 346 of file gEvGenLArDM.cxx.

kDefOptGeomDUnits

string kDefOptGeomDUnits = "g_cm3"

Definition at line 344 of file gEvGenLArDM.cxx.

kDefOptGeomLUnits

string kDefOptGeomLUnits = "mm"

Definition at line 343 of file gEvGenLArDM.cxx.

kDefOptNtpFormat

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 345 of file gEvGenLArDM.cxx.