#include <cassert>
#include <cstdlib>
#include <csignal>
#include <string>
#include <sstream>
#include <vector>
#include <map>
#include <algorithm>
#include <fstream>
#include <TSystem.h>
#include <TError.h>
#include <TTree.h>
#include <TFile.h>
#include <TH1D.h>
#include <TMath.h>
#include <TGeoVolume.h>
#include <TGeoShape.h>
#include "Framework/Conventions/Units.h"
#include "Framework/EventGen/EventRecord.h"
#include "Framework/EventGen/GFluxI.h"
#include "Framework/EventGen/GMCJDriver.h"
#include "Framework/EventGen/GMCJMonitor.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/Ntuple/NtpWriter.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGCodeList.h"
#include "Framework/Ntuple/NtpMCFormat.h"
#include "Framework/Utils/XSecSplineList.h"
#include "Framework/Utils/StringUtils.h"
#include "Framework/Utils/UnitUtils.h"
#include "Framework/Utils/AppInit.h"
#include "Framework/Utils/RunOpt.h"
#include "Framework/Utils/CmdLnArgParser.h"
#include "Framework/Utils/PrintUtils.h"
#include "Framework/Utils/SystemUtils.h"

Include dependency graph for gFNALExptEvGen.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"

Long_t	gOptRunNu

bool	gOptUsingRootGeom = false

bool	gOptUsingHistFlux = false

string	gOptFluxDriver = ""

map< string, string >	gOptFluxShortNames

PDGCodeList	gOptFluxPdg

map< int, double >	gOptTgtMix

map< int, TH1D * >	gOptFluxHst

TVector3	gOptBeamDir = TVector3(0,0,1)

TVector3	gOptBeamSpot = TVector3(0,0,0)

Double_t	gOptBeamRadius = -1

string	gOptBeamRtDependence = "x"

string	gOptRootGeom

string	gOptRootGeomTopVol = ""

string	gOptRootGeomMasterVol = ""

double	gOptGeomLUnits = 0

double	gOptGeomDUnits = 0

string	gOptExtMaxPlXml = ""

bool	gOptWriteMaxPlXml = false

string	gOptFluxFile

string	gOptDetectorLocation

int	gOptNev

double	gOptPOT

string	gOptFidCut

int	gOptNScan = 0

double	gOptZmin = -2.0e30

string	gOptEvFilePrefix

int	gOptDebug = 0

long int	gOptRanSeed

string	gOptInpXSecFile

bool	gSigTERM = false

Function Documentation

void CreateFidSelection	(	string	fidcut,
		GeomAnalyzerI *	geom_driver
	)

User defined fiducial volume cut [0][M]<SHAPE>:val1,val2,... "0" means reverse the cut (i.e. exclude the volume) "M" means the coordinates are given in the ROOT geometry "master" system and need to be transformed to "top vol" system <SHAPE> can be any of "zcyl" "box" "zpoly" "sphere" [each takes different # of args] This must be followed by a ":" and a list of values separated by punctuation (allowed separators: commas , parentheses () braces {} or brackets [] ) Value mapping: zcyl:x0,y0,radius,zmin,zmax - cylinder along z at (x0,y0) capped at z's xcyl:y0,z0,radius,xmin,ymax - cylinder along x ycyl:x0,z0,radius,ymin,ymax - cylinder along y gcyl:{x0,y0,z0}{dx,dy,dz},radius,{plane1}{plane2} – generic cylinder w/ arbitrary orientation and caps {planeX} = 4 values to define plane and orientation box:xmin,ymin,zmin,xmax,ymax,zmax - box w/ upper & lower extremes zpoly:nfaces,x0,y0,r_in,phi,zmin,zmax - nfaces sided polygon in x-y plane

Examples: 1) 0mbox:0,0,0.25,1,1,8.75 exclude (i.e. reverse) a box in master coordinates w/ corners (0,0,0.25) (1,1,8.75) 2) mzpoly:6,(2,-1),1.75,0,{0.25,8.75} six sided polygon in x-y plane, centered at x,y=(2,-1) w/ inscribed radius 1.75 no rotation (so first face is in y-z plane +r from center, i.e. hex sits on point) limited to the z range of {0.25,8.75} in the master ROOT geom coordinates 3) zcyl:(3,4),5.5,-2,10 a cylinder oriented parallel to the z axis in the "top vol" coordinates at x,y=(3,4) with radius 5.5 and z range of {-2,10}

void CreateRockBoxSelection	(	string	fidcut,
		GeomAnalyzerI *	geom_driver
	)

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 419 of file gFNALExptEvGen.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"

int main	(	int	argc,
		char **	argv
	)

Definition at line 426 of file gFNALExptEvGen.cxx.

 {
   LoadExtraOptions();
   GetCommandLineArgs(argc,argv);
 
   if ( ! RunOpt::Instance()->Tune() ) {
     LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
     exit(-1);
   }
   RunOpt::Instance()->BuildTune();
 
   // Initialization of random number generators, cross-section table,
   // messenger thresholds, cache file
   utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
   utils::app_init::CacheFile(RunOpt::Instance()->CacheFile());
   utils::app_init::RandGen(gOptRanSeed);
   utils::app_init::XSecTable(gOptInpXSecFile, false);
 
   // Set GHEP print level
   int print_level = RunOpt::Instance()->EventRecordPrintLevel();
   GHepRecord::SetPrintLevel(print_level);
 
   // *************************************************************************
   // * Create / configure the geometry driver
   // *************************************************************************
   GeomAnalyzerI * geom_driver = 0;
 
   if(gOptUsingRootGeom) {
     //
     // *** Using a realistic root-based detector geometry description
     //
 
     // creating & configuring a root geometry driver
     geometry::ROOTGeomAnalyzer * rgeom =
             new geometry::ROOTGeomAnalyzer(gOptRootGeom);
     TGeoVolume * topvol = rgeom->GetGeometry()->GetTopVolume();
     if ( ! topvol ) {
       LOG("gevgen_fnal", pFATAL) << "Null top ROOT geometry volume!";
       exit(1);
     }
     // retrieve the master volume name
     gOptRootGeomMasterVol = topvol->GetName();
 
     rgeom -> SetLengthUnits  (gOptGeomLUnits);
     rgeom -> SetDensityUnits (gOptGeomDUnits);
     rgeom -> SetTopVolName   (gOptRootGeomTopVol); // set user defined "topvol"
 
     // getting the bounding box dimensions along z so as to set the
     // appropriate upstream generation surface for the NuMI flux driver
 
     // RWH 2010-07-16:  do not try to automatically get zmin from geometry, rather
     // by default let the flux start from the window.  If the user wants to
     // override this then they need to explicitly set a "zmin".   Trying to use
     // the geometry is fraught with problems in local vs. global coordinates and
     // units where it can appear to work in some cases but it actually isn't really
     // universally correct.
     //was// TGeoShape * bounding_box = topvol->GetShape();
     //was// bounding_box->GetAxisRange(3, zmin, zmax);
     //was// zmin *= rgeom->LengthUnits();
     //was// zmax *= rgeom->LengthUnits();
 
     // switch on/off volumes as requested
     if ( (gOptRootGeomTopVol[0] == '+') || (gOptRootGeomTopVol[0] == '-') ) {
       bool exhaust = (*gOptRootGeomTopVol.c_str() == '+');
       utils::geometry::RecursiveExhaust(topvol, gOptRootGeomTopVol, exhaust);
     }
 
     // casting to the GENIE geometry driver interface
     geom_driver = dynamic_cast<GeomAnalyzerI *> (rgeom);
 
     // user specifid a fiducial volume cut ... parse that out
     if ( gOptFidCut.find("rock") != std::string::npos )
                                  CreateRockBoxSelection(gOptFidCut,rgeom);
     else if ( gOptFidCut != "" ) CreateFidSelection(gOptFidCut,rgeom);
 
   }
   else {
     //
     // *** Using a 'point' geometry with the specified target mix
     // *** ( = a list of targets with their corresponding weight fraction)
     //
 
     // creating & configuring a point geometry driver
     geometry::PointGeomAnalyzer * pgeom =
               new geometry::PointGeomAnalyzer(gOptTgtMix);
     // casting to the GENIE geometry driver interface
     geom_driver = dynamic_cast<GeomAnalyzerI *> (pgeom);
   }
 
   // *************************************************************************
   // * Create / configure the flux driver
   // *************************************************************************
   GFluxI * flux_driver =
     genie::flux::GFluxDriverFactory::Instance().GetFluxDriver(gOptFluxDriver);
   if ( ! flux_driver ) {
     // couldn't get the requested flux driver ?
     std::ostringstream s;
     s << "Known FluxDrivers:\n";
     const std::vector<std::string>& known =
       genie::flux::GFluxDriverFactory::Instance().AvailableFluxDrivers();
     std::vector<std::string>::const_iterator itr = known.begin();
     for ( ; itr != known.end(); ++itr ) s << "  " << (*itr) << "\n";
     LOG("gevgen_fnal", pFATAL)
       << "Failed to get any flux driver from GFluxDriverFactory\n"
       << "when using \"" << gOptFluxDriver << "\"\n" << s.str();
     exit(1);
   }
 
   if ( ! gOptUsingHistFlux ) {
     genie::flux::GFluxFileConfigI* flux_file_config =
       dynamic_cast<genie::flux::GFluxFileConfigI*>(flux_driver);
     if ( ! flux_file_config ) {
       LOG("gevgen_fnal", pFATAL)
         << "Failed to get GFluxFileConfigI driver from GFluxDriverFactory\n"
         << "when using \"" << gOptFluxDriver << "\"";
       exit(1);
     }
 
     //
     // *** Using the detailed ntuple neutrino flux description
     //
     flux_file_config->LoadBeamSimData(gOptFluxFile, gOptDetectorLocation);
     flux_file_config->SetUpstreamZ(gOptZmin);  // was "zmin" from bounding_box
     flux_file_config->SetNumOfCycles(0);
 
     if ( gOptFluxPdg.size() > 0 ) {
       // user specified list of neutrino PDGs
       flux_file_config->SetFluxParticles(gOptFluxPdg);
       std::ostringstream s;
       PDGCodeList::const_iterator itr = gOptFluxPdg.begin();
       for ( ; itr != gOptFluxPdg.end(); ++itr) s << (*itr) << " ";
       LOG("gevgen_fnal", pNOTICE)
         << "Limiting to nu PDGs: " << s.str();
     }
   }
   else {
     //
     // *** Using fluxes from histograms (for all specified neutrino species)
     //
     genie::flux::GCylindTH1Flux * hist_flux_driver =
       dynamic_cast<genie::flux::GCylindTH1Flux*>(flux_driver);
     if ( ! hist_flux_driver ) {
       LOG("gevgen_fnal", pFATAL)
         << "Failed to get GCylinderTH1Flux driver from GFluxDriverFactory\n"
         << "when using " << gOptFluxDriver;
       exit(1);
     }
 
     // creating & configuring a generic GCylindTH1Flux flux driver
     hist_flux_driver->SetNuDirection      (gOptBeamDir);
     hist_flux_driver->SetBeamSpot         (gOptBeamSpot);
     hist_flux_driver->SetTransverseRadius (gOptBeamRadius);
     //hist_flux_driver->SetRtDependence   (gOptBeamRtDependence); // "x");
 
     map<int,TH1D*>::iterator it = gOptFluxHst.begin();
     for( ; it != gOptFluxHst.end(); ++it) {
         int    pdg_code = it->first;
         TH1D * spectrum = it->second;
         hist_flux_driver->AddEnergySpectrum(pdg_code, spectrum);
         // once the histogram has been added to the GCylindTH1Flux driver
         // it is owned by the driver and it is up to the the driver
         // to clean up (i.e. delete it).
         // remove it from this map to avoid double deletion.
         it->second = 0;
     }
   }
 
   // these might come in handy ... avoid repeated dynamic_cast<> calls
   genie::flux::GFluxExposureI* fluxExposureI =
     dynamic_cast<genie::flux::GFluxExposureI*>(flux_driver);
   genie::flux::GFluxFileConfigI* fluxFileConfigI =
     dynamic_cast<genie::flux::GFluxFileConfigI*>(flux_driver);
 
 
   // *************************************************************************
   // * Handle chicken/egg problem: geom analyzer vs. flux.
   // * Need both at this point change geom scan defaults.
   // *************************************************************************
   if ( gOptUsingRootGeom && ! gOptUsingHistFlux ) {
 
     geometry::ROOTGeomAnalyzer * rgeom =
       dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
     if ( ! rgeom ) assert(0);
 
     rgeom -> SetDebugFlags(gOptDebug);
 
     // even if user doesn't specify gOptNScan configure to scan using flux
     if ( gOptNScan >= 0 ) {
       LOG("gevgen_fnal", pNOTICE)
         << "Using ROOTGeomAnalyzer: geom scan using flux: nparticles=" << gOptNScan;
       rgeom->SetScannerFlux(flux_driver);
       if ( gOptNScan > 0 ) rgeom->SetScannerNParticles(gOptNScan);
     } else {
       int nabs = TMath::Abs(gOptNScan);
       LOG("gevgen_fnal", pNOTICE)
         << "Using ROOTGeomAnalyzer: geom scan using box: npoints=nrays=" << nabs;
       rgeom->SetScannerNPoints(nabs);
       rgeom->SetScannerNRays(nabs);
     }
   }
 
   // *************************************************************************
   // * Create/configure the event generation driver
   // *************************************************************************
   GMCJDriver * mcj_driver = new GMCJDriver;
   mcj_driver->SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
   mcj_driver->UseFluxDriver(flux_driver);
   mcj_driver->UseGeomAnalyzer(geom_driver);
   if ( ( gOptExtMaxPlXml != "" ) && ! gOptWriteMaxPlXml ) {
     mcj_driver->UseMaxPathLengths(gOptExtMaxPlXml);
   }
   mcj_driver->Configure();
   mcj_driver->UseSplines();
   mcj_driver->ForceSingleProbScale();
 
   if ( ( gOptExtMaxPlXml != "" ) && gOptWriteMaxPlXml ) {
     geometry::ROOTGeomAnalyzer * rgeom =
       dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
     if ( rgeom ) {
       const genie::PathLengthList& maxpath = rgeom->GetMaxPathLengths();
       std::string maxplfile = gOptExtMaxPlXml;
       maxpath.SaveAsXml(maxplfile);
       // append extra info to file
       std::ofstream mpfile(maxplfile.c_str(), std::ios_base::app);
       mpfile
         << std::endl
         << "<!-- this file is only relevant for a setup compatible with:"
         << std::endl
         << "geom: " << gOptRootGeom << " top: \"" << gOptRootGeomTopVol << "\""
         << std::endl
         << "flux: " << gOptFluxFile
         << std::endl
         << "location: " << gOptDetectorLocation
         << std::endl
         << "fidcut: " << gOptFidCut
         << std::endl
         << "nscan: " << gOptNScan << " (0>= use flux, <0 use box |nscan| points/rays)"
         << std::endl
         << "zmin: " << gOptZmin << " (if |zmin| > 1e30, leave ray on flux window)"
         << std::endl
         << "-->"
         << std::endl;
       mpfile.close();
     }
   }
 
   // *************************************************************************
   // * Prepare for writing the output event tree & status file
   // *************************************************************************
 
   // Initialize an Ntuple Writer to save GHEP records into a TTree
   NtpWriter ntpw(kDefOptNtpFormat, gOptRunNu, gOptRanSeed);
   ntpw.CustomizeFilenamePrefix(gOptEvFilePrefix);
   ntpw.Initialize();
 
 
   std::vector<TBranch*>    extraBranches;
   std::vector<std::string> branchNames;
   std::vector<std::string> branchClassNames;
   std::vector<void**>      branchObjPointers;
 
   // Add custom branch(s) to the standard GENIE event tree so that
   // info on the flux neutrino parent particle can be passed-through
   if ( fluxFileConfigI ) {
     fluxFileConfigI->GetBranchInfo(branchNames,branchClassNames,
                                     branchObjPointers);
     size_t nn = branchNames.size();
     size_t nc = branchClassNames.size();
     size_t np = branchObjPointers.size();
     if ( nn != nc || nc != np ) {
      LOG("gevgen_fnal", pERROR)
        << "Inconsistent info back from \"" << gOptFluxDriver << "\" "
        << "for branch info:  " << nn << " " << nc << " " << np;
     } else {
       for (size_t ii = 0; ii < nn; ++ii) {
         const char* bname = branchNames[ii].c_str();
         const char* cname = branchClassNames[ii].c_str();
         void**&     optr  = branchObjPointers[ii];  // note critical '&' !
         if ( ! optr || ! *optr ) continue;  // no pointer supplied, skip it
         int split = 99;  // 1
         LOG("gevgen_fnal", pNOTICE)
           << "Adding extra branch \"" << bname << "\" of type \""
           << cname << "\" (" << optr << ") to output tree";
         TBranch* bptr = ntpw.EventTree()->Branch(bname,cname,optr,32000,split);
         extraBranches.push_back(bptr);
 
         if ( bptr ) {
           // don't delete this !!! we're sharing
           bptr->SetAutoDelete(false);
         } else {
           LOG("gevgen_fnal", pERROR)
           << "FAILED to add extra branch \"" << bname << "\" of type \""
           << cname << "\" to output tree";
         }
       } // loop over additions
     } // same # of entries
   } // of genie::flux::GFluxFileConfigI type
 
   // Create a MC job monitor for a periodically updated status file
   GMCJMonitor mcjmonitor(gOptRunNu);
   mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());
 
   // *************************************************************************
   // * Event generation loop
   // *************************************************************************
 
   // define handler to allow signal to end job gracefully
   signal(SIGTERM,gsSIGTERMhandler);
 
   int ievent = 0;
   while ( ! gSigTERM )
   {
      LOG("gevgen_fnal", pINFO)
           << " *** Generating event............ " << ievent;
 
      // In case the required statistics was expressed as 'number of events'
      // then quit if that number has been generated
      if ( ievent == gOptNev ) break;
 
      // In case the required statistics was expressed as 'number of POT'
      // then exit the event loop if the requested POT has been generated.
      if ( gOptPOT > 0 && fluxExposureI ) {
         double fpot = fluxExposureI->GetTotalExposure(); // current POTs used
         double psc  = mcj_driver->GlobProbScale();  // interaction prob. scale
         double pot  = fpot / psc;                   // POTs for generated sample
         if ( pot >= gOptPOT ) break;
      }
 
      // Generate a single event using neutrinos coming from the specified flux
      // and hitting the specified geometry or target mix
      EventRecord * event = mcj_driver->GenerateEvent();
 
      // Check whether a null event was returned due to the flux driver reaching
      // the end of the input flux ntuple - exit the event generation loop
      if ( ! event && flux_driver->End() ) {
         LOG("gevgen_fnal", pWARN)
           << "** The flux driver read all the input flux entries: End()==true";
         break;
      }
      if ( ! event ) {
          LOG("gevgen_fnal", pERROR)
              << "Got a null generated neutino event! Retrying ...";
          continue;
      }
      if ( print_level >= 0 ) {
        LOG("gevgen_fnal", pINFO)
          << "Generated event: " << *event;
      }
      // A valid event was generated: flux info (parent decay/prod
      // position/kinematics) for that simulated event should already
      // be connected to the right output tree branch
 
      // Add event at the output ntuple, refresh the mc job monitor & clean-up
      ntpw.AddEventRecord(ievent, event);
      mcjmonitor.Update(ievent,event);
      delete event;
      ievent++;
 
   } //1
 
   // Copy metadata tree, if available
   if ( fluxFileConfigI ) {
     TTree* t1 = fluxFileConfigI->GetMetaDataTree();
     if ( t1 ) {
       TTree* t2 = (TTree*)t1->CloneTree();
       t2->Write();
     }
   }
 
   LOG("gevgen_fnal", pINFO)
     << "The GENIE MC job is done generating events - Cleaning up & exiting...";
 
   // *************************************************************************
   // * Print job statistics &
   // * calculate normalization factor for the generated sample
   // *************************************************************************
   if ( ! gOptUsingHistFlux && gOptUsingRootGeom ) {
     // POT normalization will only be calculated if event generation was based
     // on beam simulation ntuples (not just histograms) & a detailed detector
     // geometry description.
     // Get nunber of flux neutrinos read-in by flux driver, number of flux
     // neutrinos actually thrown to the event generation driver and number
     // of neutrino interactions actually generated
     long int nflx     = 0;
     long int nflx_evg = mcj_driver-> NFluxNeutrinos();
     double   fpot     = 0;
     const char* exposureUnits = "(unknown units)";
     if ( fluxExposureI ) {
       fpot = fluxExposureI->GetTotalExposure(); // POTs used so far
       nflx = fluxExposureI->NFluxNeutrinos();
       //genie::flux::Exposure_t etype = fluxExposureI->GetExposureType();
       //exposureUnits = genie::flux::GFluxExposureI::AsString(etype);
       exposureUnits = fluxExposureI->GetExposureUnits();
     }
     if ( fluxFileConfigI ) {
       fluxFileConfigI->PrintConfig();
     }
     double psc   = mcj_driver->GlobProbScale();      // interaction prob. scale
     if ( psc <= 0.0 ) {
        LOG("gevgen_fnal", pFATAL) << "MCJobDriver GlobalProbScale was " << psc;
     }
     double pot   = fpot / psc;                       // POT for generated sample
     long int nev = ievent;
 
     LOG("gevgen_fnal", pNOTICE)
         << "\n >> Interaction probability scaling factor:  " << psc
         << "\n >> using: " << gOptFluxDriver
         << "\n >> N of flux v read-in by flux driver:      " << nflx
         << "\n >> N of flux v thrown to event gen driver:  " << nflx_evg
         << "\n >> N of generated v interactions:           " << nev
         << "\n ** Normalization for generated sample:      " << pot
         << " " << exposureUnits << " * detector";
 
     ntpw.EventTree()->SetWeight(pot); // store POT
 
   }
 
   // *************************************************************************
   // * Save & clean-up
   // *************************************************************************
 
   // Save the generated event tree & close the output file
   ntpw.Save();
 
   // Clean-up
   delete geom_driver;
   delete flux_driver;
   delete mcj_driver;
   // this list should only be histograms that have (for some reason)
   // not been handed over to the GCylindTH1Flux driver.
   map<int,TH1D*>::iterator it = gOptFluxHst.begin();
   for( ; it != gOptFluxHst.end(); ++it) {
     TH1D * spectrum = it->second;
     if(spectrum) delete spectrum;
   }
   gOptFluxHst.clear();
 
   LOG("gevgen_fnal", pNOTICE) << "Done!";
 
   return 0;
 }

void ParseFluxFileConfig ( string fopt )

Definition at line 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

TVector3 gOptBeamDir = TVector3(0,0,1)

Definition at line 392 of file gFNALExptEvGen.cxx.

Referenced by main(), and ParseFluxHst().

Double_t gOptBeamRadius = -1

Definition at line 394 of file gFNALExptEvGen.cxx.

Referenced by main(), and ParseFluxHst().

string gOptBeamRtDependence = "x"

Definition at line 395 of file gFNALExptEvGen.cxx.

TVector3 gOptBeamSpot = TVector3(0,0,0)

Definition at line 393 of file gFNALExptEvGen.cxx.

Referenced by main(), and ParseFluxHst().

int gOptDebug = 0

Definition at line 413 of file gFNALExptEvGen.cxx.

string gOptDetectorLocation

Definition at line 406 of file gFNALExptEvGen.cxx.

Referenced by main(), and ParseFluxFileConfig().

string gOptEvFilePrefix

Definition at line 412 of file gFNALExptEvGen.cxx.

string gOptExtMaxPlXml = ""

Definition at line 402 of file gFNALExptEvGen.cxx.

string gOptFidCut

Definition at line 409 of file gFNALExptEvGen.cxx.

string gOptFluxDriver = ""

Definition at line 386 of file gFNALExptEvGen.cxx.

Referenced by DetermineFluxDriver(), and main().

string gOptFluxFile

Definition at line 405 of file gFNALExptEvGen.cxx.

map<int,TH1D*> gOptFluxHst

Definition at line 391 of file gFNALExptEvGen.cxx.

Referenced by main(), and ParseFluxHst().

PDGCodeList gOptFluxPdg

Definition at line 388 of file gFNALExptEvGen.cxx.

Referenced by main(), and ParseFluxFileConfig().

map<string,string> gOptFluxShortNames

Definition at line 387 of file gFNALExptEvGen.cxx.

Referenced by DetermineFluxDriver().

double gOptGeomDUnits = 0

Definition at line 401 of file gFNALExptEvGen.cxx.

double gOptGeomLUnits = 0

Definition at line 400 of file gFNALExptEvGen.cxx.

string gOptInpXSecFile

Definition at line 415 of file gFNALExptEvGen.cxx.

int gOptNev

Definition at line 407 of file gFNALExptEvGen.cxx.

int gOptNScan = 0

Definition at line 410 of file gFNALExptEvGen.cxx.

double gOptPOT

Definition at line 408 of file gFNALExptEvGen.cxx.

Referenced by main().

long int gOptRanSeed

Definition at line 414 of file gFNALExptEvGen.cxx.

string gOptRootGeom

Definition at line 397 of file gFNALExptEvGen.cxx.

string gOptRootGeomMasterVol = ""

Definition at line 399 of file gFNALExptEvGen.cxx.

string gOptRootGeomTopVol = ""

Definition at line 398 of file gFNALExptEvGen.cxx.

Long_t gOptRunNu

Definition at line 383 of file gFNALExptEvGen.cxx.

map<int,double> gOptTgtMix

Definition at line 389 of file gFNALExptEvGen.cxx.

bool gOptUsingHistFlux = false

Definition at line 385 of file gFNALExptEvGen.cxx.

Referenced by DetermineFluxDriver(), and main().

bool gOptUsingRootGeom = false

Definition at line 384 of file gFNALExptEvGen.cxx.

bool gOptWriteMaxPlXml = false

Definition at line 403 of file gFNALExptEvGen.cxx.

double gOptZmin = -2.0e30

Definition at line 411 of file gFNALExptEvGen.cxx.

bool gSigTERM = false

Definition at line 417 of file gFNALExptEvGen.cxx.

string kDefOptEvFilePrefix = "gntp"

Definition at line 379 of file gFNALExptEvGen.cxx.

string kDefOptGeomDUnits = "g_cm3"

Definition at line 377 of file gFNALExptEvGen.cxx.

string kDefOptGeomLUnits = "mm"

Definition at line 376 of file gFNALExptEvGen.cxx.

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 378 of file gFNALExptEvGen.cxx.

Functions

Variables

Function Documentation

Variable Documentation