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

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

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.

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().

 {
   ///
   /// User defined fiducial volume cut
   ///      [0][M]<SHAPE>:val1,val2,...
   ///   "0" means reverse the cut (i.e. exclude the volume)
   ///   "M" means the coordinates are given in the ROOT geometry
   ///       "master" system and need to be transformed to "top vol" system
   ///   <SHAPE> can be any of "zcyl" "box" "zpoly" "sphere"
   ///       [each takes different # of args]
   ///   This must be followed by a ":" and a list of values separated by punctuation
   ///       (allowed separators: commas , parentheses () braces {} or brackets [] )
   ///   Value mapping:
   ///      zcly:x0,y0,radius,zmin,zmax           - cylinder along z at (x0,y0) capped at z's
   ///      box:xmin,ymin,zmin,xmax,ymax,zmax     - box w/ upper & lower extremes
   ///      zpoly:nfaces,x0,y0,r_in,phi,zmin,zmax - nfaces sided polygon in x-y plane
   //       sphere:x0,y0,z0,radius                - sphere of fixed radius at (x0,y0,z0)
   ///   Examples:
   ///      1) 0mbox:0,0,0.25,1,1,8.75
   ///         exclude (i.e. reverse) a box in master coordinates w/ corners (0,0,0.25) (1,1,8.75)
   ///      2) mzpoly:6,(2,-1),1.75,0,{0.25,8.75}
   ///         six sided polygon in x-y plane, centered at x,y=(2,-1) w/ inscribed radius 1.75
   ///         no rotation (so first face is in y-z plane +r from center, i.e. hex sits on point)
   ///         limited to the z range of {0.25,8.75} in the master ROOT geom coordinates
   ///      3) zcly:(3,4),5.5,-2,10
   ///         a cylinder oriented parallel to the z axis in the "top vol" coordinates
   ///         at x,y=(3,4) with radius 5.5 and z range of {-2,10}
   ///
   geometry::ROOTGeomAnalyzer * rgeom =
     dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
   if ( ! rgeom ) {
     LOG("gevgen_lardm", pWARN)
       << "Can not create GeomVolSelectorFiduction,"
       << " geometry driver is not ROOTGeomAnalyzer";
     return;
   }
 
   LOG("gevgen_lardm", pNOTICE) << "-F " << fidcut;
 
   genie::geometry::GeomVolSelectorFiducial* fidsel =
     new genie::geometry::GeomVolSelectorFiducial();
 
   fidsel->SetRemoveEntries(true);  // drop segments that won't be considered
 
   // convert string to lowercase
   std::transform(fidcut.begin(),fidcut.end(),fidcut.begin(),::tolower);
 
   vector<string> strtok = genie::utils::str::Split(fidcut,":");
   if ( strtok.size() != 2 ) {
     LOG("gevgen_lardm", pWARN)
       << "Can not create GeomVolSelectorFiduction,"
       << " no \":\" separating type from values.  nsplit=" << strtok.size();
     for ( unsigned int i=0; i < strtok.size(); ++i )
       LOG("gevgen_lardm",pNOTICE)
         << "strtok[" << i << "] = \"" << strtok[i] << "\"";
     return;
   }
 
   // parse out optional "x" and "m"
   string stype = strtok[0];
   bool reverse = ( stype.find("0") != string::npos );
   bool master  = ( stype.find("m") != string::npos );  // action after values are set
 
   // parse out values
   vector<double> vals;
   vector<string> valstrs = genie::utils::str::Split(strtok[1]," ,;(){}[]");
   vector<string>::const_iterator iter = valstrs.begin();
   for ( ; iter != valstrs.end(); ++iter ) {
     const string& valstr1 = *iter;
     if ( valstr1 != "" ) vals.push_back(atof(valstr1.c_str()));
   }
   size_t nvals = vals.size();
 
   std::cout << "ivals = [";
   for (unsigned int i=0; i < nvals; ++i) {
     if (i>0) cout << ",";
     std::cout << vals[i];
   }
   std::cout << "]" << std::endl;
 
   // std::vector elements are required to be adjacent so we can treat address as ptr
 
   if        ( stype.find("zcyl")   != string::npos ) {
     // cylinder along z direction at (x0,y0) radius zmin zmax
     if ( nvals < 5 )
       LOG("gevgen_lardm", pFATAL) << "MakeZCylinder needs 5 values, not " << nvals
                                 << " fidcut=\"" << fidcut << "\"";
     fidsel->MakeZCylinder(vals[0],vals[1],vals[2],vals[3],vals[4]);
 
   } else if ( stype.find("box")    != string::npos ) {
     // box (xmin,ymin,zmin) (xmax,ymax,zmax)
     if ( nvals < 6 )
       LOG("gevgen_lardm", pFATAL) << "MakeBox needs 6 values, not " << nvals
                                 << " fidcut=\"" << fidcut << "\"";
     double xyzmin[3] = { vals[0], vals[1], vals[2] };
     double xyzmax[3] = { vals[3], vals[4], vals[5] };
     fidsel->MakeBox(xyzmin,xyzmax);
 
   } else if ( stype.find("zpoly")  != string::npos ) {
     // polygon along z direction nfaces at (x0,y0) radius phi zmin zmax
     if ( nvals < 7 )
       LOG("gevgen_lardm", pFATAL) << "MakeZPolygon needs 7 values, not " << nvals
                                 << " fidcut=\"" << fidcut << "\"";
     int nfaces = (int)vals[0];
     if ( nfaces < 3 )
       LOG("gevgen_lardm", pFATAL) << "MakeZPolygon needs nfaces>=3, not " << nfaces
                                 << " fidcut=\"" << fidcut << "\"";
     fidsel->MakeZPolygon(nfaces,vals[1],vals[2],vals[3],vals[4],vals[5],vals[6]);
 
   } else if ( stype.find("sphere") != string::npos ) {
     // sphere at (x0,y0,z0) radius
     if ( nvals < 4 )
       LOG("gevgen_lardm", pFATAL) << "MakeZSphere needs 4 values, not " << nvals
                                 << " fidcut=\"" << fidcut << "\"";
     fidsel->MakeSphere(vals[0],vals[1],vals[2],vals[3]);
 
   } else {
     LOG("gevgen_lardm", pFATAL)
       << "Can not create GeomVolSelectorFiduction for shape \"" << stype << "\"";
   }
 
   if ( master  ) {
     fidsel->ConvertShapeMaster2Top(rgeom);
     LOG("gevgen_lardm", pNOTICE) << "Convert fiducial volume from master to topvol coords";
   }
   if ( reverse ) {
     fidsel->SetReverseFiducial(true);
     LOG("gevgen_lardm", pNOTICE) << "Reverse sense of fiducial volume cut";
   }
   rgeom->AdoptGeomVolSelector(fidsel);
 
 }

void CreateRockBoxSelection	(	string	fidcut,
		GeomAnalyzerI *	geom_driver
	)

Definition at line 1250 of file gEvGenLArDM.cxx.

References e, 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(), and genie::utils::str::Split().

Referenced by main().

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

void DetermineFluxDriver ( string fopt )

Definition at line 1574 of file gFNALExptEvGen.cxx.

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

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

static void gsSIGTERMhandler ( int )

static

Definition at line 376 of file gEvGenLArDM.cxx.

References gSigTERM.

Referenced by main().

 {
   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 730 of file gEvGenLArDM.cxx.

References LOG, and pNOTICE.

Referenced by main().

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

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

void ParseFluxFileConfig ( string fopt )

Definition at line 1813 of file gFNALExptEvGen.cxx.

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

Referenced by DetermineFluxDriver().

 {
   // 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 1641 of file gFNALExptEvGen.cxx.

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().

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

void PrintSyntax ( void )

Variable Documentation

int gOptDebug = 0

Definition at line 370 of file gEvGenLArDM.cxx.

Referenced by main().

double gOptDMMass

Definition at line 351 of file gEvGenLArDM.cxx.

string gOptEvFilePrefix

Definition at line 369 of file gEvGenLArDM.cxx.

string gOptExtMaxPlXml = ""

Definition at line 361 of file gEvGenLArDM.cxx.

string gOptFidCut

Definition at line 366 of file gEvGenLArDM.cxx.

Referenced by main().

string gOptFluxFile

Definition at line 364 of file gEvGenLArDM.cxx.

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

double gOptGeomDUnits = 0

Definition at line 360 of file gEvGenLArDM.cxx.

double gOptGeomLUnits = 0

Definition at line 359 of file gEvGenLArDM.cxx.

string gOptInpXSecFile

Definition at line 372 of file gEvGenLArDM.cxx.

double gOptMedRatio

Definition at line 352 of file gEvGenLArDM.cxx.

int gOptNev

Definition at line 365 of file gEvGenLArDM.cxx.

int gOptNScan = 0

Definition at line 367 of file gEvGenLArDM.cxx.

Referenced by main().

long int gOptRanSeed

Definition at line 371 of file gEvGenLArDM.cxx.

string gOptRootGeom

Definition at line 356 of file gEvGenLArDM.cxx.

string gOptRootGeomMasterVol = ""

Definition at line 358 of file gEvGenLArDM.cxx.

Referenced by CreateRockBoxSelection(), and main().

string gOptRootGeomTopVol = ""

Definition at line 357 of file gEvGenLArDM.cxx.

Long_t gOptRunNu

Definition at line 353 of file gEvGenLArDM.cxx.

map<int,double> gOptTgtMix

Definition at line 355 of file gEvGenLArDM.cxx.

bool gOptUsingRootGeom = false

Definition at line 354 of file gEvGenLArDM.cxx.

bool gOptWriteMaxPlXml = false

Definition at line 362 of file gEvGenLArDM.cxx.

Referenced by main().

double gOptZmin = -2.0e30

Definition at line 368 of file gEvGenLArDM.cxx.

Referenced by main().

double gOptZpCoupling

Definition at line 350 of file gEvGenLArDM.cxx.

bool gSigTERM = false

Definition at line 374 of file gEvGenLArDM.cxx.

Referenced by gsSIGTERMhandler(), and main().

string kDefOptEvFilePrefix = "gntp"

Definition at line 346 of file gEvGenLArDM.cxx.

string kDefOptGeomDUnits = "g_cm3"

Definition at line 344 of file gEvGenLArDM.cxx.

string kDefOptGeomLUnits = "mm"

Definition at line 343 of file gEvGenLArDM.cxx.

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 345 of file gEvGenLArDM.cxx.

Functions

Variables

Function Documentation

Variable Documentation