Functions | Variables
gT2KEvGen.cxx File Reference
#include <cassert>
#include <cstdlib>
#include <string>
#include <sstream>
#include <vector>
#include <map>
#include <TSystem.h>
#include <TTree.h>
#include <TFile.h>
#include <TH1D.h>
#include <TMath.h>
#include <TGeoVolume.h>
#include <TGeoShape.h>
#include <TList.h>
#include <TObject.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/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/Numerical/RandomGen.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/T2KEvGenMetaData.h"
#include "Framework/Utils/SystemUtils.h"
#include "Framework/Utils/PrintUtils.h"

Go to the source code of this file.

Functions

void GetCommandLineArgs (int argc, char **argv)
 
void PrintSyntax (void)
 
int main (int argc, char **argv)
 

Variables

string kDefOptGeomLUnits = "mm"
 
string kDefOptGeomDUnits = "g_cm3"
 
NtpMCFormat_t kDefOptNtpFormat = kNFGHEP
 
double kDefOptFluxNorm = 1E+21
 
string kDefOptEvFilePrefix = "gntp"
 
Long_t gOptRunNu
 
bool gOptUsingRootGeom = false
 
bool gOptUsingHistFlux = false
 
map< int, double > gOptTgtMix
 
map< int, TH1D * > gOptFluxHst
 
string gOptRootGeom
 
string gOptRootGeomTopVol = ""
 
double gOptGeomLUnits = 0
 
double gOptGeomDUnits = 0
 
string gOptExtMaxPlXml
 
string gOptFluxFile
 
string gOptDetectorLocation
 
double gOptFluxNorm
 
PDGCodeList gOptFluxNtpNuList (false)
 
int gOptFluxNCycles
 
int gOptNev
 
double gOptPOT
 
bool gOptExitAtEndOfFullFluxCycles
 
string gOptEvFilePrefix
 
bool gOptUseFluxProbs = false
 
bool gOptSaveFluxProbsFile = false
 
string gOptFluxProbFileName
 
string gOptSaveFluxProbsFileName
 
bool gOptRandomFluxOffset = false
 
long int gOptRanSeed
 
string gOptInpXSecFile
 

Function Documentation

void GetCommandLineArgs ( int  argc,
char **  argv 
)

Definition at line 936 of file gT2KEvGen.cxx.

References genie::CmdLnArgParser::ArgAsDouble(), genie::CmdLnArgParser::ArgAsInt(), genie::CmdLnArgParser::ArgAsLong(), genie::CmdLnArgParser::ArgAsString(), genie::RunOpt::EnableBareXSecPreCalc(), exit(), genie::PDGLibrary::Find(), flux, geom(), gOptDetectorLocation, gOptEvFilePrefix, gOptExitAtEndOfFullFluxCycles, gOptExtMaxPlXml, gOptFluxFile, gOptFluxHst, gOptFluxNCycles, gOptFluxNorm, gOptFluxNtpNuList, gOptFluxProbFileName, gOptGeomDUnits, gOptGeomLUnits, gOptInpXSecFile, gOptNev, gOptPOT, gOptRandomFluxOffset, gOptRanSeed, gOptRootGeom, gOptRootGeomTopVol, gOptRunNu, gOptSaveFluxProbsFile, gOptSaveFluxProbsFileName, gOptTgtMix, gOptUseFluxProbs, gOptUsingHistFlux, gOptUsingRootGeom, plot_validation_datamc::help, histo, make_syst_table_plots::ibin, genie::PDGLibrary::Instance(), genie::RunOpt::Instance(), genie::pdg::IsAntiNeutrino(), genie::pdg::IsNeutrino(), kDefOptEvFilePrefix, kDefOptFluxNorm, kDefOptGeomDUnits, kDefOptGeomLUnits, LOG, nutype(), genie::CmdLnArgParser::OptionExists(), plot_validation_datamc::parser, pDEBUG, make_root_from_grid_output::pdg, pFATAL, pINFO, pNOTICE, genie::utils::print::PrintFramedMesg(), PrintSyntax(), genie::PDGCodeList::push_back(), pWARN, genie::RunOpt::ReadFromCommandLine(), genie::utils::str::Split(), genie::utils::units::UnitFromString(), and wgt.

Referenced by main().

937 {
938  LOG("gevgen_t2k", pINFO) << "Parsing command line arguments";
939 
940  // Common run options. Set defaults and read.
941  RunOpt::Instance()->EnableBareXSecPreCalc(true);
942  RunOpt::Instance()->ReadFromCommandLine(argc,argv);
943 
944  // Parse run options for this app
945 
946  CmdLnArgParser parser(argc,argv);
947 
948  // help?
949  bool help = parser.OptionExists('h');
950  if(help) {
951  PrintSyntax();
952  exit(0);
953  }
954 
955  // run number:
956  if( parser.OptionExists('r') ) {
957  LOG("gevgen_t2k", pDEBUG) << "Reading MC run number";
958  gOptRunNu = parser.ArgAsLong('r');
959  } else {
960  LOG("gevgen_t2k", pDEBUG) << "Unspecified run number - Using default";
961  gOptRunNu = 0;
962  } //-r
963 
964  //
965  // *** geometry
966  //
967 
968  string geom = "";
969  string lunits, dunits;
970  if( parser.OptionExists('g') ) {
971  LOG("gevgen_t2k", pDEBUG) << "Getting input geometry";
972  geom = parser.ArgAsString('g');
973 
974  // is it a ROOT file that contains a ROOT geometry?
975  bool accessible_geom_file =
976  ! (gSystem->AccessPathName(geom.c_str()));
977  if (accessible_geom_file) {
978  gOptRootGeom = geom;
979  gOptUsingRootGeom = true;
980  }
981  } else {
982  LOG("gevgen_t2k", pFATAL)
983  << "No geometry option specified - Exiting";
984  PrintSyntax();
985  exit(1);
986  } //-g
987 
988  if(gOptUsingRootGeom) {
989  // using a ROOT geometry - get requested geometry units
990 
991  // legth units:
992  if( parser.OptionExists('L') ) {
993  LOG("gevgen_t2k", pDEBUG)
994  << "Checking for input geometry length units";
995  lunits = parser.ArgAsString('L');
996  } else {
997  LOG("gevgen_t2k", pDEBUG) << "Using default geometry length units";
998  lunits = kDefOptGeomLUnits;
999  } // -L
1000  // density units:
1001  if( parser.OptionExists('D') ) {
1002  LOG("gevgen_t2k", pDEBUG)
1003  << "Checking for input geometry density units";
1004  dunits = parser.ArgAsString('D');
1005  } else {
1006  LOG("gevgen_t2k", pDEBUG) << "Using default geometry density units";
1007  dunits = kDefOptGeomDUnits;
1008  } // -D
1011 
1012  // check whether an event generation volume name has been
1013  // specified -- default is the 'top volume'
1014  if( parser.OptionExists('t') ) {
1015  LOG("gevgen_t2k", pDEBUG) << "Checking for input volume name";
1016  gOptRootGeomTopVol = parser.ArgAsString('t');
1017  } else {
1018  LOG("gevgen_t2k", pDEBUG) << "Using the <master volume>";
1019  } // -t
1020 
1021  // check whether an XML file with the maximum (density weighted)
1022  // path lengths for each detector material is specified -
1023  // otherwise will compute the max path lengths at job init
1024  if( parser.OptionExists('m') ) {
1025  LOG("gevgen_t2k", pDEBUG)
1026  << "Checking for maximum path lengths XML file";
1027  gOptExtMaxPlXml = parser.ArgAsString('m');
1028  } else {
1029  LOG("gevgen_t2k", pDEBUG)
1030  << "Will compute the maximum path lengths at job init";
1031  gOptExtMaxPlXml = "";
1032  } // -m
1033  } // using root geom?
1034 
1035  else {
1036  // User has specified a target mix.
1037  // Decode the list of target pdf codes & their corresponding weight fraction
1038  // (specified as 'pdg_code_1[fraction_1],pdg_code_2[fraction_2],...')
1039  // See documentation on top section of this file.
1040  //
1041  gOptTgtMix.clear();
1042  vector<string> tgtmix = utils::str::Split(geom,",");
1043  if(tgtmix.size()==1) {
1044  int pdg = atoi(tgtmix[0].c_str());
1045  double wgt = 1.0;
1046  gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
1047  } else {
1048  vector<string>::const_iterator tgtmix_iter = tgtmix.begin();
1049  for( ; tgtmix_iter != tgtmix.end(); ++tgtmix_iter) {
1050  string tgt_with_wgt = *tgtmix_iter;
1051  string::size_type open_bracket = tgt_with_wgt.find("[");
1052  string::size_type close_bracket = tgt_with_wgt.find("]");
1053  if (open_bracket ==string::npos ||
1054  close_bracket==string::npos)
1055  {
1056  LOG("gevgen_t2k", pFATAL)
1057  << "You made an error in specifying the target mix";
1058  PrintSyntax();
1059  exit(1);
1060  }
1061  string::size_type ibeg = 0;
1062  string::size_type iend = open_bracket;
1063  string::size_type jbeg = open_bracket+1;
1064  string::size_type jend = close_bracket;
1065  int pdg = atoi(tgt_with_wgt.substr(ibeg,iend-ibeg).c_str());
1066  double wgt = atof(tgt_with_wgt.substr(jbeg,jend-jbeg).c_str());
1067  LOG("gevgen_t2k", pDEBUG)
1068  << "Adding to target mix: pdg = " << pdg << ", wgt = " << wgt;
1069  gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
1070 
1071  }// tgtmix_iter
1072  } // >1 materials in mix
1073  } // using tgt mix?
1074 
1075  //
1076  // *** flux
1077  //
1078 
1079  if( parser.OptionExists('f') ) {
1080  LOG("gevgen_t2k", pDEBUG) << "Getting input flux";
1081  string flux = parser.ArgAsString('f');
1082  gOptUsingHistFlux = (flux.find("[") != string::npos);
1083 
1084  if(!gOptUsingHistFlux) {
1085  // Using JNUBEAM flux ntuples
1086  // Extract JNUBEAM flux (root) file name & detector location.
1087  // Also extract the list of JNUBEAM neutrinos to consider (if none
1088  // is specified here then I will consider all flavours)
1089  //
1090  vector<string> fluxv = utils::str::Split(flux,",");
1091  if(fluxv.size()<2) {
1092  LOG("gevgen_t2k", pFATAL)
1093  << "You need to specify both a flux ntuple ROOT file "
1094  << " _AND_ a detector location";
1095  PrintSyntax();
1096  exit(1);
1097  }
1098  gOptFluxFile = fluxv[0];
1099  gOptDetectorLocation = fluxv[1];
1100 
1101  // Extract the list of neutrinos to consider (if any).
1102  //
1103  for(unsigned int inu = 2; inu < fluxv.size(); inu++)
1104  {
1105  gOptFluxNtpNuList.push_back( atoi(fluxv[inu].c_str()) );
1106  }
1107 
1108  } else {
1109  // Using flux from histograms
1110  // Extract the root file name & the list of histogram names & neutrino
1111  // species (specified as 'filename,histo1[species1],histo2[species2],...')
1112  // See documentation on top section of this file.
1113  //
1114  vector<string> fluxv = utils::str::Split(flux,",");
1115  if(fluxv.size()<2) {
1116  LOG("gevgen_t2k", pFATAL)
1117  << "You need to specify both a flux ntuple ROOT file "
1118  << " _AND_ a detector location";
1119  PrintSyntax();
1120  exit(1);
1121  }
1122  gOptFluxFile = fluxv[0];
1123  bool accessible_flux_file =
1124  !(gSystem->AccessPathName(gOptFluxFile.c_str()));
1125  if (!accessible_flux_file) {
1126  LOG("gevgen_t2k", pFATAL)
1127  << "Can not access flux file: " << gOptFluxFile;
1128  PrintSyntax();
1129  exit(1);
1130  }
1131  // Extract energy spectra for all specified neutrino species
1132  TFile flux_file(gOptFluxFile.c_str(), "read");
1133  for(unsigned int inu=1; inu<fluxv.size(); inu++) {
1134  string nutype_and_histo = fluxv[inu];
1135  string::size_type open_bracket = nutype_and_histo.find("[");
1136  string::size_type close_bracket = nutype_and_histo.find("]");
1137  if (open_bracket ==string::npos ||
1138  close_bracket==string::npos)
1139  {
1140  LOG("gevgen_t2k", pFATAL)
1141  << "You made an error in specifying the flux histograms";
1142  PrintSyntax();
1143  exit(1);
1144  }
1145  string::size_type ibeg = 0;
1146  string::size_type iend = open_bracket;
1147  string::size_type jbeg = open_bracket+1;
1148  string::size_type jend = close_bracket;
1149  string nutype = nutype_and_histo.substr(ibeg,iend-ibeg);
1150  string histo = nutype_and_histo.substr(jbeg,jend-jbeg);
1151  // access specified histogram from the input root file
1152  TH1D * ihst = (TH1D*) flux_file.Get(histo.c_str());
1153  if(!ihst) {
1154  LOG("gevgen_t2k", pFATAL)
1155  << "Can not find histogram: " << histo
1156  << " in flux file: " << gOptFluxFile;
1157  PrintSyntax();
1158  exit(1);
1159  }
1160  // create a local copy of the input histogram
1161  TString origname = ihst->GetName();
1162  TString tmpname; tmpname.Form("%s_", origname.Data());
1163  TH1D * spectrum = new TH1D(
1164  tmpname.Data(), ihst->GetName(), ihst->GetNbinsX(),
1165  ihst->GetXaxis()->GetXmin(), ihst->GetXaxis()->GetXmax());
1166  spectrum->SetDirectory(0);
1167  for(int ibin = 1; ibin <= ihst->GetNbinsX(); ibin++) {
1168  spectrum->SetBinContent(ibin, ihst->GetBinContent(ibin));
1169  }
1170  // get rid of original
1171  delete ihst;
1172  // rename copy
1173  spectrum->SetName(origname.Data());
1174 
1175  // convert neutrino name -> pdg code
1176  int pdg = atoi(nutype.c_str());
1177  if(!pdg::IsNeutrino(pdg) && !pdg::IsAntiNeutrino(pdg)) {
1178  LOG("gevgen_t2k", pFATAL)
1179  << "Unknown neutrino type: " << nutype;
1180  PrintSyntax();
1181  exit(1);
1182  }
1183  // store flux neutrino code / energy spectrum
1184  LOG("gevgen_t2k", pDEBUG)
1185  << "Adding energy spectrum for flux neutrino: pdg = " << pdg;
1186  gOptFluxHst.insert(map<int, TH1D*>::value_type(pdg, spectrum));
1187  }//inu
1188  if(gOptFluxHst.size()<1) {
1189  LOG("gevgen_t2k", pFATAL)
1190  << "You have not specified any flux histogram!";
1191  PrintSyntax();
1192  exit(1);
1193  }
1194  flux_file.Close();
1195  } // flux from histograms or from JNUBEAM ntuples?
1196 
1197  } else {
1198  LOG("gevgen_t2k", pFATAL) << "No flux info was specified - Exiting";
1199  PrintSyntax();
1200  exit(1);
1201  }
1202 
1203  // Use a random offset when looping over flux entries
1204  if(parser.OptionExists('R')) gOptRandomFluxOffset = true;
1205 
1206  //
1207  // *** pre-calculated flux interaction probabilities
1208  //
1209 
1210  // using pre-calculated flux interaction probabilities
1211  if( parser.OptionExists('P') ){
1212  gOptFluxProbFileName = parser.ArgAsString('P');
1213  if(gOptFluxProbFileName.length() > 0){
1214  gOptUseFluxProbs = true;
1215  bool accessible =
1216  !(gSystem->AccessPathName(gOptFluxProbFileName.c_str()));
1217  if(!accessible){
1218  LOG("gevgen_t2k", pFATAL)
1219  << "Can not access pre-calculated flux probabilities file: " << gOptFluxProbFileName;
1220  PrintSyntax();
1221  exit(1);
1222  }
1223  }
1224  else {
1225  LOG("gevgen_t2k", pFATAL)
1226  << "No flux interaction probabilites were specified - exiting";
1227  PrintSyntax();
1228  exit(1);
1229  }
1230  }
1231 
1232  // pre-generating interaction probs and saving to output file
1233  if( parser.OptionExists('S') ){
1234  gOptSaveFluxProbsFile = true;
1235  gOptSaveFluxProbsFileName = parser.ArgAsString('S');
1236  }
1237 
1238  // cannot save and run at the same time
1240  LOG("gevgen_t2k", pFATAL)
1241  << "Cannot specify both the -P and -S options at the same time!";
1242  exit(1);
1243  }
1244 
1245  // only makes sense to be setting these options for a realistic flux
1247  LOG("gevgen_t2k", pFATAL)
1248  << "Using pre-calculated flux interaction probabilities only makes "
1249  << "sense when using JNUBEAM flux option!";
1250  exit(1);
1251  }
1252 
1253  // flux file POT normalization
1254  // only relevant when using the JNUBEAM flux ntuples
1255  if( parser.OptionExists('p') ) {
1256  LOG("gevgen_t2k", pDEBUG) << "Reading flux file normalization";
1257  gOptFluxNorm = parser.ArgAsDouble('p');
1258  } else {
1259  LOG("gevgen_t2k", pDEBUG)
1260  << "Setting standard normalization for JNUBEAM flux ntuples";
1262  } //-p
1263 
1264  // number of times to cycle through the JNUBEAM flux ntuple contents
1265  if( parser.OptionExists('c') ) {
1266  LOG("gevgen_t2k", pDEBUG) << "Reading number of flux ntuple cycles";
1267  gOptFluxNCycles = parser.ArgAsInt('c');
1268  } else {
1269  LOG("gevgen_t2k", pDEBUG)
1270  << "Setting standard number of cycles for JNUBEAM flux ntuples";
1271  gOptFluxNCycles = -1;
1272  } //-c
1273 
1274  // limit on max number of events that can be generated
1275  if( parser.OptionExists('n') ) {
1276  LOG("gevgen_t2k", pDEBUG)
1277  << "Reading limit on number of events to generate";
1278  gOptNev = parser.ArgAsInt('n');
1279  } else {
1280  LOG("gevgen_t2k", pDEBUG)
1281  << "Will keep on generating events till the flux driver stops";
1282  gOptNev = -1;
1283  } //-n
1284 
1285  // exposure (in POT)
1286  bool uppc_e = parser.OptionExists('E');
1287  char pot_args = 'e';
1288  bool pot_exit = true;
1289  if(uppc_e) {
1290  pot_args = 'E';
1291  pot_exit = false;
1292  }
1293  gOptExitAtEndOfFullFluxCycles = pot_exit;
1294  if( parser.OptionExists(pot_args) ) {
1295  LOG("gevgen_t2k", pDEBUG) << "Reading requested exposure in POT";
1296  gOptPOT = parser.ArgAsDouble(pot_args);
1297  } else {
1298  LOG("gevgen_t2k", pDEBUG) << "No POT exposure was requested";
1299  gOptPOT = -1;
1300  } //-e, -E
1301 
1302  // event file prefix
1303  if( parser.OptionExists('o') ) {
1304  LOG("gevgen_t2k", pDEBUG) << "Reading the event filename prefix";
1305  gOptEvFilePrefix = parser.ArgAsString('o');
1306  } else {
1307  LOG("gevgen_t2k", pDEBUG)
1308  << "Will set the default event filename prefix";
1310  } //-o
1311 
1312 
1313  // random number seed
1314  if( parser.OptionExists("seed") ) {
1315  LOG("gevgen_t2k", pINFO) << "Reading random number seed";
1316  gOptRanSeed = parser.ArgAsLong("seed");
1317  } else {
1318  LOG("gevgen_t2k", pINFO) << "Unspecified random number seed - Using default";
1319  gOptRanSeed = -1;
1320  }
1321 
1322  // input cross-section file
1323  if( parser.OptionExists("cross-sections") ) {
1324  LOG("gevgen_t2k", pINFO) << "Reading cross-section file";
1325  gOptInpXSecFile = parser.ArgAsString("cross-sections");
1326  } else {
1327  LOG("gevgen_t2k", pINFO) << "Unspecified cross-section file";
1328  gOptInpXSecFile = "";
1329  }
1330 
1331  //
1332  // >>> perform 'sanity' checks on command line arguments
1333  //
1334 
1335  // If we use a JNUBEAM flux ntuple, the 'exposure' may be set
1336  // either as:
1337  // - a number of POTs (whichever number of flux ntuple cycles that corresponds to)
1338  // - a number of generated events (whichever number of POTs that corresponds to)
1339  // - a number of flux ntuple cycles (whichever number of POTs that corresponds to)
1340  // Only one of those options can be set.
1341  if(!gOptUsingHistFlux) {
1342  int nset=0;
1343  if(gOptPOT > 0) nset++;
1344  if(gOptFluxNCycles > 0) nset++;
1345  if(gOptNev > 0) nset++;
1346  if(nset==0) {
1347  LOG("gevgen_t2k", pWARN)
1348  << "** To use a JNUBEAM flux ntuple you need to specify an exposure, "
1349  << "either via the -c, -e or -n options";
1350  LOG("gevgen_t2k", pWARN)
1351  << "** gevgen_t2k automatically sets the exposure via '-c 1'";
1352  gOptFluxNCycles = 1;
1353  }
1354  if(nset>1) {
1355  LOG("gevgen_t2k", pFATAL)
1356  << "You can not specify more than one of the -c, -e or -n options";
1357  PrintSyntax();
1358  exit(1);
1359  }
1360  }
1361  // If we use a flux histograms (not JNUBEAM flux ntuples) then -currently- the
1362  // only way to control exposure is via a number of events
1363  if(gOptUsingHistFlux) {
1364  if(gOptNev < 0) {
1365  LOG("gevgen_t2k", pFATAL)
1366  << "If you're using flux from histograms you need to specify the -n option";
1367  PrintSyntax();
1368  exit(1);
1369  }
1370  }
1371  // If we don't use a detailed ROOT detector geometry (but just a target mix) then
1372  // don't accept POT as a way to control job statistics (not enough info is passed
1373  // in the target mix to compute POT & the calculation can be easily done offline)
1374  if(!gOptUsingRootGeom) {
1375  if(gOptPOT > 0) {
1376  LOG("gevgen_t2k", pFATAL)
1377  << "You may not use the -e, -E options "
1378  << "without a detailed detector geometry description input";
1379  exit(1);
1380  }
1381  }
1382 
1383  //
1384  // >>> print the command line options
1385  //
1386  PDGLibrary * pdglib = PDGLibrary::Instance();
1387 
1388  ostringstream gminfo;
1389  if (gOptUsingRootGeom) {
1390  gminfo << "Using ROOT geometry - file: " << gOptRootGeom
1391  << ", top volume: "
1392  << ((gOptRootGeomTopVol.size()==0) ? "<master volume>" : gOptRootGeomTopVol)
1393  << ", max{PL} file: "
1394  << ((gOptExtMaxPlXml.size()==0) ? "<none>" : gOptExtMaxPlXml)
1395  << ", length units: " << lunits
1396  << ", density units: " << dunits;
1397  } else {
1398  gminfo << "Using target mix - ";
1399  map<int,double>::const_iterator iter;
1400  for(iter = gOptTgtMix.begin(); iter != gOptTgtMix.end(); ++iter) {
1401  int pdg_code = iter->first;
1402  double wgt = iter->second;
1403  TParticlePDG * p = pdglib->Find(pdg_code);
1404  if(p) {
1405  string name = p->GetName();
1406  gminfo << "(" << name << ") -> " << 100*wgt << "% / ";
1407  }//p?
1408  }
1409  }
1410 
1411  ostringstream fluxinfo;
1412  if(gOptUsingHistFlux) {
1413  fluxinfo << "Using flux histograms - ";
1414  map<int,TH1D*>::const_iterator iter;
1415  for(iter = gOptFluxHst.begin(); iter != gOptFluxHst.end(); ++iter) {
1416  int pdg_code = iter->first;
1417  TH1D * spectrum = iter->second;
1418  TParticlePDG * p = pdglib->Find(pdg_code);
1419  if(p) {
1420  string name = p->GetName();
1421  fluxinfo << "(" << name << ") -> " << spectrum->GetName() << " / ";
1422  }//p?
1423  }
1424  } else {
1425  fluxinfo << "Using JNUBEAM flux ntuple - "
1426  << "file: " << gOptFluxFile
1427  << ", location: " << gOptDetectorLocation
1428  << ", pot norm: " << gOptFluxNorm;
1429 
1430  if( gOptFluxNtpNuList.size() > 0 ) {
1431  fluxinfo << ", ** this job is considering only: ";
1432  for(unsigned int inu = 0; inu < gOptFluxNtpNuList.size(); inu++) {
1433  fluxinfo << gOptFluxNtpNuList[inu];
1434  if(inu < gOptFluxNtpNuList.size()-1) fluxinfo << ",";
1435  }
1436  }
1437 
1438  }
1439 
1440  ostringstream exposure;
1441  if(gOptPOT > 0)
1442  exposure << "Number of POTs = " << gOptPOT;
1443  if(gOptFluxNCycles > 0)
1444  exposure << "Number of flux cycles = " << gOptFluxNCycles;
1445  if(gOptNev > 0)
1446  exposure << "Number of events = " << gOptNev;
1447 
1448  LOG("gevgen_t2k", pNOTICE)
1449  << "\n\n"
1450  << utils::print::PrintFramedMesg("T2K event generation job configuration");
1451 
1452  LOG("gevgen_t2k", pNOTICE)
1453  << "\n - Run number: " << gOptRunNu
1454  << "\n - Random number seed: " << gOptRanSeed
1455  << "\n - Using cross-section file: " << gOptInpXSecFile
1456  << "\n - Flux @ " << fluxinfo.str()
1457  << "\n - Geometry @ " << gminfo.str()
1458  << "\n - Exposure @ " << exposure.str();
1459 
1460  LOG("gevgen_t2k", pNOTICE) << *RunOpt::Instance();
1461 }
map< int, TH1D * > gOptFluxHst
Definition: gT2KEvGen.cxx:488
string gOptRootGeomTopVol
Definition: gT2KEvGen.cxx:490
string gOptRootGeom
Definition: gT2KEvGen.cxx:489
const XML_Char * name
Definition: expat.h:151
double gOptFluxNorm
Definition: gT2KEvGen.cxx:496
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
string kDefOptGeomDUnits
Definition: gT2KEvGen.cxx:477
map< int, double > gOptTgtMix
Definition: gT2KEvGen.cxx:487
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
bool gOptUsingHistFlux
Definition: gT2KEvGen.cxx:486
string gOptDetectorLocation
Definition: gT2KEvGen.cxx:495
double gOptPOT
Definition: gT2KEvGen.cxx:500
Long_t gOptRunNu
Definition: gT2KEvGen.cxx:484
Loaders::FluxType flux
string gOptFluxProbFileName
Definition: gT2KEvGen.cxx:505
string gOptFluxFile
Definition: gT2KEvGen.cxx:494
int gOptNev
Definition: gT2KEvGen.cxx:499
string gOptInpXSecFile
Definition: gT2KEvGen.cxx:509
bool gOptRandomFluxOffset
Definition: gT2KEvGen.cxx:507
unsigned int nutype(unsigned int u)
Definition: runWimpSim.h:122
string kDefOptGeomLUnits
Definition: gT2KEvGen.cxx:476
string gOptExtMaxPlXml
Definition: gT2KEvGen.cxx:493
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Definition: typedefs.hpp:11
void PrintSyntax(void)
Definition: gT2KEvGen.cxx:1463
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
bool gOptUseFluxProbs
Definition: gT2KEvGen.cxx:503
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
double UnitFromString(string u)
Definition: UnitUtils.cxx:26
string kDefOptEvFilePrefix
Definition: gT2KEvGen.cxx:480
double kDefOptFluxNorm
Definition: gT2KEvGen.cxx:479
TH2D * histo
#define pINFO
Definition: Messenger.h:63
const ana::Var wgt
bool gOptUsingRootGeom
Definition: gT2KEvGen.cxx:485
#define pWARN
Definition: Messenger.h:61
PDGCodeList gOptFluxNtpNuList(false)
string gOptSaveFluxProbsFileName
Definition: gT2KEvGen.cxx:506
long int gOptRanSeed
Definition: gT2KEvGen.cxx:508
double gOptGeomDUnits
Definition: gT2KEvGen.cxx:492
bool gOptSaveFluxProbsFile
Definition: gT2KEvGen.cxx:504
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:42
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
string gOptEvFilePrefix
Definition: gT2KEvGen.cxx:502
int gOptFluxNCycles
Definition: gT2KEvGen.cxx:498
double gOptGeomLUnits
Definition: gT2KEvGen.cxx:491
void geom(int which=0)
Definition: geom.C:163
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
exit(0)
string PrintFramedMesg(string mesg, unsigned int nl=1, const char f='*')
Definition: PrintUtils.cxx:171
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:62
bool gOptExitAtEndOfFullFluxCycles
Definition: gT2KEvGen.cxx:501
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
#define pDEBUG
Definition: Messenger.h:64
int main ( int  argc,
char **  argv 
)

Definition at line 512 of file gT2KEvGen.cxx.

References genie::flux::GCylindTH1Flux::AddEnergySpectrum(), genie::NtpWriter::AddEventRecord(), ana::assert(), genie::RunOpt::BuildTune(), genie::utils::app_init::CacheFile(), genie::GMCJDriver::Configure(), genie::NtpWriter::CustomizeFilenamePrefix(), genie::flux::GJPARCNuFlux::DisableOffset(), genie::flux::GJPARCNuFlux::End(), genie::NtpWriter::EventTree(), exit(), flux, genie::GMCJDriver::ForceSingleProbScale(), genie::GMCJDriver::GenerateEvent(), GetCommandLineArgs(), genie::geometry::ROOTGeomAnalyzer::GetGeometry(), genie::GMCJDriver::GlobProbScale(), gOptDetectorLocation, gOptEvFilePrefix, gOptExitAtEndOfFullFluxCycles, gOptExtMaxPlXml, gOptFluxFile, gOptFluxHst, gOptFluxNCycles, gOptFluxNorm, gOptFluxNtpNuList, gOptFluxProbFileName, gOptGeomDUnits, gOptGeomLUnits, gOptInpXSecFile, gOptNev, gOptPOT, gOptRandomFluxOffset, gOptRanSeed, gOptRootGeom, gOptRootGeomTopVol, gOptRunNu, gOptSaveFluxProbsFile, gOptSaveFluxProbsFileName, gOptTgtMix, gOptUseFluxProbs, gOptUsingHistFlux, gOptUsingRootGeom, genie::NtpWriter::Initialize(), genie::RunOpt::Instance(), makeTrainCVSamples::int, it, kDefOptNtpFormat, genie::geometry::ROOTGeomAnalyzer::LengthUnits(), genie::flux::GJPARCNuFlux::LoadBeamSimData(), genie::GMCJDriver::LoadFluxProbabilities(), LOG, genie::utils::app_init::MesgThresholds(), runNovaSAM::metadata, genie::flux::GJPARCNuFlux::PassThroughInfo(), pERROR, pFATAL, pINFO, pNOTICE, pot, genie::flux::GJPARCNuFlux::POT_1cycle(), genie::flux::GJPARCNuFlux::POT_curravg(), genie::GMCJDriver::PreCalcFluxProbabilities(), genie::utils::app_init::RandGen(), genie::utils::geometry::RecursiveExhaust(), genie::NtpWriter::Save(), genie::GMCJDriver::SaveFluxProbabilities(), genie::flux::GCylindTH1Flux::SetBeamSpot(), genie::GMCJDriver::SetEventGeneratorList(), genie::flux::GJPARCNuFlux::SetFilePOT(), genie::flux::GJPARCNuFlux::SetFluxParticles(), genie::flux::GCylindTH1Flux::SetNuDirection(), genie::flux::GJPARCNuFlux::SetNumOfCycles(), genie::GHepRecord::SetPrintLevel(), genie::GMCJMonitor::SetRefreshRate(), genie::flux::GCylindTH1Flux::SetTransverseRadius(), genie::flux::GJPARCNuFlux::SetUpstreamZ(), log::success(), genie::GMCJDriver::SumFluxIntProbs(), genie::GMCJMonitor::Update(), genie::GMCJDriver::UseFluxDriver(), genie::GMCJDriver::UseGeomAnalyzer(), genie::GMCJDriver::UseMaxPathLengths(), genie::GMCJDriver::UseSplines(), genie::utils::app_init::XSecTable(), and make_true_q0q3_plots::zmax.

513 {
514  // Parse command line arguments
515  GetCommandLineArgs(argc,argv);
516 
517 
518  if ( ! RunOpt::Instance()->Tune() ) {
519  LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
520  exit(-1);
521  }
522  RunOpt::Instance()->BuildTune();
523 
524  // Initialization of random number generators, cross-section table,
525  // messenger thresholds, cache file
526  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
527  utils::app_init::CacheFile(RunOpt::Instance()->CacheFile());
530 
531  // Set GHEP print level
532  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
533 
534  // *************************************************************************
535  // * Create / configure the geometry driver
536  // *************************************************************************
537  GeomAnalyzerI * geom_driver = 0;
538  double zmin=0, zmax=0;
539 
540  if(gOptUsingRootGeom) {
541  //
542  // *** Using a realistic root-based detector geometry description
543  //
544 
545  // creating & configuring a root geometry driver
548  rgeom -> SetLengthUnits (gOptGeomLUnits);
549  rgeom -> SetDensityUnits (gOptGeomDUnits);
550  // getting the bounding box dimensions, before setting topvolume,
551  // along z so as to set the appropriate upstream generation surface
552  // for the JPARC flux driver
553  TGeoVolume * topvol = rgeom->GetGeometry()->GetTopVolume();
554  TGeoShape * bounding_box = topvol->GetShape();
555  bounding_box->GetAxisRange(3, zmin, zmax);
556  zmin *= rgeom->LengthUnits();
557  zmax *= rgeom->LengthUnits();
558  // now update to the requested topvolume for use in recursive exhaust method
559  rgeom -> SetTopVolName (gOptRootGeomTopVol);
560  topvol = rgeom->GetGeometry()->GetTopVolume();
561  if(!topvol) {
562  LOG("gevgen_t2k", pFATAL) << "Null top ROOT geometry volume!";
563  exit(1);
564  }
565  // switch on/off volumes as requested
566  if ( (gOptRootGeomTopVol[0] == '+') || (gOptRootGeomTopVol[0] == '-') ) {
567  bool exhaust = (*gOptRootGeomTopVol.c_str() == '+');
569  }
570 
571  // casting to the GENIE geometry driver interface
572  geom_driver = dynamic_cast<GeomAnalyzerI *> (rgeom);
573  }
574  else {
575  //
576  // *** Using a 'point' geometry with the specified target mix
577  // *** ( = a list of targets with their corresponding weight fraction)
578  //
579 
580  // creating & configuring a point geometry driver
583  // casting to the GENIE geometry driver interface
584  geom_driver = dynamic_cast<GeomAnalyzerI *> (pgeom);
585  }
586 
587  // *************************************************************************
588  // * Create / configure the flux driver
589  // *************************************************************************
590  GFluxI * flux_driver = 0;
591 
592  flux::GJPARCNuFlux * jparc_flux_driver = 0;
593  flux::GCylindTH1Flux * hst_flux_driver = 0;
594 
595  if(!gOptUsingHistFlux) {
596  //
597  // *** Using the detailed JPARC neutrino flux desription by feeding-in
598  // *** the JNUBEAM flux simulation ntuples
599  //
600 
601  // creating JPARC neutrino flux driver
602  jparc_flux_driver = new flux::GJPARCNuFlux;
603  // before loading the beam sim data set whether to use a random offset when looping
604  if(gOptRandomFluxOffset == false) jparc_flux_driver->DisableOffset();
605  // specify input JNUBEAM file & detector location
606  bool beam_sim_data_success = jparc_flux_driver->LoadBeamSimData(gOptFluxFile, gOptDetectorLocation);
607  if(!beam_sim_data_success) {
608  LOG("gevgen_t2k", pFATAL)
609  << "The flux driver has not been properly configured. Exiting";
610  exit(1);
611  }
612  // specify JNUBEAM normalization
613  jparc_flux_driver->SetFilePOT(gOptFluxNorm);
614  // specify upstream generation surface
615  jparc_flux_driver->SetUpstreamZ(zmin);
616  // specify which flavours to consider -
617  // if no neutrino list was specified then the MC job will consider all flavours
618  if( gOptFluxNtpNuList.size() > 0 ) {
619  jparc_flux_driver->SetFluxParticles(gOptFluxNtpNuList);
620  }
621 
622  // casting to the GENIE flux driver interface
623  flux_driver = dynamic_cast<GFluxI *> (jparc_flux_driver);
624  }
625  else {
626  //
627  // *** Using fluxes from histograms (for all specified neutrino species)
628  //
629 
630  // creating & configuring a generic GCylindTH1Flux flux driver
631  TVector3 bdir (0,0,1); // dir along +z
632  TVector3 bspot(0,0,0);
633  hst_flux_driver = new flux::GCylindTH1Flux;
634  hst_flux_driver->SetNuDirection (bdir);
635  hst_flux_driver->SetBeamSpot (bspot);
636  hst_flux_driver->SetTransverseRadius (-1);
637  map<int,TH1D*>::iterator it = gOptFluxHst.begin();
638  for( ; it != gOptFluxHst.end(); ++it) {
639  int pdg_code = it->first;
640  TH1D * spectrum = new TH1D(*(it->second));
641  hst_flux_driver->AddEnergySpectrum(pdg_code, spectrum);
642  }
643  // casting to the GENIE flux driver interface
644  flux_driver = dynamic_cast<GFluxI *> (hst_flux_driver);
645  }
646 
647  // *************************************************************************
648  // * Create/configure the event generation driver
649  // *************************************************************************
650  GMCJDriver * mcj_driver = new GMCJDriver;
651  mcj_driver->SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
652  mcj_driver->UseFluxDriver(flux_driver);
653  mcj_driver->UseGeomAnalyzer(geom_driver);
654  mcj_driver->UseMaxPathLengths(gOptExtMaxPlXml);
655  // do not calculate probability scales if using pre-generated flux probs
656  bool calc_prob_scales = (gOptSaveFluxProbsFile || gOptUseFluxProbs) ? false : true;
657  mcj_driver->Configure(calc_prob_scales);
658  mcj_driver->UseSplines();
659  mcj_driver->ForceSingleProbScale();
660 
661  // *************************************************************************
662  // * If specified use pre-calculated flux interaction probabilities instead
663  // * of preselecting based on max path lengths
664  // *************************************************************************
665 
667 
668  bool success = false;
669 
670  // set flux probs output file name
672  // default output name is ${FLUFILENAME}.${TOPVOL}.flxprobs.root
673  string basename = gOptFluxFile.substr(gOptFluxFile.rfind("/")+1);
674  string name = basename.substr(0, basename.rfind("."));
675  if(gOptRootGeomTopVol.length()>0)
676  name += "."+gOptRootGeomTopVol+".flxprobs.root";
677  else
678  name += ".master.flxprobs.root";
679  // if specified override with cmd line option
681  // Tell the driver save pre-generated probabilities to an output file
682  mcj_driver->SaveFluxProbabilities(name);
683  }
684 
685  // Either load pre-generated flux probabilities
686  if(gOptFluxProbFileName.size() > 0){
687  success = mcj_driver->LoadFluxProbabilities(gOptFluxProbFileName);
688  }
689  // Or pre-calculate them
690  else success = mcj_driver->PreCalcFluxProbabilities();
691 
692  if(success){
693  LOG("gevgen_t2k", pNOTICE)
694  << "Successfully calculated/loaded flux interaction probabilities!";
695  // Print out a list of expected number of events per POT and per cycle
696  // based on the pre-generated flux interaction probabilities
697  map<int, double> sum_probs_map = mcj_driver->SumFluxIntProbs();
698  map<int, double>::const_iterator sum_probs_it = sum_probs_map.begin();
699  double ntot_per_pot = 0.0;
700  double ntot_per_cycle = 0.0;
701  double pscale = mcj_driver->GlobProbScale();
702  double pot_1cycle = jparc_flux_driver->POT_1cycle();
703  LOG("T2KProdInfo", pNOTICE) <<
704  "Expected event rates based on flux interaction probabilities:";
705  for(; sum_probs_it != sum_probs_map.end(); sum_probs_it++){
706  double sum_probs = sum_probs_it->second;
707  double nevts_per_cycle = sum_probs / pscale; // take into account rescale
708  double nevts_per_pot = sum_probs/pot_1cycle; // == (sum_probs*pscale)/(pot_1cycle*pscale)
709  ntot_per_pot += nevts_per_pot;
710  ntot_per_cycle += nevts_per_cycle;
711  LOG("T2KProdInfo", pNOTICE) <<
712  " PDG "<< sum_probs_it->first << ": " << nevts_per_cycle <<
713  " Events/Cycle, "<< nevts_per_pot << " Events/POT";
714  }
715  LOG("T2KProdInfo", pNOTICE) << " -----------";
716  LOG("T2KProdInfo", pNOTICE) <<
717  " All neutrino species: " << ntot_per_cycle <<
718  " Events/Cycle, "<< ntot_per_pot << " Events/POT";
719  LOG("T2KProdInfo", pNOTICE) <<
720  "N.B. This assumes input flux file corresponds to "<< pot_1cycle <<
721  "POT, ensure this is correct if using these numbers!";
722  }
723  else {
724  LOG("gevgen_t2k", pFATAL)
725  << "Failed to calculated/loaded flux interaction probabilities!";
726  return 1;
727  }
728 
729  // Exit now if just pre-generating interaction probabilities
731  LOG("gevgen_t2k", pNOTICE)
732  << "Will not generate events - just pre-calculating flux interaction"
733  << "probabilities";
734  return 0;
735  }
736  } // Pre-calculated flux interaction probabilities
737 
738  // *************************************************************************
739  // * Work out number of cycles for current exposure settings
740  // *************************************************************************
741 
742  if(!gOptUsingHistFlux) {
743  // If a number of POT was requested, then work out how many flux ntuple
744  // cycles are required for accumulating those statistics
745  if(gOptPOT>0) {
746  double fpot_1c = jparc_flux_driver->POT_1cycle(); // flux POT / cycle
747  double psc = mcj_driver->GlobProbScale(); // interaction prob. scale
748  double pot_1c = fpot_1c / psc; // actual POT / cycle
749  int ncycles = (int) TMath::Max(1., TMath::Ceil(gOptPOT/pot_1c));
750 
751  LOG("gevgen_t2k", pNOTICE)
752  << " *** POT/cycle: " << pot_1c;
753  LOG("gevgen_t2k", pNOTICE)
754  << " *** Requested POT will be accumulated in: "
755  << ncycles << " flux ntuple cycles";
756 
757  jparc_flux_driver->SetNumOfCycles(ncycles);
758  }
759  // If a number of events was requested, then set the number of flux
760  // ntuple cycles to 'infinite'
761  else if(gOptNev>0) {
762  jparc_flux_driver->SetNumOfCycles(0);
763  }
764  // Just set the number of cycles to the requested value
765  else {
766  jparc_flux_driver->SetNumOfCycles(gOptFluxNCycles);
767  }
768  }
769 
770  // *************************************************************************
771  // * Prepare for writing the output event tree & status file
772  // *************************************************************************
773 
774  // Initialize an Ntuple Writer to save GHEP records into a TTree
776  ntpw.CustomizeFilenamePrefix(gOptEvFilePrefix);
777  ntpw.Initialize();
778 
779  // Add a custom-branch at the standard GENIE event tree so that
780  // info on the flux neutrino parent particle can be passed-through
781  flux::GJPARCNuFluxPassThroughInfo * flux_info = 0;
782  if(!gOptUsingHistFlux) {
783  TBranch * flux = ntpw.EventTree()->Branch("flux",
784  "genie::flux::GJPARCNuFluxPassThroughInfo", &flux_info, 32000, 1);
785  assert(flux);
786  flux->SetAutoDelete(kFALSE);
787  }
788 
789  // Create a MC job monitor for a periodically updated status file
790  GMCJMonitor mcjmonitor(gOptRunNu);
791  mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());
792 
793  // *************************************************************************
794  // * Event generation loop
795  // *************************************************************************
796 
797  int ievent = 0;
798  while (1)
799  {
800  LOG("gevgen_t2k", pNOTICE)
801  << " *** Generating event............ " << ievent;
802 
803  // In case the required statistics was expressed as 'number of events'
804  // then quit if that number has been generated
805  if(ievent == gOptNev) break;
806 
807  // In case the required statistics was expressed as 'number of POT' and
808  // the user does not want to wait till the end of the flux cycle to exit
809  // the event loop, then quit if the requested POT has been generated.
810  // In this case the computed POT may not be as accurate as if the program
811  // was waiting for the current flux cycle to be completed.
813  double fpot = jparc_flux_driver->POT_curravg(); // current POT in flux file
814  double psc = mcj_driver->GlobProbScale(); // interaction prob. scale
815  double pot = fpot / psc; // POT for generated sample
816  if(pot >= gOptPOT) break;
817  }
818 
819  // Generate a single event using neutrinos coming from the specified flux
820  // and hitting the specified geometry or target mix
821  EventRecord * event = mcj_driver->GenerateEvent();
822 
823  // Check whether a null event was returned due to the flux driver reaching
824  // the end of the input flux ntuple - exit the event generation loop
825  if(!event && jparc_flux_driver->End()) {
826  LOG("gevgen_t2k", pNOTICE)
827  << "** The JPARC flux driver read all the input flux ntuple entries";
828  break;
829  }
830  if(!event) {
831  LOG("gevgen_t2k", pERROR)
832  << "Got a null generated neutino event! Retrying ...";
833  continue;
834  }
835  LOG("gevgen_t2k", pINFO)
836  << "Generated event: " << *event;
837 
838  // A valid event was generated: extract flux info (parent decay/prod
839  // position/kinematics) for that simulated event so that it can be
840  // passed-through.
841  // Can only do so if I am generating events using the JNUBEAM flux
842  // ntuples, not simple histograms
843  if(!gOptUsingHistFlux) {
844  flux_info = new flux::GJPARCNuFluxPassThroughInfo(
845  jparc_flux_driver->PassThroughInfo());
846  LOG("gevgen_t2k", pINFO)
847  << "Pass-through flux info associated with generated event: "
848  << *flux_info;
849  }
850 
851  // Add event at the output ntuple, refresh the mc job monitor & clean-up
852  ntpw.AddEventRecord(ievent, event);
853  mcjmonitor.Update(ievent,event);
854  delete event;
855  if(flux_info) delete flux_info;
856  ievent++;
857  } //1
858 
859  LOG("gevgen_t2k", pNOTICE)
860  << "The GENIE MC job is done generaing events - Cleaning up & exiting...";
861 
862  // *************************************************************************
863  // * Print job statistics &
864  // * calculate normalization factor for the generated sample
865  // *************************************************************************
867  {
868  // POT normalization will only be calculated if event generation was based
869  // on beam simulation outputs (not just histograms) & a detailed detector
870  // geometry description.
871  double fpot = jparc_flux_driver->POT_curravg(); // current POT in flux file
872  double psc = mcj_driver->GlobProbScale(); // interaction prob. scale
873  double pot = fpot / psc; // POT for generated sample
874  // Get nunber of flux neutrinos read-in by flux friver, number of flux
875  // neutrinos actually thrown to the event generation driver and number
876  // of neutrino interactions actually generated
877  long int nflx_evg = mcj_driver -> NFluxNeutrinos();
878  long int nflx = jparc_flux_driver -> NFluxNeutrinos();
879  long int nev = ievent;
880 
881  LOG("gevgen_t2k", pNOTICE)
882  << "\n >> Actual JNUBEAM flux file normalization: " << fpot
883  << " POT * " << ((gOptDetectorLocation == "sk") ? "cm^2" : "det")
884  << "\n >> Interaction probability scaling factor: " << psc
885  << "\n >> N of flux v read-in by flux driver: " << nflx
886  << "\n >> N of flux v thrown to event gen driver: " << nflx_evg
887  << "\n >> N of generated v interactions: " << nev
888  << "\n ** Normalization for generated sample: " << pot
889  << " POT * " << ((gOptDetectorLocation == "sk") ? "cm^2" : "det");
890 
891  ntpw.EventTree()->SetWeight(pot); // POT
892  }
893 
894  // *************************************************************************
895  // * MC job meta-data
896  // *************************************************************************
897 
899 
900  metadata -> jnubeam_file = gOptFluxFile;
901  metadata -> detector_location = gOptDetectorLocation;
902  metadata -> geom_file = gOptRootGeom;
903  metadata -> geom_top_volume = gOptRootGeomTopVol;
904  metadata -> geom_length_units = gOptGeomLUnits;
905  metadata -> geom_density_units = gOptGeomDUnits;
906  metadata -> using_root_geom = gOptUsingRootGeom;
907  metadata -> target_mix = gOptTgtMix;
908  metadata -> using_hist_flux = gOptUsingHistFlux;
909  metadata -> flux_hists = gOptFluxHst;
910 
911  ntpw.EventTree()->GetUserInfo()->Add(metadata);
912 
913  // *************************************************************************
914  // * Save & clean-up
915  // *************************************************************************
916 
917  // Save the generated event tree & close the output file
918  ntpw.Save();
919 
920  // Clean-up
921  delete geom_driver;
922  delete flux_driver;
923  delete mcj_driver;
924  map<int,TH1D*>::iterator it = gOptFluxHst.begin();
925  for( ; it != gOptFluxHst.end(); ++it) {
926  TH1D * spectrum = it->second;
927  if(spectrum) delete spectrum;
928  }
929  gOptFluxHst.clear();
930 
931  LOG("gevgen_t2k", pNOTICE) << "Done!";
932 
933  return 0;
934 }
map< int, TH1D * > gOptFluxHst
Definition: gT2KEvGen.cxx:488
string gOptRootGeomTopVol
Definition: gT2KEvGen.cxx:490
string gOptRootGeom
Definition: gT2KEvGen.cxx:489
void RandGen(long int seed)
Definition: AppInit.cxx:31
const XML_Char * name
Definition: expat.h:151
double gOptFluxNorm
Definition: gT2KEvGen.cxx:496
double POT_curravg(void)
current average POT
void XSecTable(string inpfile, bool require_table)
Definition: AppInit.cxx:39
bool PreCalcFluxProbabilities(void)
Definition: GMCJDriver.cxx:239
map< int, double > SumFluxIntProbs(void) const
Definition: GMCJDriver.h:74
set< int >::iterator it
bool LoadBeamSimData(string filename, string det_loc)
load a jnubeam root flux ntuple
bool End(void)
true if no more flux nu&#39;s can be thrown (eg reaching end of beam sim ntuples)
Definition: GJPARCNuFlux.h:67
#define pERROR
Definition: Messenger.h:60
void SetUpstreamZ(double z0)
set flux neutrino initial z position (upstream of the detector)
void SetEventGeneratorList(string listname)
Definition: GMCJDriver.cxx:157
double POT_1cycle(void)
flux POT per cycle
map< int, double > gOptTgtMix
Definition: gT2KEvGen.cxx:487
#define pFATAL
Definition: Messenger.h:57
bool gOptUsingHistFlux
Definition: gT2KEvGen.cxx:486
string gOptDetectorLocation
Definition: gT2KEvGen.cxx:495
NtpMCFormat_t kDefOptNtpFormat
Definition: gT2KEvGen.cxx:478
double gOptPOT
Definition: gT2KEvGen.cxx:500
void SetFluxParticles(const PDGCodeList &particles)
specify list of flux neutrino species
Long_t gOptRunNu
Definition: gT2KEvGen.cxx:484
Loaders::FluxType flux
string gOptFluxProbFileName
Definition: gT2KEvGen.cxx:505
string gOptFluxFile
Definition: gT2KEvGen.cxx:494
int gOptNev
Definition: gT2KEvGen.cxx:499
string gOptInpXSecFile
Definition: gT2KEvGen.cxx:509
A GENIE flux driver encapsulating the JPARC neutrino flux. It reads-in the official JPARC neutrino fl...
Definition: GJPARCNuFlux.h:51
void UseFluxDriver(GFluxI *flux)
Definition: GMCJDriver.cxx:174
bool gOptRandomFluxOffset
Definition: gT2KEvGen.cxx:507
Simple class to create & update MC job status files and env. vars. This is used to be able to keep tr...
Definition: GMCJMonitor.h:30
A GENIE `MC Job Driver&#39;. Can be used for setting up complicated event generation cases involving deta...
Definition: GMCJDriver.h:47
double GlobProbScale(void) const
Definition: GMCJDriver.h:72
string gOptExtMaxPlXml
Definition: gT2KEvGen.cxx:493
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
void ForceSingleProbScale(void)
Definition: GMCJDriver.cxx:219
void SetFilePOT(double pot)
flux file norm is in /N POT/det [ND] or /N POT/cm^2 [FD]. Specify N (typically 1E+21) ...
virtual double LengthUnits(void) const
A generic GENIE flux driver. Generates a &#39;cylindrical&#39; neutrino beam along the input direction...
bool gOptUseFluxProbs
Definition: gT2KEvGen.cxx:503
bool UseMaxPathLengths(string xml_filename)
Definition: GMCJDriver.cxx:190
#define pot
void Configure(bool calc_prob_scales=true)
Definition: GMCJDriver.cxx:446
A ROOT/GEANT4 geometry driver.
def success(message)
Definition: log.py:5
#define pINFO
Definition: Messenger.h:63
void GetCommandLineArgs(int argc, char **argv)
Definition: gT2KEvGen.cxx:936
void SetTransverseRadius(double Rt)
void DisableOffset(void)
switch off random offset, must be called before LoadBeamSimData to have any effect ...
Definition: GJPARCNuFlux.h:84
void SetNuDirection(const TVector3 &direction)
bool gOptUsingRootGeom
Definition: gT2KEvGen.cxx:485
PDGCodeList gOptFluxNtpNuList(false)
string gOptSaveFluxProbsFileName
Definition: gT2KEvGen.cxx:506
EventRecord * GenerateEvent(void)
Definition: GMCJDriver.cxx:838
long int gOptRanSeed
Definition: gT2KEvGen.cxx:508
double gOptGeomDUnits
Definition: gT2KEvGen.cxx:492
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
bool gOptSaveFluxProbsFile
Definition: gT2KEvGen.cxx:504
A utility class to facilitate creating the GENIE MC Ntuple from the output GENIE GHEP event records...
Definition: NtpWriter.h:40
string gOptEvFilePrefix
Definition: gT2KEvGen.cxx:502
int gOptFluxNCycles
Definition: gT2KEvGen.cxx:498
double gOptGeomLUnits
Definition: gT2KEvGen.cxx:491
exit(0)
bool LoadFluxProbabilities(string filename)
Definition: GMCJDriver.cxx:390
void SetNumOfCycles(int n)
set how many times to cycle through the ntuple (default: 1 / n=0 means &#39;infinite&#39;) ...
assert(nhit_max >=nhit_nbins)
A vector of EventGeneratorI objects.
void SetBeamSpot(const TVector3 &spot)
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:100
void AddEnergySpectrum(int nu_pdgc, TH1D *spectrum)
Utility class to store MC job meta-data.
#define pNOTICE
Definition: Messenger.h:62
bool gOptExitAtEndOfFullFluxCycles
Definition: gT2KEvGen.cxx:501
Defines the GENIE Geometry Analyzer Interface.
Definition: GeomAnalyzerI.h:30
void RecursiveExhaust(TGeoVolume *topvol, string volnames, bool exhaust)
Definition: GeoUtils.cxx:26
void UseGeomAnalyzer(GeomAnalyzerI *geom)
Definition: GMCJDriver.cxx:179
void UseSplines(bool useLogE=true)
Definition: GMCJDriver.cxx:184
The PointGeomAnalyzer class is the simplest implementation of the GeomAnalyserI interface and defines...
void CacheFile(string inpfile)
Definition: AppInit.cxx:118
void SaveFluxProbabilities(string outfilename)
Definition: GMCJDriver.cxx:437
const GJPARCNuFluxPassThroughInfo & PassThroughInfo(void)
GJPARCNuFluxPassThroughInfo.
Definition: GJPARCNuFlux.h:93
virtual TGeoManager * GetGeometry(void) const
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
void PrintSyntax ( void  )

Definition at line 1463 of file gT2KEvGen.cxx.

References LOG, and pFATAL.

Referenced by GetCommandLineArgs().

1464 {
1465  LOG("gevgen_t2k", pFATAL)
1466  << "\n **Syntax**"
1467  << "\n gevgen_t2k [-h] "
1468  << "\n [-r run#]"
1469  << "\n -f flux"
1470  << "\n -g geometry"
1471  << "\n [-p pot_normalization_of_flux_file]"
1472  << "\n [-t top_volume_name_at_geom]"
1473  << "\n [-P pre_gen_prob_file]"
1474  << "\n [-S] [output_name]"
1475  << "\n [-m max_path_lengths_xml_file]"
1476  << "\n [-L length_units_at_geom]"
1477  << "\n [-D density_units_at_geom]"
1478  << "\n [-n n_of_events]"
1479  << "\n [-c flux_cycles]"
1480  << "\n [-e, -E exposure_in_POTs]"
1481  << "\n [-o output_event_file_prefix]"
1482  << "\n [-R]"
1483  << "\n [--seed random_number_seed]"
1484  << "\n --cross-sections xml_file"
1485  << "\n [--event-generator-list list_name]"
1486  << "\n [--message-thresholds xml_file]"
1487  << "\n [--unphysical-event-mask mask]"
1488  << "\n [--event-record-print-level level]"
1489  << "\n [--mc-job-status-refresh-rate rate]"
1490  << "\n [--cache-file root_file]"
1491  << "\n"
1492  << " Please also read the detailed documentation at http://www.genie-mc.org"
1493  << " or look at the source code: $GENIE/src/Apps/gT2KEvGen.cxx"
1494  << "\n";
1495 }
#define pFATAL
Definition: Messenger.h:57
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97

Variable Documentation

string gOptDetectorLocation

Definition at line 495 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

string gOptEvFilePrefix

Definition at line 502 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

bool gOptExitAtEndOfFullFluxCycles

Definition at line 501 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

string gOptExtMaxPlXml

Definition at line 493 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

string gOptFluxFile

Definition at line 494 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

map<int,TH1D*> gOptFluxHst

Definition at line 488 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

int gOptFluxNCycles

Definition at line 498 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

double gOptFluxNorm

Definition at line 496 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

PDGCodeList gOptFluxNtpNuList(false)

Referenced by GetCommandLineArgs(), and main().

string gOptFluxProbFileName

Definition at line 505 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

double gOptGeomDUnits = 0

Definition at line 492 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

double gOptGeomLUnits = 0

Definition at line 491 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

string gOptInpXSecFile

Definition at line 509 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

int gOptNev

Definition at line 499 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

double gOptPOT

Definition at line 500 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

bool gOptRandomFluxOffset = false

Definition at line 507 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

long int gOptRanSeed

Definition at line 508 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

string gOptRootGeom

Definition at line 489 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

string gOptRootGeomTopVol = ""

Definition at line 490 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

Long_t gOptRunNu

Definition at line 484 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

bool gOptSaveFluxProbsFile = false

Definition at line 504 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

string gOptSaveFluxProbsFileName

Definition at line 506 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

map<int,double> gOptTgtMix

Definition at line 487 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

bool gOptUseFluxProbs = false

Definition at line 503 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

bool gOptUsingHistFlux = false

Definition at line 486 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

bool gOptUsingRootGeom = false

Definition at line 485 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

string kDefOptEvFilePrefix = "gntp"

Definition at line 480 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs().

double kDefOptFluxNorm = 1E+21

Definition at line 479 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs().

string kDefOptGeomDUnits = "g_cm3"

Definition at line 477 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs().

string kDefOptGeomLUnits = "mm"

Definition at line 476 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs().

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 478 of file gT2KEvGen.cxx.

Referenced by main().