GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions | Variables
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

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

1118 {
1119  ///
1120  /// User defined fiducial volume cut
1121  /// [0][M]<SHAPE>:val1,val2,...
1122  /// "0" means reverse the cut (i.e. exclude the volume)
1123  /// "M" means the coordinates are given in the ROOT geometry
1124  /// "master" system and need to be transformed to "top vol" system
1125  /// <SHAPE> can be any of "zcyl" "box" "zpoly" "sphere"
1126  /// [each takes different # of args]
1127  /// This must be followed by a ":" and a list of values separated by punctuation
1128  /// (allowed separators: commas , parentheses () braces {} or brackets [] )
1129  /// Value mapping:
1130  /// zcly:x0,y0,radius,zmin,zmax - cylinder along z at (x0,y0) capped at z's
1131  /// box:xmin,ymin,zmin,xmax,ymax,zmax - box w/ upper & lower extremes
1132  /// zpoly:nfaces,x0,y0,r_in,phi,zmin,zmax - nfaces sided polygon in x-y plane
1133  // sphere:x0,y0,z0,radius - sphere of fixed radius at (x0,y0,z0)
1134  /// Examples:
1135  /// 1) 0mbox:0,0,0.25,1,1,8.75
1136  /// exclude (i.e. reverse) a box in master coordinates w/ corners (0,0,0.25) (1,1,8.75)
1137  /// 2) mzpoly:6,(2,-1),1.75,0,{0.25,8.75}
1138  /// six sided polygon in x-y plane, centered at x,y=(2,-1) w/ inscribed radius 1.75
1139  /// no rotation (so first face is in y-z plane +r from center, i.e. hex sits on point)
1140  /// limited to the z range of {0.25,8.75} in the master ROOT geom coordinates
1141  /// 3) zcly:(3,4),5.5,-2,10
1142  /// a cylinder oriented parallel to the z axis in the "top vol" coordinates
1143  /// at x,y=(3,4) with radius 5.5 and z range of {-2,10}
1144  ///
1145  geometry::ROOTGeomAnalyzer * rgeom =
1146  dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
1147  if ( ! rgeom ) {
1148  LOG("gevgen_lardm", pWARN)
1149  << "Can not create GeomVolSelectorFiduction,"
1150  << " geometry driver is not ROOTGeomAnalyzer";
1151  return;
1152  }
1153 
1154  LOG("gevgen_lardm", pNOTICE) << "-F " << fidcut;
1155 
1158 
1159  fidsel->SetRemoveEntries(true); // drop segments that won't be considered
1160 
1161  // convert string to lowercase
1162  std::transform(fidcut.begin(),fidcut.end(),fidcut.begin(),::tolower);
1163 
1164  vector<string> strtok = genie::utils::str::Split(fidcut,":");
1165  if ( strtok.size() != 2 ) {
1166  LOG("gevgen_lardm", pWARN)
1167  << "Can not create GeomVolSelectorFiduction,"
1168  << " no \":\" separating type from values. nsplit=" << strtok.size();
1169  for ( unsigned int i=0; i < strtok.size(); ++i )
1170  LOG("gevgen_lardm",pNOTICE)
1171  << "strtok[" << i << "] = \"" << strtok[i] << "\"";
1172  return;
1173  }
1174 
1175  // parse out optional "x" and "m"
1176  string stype = strtok[0];
1177  bool reverse = ( stype.find("0") != string::npos );
1178  bool master = ( stype.find("m") != string::npos ); // action after values are set
1179 
1180  // parse out values
1181  vector<double> vals;
1182  vector<string> valstrs = genie::utils::str::Split(strtok[1]," ,;(){}[]");
1183  vector<string>::const_iterator iter = valstrs.begin();
1184  for ( ; iter != valstrs.end(); ++iter ) {
1185  const string& valstr1 = *iter;
1186  if ( valstr1 != "" ) vals.push_back(atof(valstr1.c_str()));
1187  }
1188  size_t nvals = vals.size();
1189 
1190  std::cout << "ivals = [";
1191  for (unsigned int i=0; i < nvals; ++i) {
1192  if (i>0) cout << ",";
1193  std::cout << vals[i];
1194  }
1195  std::cout << "]" << std::endl;
1196 
1197  // std::vector elements are required to be adjacent so we can treat address as ptr
1198 
1199  if ( stype.find("zcyl") != string::npos ) {
1200  // cylinder along z direction at (x0,y0) radius zmin zmax
1201  if ( nvals < 5 )
1202  LOG("gevgen_lardm", pFATAL) << "MakeZCylinder needs 5 values, not " << nvals
1203  << " fidcut=\"" << fidcut << "\"";
1204  fidsel->MakeZCylinder(vals[0],vals[1],vals[2],vals[3],vals[4]);
1205 
1206  } else if ( stype.find("box") != string::npos ) {
1207  // box (xmin,ymin,zmin) (xmax,ymax,zmax)
1208  if ( nvals < 6 )
1209  LOG("gevgen_lardm", pFATAL) << "MakeBox needs 6 values, not " << nvals
1210  << " fidcut=\"" << fidcut << "\"";
1211  double xyzmin[3] = { vals[0], vals[1], vals[2] };
1212  double xyzmax[3] = { vals[3], vals[4], vals[5] };
1213  fidsel->MakeBox(xyzmin,xyzmax);
1214 
1215  } else if ( stype.find("zpoly") != string::npos ) {
1216  // polygon along z direction nfaces at (x0,y0) radius phi zmin zmax
1217  if ( nvals < 7 )
1218  LOG("gevgen_lardm", pFATAL) << "MakeZPolygon needs 7 values, not " << nvals
1219  << " fidcut=\"" << fidcut << "\"";
1220  int nfaces = (int)vals[0];
1221  if ( nfaces < 3 )
1222  LOG("gevgen_lardm", pFATAL) << "MakeZPolygon needs nfaces>=3, not " << nfaces
1223  << " fidcut=\"" << fidcut << "\"";
1224  fidsel->MakeZPolygon(nfaces,vals[1],vals[2],vals[3],vals[4],vals[5],vals[6]);
1225 
1226  } else if ( stype.find("sphere") != string::npos ) {
1227  // sphere at (x0,y0,z0) radius
1228  if ( nvals < 4 )
1229  LOG("gevgen_lardm", pFATAL) << "MakeZSphere needs 4 values, not " << nvals
1230  << " fidcut=\"" << fidcut << "\"";
1231  fidsel->MakeSphere(vals[0],vals[1],vals[2],vals[3]);
1232 
1233  } else {
1234  LOG("gevgen_lardm", pFATAL)
1235  << "Can not create GeomVolSelectorFiduction for shape \"" << stype << "\"";
1236  }
1237 
1238  if ( master ) {
1239  fidsel->ConvertShapeMaster2Top(rgeom);
1240  LOG("gevgen_lardm", pNOTICE) << "Convert fiducial volume from master to topvol coords";
1241  }
1242  if ( reverse ) {
1243  fidsel->SetReverseFiducial(true);
1244  LOG("gevgen_lardm", pNOTICE) << "Reverse sense of fiducial volume cut";
1245  }
1246  rgeom->AdoptGeomVolSelector(fidsel);
1247 
1248 }
void MakeSphere(Double_t x0, Double_t y0, Double_t z0, Double_t radius)
GENIE Interface for user-defined volume selector functors Trim path segments based on the intersectio...
virtual GeomVolSelectorI * AdoptGeomVolSelector(GeomVolSelectorI *selector)
configure processing to perform path segment trimming
void MakeZCylinder(Double_t x0, Double_t y0, Double_t radius, Double_t zmin, Double_t zmax)
#define pFATAL
Definition: Messenger.h:56
virtual void ConvertShapeMaster2Top(const ROOTGeomAnalyzer *rgeom)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void MakeZPolygon(Int_t n, Double_t x0, Double_t y0, Double_t inradius, Double_t phi0deg, Double_t zmin, Double_t zmax)
A ROOT/GEANT4 geometry driver.
#define pWARN
Definition: Messenger.h:60
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
#define pNOTICE
Definition: Messenger.h:61
void MakeBox(Double_t *xyzmin, Double_t *xyzmax)
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().

1251 {
1252 
1253  if( fidcut.find_first_not_of(" \t\n") != 0) // trim any leading whitespace
1254  fidcut.erase( 0, fidcut.find_first_not_of(" \t\n") );
1255 
1256  // convert string to lowercase
1257  std::transform(fidcut.begin(),fidcut.end(),fidcut.begin(),::tolower);
1258 
1260  dynamic_cast<genie::geometry::ROOTGeomAnalyzer *>(geom_driver);
1261  if ( ! rgeom ) {
1262  LOG("gevgen_numi", pWARN)
1263  << "Can not create GeomVolSelectorRockBox,"
1264  << " geometry driver is not ROOTGeomAnalyzer";
1265  return;
1266  }
1267 
1268  LOG("gevgen_numi", pWARN) << "fiducial (rock) cut: " << fidcut;
1269 
1270  // for now, only fiducial no "rock box"
1273 
1274  vector<string> strtok = genie::utils::str::Split(fidcut,":");
1275  if ( strtok.size() != 2 ) {
1276  LOG("gevgen_numi", pWARN)
1277  << "Can not create GeomVolSelectorRockBox,"
1278  << " no \":\" separating type from values. nsplit=" << strtok.size();
1279  for ( unsigned int i=0; i < strtok.size(); ++i )
1280  LOG("gevgen_numi", pWARN)
1281  << "strtok[" << i << "] = \"" << strtok[i] << "\"";
1282  return;
1283  }
1284 
1285  string stype = strtok[0];
1286 
1287  // parse out values
1288  vector<double> vals;
1289  vector<string> valstrs = genie::utils::str::Split(strtok[1]," ,;(){}[]\t\n\r");
1290  vector<string>::const_iterator iter = valstrs.begin();
1291  for ( ; iter != valstrs.end(); ++iter ) {
1292  const string& valstr1 = *iter;
1293  if ( valstr1 != "" ) {
1294  double aval = atof(valstr1.c_str());
1295  LOG("gevgen_numi", pWARN) << "rock value [" << vals.size() << "] "
1296  << aval;
1297  vals.push_back(aval);
1298  }
1299  }
1300  size_t nvals = vals.size();
1301 
1302  rocksel->SetRemoveEntries(true); // drop segments that won't be considered
1303 
1304  // assume coordinates are in the *master* (not "top volume") system
1305  // need to set fTopVolume to fWorldVolume
1306  //fTopVolume = fWorldVolume;
1307  //rgeom->SetTopVolName(fTopVolume.c_str());
1309  rgeom->SetTopVolName(gOptRootGeomMasterVol);
1310 
1311  if ( nvals < 6 ) {
1312  LOG("gevgen_numi", pFATAL) << "rockbox needs at "
1313  << "least 6 values, found "
1314  << nvals << "in \""
1315  << strtok[1] << "\"";
1316  exit(1);
1317 
1318  }
1319  double xyzmin[3] = { vals[0], vals[1], vals[2] };
1320  double xyzmax[3] = { vals[3], vals[4], vals[5] };
1321 
1322  bool rockonly = true;
1323  double wallmin = 800.; // geometry in cm, ( 8 meter buffer)
1324  double dedx = 2.5 * 1.7e-3; // GeV/cm, rho=2.5, 1.7e-3 ~ rock like loss
1325  double fudge = 1.05;
1326 
1327  if ( nvals >= 7 ) rockonly = vals[6];
1328  if ( nvals >= 8 ) wallmin = vals[7];
1329  if ( nvals >= 9 ) dedx = vals[8];
1330  if ( nvals >= 10 ) fudge = vals[9];
1331 
1332  rocksel->SetRockBoxMinimal(xyzmin,xyzmax);
1333  rocksel->SetMinimumWall(wallmin);
1334  rocksel->SetDeDx(dedx/fudge);
1335 
1336  if ( nvals >= 11 ) rocksel->SetExpandFromInclusion((int)vals[10]);
1337 
1338  // if not rock-only then make a tiny exclusion bubble
1339  // call to MakeBox shouldn't be necessary
1340  // should be done by SetRockBoxMinimal but for some GENIE versions isn't
1341  if ( ! rockonly ) rocksel->MakeSphere(0,0,0,1.0e-10);
1342  else rocksel->MakeBox(xyzmin,xyzmax);
1343 
1344  rgeom->AdoptGeomVolSelector(rocksel);
1345 
1346 }
void MakeSphere(Double_t x0, Double_t y0, Double_t z0, Double_t radius)
string gOptRootGeomMasterVol
#define pFATAL
Definition: Messenger.h:56
const double e
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
string gOptRootGeomTopVol
Definition: gAtmoEvGen.cxx:301
GENIE Interface for limiting vertex selection in the rock to a volume that depends (in part) on the n...
A ROOT/GEANT4 geometry driver.
#define pWARN
Definition: Messenger.h:60
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
void SetRockBoxMinimal(Double_t *xyzmin, Double_t *xyzmax)
void MakeBox(Double_t *xyzmin, Double_t *xyzmax)
void DetermineFluxDriver ( string  fopt)

Definition at line 1574 of file gFNALExptEvGen.cxx.

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

1575 {
1576  // based on the -f option string determine which flux driver to use
1577  // this may take some guessing
1578 
1579  // first look for strings that look like "<proto>:..."
1580  // or ....<proto>.root,....
1581  // where "<proto>" is a key the gOptFluxShortNames map
1582 
1583  map<string,string>::const_iterator mitr = gOptFluxShortNames.begin();
1584  map<string,string>::const_iterator mitr_end = gOptFluxShortNames.end();
1585  for ( ; mitr != mitr_end; ++mitr ) {
1586  string proto = mitr->first + string(":");
1587  string gproto = string("g") + proto;
1588  string protor = proto + ".root,";
1589  string full = mitr->second;
1590  if ( fopt.find(proto) == 0 ) {
1591  fopt.erase(0,proto.size());
1592  gOptFluxDriver = full;
1593  break;
1594  } else if ( fopt.find(gproto) == 0 ) {
1595  fopt.erase(0,gproto.size());
1596  gOptFluxDriver = full;
1597  break;
1598  } else if ( fopt.find(protor) != std::string::npos ) {
1599  gOptFluxDriver = full;
1600  break;
1601  }
1602  }
1603  // tested all cases where user might have specified explicitly
1604  // or been part of an extended file extension
1605  // this is where it gets messy
1606  if ( gOptFluxDriver == "" ) {
1607 
1608  // not specified ? guess from file name itself
1609  if ( fopt.find("gsimple") != std::string::npos ) {
1610  // put dk2nu after gsimple in case simple files are derived from dk2nu
1611  // then both are in name we should choose gsimple
1612  gOptFluxDriver = "genie::flux::GSimpleNtpFlux";
1613  } else if ( fopt.find("dk2nu") != std::string::npos ) {
1614  gOptFluxDriver = "genie::flux::GDk2NuFlux";
1615  } else {
1616  // does it look like the histogram format
1617  const char* hstrings[] = { ",12[", ",+12[", ",-12[",
1618  ",14[", ",+14[", ",-14[",
1619  ",16[", ",+16[", ",-16[" };
1620  size_t nh = sizeof(hstrings)/sizeof(const char*);
1621  for (size_t ih=0; ih<nh; ++ih) {
1622  if ( fopt.find(hstrings[ih]) != std::string::npos ) {
1623  // hey!
1624  gOptFluxDriver = "genie::flux::GCylindTH1Flux";
1625  break;
1626  }
1627  } // loop over possible histogram specifiers
1628  }
1629 
1630  // fall through default ... hope it works
1631  if ( gOptFluxDriver == "" ) {
1632  gOptFluxDriver = "genie::flux::GNuMIFlux";
1633  }
1634  }
1635 
1636  gOptUsingHistFlux = ( gOptFluxDriver == "genie::flux::GCylindTH1Flux" );
1637  if ( gOptUsingHistFlux ) ParseFluxHst(fopt);
1638  else ParseFluxFileConfig(fopt);
1639 }
void ParseFluxFileConfig(string fopt)
map< string, string > gOptFluxShortNames
void ParseFluxHst(string fopt)
string gOptFluxDriver
bool gOptUsingHistFlux
void GetCommandLineArgs ( int  argc,
char **  argv 
)
static void gsSIGTERMhandler ( int  )
static

Definition at line 376 of file gEvGenLArDM.cxx.

References gSigTERM.

Referenced by main().

377 {
378  gSigTERM = true;
379  std::cerr << "Caught SIGTERM" << std::endl;
380 }
bool gSigTERM
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().

731 {
732  /// potentially load extra libraries that might extend the list of
733  /// potential flux drivers, and how to map short names to classes ...
734  // we really should at this point read some external file to get
735  // an expandable list of libraries ... but for now just hard code it
736 
737  vector<string> extraLibs;
738 
739  ///***** this part should come from reading an external file
740  /// placeholder file $GENIE/config/FluxDriverExpansion.xml
741 
742  extraLibs.push_back("libdk2nuTree");
743  extraLibs.push_back("libdk2nuGenie");
744 
745  ///******* done with fake "read"
746 
747  // see if there are any libraries to load
748  // just let ROOT do it ... check error code on return
749  // but tweak ROOT's ERROR message output so we don't see huge messages
750  // for failures
751  // gErrorIgnoreLevel to kError (from TError.h)
752 
753  Int_t prevErrorLevel = gErrorIgnoreLevel;
754  gErrorIgnoreLevel = kFatal;
755  for (size_t ilib=0; ilib < extraLibs.size(); ++ilib) {
756  // library names should be like libXYZZY without extension [e.g. .so]
757  // but with the leading "lib"
758  int loadStatus = gSystem->Load(extraLibs[ilib].c_str());
759  const char* statWords = "failed to load";
760  if ( loadStatus == 0 ) { statWords = "successfully loaded"; }
761  else if ( loadStatus == 1 ) { statWords = "already had loaded"; }
762  else if ( loadStatus == -1 ) { statWords = "couldn't find library"; }
763  else if ( loadStatus == -2 ) { statWords = "mismatched version"; }
764 
765  LOG("gevgen_lardm",pNOTICE)
766  << statWords << " (" << loadStatus << ") " << extraLibs[ilib];
767  }
768  // restore the ROOT error message level
769  gErrorIgnoreLevel = prevErrorLevel;
770 
771 }
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pNOTICE
Definition: Messenger.h:61
int main ( int  argc,
char **  argv 
)

Definition at line 383 of file gEvGenLArDM.cxx.

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::RunOpt::Instance(), genie::PDGLibrary::Instance(), genie::AlgConfigPool::Instance(), genie::flux::GFluxDriverFactory::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::units::s, 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().

384 {
386  GetCommandLineArgs(argc,argv);
387  PDGLibrary::Instance()->AddDarkMatter(gOptDMMass,gOptMedRatio);
388  if (gOptZpCoupling > 0.) {
389  Registry * r = AlgConfigPool::Instance()->CommonList("Param", "BoostedDarkMatter");
390  r->UnLock();
391  r->Set("ZpCoupling", gOptZpCoupling);
392  r->Lock();
393  }
394 
395  if ( ! RunOpt::Instance()->Tune() ) {
396  LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
397  exit(-1);
398  }
399  RunOpt::Instance()->BuildTune();
400 
401  // Initialization of random number generators, cross-section table,
402  // messenger thresholds, cache file
403  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
404  utils::app_init::CacheFile(RunOpt::Instance()->CacheFile());
407 
408  // Set GHEP print level
409  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
410 
411  // *************************************************************************
412  // * Create / configure the geometry driver
413  // *************************************************************************
414  GeomAnalyzerI * geom_driver = 0;
415 
416  if(gOptUsingRootGeom) {
417  //
418  // *** Using a realistic root-based detector geometry description
419  //
420 
421  // creating & configuring a root geometry driver
424  TGeoVolume * topvol = rgeom->GetGeometry()->GetTopVolume();
425  if ( ! topvol ) {
426  LOG("gevgen_numi", pFATAL) << "Null top ROOT geometry volume!";
427  exit(1);
428  }
429  // retrieve the master volume name
430  gOptRootGeomMasterVol = topvol->GetName();
431 
432  rgeom -> SetLengthUnits (gOptGeomLUnits);
433  rgeom -> SetDensityUnits (gOptGeomDUnits);
434  rgeom -> SetTopVolName (gOptRootGeomTopVol); // set user defined "topvol"
435 
436  // getting the bounding box dimensions along z so as to set the
437  // appropriate upstream generation surface for the NuMI flux driver
438 
439  // RWH 2010-07-16: do not try to automatically get zmin from geometry, rather
440  // by default let the flux start from the window. If the user wants to
441  // override this then they need to explicitly set a "zmin". Trying to use
442  // the geometry is fraught with problems in local vs. global coordinates and
443  // units where it can appear to work in some cases but it actually isn't really
444  // universally correct.
445  //was// TGeoShape * bounding_box = topvol->GetShape();
446  //was// bounding_box->GetAxisRange(3, zmin, zmax);
447  //was// zmin *= rgeom->LengthUnits();
448  //was// zmax *= rgeom->LengthUnits();
449 
450  // switch on/off volumes as requested
451  if ( (gOptRootGeomTopVol[0] == '+') || (gOptRootGeomTopVol[0] == '-') ) {
452  bool exhaust = (*gOptRootGeomTopVol.c_str() == '+');
454  }
455 
456  // casting to the GENIE geometry driver interface
457  geom_driver = dynamic_cast<GeomAnalyzerI *> (rgeom);
458 
459  // user specifid a fiducial volume cut ... parse that out
460  if ( gOptFidCut.find("rock") != std::string::npos )
462  else if ( gOptFidCut != "" ) CreateFidSelection(gOptFidCut,rgeom);
463 
464  }
465  else {
466  //
467  // *** Using a 'point' geometry with the specified target mix
468  // *** ( = a list of targets with their corresponding weight fraction)
469  //
470 
471  // creating & configuring a point geometry driver
474  // casting to the GENIE geometry driver interface
475  geom_driver = dynamic_cast<GeomAnalyzerI *> (pgeom);
476  }
477 
478  // *************************************************************************
479  // * Create / configure the flux driver
480  // *************************************************************************
481  GFluxI * flux_driver =
482  genie::flux::GFluxDriverFactory::Instance().GetFluxDriver("genie::flux::GSimpleNtpFlux");
483  if ( ! flux_driver ) {
484  // couldn't get the requested flux driver ?
485  std::ostringstream s;
486  s << "Known FluxDrivers:\n";
487  const std::vector<std::string>& known =
489  std::vector<std::string>::const_iterator itr = known.begin();
490  for ( ; itr != known.end(); ++itr ) s << " " << (*itr) << "\n";
491  LOG("gevgen_lardm", pFATAL)
492  << "Failed to get any flux driver from GFluxDriverFactory";
493  exit(1);
494  }
495 
496  genie::flux::GFluxFileConfigI* flux_file_config =
497  dynamic_cast<genie::flux::GFluxFileConfigI*>(flux_driver);
498  if ( ! flux_file_config ) {
499  LOG("gevgen_lardm", pFATAL)
500  << "Failed to get GFluxFileConfigI driver from GFluxDriverFactory";
501  exit(1);
502  }
503 
504  //
505  // *** Using the detailed ntuple neutrino flux description
506  //
507  flux_file_config->LoadBeamSimData(gOptFluxFile, "");
508  flux_file_config->SetUpstreamZ(gOptZmin); // was "zmin" from bounding_box
509  flux_file_config->SetNumOfCycles(0);
510  PDGCodeList dm_pdg;
511  dm_pdg.push_back(kPdgDarkMatter);
512  flux_file_config->SetFluxParticles(dm_pdg);
513 
514  // these might come in handy ... avoid repeated dynamic_cast<> calls
515  genie::flux::GFluxFileConfigI* fluxFileConfigI =
516  dynamic_cast<genie::flux::GFluxFileConfigI*>(flux_driver);
517 
518 
519  // *************************************************************************
520  // * Handle chicken/egg problem: geom analyzer vs. flux.
521  // * Need both at this point change geom scan defaults.
522  // *************************************************************************
523  if ( gOptUsingRootGeom ) {
524 
526  dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
527  if ( ! rgeom ) assert(0);
528 
529  rgeom -> SetDebugFlags(gOptDebug);
530 
531  // even if user doesn't specify gOptNScan configure to scan using flux
532  if ( gOptNScan >= 0 ) {
533  LOG("gevgen_lardm", pNOTICE)
534  << "Using ROOTGeomAnalyzer: geom scan using flux: nparticles=" << gOptNScan;
535  rgeom->SetScannerFlux(flux_driver);
536  if ( gOptNScan > 0 ) rgeom->SetScannerNParticles(gOptNScan);
537  } else {
538  int nabs = TMath::Abs(gOptNScan);
539  LOG("gevgen_lardm", pNOTICE)
540  << "Using ROOTGeomAnalyzer: geom scan using box: npoints=nrays=" << nabs;
541  rgeom->SetScannerNPoints(nabs);
542  rgeom->SetScannerNRays(nabs);
543  }
544  }
545 
546  // *************************************************************************
547  // * Create/configure the event generation driver
548  // *************************************************************************
549  GMCJDriver * mcj_driver = new GMCJDriver;
550  mcj_driver->SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
551  mcj_driver->UseFluxDriver(flux_driver);
552  mcj_driver->UseGeomAnalyzer(geom_driver);
553  if ( ( gOptExtMaxPlXml != "" ) && ! gOptWriteMaxPlXml ) {
554  mcj_driver->UseMaxPathLengths(gOptExtMaxPlXml);
555  }
556  mcj_driver->Configure();
557  mcj_driver->UseSplines();
558  mcj_driver->ForceSingleProbScale();
559 
560  if ( ( gOptExtMaxPlXml != "" ) && gOptWriteMaxPlXml ) {
562  dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
563  if ( rgeom ) {
564  const genie::PathLengthList& maxpath = rgeom->GetMaxPathLengths();
565  std::string maxplfile = gOptExtMaxPlXml;
566  maxpath.SaveAsXml(maxplfile);
567  // append extra info to file
568  std::ofstream mpfile(maxplfile.c_str(), std::ios_base::app);
569  mpfile
570  << std::endl
571  << "<!-- this file is only relevant for a setup compatible with:"
572  << std::endl
573  << "geom: " << gOptRootGeom << " top: \"" << gOptRootGeomTopVol << "\""
574  << std::endl
575  << "flux: " << gOptFluxFile
576  << std::endl
577  << "fidcut: " << gOptFidCut
578  << std::endl
579  << "nscan: " << gOptNScan << " (0>= use flux, <0 use box |nscan| points/rays)"
580  << std::endl
581  << "zmin: " << gOptZmin << " (if |zmin| > 1e30, leave ray on flux window)"
582  << std::endl
583  << "-->"
584  << std::endl;
585  mpfile.close();
586  }
587  }
588 
589  // *************************************************************************
590  // * Prepare for writing the output event tree & status file
591  // *************************************************************************
592 
593  // Initialize an Ntuple Writer to save GHEP records into a TTree
595  ntpw.CustomizeFilenamePrefix(gOptEvFilePrefix);
596  ntpw.Initialize();
597 
598 
599  std::vector<TBranch*> extraBranches;
600  std::vector<std::string> branchNames;
601  std::vector<std::string> branchClassNames;
602  std::vector<void**> branchObjPointers;
603 
604  // Add custom branch(s) to the standard GENIE event tree so that
605  // info on the flux neutrino parent particle can be passed-through
606  if ( fluxFileConfigI ) {
607  fluxFileConfigI->GetBranchInfo(branchNames,branchClassNames,
608  branchObjPointers);
609  size_t nn = branchNames.size();
610  size_t nc = branchClassNames.size();
611  size_t np = branchObjPointers.size();
612  if ( nn != nc || nc != np ) {
613  LOG("gevgen_lardm", pERROR)
614  << "Inconsistent info back from flux driver "
615  << "for branch info: " << nn << " " << nc << " " << np;
616  } else {
617  for (size_t ii = 0; ii < nn; ++ii) {
618  const char* bname = branchNames[ii].c_str();
619  const char* cname = branchClassNames[ii].c_str();
620  void**& optr = branchObjPointers[ii]; // note critical '&' !
621  if ( ! optr || ! *optr ) continue; // no pointer supplied, skip it
622  int split = 99; // 1
623  LOG("gevgen_lardm", pNOTICE)
624  << "Adding extra branch \"" << bname << "\" of type \""
625  << cname << "\" (" << optr << ") to output tree";
626  TBranch* bptr = ntpw.EventTree()->Branch(bname,cname,optr,32000,split);
627  extraBranches.push_back(bptr);
628 
629  if ( bptr ) {
630  // don't delete this !!! we're sharing
631  bptr->SetAutoDelete(false);
632  } else {
633  LOG("gevgen_lardm", pERROR)
634  << "FAILED to add extra branch \"" << bname << "\" of type \""
635  << cname << "\" to output tree";
636  }
637  } // loop over additions
638  } // same # of entries
639  } // of genie::flux::GFluxFileConfigI type
640 
641  // Create a MC job monitor for a periodically updated status file
642  GMCJMonitor mcjmonitor(gOptRunNu);
643  mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());
644 
645  // *************************************************************************
646  // * Event generation loop
647  // *************************************************************************
648 
649  // define handler to allow signal to end job gracefully
650  signal(SIGTERM,gsSIGTERMhandler);
651 
652  int ievent = 0;
653  while ( ! gSigTERM )
654  {
655  LOG("gevgen_lardm", pINFO)
656  << " *** Generating event............ " << ievent;
657 
658  // In case the required statistics was expressed as 'number of events'
659  // then quit if that number has been generated
660  if ( ievent == gOptNev ) break;
661 
662  // Generate a single event using neutrinos coming from the specified flux
663  // and hitting the specified geometry or target mix
664  EventRecord * event = mcj_driver->GenerateEvent();
665 
666  // Check whether a null event was returned due to the flux driver reaching
667  // the end of the input flux ntuple - exit the event generation loop
668  if ( ! event && flux_driver->End() ) {
669  LOG("gevgen_lardm", pWARN)
670  << "** The flux driver read all the input flux entries: End()==true";
671  break;
672  }
673  if ( ! event ) {
674  LOG("gevgen_lardm", pERROR)
675  << "Got a null generated neutino event! Retrying ...";
676  continue;
677  }
678  LOG("gevgen_lardm", pINFO)
679  << "Generated event: " << *event;
680 
681  // A valid event was generated: flux info (parent decay/prod
682  // position/kinematics) for that simulated event should already
683  // be connected to the right output tree branch
684 
685  // Add event at the output ntuple, refresh the mc job monitor & clean-up
686  ntpw.AddEventRecord(ievent, event);
687  mcjmonitor.Update(ievent,event);
688  delete event;
689  ievent++;
690 
691  } //1
692 
693  // Copy metadata tree, if available
694  if ( fluxFileConfigI ) {
695  TTree* t1 = fluxFileConfigI->GetMetaDataTree();
696  if ( t1 ) {
697  size_t nmeta = t1->GetEntries();
698  TTree* t2 = (TTree*)t1->Clone(0);
699  for (size_t i = 0; i < nmeta; ++i) {
700  t1->GetEntry(i);
701  t2->Fill();
702  }
703  t2->Write();
704  }
705  }
706 
707  LOG("gevgen_lardm", pINFO)
708  << "The GENIE MC job is done generating events - Cleaning up & exiting...";
709 
710  // *************************************************************************
711  // * Save & clean-up
712  // *************************************************************************
713 
714  // Save the generated event tree & close the output file
715  ntpw.Save();
716 
717  // Clean-up
718  delete geom_driver;
719  delete flux_driver;
720  delete mcj_driver;
721  // this list should only be histograms that have (for some reason)
722  // not been handed over to the GCylindTH1Flux driver.
723 
724  LOG("gevgen_lardm", pNOTICE) << "Done!";
725 
726  return 0;
727 }
void RandGen(long int seed)
Definition: AppInit.cxx:30
void XSecTable(string inpfile, bool require_table)
Definition: AppInit.cxx:38
void CreateRockBoxSelection(string fidcut, GeomAnalyzerI *geom_driver)
virtual void LoadBeamSimData(const std::vector< std::string > &filenames, const std::string &det_loc)=0
string gOptEvFilePrefix
Definition: gAtmoEvGen.cxx:310
#define pERROR
Definition: Messenger.h:59
double gOptDMMass
Definition: gEvGenDM.cxx:221
int gOptDebug
string gOptRootGeomMasterVol
void SetEventGeneratorList(string listname)
Definition: GMCJDriver.cxx:66
double gOptMedRatio
Definition: gEvGenDM.cxx:224
#define pFATAL
Definition: Messenger.h:56
double gOptGeomDUnits
Definition: gAtmoEvGen.cxx:303
const int kPdgDarkMatter
Definition: PDGCodes.h:218
virtual const PathLengthList & GetMaxPathLengths(void) const
static constexpr double s
Definition: Units.h:95
double gOptZpCoupling
Definition: gEvGenDM.cxx:222
double gOptZmin
const std::vector< std::string > & AvailableFluxDrivers() const
int gOptNScan
map< int, double > gOptTgtMix
Definition: gAtmoEvGen.cxx:299
string gOptFidCut
virtual void SetNumOfCycles(long int ncycle)
limit cycling through input files
virtual void SetScannerNParticles(int np)
A list of PDG codes.
Definition: PDGCodeList.h:32
bool gOptWriteMaxPlXml
void CreateFidSelection(string fidcut, GeomAnalyzerI *geom_driver)
virtual void SetUpstreamZ(double z0)
void UseFluxDriver(GFluxI *flux)
Definition: GMCJDriver.cxx:83
Object to be filled with the neutrino path-length, for all detector geometry materials, when starting from a position x and travelling along the direction of the neutrino 4-momentum.
Simple class to create &amp; update MC job status files and env. vars. This is used to be able to keep tr...
Definition: GMCJMonitor.h:31
A GENIE `MC Job Driver&#39;. Can be used for setting up complicated event generation cases involving deta...
Definition: GMCJDriver.h:46
void SaveAsXml(string filename) const
string gOptInpXSecFile
Definition: gAtmoEvGen.cxx:313
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
string gOptRootGeomTopVol
Definition: gAtmoEvGen.cxx:301
void ForceSingleProbScale(void)
Definition: GMCJDriver.cxx:172
virtual void SetScannerNRays(int nr)
bool UseMaxPathLengths(string xml_filename)
Definition: GMCJDriver.cxx:99
bool gSigTERM
int gOptNev
Definition: gAtmoEvGen.cxx:305
static void gsSIGTERMhandler(int)
void Configure(bool calc_prob_scales=true)
Definition: GMCJDriver.cxx:399
A ROOT/GEANT4 geometry driver.
#define pINFO
Definition: Messenger.h:62
void Lock(void)
locks the registry
Definition: Registry.cxx:148
string gOptRootGeom
Definition: gAtmoEvGen.cxx:300
#define pWARN
Definition: Messenger.h:60
EventRecord * GenerateEvent(void)
Definition: GMCJDriver.cxx:815
virtual void GetBranchInfo(std::vector< std::string > &branchNames, std::vector< std::string > &branchClassNames, std::vector< void ** > &branchObjPointers)
Long_t gOptRunNu
Definition: gAtmoEvGen.cxx:295
void UnLock(void)
unlocks the registry (doesn&#39;t unlock items)
Definition: Registry.cxx:153
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:37
void LoadExtraOptions(void)
A utility class to facilitate creating the GENIE MC Ntuple from the output GENIE GHEP event records...
Definition: NtpWriter.h:39
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
virtual void SetScannerNPoints(int np)
set geometry driver&#39;s configuration options
virtual void SetScannerFlux(GFluxI *f)
static GFluxDriverFactory & Instance()
string gOptFluxFile
virtual bool End(void)=0
true if no more flux nu&#39;s can be thrown (eg reaching end of beam sim ntuples)
NtpMCFormat_t kDefOptNtpFormat
Definition: gAtmoEvGen.cxx:319
A vector of EventGeneratorI objects.
virtual void SetFluxParticles(const PDGCodeList &particles)
specify list of flux neutrino species
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:99
#define pNOTICE
Definition: Messenger.h:61
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:563
Defines the GENIE Geometry Analyzer Interface.
Definition: GeomAnalyzerI.h:29
void Set(RgIMapPair entry)
Definition: Registry.cxx:267
void RecursiveExhaust(TGeoVolume *topvol, string volnames, bool exhaust)
Definition: GeoUtils.cxx:16
double gOptGeomLUnits
Definition: gAtmoEvGen.cxx:302
void UseGeomAnalyzer(GeomAnalyzerI *geom)
Definition: GMCJDriver.cxx:88
void UseSplines(bool useLogE=true)
Definition: GMCJDriver.cxx:93
The PointGeomAnalyzer class is the simplest implementation of the GeomAnalyserI interface and defines...
void CacheFile(string inpfile)
Definition: AppInit.cxx:117
string gOptExtMaxPlXml
Definition: gAtmoEvGen.cxx:304
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
bool gOptUsingRootGeom
Definition: gAtmoEvGen.cxx:298
virtual TGeoManager * GetGeometry(void) const
genie::GFluxI * GetFluxDriver(const std::string &)
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:29
long int gOptRanSeed
Definition: gAtmoEvGen.cxx:312
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().

1814 {
1815  // Using gnumi/gsimple/dk2nu beam flux ntuples
1816  // Extract beam flux (root) file name & detector location
1817  //
1818  vector<string> fluxv = utils::str::Split(flux,",");
1819  if(fluxv.size()<2) {
1820  LOG("gevgen_fnal", pFATAL)
1821  << "You need to specify both a flux ntuple ROOT file "
1822  << " _AND_ a detector location";
1823  PrintSyntax();
1824  exit(1);
1825  }
1826  gOptFluxFile = fluxv[0];
1827  gOptDetectorLocation = fluxv[1];
1828 
1829  for ( size_t j = 2; j < fluxv.size(); ++j ) {
1830  int ipdg = atoi(fluxv[j].c_str());
1831  gOptFluxPdg.push_back(ipdg);
1832  }
1833 
1834 }
PDGCodeList gOptFluxPdg
#define pFATAL
Definition: Messenger.h:56
string gOptDetectorLocation
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
string gOptFluxFile
void PrintSyntax(void)
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:58
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().

1642 {
1643  // Using flux from histograms
1644  // Extract the root file name & the list of histogram names & neutrino
1645  // species (specified as 'filename,histo1[species1],histo2[species2],...')
1646  // for variable width histograms, default to multiply in the width
1647  // histo1[species1]WIDTH = multiply in the width
1648  // histo1[species1]NOWIDTH = don't multiply in the width
1649  // possibly with configuration ";r=1.1,dir=(0.1,0.2,0.3),spot=(-1,-2,-3)"
1650  // appended
1651  // See documentation on top section of this file.
1652 
1653  vector<string> fluxtop = utils::str::Split(flux,";");
1654  string beamsettings;
1655  string histosettings;
1656  if ( fluxtop.size() == 1 ) {
1657  histosettings = fluxtop[0];
1658  LOG("gevgen_fnal", pNOTICE)
1659  << "ParseFluxHst: no settings for radius, direction, spot position"
1660  << " using defaults, r=" << gOptBeamRadius;
1661  } else if ( fluxtop.size() > 2 ) {
1662  LOG("gevgen_fnal", pFATAL)
1663  << "ParseFluxHst: too many ';' separated fields";
1664  PrintSyntax();
1665  exit(1);
1666  } else {
1667  // decide which is which depending on which has "=" in it
1668  string::size_type eqpos0 = fluxtop[0].find("=");
1669  string::size_type eqpos1 = fluxtop[1].find("=");
1670  if (eqpos0 == string::npos && eqpos1 != string::npos ) {
1671  histosettings = fluxtop[0];
1672  beamsettings = fluxtop[1];
1673  } else if (eqpos0 != string::npos && eqpos1 == string::npos ) {
1674  beamsettings = fluxtop[0];
1675  histosettings = fluxtop[1];
1676  } else {
1677  LOG("gevgen_fnal", pFATAL)
1678  << "ParseFluxHst: too many / not enough fields with '=' "
1679  << "\n" << fluxtop[0] << "\n" << fluxtop[1];
1680  PrintSyntax();
1681  exit(1);
1682  }
1683  // now parse the string we have
1684  double r=-1, dx=0, dy=0, dz=1, x=0, y=0, z=0;
1685  int nscan = sscanf(beamsettings.c_str(),
1686  "r=%lf,dir=(%lf,%lf,%lf),spot=(%lf,%lf,%lf)",
1687  &r,&dx,&dy,&dz,&x,&y,&z);
1688  cout << "nscan = " << nscan << endl;
1689  gOptBeamRadius = r;
1690  gOptBeamDir = TVector3(dx,dy,dz);
1691  gOptBeamSpot = TVector3(x,y,z);
1692 
1693  }
1694  LOG("gevgen_fnal", pNOTICE)
1695  << "ParseFluxHst: "
1696  << " using r=" << gOptBeamRadius
1697  << ",dir=("
1698  << gOptBeamDir.Px() << ","
1699  << gOptBeamDir.Py() << ","
1700  << gOptBeamDir.Pz() << ")"
1701  << ",spot=("
1702  << gOptBeamSpot.X() << ","
1703  << gOptBeamSpot.Y() << ","
1704  << gOptBeamSpot.Z() << ")";
1705 
1706 
1707  vector<string> fluxv = utils::str::Split(histosettings,",");
1708  if(fluxv.size()<2) {
1709  LOG("gevgen_fnal", pFATAL)
1710  << "You need to specify both a flux histogram ROOT file "
1711  << " _AND_ at least one histogram[pdg] mapping";
1712  PrintSyntax();
1713  exit(1);
1714  }
1715  gOptFluxFile = fluxv[0];
1716  bool accessible_flux_file = !(gSystem->AccessPathName(gOptFluxFile.c_str()));
1717  if (!accessible_flux_file) {
1718  LOG("gevgen_fnal", pFATAL)
1719  << "Can not access flux file: " << gOptFluxFile;
1720  PrintSyntax();
1721  exit(1);
1722  }
1723  // Extract energy spectra for all specified neutrino species
1724  TFile flux_file(gOptFluxFile.c_str(), "read");
1725  for(unsigned int inu=1; inu<fluxv.size(); inu++) {
1726  string nutype_and_histo = fluxv[inu];
1727  string::size_type open_bracket = nutype_and_histo.find("[");
1728  string::size_type close_bracket = nutype_and_histo.find("]");
1729  if (open_bracket ==string::npos ||
1730  close_bracket==string::npos)
1731  {
1732  LOG("gevgen_fnal", pFATAL)
1733  << "You made an error in specifying the flux histograms";
1734  PrintSyntax();
1735  exit(1);
1736  }
1737  string::size_type ibeg = 0;
1738  string::size_type iend = open_bracket;
1739  string::size_type jbeg = open_bracket+1;
1740  string::size_type jend = close_bracket;
1741  string nutype = nutype_and_histo.substr(ibeg,iend-ibeg);
1742  string histo = nutype_and_histo.substr(jbeg,jend-jbeg);
1743  string extra = nutype_and_histo.substr(jend+1,string::npos);
1744  std::transform(extra.begin(),extra.end(),extra.begin(),::toupper);
1745  LOG("gevgen_fnal", pNOTICE) // pDEBUG
1746  << " =======> nutype " << nutype << " histo " << histo << " extra " << extra;
1747  // access specified histogram from the input root file
1748  TH1D * ihst = (TH1D*) flux_file.Get(histo.c_str());
1749  if(!ihst) {
1750  LOG("gevgen_fnal", pFATAL)
1751  << "Can not find histogram: " << histo
1752  << " in flux file: " << gOptFluxFile;
1753  PrintSyntax();
1754  exit(1);
1755  }
1756 
1757  // Copy in the flux histogram from the root file
1758  // use Clone rather than assuming fix bin widths and rebooking
1759  TH1D* spectrum = (TH1D*)ihst->Clone();
1760  spectrum->SetNameTitle("spectrum","neutrino_flux");
1761  spectrum->SetDirectory(0);
1762  // get rid of original
1763  delete ihst;
1764 
1765  bool force_binwidth = false;
1766 #if ROOT_VERSION_CODE <= ROOT_VERSION(9,99,99)
1767  // GetRandom() sampling on variable bin width histograms does not
1768  // correctly account for bin widths for all versions of ROOT prior
1769  // to (currently forever). At some point this might change and
1770  // the necessity of this code snippet will go away
1771  TAxis* xaxis = spectrum->GetXaxis();
1772  if (xaxis->IsVariableBinSize()) force_binwidth = true;
1773 #endif
1774  if ( extra == "WIDTH" ) force_binwidth = true;
1775  if ( extra == "NOWIDTH" ) force_binwidth = false;
1776  if ( force_binwidth ) {
1777  LOG("gevgen_fnal", pNOTICE)
1778  << "multiplying by bin width for histogram " << histo
1779  << " as " << spectrum->GetName() << " for nutype " << nutype
1780  << " from " << gOptFluxFile;
1781  for(int ibin = 1; ibin <= spectrum->GetNbinsX(); ++ibin) {
1782  double data = spectrum->GetBinContent(ibin);
1783  double width = spectrum->GetBinWidth(ibin);
1784  spectrum->SetBinContent(ibin,data*width);
1785  }
1786  }
1787 
1788  // convert neutrino name -> pdg code
1789  int pdg = atoi(nutype.c_str());
1790  if(!pdg::IsNeutrino(pdg) && !pdg::IsAntiNeutrino(pdg)) {
1791  LOG("gevgen_fnal", pFATAL)
1792  << "Unknown neutrino type: " << nutype;
1793  PrintSyntax();
1794  exit(1);
1795  }
1796  // store flux neutrino code / energy spectrum
1797  LOG("gevgen_fnal", pDEBUG)
1798  << "Adding energy spectrum for flux neutrino: pdg = " << pdg;
1799  gOptFluxHst.insert(map<int, TH1D*>::value_type(pdg, spectrum));
1800  }//inu
1801 
1802  if(gOptFluxHst.size()<1) {
1803  LOG("gevgen_fnal", pFATAL)
1804  << "You have not specified any flux histogram!";
1805  PrintSyntax();
1806  exit(1);
1807  }
1808 
1809  flux_file.Close();
1810 }
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
#define pFATAL
Definition: Messenger.h:56
TVector3 gOptBeamDir
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
TVector3 gOptBeamSpot
Double_t gOptBeamRadius
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:36
string gOptFluxFile
#define pNOTICE
Definition: Messenger.h:61
void PrintSyntax(void)
map< int, TH1D * > gOptFluxHst
#define pDEBUG
Definition: Messenger.h:63
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.