gEvGenLArDM.cxx
Go to the documentation of this file.
1 //______________________________________________________________________________
2 /*!
3 
4 \program gevgen_lardm
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  *** Synopsis :
11 
12  gevgen_lardm [-h]
13  [-r run#]
14  -f flux
15  -g geometry
16  -M dm mass
17  [-c zp_coupling]
18  [-v med_ratio]
19  [-t top_volume_name_at_geom || -t +Vol1-Vol2...]
20  [-m max_path_lengths_xml_file]
21  [-L length_units_at_geom]
22  [-D density_units_at_geom]
23  [-n n_of_events]
24  [-o output_event_file_prefix]
25  [-F fid_cut_string]
26  [-S nrays]
27  [-z zmin]
28  [-d debug flags]
29  [--seed random_number_seed]
30  [ --cross-sections xml_file]
31  [--event-generator-list list_name]
32  [--tune genie_tune]
33  [--message-thresholds xml_file]
34  [--unphysical-event-mask mask]
35  [--event-record-print-level level]
36  [--mc-job-status-refresh-rate rate]
37  [--cache-file root_file]
38 
39  *** Options :
40 
41  [] Denotes an optional argument
42 
43  -h
44  Prints out the gevgen_lardm syntax and exits.
45  -r
46  Specifies the MC run number [default: 1000]
47  -g
48  Input 'geometry'.
49  This option can be used to specify any of:
50  1 > A ROOT file containing a ROOT/GEANT geometry description
51  [Examples]
52  - To use the master volume from the ROOT geometry stored
53  in the /some/path/nova-geom.root file, type:
54  '-g /some/path/nova-geom.root'
55  2 > A mix of target materials, each with its corresponding weight,
56  typed as a comma-separated list of nuclear PDG codes (in the
57  std PDG2006 convention: 10LZZZAAAI) with the weight fractions
58  in brackets, eg code1[fraction1],code2[fraction2],...
59  If that option is used (no detailed input geometry description)
60  then the interaction vertices are distributed in the detector
61  by the detector MC.
62  [Examples]
63  - To use a target mix of 95% O16 and 5% H type:
64  '-g 1000080160[0.95],1000010010[0.05]'
65  - To use a target which is 100% C12, type:
66  '-g 1000060120'
67  -M
68  Specifies the DM mass
69  -c
70  Specifies the Z' coupling constant
71  Default: Value in UserPhysicsOptions.xml
72  -t
73  Input 'top volume' for event generation -
74  can be used to force event generation in given sub-detector.
75  [default: the 'master volume' of the input geometry]
76  You can also use the -t option to switch generation on/off at
77  multiple volumes as, for example, in:
78  `-t +Vol1-Vol2+Vol3-Vol4',
79  `-t "+Vol1 -Vol2 +Vol3 -Vol4"',
80  `-t -Vol2-Vol4+Vol1+Vol3',
81  `-t "-Vol2 -Vol4 +Vol1 +Vol3"'m
82  where:
83  "+Vol1" and "+Vol3" tells GENIE to `switch on' Vol1 and Vol3, while
84  "-Vol2" and "-Vol4" tells GENIE to `switch off' Vol2 and Vol4.
85  If the very first character is a '+', GENIE will neglect all volumes
86  except the ones explicitly turned on. Vice versa, if the very first
87  character is a `-', GENIE will keep all volumes except the ones
88  explicitly turned off (feature contributed by J.Holeczek).
89  -m
90  An XML file (generated by gmxpl) with the max (density weighted)
91  path-lengths for each target material in the input ROOT geometry.
92  If no file is input, then the geometry will be scanned at MC job
93  initialization to determine those max path lengths.
94  Supplying this file can speed-up the MC job initialization.
95  -L
96  Input geometry length units, eg "m", "cm", "mm", ...
97  [default: mm]
98  -D
99  Input geometry density units, eg "g_cm3", "clhep_def_density_unit",...
100  [default: g_cm3]
101  -f
102  Input 'neutrino flux'.
103  This option can be used to specify any of:
104  1 > A g[4]numi[_flugg] beam simulation output file
105  and the detector location
106  The general sytax is:
107  -f /full/path/flux_file.root,detector,flavor1,flavor2...
108  [Notes]
109  - For more information on the flux ntuples, see the gNuMI doc.
110  - The original hbook ntuples need to be converted to a ROOT
111  format using the h2root ROOT utility.
112  - If flavors aren't specified then use default (12,-12,14,-14)
113  - See GNuMIFlux.xml or Dk2Nu.xml for all supported detector
114  locations
115  - The g[4]NuMI[_flugg] flux ntuples are read via GENIE's
116  GNuMIFlux driver, and dk2nu file via an external
117  product w/ the GDk2NuFlux driver (if it can be found).
118  This customized GENIE event generation driver passes-through
119  the complete gNuMI input flux information (eg parent decay
120  kinematics / position etc) for each neutrino event it
121  generates (as an additional 'flux' branch is added on the
122  output event tree).
123  [Examples]
124  - To use the gNuMI flux ntuple flux.root at MINOS near
125  detector location, type:
126  '-f /path/flux.root,MINOS-NearDet'
127  1a> Similar to 1 above, but filename contains "dk2nu", then use
128  the GDk2NuFlux driver
129  1b> Similar to 1 above, but filename contains "gsimple", then
130  use GSimpleNtpFlux driver
131  2 > A set of histograms stored in a ROOT file.
132  The general syntax is:
133  -f /path/histogram_file.root,neutrino_code[histo_name],...
134  [Notes]
135  - The neutrino codes are the PDG ones.
136  - The 'neutrino_code[histogram_name]' part of the option can
137  be repeated multiple times (separated by commas), once for
138  each flux neutrino species you want to consider, eg
139  '-f somefile.root,12[nuehst],-12[nuebarhst],14[numuhst]'
140  - When using flux from histograms then there is no point in
141  using a 'detailed detector geometry description' as your
142  flux input contains no directional information for those
143  flux neutrinos.
144  The neutrino direction is conventionally set to be
145  +z {x=0,y=0}.
146  So, when using this option you must be using a simple
147  'target mix'
148  See the -g option for possible geometry settings.
149  If you want to use the detailed detector geometry
150  description then you should be using a driver that
151  supplies a detailed simulated beam flux.
152  - When using flux from histograms there is no branch with
153  neutrino parent information added in the output event
154  tree as your flux input contains no such information.
155  - Note that the relative normalization of the flux histograms
156  is taken into account and is reflected in the relative
157  frequency of flux neutrinos thrown by the flux driver
158  [Examples]
159  - To use the histogram 'h100' (representing the nu_mu flux)
160  and the histogram 'h300' (representing the nu_e flux) and
161  the histogram 'h301' (representing the nu_e_bar flux) from
162  the flux.root file in /path/ type:
163  '-f /path/flux.root,14[h100],12[h300],-12[h301]
164 
165  -n
166  Specifies how many events to generate.
167 
168  -F
169  Apply a fiducial cut (for now hard coded ... generalize)
170  Only used with ROOTGeomAnalyzer
171  if string starts with "-" then reverses sense (ie. anti-fiducial)
172  -S
173  Number of rays to use to scan geometry for max path length
174  Only used with ROOTGeomAnalyzer & GNuMIFlux
175  +N Use flux to scan geometry for max path length
176  -N Use N rays x N points on each face of a box
177  -z
178  Z from which to start flux ray in user world coordinates
179  Only use with ROOTGeomAnalyzer & GNuMIFlux
180  If left unset then flux originates on the flux window
181  [No longer attempts to determine z from geometry, generally
182  got this wrong]
183  -o
184  Sets the prefix of the output event file.
185  The output filename is built as:
186  [prefix].[run_number].[event_tree_format].[file_format]
187  The default output filename is:
188  gntp.[run_number].ghep.root
189  This cmd line arguments lets you override 'gntp'
190  --seed
191  Random number seed.
192  --cross-sections
193  Name (incl. full path) of an XML file with pre-computed
194  cross-section values used for constructing splines.
195  --tune
196  Specifies a GENIE comprehensive neutrino interaction model tune.
197  [default: "Default"].
198  --message-thresholds
199  Allows users to customize the message stream thresholds.
200  The thresholds are specified using an XML file.
201  See $GENIE/config/Messenger.xml for the XML schema.
202  --unphysical-event-mask
203  Allows users to specify a 16-bit mask to allow certain types of
204  unphysical events to be written in the output file.
205  [default: all unphysical events are rejected]
206  --event-record-print-level
207  Allows users to set the level of information shown when the event
208  record is printed in the screen. See GHepRecord::Print().
209  --mc-job-status-refresh-rate
210  Allows users to customize the refresh rate of the status file.
211  --cache-file
212  Allows users to specify a cache file so that the cache can be
213  re-used in subsequent MC jobs.
214 
215  *** Examples:
216 
217  (1) shell% gevgen_lardm
218  -r 1001
219  -f /data/mc_inputs/flux/flux_00001.root,MINOS-NearDet,12,-12
220  -g /data/mc_inputs/geom/minos.root
221  -L mm -D g_cm3
222  -e 5E+17
223  --cross-sections /data/xsec.xml
224 
225  Generate events (run number 1001) using the gNuMI flux ntuple in
226  /data/mc_inputs/flux/v1/flux_00001.root
227  The job will load the MINOS near detector detector geometry
228  description from /data/mc_inputs/geom/minos.root and interpret it
229  using 'mm' as the length unit and 'g/cm^3' as the density unit.
230  The job will stop as it accumulates a sample corresponding to
231  5E+17 POT.
232  Pre-computed cross-section data are loaded from /data/xsec.xml
233 
234  (2) shell% gevgen_lardm
235  -r 1001
236  -f /data/t2k/flux/hst/flux.root,12[h100],-12[h101],14[h200]
237  -g 1000080160[0.95],1000010010[0.05]
238  -n 50000
239  --cross-sections /data/xsec.xml
240 
241  Please read the GENIE user manual for more information.
242 
243 \author Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
244  University of Liverpool & STFC Rutherford Appleton Lab
245 
246  Robert Hatcher <rhatcher \at fnal.gov>
247  Fermi National Accelerator Laboratory
248 
249 \created August 20, 2008
250 
251 \cpright Copyright (c) 2003-2019, The GENIE Collaboration
252  For the full text of the license visit http://copyright.genie-mc.org
253  or see $GENIE/LICENSE
254 */
255 //_________________________________________________________________________________________
256 
257 #include <cassert>
258 #include <cstdlib>
259 #include <csignal>
260 
261 #include <string>
262 #include <sstream>
263 #include <vector>
264 #include <map>
265 #include <algorithm> // for transform()
266 #include <fstream>
267 
268 #include <TSystem.h>
269 #include <TError.h> // for gErrorIgnoreLevel
270 #include <TTree.h>
271 #include <TFile.h>
272 #include <TH1D.h>
273 #include <TMath.h>
274 #include <TGeoVolume.h>
275 #include <TGeoShape.h>
276 
293 #include "Framework/Utils/AppInit.h"
294 #include "Framework/Utils/RunOpt.h"
298 
299 #ifdef __GENIE_FLUX_DRIVERS_ENABLED__
303 
304 //#include "Tools/Flux/GNuMIFlux.h"
305 //#include "Tools/Flux/GSimpleNtpFlux.h"
306 // #ifdef __DK2NU_FLUX_DRIVER_AVAILABLE__
307 // #include "dk2nu/tree/dk2nu.h"
308 // #include "dk2nu/tree/dkmeta.h"
309 // #include "dk2nu/tree/NuChoice.h"
310 // #include "dk2nu/genie/GDk2NuFlux.h"
311 // #endif
312 
313 #endif
314 
315 #ifdef __GENIE_GEOM_DRIVERS_ENABLED__
316 #include "Tools/Geometry/GeoUtils.h"
321 #endif
322 
323 using std::string;
324 using std::vector;
325 using std::map;
326 using std::ostringstream;
327 
328 using namespace genie;
329 
330 // Forward function declarations
331 //
332 void LoadExtraOptions (void);
333 void GetCommandLineArgs (int argc, char ** argv);
334 void PrintSyntax (void);
335 void CreateFidSelection (string fidcut, GeomAnalyzerI* geom_driver);
336 void CreateRockBoxSelection (string fidcut, GeomAnalyzerI* geom_driver);
337 void DetermineFluxDriver(string fopt);
338 void ParseFluxHst (string fopt);
339 void ParseFluxFileConfig(string fopt);
340 
341 // Default options (override them using the command line arguments):
342 //
343 string kDefOptGeomLUnits = "mm"; // default geometry length units
344 string kDefOptGeomDUnits = "g_cm3"; // default geometry density units
345 NtpMCFormat_t kDefOptNtpFormat = kNFGHEP; // default event tree format
346 string kDefOptEvFilePrefix = "gntp";
347 
348 // User-specified options:
349 //
350 double gOptZpCoupling; // mediator coupling
351 double gOptDMMass; // DM Mass
352 double gOptMedRatio; // ratio of mediator to DM mass
353 Long_t gOptRunNu; // run number
354 bool gOptUsingRootGeom = false; // using root geom or target mix?
355 map<int,double> gOptTgtMix; // target mix (tgt pdg -> wght frac) / if not using detailed root geom
356 string gOptRootGeom; // input ROOT file with realistic detector geometry
357 string gOptRootGeomTopVol = ""; // input geometry top event generation volume
358 string gOptRootGeomMasterVol = ""; // (highest level of geometry)
359 double gOptGeomLUnits = 0; // input geometry length units
360 double gOptGeomDUnits = 0; // input geometry density units
361 string gOptExtMaxPlXml = ""; // max path lengths XML file for input geometry
362 bool gOptWriteMaxPlXml = false; // rather than read file, write the file
363  // triggered by leading '+' on given ofilename
364 string gOptFluxFile; // ROOT file with gnumi beam flux ntuple
365 int gOptNev; // number of events to generate
366 string gOptFidCut; // fiducial cut selection
367 int gOptNScan = 0; // # of geometry scan rays
368 double gOptZmin = -2.0e30; // starting z position [ if abs() < 1e30 ]
369 string gOptEvFilePrefix; // event file prefix
370 int gOptDebug = 0; // debug flags
371 long int gOptRanSeed; // random number seed
372 string gOptInpXSecFile; // cross-section splines
373 
374 bool gSigTERM = false; // was TERM signal sent?
375 
376 static void gsSIGTERMhandler(int /* s */)
377 {
378  gSigTERM = true;
379  std::cerr << "Caught SIGTERM" << std::endl;
380 }
381 
382 //____________________________________________________________________________
383 int main(int argc, char ** argv)
384 {
386  GetCommandLineArgs(argc,argv);
388  if (gOptZpCoupling > 0.) {
389  Registry * r = AlgConfigPool::Instance()->CommonList("Param", "BoostedDarkMatter");
390  r->UnLock();
391  r->Set("ZpCoupling", gOptZpCoupling);
392  r->Lock();
393  }
394 
395  if ( ! RunOpt::Instance()->Tune() ) {
396  LOG("gmkspl", pFATAL) << " No TuneId in RunOption";
397  exit(-1);
398  }
400 
401  // Initialization of random number generators, cross-section table,
402  // messenger thresholds, cache file
403  utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles());
407 
408  // Set GHEP print level
409  GHepRecord::SetPrintLevel(RunOpt::Instance()->EventRecordPrintLevel());
410 
411  // *************************************************************************
412  // * Create / configure the geometry driver
413  // *************************************************************************
414  GeomAnalyzerI * geom_driver = 0;
415 
416  if(gOptUsingRootGeom) {
417  //
418  // *** Using a realistic root-based detector geometry description
419  //
420 
421  // creating & configuring a root geometry driver
424  TGeoVolume * topvol = rgeom->GetGeometry()->GetTopVolume();
425  if ( ! topvol ) {
426  LOG("gevgen_numi", pFATAL) << "Null top ROOT geometry volume!";
427  exit(1);
428  }
429  // retrieve the master volume name
430  gOptRootGeomMasterVol = topvol->GetName();
431 
432  rgeom -> SetLengthUnits (gOptGeomLUnits);
433  rgeom -> SetDensityUnits (gOptGeomDUnits);
434  rgeom -> SetTopVolName (gOptRootGeomTopVol); // set user defined "topvol"
435 
436  // getting the bounding box dimensions along z so as to set the
437  // appropriate upstream generation surface for the NuMI flux driver
438 
439  // RWH 2010-07-16: do not try to automatically get zmin from geometry, rather
440  // by default let the flux start from the window. If the user wants to
441  // override this then they need to explicitly set a "zmin". Trying to use
442  // the geometry is fraught with problems in local vs. global coordinates and
443  // units where it can appear to work in some cases but it actually isn't really
444  // universally correct.
445  //was// TGeoShape * bounding_box = topvol->GetShape();
446  //was// bounding_box->GetAxisRange(3, zmin, zmax);
447  //was// zmin *= rgeom->LengthUnits();
448  //was// zmax *= rgeom->LengthUnits();
449 
450  // switch on/off volumes as requested
451  if ( (gOptRootGeomTopVol[0] == '+') || (gOptRootGeomTopVol[0] == '-') ) {
452  bool exhaust = (*gOptRootGeomTopVol.c_str() == '+');
454  }
455 
456  // casting to the GENIE geometry driver interface
457  geom_driver = dynamic_cast<GeomAnalyzerI *> (rgeom);
458 
459  // user specifid a fiducial volume cut ... parse that out
460  if ( gOptFidCut.find("rock") != std::string::npos )
462  else if ( gOptFidCut != "" ) CreateFidSelection(gOptFidCut,rgeom);
463 
464  }
465  else {
466  //
467  // *** Using a 'point' geometry with the specified target mix
468  // *** ( = a list of targets with their corresponding weight fraction)
469  //
470 
471  // creating & configuring a point geometry driver
474  // casting to the GENIE geometry driver interface
475  geom_driver = dynamic_cast<GeomAnalyzerI *> (pgeom);
476  }
477 
478  // *************************************************************************
479  // * Create / configure the flux driver
480  // *************************************************************************
481  GFluxI * flux_driver =
482  genie::flux::GFluxDriverFactory::Instance().GetFluxDriver("genie::flux::GSimpleNtpFlux");
483  if ( ! flux_driver ) {
484  // couldn't get the requested flux driver ?
485  std::ostringstream s;
486  s << "Known FluxDrivers:\n";
487  const std::vector<std::string>& known =
489  std::vector<std::string>::const_iterator itr = known.begin();
490  for ( ; itr != known.end(); ++itr ) s << " " << (*itr) << "\n";
491  LOG("gevgen_lardm", pFATAL)
492  << "Failed to get any flux driver from GFluxDriverFactory";
493  exit(1);
494  }
495 
496  genie::flux::GFluxFileConfigI* flux_file_config =
497  dynamic_cast<genie::flux::GFluxFileConfigI*>(flux_driver);
498  if ( ! flux_file_config ) {
499  LOG("gevgen_lardm", pFATAL)
500  << "Failed to get GFluxFileConfigI driver from GFluxDriverFactory";
501  exit(1);
502  }
503 
504  //
505  // *** Using the detailed ntuple neutrino flux description
506  //
507  flux_file_config->LoadBeamSimData(gOptFluxFile, "");
508  flux_file_config->SetUpstreamZ(gOptZmin); // was "zmin" from bounding_box
509  flux_file_config->SetNumOfCycles(0);
510  PDGCodeList dm_pdg;
511  dm_pdg.push_back(kPdgDarkMatter);
512  flux_file_config->SetFluxParticles(dm_pdg);
513 
514  // these might come in handy ... avoid repeated dynamic_cast<> calls
515  genie::flux::GFluxFileConfigI* fluxFileConfigI =
516  dynamic_cast<genie::flux::GFluxFileConfigI*>(flux_driver);
517 
518 
519  // *************************************************************************
520  // * Handle chicken/egg problem: geom analyzer vs. flux.
521  // * Need both at this point change geom scan defaults.
522  // *************************************************************************
523  if ( gOptUsingRootGeom ) {
524 
526  dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
527  if ( ! rgeom ) assert(0);
528 
529  rgeom -> SetDebugFlags(gOptDebug);
530 
531  // even if user doesn't specify gOptNScan configure to scan using flux
532  if ( gOptNScan >= 0 ) {
533  LOG("gevgen_lardm", pNOTICE)
534  << "Using ROOTGeomAnalyzer: geom scan using flux: nparticles=" << gOptNScan;
535  rgeom->SetScannerFlux(flux_driver);
536  if ( gOptNScan > 0 ) rgeom->SetScannerNParticles(gOptNScan);
537  } else {
538  int nabs = TMath::Abs(gOptNScan);
539  LOG("gevgen_lardm", pNOTICE)
540  << "Using ROOTGeomAnalyzer: geom scan using box: npoints=nrays=" << nabs;
541  rgeom->SetScannerNPoints(nabs);
542  rgeom->SetScannerNRays(nabs);
543  }
544  }
545 
546  // *************************************************************************
547  // * Create/configure the event generation driver
548  // *************************************************************************
549  GMCJDriver * mcj_driver = new GMCJDriver;
551  mcj_driver->UseFluxDriver(flux_driver);
552  mcj_driver->UseGeomAnalyzer(geom_driver);
553  if ( ( gOptExtMaxPlXml != "" ) && ! gOptWriteMaxPlXml ) {
554  mcj_driver->UseMaxPathLengths(gOptExtMaxPlXml);
555  }
556  mcj_driver->Configure();
557  mcj_driver->UseSplines();
558  mcj_driver->ForceSingleProbScale();
559 
560  if ( ( gOptExtMaxPlXml != "" ) && gOptWriteMaxPlXml ) {
562  dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
563  if ( rgeom ) {
564  const genie::PathLengthList& maxpath = rgeom->GetMaxPathLengths();
565  std::string maxplfile = gOptExtMaxPlXml;
566  maxpath.SaveAsXml(maxplfile);
567  // append extra info to file
568  std::ofstream mpfile(maxplfile.c_str(), std::ios_base::app);
569  mpfile
570  << std::endl
571  << "<!-- this file is only relevant for a setup compatible with:"
572  << std::endl
573  << "geom: " << gOptRootGeom << " top: \"" << gOptRootGeomTopVol << "\""
574  << std::endl
575  << "flux: " << gOptFluxFile
576  << std::endl
577  << "fidcut: " << gOptFidCut
578  << std::endl
579  << "nscan: " << gOptNScan << " (0>= use flux, <0 use box |nscan| points/rays)"
580  << std::endl
581  << "zmin: " << gOptZmin << " (if |zmin| > 1e30, leave ray on flux window)"
582  << std::endl
583  << "-->"
584  << std::endl;
585  mpfile.close();
586  }
587  }
588 
589  // *************************************************************************
590  // * Prepare for writing the output event tree & status file
591  // *************************************************************************
592 
593  // Initialize an Ntuple Writer to save GHEP records into a TTree
596  ntpw.Initialize();
597 
598 
599  std::vector<TBranch*> extraBranches;
600  std::vector<std::string> branchNames;
601  std::vector<std::string> branchClassNames;
602  std::vector<void**> branchObjPointers;
603 
604  // Add custom branch(s) to the standard GENIE event tree so that
605  // info on the flux neutrino parent particle can be passed-through
606  if ( fluxFileConfigI ) {
607  fluxFileConfigI->GetBranchInfo(branchNames,branchClassNames,
608  branchObjPointers);
609  size_t nn = branchNames.size();
610  size_t nc = branchClassNames.size();
611  size_t np = branchObjPointers.size();
612  if ( nn != nc || nc != np ) {
613  LOG("gevgen_lardm", pERROR)
614  << "Inconsistent info back from flux driver "
615  << "for branch info: " << nn << " " << nc << " " << np;
616  } else {
617  for (size_t ii = 0; ii < nn; ++ii) {
618  const char* bname = branchNames[ii].c_str();
619  const char* cname = branchClassNames[ii].c_str();
620  void**& optr = branchObjPointers[ii]; // note critical '&' !
621  if ( ! optr || ! *optr ) continue; // no pointer supplied, skip it
622  int split = 99; // 1
623  LOG("gevgen_lardm", pNOTICE)
624  << "Adding extra branch \"" << bname << "\" of type \""
625  << cname << "\" (" << optr << ") to output tree";
626  TBranch* bptr = ntpw.EventTree()->Branch(bname,cname,optr,32000,split);
627  extraBranches.push_back(bptr);
628 
629  if ( bptr ) {
630  // don't delete this !!! we're sharing
631  bptr->SetAutoDelete(false);
632  } else {
633  LOG("gevgen_lardm", pERROR)
634  << "FAILED to add extra branch \"" << bname << "\" of type \""
635  << cname << "\" to output tree";
636  }
637  } // loop over additions
638  } // same # of entries
639  } // of genie::flux::GFluxFileConfigI type
640 
641  // Create a MC job monitor for a periodically updated status file
642  GMCJMonitor mcjmonitor(gOptRunNu);
643  mcjmonitor.SetRefreshRate(RunOpt::Instance()->MCJobStatusRefreshRate());
644 
645  // *************************************************************************
646  // * Event generation loop
647  // *************************************************************************
648 
649  // define handler to allow signal to end job gracefully
650  signal(SIGTERM,gsSIGTERMhandler);
651 
652  int ievent = 0;
653  while ( ! gSigTERM )
654  {
655  LOG("gevgen_lardm", pINFO)
656  << " *** Generating event............ " << ievent;
657 
658  // In case the required statistics was expressed as 'number of events'
659  // then quit if that number has been generated
660  if ( ievent == gOptNev ) break;
661 
662  // Generate a single event using neutrinos coming from the specified flux
663  // and hitting the specified geometry or target mix
664  EventRecord * event = mcj_driver->GenerateEvent();
665 
666  // Check whether a null event was returned due to the flux driver reaching
667  // the end of the input flux ntuple - exit the event generation loop
668  if ( ! event && flux_driver->End() ) {
669  LOG("gevgen_lardm", pWARN)
670  << "** The flux driver read all the input flux entries: End()==true";
671  break;
672  }
673  if ( ! event ) {
674  LOG("gevgen_lardm", pERROR)
675  << "Got a null generated neutino event! Retrying ...";
676  continue;
677  }
678  LOG("gevgen_lardm", pINFO)
679  << "Generated event: " << *event;
680 
681  // A valid event was generated: flux info (parent decay/prod
682  // position/kinematics) for that simulated event should already
683  // be connected to the right output tree branch
684 
685  // Add event at the output ntuple, refresh the mc job monitor & clean-up
686  ntpw.AddEventRecord(ievent, event);
687  mcjmonitor.Update(ievent,event);
688  delete event;
689  ievent++;
690 
691  } //1
692 
693  // Copy metadata tree, if available
694  if ( fluxFileConfigI ) {
695  TTree* t1 = fluxFileConfigI->GetMetaDataTree();
696  if ( t1 ) {
697  size_t nmeta = t1->GetEntries();
698  TTree* t2 = (TTree*)t1->Clone(0);
699  for (size_t i = 0; i < nmeta; ++i) {
700  t1->GetEntry(i);
701  t2->Fill();
702  }
703  t2->Write();
704  }
705  }
706 
707  LOG("gevgen_lardm", pINFO)
708  << "The GENIE MC job is done generating events - Cleaning up & exiting...";
709 
710  // *************************************************************************
711  // * Save & clean-up
712  // *************************************************************************
713 
714  // Save the generated event tree & close the output file
715  ntpw.Save();
716 
717  // Clean-up
718  delete geom_driver;
719  delete flux_driver;
720  delete mcj_driver;
721  // this list should only be histograms that have (for some reason)
722  // not been handed over to the GCylindTH1Flux driver.
723 
724  LOG("gevgen_lardm", pNOTICE) << "Done!";
725 
726  return 0;
727 }
728 
729 //____________________________________________________________________________
731 {
732  /// potentially load extra libraries that might extend the list of
733  /// potential flux drivers, and how to map short names to classes ...
734  // we really should at this point read some external file to get
735  // an expandable list of libraries ... but for now just hard code it
736 
737  vector<string> extraLibs;
738 
739  ///***** this part should come from reading an external file
740  /// placeholder file $GENIE/config/FluxDriverExpansion.xml
741 
742  extraLibs.push_back("libdk2nuTree");
743  extraLibs.push_back("libdk2nuGenie");
744 
745  ///******* done with fake "read"
746 
747  // see if there are any libraries to load
748  // just let ROOT do it ... check error code on return
749  // but tweak ROOT's ERROR message output so we don't see huge messages
750  // for failures
751  // gErrorIgnoreLevel to kError (from TError.h)
752 
753  Int_t prevErrorLevel = gErrorIgnoreLevel;
754  gErrorIgnoreLevel = kFatal;
755  for (size_t ilib=0; ilib < extraLibs.size(); ++ilib) {
756  // library names should be like libXYZZY without extension [e.g. .so]
757  // but with the leading "lib"
758  int loadStatus = gSystem->Load(extraLibs[ilib].c_str());
759  const char* statWords = "failed to load";
760  if ( loadStatus == 0 ) { statWords = "successfully loaded"; }
761  else if ( loadStatus == 1 ) { statWords = "already had loaded"; }
762  else if ( loadStatus == -1 ) { statWords = "couldn't find library"; }
763  else if ( loadStatus == -2 ) { statWords = "mismatched version"; }
764 
765  LOG("gevgen_lardm",pNOTICE)
766  << statWords << " (" << loadStatus << ") " << extraLibs[ilib];
767  }
768  // restore the ROOT error message level
769  gErrorIgnoreLevel = prevErrorLevel;
770 
771 }
772 
773 //____________________________________________________________________________
774 void GetCommandLineArgs(int argc, char ** argv)
775 {
776  LOG("gevgen_lardm", pINFO) << "Parsing command line arguments";
777 
778  // Common run options. Set defaults and read.
781 
782  // Parse run options for this app
783 
784  CmdLnArgParser parser(argc,argv);
785 
786  // help?
787  bool help = parser.OptionExists('h');
788  if(help) {
789  PrintSyntax();
790  exit(0);
791  }
792 
793  // run number:
794  if ( parser.OptionExists('r') ) {
795  LOG("gevgen_lardm", pDEBUG) << "Reading MC run number";
796  gOptRunNu = parser.ArgAsLong('r');
797  } else {
798  LOG("gevgen_lardm", pDEBUG)
799  << "Unspecified run number - Using default";
800  gOptRunNu = 0;
801  } //-r
802 
803  //
804  // *** geometry
805  //
806 
807  string geom = "";
808  string lunits, dunits;
809  if( parser.OptionExists('g') ) {
810  LOG("gevgen_lardm", pDEBUG) << "Getting input geometry";
811  geom = parser.ArgAsString('g');
812 
813  // is it a ROOT file that contains a ROOT geometry?
814  bool accessible_geom_file =
815  ! (gSystem->AccessPathName(geom.c_str()));
816  if (accessible_geom_file) {
817  gOptRootGeom = geom;
818  gOptUsingRootGeom = true;
819  }
820  } else {
821  LOG("gevgen_lardm", pFATAL)
822  << "No geometry option specified - Exiting";
823  PrintSyntax();
824  exit(1);
825  } //-g
826 
827  // dark matter mass
828  if( parser.OptionExists('M') ) {
829  LOG("gevgen_lardm", pINFO) << "Reading dark matter mass";
830  gOptDMMass = parser.ArgAsDouble('M');
831  } else {
832  LOG("gevgen_lardm", pFATAL) << "Unspecified dark matter mass - Exiting";
833  PrintSyntax();
834  exit(1);
835  } // -M
836 
837  // mediator coupling
838  if( parser.OptionExists('c') ) {
839  LOG("gevgen_lardm", pINFO) << "Reading mediator coupling";
840  gOptZpCoupling = parser.ArgAsDouble('c');
841  } else {
842  LOG("gevgen_lardm", pINFO) << "Unspecified mediator coupling - Using value from config file";
843  gOptZpCoupling = -1.;
844  } // -c
845 
846  // mediator mass ratio
847  if( parser.OptionExists('v') ) {
848  LOG("gevgen_lardm", pINFO) << "Reading mediator mass ratio";
849  gOptMedRatio = parser.ArgAsDouble('v');
850  } else {
851  LOG("gevgen_lardm", pINFO) << "Unspecified mediator mass ratio - Using default";
852  gOptMedRatio = 0.5;
853  } // -v
854 
855  if(gOptUsingRootGeom) {
856  // using a ROOT geometry - get requested geometry units
857 
858  // legth units:
859  if( parser.OptionExists('L') ) {
860  LOG("gevgen_lardm", pDEBUG)
861  << "Checking for input geometry length units";
862  lunits = parser.ArgAsString('L');
863  } else {
864  LOG("gevgen_lardm", pDEBUG) << "Using default geometry length units";
865  lunits = kDefOptGeomLUnits;
866  } // -L
867  // density units:
868  if( parser.OptionExists('D') ) {
869  LOG("gevgen_lardm", pDEBUG)
870  << "Checking for input geometry density units";
871  dunits = parser.ArgAsString('D');
872  } else {
873  LOG("gevgen_lardm", pDEBUG) << "Using default geometry density units";
874  dunits = kDefOptGeomDUnits;
875  } // -D
878 
879  // check whether an event generation volume name has been
880  // specified -- default is the 'top volume'
881  if( parser.OptionExists('t') ) {
882  LOG("gevgen_lardm", pDEBUG) << "Checking for input volume name";
883  gOptRootGeomTopVol = parser.ArgAsString('t');
884  } else {
885  LOG("gevgen_lardm", pDEBUG) << "Using the <master volume>";
886  } // -t
887 
888  // check whether an XML file with the maximum (density weighted)
889  // path lengths for each detector material is specified -
890  // otherwise will compute the max path lengths at job init
891  // if passed name starts with '+', then compute max at job init, but write out the result
892  if ( parser.OptionExists('m') ) {
893  LOG("gevgen_lardm", pDEBUG)
894  << "Checking for maximum path lengths XML file";
895  gOptExtMaxPlXml = parser.ArgAsString('m');
896  gOptWriteMaxPlXml = false;
897  if ( gOptExtMaxPlXml[0] == '+' ) {
898  gOptExtMaxPlXml = gOptExtMaxPlXml.substr(1,std::string::npos);
899  gOptWriteMaxPlXml = true;
900  LOG("gevgen_lardm", pINFO)
901  << "Will write maximum path lengths XML file: " << gOptExtMaxPlXml;
902  }
903  } else {
904  LOG("gevgen_lardm", pDEBUG)
905  << "Will compute the maximum path lengths at job init";
906  gOptExtMaxPlXml = "";
907  } // -m
908 
909  // fidcut:
910  if( parser.OptionExists('F') ) {
911  LOG("gevgen_lardm", pDEBUG) << "Using Fiducial cut?";
912  gOptFidCut = parser.ArgAsString('F');
913  } else {
914  LOG("gevgen_lardm", pDEBUG) << "No fiducial volume cut";
915  gOptFidCut = "";
916  } //-F
917 
918  // how to scan the geometry (if relevant)
919  if( parser.OptionExists('S') ) {
920  LOG("gevgen_lardm", pDEBUG) << "Reading requested geom scan count";
921  gOptNScan = parser.ArgAsInt('S');
922  } else {
923  LOG("gevgen_lardm", pDEBUG) << "No geom scan count was requested";
924  gOptNScan = 0;
925  } //-S
926 
927  // z for flux rays to start
928  if( parser.OptionExists('z') ) {
929  LOG("gevgen_lardm", pDEBUG) << "Reading requested zmin";
930  gOptZmin = parser.ArgAsDouble('z');
931  } else {
932  LOG("gevgen_lardm", pDEBUG) << "No zmin was requested";
933  gOptZmin = -2.0e30; // < -1.0e30 ==> leave it on flux window
934  } //-z
935 
936  // debug flags
937  if ( parser.OptionExists('d') ) {
938  LOG("gevgen_lardm", pDEBUG) << "Reading debug flag value";
939  gOptDebug = parser.ArgAsInt('d');
940  } else {
941  LOG("gevgen_lardm", pDEBUG) << "Unspecified debug flags - Using default";
942  gOptDebug = 0;
943  } //-d
944 
945  } // using root geom?
946 
947  else {
948  // User has specified a target mix.
949  // Decode the list of target pdf codes & their corresponding weight fraction
950  // (specified as 'pdg_code_1[fraction_1],pdg_code_2[fraction_2],...')
951  // See documentation on top section of this file.
952  //
953  gOptTgtMix.clear();
954  vector<string> tgtmix = utils::str::Split(geom,",");
955  if(tgtmix.size()==1) {
956  int pdg = atoi(tgtmix[0].c_str());
957  double wgt = 1.0;
958  gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
959  } else {
960  vector<string>::const_iterator tgtmix_iter = tgtmix.begin();
961  for( ; tgtmix_iter != tgtmix.end(); ++tgtmix_iter) {
962  string tgt_with_wgt = *tgtmix_iter;
963  string::size_type open_bracket = tgt_with_wgt.find("[");
964  string::size_type close_bracket = tgt_with_wgt.find("]");
965  if (open_bracket ==string::npos ||
966  close_bracket==string::npos)
967  {
968  LOG("gevgen_lardm", pFATAL)
969  << "You made an error in specifying the target mix";
970  PrintSyntax();
971  exit(1);
972  }
973  string::size_type ibeg = 0;
974  string::size_type iend = open_bracket;
975  string::size_type jbeg = open_bracket+1;
976  string::size_type jend = close_bracket;
977  int pdg = atoi(tgt_with_wgt.substr(ibeg,iend-ibeg).c_str());
978  double wgt = atof(tgt_with_wgt.substr(jbeg,jend-jbeg).c_str());
979  LOG("gevgen_lardm", pDEBUG)
980  << "Adding to target mix: pdg = " << pdg << ", wgt = " << wgt;
981  gOptTgtMix.insert(map<int, double>::value_type(pdg, wgt));
982 
983  }// tgtmix_iter
984  } // >1 materials in mix
985  } // using tgt mix?
986 
987  //
988  // *** flux
989  //
990  if ( parser.OptionExists('f') ) {
991  LOG("gevgen_lardm", pDEBUG) << "Getting input flux";
992  gOptFluxFile = parser.ArgAsString('f');
993  } else {
994  LOG("gevgen_lardm", pFATAL) << "No flux info was specified - Exiting";
995  PrintSyntax();
996  exit(1);
997  }
998 
999  // number of events to generate
1000  if( parser.OptionExists('n') ) {
1001  LOG("gevgen_lardm", pDEBUG)
1002  << "Reading limit on number of events to generate";
1003  gOptNev = parser.ArgAsInt('n');
1004  } else {
1005  LOG("gevgen_lardm", pDEBUG)
1006  << "Will keep on generating events till the flux driver stops";
1007  gOptNev = -1;
1008  } //-n
1009 
1010  // event file prefix
1011  if( parser.OptionExists('o') ) {
1012  LOG("gevgen_lardm", pDEBUG) << "Reading the event filename prefix";
1013  gOptEvFilePrefix = parser.ArgAsString('o');
1014  } else {
1015  LOG("gevgen_lardm", pDEBUG)
1016  << "Will set the default event filename prefix";
1018  } //-o
1019 
1020 
1021  // random number seed
1022  if( parser.OptionExists("seed") ) {
1023  LOG("gevgen_lardm", pINFO) << "Reading random number seed";
1024  gOptRanSeed = parser.ArgAsLong("seed");
1025  } else {
1026  LOG("gevgen_lardm", pINFO) << "Unspecified random number seed - Using default";
1027  gOptRanSeed = -1;
1028  }
1029 
1030  // input cross-section file
1031  if( parser.OptionExists("cross-sections") ) {
1032  LOG("gevgen_lardm", pINFO) << "Reading cross-section file";
1033  gOptInpXSecFile = parser.ArgAsString("cross-sections");
1034  } else {
1035  LOG("gevgen_lardm", pINFO) << "Unspecified cross-section file";
1036  gOptInpXSecFile = "";
1037  }
1038 
1039  //
1040  // >>> print the command line options
1041  //
1042 
1043  PDGLibrary * pdglib = PDGLibrary::Instance();
1044 
1045  ostringstream gminfo;
1046  if (gOptUsingRootGeom) {
1047  gminfo << "Using ROOT geometry - file = " << gOptRootGeom
1048  << ", top volume = "
1049  << ((gOptRootGeomTopVol.size()==0) ? "<master volume>" : gOptRootGeomTopVol)
1050  << ", max{PL} file = "
1051  << ((gOptExtMaxPlXml.size()==0) ? "<none>" : gOptExtMaxPlXml)
1052  << ", length units = " << lunits
1053  << ", density units = " << dunits;
1054  } else {
1055  gminfo << "Using target mix: ";
1056  map<int,double>::const_iterator iter;
1057  for(iter = gOptTgtMix.begin(); iter != gOptTgtMix.end(); ++iter) {
1058  int pdg_code = iter->first;
1059  double wgt = iter->second;
1060  TParticlePDG * p = pdglib->Find(pdg_code);
1061  if(p) {
1062  string name = p->GetName();
1063  gminfo << "(" << name << ") -> " << 100*wgt << "% / ";
1064  }//p?
1065  }
1066  }
1067 
1068  ostringstream fluxinfo;
1069  fluxinfo << "file = " << gOptFluxFile;
1070 
1071  ostringstream exposure;
1072  if(gOptNev > 0)
1073  exposure << "Number of events = " << gOptNev;
1074 
1075  LOG("gevgen_lardm", pNOTICE)
1076  << "\n\n"
1077  << utils::print::PrintFramedMesg("NuMI expt. event generation job configuration");
1078 
1079  LOG("gevgen_lardm", pNOTICE)
1080  << "\n - Run number: " << gOptRunNu
1081  << "\n - Random number seed: " << gOptRanSeed
1082  << "\n - Using cross-section file: " << gOptInpXSecFile
1083  << "\n - Flux @ " << fluxinfo.str()
1084  << "\n - Geometry @ " << gminfo.str()
1085  << "\n - Exposure @ " << exposure.str();
1086 
1087  LOG("gevgen_lardm", pNOTICE) << *RunOpt::Instance();
1088 }
1089 //____________________________________________________________________________
1090 void PrintSyntax(void)
1091 {
1092  LOG("gevgen_lardm", pFATAL)
1093  << "\n **Syntax**"
1094  << "\n gevgen_lardm [-h] [-r run#]"
1095  << "\n -f flux -g geometry -M dm_mass"
1096  << "\n [-c zp_coupling] [-v med_ratio]"
1097  << "\n [-t top_volume_name_at_geom] [-m max_path_lengths_xml_file]"
1098  << "\n [-L length_units_at_geom] [-D density_units_at_geom]"
1099  << "\n [-n n_of_events] "
1100  << "\n [-o output_event_file_prefix]"
1101  << "\n [-F fid_cut_string] [-S nrays_scan]"
1102  << "\n [-z zmin_start]"
1103  << "\n [--seed random_number_seed]"
1104  << "\n --cross-sections xml_file"
1105  << "\n [--event-generator-list list_name]"
1106  << "\n [--message-thresholds xml_file]"
1107  << "\n [--unphysical-event-mask mask]"
1108  << "\n [--event-record-print-level level]"
1109  << "\n [--mc-job-status-refresh-rate rate]"
1110  << "\n [--cache-file root_file]"
1111  << "\n"
1112  << " Please also read the detailed documentation at "
1113  << "$GENIE/src/Apps/gFNALExptEvGen.cxx"
1114  << "\n";
1115 }
1116 //____________________________________________________________________________
1117 void CreateFidSelection (string fidcut, GeomAnalyzerI* geom_driver)
1118 {
1119  ///
1120  /// User defined fiducial volume cut
1121  /// [0][M]<SHAPE>:val1,val2,...
1122  /// "0" means reverse the cut (i.e. exclude the volume)
1123  /// "M" means the coordinates are given in the ROOT geometry
1124  /// "master" system and need to be transformed to "top vol" system
1125  /// <SHAPE> can be any of "zcyl" "box" "zpoly" "sphere"
1126  /// [each takes different # of args]
1127  /// This must be followed by a ":" and a list of values separated by punctuation
1128  /// (allowed separators: commas , parentheses () braces {} or brackets [] )
1129  /// Value mapping:
1130  /// zcly:x0,y0,radius,zmin,zmax - cylinder along z at (x0,y0) capped at z's
1131  /// box:xmin,ymin,zmin,xmax,ymax,zmax - box w/ upper & lower extremes
1132  /// zpoly:nfaces,x0,y0,r_in,phi,zmin,zmax - nfaces sided polygon in x-y plane
1133  // sphere:x0,y0,z0,radius - sphere of fixed radius at (x0,y0,z0)
1134  /// Examples:
1135  /// 1) 0mbox:0,0,0.25,1,1,8.75
1136  /// exclude (i.e. reverse) a box in master coordinates w/ corners (0,0,0.25) (1,1,8.75)
1137  /// 2) mzpoly:6,(2,-1),1.75,0,{0.25,8.75}
1138  /// six sided polygon in x-y plane, centered at x,y=(2,-1) w/ inscribed radius 1.75
1139  /// no rotation (so first face is in y-z plane +r from center, i.e. hex sits on point)
1140  /// limited to the z range of {0.25,8.75} in the master ROOT geom coordinates
1141  /// 3) zcly:(3,4),5.5,-2,10
1142  /// a cylinder oriented parallel to the z axis in the "top vol" coordinates
1143  /// at x,y=(3,4) with radius 5.5 and z range of {-2,10}
1144  ///
1145  geometry::ROOTGeomAnalyzer * rgeom =
1146  dynamic_cast<geometry::ROOTGeomAnalyzer *>(geom_driver);
1147  if ( ! rgeom ) {
1148  LOG("gevgen_lardm", pWARN)
1149  << "Can not create GeomVolSelectorFiduction,"
1150  << " geometry driver is not ROOTGeomAnalyzer";
1151  return;
1152  }
1153 
1154  LOG("gevgen_lardm", pNOTICE) << "-F " << fidcut;
1155 
1158 
1159  fidsel->SetRemoveEntries(true); // drop segments that won't be considered
1160 
1161  // convert string to lowercase
1162  std::transform(fidcut.begin(),fidcut.end(),fidcut.begin(),::tolower);
1163 
1164  vector<string> strtok = genie::utils::str::Split(fidcut,":");
1165  if ( strtok.size() != 2 ) {
1166  LOG("gevgen_lardm", pWARN)
1167  << "Can not create GeomVolSelectorFiduction,"
1168  << " no \":\" separating type from values. nsplit=" << strtok.size();
1169  for ( unsigned int i=0; i < strtok.size(); ++i )
1170  LOG("gevgen_lardm",pNOTICE)
1171  << "strtok[" << i << "] = \"" << strtok[i] << "\"";
1172  return;
1173  }
1174 
1175  // parse out optional "x" and "m"
1176  string stype = strtok[0];
1177  bool reverse = ( stype.find("0") != string::npos );
1178  bool master = ( stype.find("m") != string::npos ); // action after values are set
1179 
1180  // parse out values
1181  vector<double> vals;
1182  vector<string> valstrs = genie::utils::str::Split(strtok[1]," ,;(){}[]");
1183  vector<string>::const_iterator iter = valstrs.begin();
1184  for ( ; iter != valstrs.end(); ++iter ) {
1185  const string& valstr1 = *iter;
1186  if ( valstr1 != "" ) vals.push_back(atof(valstr1.c_str()));
1187  }
1188  size_t nvals = vals.size();
1189 
1190  std::cout << "ivals = [";
1191  for (unsigned int i=0; i < nvals; ++i) {
1192  if (i>0) cout << ",";
1193  std::cout << vals[i];
1194  }
1195  std::cout << "]" << std::endl;
1196 
1197  // std::vector elements are required to be adjacent so we can treat address as ptr
1198 
1199  if ( stype.find("zcyl") != string::npos ) {
1200  // cylinder along z direction at (x0,y0) radius zmin zmax
1201  if ( nvals < 5 )
1202  LOG("gevgen_lardm", pFATAL) << "MakeZCylinder needs 5 values, not " << nvals
1203  << " fidcut=\"" << fidcut << "\"";
1204  fidsel->MakeZCylinder(vals[0],vals[1],vals[2],vals[3],vals[4]);
1205 
1206  } else if ( stype.find("box") != string::npos ) {
1207  // box (xmin,ymin,zmin) (xmax,ymax,zmax)
1208  if ( nvals < 6 )
1209  LOG("gevgen_lardm", pFATAL) << "MakeBox needs 6 values, not " << nvals
1210  << " fidcut=\"" << fidcut << "\"";
1211  double xyzmin[3] = { vals[0], vals[1], vals[2] };
1212  double xyzmax[3] = { vals[3], vals[4], vals[5] };
1213  fidsel->MakeBox(xyzmin,xyzmax);
1214 
1215  } else if ( stype.find("zpoly") != string::npos ) {
1216  // polygon along z direction nfaces at (x0,y0) radius phi zmin zmax
1217  if ( nvals < 7 )
1218  LOG("gevgen_lardm", pFATAL) << "MakeZPolygon needs 7 values, not " << nvals
1219  << " fidcut=\"" << fidcut << "\"";
1220  int nfaces = (int)vals[0];
1221  if ( nfaces < 3 )
1222  LOG("gevgen_lardm", pFATAL) << "MakeZPolygon needs nfaces>=3, not " << nfaces
1223  << " fidcut=\"" << fidcut << "\"";
1224  fidsel->MakeZPolygon(nfaces,vals[1],vals[2],vals[3],vals[4],vals[5],vals[6]);
1225 
1226  } else if ( stype.find("sphere") != string::npos ) {
1227  // sphere at (x0,y0,z0) radius
1228  if ( nvals < 4 )
1229  LOG("gevgen_lardm", pFATAL) << "MakeZSphere needs 4 values, not " << nvals
1230  << " fidcut=\"" << fidcut << "\"";
1231  fidsel->MakeSphere(vals[0],vals[1],vals[2],vals[3]);
1232 
1233  } else {
1234  LOG("gevgen_lardm", pFATAL)
1235  << "Can not create GeomVolSelectorFiduction for shape \"" << stype << "\"";
1236  }
1237 
1238  if ( master ) {
1239  fidsel->ConvertShapeMaster2Top(rgeom);
1240  LOG("gevgen_lardm", pNOTICE) << "Convert fiducial volume from master to topvol coords";
1241  }
1242  if ( reverse ) {
1243  fidsel->SetReverseFiducial(true);
1244  LOG("gevgen_lardm", pNOTICE) << "Reverse sense of fiducial volume cut";
1245  }
1246  rgeom->AdoptGeomVolSelector(fidsel);
1247 
1248 }
1249 //____________________________________________________________________________
1250 void CreateRockBoxSelection (string fidcut, GeomAnalyzerI* geom_driver)
1251 {
1252 
1253  if( fidcut.find_first_not_of(" \t\n") != 0) // trim any leading whitespace
1254  fidcut.erase( 0, fidcut.find_first_not_of(" \t\n") );
1255 
1256  // convert string to lowercase
1257  std::transform(fidcut.begin(),fidcut.end(),fidcut.begin(),::tolower);
1258 
1260  dynamic_cast<genie::geometry::ROOTGeomAnalyzer *>(geom_driver);
1261  if ( ! rgeom ) {
1262  LOG("gevgen_numi", pWARN)
1263  << "Can not create GeomVolSelectorRockBox,"
1264  << " geometry driver is not ROOTGeomAnalyzer";
1265  return;
1266  }
1267 
1268  LOG("gevgen_numi", pWARN) << "fiducial (rock) cut: " << fidcut;
1269 
1270  // for now, only fiducial no "rock box"
1273 
1274  vector<string> strtok = genie::utils::str::Split(fidcut,":");
1275  if ( strtok.size() != 2 ) {
1276  LOG("gevgen_numi", pWARN)
1277  << "Can not create GeomVolSelectorRockBox,"
1278  << " no \":\" separating type from values. nsplit=" << strtok.size();
1279  for ( unsigned int i=0; i < strtok.size(); ++i )
1280  LOG("gevgen_numi", pWARN)
1281  << "strtok[" << i << "] = \"" << strtok[i] << "\"";
1282  return;
1283  }
1284 
1285  string stype = strtok[0];
1286 
1287  // parse out values
1288  vector<double> vals;
1289  vector<string> valstrs = genie::utils::str::Split(strtok[1]," ,;(){}[]\t\n\r");
1290  vector<string>::const_iterator iter = valstrs.begin();
1291  for ( ; iter != valstrs.end(); ++iter ) {
1292  const string& valstr1 = *iter;
1293  if ( valstr1 != "" ) {
1294  double aval = atof(valstr1.c_str());
1295  LOG("gevgen_numi", pWARN) << "rock value [" << vals.size() << "] "
1296  << aval;
1297  vals.push_back(aval);
1298  }
1299  }
1300  size_t nvals = vals.size();
1301 
1302  rocksel->SetRemoveEntries(true); // drop segments that won't be considered
1303 
1304  // assume coordinates are in the *master* (not "top volume") system
1305  // need to set fTopVolume to fWorldVolume
1306  //fTopVolume = fWorldVolume;
1307  //rgeom->SetTopVolName(fTopVolume.c_str());
1309  rgeom->SetTopVolName(gOptRootGeomMasterVol);
1310 
1311  if ( nvals < 6 ) {
1312  LOG("gevgen_numi", pFATAL) << "rockbox needs at "
1313  << "least 6 values, found "
1314  << nvals << "in \""
1315  << strtok[1] << "\"";
1316  exit(1);
1317 
1318  }
1319  double xyzmin[3] = { vals[0], vals[1], vals[2] };
1320  double xyzmax[3] = { vals[3], vals[4], vals[5] };
1321 
1322  bool rockonly = true;
1323  double wallmin = 800.; // geometry in cm, ( 8 meter buffer)
1324  double dedx = 2.5 * 1.7e-3; // GeV/cm, rho=2.5, 1.7e-3 ~ rock like loss
1325  double fudge = 1.05;
1326 
1327  if ( nvals >= 7 ) rockonly = vals[6];
1328  if ( nvals >= 8 ) wallmin = vals[7];
1329  if ( nvals >= 9 ) dedx = vals[8];
1330  if ( nvals >= 10 ) fudge = vals[9];
1331 
1332  rocksel->SetRockBoxMinimal(xyzmin,xyzmax);
1333  rocksel->SetMinimumWall(wallmin);
1334  rocksel->SetDeDx(dedx/fudge);
1335 
1336  if ( nvals >= 11 ) rocksel->SetExpandFromInclusion((int)vals[10]);
1337 
1338  // if not rock-only then make a tiny exclusion bubble
1339  // call to MakeBox shouldn't be necessary
1340  // should be done by SetRockBoxMinimal but for some GENIE versions isn't
1341  if ( ! rockonly ) rocksel->MakeSphere(0,0,0,1.0e-10);
1342  else rocksel->MakeBox(xyzmin,xyzmax);
1343 
1344  rgeom->AdoptGeomVolSelector(rocksel);
1345 
1346 }
1347 //____________________________________________________________________________
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)
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)
void CreateRockBoxSelection(string fidcut, GeomAnalyzerI *geom_driver)
double ArgAsDouble(char opt)
virtual void LoadBeamSimData(const std::vector< std::string > &filenames, const std::string &det_loc)=0
string ArgAsString(char opt)
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
#define pERROR
Definition: Messenger.h:60
void AddDarkMatter(double mass, double med_ratio)
Definition: PDGLibrary.cxx:113
int gOptDebug
void MakeZCylinder(Double_t x0, Double_t y0, Double_t radius, Double_t zmin, Double_t zmax)
string gOptRootGeomMasterVol
void SetEventGeneratorList(string listname)
Definition: GMCJDriver.cxx:157
void DetermineFluxDriver(string fopt)
void ReadFromCommandLine(int argc, char **argv)
Definition: RunOpt.cxx:107
const char * p
Definition: xmltok.h:285
#define pFATAL
Definition: Messenger.h:57
const int kPdgDarkMatter
Definition: PDGCodes.h:195
virtual const PathLengthList & GetMaxPathLengths(void) const
int gOptNev
string gOptExtMaxPlXml
string kDefOptEvFilePrefix
string gOptInpXSecFile
virtual void ConvertShapeMaster2Top(const ROOTGeomAnalyzer *rgeom)
void Update(int iev, const EventRecord *event)
Definition: GMCJMonitor.cxx:58
OStream cerr
Definition: OStream.cxx:7
double gOptZmin
const std::vector< std::string > & AvailableFluxDrivers() const
int gOptNScan
string gOptFidCut
virtual void SetNumOfCycles(long int ncycle)
limit cycling through input files
virtual void SetScannerNParticles(int np)
Registry * CommonList(const string &file_id, const string &set_name) const
A list of PDG codes.
Definition: PDGCodeList.h:33
bool gOptWriteMaxPlXml
void CreateFidSelection(string fidcut, GeomAnalyzerI *geom_driver)
virtual void SetUpstreamZ(double z0)
void SetRefreshRate(int rate)
Definition: GMCJMonitor.cxx:53
int main(int argc, char **argv)
void UseFluxDriver(GFluxI *flux)
Definition: GMCJDriver.cxx:174
double gOptZpCoupling
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
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
void SaveAsXml(string filename) const
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Index size_type
Definition: typedefs.hpp:11
void ParseFluxFileConfig(string fopt)
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 kDefOptGeomLUnits
map< int, double > gOptTgtMix
void ForceSingleProbScale(void)
Definition: GMCJDriver.cxx:219
virtual void SetScannerNRays(int nr)
bool UseMaxPathLengths(string xml_filename)
Definition: GMCJDriver.cxx:190
double UnitFromString(string u)
Definition: UnitUtils.cxx:26
bool gSigTERM
GENIE Interface for limiting vertex selection in the rock to a volume that depends (in part) on the n...
static void gsSIGTERMhandler(int)
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.
void Save(void)
get the even tree
Definition: NtpWriter.cxx:214
#define pINFO
Definition: Messenger.h:63
void Lock(void)
locks the registry
Definition: Registry.cxx:163
double t2
Long_t gOptRunNu
void AddEventRecord(int ievent, const EventRecord *ev_rec)
save the event tree
Definition: NtpWriter.cxx:69
const ana::Var wgt
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.
void UnLock(void)
unlocks the registry (doesn&#39;t unlock items)
Definition: Registry.cxx:168
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:38
void CustomizeFilenamePrefix(string prefix)
Definition: NtpWriter.cxx:123
string gOptEvFilePrefix
void LoadExtraOptions(void)
void Initialize(void)
add event
Definition: NtpWriter.cxx:95
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
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
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
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
static GFluxDriverFactory & Instance()
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
exit(0)
string gOptFluxFile
enum BeamMode nc
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)
TRandom3 r(0)
long int gOptRanSeed
void PrintSyntax(void)
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 kDefOptGeomDUnits
double gOptGeomDUnits
void MesgThresholds(string inpfile)
Definition: AppInit.cxx:100
double gOptGeomLUnits
void GetCommandLineArgs(int argc, char **argv)
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)
Defines the GENIE Geometry Analyzer Interface.
Definition: GeomAnalyzerI.h:30
string gOptRootGeom
Float_t e
Definition: plot.C:35
void Set(RgIMapPair entry)
Definition: Registry.cxx:282
void RecursiveExhaust(TGeoVolume *topvol, string volnames, bool exhaust)
Definition: GeoUtils.cxx:26
bool gOptUsingRootGeom
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 ParseFluxHst(string fopt)
void CacheFile(string inpfile)
Definition: AppInit.cxx:118
string gOptRootGeomTopVol
void EnableBareXSecPreCalc(bool flag)
Definition: RunOpt.h:61
double gOptDMMass
void push_back(int pdg_code)
Definition: PDGCodeList.cxx:67
virtual TGeoManager * GetGeometry(void) const
static AlgConfigPool * Instance()
genie::GFluxI * GetFluxDriver(const std::string &)
double gOptMedRatio
#define pDEBUG
Definition: Messenger.h:64
GENIE Interface for user-defined flux classes.
Definition: GFluxI.h:37
NtpMCFormat_t kDefOptNtpFormat
enum BeamMode string