gFNALExptEvGen.cxx
Go to the documentation of this file.
1 //______________________________________________________________________________
2 /*!
3 
4 \program gevgen_fnal
5 
6 \brief A GENIE event generation driver `customized' for the FNAL neutrino
7  experiments using flux drivers for file types used by those expts.
8  This program was adapted from gevgen_t2k used at T2K.
9 
10  This driver can use either the FNAL-supported neutrino flux ntuples
11  (generated by gNuMI, g4numi, g4numi_flugg), or "dk2nu" format,
12  or plain flux histograms for all neutrino species you are considering.
13  You can specify either a ROOT-based detailed detector geometry
14  description or a simple target mix. See below for further details.
15 
16  Users should note that the generic GENIE event generation driver
17  (gevgen) may still be a more appropriate tool to use for the simpler
18  event generation casesrequired for many 4-vector level / systematic
19  studies.
20  Please see the GENIE documentation (http://www.genie-mc.org) and
21  contact me <costas.andreopoulos \at stfc.ac.uk> if in doubt.
22 
23  *** Synopsis :
24 
25  gevgen_fnal [-h]
26  [-r run#]
27  -f flux
28  -g geometry
29  [-t top_volume_name_at_geom || -t +Vol1-Vol2...]
30  [-m max_path_lengths_xml_file]
31  [-L length_units_at_geom]
32  [-D density_units_at_geom]
33  [-n n_of_events]
34  [-e exposure_in_POTs]
35  [-o output_event_file_prefix]
36  [-F fid_cut_string]
37  [-S nrays]
38  [-z zmin]
39  [-d debug flags]
40  [--seed random_number_seed]
41  --cross-sections xml_file
42  [--event-generator-list list_name]
43  [--tune genie_tune]
44  [--message-thresholds xml_file]
45  [--unphysical-event-mask mask]
46  [--event-record-print-level level]
47  [--mc-job-status-refresh-rate rate]
48  [--cache-file root_file]
49 
50  *** Options :
51 
52  [] Denotes an optional argument
53 
54  -h
55  Prints out the gevgen_fnal syntax and exits.
56  -r
57  Specifies the MC run number [default: 1000]
58  -g
59  Input 'geometry'.
60  This option can be used to specify any of:
61  1 > A ROOT file containing a ROOT/GEANT geometry description
62  [Examples]
63  - To use the master volume from the ROOT geometry stored
64  in the /some/path/nova-geom.root file, type:
65  '-g /some/path/nova-geom.root'
66  2 > A mix of target materials, each with its corresponding weight,
67  typed as a comma-separated list of nuclear PDG codes (in the
68  std PDG2006 convention: 10LZZZAAAI) with the weight fractions
69  in brackets, eg code1[fraction1],code2[fraction2],...
70  If that option is used (no detailed input geometry description)
71  then the interaction vertices are distributed in the detector
72  by the detector MC.
73  [Examples]
74  - To use a target mix of 95% O16 and 5% H type:
75  '-g 1000080160[0.95],1000010010[0.05]'
76  - To use a target which is 100% C12, type:
77  '-g 1000060120'
78  -t
79  Input 'top volume' for event generation -
80  can be used to force event generation in given sub-detector.
81  [default: the 'master volume' of the input geometry]
82  You can also use the -t option to switch generation on/off at
83  multiple volumes as, for example, in:
84  `-t +Vol1-Vol2+Vol3-Vol4',
85  `-t "+Vol1 -Vol2 +Vol3 -Vol4"',
86  `-t -Vol2-Vol4+Vol1+Vol3',
87  `-t "-Vol2 -Vol4 +Vol1 +Vol3"'m
88  where:
89  "+Vol1" and "+Vol3" tells GENIE to `switch on' Vol1 and Vol3, while
90  "-Vol2" and "-Vol4" tells GENIE to `switch off' Vol2 and Vol4.
91  If the very first character is a '+', GENIE will neglect all volumes
92  except the ones explicitly turned on. Vice versa, if the very first
93  character is a `-', GENIE will keep all volumes except the ones
94  explicitly turned off (feature contributed by J.Holeczek).
95  -m
96  An XML file (generated by gmxpl) with the max (density weighted)
97  path-lengths for each target material in the input ROOT geometry.
98  If no file is input, then the geometry will be scanned at MC job
99  initialization to determine those max path lengths.
100  Supplying this file can speed-up the MC job initialization.
101  -L
102  Input geometry length units, eg "m", "cm", "mm", ...
103  [default: mm]
104  -D
105  Input geometry density units, eg "g_cm3", "clhep_def_density_unit",...
106  [default: g_cm3]
107  -f
108  Input 'neutrino flux'.
109  This option can be used to specify any of:
110  1 > A g[4][numi|lbne][_flugg] beam simulation output file
111  and the detector location
112  The general sytax is:
113  -f /full/path/flux_file.root,detector,flavor1,flavor2...
114  [Notes]
115  - For more information on the flux ntuples, see the gNuMI doc.
116  - The original hbook ntuples need to be converted to a ROOT
117  format using the h2root ROOT utility.
118  - If flavors aren't specified then use default (12,-12,14,-14)
119  - See GNuMIFlux.xml or Dk2Nu.xml for all supported detector
120  locations
121  - The g[4][NuMI|LBNE][_flugg] flux ntuples are read via GENIE's
122  GNuMIFlux driver, and dk2nu file via an external
123  product w/ the GDk2NuFlux driver (if it can be found).
124  This customized GENIE event generation driver passes-through
125  the complete gNuMI input flux information (eg parent decay
126  kinematics / position etc) for each neutrino event it
127  generates (as an additional 'flux' branch is added on the
128  output event tree).
129  [Examples]
130  - To use the gNuMI flux ntuple flux.root at MINOS near
131  detector location, type:
132  '-f /path/flux.root,MINOS-NearDet'
133  1a> Similar to 1 above, but filename contains "dk2nu", then use
134  the GDk2NuFlux driver
135  1b> Similar to 1 above, but filename contains "gsimple", then
136  use GSimpleNtpFlux driver
137  2 > A set of histograms stored in a ROOT file.
138  The general syntax is:
139  -f /path/histogram_file.root,neutrino_code[histo_name],...
140  [Notes]
141  - The neutrino codes are the PDG ones.
142  - The 'neutrino_code[histogram_name]' part of the option can
143  be repeated multiple times (separated by commas), once for
144  each flux neutrino species you want to consider, eg
145  '-f somefile.root,12[nuehst],-12[nuebarhst],14[numuhst]'
146  - When using flux from histograms then there is no point in
147  using a 'detailed detector geometry description' as your
148  flux input contains no directional information for those
149  flux neutrinos.
150  The neutrino direction is conventionally set to be
151  +z {x=0,y=0}.
152  So, when using this option you must be using a simple
153  'target mix'
154  See the -g option for possible geometry settings.
155  If you want to use the detailed detector geometry
156  description then you should be using a driver that
157  supplies a detailed simulated beam flux.
158  - When using flux from histograms there is no branch with
159  neutrino parent information added in the output event
160  tree as your flux input contains no such information.
161  - Note that the relative normalization of the flux histograms
162  is taken into account and is reflected in the relative
163  frequency of flux neutrinos thrown by the flux driver
164  [Examples]
165  - To use the histogram 'h100' (representing the nu_mu flux)
166  and the histogram 'h300' (representing the nu_e flux) and
167  the histogram 'h301' (representing the nu_e_bar flux) from
168  the flux.root file in /path/ type:
169  '-f /path/flux.root,14[h100],12[h300],-12[h301]
170 
171  -e
172  Specifies how many POTs to generate.
173  -n
174  Specifies how many events to generate.
175 
176  -------
177  [Note on exposure / statistics]
178  Both -e and -n options can be used to set the exposure.
179  - If the input flux is a non-histogram driver then any of these
180  options can be used (one at a time).
181  - If the input flux is described with histograms then only the -n
182  option is available.
183  -------
184 
185  -F
186  Apply a fiducial cut (for now hard coded ... generalize)
187  Only used with ROOTGeomAnalyzer
188  if string starts with "-" then reverses sense (ie. anti-fiducial)
189  -S
190  Number of rays to use to scan geometry for max path length
191  Only used with ROOTGeomAnalyzer & { GNuMIFlux, GSimpleNtpFlux, GDk2NuFlux }
192  +N Use flux to scan geometry for max path length
193  -N Use N rays x N points on each face of a box
194  -z
195  Z from which to start flux ray in user world coordinates
196  Only use with ROOTGeomAnalyzer & { GNuMIFlux, GSimpleNtpFlux, GDk2NuFlux }
197  If left unset then flux originates on the flux window
198  [No longer attempts to determine z from geometry, generally
199  got this wrong]
200  -o
201  Sets the prefix of the output event file.
202  The output filename is built as:
203  [prefix].[run_number].[event_tree_format].[file_format]
204  The default output filename is:
205  gntp.[run_number].ghep.root
206  This cmd line arguments lets you override 'gntp'
207  --seed
208  Random number seed.
209  --cross-sections
210  Name (incl. full path) of an XML file with pre-computed
211  cross-section values used for constructing splines.
212  --tune
213  Specifies a GENIE comprehensive neutrino interaction model tune.
214  [default: "Default"].
215  --message-thresholds
216  Allows users to customize the message stream thresholds.
217  The thresholds are specified using an XML file.
218  See $GENIE/config/Messenger.xml for the XML schema.
219  --unphysical-event-mask
220  Allows users to specify a 16-bit mask to allow certain types of
221  unphysical events to be written in the output file.
222  [default: all unphysical events are rejected]
223  --event-record-print-level
224  Allows users to set the level of information shown when the event
225  record is printed in the screen. See GHepRecord::Print().
226  --mc-job-status-refresh-rate
227  Allows users to customize the refresh rate of the status file.
228  --cache-file
229  Allows users to specify a cache file so that the cache can be
230  re-used in subsequent MC jobs.
231 
232  *** Examples:
233 
234  (1) shell% gevgen_fnal
235  -r 1001
236  -f /data/mc_inputs/flux/flux_00001.root,MINOS-NearDet,12,-12
237  -g /data/mc_inputs/geom/minos.root
238  -L mm -D g_cm3
239  -e 5E+17
240  --cross-sections /data/xsec.xml
241 
242  Generate events (run number 1001) using the flux ntuple in
243  /data/mc_inputs/flux/v1/flux_00001.root
244  The job will load the MINOS near detector detector geometry
245  description from /data/mc_inputs/geom/minos.root and interpret it
246  using 'mm' as the length unit and 'g/cm^3' as the density unit.
247  The job will stop as it accumulates a sample corresponding to
248  5E+17 POT.
249  Pre-computed cross-section data are loaded from /data/xsec.xml
250 
251  (2) shell% gevgen_fnal
252  -r 1001
253  -f /data/t2k/flux/hst/flux.root,12[h100],-12[h101],14[h200]
254  -g 1000080160[0.95],1000010010[0.05]
255  -n 50000
256  --cross-sections /data/xsec.xml
257 
258  Please read the GENIE user manual for more information.
259 
260 \author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
261  University of Liverpool & STFC Rutherford Appleton Lab
262 
263  Robert Hatcher <rhatcher \at fnal.gov>
264  Fermi National Accelerator Laboratory
265 
266 \created August 20, 2008
267 
268 \cpright Copyright (c) 2003-2019, The GENIE Collaboration
269  For the full text of the license visit http://copyright.genie-mc.org
270  or see $GENIE/LICENSE
271 */
272 //_________________________________________________________________________________________
273 
274 #include <cassert>
275 #include <cstdlib>
276 #include <csignal>
277 
278 #include <string>
279 #include <sstream>
280 #include <vector>
281 #include <map>
282 #include <algorithm> // for transform()
283 #include <fstream>
284 
285 #include <TSystem.h>
286 #include <TError.h> // for gErrorIgnoreLevel
287 #include <TTree.h>
288 #include <TFile.h>
289 #include <TH1D.h>
290 #include <TMath.h>
291 #include <TGeoVolume.h>
292 #include <TGeoShape.h>
293 
309 #include "Framework/Utils/AppInit.h"
310 #include "Framework/Utils/RunOpt.h"
314 
315 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
320 
321 //#include "Tools/Flux/GNuMIFlux.h"
322 //#include "Tools/Flux/GSimpleNtpFlux.h"
323 // #ifdef __DK2NU_FLUX_DRIVER_AVAILABLE__
324 // #include "dk2nu/tree/dk2nu.h"
325 // #include "dk2nu/tree/dkmeta.h"
326 // #include "dk2nu/tree/NuChoice.h"
327 // #include "dk2nu/genie/GDk2NuFlux.h"
328 // #endif
329 
330 #endif
331 
332 #ifdef __GENIE_GEOM_DRIVERS_ENABLED__
333 #include "Tools/Geometry/GeoUtils.h"
338 #endif
339 
340 using std::string;
341 using std::vector;
342 using std::map;
343 using std::ostringstream;
344 
345 using namespace genie;
346 
347 // Forward function declarations
348 //
349 void LoadExtraOptions (void);
350 void GetCommandLineArgs (int argc, char ** argv);
351 void PrintSyntax (void);
352 void CreateFidSelection (string fidcut, GeomAnalyzerI* geom_driver);
353 void CreateRockBoxSelection (string fidcut, GeomAnalyzerI* geom_driver);
354 void DetermineFluxDriver(string fopt);
355 void ParseFluxHst (string fopt);
356 void ParseFluxFileConfig(string fopt);
357 
358 // Default options (override them using the command line arguments):
359 //
360 string kDefOptGeomLUnits = "mm"; // default geometry length units
361 string kDefOptGeomDUnits = "g_cm3"; // default geometry density units
362 NtpMCFormat_t kDefOptNtpFormat = kNFGHEP; // default event tree format
363 string kDefOptEvFilePrefix = "gntp";
364 
365 // User-specified options:
366 //
367 Long_t gOptRunNu; // run number
368 bool gOptUsingRootGeom = false; // using root geom or target mix?
369 bool gOptUsingHistFlux = false; // using beam flux ntuples or flux from histograms?
370 string gOptFluxDriver = ""; // flux driver class to use
371 map<string,string> gOptFluxShortNames;
372 PDGCodeList gOptFluxPdg; // list of neutrino flavors to accept
373 map<int,double> gOptTgtMix; // target mix (tgt pdg -> wght frac) / if not using detailed root geom
374 map<int,TH1D*> gOptFluxHst; // flux histos (nu pdg -> spectrum) / if not using beam sim ntuples
375 string gOptRootGeom; // input ROOT file with realistic detector geometry
376 string gOptRootGeomTopVol = ""; // input geometry top event generation volume
377 string gOptRootGeomMasterVol = ""; // (highest level of geometry)
378 double gOptGeomLUnits = 0; // input geometry length units
379 double gOptGeomDUnits = 0; // input geometry density units
380 string gOptExtMaxPlXml = ""; // max path lengths XML file for input geometry
381 bool gOptWriteMaxPlXml = false; // rather than read file, write the file
382  // triggered by leading '+' on given ofilename
383 string gOptFluxFile; // ROOT file with beam flux ntuple
384 string gOptDetectorLocation; // detector location (see GNuMIFlux.xml for supported locations))
385 int gOptNev; // number of events to generate
386 double gOptPOT; // exposure (in POT)
387 string gOptFidCut; // fiducial cut selection
388 int gOptNScan = 0; // # of geometry scan rays
389 double gOptZmin = -2.0e30; // starting z position [ if abs() < 1e30 ]
390 string gOptEvFilePrefix; // event file prefix
391 int gOptDebug = 0; // debug flags
392 long int gOptRanSeed; // random number seed
393 string gOptInpXSecFile; // cross-section splines
394 
395 bool gSigTERM = false; // was TERM signal sent?
396 
397 static void gsSIGTERMhandler(int /* s */)
398 {
399  gSigTERM = true;
400  std::cerr << "Caught SIGTERM" << std::endl;
401 }
402 
403 //____________________________________________________________________________
404 int main(int argc, char ** argv)
405 {
407  GetCommandLineArgs(argc,argv);
408 
409  if ( ! RunOpt::Instance()->Tune() ) {
410  LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
411  exit(-1);
412  }
414 
415  // Initialization of random number generators, cross-section table,
416  // messenger thresholds, cache file
417  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
421 
422  // Set GHEP print level
423  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
424 
425  // *************************************************************************
426  // * Create / configure the geometry driver
427  // *************************************************************************
428  GeomAnalyzerI * geom_driver = 0;
429 
430  if(gOptUsingRootGeom) {
431  //
432  // *** Using a realistic root-based detector geometry description
433  //
434 
435  // creating & configuring a root geometry driver
438  TGeoVolume * topvol = rgeom->GetGeometry()->GetTopVolume();
439  if ( ! topvol ) {
440  LOG("gevgen_fnal", pFATAL) << "Null top ROOT geometry volume!";
441  exit(1);
442  }
443  // retrieve the master volume name
444  gOptRootGeomMasterVol = topvol->GetName();
445 
446  rgeom -> SetLengthUnits (gOptGeomLUnits);
447  rgeom -> SetDensityUnits (gOptGeomDUnits);
448  rgeom -> SetTopVolName (gOptRootGeomTopVol); // set user defined "topvol"
449 
450  // getting the bounding box dimensions along z so as to set the
451  // appropriate upstream generation surface for the NuMI flux driver
452 
453  // RWH 2010-07-16: do not try to automatically get zmin from geometry, rather
454  // by default let the flux start from the window. If the user wants to
455  // override this then they need to explicitly set a "zmin". Trying to use
456  // the geometry is fraught with problems in local vs. global coordinates and
457  // units where it can appear to work in some cases but it actually isn't really
458  // universally correct.
459  //was// TGeoShape * bounding_box = topvol->GetShape();
460  //was// bounding_box->GetAxisRange(3, zmin, zmax);
461  //was// zmin *= rgeom->LengthUnits();
462  //was// zmax *= rgeom->LengthUnits();
463 
464  // switch on/off volumes as requested
465  if ( (gOptRootGeomTopVol[0] == '+') || (gOptRootGeomTopVol[0] == '-') ) {
466  bool exhaust = (*gOptRootGeomTopVol.c_str() == '+');
468  }
469 
470  // casting to the GENIE geometry driver interface
471  geom_driver = dynamic_cast<GeomAnalyzerI *> (rgeom);
472 
473  // user specifid a fiducial volume cut ... parse that out
474  if ( gOptFidCut.find("rock") != std::string::npos )
476  else if ( gOptFidCut != "" ) CreateFidSelection(gOptFidCut,rgeom);
477 
478  }
479  else {
480  //
481  // *** Using a 'point' geometry with the specified target mix
482  // *** ( = a list of targets with their corresponding weight fraction)
483  //
484 
485  // creating & configuring a point geometry driver
488  // casting to the GENIE geometry driver interface
489  geom_driver = dynamic_cast<GeomAnalyzerI *> (pgeom);
490  }
491 
492  // *************************************************************************
493  // * Create / configure the flux driver
494  // *************************************************************************
495  GFluxI * flux_driver =
497  if ( ! flux_driver ) {
498  // couldn't get the requested flux driver ?
499  std::ostringstream s;
500  s << "Known FluxDrivers:\n";
501  const std::vector<std::string>& known =
503  std::vector<std::string>::const_iterator itr = known.begin();
504  for ( ; itr != known.end(); ++itr ) s << " " << (*itr) << "\n";
505  LOG("gevgen_fnal", pFATAL)
506  << "Failed to get any flux driver from GFluxDriverFactory\n"
507  << "when using \"" << gOptFluxDriver << "\"\n" << s.str();
508  exit(1);
509  }
510 
511  if ( ! gOptUsingHistFlux ) {
512  genie::flux::GFluxFileConfigI* flux_file_config =
513  dynamic_cast<genie::flux::GFluxFileConfigI*>(flux_driver);
514  if ( ! flux_file_config ) {
515  LOG("gevgen_fnal", pFATAL)
516  << "Failed to get GFluxFileConfigI driver from GFluxDriverFactory\n"
517  << "when using \"" << gOptFluxDriver << "\"";
518  exit(1);
519  }
520 
521  //
522  // *** Using the detailed ntuple neutrino flux description
523  //
525  flux_file_config->SetUpstreamZ(gOptZmin); // was "zmin" from bounding_box
526  flux_file_config->SetNumOfCycles(0);
527 
528  if ( gOptFluxPdg.size() > 0 ) {
529  // user specified list of neutrino PDGs
530  flux_file_config->SetFluxParticles(gOptFluxPdg);
531  std::ostringstream s;
532  PDGCodeList::const_iterator itr = gOptFluxPdg.begin();
533  for ( ; itr != gOptFluxPdg.end(); ++itr) s << (*itr) << " ";
534  LOG("gevgen_fnal", pNOTICE)
535  << "Limiting to nu PDGs: " << s.str();
536  }
537  }
538  else {
539  //
540  // *** Using fluxes from histograms (for all specified neutrino species)
541  //
542  genie::flux::GCylindTH1Flux * hist_flux_driver =
543  dynamic_cast<genie::flux::GCylindTH1Flux*>(flux_driver);
544  if ( ! hist_flux_driver ) {
545  LOG("gevgen_fnal", pFATAL)
546  << "Failed to get GCylinderTH1Flux driver from GFluxDriverFactory\n"
547  << "when using " << gOptFluxDriver;
548  exit(1);
549  }
550 
551  // creating & configuring a generic GCylindTH1Flux flux driver
552  TVector3 bdir (0,0,1); // dir along +z
553  TVector3 bspot(0,0,0);
554 
555  hist_flux_driver->SetNuDirection (bdir);
556  hist_flux_driver->SetBeamSpot (bspot);
557  hist_flux_driver->SetTransverseRadius (-1);
558  map<int,TH1D*>::iterator it = gOptFluxHst.begin();
559  for( ; it != gOptFluxHst.end(); ++it) {
560  int pdg_code = it->first;
561  TH1D * spectrum = it->second;
562  hist_flux_driver->AddEnergySpectrum(pdg_code, spectrum);
563  // once the histogram has been added to the GCylindTH1Flux driver
564  // it is owned by the driver and it is up to the the driver
565  // to clean up (i.e. delete it).
566  // remove it from this map to avoid double deletion.
567  it->second = 0;
568  }
569  }
570 
571  // these might come in handy ... avoid repeated dynamic_cast<> calls
572  genie::flux::GFluxExposureI* fluxExposureI =
573  dynamic_cast<genie::flux::GFluxExposureI*>(flux_driver);
574  genie::flux::GFluxFileConfigI* fluxFileConfigI =
575  dynamic_cast<genie::flux::GFluxFileConfigI*>(flux_driver);
576 
577 
578  // *************************************************************************
579  // * Handle chicken/egg problem: geom analyzer vs. flux.
580  // * Need both at this point change geom scan defaults.
581  // *************************************************************************
583 
585  dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
586  if ( ! rgeom ) assert(0);
587 
588  rgeom -> SetDebugFlags(gOptDebug);
589 
590  // even if user doesn't specify gOptNScan configure to scan using flux
591  if ( gOptNScan >= 0 ) {
592  LOG("gevgen_fnal", pNOTICE)
593  << "Using ROOTGeomAnalyzer: geom scan using flux: nparticles=" << gOptNScan;
594  rgeom->SetScannerFlux(flux_driver);
595  if ( gOptNScan > 0 ) rgeom->SetScannerNParticles(gOptNScan);
596  } else {
597  int nabs = TMath::Abs(gOptNScan);
598  LOG("gevgen_fnal", pNOTICE)
599  << "Using ROOTGeomAnalyzer: geom scan using box: npoints=nrays=" << nabs;
600  rgeom->SetScannerNPoints(nabs);
601  rgeom->SetScannerNRays(nabs);
602  }
603  }
604 
605  // *************************************************************************
606  // * Create/configure the event generation driver
607  // *************************************************************************
608  GMCJDriver * mcj_driver = new GMCJDriver;
610  mcj_driver->UseFluxDriver(flux_driver);
611  mcj_driver->UseGeomAnalyzer(geom_driver);
612  if ( ( gOptExtMaxPlXml != "" ) && ! gOptWriteMaxPlXml ) {
613  mcj_driver->UseMaxPathLengths(gOptExtMaxPlXml);
614  }
615  mcj_driver->Configure();
616  mcj_driver->UseSplines();
617  mcj_driver->ForceSingleProbScale();
618 
619  if ( ( gOptExtMaxPlXml != "" ) && gOptWriteMaxPlXml ) {
621  dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
622  if ( rgeom ) {
623  const genie::PathLengthList& maxpath = rgeom->GetMaxPathLengths();
624  std::string maxplfile = gOptExtMaxPlXml;
625  maxpath.SaveAsXml(maxplfile);
626  // append extra info to file
627  std::ofstream mpfile(maxplfile.c_str(), std::ios_base::app);
628  mpfile
629  << std::endl
630  << "<!-- this file is only relevant for a setup compatible with:"
631  << std::endl
632  << "geom: " << gOptRootGeom << " top: \"" << gOptRootGeomTopVol << "\""
633  << std::endl
634  << "flux: " << gOptFluxFile
635  << std::endl
636  << "location: " << gOptDetectorLocation
637  << std::endl
638  << "fidcut: " << gOptFidCut
639  << std::endl
640  << "nscan: " << gOptNScan << " (0>= use flux, <0 use box |nscan| points/rays)"
641  << std::endl
642  << "zmin: " << gOptZmin << " (if |zmin| > 1e30, leave ray on flux window)"
643  << std::endl
644  << "-->"
645  << std::endl;
646  mpfile.close();
647  }
648  }
649 
650  // *************************************************************************
651  // * Prepare for writing the output event tree & status file
652  // *************************************************************************
653 
654  // Initialize an Ntuple Writer to save GHEP records into a TTree
657  ntpw.Initialize();
658 
659 
660  std::vector<TBranch*> extraBranches;
661  std::vector<std::string> branchNames;
662  std::vector<std::string> branchClassNames;
663  std::vector<void**> branchObjPointers;
664 
665  // Add custom branch(s) to the standard GENIE event tree so that
666  // info on the flux neutrino parent particle can be passed-through
667  if ( fluxFileConfigI ) {
668  fluxFileConfigI->GetBranchInfo(branchNames,branchClassNames,
669  branchObjPointers);
670  size_t nn = branchNames.size();
671  size_t nc = branchClassNames.size();
672  size_t np = branchObjPointers.size();
673  if ( nn != nc || nc != np ) {
674  LOG("gevgen_fnal", pERROR)
675  << "Inconsistent info back from \"" << gOptFluxDriver << "\" "
676  << "for branch info: " << nn << " " << nc << " " << np;
677  } else {
678  for (size_t ii = 0; ii < nn; ++ii) {
679  const char* bname = branchNames[ii].c_str();
680  const char* cname = branchClassNames[ii].c_str();
681  void**& optr = branchObjPointers[ii]; // note critical '&' !
682  if ( ! optr || ! *optr ) continue; // no pointer supplied, skip it
683  int split = 99; // 1
684  LOG("gevgen_fnal", pNOTICE)
685  << "Adding extra branch \"" << bname << "\" of type \""
686  << cname << "\" (" << optr << ") to output tree";
687  TBranch* bptr = ntpw.EventTree()->Branch(bname,cname,optr,32000,split);
688  extraBranches.push_back(bptr);
689 
690  if ( bptr ) {
691  // don't delete this !!! we're sharing
692  bptr->SetAutoDelete(false);
693  } else {
694  LOG("gevgen_fnal", pERROR)
695  << "FAILED to add extra branch \"" << bname << "\" of type \""
696  << cname << "\" to output tree";
697  }
698  } // loop over additions
699  } // same # of entries
700  } // of genie::flux::GFluxFileConfigI type
701 
702  // Create a MC job monitor for a periodically updated status file
703  GMCJMonitor mcjmonitor(gOptRunNu);
704  mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());
705 
706  // *************************************************************************
707  // * Event generation loop
708  // *************************************************************************
709 
710  // define handler to allow signal to end job gracefully
711  signal(SIGTERM,gsSIGTERMhandler);
712 
713  int ievent = 0;
714  while ( ! gSigTERM )
715  {
716  LOG("gevgen_fnal", pINFO)
717  << " *** Generating event............ " << ievent;
718 
719  // In case the required statistics was expressed as 'number of events'
720  // then quit if that number has been generated
721  if ( ievent == gOptNev ) break;
722 
723  // In case the required statistics was expressed as 'number of POT'
724  // then exit the event loop if the requested POT has been generated.
725  if ( gOptPOT > 0 && fluxExposureI ) {
726  double fpot = fluxExposureI->GetTotalExposure(); // current POTs used
727  double psc = mcj_driver->GlobProbScale(); // interaction prob. scale
728  double pot = fpot / psc; // POTs for generated sample
729  if ( pot >= gOptPOT ) break;
730  }
731 
732  // Generate a single event using neutrinos coming from the specified flux
733  // and hitting the specified geometry or target mix
734  EventRecord * event = mcj_driver->GenerateEvent();
735 
736  // Check whether a null event was returned due to the flux driver reaching
737  // the end of the input flux ntuple - exit the event generation loop
738  if ( ! event && flux_driver->End() ) {
739  LOG("gevgen_fnal", pWARN)
740  << "** The flux driver read all the input flux entries: End()==true";
741  break;
742  }
743  if ( ! event ) {
744  LOG("gevgen_fnal", pERROR)
745  << "Got a null generated neutino event! Retrying ...";
746  continue;
747  }
748  LOG("gevgen_fnal", pINFO)
749  << "Generated event: " << *event;
750 
751  // A valid event was generated: flux info (parent decay/prod
752  // position/kinematics) for that simulated event should already
753  // be connected to the right output tree branch
754 
755  // Add event at the output ntuple, refresh the mc job monitor & clean-up
756  ntpw.AddEventRecord(ievent, event);
757  mcjmonitor.Update(ievent,event);
758  delete event;
759  ievent++;
760 
761  } //1
762 
763  // Copy metadata tree, if available
764  if ( fluxFileConfigI ) {
765  TTree* t1 = fluxFileConfigI->GetMetaDataTree();
766  if ( t1 ) {
767  size_t nmeta = t1->GetEntries();
768  TTree* t2 = (TTree*)t1->Clone(0);
769  for (size_t i = 0; i < nmeta; ++i) {
770  t1->GetEntry(i);
771  t2->Fill();
772  }
773  t2->Write();
774  }
775  }
776 
777  LOG("gevgen_fnal", pINFO)
778  << "The GENIE MC job is done generating events - Cleaning up & exiting...";
779 
780  // *************************************************************************
781  // * Print job statistics &
782  // * calculate normalization factor for the generated sample
783  // *************************************************************************
785  // POT normalization will only be calculated if event generation was based
786  // on beam simulation ntuples (not just histograms) & a detailed detector
787  // geometry description.
788  // Get nunber of flux neutrinos read-in by flux driver, number of flux
789  // neutrinos actually thrown to the event generation driver and number
790  // of neutrino interactions actually generated
791  long int nflx = 0;
792  long int nflx_evg = mcj_driver-> NFluxNeutrinos();
793  double fpot = 0;
794  const char* exposureUnits = "(unknown units)";
795  if ( fluxExposureI ) {
796  fpot = fluxExposureI->GetTotalExposure(); // POTs used so far
797  nflx = fluxExposureI->NFluxNeutrinos();
798  //genie::flux::Exposure_t etype = fluxExposureI->GetExposureType();
799  //exposureUnits = genie::flux::GFluxExposureI::AsString(etype);
800  exposureUnits = fluxExposureI->GetExposureUnits();
801  }
802  if ( fluxFileConfigI ) {
803  fluxFileConfigI->PrintConfig();
804  }
805  double psc = mcj_driver->GlobProbScale(); // interaction prob. scale
806  if ( psc <= 0.0 ) {
807  LOG("gevgen_fnal", pFATAL) << "MCJobDriver GlobalProbScale was " << psc;
808  }
809  double pot = fpot / psc; // POT for generated sample
810  long int nev = ievent;
811 
812  LOG("gevgen_fnal", pNOTICE)
813  << "\n >> Interaction probability scaling factor: " << psc
814  << "\n >> using: " << gOptFluxDriver
815  << "\n >> N of flux v read-in by flux driver: " << nflx
816  << "\n >> N of flux v thrown to event gen driver: " << nflx_evg
817  << "\n >> N of generated v interactions: " << nev
818  << "\n ** Normalization for generated sample: " << pot
819  << " " << exposureUnits << " * detector";
820 
821  ntpw.EventTree()->SetWeight(pot); // store POT
822 
823  }
824 
825  // *************************************************************************
826  // * Save & clean-up
827  // *************************************************************************
828 
829  // Save the generated event tree & close the output file
830  ntpw.Save();
831 
832  // Clean-up
833  delete geom_driver;
834  delete flux_driver;
835  delete mcj_driver;
836  // this list should only be histograms that have (for some reason)
837  // not been handed over to the GCylindTH1Flux driver.
838  map<int,TH1D*>::iterator it = gOptFluxHst.begin();
839  for( ; it != gOptFluxHst.end(); ++it) {
840  TH1D * spectrum = it->second;
841  if(spectrum) delete spectrum;
842  }
843  gOptFluxHst.clear();
844 
845  LOG("gevgen_fnal", pNOTICE) << "Done!";
846 
847  return 0;
848 }
849 
850 //____________________________________________________________________________
852 {
853  /// potentially load extra libraries that might extend the list of
854  /// potential flux drivers, and how to map short names to classes ...
855  // we really should at this point read some external file to get
856  // an expandable list of libraries ... but for now just hard code it
857 
858  vector<string> extraLibs;
859 
860  ///***** this part should come from reading an external file
861  /// placeholder file $GENIE/config/FluxDriverExpansion.xml
862 
863  extraLibs.push_back("libdk2nuTree");
864  extraLibs.push_back("libdk2nuGenie");
865 
866  // what one might expect to find at the beginning of -f <arg>
867 
868  gOptFluxShortNames["histo"] = "genie::flux::GCylindTH1Flux";
869  gOptFluxShortNames["hist"] = "genie::flux::GCylindTH1Flux";
870 
871  gOptFluxShortNames["simple"] = "genie::flux::GSimpleNtpFlux";
872  gOptFluxShortNames["numi"] = "genie::flux::GNuMIFlux";
873  gOptFluxShortNames["dk2nu"] = "genie::flux::GDk2NuFlux";
874 
875  ///******* done with fake "read"
876 
877  // see if there are any libraries to load
878  // just let ROOT do it ... check error code on return
879  // but tweak ROOT's ERROR message output so we don't see huge messages
880  // for failures
881  // gErrorIgnoreLevel to kError (from TError.h)
882 
883  Int_t prevErrorLevel = gErrorIgnoreLevel;
884  gErrorIgnoreLevel = kFatal;
885  for (size_t ilib=0; ilib < extraLibs.size(); ++ilib) {
886  // library names should be like libXYZZY without extension [e.g. .so]
887  // but with the leading "lib"
888  int loadStatus = gSystem->Load(extraLibs[ilib].c_str());
889  const char* statWords = "failed to load";
890  if ( loadStatus == 0 ) { statWords = "successfully loaded"; }
891  else if ( loadStatus == 1 ) { statWords = "already had loaded"; }
892  else if ( loadStatus == -1 ) { statWords = "couldn't find library"; }
893  else if ( loadStatus == -2 ) { statWords = "mismatched version"; }
894 
895  LOG("gevgen_fnal",pNOTICE)
896  << statWords << " (" << loadStatus << ") " << extraLibs[ilib];
897  }
898  // restore the ROOT error message level
899  gErrorIgnoreLevel = prevErrorLevel;
900 
901 }
902 
903 //____________________________________________________________________________
904 void GetCommandLineArgs(int argc, char ** argv)
905 {
906  LOG("gevgen_fnal", pINFO) << "Parsing command line arguments";
907 
908  // Common run options. Set defaults and read.
911 
912  // Parse run options for this app
913 
914  CmdLnArgParser parser(argc,argv);
915 
916  // help?
917  bool help = parser.OptionExists('h');
918  if(help) {
919  PrintSyntax();
920  exit(0);
921  }
922 
923  // run number:
924  if ( parser.OptionExists('r') ) {
925  LOG("gevgen_fnal", pDEBUG) << "Reading MC run number";
926  gOptRunNu = parser.ArgAsLong('r');
927  } else {
928  LOG("gevgen_fnal", pDEBUG)
929  << "Unspecified run number - Using default";
930  gOptRunNu = 0;
931  } //-r
932 
933  //
934  // *** geometry
935  //
936 
937  string geom = "";
938  string lunits, dunits;
939  if( parser.OptionExists('g') ) {
940  LOG("gevgen_fnal", pDEBUG) << "Getting input geometry";
941  geom = parser.ArgAsString('g');
942 
943  // is it a ROOT file that contains a ROOT geometry?
944  bool accessible_geom_file =
945  ! (gSystem->AccessPathName(geom.c_str()));
946  if (accessible_geom_file) {
947  gOptRootGeom = geom;
948  gOptUsingRootGeom = true;
949  }
950  } else {
951  LOG("gevgen_fnal", pFATAL)
952  << "No geometry option specified - Exiting";
953  PrintSyntax();
954  exit(1);
955  } //-g
956 
957  if(gOptUsingRootGeom) {
958  // using a ROOT geometry - get requested geometry units
959 
960  // legth units:
961  if( parser.OptionExists('L') ) {
962  LOG("gevgen_fnal", pDEBUG)
963  << "Checking for input geometry length units";
964  lunits = parser.ArgAsString('L');
965  } else {
966  LOG("gevgen_fnal", pDEBUG) << "Using default geometry length units";
967  lunits = kDefOptGeomLUnits;
968  } // -L
969  // density units:
970  if( parser.OptionExists('D') ) {
971  LOG("gevgen_fnal", pDEBUG)
972  << "Checking for input geometry density units";
973  dunits = parser.ArgAsString('D');
974  } else {
975  LOG("gevgen_fnal", pDEBUG) << "Using default geometry density units";
976  dunits = kDefOptGeomDUnits;
977  } // -D
980 
981  // check whether an event generation volume name has been
982  // specified -- default is the 'top volume'
983  if( parser.OptionExists('t') ) {
984  LOG("gevgen_fnal", pDEBUG) << "Checking for input volume name";
985  gOptRootGeomTopVol = parser.ArgAsString('t');
986  } else {
987  LOG("gevgen_fnal", pDEBUG) << "Using the <master volume>";
988  } // -t
989 
990  // check whether an XML file with the maximum (density weighted)
991  // path lengths for each detector material is specified -
992  // otherwise will compute the max path lengths at job init
993  // if passed name starts with '+', then compute max at job init, but write out the result
994  if ( parser.OptionExists('m') ) {
995  LOG("gevgen_fnal", pDEBUG)
996  << "Checking for maximum path lengths XML file";
997  gOptExtMaxPlXml = parser.ArgAsString('m');
998  gOptWriteMaxPlXml = false;
999  if ( gOptExtMaxPlXml[0] == '+' ) {
1000  gOptExtMaxPlXml = gOptExtMaxPlXml.substr(1,std::string::npos);
1001  gOptWriteMaxPlXml = true;
1002  LOG("gevgen_fnal", pINFO)
1003  << "Will write maximum path lengths XML file: " << gOptExtMaxPlXml;
1004  }
1005  } else {
1006  LOG("gevgen_fnal", pDEBUG)
1007  << "Will compute the maximum path lengths at job init";
1008  gOptExtMaxPlXml = "";
1009  } // -m
1010 
1011  // fidcut:
1012  if( parser.OptionExists('F') ) {
1013  LOG("gevgen_fnal", pDEBUG) << "Using Fiducial cut?";
1014  gOptFidCut = parser.ArgAsString('F');
1015  } else {
1016  LOG("gevgen_fnal", pDEBUG) << "No fiducial volume cut";
1017  gOptFidCut = "";
1018  } //-F
1019 
1020  if(!gOptUsingHistFlux) {
1021  // how to scan the geometry (if relevant)
1022  if( parser.OptionExists('S') ) {
1023  LOG("gevgen_fnal", pDEBUG) << "Reading requested geom scan count";
1024  gOptNScan = parser.ArgAsInt('S');
1025  } else {
1026  LOG("gevgen_fnal", pDEBUG) << "No geom scan count was requested";
1027  gOptNScan = 0;
1028  } //-S
1029 
1030  // z for flux rays to start
1031  if( parser.OptionExists('z') ) {
1032  LOG("gevgen_fnal", pDEBUG) << "Reading requested zmin";
1033  gOptZmin = parser.ArgAsDouble('z');
1034  } else {
1035  LOG("gevgen_fnal", pDEBUG) << "No zmin was requested";
1036  gOptZmin = -2.0e30; // < -1.0e30 ==> leave it on flux window
1037  } //-z
1038 
1039  // debug flags
1040  if ( parser.OptionExists('d') ) {
1041  LOG("gevgen_fnal", pDEBUG) << "Reading debug flag value";
1042  gOptDebug = parser.ArgAsInt('d');
1043  } else {
1044  LOG("gevgen_fnal", pDEBUG) << "Unspecified debug flags - Using default";
1045  gOptDebug = 0;
1046  } //-d
1047 
1048  } // root geom && gnumi flux
1049 
1050  } // using root geom?
1051 
1052  else {
1053  // User has specified a target mix.
1054  // Decode the list of target pdf codes & their corresponding weight fraction
1055  // (specified as 'pdg_code_1[fraction_1],pdg_code_2[fraction_2],...')
1056  // See documentation on top section of this file.
1057  //
1058  gOptTgtMix.clear();
1059  vector<string> tgtmix = utils::str::Split(geom,",");
1060  if(tgtmix.size()==1) {
1061  int pdg = atoi(tgtmix[0].c_str());
1062  double wgt = 1.0;
1063  gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
1064  } else {
1065  vector<string>::const_iterator tgtmix_iter = tgtmix.begin();
1066  for( ; tgtmix_iter != tgtmix.end(); ++tgtmix_iter) {
1067  string tgt_with_wgt = *tgtmix_iter;
1068  string::size_type open_bracket = tgt_with_wgt.find("[");
1069  string::size_type close_bracket = tgt_with_wgt.find("]");
1070  if (open_bracket ==string::npos ||
1071  close_bracket==string::npos)
1072  {
1073  LOG("gevgen_fnal", pFATAL)
1074  << "You made an error in specifying the target mix";
1075  PrintSyntax();
1076  exit(1);
1077  }
1078  string::size_type ibeg = 0;
1079  string::size_type iend = open_bracket;
1080  string::size_type jbeg = open_bracket+1;
1081  string::size_type jend = close_bracket;
1082  int pdg = atoi(tgt_with_wgt.substr(ibeg,iend-ibeg).c_str());
1083  double wgt = atof(tgt_with_wgt.substr(jbeg,jend-jbeg).c_str());
1084  LOG("gevgen_fnal", pDEBUG)
1085  << "Adding to target mix: pdg = " << pdg << ", wgt = " << wgt;
1086  gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
1087 
1088  }// tgtmix_iter
1089  } // >1 materials in mix
1090  } // using tgt mix?
1091 
1092  //
1093  // *** flux
1094  //
1095  if ( parser.OptionExists('f') ) {
1096  LOG("gevgen_fnal", pDEBUG) << "Getting input flux";
1097  DetermineFluxDriver(parser.ArgAsString('f'));
1098  } else {
1099  LOG("gevgen_fnal", pFATAL) << "No flux info was specified - Exiting";
1100  PrintSyntax();
1101  exit(1);
1102  }
1103 
1104  // number of events to generate
1105  if( parser.OptionExists('n') ) {
1106  LOG("gevgen_fnal", pDEBUG)
1107  << "Reading limit on number of events to generate";
1108  gOptNev = parser.ArgAsInt('n');
1109  } else {
1110  LOG("gevgen_fnal", pDEBUG)
1111  << "Will keep on generating events till the flux driver stops";
1112  gOptNev = -1;
1113  } //-n
1114 
1115  // statistics to generate in terms of POT
1116  if( parser.OptionExists('e') ) {
1117  LOG("gevgen_fnal", pDEBUG) << "Reading requested exposure in POT";
1118  gOptPOT = parser.ArgAsDouble('e');
1119  } else {
1120  LOG("gevgen_fnal", pDEBUG) << "No POT exposure was requested";
1121  gOptPOT = -1;
1122  } //-e
1123 
1124  // event file prefix
1125  if( parser.OptionExists('o') ) {
1126  LOG("gevgen_fnal", pDEBUG) << "Reading the event filename prefix";
1127  gOptEvFilePrefix = parser.ArgAsString('o');
1128  } else {
1129  LOG("gevgen_fnal", pDEBUG)
1130  << "Will set the default event filename prefix";
1132  } //-o
1133 
1134 
1135  // random number seed
1136  if( parser.OptionExists("seed") ) {
1137  LOG("gevgen_fnal", pINFO) << "Reading random number seed";
1138  gOptRanSeed = parser.ArgAsLong("seed");
1139  } else {
1140  LOG("gevgen_fnal", pINFO) << "Unspecified random number seed - Using default";
1141  gOptRanSeed = -1;
1142  }
1143 
1144  // input cross-section file
1145  if( parser.OptionExists("cross-sections") ) {
1146  LOG("gevgen_fnal", pINFO) << "Reading cross-section file";
1147  gOptInpXSecFile = parser.ArgAsString("cross-sections");
1148  } else {
1149  LOG("gevgen_fnal", pINFO) << "Unspecified cross-section file";
1150  gOptInpXSecFile = "";
1151  }
1152 
1153 
1154  //
1155  // >>> perform 'sanity' checks on command line arguments
1156  //
1157 
1158  // Tthe 'exposure' may be set either as:
1159  // - a number of POTs
1160  // - a number of generated events
1161  // Only one of those options can be set.
1162  if(!gOptUsingHistFlux) {
1163  int nset=0;
1164  if(gOptPOT > 0) nset++;
1165  if(gOptNev > 0) nset++;
1166  if(nset==0) {
1167  LOG("gevgen_fnal", pFATAL)
1168  << "** To use a gNuMI flux ntuple you need to specify an exposure, "
1169  << "either via the -e or -n options";
1170  PrintSyntax();
1171  exit(1);
1172  }
1173  if(nset>1) {
1174  LOG("gevgen_fnal", pFATAL)
1175  << "You can not specify more than one of the -e or -n options";
1176  PrintSyntax();
1177  exit(1);
1178  }
1179  }
1180  // If we use a flux histograms (not flux ntuples) then -currently- the
1181  // only way to control exposure is via a number of events
1182  if(gOptUsingHistFlux) {
1183  if(gOptNev < 0) {
1184  LOG("gevgen_fnal", pFATAL)
1185  << "If you're using flux from histograms you need to specify the -n option";
1186  PrintSyntax();
1187  exit(1);
1188  }
1189  }
1190  // If we don't use a detailed ROOT detector geometry (but just a target mix)
1191  // then don't accept POT as a way to control job statistics (not enough info
1192  // is passed in the target mix to compute POT & the calculation can be easily
1193  // done offline)
1194  if(!gOptUsingRootGeom) {
1195  if(gOptPOT > 0) {
1196  LOG("gevgen_fnal", pFATAL)
1197  << "You may not use the -e option without a detector geometry description";
1198  exit(1);
1199  }
1200  }
1201 
1202  //
1203  // >>> print the command line options
1204  //
1205 
1206  PDGLibrary * pdglib = PDGLibrary::Instance();
1207 
1208  ostringstream gminfo;
1209  if (gOptUsingRootGeom) {
1210  gminfo << "Using ROOT geometry - file = " << gOptRootGeom
1211  << ", top volume = "
1212  << ((gOptRootGeomTopVol.size()==0) ? "<master volume>" : gOptRootGeomTopVol)
1213  << ", max{PL} file = "
1214  << ((gOptExtMaxPlXml.size()==0) ? "<none>" : gOptExtMaxPlXml)
1215  << ", length units = " << lunits
1216  << ", density units = " << dunits;
1217  } else {
1218  gminfo << "Using target mix: ";
1219  map<int,double>::const_iterator iter;
1220  for(iter = gOptTgtMix.begin(); iter != gOptTgtMix.end(); ++iter) {
1221  int pdg_code = iter->first;
1222  double wgt = iter->second;
1223  TParticlePDG * p = pdglib->Find(pdg_code);
1224  if(p) {
1225  string name = p->GetName();
1226  gminfo << "(" << name << ") -> " << 100*wgt << "% / ";
1227  }//p?
1228  }
1229  }
1230 
1231  ostringstream fluxinfo;
1232  if (gOptUsingHistFlux) {
1233  fluxinfo << "Using histograms: ";
1234  map<int,TH1D*>::const_iterator iter;
1235  for(iter = gOptFluxHst.begin(); iter != gOptFluxHst.end(); ++iter) {
1236  int pdg_code = iter->first;
1237  TH1D * spectrum = iter->second;
1238  TParticlePDG * p = pdglib->Find(pdg_code);
1239  if(p) {
1240  string name = p->GetName();
1241  fluxinfo << "(" << name << ") -> " << spectrum->GetName() << " / ";
1242  }//p?
1243  }
1244  } else {
1245  fluxinfo << "Using " << gOptFluxDriver << " flux driver- "
1246  << "file = " << gOptFluxFile
1247  << ", location = " << gOptDetectorLocation;
1248  }
1249 
1250  ostringstream exposure;
1251  if(gOptPOT > 0)
1252  exposure << "Number of POTs = " << gOptPOT;
1253  if(gOptNev > 0)
1254  exposure << "Number of events = " << gOptNev;
1255 
1256 
1257  LOG("gevgen_fnal", pNOTICE)
1258  << "\n\n"
1259  << utils::print::PrintFramedMesg("FNAL expt. event generation job configuration");
1260 
1261  LOG("gevgen_fnal", pNOTICE)
1262  << "\n - Run number: " << gOptRunNu
1263  << "\n - Random number seed: " << gOptRanSeed
1264  << "\n - Using cross-section file: " << gOptInpXSecFile
1265  << "\n - Flux @ " << fluxinfo.str()
1266  << "\n - Geometry @ " << gminfo.str()
1267  << "\n - Exposure @ " << exposure.str();
1268 
1269  LOG("gevgen_fnal", pNOTICE) << *RunOpt::Instance();
1270 }
1271 //____________________________________________________________________________
1272 void PrintSyntax(void)
1273 {
1274  LOG("gevgen_fnal", pFATAL)
1275  << "\n **Syntax**"
1276  << "\n gevgen_fnal [-h] [-r run#]"
1277  << "\n -f flux -g geometry"
1278  << "\n [-t top_volume_name_at_geom] [-m max_path_lengths_xml_file]"
1279  << "\n [-L length_units_at_geom] [-D density_units_at_geom]"
1280  << "\n [-n n_of_events] [-e exposure_in_POTs]"
1281  << "\n [-o output_event_file_prefix]"
1282  << "\n [-F fid_cut_string] [-S nrays_scan]"
1283  << "\n [-z zmin_start]"
1284  << "\n [--seed random_number_seed]"
1285  << "\n --cross-sections xml_file"
1286  << "\n [--event-generator-list list_name]"
1287  << "\n [--message-thresholds xml_file]"
1288  << "\n [--unphysical-event-mask mask]"
1289  << "\n [--event-record-print-level level]"
1290  << "\n [--mc-job-status-refresh-rate rate]"
1291  << "\n [--cache-file root_file]"
1292  << "\n"
1293  << " Please also read the detailed documentation at "
1294  << "$GENIE/src/Apps/gFNALExptEvGen.cxx"
1295  << "\n";
1296 }
1297 //____________________________________________________________________________
1298 void CreateFidSelection (string fidcut, GeomAnalyzerI* geom_driver)
1299 {
1300  ///
1301  /// User defined fiducial volume cut
1302  /// [0][M]<SHAPE>:val1,val2,...
1303  /// "0" means reverse the cut (i.e. exclude the volume)
1304  /// "M" means the coordinates are given in the ROOT geometry
1305  /// "master" system and need to be transformed to "top vol" system
1306  /// <SHAPE> can be any of "zcyl" "box" "zpoly" "sphere"
1307  /// [each takes different # of args]
1308  /// This must be followed by a ":" and a list of values separated by punctuation
1309  /// (allowed separators: commas , parentheses () braces {} or brackets [] )
1310  /// Value mapping:
1311  /// zcly:x0,y0,radius,zmin,zmax - cylinder along z at (x0,y0) capped at z's
1312  /// box:xmin,ymin,zmin,xmax,ymax,zmax - box w/ upper & lower extremes
1313  /// zpoly:nfaces,x0,y0,r_in,phi,zmin,zmax - nfaces sided polygon in x-y plane
1314  // sphere:x0,y0,z0,radius - sphere of fixed radius at (x0,y0,z0)
1315  /// Examples:
1316  /// 1) 0mbox:0,0,0.25,1,1,8.75
1317  /// exclude (i.e. reverse) a box in master coordinates w/ corners (0,0,0.25) (1,1,8.75)
1318  /// 2) mzpoly:6,(2,-1),1.75,0,{0.25,8.75}
1319  /// six sided polygon in x-y plane, centered at x,y=(2,-1) w/ inscribed radius 1.75
1320  /// no rotation (so first face is in y-z plane +r from center, i.e. hex sits on point)
1321  /// limited to the z range of {0.25,8.75} in the master ROOT geom coordinates
1322  /// 3) zcly:(3,4),5.5,-2,10
1323  /// a cylinder oriented parallel to the z axis in the "top vol" coordinates
1324  /// at x,y=(3,4) with radius 5.5 and z range of {-2,10}
1325  ///
1326  geometry::ROOTGeomAnalyzer * rgeom =
1327  dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
1328  if ( ! rgeom ) {
1329  LOG("gevgen_fnal", pWARN)
1330  << "Can not create GeomVolSelectorFiduction,"
1331  << " geometry driver is not ROOTGeomAnalyzer";
1332  return;
1333  }
1334 
1335  LOG("gevgen_fnal", pNOTICE) << "-F " << fidcut;
1336 
1339 
1340  fidsel->SetRemoveEntries(true); // drop segments that won't be considered
1341 
1342  // convert string to lowercase
1343  std::transform(fidcut.begin(),fidcut.end(),fidcut.begin(),::tolower);
1344 
1345  vector<string> strtok = genie::utils::str::Split(fidcut,":");
1346  if ( strtok.size() != 2 ) {
1347  LOG("gevgen_fnal", pWARN)
1348  << "Can not create GeomVolSelectorFiduction,"
1349  << " no \":\" separating type from values. nsplit=" << strtok.size();
1350  for ( unsigned int i=0; i < strtok.size(); ++i )
1351  LOG("gevgen_fnal",pNOTICE)
1352  << "strtok[" << i << "] = \"" << strtok[i] << "\"";
1353  return;
1354  }
1355 
1356  // parse out optional "x" and "m"
1357  string stype = strtok[0];
1358  bool reverse = ( stype.find("0") != string::npos );
1359  bool master = ( stype.find("m") != string::npos ); // action after values are set
1360 
1361  // parse out values
1362  vector<double> vals;
1363  vector<string> valstrs = genie::utils::str::Split(strtok[1]," ,;(){}[]");
1364  vector<string>::const_iterator iter = valstrs.begin();
1365  for ( ; iter != valstrs.end(); ++iter ) {
1366  const string& valstr1 = *iter;
1367  if ( valstr1 != "" ) vals.push_back(atof(valstr1.c_str()));
1368  }
1369  size_t nvals = vals.size();
1370 
1371  std::cout << "ivals = [";
1372  for (unsigned int i=0; i < nvals; ++i) {
1373  if (i>0) cout << ",";
1374  std::cout << vals[i];
1375  }
1376  std::cout << "]" << std::endl;
1377 
1378  // std::vector elements are required to be adjacent so we can treat address as ptr
1379 
1380  if ( stype.find("zcyl") != string::npos ) {
1381  // cylinder along z direction at (x0,y0) radius zmin zmax
1382  if ( nvals < 5 )
1383  LOG("gevgen_fnal", pFATAL) << "MakeZCylinder needs 5 values, not " << nvals
1384  << " fidcut=\"" << fidcut << "\"";
1385  fidsel->MakeZCylinder(vals[0],vals[1],vals[2],vals[3],vals[4]);
1386 
1387  } else if ( stype.find("box") != string::npos ) {
1388  // box (xmin,ymin,zmin) (xmax,ymax,zmax)
1389  if ( nvals < 6 )
1390  LOG("gevgen_fnal", pFATAL) << "MakeBox needs 6 values, not " << nvals
1391  << " fidcut=\"" << fidcut << "\"";
1392  double xyzmin[3] = { vals[0], vals[1], vals[2] };
1393  double xyzmax[3] = { vals[3], vals[4], vals[5] };
1394  fidsel->MakeBox(xyzmin,xyzmax);
1395 
1396  } else if ( stype.find("zpoly") != string::npos ) {
1397  // polygon along z direction nfaces at (x0,y0) radius phi zmin zmax
1398  if ( nvals < 7 )
1399  LOG("gevgen_fnal", pFATAL) << "MakeZPolygon needs 7 values, not " << nvals
1400  << " fidcut=\"" << fidcut << "\"";
1401  int nfaces = (int)vals[0];
1402  if ( nfaces < 3 )
1403  LOG("gevgen_fnal", pFATAL) << "MakeZPolygon needs nfaces>=3, not " << nfaces
1404  << " fidcut=\"" << fidcut << "\"";
1405  fidsel->MakeZPolygon(nfaces,vals[1],vals[2],vals[3],vals[4],vals[5],vals[6]);
1406 
1407  } else if ( stype.find("sphere") != string::npos ) {
1408  // sphere at (x0,y0,z0) radius
1409  if ( nvals < 4 )
1410  LOG("gevgen_fnal", pFATAL) << "MakeZSphere needs 4 values, not " << nvals
1411  << " fidcut=\"" << fidcut << "\"";
1412  fidsel->MakeSphere(vals[0],vals[1],vals[2],vals[3]);
1413 
1414  } else {
1415  LOG("gevgen_fnal", pFATAL)
1416  << "Can not create GeomVolSelectorFiduction for shape \"" << stype << "\"";
1417  }
1418 
1419  if ( master ) {
1420  fidsel->ConvertShapeMaster2Top(rgeom);
1421  LOG("gevgen_fnal", pNOTICE) << "Convert fiducial volume from master to topvol coords";
1422  }
1423  if ( reverse ) {
1424  fidsel->SetReverseFiducial(true);
1425  LOG("gevgen_fnal", pNOTICE) << "Reverse sense of fiducial volume cut";
1426  }
1427  rgeom->AdoptGeomVolSelector(fidsel);
1428 
1429 }
1430 //____________________________________________________________________________
1431 void CreateRockBoxSelection (string fidcut, GeomAnalyzerI* geom_driver)
1432 {
1433 
1434  if( fidcut.find_first_not_of(" \t\n") != 0) // trim any leading whitespace
1435  fidcut.erase( 0, fidcut.find_first_not_of(" \t\n") );
1436 
1437  // convert string to lowercase
1438  std::transform(fidcut.begin(),fidcut.end(),fidcut.begin(),::tolower);
1439 
1441  dynamic_cast<genie::geometry::ROOTGeomAnalyzer *>(geom_driver);
1442  if ( ! rgeom ) {
1443  LOG("gevgen_fnal", pWARN)
1444  << "Can not create GeomVolSelectorRockBox,"
1445  << " geometry driver is not ROOTGeomAnalyzer";
1446  return;
1447  }
1448 
1449  LOG("gevgen_fnal", pWARN) << "fiducial (rock) cut: " << fidcut;
1450 
1451  // for now, only fiducial no "rock box"
1454 
1455  vector<string> strtok = genie::utils::str::Split(fidcut,":");
1456  if ( strtok.size() != 2 ) {
1457  LOG("gevgen_fnal", pWARN)
1458  << "Can not create GeomVolSelectorRockBox,"
1459  << " no \":\" separating type from values. nsplit=" << strtok.size();
1460  for ( unsigned int i=0; i < strtok.size(); ++i )
1461  LOG("gevgen_fnal", pWARN)
1462  << "strtok[" << i << "] = \"" << strtok[i] << "\"";
1463  return;
1464  }
1465 
1466  string stype = strtok[0];
1467 
1468  // parse out values
1469  vector<double> vals;
1470  vector<string> valstrs = genie::utils::str::Split(strtok[1]," ,;(){}[]\t\n\r");
1471  vector<string>::const_iterator iter = valstrs.begin();
1472  for ( ; iter != valstrs.end(); ++iter ) {
1473  const string& valstr1 = *iter;
1474  if ( valstr1 != "" ) {
1475  double aval = atof(valstr1.c_str());
1476  LOG("gevgen_fnal", pWARN) << "rock value [" << vals.size() << "] "
1477  << aval;
1478  vals.push_back(aval);
1479  }
1480  }
1481  size_t nvals = vals.size();
1482 
1483  rocksel->SetRemoveEntries(true); // drop segments that won't be considered
1484 
1485  // assume coordinates are in the *master* (not "top volume") system
1486  // need to set fTopVolume to fWorldVolume
1487  //fTopVolume = fWorldVolume;
1488  //rgeom->SetTopVolName(fTopVolume.c_str());
1490  rgeom->SetTopVolName(gOptRootGeomMasterVol);
1491 
1492  if ( nvals < 6 ) {
1493  LOG("gevgen_fnal", pFATAL) << "rockbox needs at "
1494  << "least 6 values, found "
1495  << nvals << "in \""
1496  << strtok[1] << "\"";
1497  exit(1);
1498 
1499  }
1500  double xyzmin[3] = { vals[0], vals[1], vals[2] };
1501  double xyzmax[3] = { vals[3], vals[4], vals[5] };
1502 
1503  bool rockonly = true;
1504  double wallmin = 800.; // geometry in cm, ( 8 meter buffer)
1505  double dedx = 2.5 * 1.7e-3; // GeV/cm, rho=2.5, 1.7e-3 ~ rock like loss
1506  double fudge = 1.05;
1507 
1508  if ( nvals >= 7 ) rockonly = vals[6];
1509  if ( nvals >= 8 ) wallmin = vals[7];
1510  if ( nvals >= 9 ) dedx = vals[8];
1511  if ( nvals >= 10 ) fudge = vals[9];
1512 
1513  rocksel->SetRockBoxMinimal(xyzmin,xyzmax);
1514  rocksel->SetMinimumWall(wallmin);
1515  rocksel->SetDeDx(dedx/fudge);
1516 
1517  if ( nvals >= 11 ) rocksel->SetExpandFromInclusion((int)vals[10]);
1518 
1519  // if not rock-only then make a tiny exclusion bubble
1520  // call to MakeBox shouldn't be necessary
1521  // should be done by SetRockBoxMinimal but for some GENIE versions isn't
1522  if ( ! rockonly ) rocksel->MakeSphere(0,0,0,1.0e-10);
1523  else rocksel->MakeBox(xyzmin,xyzmax);
1524 
1525  rgeom->AdoptGeomVolSelector(rocksel);
1526 
1527 }
1528 
1529 //____________________________________________________________________________
1530 void DetermineFluxDriver(string fopt)
1531 {
1532  // based on the -f option string determine which flux driver to use
1533  // this may take some guessing
1534 
1535  // first look for strings that look like "<proto>:..."
1536  // or ....<proto>.root,....
1537  // where "<proto>" is a key the gOptFluxShortNames map
1538 
1539  map<string,string>::const_iterator mitr = gOptFluxShortNames.begin();
1540  map<string,string>::const_iterator mitr_end = gOptFluxShortNames.end();
1541  for ( ; mitr != mitr_end; ++mitr ) {
1542  string proto = mitr->first + string(":");
1543  string gproto = string("g") + proto;
1544  string protor = proto + ".root,";
1545  string full = mitr->second;
1546  if ( fopt.find(proto) == 0 ) {
1547  fopt.erase(0,proto.size());
1548  gOptFluxDriver = full;
1549  break;
1550  } else if ( fopt.find(gproto) == 0 ) {
1551  fopt.erase(0,gproto.size());
1552  gOptFluxDriver = full;
1553  break;
1554  } else if ( fopt.find(protor) != std::string::npos ) {
1555  gOptFluxDriver = full;
1556  break;
1557  }
1558  }
1559  // tested all cases where user might have specified explicitly
1560  // or been part of an extended file extension
1561  // this is where it gets messy
1562  if ( gOptFluxDriver == "" ) {
1563 
1564  // not specified ? guess from file name itself
1565  if ( fopt.find("gsimple") != std::string::npos ) {
1566  // put dk2nu after gsimple in case simple files are derived from dk2nu
1567  // then both are in name we should choose gsimple
1568  gOptFluxDriver = "genie::flux::GSimpleNtpFlux";
1569  } else if ( fopt.find("dk2nu") != std::string::npos ) {
1570  gOptFluxDriver = "genie::flux::GDk2NuFlux";
1571  } else {
1572  // does it look like the histogram format
1573  const char* hstrings[] = { ",12[", ",+12[", ",-12[",
1574  ",14[", ",+14[", ",-14[",
1575  ",16[", ",+16[", ",-16[" };
1576  size_t nh = sizeof(hstrings)/sizeof(const char*);
1577  for (size_t ih=0; ih<nh; ++ih) {
1578  if ( fopt.find(hstrings[ih]) != std::string::npos ) {
1579  // hey!
1580  gOptFluxDriver = "genie::flux::GCylindTH1Flux";
1581  break;
1582  }
1583  } // loop over possible histogram specifiers
1584  }
1585 
1586  // fall through default ... hope it works
1587  if ( gOptFluxDriver == "" ) {
1588  gOptFluxDriver = "genie::flux::GNuMIFlux";
1589  }
1590  }
1591 
1592  gOptUsingHistFlux = ( gOptFluxDriver == "genie::flux::GCylindTH1Flux" );
1593  if ( gOptUsingHistFlux ) ParseFluxHst(fopt);
1594  else ParseFluxFileConfig(fopt);
1595 }
1596 //____________________________________________________________________________
1597 void ParseFluxHst(string flux)
1598 {
1599  // Using flux from histograms
1600  // Extract the root file name & the list of histogram names & neutrino
1601  // species (specified as 'filename,histo1[species1],histo2[species2],...')
1602  // See documentation on top section of this file.
1603  //
1604 
1605  vector<string> fluxv = utils::str::Split(flux,",");
1606  if(fluxv.size()<2) {
1607  LOG("gevgen_fnal", pFATAL)
1608  << "You need to specify both a flux ntuple ROOT file "
1609  << " _AND_ a detector location";
1610  PrintSyntax();
1611  exit(1);
1612  }
1613  gOptFluxFile = fluxv[0];
1614  bool accessible_flux_file = !(gSystem->AccessPathName(gOptFluxFile.c_str()));
1615  if (!accessible_flux_file) {
1616  LOG("gevgen_fnal", pFATAL)
1617  << "Can not access flux file: " << gOptFluxFile;
1618  PrintSyntax();
1619  exit(1);
1620  }
1621  // Extract energy spectra for all specified neutrino species
1622  TFile flux_file(gOptFluxFile.c_str(), "read");
1623  for(unsigned int inu=1; inu<fluxv.size(); inu++) {
1624  string nutype_and_histo = fluxv[inu];
1625  string::size_type open_bracket = nutype_and_histo.find("[");
1626  string::size_type close_bracket = nutype_and_histo.find("]");
1627  if (open_bracket ==string::npos ||
1628  close_bracket==string::npos)
1629  {
1630  LOG("gevgen_fnal", pFATAL)
1631  << "You made an error in specifying the flux histograms";
1632  PrintSyntax();
1633  exit(1);
1634  }
1635  string::size_type ibeg = 0;
1636  string::size_type iend = open_bracket;
1637  string::size_type jbeg = open_bracket+1;
1638  string::size_type jend = close_bracket;
1639  string nutype = nutype_and_histo.substr(ibeg,iend-ibeg);
1640  string histo = nutype_and_histo.substr(jbeg,jend-jbeg);
1641  // access specified histogram from the input root file
1642  TH1D * ihst = (TH1D*) flux_file.Get(histo.c_str());
1643  if(!ihst) {
1644  LOG("gevgen_fnal", pFATAL)
1645  << "Can not find histogram: " << histo
1646  << " in flux file: " << gOptFluxFile;
1647  PrintSyntax();
1648  exit(1);
1649  }
1650  // create a local copy of the input histogram
1651  TH1D * spectrum = new TH1D(
1652  histo.c_str(), histo.c_str(), ihst->GetNbinsX(),
1653  ihst->GetXaxis()->GetXmin(), ihst->GetXaxis()->GetXmax());
1654  spectrum->SetDirectory(0);
1655  for(int ibin = 1; ibin <= ihst->GetNbinsX(); ibin++) {
1656  spectrum->SetBinContent(ibin, ihst->GetBinContent(ibin));
1657  }
1658  // convert neutrino name -> pdg code
1659  int pdg = atoi(nutype.c_str());
1660  if(!pdg::IsNeutrino(pdg) && !pdg::IsAntiNeutrino(pdg)) {
1661  LOG("gevgen_fnal", pFATAL)
1662  << "Unknown neutrino type: " << nutype;
1663  PrintSyntax();
1664  exit(1);
1665  }
1666  // store flux neutrino code / energy spectrum
1667  LOG("gevgen_fnal", pDEBUG)
1668  << "Adding energy spectrum for flux neutrino: pdg = " << pdg;
1669  gOptFluxHst.insert(map<int, TH1D*>::value_type(pdg, spectrum));
1670  }//inu
1671 
1672  if(gOptFluxHst.size()<1) {
1673  LOG("gevgen_fnal", pFATAL)
1674  << "You have not specified any flux histogram!";
1675  PrintSyntax();
1676  exit(1);
1677  }
1678 
1679  flux_file.Close();
1680 }
1681 
1682 //____________________________________________________________________________
1684 {
1685  // Using gnumi/gsimple/dk2nu beam flux ntuples
1686  // Extract beam flux (root) file name & detector location
1687  //
1688  vector<string> fluxv = utils::str::Split(flux,",");
1689  if(fluxv.size()<2) {
1690  LOG("gevgen_fnal", pFATAL)
1691  << "You need to specify both a flux ntuple ROOT file "
1692  << " _AND_ a detector location";
1693  PrintSyntax();
1694  exit(1);
1695  }
1696  gOptFluxFile = fluxv[0];
1697  gOptDetectorLocation = fluxv[1];
1698 
1699  for ( size_t j = 2; j < fluxv.size(); ++j ) {
1700  int ipdg = atoi(fluxv[j].c_str());
1701  gOptFluxPdg.push_back(ipdg);
1702  }
1703 
1704 }
1705 
1706 //____________________________________________________________________________
string gOptFluxFile
void MakeSphere(Double_t x0, Double_t y0, Double_t z0, Double_t radius)
void RandGen(long int seed)
Definition: AppInit.cxx:31
static void SetPrintLevel(int print_level)
Definition: GHepRecord.cxx:992
void split(double tt, double *fr)
const char * GetExposureUnits() const
what units are returned by GetTotalExposure?
GENIE Interface for user-defined volume selector functors Trim path segments based on the intersectio...
const XML_Char * name
Definition: expat.h:151
A class for generating concrete GFluxI derived classes based on the factory pattern. This code supplies a CPP macro which allows the classes to self-register and thus no modification of this class is needed in order to expand the list of classes it knows about.
void XSecTable(string inpfile, bool require_table)
Definition: AppInit.cxx:39
virtual GeomVolSelectorI * AdoptGeomVolSelector(GeomVolSelectorI *selector)
configure processing to perform path segment trimming
long ArgAsLong(char opt)
PDGCodeList gOptFluxPdg
double ArgAsDouble(char opt)
void ParseFluxFileConfig(string fopt)
string kDefOptEvFilePrefix
virtual void LoadBeamSimData(const std::vector< std::string > &filenames, const std::string &det_loc)=0
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:108
string ArgAsString(char opt)
set< int >::iterator it
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
void LoadExtraOptions(void)
void MakeZCylinder(Double_t x0, Double_t y0, Double_t radius, Double_t zmin, Double_t zmax)
void SetEventGeneratorList(string listname)
Definition: GMCJDriver.cxx:157
void GetCommandLineArgs(int argc, char **argv)
virtual void PrintConfig()=0
print the current configuration
bool gOptUsingRootGeom
void ReadFromCommandLine(int argc, char **argv)
Definition: RunOpt.cxx:107
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
string gOptDetectorLocation
string gOptFidCut
virtual const PathLengthList & GetMaxPathLengths(void) const
void CreateRockBoxSelection(string fidcut, GeomAnalyzerI *geom_driver)
virtual void ConvertShapeMaster2Top(const ROOTGeomAnalyzer *rgeom)
void Update(int iev, const EventRecord *event)
Definition: GMCJMonitor.cxx:58
OStream cerr
Definition: OStream.cxx:7
double gOptGeomDUnits
const std::vector< std::string > & AvailableFluxDrivers() const
bool gSigTERM
int gOptNev
static void gsSIGTERMhandler(int)
virtual void SetNumOfCycles(long int ncycle)
limit cycling through input files
string gOptRootGeomMasterVol
virtual void SetScannerNParticles(int np)
bool gOptWriteMaxPlXml
void DetermineFluxDriver(string fopt)
virtual long int NFluxNeutrinos() const =0
of rays generated
Loaders::FluxType flux
long int gOptRanSeed
int gOptNScan
A list of PDG codes.
Definition: PDGCodeList.h:33
NtpMCFormat_t kDefOptNtpFormat
virtual void SetUpstreamZ(double z0)
void SetRefreshRate(int rate)
Definition: GMCJMonitor.cxx:53
void UseFluxDriver(GFluxI *flux)
Definition: GMCJDriver.cxx:174
Object to be filled with the neutrino path-length, for all detector geometry materials, when starting from a position x and travelling along the direction of the neutrino 4-momentum.
const XML_Char * s
Definition: expat.h:262
unsigned int nutype(unsigned int u)
Definition: runWimpSim.h:122
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
void SaveAsXml(string filename) const
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Definition: typedefs.hpp:11
TTree * EventTree(void)
Definition: NtpWriter.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
string kDefOptGeomDUnits
Definition: lz4.cxx:387
void ForceSingleProbScale(void)
Definition: GMCJDriver.cxx:219
void ParseFluxHst(string fopt)
virtual void SetScannerNRays(int nr)
A generic GENIE flux driver. Generates a &#39;cylindrical&#39; neutrino beam along the input direction...
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:116
bool UseMaxPathLengths(string xml_filename)
Definition: GMCJDriver.cxx:190
double UnitFromString(string u)
Definition: UnitUtils.cxx:26
double gOptGeomLUnits
GENIE interface for uniform flux exposure iterface.
GENIE Interface for limiting vertex selection in the rock to a volume that depends (in part) on the n...
#define pot
void Configure(bool calc_prob_scales=true)
Definition: GMCJDriver.cxx:446
void MakeZPolygon(Int_t n, Double_t x0, Double_t y0, Double_t inradius, Double_t phi0deg, Double_t zmin, Double_t zmax)
A ROOT/GEANT4 geometry driver.
map< string, string > gOptFluxShortNames
const double j
Definition: BetheBloch.cxx:29
void Save(void)
get the even tree
Definition: NtpWriter.cxx:214
TH2D * histo
#define pINFO
Definition: Messenger.h:63
map< int, double > gOptTgtMix
string gOptRootGeomTopVol
double t2
string gOptExtMaxPlXml
void AddEventRecord(int ievent, const EventRecord *ev_rec)
save the event tree
Definition: NtpWriter.cxx:69
void SetTransverseRadius(double Rt)
const ana::Var wgt
void SetNuDirection(const TVector3 &direction)
void BuildTune()
build tune and inform XSecSplineList
Definition: RunOpt.cxx:100
#define pWARN
Definition: Messenger.h:61
EventRecord * GenerateEvent(void)
Definition: GMCJDriver.cxx:838
OStream cout
Definition: OStream.cxx:6
virtual void GetBranchInfo(std::vector< std::string > &branchNames, std::vector< std::string > &branchClassNames, std::vector< void ** > &branchObjPointers)
GENIE interface for uniform flux exposure iterface.
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
Long_t gOptRunNu
void CustomizeFilenamePrefix(string prefix)
Definition: NtpWriter.cxx:123
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void PrintSyntax(void)
string gOptInpXSecFile
void Initialize(void)
add event
Definition: NtpWriter.cxx:95
int main(int argc, char **argv)
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
static RunOpt * Instance(void)
Definition: RunOpt.cxx:62
vector< string > Split(string input, string delim)
Definition: StringUtils.cxx:42
int gOptDebug
A utility class to facilitate creating the GENIE MC Ntuple from the output GENIE GHEP event records...
Definition: NtpWriter.h:40
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
string gOptEvFilePrefix
string gOptRootGeom
virtual void SetScannerNPoints(int np)
set geometry driver&#39;s configuration options
virtual void SetScannerFlux(GFluxI *f)
app
Definition: demo.py:189
void geom(int which=0)
Definition: geom.C:163
virtual double GetTotalExposure() const =0
static GFluxDriverFactory & Instance()
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
exit(0)
assert(nhit_max >=nhit_nbins)
virtual bool End(void)=0
true if no more flux nu&#39;s can be thrown (eg reaching end of beam sim ntuples)
string PrintFramedMesg(string mesg, unsigned int nl=1, const char f='*')
Definition: PrintUtils.cxx:171
A vector of EventGeneratorI objects.
virtual void SetFluxParticles(const PDGCodeList &particles)
specify list of flux neutrino species
string kDefOptGeomLUnits
void SetBeamSpot(const TVector3 &spot)
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:100
void AddEnergySpectrum(int nu_pdgc, TH1D *spectrum)
Command line argument parser.
enum genie::ENtpMCFormat NtpMCFormat_t
#define pNOTICE
Definition: Messenger.h:62
void SetRockBoxMinimal(Double_t *xyzmin, Double_t *xyzmax)
void MakeBox(Double_t *xyzmin, Double_t *xyzmax)
double gOptPOT
Defines the GENIE Geometry Analyzer Interface.
Definition: GeomAnalyzerI.h:30
Float_t e
Definition: plot.C:35
void CreateFidSelection(string fidcut, GeomAnalyzerI *geom_driver)
void RecursiveExhaust(TGeoVolume *topvol, string volnames, bool exhaust)
Definition: GeoUtils.cxx:26
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
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...
bool OptionExists(char opt)
was option set?
void CacheFile(string inpfile)
Definition: AppInit.cxx:118
double gOptZmin
void EnableBareXSecPreCalc(bool flag)
Definition: RunOpt.h:61
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
map< int, TH1D * > gOptFluxHst
virtual TGeoManager * GetGeometry(void) const
genie::GFluxI * GetFluxDriver(const std::string &)
#define pDEBUG
Definition: Messenger.h:64
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
string gOptFluxDriver
bool gOptUsingHistFlux