GENIEGenerator
Loading...
Searching...
No Matches
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"
Include dependency graph for gT2KEvGen.cxx:

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

◆ GetCommandLineArgs()

void GetCommandLineArgs ( int argc,
char ** argv )

Definition at line 948 of file gT2KEvGen.cxx.

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

References genie::CmdLnArgParser::ArgAsDouble(), genie::CmdLnArgParser::ArgAsInt(), genie::CmdLnArgParser::ArgAsLong(), genie::CmdLnArgParser::ArgAsString(), genie::RunOpt::EnableBareXSecPreCalc(), genie::PDGLibrary::Find(), 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, genie::PDGLibrary::Instance(), genie::RunOpt::Instance(), genie::pdg::IsAntiNeutrino(), genie::pdg::IsNeutrino(), kDefOptEvFilePrefix, kDefOptFluxNorm, kDefOptGeomDUnits, kDefOptGeomLUnits, LOG, lunits, genie::CmdLnArgParser::OptionExists(), pDEBUG, pFATAL, pINFO, pNOTICE, genie::utils::print::PrintFramedMesg(), PrintSyntax(), pWARN, genie::RunOpt::ReadFromCommandLine(), genie::utils::str::Split(), and genie::utils::units::UnitFromString().

Referenced by main().

◆ main()

int main ( int argc,
char ** argv )

Definition at line 524 of file gT2KEvGen.cxx.

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

References genie::flux::GCylindTH1Flux::AddEnergySpectrum(), genie::NtpWriter::AddEventRecord(), 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(), 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(), kDefOptNtpFormat, genie::geometry::ROOTGeomAnalyzer::LengthUnits(), genie::flux::GJPARCNuFlux::LoadBeamSimData(), genie::GMCJDriver::LoadFluxProbabilities(), LOG, genie::utils::app_init::MesgThresholds(), genie::flux::GJPARCNuFlux::PassThroughInfo(), pERROR, pFATAL, pINFO, pNOTICE, 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(), genie::GMCJDriver::SumFluxIntProbs(), genie::GMCJMonitor::Update(), genie::GMCJDriver::UseFluxDriver(), genie::GMCJDriver::UseGeomAnalyzer(), genie::GMCJDriver::UseMaxPathLengths(), genie::GMCJDriver::UseSplines(), and genie::utils::app_init::XSecTable().

◆ PrintSyntax()

void PrintSyntax ( void )

Definition at line 1502 of file gT2KEvGen.cxx.

1503{
1504 LOG("gevgen_t2k", pFATAL)
1505 << "\n **Syntax**"
1506 << "\n gevgen_t2k [-h] "
1507 << "\n [-r run#]"
1508 << "\n -f flux"
1509 << "\n -g geometry"
1510 << "\n [-p pot_normalization_of_flux_file]"
1511 << "\n [-t top_volume_name_at_geom]"
1512 << "\n [-P pre_gen_prob_file]"
1513 << "\n [-S] [output_name]"
1514 << "\n [-m max_path_lengths_xml_file]"
1515 << "\n [-L length_units_at_geom]"
1516 << "\n [-D density_units_at_geom]"
1517 << "\n [-n n_of_events]"
1518 << "\n [-c flux_cycles]"
1519 << "\n [-e, -E exposure_in_POTs]"
1520 << "\n [-o output_event_file_prefix]"
1521 << "\n [-R]"
1522 << "\n [--seed random_number_seed]"
1523 << "\n --cross-sections xml_file"
1524 << "\n [--event-generator-list list_name]"
1525 << "\n [--message-thresholds xml_file]"
1526 << "\n [--unphysical-event-mask mask]"
1527 << "\n [--event-record-print-level level]"
1528 << "\n [--mc-job-status-refresh-rate rate]"
1529 << "\n [--cache-file root_file]"
1530 << "\n"
1531 << " Please also read the detailed documentation at http://www.genie-mc.org"
1532 << " or look at the source code: $GENIE/src/Apps/gT2KEvGen.cxx"
1533 << "\n";
1534}

References LOG, and pFATAL.

Referenced by GetCommandLineArgs().

Variable Documentation

◆ gOptDetectorLocation

string gOptDetectorLocation

Definition at line 507 of file gT2KEvGen.cxx.

◆ gOptEvFilePrefix

string gOptEvFilePrefix

Definition at line 514 of file gT2KEvGen.cxx.

◆ gOptExitAtEndOfFullFluxCycles

bool gOptExitAtEndOfFullFluxCycles

Definition at line 513 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

◆ gOptExtMaxPlXml

string gOptExtMaxPlXml

Definition at line 505 of file gT2KEvGen.cxx.

◆ gOptFluxFile

string gOptFluxFile

Definition at line 506 of file gT2KEvGen.cxx.

◆ gOptFluxHst

map<int,TH1D*> gOptFluxHst

Definition at line 500 of file gT2KEvGen.cxx.

◆ gOptFluxNCycles

int gOptFluxNCycles

Definition at line 510 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

◆ gOptFluxNorm

double gOptFluxNorm

Definition at line 508 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

◆ gOptFluxNtpNuList

PDGCodeList gOptFluxNtpNuList(false) ( false )

Referenced by GetCommandLineArgs(), and main().

◆ gOptFluxProbFileName

string gOptFluxProbFileName

Definition at line 517 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

◆ gOptGeomDUnits

double gOptGeomDUnits = 0

Definition at line 504 of file gT2KEvGen.cxx.

◆ gOptGeomLUnits

double gOptGeomLUnits = 0

Definition at line 503 of file gT2KEvGen.cxx.

◆ gOptInpXSecFile

string gOptInpXSecFile

Definition at line 521 of file gT2KEvGen.cxx.

◆ gOptNev

int gOptNev

Definition at line 511 of file gT2KEvGen.cxx.

◆ gOptPOT

double gOptPOT

Definition at line 512 of file gT2KEvGen.cxx.

◆ gOptRandomFluxOffset

bool gOptRandomFluxOffset = false

Definition at line 519 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

◆ gOptRanSeed

long int gOptRanSeed

Definition at line 520 of file gT2KEvGen.cxx.

◆ gOptRootGeom

string gOptRootGeom

Definition at line 501 of file gT2KEvGen.cxx.

◆ gOptRootGeomTopVol

string gOptRootGeomTopVol = ""

Definition at line 502 of file gT2KEvGen.cxx.

◆ gOptRunNu

Long_t gOptRunNu

Definition at line 496 of file gT2KEvGen.cxx.

◆ gOptSaveFluxProbsFile

bool gOptSaveFluxProbsFile = false

Definition at line 516 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

◆ gOptSaveFluxProbsFileName

string gOptSaveFluxProbsFileName

Definition at line 518 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

◆ gOptTgtMix

map<int,double> gOptTgtMix

Definition at line 499 of file gT2KEvGen.cxx.

◆ gOptUseFluxProbs

bool gOptUseFluxProbs = false

Definition at line 515 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs(), and main().

◆ gOptUsingHistFlux

bool gOptUsingHistFlux = false

Definition at line 498 of file gT2KEvGen.cxx.

◆ gOptUsingRootGeom

bool gOptUsingRootGeom = false

Definition at line 497 of file gT2KEvGen.cxx.

◆ kDefOptEvFilePrefix

string kDefOptEvFilePrefix = "gntp"

Definition at line 492 of file gT2KEvGen.cxx.

◆ kDefOptFluxNorm

double kDefOptFluxNorm = 1E+21

Definition at line 491 of file gT2KEvGen.cxx.

Referenced by GetCommandLineArgs().

◆ kDefOptGeomDUnits

string kDefOptGeomDUnits = "g_cm3"

Definition at line 489 of file gT2KEvGen.cxx.

◆ kDefOptGeomLUnits

string kDefOptGeomLUnits = "mm"

Definition at line 488 of file gT2KEvGen.cxx.

◆ kDefOptNtpFormat

NtpMCFormat_t kDefOptNtpFormat = kNFGHEP

Definition at line 490 of file gT2KEvGen.cxx.