LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
GENIEHelper.cxx
Go to the documentation of this file.
1 
11 // C/C++ includes
12 #include <math.h>
13 #include <map>
14 #include <fstream>
15 #include <iomanip>
16 #include <algorithm>
17 #include <sstream>
18 #include <glob.h>
19 #include <cstdlib> // for unsetenv()
20 
21 //ROOT includes
22 #include "TH1.h"
23 #include "TH2.h" //used by GAtmoFlux
24 #include "TF1.h"
25 #include "TFile.h"
26 #include "TDirectory.h"
27 #include "TVector3.h"
28 #include "TLorentzVector.h"
29 #include "TCollection.h"
30 #include "TSystem.h"
31 #include "TString.h"
32 #include "TRandom.h" //needed for gRandom to be defined
33 #include "TRegexp.h"
34 #include "TMath.h"
35 #include "TStopwatch.h"
36 #include "TRotation.h"
37 
38 //GENIE includes
39 #ifdef GENIE_PRE_R3
40  #include "GENIE/Messenger/Messenger.h"
41  #include "GENIE/Conventions/GVersion.h"
42  #include "GENIE/Conventions/Units.h"
43  #include "GENIE/EVGCore/EventRecord.h"
44  #include "GENIE/EVGDrivers/GMCJDriver.h"
45  #include "GENIE/GHEP/GHepUtils.h"
46  #include "GENIE/FluxDrivers/GCylindTH1Flux.h"
47  #include "GENIE/FluxDrivers/GMonoEnergeticFlux.h"
48  #include "GENIE/FluxDrivers/GNuMIFlux.h"
49  #include "GENIE/FluxDrivers/GSimpleNtpFlux.h"
50  #include "GENIE/FluxDrivers/GFluxDriverFactory.h"
51  #include "GENIE/FluxDrivers/GBGLRSAtmoFlux.h" //for atmo nu generation
52  #include "GENIE/FluxDrivers/GFLUKAAtmoFlux.h" //for atmo nu generation
53  #if __GENIE_RELEASE_CODE__ >= GRELCODE(2,12,2)
54  #include "GENIE/FluxDrivers/GHAKKMAtmoFlux.h" // for atmo nu generation
55  #endif
56  #include "GENIE/FluxDrivers/GAtmoFlux.h" //for atmo nu generation
57 
58  #pragma GCC diagnostic push
59  #pragma GCC diagnostic ignored "-Wunused-variable"
60  #include "GENIE/Conventions/Constants.h" //for calculating event kinematics
61  #pragma GCC diagnostic pop
62 
63  #include "GENIE/PDG/PDGLibrary.h"
64  #include "GENIE/PDG/PDGCodes.h"
65  #include "GENIE/Utils/AppInit.h"
66  #include "GENIE/Utils/RunOpt.h"
67 
68  #include "GENIE/Geo/ROOTGeomAnalyzer.h"
69  #include "GENIE/Geo/GeomVolSelectorFiducial.h"
70  #include "GENIE/Geo/GeomVolSelectorRockBox.h"
71  #include "GENIE/Utils/StringUtils.h"
72  #include "GENIE/Utils/XmlParserUtils.h"
73  #include "GENIE/Interaction/InitialState.h"
74  #include "GENIE/Interaction/Interaction.h"
75  #include "GENIE/Interaction/Kinematics.h"
76  #include "GENIE/Interaction/KPhaseSpace.h"
77  #include "GENIE/Interaction/ProcessInfo.h"
78  #include "GENIE/Interaction/XclsTag.h"
79  #include "GENIE/GHEP/GHepParticle.h"
80  #include "GENIE/PDG/PDGCodeList.h"
81 
82  #include "GENIE/FluxDrivers/GFluxBlender.h"
83  #include "GENIE/FluxDrivers/GFlavorMixerI.h"
84  #include "GENIE/FluxDrivers/GFlavorMap.h"
85  #include "GENIE/FluxDrivers/GFlavorMixerFactory.h"
86 
87 #else
88  // GENIE R-3 reorganized headers
89  #include "GENIE/Framework/Messenger/Messenger.h"
90  #include "GENIE/Framework/Conventions/GVersion.h"
91  #include "GENIE/Framework/Utils/StringUtils.h"
92  #include "GENIE/Framework/Utils/XmlParserUtils.h"
93  #include "GENIE/Framework/Messenger/Messenger.h"
94 
95  #pragma GCC diagnostic push
96  #pragma GCC diagnostic ignored "-Wunused-variable"
97  // constants for calculating event kinematics
98  #include "GENIE/Framework/Conventions/Constants.h"
99  #pragma GCC diagnostic pop
100 
101  #include "GENIE/Framework/Conventions/Units.h"
102  #include "GENIE/Framework/ParticleData/PDGCodes.h"
103  #include "GENIE/Framework/ParticleData/PDGCodeList.h"
104  #include "GENIE/Framework/ParticleData/PDGLibrary.h"
105 
106  #include "GENIE/Framework/GHEP/GHepParticle.h"
107  #include "GENIE/Framework/GHEP/GHepUtils.h"
108 
109  #include "GENIE/Framework/Interaction/InitialState.h"
110  #include "GENIE/Framework/Interaction/Interaction.h"
111  #include "GENIE/Framework/Interaction/Kinematics.h"
112  #include "GENIE/Framework/Interaction/KPhaseSpace.h"
113  #include "GENIE/Framework/Interaction/ProcessInfo.h"
114  #include "GENIE/Framework/Interaction/XclsTag.h"
115 
116  #include "GENIE/Framework/EventGen/GFluxI.h"
117  #include "GENIE/Framework/EventGen/EventRecord.h"
118  #include "GENIE/Framework/EventGen/GMCJDriver.h"
119 
120  #include "GENIE/Framework/Utils/AppInit.h"
121  #include "GENIE/Framework/Utils/RunOpt.h"
122 
123  #include "GENIE/Tools/Geometry/ROOTGeomAnalyzer.h"
124  #include "GENIE/Tools/Geometry/GeomVolSelectorFiducial.h"
125  #include "GENIE/Tools/Geometry/GeomVolSelectorRockBox.h"
126 
127  #include "GENIE/Tools/Flux/GFluxBlender.h"
128  #include "GENIE/Tools/Flux/GFlavorMixerI.h"
129  #include "GENIE/Tools/Flux/GFlavorMap.h"
130  #include "GENIE/Tools/Flux/GFlavorMixerFactory.h"
131  #include "GENIE/Tools/Flux/GFluxDriverFactory.h"
132 
133  #include "GENIE/Tools/Flux/GCylindTH1Flux.h"
134  #include "GENIE/Tools/Flux/GMonoEnergeticFlux.h"
135  #include "GENIE/Tools/Flux/GNuMIFlux.h"
136  #include "GENIE/Tools/Flux/GSimpleNtpFlux.h"
137  #include "GENIE/Tools/Flux/GAtmoFlux.h"
138  #include "GENIE/Tools/Flux/GBGLRSAtmoFlux.h" // was GBartolAtmoFlux.h
139  #include "GENIE/Tools/Flux/GFLUKAAtmoFlux.h" // was GFlukaAtmo3DFlux.h
140  #include "GENIE/Tools/Flux/GHAKKMAtmoFlux.h" //
141 
142 #endif
143 
144 
145 // dk2nu flux
146 #include "dk2nu/tree/dk2nu.h"
147 #include "dk2nu/tree/dkmeta.h"
148 #include "dk2nu/tree/NuChoice.h"
149 #include "dk2nu/genie/GDk2NuFlux.h"
150 
151 // NuGen includes
152 #include "nugen/EventGeneratorBase/evgenbase.h"
153 #include "nugen/EventGeneratorBase/GENIE/GENIEHelper.h"
154 
155 #include "nugen/EventGeneratorBase/GENIE/GENIE2ART.h"
156 #include "nugen/EventGeneratorBase/GENIE/EvtTimeShiftFactory.h"
157 #include "nugen/EventGeneratorBase/GENIE/EvtTimeShiftI.h"
158 
159 #include "nugen/EventGeneratorBase/GENIE/GPowerSpectrumAtmoFlux.h"
160 
161 // nusimdata includes
162 #include "nusimdata/SimulationBase/MCTruth.h"
163 #include "nusimdata/SimulationBase/MCFlux.h"
164 #include "nusimdata/SimulationBase/GTruth.h"
165 #include "nusimdata/SimulationBase/MCParticle.h"
166 #include "nusimdata/SimulationBase/MCNeutrino.h"
167 
168 // Framework includes
169 #include "art/Framework/Services/Registry/ServiceHandle.h"
170 #include "cetlib/search_path.h"
171 #include "cetlib/getenv.h"
172 #include "cetlib/split_path.h"
173 #include "cetlib_except/exception.h"
174 #include "fhiclcpp/ParameterSet.h"
175 #include "messagefacility/MessageLogger/MessageLogger.h"
176 
177 #include "cetlib_except/exception.h"
178 
179 // can't find IFDH_service.h header ... unless ups depends on ifdh_art
180 //
181 // #undef USE_IFDH_SERVICE
182 
183 #ifndef NO_IFDH_LIB
184  #define USE_IFDH_SERVICE 1
185  // IFDHC
186  #ifdef USE_IFDH_SERVICE
187  #include "ifdh_art/IFDHService/IFDH_service.h"
188  #else
189  // bare IFDHC
190  #include "ifdh.h"
191  #endif
192 #else
193  #undef USE_IFDH_SERVICE
194  // nothing doing ... use ifdef to hide any reference that might need header
195  #include <cassert>
196 #endif
197 
198 // The GENIE logging macros require the name 'Messenger' to be
199 // accessible during the macro expansion. This is a bug in GENIE, but
200 // for now, the easiest way to do this is to bring genie::Messenger
201 // into the global namespace by a using declaration.
202 using genie::Messenger;
203 
204 namespace evgb {
205 
206  static const int kNue = 0;
207  static const int kNueBar = 1;
208  static const int kNuMu = 2;
209  static const int kNuMuBar = 3;
210  static const int kNuTau = 4;
211  static const int kNuTauBar = 5;
212 
213  //--------------------------------------------------
214  GENIEHelper::GENIEHelper(fhicl::ParameterSet const& pset,
215  TGeoManager* geoManager,
216  std::string const& rootFile,
217  double const& detectorMass)
218  : fGeoManager (geoManager)
219  , fGeoFile (rootFile)
220  , fGenieEventRecord (0)
221  , fGeomD (0)
222  , fFluxD (0)
223  , fFluxD2GMCJD (0)
224  , fDriver (0)
225  , fIFDH (0)
226  , fHelperRandom (0)
227  , fUseHelperRndGen4GENIE(pset.get< bool >("UseHelperRndGen4GENIE",true))
228  , fTimeShifter (0)
229  , fFluxType (pset.get< std::string >("FluxType") )
230  , fFluxSearchPaths (pset.get< std::string >("FluxSearchPaths","") )
231  , fFluxFilePatterns (pset.get< std::vector<std::string> >("FluxFiles") )
232  , fMaxFluxFileMB (pset.get< int >("MaxFluxFileMB", 2000) ) // 2GB max default
233  , fMaxFluxFileNumber (pset.get< int >("MaxFluxFileNumber",9999) ) // at most 9999 files
234  , fFluxCopyMethod (pset.get< std::string >("FluxCopyMethod","DIRECT")) // "DIRECT" = old direct access method
235  , fFluxCleanup (pset.get< std::string >("FluxCleanup","/var/tmp") ) // "ALWAYS", "NEVER", "/var/tmp"
236  , fBeamName (pset.get< std::string >("BeamName") )
237  , fFluxRotCfg (pset.get< std::string >("FluxRotCfg","none") )
238  , fFluxRotValues (pset.get< std::vector<double> >("FluxRotValues", {} ) ) // default empty vector
239  , fFluxRotation (0)
240  , fTopVolume (pset.get< std::string >("TopVolume") )
241  , fWorldVolume ("volWorld")
242  , fDetLocation (pset.get< std::string >("DetectorLocation") )
243  , fFluxUpstreamZ (pset.get< double >("FluxUpstreamZ", -2.e30) )
244  , fEventsPerSpill (pset.get< double >("EventsPerSpill", 0) )
245  , fPOTPerSpill (pset.get< double >("POTPerSpill", 0.0) )
246  , fHistEventsPerSpill(0.)
247  , fSpillEvents (0)
248  , fSpillExposure (0.)
249  , fTotalExposure (0.)
250  , fMonoEnergy (pset.get< double >("MonoEnergy", 2.0) )
251  , fFunctionalFlux (pset.get< std::string >("FunctionalFlux", "x") )
252  , fFunctionalBinning (pset.get< int >("FunctionalBinning", 10000) )
253  , fEmin (pset.get< double >("FluxEmin", 0) )
254  , fEmax (pset.get< double >("FluxEmax", 10) )
255  , fBeamRadius (pset.get< double >("BeamRadius", 3.0) )
256  , fDetectorMass (detectorMass)
257  , fSurroundingMass (pset.get< double >("SurroundingMass", 0.) )
258  , fGlobalTimeOffset (pset.get< double >("GlobalTimeOffset", 1.e4) )
259  , fRandomTimeOffset (pset.get< double >("RandomTimeOffset", 1.e4) )
260  , fSpillTimeConfig (pset.get< std::string >("SpillTimeConfig", "") )
261  , fAddGenieVtxTime (pset.get< bool >("AddGenieVtxTime", false) )
262  , fForceApplyFlxWgt (pset.get< bool >("ForceApplyFlxWgt", true) )
263  , fGenFlavors (pset.get< std::vector<int> >("GenFlavors") )
264  , fAtmoEmin (pset.get< double >("AtmoEmin", 0.1) )
265  , fAtmoEmax (pset.get< double >("AtmoEmax", 10.0) )
266  , fAtmoRl (pset.get< double >("Rl", 20.0) )
267  , fAtmoRt (pset.get< double >("Rt", 20.0) )
268  , fAtmoSpectralIndex (pset.get< double >("SpectralIndex", 2.0) )
269  , fEnvironment (pset.get< std::vector<std::string> >("Environment") )
270  , fXSecTable (pset.get< std::string >("XSecTable", "") ) //e.g. "gxspl-FNALsmall.xml"
271  , fTuneName (pset.get< std::string >("TuneName","${GENIE_XSEC_TUNE}") )
272  , fEventGeneratorList(pset.get< std::string >("EventGeneratorList", "Default") )
273  , fGXMLPATH (pset.get< std::string >("GXMLPATH", "") )
274  , fGMSGLAYOUT (pset.get< std::string >("GMSGLAYOUT", "") ) // [BASIC] or SIMPLE
275  , fGENIEMsgThresholds(pset.get< std::string >("GENIEMsgThresholds", "") ) // : separate list of files
276  , fGHepPrintLevel (pset.get< int >("GHepPrintLevel", -1) ) // see GHepRecord::SetPrintLevel() -1=no-print
277  , fMixerConfig (pset.get< std::string >("MixerConfig", "none") )
278  , fMixerBaseline (pset.get< double >("MixerBaseline", 0.) )
279  , fUseBlenderDist (pset.get< bool >("UseBlenderDist", true) )
280  , fFiducialCut (pset.get< std::string >("FiducialCut", "none") )
281  , fGeomScan (pset.get< std::string >("GeomScan", "default") )
282  , fDebugFlags (pset.get< unsigned int >("DebugFlags", 0) )
283  {
284 
285  // fEnvironment is (generally) deprecated ... print out any settings
286  // but for some make it fatal because it's no longer supported and
287  // can lead to unexpected behaviour
288  if ( fEnvironment.size() > 0 ) {
289  bool fatal = false;
290  std::ostringstream fenvout, fenvfatal;
291  fenvout << " fEnviroment.size() = " << fEnvironment.size();
292  for (size_t i = 0; i < fEnvironment.size(); i += 2) {
293  fenvout << std::endl << " [" << std::setw(20) << std::left
294  << fEnvironment[i] << "] ==> "
295  << fEnvironment[i+1];
296  // fatal: GEVGL, GSPLOAD
297  // harmless: GXMLPATH, GMSGLAYOUT, GMSGCONF, GPRODMODE
298  if ( fEnvironment[i] == "GEVGL" ||
299  fEnvironment[i] == "GSPLOAD" ) {
300  fatal = true;
301  std::string altparam = "";
302  if ( fEnvironment[i] == "GEVGL" ) altparam = "EventGeneratorList";
303  if ( fEnvironment[i] == "GSPLOAD" ) altparam = "XSecTable";
304 
305  fenvfatal << std::endl << "use of \"" << fEnvironment[i] << "\""
306  << " is no longer supported as a Environment fcl parameter pair key,"
307  << "\n"
308  << " please remove it from from the list and use"
309  << " the appropriate direct fcl parameter \""
310  << altparam << "\"" << std::endl;
311  }
312  }
313  mf::LogInfo("GENIEHelper")
314  << " Use of fEnvironment parameters is deprecated:\n"
315  << fenvout.str();
316  if ( fatal ) {
317  mf::LogError("GENIEHelper") << fenvfatal.str();
318  throw cet::exception("GENIEHelper")
319  << " bad use of Environment fcl parameter";
320  }
321  }
322 
323  // for histogram, mono-energetic, functional form fluxes ...
324  std::vector<double> beamCenter (pset.get< std::vector<double> >("BeamCenter") );
325  std::vector<double> beamDirection(pset.get< std::vector<double> >("BeamDirection"));
326  fBeamCenter.SetXYZ(beamCenter[0], beamCenter[1], beamCenter[2]);
327  fBeamDirection.SetXYZ(beamDirection[0], beamDirection[1], beamDirection[2]);
328 
329  // special processing of GSEED (GENIE's random seed)... priority:
330  // if set in .fcl file RandomSeed variable, use that
331  // else if already set in environment use that
332  // else use evgb::GetRandomNumberSeed()
333  int dfltseed;
334  const char* gseedstr = std::getenv("GSEED");
335  if ( gseedstr ) {
336  dfltseed = strtol(gseedstr,NULL,0);
337  } else {
338  dfltseed = evgb::GetRandomNumberSeed();
339  }
340  int seedval = pset.get< int >("RandomSeed", dfltseed);
341  // initialize random # generator for use within GENIEHelper
342  mf::LogInfo("GENIEHelper") << "Init HelperRandom with seed " << seedval;
343  fHelperRandom = new TRandom3(seedval);
344 
345  // clean up user input
346  // also classifies flux type to simplify tests
347  // e.g. atmo_ tree_
348  RegularizeFluxType();
349 
352 
353  // for tree-based fluxes
354  // (e.g. "tree_numi" (nee "ntuple"), "tree_simple" and "tree_dk2nu")
355  // we don't care about order and don't want duplicates
356  // for others order might matter
357  if ( fFluxType.find("tree_") == 0 ) SqueezeFilePatterns();
358 
359  ExpandFluxPaths();
360  if ( fFluxCopyMethod == "DIRECT" ) ExpandFluxFilePatternsDirect();
361  else ExpandFluxFilePatternsIFDH();
362 
365  if ( fFluxType.find("atmo_") == 0 ||
366  fFluxType.find("astro_") == 0 ) {
367  BuildFluxRotation();
368  }
369 
372  // they should come in pairs of variable name key, then value
373 
374  // Process GXMLPATH extensions first, so they are available
375  // when GENIE starts to get initialized; these might be
376  // alternative locations for configurations (including
377  // the GENIE Messenger system).
378  SetGXMLPATH();
379 
380  // Also set GENIE log4cpp Messenger layout format before
381  // initializing GENIE (can't be changed after singleton is created)
382  SetGMSGLAYOUT();
383 
384  // Now initialize GENIE Messenger service
385  StartGENIEMessenger(pset.get<std::string>("ProductionMode","false"));
386 
387  // In case we're printing the event record, how verbose should it be
388  genie::GHepRecord::SetPrintLevel(fGHepPrintLevel);
389 
390  // Set GENIE's random # seed
391  mf::LogInfo("GENIEHelper")
392  << "Init genie::utils::app_init::RandGen() with seed " << seedval;
393  genie::utils::app_init::RandGen(seedval);
394 
395  // special things for atmos fluxes
396  if ( fFluxType.find("atmo_") == 0 ) AtmoFluxCheck();
397 
398  // make the histogram associations
399  if ( fFluxType.find("histogram") == 0 ) HistogramFluxCheck();
400 
401  std::string flvlist;
402  for ( std::vector<int>::iterator itr = fGenFlavors.begin(); itr != fGenFlavors.end(); itr++ )
403  flvlist += Form(" %d",*itr);
404 
405  if ( fFluxType.find("mono") == 0 ) {
406  fEventsPerSpill = 1;
407  mf::LogInfo("GENIEHelper")
408  << "Generating monoenergetic (" << fMonoEnergy
409  << " GeV) neutrinos with the following flavors: "
410  << flvlist;
411  } else if ( fFluxType.find("function") == 0 ) {
412  fEventsPerSpill = 1;
413  mf::LogInfo("GENIEHelper")
414  << "Generating neutrinos using the functional form: "
415  << fFunctionalFlux << " E [" << fEmin << ":" << fEmax
416  << "] GeV with " << fFunctionalBinning << " bins "
417  << "with the following flavors: " << flvlist;
418  } else if (fFluxType.find("PowerSpectrum") != std::string::npos) {
419  mf::LogInfo("GENIEHelper")
420  << "Generating neutrinos using the a power spectrum with Spectral index = "
421  << fAtmoSpectralIndex
422  << ", with the following flavors: " << flvlist
423  << "\nThe energy range is between: " << fAtmoEmin << " GeV and "
424  << fAtmoEmax << " GeV."
425  << '\n'
426  << " Generation surface of: (" << fAtmoRl << ","
427  << fAtmoRt << ")";
428  } else {
429 
430  // flux methods other than "mono" and "function" require files
431  std::string fileliststr;
432  if ( fSelectedFluxFiles.empty() ) {
433  fileliststr = "NO FLUX FILES FOUND!";
434  mf::LogWarning("GENIEHelper") << fileliststr;
435  }
436  else {
437  std::vector<std::string>::iterator sitr = fSelectedFluxFiles.begin();
438  for ( ; sitr != fSelectedFluxFiles.end(); sitr++) {
439  fileliststr += "\n\t";
440  fileliststr += *sitr;
441  }
442  }
443  mf::LogInfo("GENIEHelper")
444  << "Generating flux with the following flavors: " << flvlist
445  << "\nand these file patterns: " << fileliststr;
446 
447  }
448 
449  // disallow conflicting settings here
450  if ( ( fEventsPerSpill != 0 && fPOTPerSpill != 0 ) ||
451  ( fEventsPerSpill == 0 && fPOTPerSpill == 0 ) ) {
452  throw cet::exception("GENIEHelper")
453  << "one or the other of EventsPerSpill ("
454  << fEventsPerSpill << ") or "
455  << "POTPerSpill ("
456  << fPOTPerSpill << ") needs to be zero (but not both)";
457  }
458 
459  if ( fEventsPerSpill != 0 ) {
460  mf::LogInfo("GENIEHelper") << "Generating " << fEventsPerSpill
461  << " events for each spill";
462 
463  } else {
464  mf::LogInfo("GENIEHelper") << "Using " << fPOTPerSpill
465  << " pot for each spill";
466  }
467 
468  // how to distribute events in time
469  if ( fSpillTimeConfig != "" ) {
470  evgb::EvtTimeShiftFactory& timeShiftFactory = evgb::EvtTimeShiftFactory::Instance();
471  fTimeShifter = timeShiftFactory.GetEvtTimeShift(fSpillTimeConfig);
472  if ( fTimeShifter ) {
473  if ( ! fTimeShifter->IsRandomGeneratorSeeded() ) {
474  mf::LogInfo("GENIEHelper")
475  << "Initialize spill time random generator with seed " << seedval;
476  fTimeShifter->GetRandomGenerator()->SetSeed(seedval);
477  }
478  fTimeShifter->PrintConfig();
479  } else {
480  timeShiftFactory.Print();
481  }
482  }
483 
484  return;
485  }
486 
487  //--------------------------------------------------
488  GENIEHelper::~GENIEHelper()
489  {
490  // user request writing out the scan of the geometry
491  if ( fGeomD && fMaxPathOutInfo != "" ) {
492  genie::geometry::ROOTGeomAnalyzer* rgeom =
493  dynamic_cast<genie::geometry::ROOTGeomAnalyzer*>(fGeomD);
494 
495  string filename = "maxpathlength.xml";
496  mf::LogInfo("GENIEHelper")
497  << "Saving MaxPathLengths as: \"" << filename << "\"";
498 
499  const genie::PathLengthList& maxpath = rgeom->GetMaxPathLengths();
500 
501  maxpath.SaveAsXml(filename);
502  // append extra info to file
503  std::ofstream mpfile(filename.c_str(), std::ios_base::app);
504  mpfile
505  << std::endl
506  << "<!-- this file is only relevant for a setup compatible with:"
507  << std::endl
508  << fMaxPathOutInfo
509  << std::endl
510  << "-->"
511  << std::endl;
512  mpfile.close();
513  } // finished writing max path length XML file (if requested)
514 
515  // protect against lack of driver due to not getting to Initialize()
516  // (called from module's beginRun() method)
517  if ( ! fDriver || ! fFluxD ) {
518  mf::LogInfo("GENIEHelper")
519  << "~GENIEHelper called, but previously failed to construct "
520  << ( (fDriver) ? "":" genie::GMCJDriver" )
521  << ( (fFluxD) ? "":" genie::GFluxI" );
522  } else {
523 
524  double probscale = fDriver->GlobProbScale();
525  double rawpots = 0;
526 
527  // rather than ask individual flux drivers for info
528  // use the unified interface
529 
530  genie::flux::GFluxExposureI* fexposure =
531  dynamic_cast<genie::flux::GFluxExposureI*>(fFluxD);
532  if ( fexposure ) {
533  rawpots = fexposure->GetTotalExposure();
534  }
535  genie::flux::GFluxFileConfigI* ffileconfig =
536  dynamic_cast<genie::flux::GFluxFileConfigI*>(fFluxD);
537  if ( ffileconfig ) {
538  ffileconfig->PrintConfig();
539  }
540 
541  mf::LogInfo("GENIEHelper")
542  << " Total Exposure " << fTotalExposure
543  << " GMCJDriver GlobProbScale " << probscale
544  << " FluxDriver base pots " << rawpots
545  << " corrected POTS " << rawpots/TMath::Max(probscale,1.0e-100);
546  }
547 
548  // clean up owned genie object (other genie obj are ref ptrs)
549  delete fGenieEventRecord;
550  delete fDriver;
551  delete fHelperRandom;
552 
553 #ifndef NO_IFDH_LIB
554  #ifdef USE_IFDH_SERVICE
555  art::ServiceHandle<IFDH> ifdhp;
556  if ( fFluxCleanup.find("ALWAYS") == 0 ) {
557  ifdhp->cleanup();
558  } else if ( fFluxCleanup.find("/var/tmp") == 0 ) {
559  auto ffitr = fSelectedFluxFiles.begin();
560  for ( ; ffitr != fSelectedFluxFiles.end(); ++ffitr ) {
561  std::string ff = *ffitr;
562  if ( ff.find("/var/tmp") == 0 ) {
563  mf::LogDebug("GENIEHelper") << "delete " << ff;
564  ifdhp->rm(ff);
565  }
566  }
567  }
568  #else
569  if ( fIFDH ) {
570  if ( fFluxCleanup.find("ALWAYS") == 0 ) {
571  fIFDH->cleanup();
572  } else if ( fFluxCleanup.find("/var/tmp") == 0 ) {
573  auto ffitr = fSelectedFluxFiles.begin();
574  for ( ; ffitr != fSelectedFluxFiles.end(); ++ffitr ) {
575  std::string ff = *ffitr;
576  if ( ff.find("/var/tmp") == 0 ) {
577  mf::LogDebug("GENIEHelper") << "delete " << ff;
578  fIFDH->rm(ff);
579  }
580  }
581  }
582  delete fIFDH;
583  fIFDH = 0;
584  }
585  #endif
586 #endif
587 
588  }
589 
590  //--------------------------------------------------
591  double GENIEHelper::TotalHistFlux()
592  {
593  if ( fFluxType.find("mono") == 0 ||
594  fFluxType.find("function") == 0 ||
595  fFluxType.find("tree_") == 0 ||
596  fFluxType.find("atmo_") == 0 ) {
597  // shouldn't be asking for this for any of the above
598  // perhaps not-not "function"?
599  return -999.;
600  }
601 
602  return fTotalHistFlux;
603  }
604 
605  //--------------------------------------------------
606  void GENIEHelper::Initialize()
607  {
608 #ifdef GENIE_PRE_R3
609  // trigger early initialization of PDG database & GENIE message service
610  // just to get it out of the way and not intermixed with other output
611  genie::PDGLibrary::Instance();
612 #else
613  // get the GENIE banner out of the way
614  // no longer can use genie::PDGLibrary::Instance() to do this
615  // because that must happen, in some cases in v3_02_xx, after the tune
616  // is determined
617  // banner is triggered by first use of GENIE Messenger
618  // avoid using GENIE macros (possible conflict with mf macros)
619  // LOG("GENIE",pInfo) << "Trigger GENIE banner";
620  (*genie::Messenger::Instance())("GENIE") << log4cpp::Priority::INFO
621  << "Trigger GENIE banner";
622 #endif
623 
624 
625  mf::LogInfo("GENIEHelper") << "GENIE EventGeneratorList using \""
626  << fEventGeneratorList << "\"";
627 
628  evgb::SetEventGeneratorListAndTune(fEventGeneratorList,fTuneName);
629 
630 #ifdef GENIE_PRE_R3
631  // no tunes ... so can report current setting if different from fcl
632 #else
633  genie::RunOpt* grunopt = genie::RunOpt::Instance();
634 
635  // RunOpt's SetEventGeneratorList() wasn't introduced until R-3
636  // so above call could not have modified it in RunOpt
637  std::string currentEvtGenListName = grunopt->EventGeneratorList();
638  if ( currentEvtGenListName != fEventGeneratorList ) {
639  mf::LogInfo("GENIEHelper") << " EventGeneratorList name changed from \""
640  << fEventGeneratorList << "\" to \""
641  << currentEvtGenListName << "\""
642  << " by evgb::SetEventGeneratorListAndTune";
643  fEventGeneratorList = currentEvtGenListName;
644  }
645 
646  std::string currentTuneName = grunopt->Tune()->Name();
647  if ( currentTuneName != fTuneName ) {
648  mf::LogInfo("GENIEHelper") << " TuneName name changed from \""
649  << fTuneName << "\" to \""
650  << currentTuneName << "\""
651  << " by evgb::SetEventGeneratorListAndTune";
652  fTuneName = currentTuneName;
653  }
654 #endif
655 
656 
657  fDriver = new genie::GMCJDriver(); // this needs to be before ConfigGeomScan
658  // let the driver know which list
659  // (for R-3 this is in _addition_ to setting it in RunOpt)
660  // pre-R-3 one could not set it in RunOpt directly (ie. without parsing
661  // command line args), thus don't try to fetch it from there.
662  // post-R-3 one *-could-* use
663  // fDriver->SetEventGeneratorList(RunOpt::Instance()->EventGeneratorList());
664  fDriver->SetEventGeneratorList(fEventGeneratorList);
665 
666  // Figure out which cross section file to use
667  // post R-2_8_0 this actually triggers reading the file
668  ReadXSecTable();
669 
670  // initialize the Geometry and Flux drivers
671  InitializeGeometry();
672  InitializeFluxDriver();
673 
674  fDriver->UseFluxDriver(fFluxD2GMCJD);
675  fDriver->UseGeomAnalyzer(fGeomD);
676 
677  // must come after creation of Geom, Flux and GMCJDriver
678  ConfigGeomScan(); // could trigger fDriver->UseMaxPathLengths(*xmlfile*)
679 
680  fDriver->Configure(); // could trigger GeomDriver::ComputeMaxPathLengths()
681  fDriver->UseSplines();
682  fDriver->ForceSingleProbScale();
683 
684  if ( fFluxType.find("histogram") == 0 && fEventsPerSpill < 0.01 ) {
685  // fluxes are assumed to be given in units of neutrinos/cm^2/1e20POT/energy
686  // integral over all fluxes removes energy dependence
687  // histograms should have bin width that reflects the value of the /energy bit
688  // ie if /energy = /50MeV then the bin width should be 50 MeV
689 
690  // determine product of pot/spill, mass, and cross section
691  // events = flux * pot * 10^-38 cm^2 (xsec) * (mass detector (in kg) / nucleon mass (in kg))
692  fXSecMassPOT = 1.e-38*1.e-20;
693  fXSecMassPOT *= fPOTPerSpill*(fDetectorMass+fSurroundingMass)/(1.67262158e-27);
694 
695  mf::LogInfo("GENIEHelper") << "Number of events per spill will be based on poisson mean of "
696  << fXSecMassPOT*fTotalHistFlux;
697 
698  fHistEventsPerSpill = fHelperRandom->Poisson(fXSecMassPOT*fTotalHistFlux);
699  }
700 
701  // set the pot/event counters to zero
702  fSpillEvents = 0;
703  fSpillExposure = 0.;
704  fTotalExposure = 0.;
705 
706  // If the flux driver knows how to keep track of exposure (time,pots)
707  // reset it now as some might have been used in determining
708  // the geometry maxpathlength or internally scanning for weights.
709  // Should be happen for all cases since R-2_8_0 .. but no harm doing it ourselves
710 
711  fFluxD->Clear("CycleHistory");
712 
713  return;
714  }
715 
716  //--------------------------------------------------
717  void GENIEHelper::RegularizeFluxType()
718  {
719  // Regularize fFluxType string to sensible setting
720 
721  std::string tmpFluxType = fFluxType;
722 
723  // remove lead/trailing whitespace
724  size_t ftlead = tmpFluxType.find_first_not_of(" \t\n");
725  if ( ftlead ) tmpFluxType.erase( 0, ftlead );
726  size_t ftlen = tmpFluxType.length();
727  size_t fttrail = tmpFluxType.find_last_not_of(" \t\n");
728  if ( fttrail != ftlen ) tmpFluxType.erase( fttrail+1, ftlen );
729 
730  // strip off leading catagories ... we'll put them back later
731  // so we don't accidently allow arbitrary strings
732  if ( tmpFluxType.find("atmos_") == 0 ) tmpFluxType.erase(0,6);
733  if ( tmpFluxType.find("atmo_") == 0 ) tmpFluxType.erase(0,5);
734  if ( tmpFluxType.find("tree_") == 0 ) tmpFluxType.erase(0,5);
735 
736  // make reasonable inferences of what the user intended
737 
738  // simple fluxes
739  if ( tmpFluxType.find("hist") != std::string::npos ) tmpFluxType = "histogram";
740  if ( tmpFluxType.find("func") != std::string::npos ) tmpFluxType = "function";
741  if ( tmpFluxType.find("mono") != std::string::npos ) tmpFluxType = "mono";
742  // Atmospheric fluxes
743  // prior to R-2_11_0 BGLRS was "BARTOL" and HAKKM was "HONDA"
744  if ( tmpFluxType.find("FLUKA") != std::string::npos ) tmpFluxType = "atmo_FLUKA";
745  if ( tmpFluxType.find("BARTOL") != std::string::npos ) tmpFluxType = "atmo_BGLRS";
746  if ( tmpFluxType.find("BGLRS") != std::string::npos ) tmpFluxType = "atmo_BGLRS";
747  if ( tmpFluxType.find("HONDA") != std::string::npos ) tmpFluxType = "atmo_HAKKM";
748  if ( tmpFluxType.find("HAKKM") != std::string::npos ) tmpFluxType = "atmo_HAKKM";
749  if ( tmpFluxType.find("POWER") != std::string::npos ) tmpFluxType = "atmo_PowerSpectrum";
750  // TTree-based fluxes (old "ntuple" is really "numi")
751  // we're allowed to randomize the order here, and squeeze out duplicates
752  if ( tmpFluxType.find("simple") != std::string::npos ) tmpFluxType = "tree_simple";
753  if ( tmpFluxType.find("ntuple") != std::string::npos ) tmpFluxType = "tree_numi";
754  if ( tmpFluxType.find("numi") != std::string::npos ) tmpFluxType = "tree_numi";
755  if ( tmpFluxType.find("dk2nu") != std::string::npos ) tmpFluxType = "tree_dk2nu";
756 
757  fFluxType = tmpFluxType;
758  }
759 
760  //--------------------------------------------------
761  void GENIEHelper::SqueezeFilePatterns()
762  {
763  // for "numi" (nee "ntuple"), "simple" and "dk2nu" squeeze the patterns
764  // so there are no duplicates; for the others we want to preserve order
765 
766  // convert vector<> to a set<> and back to vector<>
767  // to avoid having duplicate patterns in the list
768  std::set<std::string> fluxpattset(fFluxFilePatterns.begin(),
769  fFluxFilePatterns.end());
771  //std::copy(fFluxFilePatterns.begin(),fFluxFilePatterns.end(),
772  // std::inserter(fluxpattset,fluxpattset.begin()));
773  fFluxFilePatterns.clear(); // clear vector, copy unique set back
774  std::copy(fluxpattset.begin(),fluxpattset.end(),
775  std::back_inserter(fFluxFilePatterns));
776  }
777 
778  //--------------------------------------------------
779  void GENIEHelper::AtmoFluxCheck()
780  {
782 
783  if ( fGenFlavors.size() != fSelectedFluxFiles.size() ) {
784  mf::LogInfo("GENIEHelper")
785  << "ERROR: The number of generated neutrino flavors ("
786  << fGenFlavors.size()
787  << ") doesn't correspond to the number of files ("
788  << fSelectedFluxFiles.size() << ")!!!";
789  throw cet::exception("GENIEHelper")
790  << "ERROR: atmo_ flavors != files";
791  } else {
792  std::ostringstream atmofluxmatch;
793  for (size_t indx=0; indx < fGenFlavors.size(); ++indx ) {
794  atmofluxmatch << " " << std::setw(3) << fGenFlavors[indx]
795  << " " << fSelectedFluxFiles[indx] << "\n";
796  }
797  mf::LogInfo("GENIEHelper")
798  << "atmo flux assignment : \n" << atmofluxmatch.str();
799  }
800 
801  if ( fEventsPerSpill != 1 ) {
802  mf::LogInfo("GENIEHelper")
803  << "ERROR: For Atmospheric Neutrino generation,"
804  << " EventPerSpill need to be 1!!";
805  throw cet::exception("GENIEHelper")
806  << "ERROR: " << fFluxType << " EventsPerSpill wasn't 1 ("
807  << fEventsPerSpill << ")";
808  }
809 
810  std::ostringstream atmofluxinfo;
811 
812  if (fFluxType.find("FLUKA") != std::string::npos ){
813  atmofluxinfo << " The fluxes are from FLUKA";
814  }
815  else if (fFluxType.find("BARTOL") != std::string::npos ||
816  fFluxType.find("BGLRS") != std::string::npos ){
817  atmofluxinfo << " The fluxes are from BARTOL/BGLRS";
818  }
819  else if (fFluxType.find("HONDA") != std::string::npos ||
820  fFluxType.find("HAKKM") != std::string::npos ){
821  atmofluxinfo << " The fluxes are from HONDA/HAKKM";
822  }
823  else if (fFluxType.find("PowerSpectrum") != std::string::npos){
824  atmofluxinfo << " The fluxes are interpolated HONDA";
825  }
826  else {
827  mf::LogInfo("GENIEHelper")
828  << "Unknown atmo_ flux simulation: " << fFluxType;
829  throw cet::exception("GENIEHelper")
830  << "ERROR: bad atmo_ flux type " << fFluxType;
831  }
832 
833  atmofluxinfo
834  << '\n'
835  << " The energy range is between: " << fAtmoEmin << " GeV and "
836  << fAtmoEmax << " GeV."
837  << '\n'
838  << " Generation surface of: (" << fAtmoRl << ","
839  << fAtmoRt << ")";
840 
841  mf::LogInfo("GENIEHelper") << atmofluxinfo.str();
842 
843  }
844 
845  //--------------------------------------------------
846  void GENIEHelper::HistogramFluxCheck()
847  {
848 
849  mf::LogInfo("GENIEHelper")
850  << "setting beam direction and center at "
851  << fBeamDirection.X() << " " << fBeamDirection.Y()
852  << " " << fBeamDirection.Z()
853  << " (" << fBeamCenter.X() << "," << fBeamCenter.Y()
854  << "," << fBeamCenter.Z()
855  << ") with radius " << fBeamRadius;
856 
857  TDirectory *savedir = gDirectory;
858 
859  fFluxHistograms.clear();
860 
861  TFile tf((*fSelectedFluxFiles.begin()).c_str());
862  tf.ls();
863 
864  for ( std::vector<int>::iterator flvitr = fGenFlavors.begin(); flvitr != fGenFlavors.end(); flvitr++){
865  const char* histname = "none";
866  switch ( *flvitr ) {
867  case 12: histname = "nue"; break;
868  case -12: histname = "nuebar"; break;
869  case 14: histname = "numu"; break;
870  case -14: histname = "numubar"; break;
871  case 16: histname = "nutau"; break;
872  case -16: histname = "nutaubar"; break;
873  default: {
874  throw cet::exception("GENIEHelper")
875  << "ERROR: no support for histogram flux with flavor PDG="
876  << *flvitr;
877  }
878  }
879  fFluxHistograms.push_back(dynamic_cast<TH1D *>(tf.Get(histname)));
880  }
881 
882  for ( unsigned int i = 0; i < fFluxHistograms.size(); ++i ) {
883  fFluxHistograms[i]->SetDirectory(savedir);
884  fTotalHistFlux += fFluxHistograms[i]->Integral();
885  }
886 
887  mf::LogInfo("GENIEHelper")
888  << "total histogram flux over desired flavors = "
889  << fTotalHistFlux;
890 
891  }
892 
893  //--------------------------------------------------
894  void GENIEHelper::InitializeGeometry()
895  {
896  genie::geometry::ROOTGeomAnalyzer *rgeom =
897  new genie::geometry::ROOTGeomAnalyzer(fGeoManager);
898 
899  // pass some of the debug flag bits on to the geometry manager
900  int geomFlags = ( fDebugFlags >> 16 ) & 0xFF ;
901  if ( geomFlags ) {
902  int keep = ( geomFlags >> 7 );
903  mf::LogInfo("GENIEHelper")
904  << "InitializeGeometry set debug 0x"
905  << std::hex << geomFlags << std::dec
906  << " keepSegPath " << keep;
907  rgeom->SetDebugFlags(geomFlags);
908  if ( keep ) rgeom->SetKeepSegPath(true);
909  }
910 
911  // get the world volume name from the geometry
912  fWorldVolume = fGeoManager->GetTopVolume()->GetName();
913 
914  // the detector geometry uses cgs units.
915  rgeom->SetLengthUnits(genie::units::centimeter);
916  rgeom->SetDensityUnits(genie::units::gram_centimeter3);
917  rgeom->SetTopVolName(fTopVolume.c_str());
918  rgeom->SetMixtureWeightsSum(1.);
919  mf::LogInfo("GENIEHelper")
920  << "GENIEHelper::InitializeGeometry using TopVolume '"
921  << fTopVolume << "'";
922  // casting to the GENIE geometry driver interface
923  fGeomD = rgeom; // dynamic_cast<genie::GeomAnalyzerI *>(rgeom);
924  InitializeFiducialSelection();
925 
926  return;
927  }
928 
929  //--------------------------------------------------
930  void GENIEHelper::InitializeFiducialSelection()
931  {
932  genie::GeomAnalyzerI* geom_driver = fGeomD; // GENIEHelper name -> gNuMIExptEvGen name
933  std::string fidcut = fFiducialCut; // ditto
934 
935  if( fidcut.find_first_not_of(" \t\n") != 0) // trim any leading whitespace
936  fidcut.erase( 0, fidcut.find_first_not_of(" \t\n") );
937 
938  // convert string to lowercase
939  std::transform(fidcut.begin(),fidcut.end(),fidcut.begin(),::tolower);
940 
941  if ( "" == fidcut || "none" == fidcut ) return;
942 
943  if ( fidcut.find("rock") != string::npos ) {
944  // deal with RockBox separately than basic shapes
945  InitializeRockBoxSelection();
946  return;
947  }
948 
949  // below is as it is in $GENIE/src/support/numi/EvGen/gNuMIExptEvGen
950  // except the change in message logger from log4cpp (GENIE) to cet's MessageLogger used by art
951 
966  // sphere:x0,y0,z0,radius - sphere of fixed radius at (x0,y0,z0)
978  genie::geometry::ROOTGeomAnalyzer * rgeom =
979  dynamic_cast<genie::geometry::ROOTGeomAnalyzer *>(geom_driver);
980  if ( ! rgeom ) {
981  mf::LogWarning("GENIEHelpler")
982  << "Can not create GeomVolSelectorFiduction,"
983  << " geometry driver is not ROOTGeomAnalyzer";
984  return;
985  }
986 
987  mf::LogInfo("GENIEHelper") << "fiducial cut: " << fidcut;
988 
989  // for now, only fiducial no "rock box"
990  genie::geometry::GeomVolSelectorFiducial* fidsel =
991  new genie::geometry::GeomVolSelectorFiducial();
992 
993  fidsel->SetRemoveEntries(true); // drop segments that won't be considered
994 
995  vector<string> strtok = genie::utils::str::Split(fidcut,":");
996  if ( strtok.size() != 2 ) {
997  mf::LogWarning("GENIEHelper")
998  << "Can not create GeomVolSelectorFiduction,"
999  << " no \":\" separating type from values. nsplit=" << strtok.size();
1000  for ( unsigned int i=0; i < strtok.size(); ++i )
1001  mf::LogWarning("GENIEHelper")
1002  << "strtok[" << i << "] = \"" << strtok[i] << "\"";
1003  return;
1004  }
1005 
1006  // parse out optional "x" and "m"
1007  string stype = strtok[0];
1008  bool reverse = ( stype.find("0") != string::npos );
1009  bool master = ( stype.find("m") != string::npos ); // action after values are set
1010 
1011  // parse out values
1012  vector<double> vals;
1013  vector<string> valstrs = genie::utils::str::Split(strtok[1]," ,;(){}[]");
1014  vector<string>::const_iterator iter = valstrs.begin();
1015  for ( ; iter != valstrs.end(); ++iter ) {
1016  const string& valstr1 = *iter;
1017  if ( valstr1 != "" ) vals.push_back(atof(valstr1.c_str()));
1018  }
1019  size_t nvals = vals.size();
1020  // pad it out to at least 7 entries to avoid index issues if used
1021  for ( size_t nadd = 0; nadd < 7-nvals; ++nadd ) vals.push_back(0);
1022 
1023  //std::cout << "ivals = [";
1024  //for (unsigned int i=0; i < nvals; ++i) {
1025  // if (i>0) cout << ",";
1026  // std::cout << vals[i];
1027  //}
1028  //std::cout << "]" << std::endl;
1029 
1030  // std::vector elements are required to be adjacent so we can treat address as ptr
1031 
1032  if ( stype.find("zcyl") != string::npos ) {
1033  // cylinder along z direction at (x0,y0) radius zmin zmax
1034  if ( nvals < 5 )
1035  mf::LogError("GENIEHelper") << "MakeZCylinder needs 5 values, not " << nvals
1036  << " fidcut=\"" << fidcut << "\"";
1037  fidsel->MakeZCylinder(vals[0],vals[1],vals[2],vals[3],vals[4]);
1038 
1039  } else if ( stype.find("box") != string::npos ) {
1040  // box (xmin,ymin,zmin) (xmax,ymax,zmax)
1041  if ( nvals < 6 )
1042  mf::LogError("GENIEHelper") << "MakeBox needs 6 values, not " << nvals
1043  << " fidcut=\"" << fidcut << "\"";
1044  double xyzmin[3] = { vals[0], vals[1], vals[2] };
1045  double xyzmax[3] = { vals[3], vals[4], vals[5] };
1046  fidsel->MakeBox(xyzmin,xyzmax);
1047 
1048  } else if ( stype.find("zpoly") != string::npos ) {
1049  // polygon along z direction nfaces at (x0,y0) radius phi zmin zmax
1050  if ( nvals < 7 )
1051  mf::LogError("GENIEHelper") << "MakeZPolygon needs 7 values, not " << nvals
1052  << " fidcut=\"" << fidcut << "\"";
1053  int nfaces = (int)vals[0];
1054  if ( nfaces < 3 )
1055  mf::LogError("GENIEHelper") << "MakeZPolygon needs nfaces>=3, not " << nfaces
1056  << " fidcut=\"" << fidcut << "\"";
1057  fidsel->MakeZPolygon(nfaces,vals[1],vals[2],vals[3],vals[4],vals[5],vals[6]);
1058 
1059  } else if ( stype.find("sphere") != string::npos ) {
1060  // sphere at (x0,y0,z0) radius
1061  if ( nvals < 4 )
1062  mf::LogError("GENIEHelper") << "MakeZSphere needs 4 values, not " << nvals
1063  << " fidcut=\"" << fidcut << "\"";
1064  fidsel->MakeSphere(vals[0],vals[1],vals[2],vals[3]);
1065 
1066  } else {
1067  mf::LogError("GENIEHelper")
1068  << "Can not create GeomVolSelectorFiduction for shape \"" << stype << "\"";
1069  }
1070 
1071  if ( master ) {
1072  fidsel->ConvertShapeMaster2Top(rgeom);
1073  mf::LogInfo("GENIEHelper") << "Convert fiducial volume from master to topvol coords";
1074  }
1075  if ( reverse ) {
1076  fidsel->SetReverseFiducial(true);
1077  mf::LogInfo("GENIEHelper") << "Reverse sense of fiducial volume cut";
1078  }
1079 
1080  rgeom->AdoptGeomVolSelector(fidsel);
1081 
1082  }
1083 
1084  //--------------------------------------------------
1085  void GENIEHelper::InitializeRockBoxSelection()
1086  {
1087  genie::GeomAnalyzerI* geom_driver = fGeomD; // GENIEHelper name -> gNuMIExptEvGen name
1088  std::string fidcut = fFiducialCut; // ditto
1089 
1090  if( fidcut.find_first_not_of(" \t\n") != 0) // trim any leading whitespace
1091  fidcut.erase( 0, fidcut.find_first_not_of(" \t\n") );
1092 
1093  // convert string to lowercase
1094  std::transform(fidcut.begin(),fidcut.end(),fidcut.begin(),::tolower);
1095 
1096  genie::geometry::ROOTGeomAnalyzer * rgeom =
1097  dynamic_cast<genie::geometry::ROOTGeomAnalyzer *>(geom_driver);
1098  if ( ! rgeom ) {
1099  mf::LogWarning("GENIEHelpler")
1100  << "Can not create GeomVolSelectorRockBox,"
1101  << " geometry driver is not ROOTGeomAnalyzer";
1102  return;
1103  }
1104 
1105  mf::LogInfo("GENIEHelper") << "fiducial (rock) cut: " << fidcut;
1106 
1107  // for now, only fiducial no "rock box"
1108  genie::geometry::GeomVolSelectorRockBox* rocksel =
1109  new genie::geometry::GeomVolSelectorRockBox();
1110 
1111  vector<string> strtok = genie::utils::str::Split(fidcut,":");
1112  if ( strtok.size() != 2 ) {
1113  mf::LogWarning("GENIEHelper")
1114  << "Can not create GeomVolSelectorRockBox,"
1115  << " no \":\" separating type from values. nsplit=" << strtok.size();
1116  for ( unsigned int i=0; i < strtok.size(); ++i )
1117  mf::LogWarning("GENIEHelper")
1118  << "strtok[" << i << "] = \"" << strtok[i] << "\"";
1119  return;
1120  }
1121 
1122  string stype = strtok[0];
1123 
1124  // parse out values
1125  vector<double> vals;
1126  vector<string> valstrs = genie::utils::str::Split(strtok[1]," ,;(){}[]\t\n\r");
1127  vector<string>::const_iterator iter = valstrs.begin();
1128  for ( ; iter != valstrs.end(); ++iter ) {
1129  const string& valstr1 = *iter;
1130  if ( valstr1 != "" ) {
1131  double aval = atof(valstr1.c_str());
1132  mf::LogDebug("GENIEHelper") << "rock value [" << vals.size() << "] "
1133  << aval;
1134  vals.push_back(aval);
1135  }
1136  }
1137  size_t nvals = vals.size();
1138 
1139  rocksel->SetRemoveEntries(true); // drop segments that won't be considered
1140 
1141  // assume coordinates are in the *master* (not "top volume") system
1142  // need to set fTopVolume to fWorldVolume as Sample() will keep setting it
1143  fTopVolume = fWorldVolume;
1144  rgeom->SetTopVolName(fTopVolume.c_str());
1145 
1146  if ( nvals < 6 ) {
1147  throw cet::exception("GENIEHelper") << "rockbox needs at "
1148  << "least 6 values, found "
1149  << nvals << "in \""
1150  << strtok[1] << "\"";
1151 
1152  }
1153  double xyzmin[3] = { vals[0], vals[1], vals[2] };
1154  double xyzmax[3] = { vals[3], vals[4], vals[5] };
1155 
1156  bool rockonly = true;
1157  double wallmin = 800.; // geometry in cm, ( 8 meter buffer)
1158  double dedx = 2.5 * 1.7e-3; // GeV/cm, rho=2.5, 1.7e-3 ~ rock like loss
1159  double fudge = 1.05;
1160 
1161  if ( nvals >= 7 ) rockonly = vals[6];
1162  if ( nvals >= 8 ) wallmin = vals[7];
1163  if ( nvals >= 9 ) dedx = vals[8];
1164  if ( nvals >= 10 ) fudge = vals[9];
1165 
1166  rocksel->SetRockBoxMinimal(xyzmin,xyzmax);
1167  rocksel->SetMinimumWall(wallmin);
1168  rocksel->SetDeDx(dedx/fudge);
1169 
1170  // if not rock-only then make a tiny exclusion bubble
1171  // call to MakeBox shouldn't be necessary
1172  // should be done by SetRockBoxMinimal but for some GENIE versions isn't
1173  if ( ! rockonly ) rocksel->MakeSphere(0,0,0,1.0e-10);
1174  else rocksel->MakeBox(xyzmin,xyzmax);
1175 
1176  rgeom->AdoptGeomVolSelector(rocksel);
1177  }
1178 
1179  //--------------------------------------------------
1180  void GENIEHelper::InitializeFluxDriver()
1181  {
1182 
1183  // simplify a lot of things ...
1184  // but for now this part only handles the 3 ntuple styles
1185  // that support the GFluxFileConfig mix-in
1186  // not the atmos, histo or mono versions
1187 
1188  std::string fluxName = "";
1189 
1190  // what looks like the start of a fully qualified class name
1191  // or one of our tree_ classes "numi" "simple" "dk2nu"
1192  // but only know how to configure those that dervie from:
1193  // genie::flux::GFluxFileConfigI*
1194  if ( fFluxType.find("genie::flux::") != std::string::npos )
1195  fluxName = fFluxType;
1196  else if ( fFluxType.find("tree_numi") == 0 )
1197  fluxName = "genie::flux::GNuMIFlux";
1198  else if ( fFluxType.find("tree_simple") == 0 )
1199  fluxName = "genie::flux::GSimpleNtpFlux";
1200  else if ( fFluxType.find("tree_dk2nu") == 0 )
1201  fluxName = "genie::flux::GDk2NuFlux";
1202 
1203  if ( fluxName != "" ) {
1204  // any fall through to hopefully be handled below ...
1205  genie::flux::GFluxDriverFactory& fluxDFactory =
1206  genie::flux::GFluxDriverFactory::Instance();
1207  fFluxD = fluxDFactory.GetFluxDriver(fluxName);
1208  if ( ! fFluxD ) {
1209  mf::LogInfo("GENIEHelper")
1210  << "genie::flux::GFluxDriverFactory had not result for '"
1211  << fluxName << "' (fFluxType was '" << fFluxType << "'";
1212  // fall through in hopes someone later picks it up
1213  } else {
1214  // got something
1215  // for the ones that support GFluxFileConfigI (numi,simple,dk2nu)
1216  // initialize them
1217  genie::flux::GFluxFileConfigI* ffileconfig =
1218  dynamic_cast<genie::flux::GFluxFileConfigI*>(fFluxD);
1219  if ( ffileconfig ) {
1220  ffileconfig->LoadBeamSimData(fSelectedFluxFiles,fDetLocation);
1221  ffileconfig->PrintConfig();
1222  // initialize to only use neutrino flavors requested by user
1223  genie::PDGCodeList probes;
1224  for ( std::vector<int>::iterator flvitr = fGenFlavors.begin(); flvitr != fGenFlavors.end(); flvitr++ )
1225  probes.push_back(*flvitr);
1226  ffileconfig->SetFluxParticles(probes);
1227  if ( TMath::Abs(fFluxUpstreamZ) < 1.0e30 ) ffileconfig->SetUpstreamZ(fFluxUpstreamZ);
1228  }
1229  }
1230  } // is genie::flux:: or tree_{numi|simple|dk2nu}
1231 
1232  if ( fFluxType.find("histogram") == 0 ) {
1233 
1234  genie::flux::GCylindTH1Flux* histFlux = new genie::flux::GCylindTH1Flux();
1235 
1236  // now add the different fluxes - fluxes were added to the vector in the same
1237  // order that the flavors appear in fGenFlavors
1238  int ctr = 0;
1239  for ( std::vector<int>::iterator i = fGenFlavors.begin(); i != fGenFlavors.end(); i++ ) {
1240  histFlux->AddEnergySpectrum(*i, fFluxHistograms[ctr]);
1241  ++ctr;
1242  } //end loop to add flux histograms to driver
1243 
1244  histFlux->SetNuDirection(fBeamDirection);
1245  histFlux->SetBeamSpot(fBeamCenter);
1246  histFlux->SetTransverseRadius(fBeamRadius);
1247 
1248  fFluxD = histFlux; // dynamic_cast<genie::GFluxI *>(histFlux);
1249 
1250  } // end if using a histogram
1251  else if ( fFluxType.find("mono") == 0 ) {
1252 
1253  // weight each species equally in the generation
1254  double weight = 1./(1.*fGenFlavors.size());
1255  //make a map of pdg to weight codes
1256  std::map<int, double> pdgwmap;
1257  for ( std::vector<int>::iterator i = fGenFlavors.begin(); i != fGenFlavors.end(); i++ )
1258  pdgwmap[*i] = weight;
1259 
1260  genie::flux::GMonoEnergeticFlux *monoflux = new genie::flux::GMonoEnergeticFlux(fMonoEnergy, pdgwmap);
1261  monoflux->SetDirectionCos(fBeamDirection.X(), fBeamDirection.Y(), fBeamDirection.Z());
1262  monoflux->SetRayOrigin(fBeamCenter.X(), fBeamCenter.Y(), fBeamCenter.Z());
1263  fFluxD = monoflux; // dynamic_cast<genie::GFluxI *>(monoflux);
1264 
1265  } //end if using monoenergetic beam
1266  else if ( fFluxType.find("function") == 0 ) {
1267 
1268  genie::flux::GCylindTH1Flux* histFlux = new genie::flux::GCylindTH1Flux();
1269  TF1* input_func = new TF1("input_func", fFunctionalFlux.c_str(), fEmin, fEmax);
1270  TH1D* spectrum = new TH1D("spectrum", "neutrino flux", fFunctionalBinning, fEmin, fEmax);
1271  spectrum->Add(input_func);
1272 
1273  for ( std::vector<int>::iterator i = fGenFlavors.begin(); i != fGenFlavors.end(); i++ ) {
1274  histFlux->AddEnergySpectrum(*i, spectrum);
1275  } //end loop to add flux histograms to driver
1276 
1277  histFlux->SetNuDirection(fBeamDirection);
1278  histFlux->SetBeamSpot(fBeamCenter);
1279  histFlux->SetTransverseRadius(fBeamRadius);
1280 
1281  fFluxD = histFlux; // dynamic_cast<genie::GFluxI *>(histFlux);
1282  delete input_func;
1283  } //end if using function beam
1284 
1285  else if(fFluxType.find("PowerSpectrum") != std::string::npos){
1286  genie::flux::GPowerSpectrumAtmoFlux *power_flux = new genie::flux::GPowerSpectrumAtmoFlux;
1287  power_flux->SetSpectralIndex(fAtmoSpectralIndex);
1288  power_flux->SetFlavors(fGenFlavors);
1289  power_flux->SetMinEnergy(fAtmoEmin);
1290  power_flux->SetMaxEnergy(fAtmoEmax);
1291 
1292  mf::LogInfo("GENIEHelper") << "Setting Emin=" << fEmin << " ; Emax=" << fEmax;
1293 
1294  for ( size_t j = 0; j < fGenFlavors.size(); ++j ) {
1295  int flavor = fGenFlavors[j];
1296  std::string flxfile = fSelectedFluxFiles[j];
1297  power_flux->AddFluxFile(flavor,flxfile);
1298  }
1299 
1300 
1301  if ( fFluxRotation ){
1302  power_flux->SetUserCoordSystem(*fFluxRotation);
1303  std::ostringstream atmoCfgText;
1304  const int w=13, p=6;
1305  auto old_p = atmoCfgText.precision(p);
1306  atmoCfgText << "\n UserCoordSystem rotation:\n"
1307  << " [ "
1308  << std::setw(w) << fFluxRotation->XX() << " "
1309  << std::setw(w) << fFluxRotation->XY() << " "
1310  << std::setw(w) << fFluxRotation->XZ() << " ]\n"
1311  << " [ "
1312  << std::setw(w) << fFluxRotation->YX() << " "
1313  << std::setw(w) << fFluxRotation->YY() << " "
1314  << std::setw(w) << fFluxRotation->YZ() << " ]\n"
1315  << " [ "
1316  << std::setw(w) << fFluxRotation->ZX() << " "
1317  << std::setw(w) << fFluxRotation->ZY() << " "
1318  << std::setw(w) << fFluxRotation->ZZ() << " ]\n";
1319  atmoCfgText.precision(old_p);
1320  mf::LogInfo("GENIEHelper") << atmoCfgText.str();
1321  }
1322 
1323  power_flux->LoadFluxData();
1324 
1325  // configure flux generation surface:
1326  power_flux->SetRadii(fAtmoRl, fAtmoRt);
1327 
1328  fFluxD = power_flux;
1329  }
1330 
1331  // Using the atmospheric fluxes
1332  else if ( fFluxType.find("atmo_") == 0 ) {
1333 
1334  // Instantiate appropriate concrete flux driver
1335  genie::flux::GAtmoFlux *atmo_flux_driver = 0;
1336 
1337  if ( fFluxType.find("FLUKA") != std::string::npos ) {
1338  genie::flux::GFLUKAAtmoFlux * fluka_flux =
1339  new genie::flux::GFLUKAAtmoFlux;
1340  atmo_flux_driver = dynamic_cast<genie::flux::GAtmoFlux *>(fluka_flux);
1341  }
1342  if ( fFluxType.find("BARTOL") != std::string::npos ||
1343  fFluxType.find("BGLRS") != std::string::npos ) {
1344  genie::flux::GBGLRSAtmoFlux * bartol_flux =
1345  new genie::flux::GBGLRSAtmoFlux;
1346  atmo_flux_driver = dynamic_cast<genie::flux::GAtmoFlux *>(bartol_flux);
1347  }
1348 #if __GENIE_RELEASE_CODE__ >= GRELCODE(2,12,2)
1349  if (fFluxType.find("atmo_HONDA") != std::string::npos ||
1350  fFluxType.find("atmo_HAKKM") != std::string::npos ) {
1351  genie::flux::GHAKKMAtmoFlux * honda_flux =
1352  new genie::flux::GHAKKMAtmoFlux;
1353  atmo_flux_driver = dynamic_cast<genie::flux::GAtmoFlux *>(honda_flux);
1354  }
1355 #endif
1356 
1357  atmo_flux_driver->ForceMinEnergy(fAtmoEmin);
1358  atmo_flux_driver->ForceMaxEnergy(fAtmoEmax);
1359  if ( fFluxRotation ) atmo_flux_driver->SetUserCoordSystem(*fFluxRotation);
1360 
1361  std::ostringstream atmoCfgText;
1362  atmoCfgText << "Configuration for " << fFluxType
1363  << ", Rl " << fAtmoRl << " Rt " << fAtmoRt;
1364  for ( size_t j = 0; j < fGenFlavors.size(); ++j ) {
1365  int flavor = fGenFlavors[j];
1366  std::string flxfile = fSelectedFluxFiles[j];
1367  atmo_flux_driver->AddFluxFile(flavor,flxfile); // pre-R-2_11_0 was SetFluxFile()
1368  atmoCfgText << "\n FLAVOR: " << std::setw(3) << flavor
1369  << " FLUX FILE: " << flxfile;
1370  }
1371  if ( fFluxRotation ) {
1372  const int w=13, p=6;
1373  auto old_p = atmoCfgText.precision(p);
1374  atmoCfgText << "\n UserCoordSystem rotation:\n"
1375  << " [ "
1376  << std::setw(w) << fFluxRotation->XX() << " "
1377  << std::setw(w) << fFluxRotation->XY() << " "
1378  << std::setw(w) << fFluxRotation->XZ() << " ]\n"
1379  << " [ "
1380  << std::setw(w) << fFluxRotation->YX() << " "
1381  << std::setw(w) << fFluxRotation->YY() << " "
1382  << std::setw(w) << fFluxRotation->YZ() << " ]\n"
1383  << " [ "
1384  << std::setw(w) << fFluxRotation->ZX() << " "
1385  << std::setw(w) << fFluxRotation->ZY() << " "
1386  << std::setw(w) << fFluxRotation->ZZ() << " ]\n";
1387  atmoCfgText.precision(old_p);
1388  }
1389  mf::LogInfo("GENIEHelper") << atmoCfgText.str();
1390 
1391  atmo_flux_driver->LoadFluxData();
1392 
1393  // configure flux generation surface:
1394  atmo_flux_driver->SetRadii(fAtmoRl, fAtmoRt);
1395 
1396  fFluxD = atmo_flux_driver;//dynamic_cast<genie::GFluxI *>(atmo_flux_driver);
1397 
1398  } //end if using atmospheric fluxes
1399 
1400  if ( ! fFluxD ) {
1401  mf::LogError("GENIEHelper")
1402  << "Failed to contruct base flux driver for FluxType '"
1403  << fFluxType << "'";
1404  throw cet::exception("GENIEHelper")
1405  << "Failed to contruct base flux driver for FluxType '"
1406  << fFluxType << "'\n"
1407  << __FILE__ << ":" << __LINE__ << "\n";
1408  }
1409 
1410  //
1411  // Is the user asking to do flavor mixing?
1412  //
1413  fFluxD2GMCJD = fFluxD; // default: genie's GMCJDriver uses the bare flux generator
1414  if( fMixerConfig.find_first_not_of(" \t\n") != 0) // trim any leading whitespace
1415  fMixerConfig.erase( 0, fMixerConfig.find_first_not_of(" \t\n") );
1416  std::string keyword = fMixerConfig.substr(0,fMixerConfig.find_first_of(" \t\n"));
1417  if ( keyword != "none" && keyword != "" ) {
1418  // Wrap the true flux driver up in the adapter to allow flavor mixing
1419  genie::flux::GFlavorMixerI* mixer = 0;
1420  // here is where we map MixerConfig string keyword to actual class
1421  // first is a special case that is part of GENIE proper
1422  if ( keyword == "map" || keyword == "swap" || keyword == "fixedfrac" )
1423  mixer = new genie::flux::GFlavorMap();
1424  // if it wasn't one of the predefined known mixers then
1425  // see if the factory knows about it and can create one
1426  // assuming the keyword (first token) is the class name
1427  if ( ! mixer ) {
1428  genie::flux::GFlavorMixerFactory& mixerFactory =
1429  genie::flux::GFlavorMixerFactory::Instance();
1430  mixer = mixerFactory.GetFlavorMixer(keyword);
1431  if ( mixer ) {
1432  // remove class name from config string
1433  fMixerConfig.erase(0,keyword.size());
1434  // trim any leading whitespace
1435  if ( fMixerConfig.find_first_not_of(" \t\n") != 0 )
1436  fMixerConfig.erase( 0, fMixerConfig.find_first_not_of(" \t\n") );
1437  } else {
1438  const std::vector<std::string>& knownMixers =
1439  mixerFactory.AvailableFlavorMixers();
1440  std::ostringstream mixers;
1441  for (unsigned int j=0; j < knownMixers.size(); ++j ) {
1442  mixers << "\n [" << std::setw(2) << j << "] " << knownMixers[j];
1443  }
1444  mf::LogWarning("GENIEHelper")
1445  << " GFlavorMixerFactory known mixers: " << mixers.str();
1446  }
1447  }
1448  // configure the mixer
1449  if ( mixer ) mixer->Config(fMixerConfig);
1450  else {
1451  mf::LogWarning("GENIEHelper")
1452  << "GENIEHelper MixerConfig keyword was \"" << keyword
1453  << "\" but that did not map to a class; " << std::endl
1454  << "GFluxBlender in use, but no mixer";
1455  }
1456 
1457  genie::GFluxI* realFluxD = fFluxD;
1458  genie::flux::GFluxBlender* blender = new genie::flux::GFluxBlender();
1459  blender->SetBaselineDist(fMixerBaseline);
1460  blender->AdoptFluxGenerator(realFluxD);
1461  blender->AdoptFlavorMixer(mixer);
1462  fFluxD2GMCJD = blender;
1463  if ( fDebugFlags & 0x01 ) {
1464  if ( mixer ) mixer->PrintConfig();
1465  blender->PrintConfig();
1466  std::cout << std::flush;
1467  }
1468  }
1469 
1470  return;
1471  }
1472 
1473  //--------------------------------------------------
1474  void GENIEHelper::ConfigGeomScan()
1475  {
1476 
1477  // trim any leading whitespace
1478  if( fGeomScan.find_first_not_of(" \t\n") != 0)
1479  fGeomScan.erase( 0, fGeomScan.find_first_not_of(" \t\n") );
1480 
1481  if ( fGeomScan.find("default") == 0 ) return;
1482 
1483  genie::geometry::ROOTGeomAnalyzer* rgeom =
1484  dynamic_cast<genie::geometry::ROOTGeomAnalyzer*>(fGeomD);
1485 
1486  if ( !rgeom ) {
1487  throw cet::exception("GENIEHelper") << "fGeomD wasn't of type "
1488  << "genie::geometry::ROOTGeomAnalyzer*";
1489  }
1490 
1491  // convert string to lowercase
1492 
1493  // parse out string
1494  vector<string> strtok = genie::utils::str::Split(fGeomScan," ");
1495  // first value is a string, others should be numbers unless "file:"
1496  string scanmethod = strtok[0];
1497 
1498  // convert key string to lowercase (but not potential file name)
1499  std::transform(scanmethod.begin(),scanmethod.end(),scanmethod.begin(),::tolower);
1500 
1501  if ( scanmethod.find("file") == 0 ) {
1502  // xml expand path before passing in
1503  string filename = strtok[1];
1504  string fullname = genie::utils::xml::GetXMLFilePath(filename);
1505  mf::LogInfo("GENIEHelper")
1506  << "ConfigGeomScan getting MaxPathLengths from \"" << fullname << "\"";
1507  fDriver->UseMaxPathLengths(fullname);
1508  return;
1509  }
1510 
1511  vector<double> vals;
1512  for ( size_t indx=1; indx < strtok.size(); ++indx ) {
1513  const string& valstr1 = strtok[indx];
1514  if ( valstr1 != "" ) vals.push_back(atof(valstr1.c_str()));
1515  }
1516  size_t nvals = vals.size();
1517  // pad it out to at least 4 entries to avoid index issues
1518  for ( size_t nadd = 0; nadd < 4-nvals; ++nadd ) vals.push_back(0);
1519 
1520  double safetyfactor = 0;
1521  int writeout = 0;
1522  if ( scanmethod.find("box") == 0 ) {
1523  // use box method
1524  int np = (int)vals[0];
1525  int nr = (int)vals[1];
1526  if ( nvals >= 3 ) safetyfactor = vals[2];
1527  if ( nvals >= 4 ) writeout = vals[3];
1528  // protect against too small values
1529  if ( np <= 10 ) np = rgeom->ScannerNPoints();
1530  if ( nr <= 10 ) nr = rgeom->ScannerNRays();
1531  mf::LogInfo("GENIEHelper")
1532  << "ConfigGeomScan scan using box " << np << " points, "
1533  << nr << " rays";
1534  rgeom->SetScannerNPoints(np);
1535  rgeom->SetScannerNRays(nr);
1536  } else if ( scanmethod.find("flux") == 0 ) {
1537  // use flux method
1538  int np = (int)vals[0];
1539  if ( nvals >= 2 ) safetyfactor = vals[1];
1540  if ( nvals >= 3 ) writeout = vals[2];
1541  // sanity check, need some number of rays to explore the geometry
1542  // negative now means adjust rays to enu_max (GENIE svn 3676)
1543  if ( abs(np) <= 100 ) {
1544  int npnew = rgeom->ScannerNParticles();
1545  if ( np < 0 ) npnew = -abs(npnew);
1546  mf::LogWarning("GENIEHelper")
1547  << "Too few rays requested for geometry scan: " << np
1548  << ", use: " << npnew << "instead";
1549  np = npnew;
1550  }
1551  mf::LogInfo("GENIEHelper")
1552  << "ConfigGeomScan scan using " << np << " flux particles"
1553  << ( (np>0) ? "" : " with ray energy pushed to flux driver maximum" );
1554  rgeom->SetScannerFlux(fFluxD);
1555  rgeom->SetScannerNParticles(np);
1556  }
1557  else{
1558  // unknown
1559  throw cet::exception("GENIEHelper") << "fGeomScan unknown method: \""
1560  << fGeomScan << "\"";
1561  }
1562  if ( safetyfactor > 0 ) {
1563  mf::LogInfo("GENIEHelper")
1564  << "ConfigGeomScan setting safety factor to " << safetyfactor;
1565  rgeom->SetMaxPlSafetyFactor(safetyfactor);
1566  }
1567  if ( writeout != 0 ) SetMaxPathOutInfo();
1568  }
1569 
1570  //--------------------------------------------------
1571  void GENIEHelper::SetMaxPathOutInfo()
1572  {
1573  // create an info string based on:
1574  // ROOT geometry, TopVolume, FiducialCut, GeomScan, Flux
1575 
1576  mf::LogInfo("GENIEHelper") << "about to create MaxPathOutInfo";
1577 
1578  fMaxPathOutInfo = "\n";
1579  fMaxPathOutInfo += " FluxType: " + fFluxType + "\n";
1580  fMaxPathOutInfo += " BeamName: " + fBeamName + "\n";
1581  fMaxPathOutInfo += " FluxFiles: ";
1582  std::vector<string>::iterator ffitr = fSelectedFluxFiles.begin();
1583  for ( ; ffitr != fSelectedFluxFiles.end() ; ++ffitr )
1584  fMaxPathOutInfo += "\n " + *ffitr;
1585  fMaxPathOutInfo += "\n";
1586  fMaxPathOutInfo += " DetLocation: " + fDetLocation + "\n";
1587  fMaxPathOutInfo += " ROOTFile: " + fGeoFile + "\n";
1588  fMaxPathOutInfo += " WorldVolume: " + fWorldVolume + "\n";
1589  fMaxPathOutInfo += " TopVolume: " + fTopVolume + "\n";
1590  fMaxPathOutInfo += " FiducialCut: " + fFiducialCut + "\n";
1591  fMaxPathOutInfo += " GeomScan: " + fGeomScan + "\n";
1592 
1593  mf::LogInfo("GENIEHelper") << "MaxPathOutInfo: \""
1594  << fMaxPathOutInfo << "\"";
1595 
1596  }
1597 
1598  //--------------------------------------------------
1599  bool GENIEHelper::Stop()
1600  {
1601  // std::cout << "in GENIEHelper::Stop(), fEventsPerSpill = " << fEventsPerSpill
1602  // << " fPOTPerSpill = " << fPOTPerSpill << " fSpillExposure = " << fSpillExposure
1603  // << " fSpillEvents = " << fSpillEvents
1604  // << " fHistEventsPerSpill = " << fHistEventsPerSpill << std::endl;
1605 
1606  // determine if we should keep throwing neutrinos for
1607  // this spill or move on
1608 
1609  if ( fEventsPerSpill > 0 ) {
1610  if ( fSpillEvents < fEventsPerSpill ) return false;
1611  } else {
1612  // exposure (POT) based
1613  if ( fFluxType.find("tree_") == 0 ) {
1614  if ( fSpillExposure < fPOTPerSpill ) {
1615  return false;
1616  }
1617  }
1618  else if ( fFluxType.find("histogram") == 0 ) {
1619  if ( fSpillEvents < fHistEventsPerSpill ) return false;
1620  else fSpillExposure = fPOTPerSpill;
1621  }
1622  }
1623 
1624  if ( fFluxType.find("atmo_") == 0 ) {
1625  // the exposure for atmo is in SECONDS. In order to get seconds,
1626  // it needs to be normalized by 1e4 to take into account the units
1627  // discrepency between AtmoFluxDriver(/m2) and Generate(/cm2)
1628  // and it need to be normalized by the generation surface area since
1629  // it's not taken into accoutn in the flux driver
1630 
1631  long int nNeutrinos;
1632 
1633  if(fFluxType.find("PowerSpectrum") != std::string::npos){
1634  nNeutrinos = dynamic_cast<genie::flux::GPowerSpectrumAtmoFlux *>(fFluxD)->NFluxNeutrinos();
1635  fTotalExposure = nNeutrinos/fGenFlavors.size()/fDriver->GlobProbScale();
1636  mf::LogInfo("GENIEHelper")
1637  << "===> Atmo Pscale*Ngen/Nflavours = " << fTotalExposure;
1638  }
1639  else{
1640  nNeutrinos = dynamic_cast<genie::flux::GAtmoFlux *>(fFluxD)->NFluxNeutrinos();
1641  fTotalExposure = nNeutrinos * 1.0e4 / (TMath::Pi() * fAtmoRt*fAtmoRt);
1642  mf::LogInfo("GENIEHelper")
1643  << "===> Atmo EXPOSURE = " << fTotalExposure << " seconds.";
1644  }
1645 
1646 
1647 
1648 
1649 
1650  } else {
1651  fTotalExposure += fSpillExposure;
1652  }
1653 
1654  // made it to here, means need to reset the counters
1655  fSpillEvents = 0;
1656  fSpillExposure = 0.;
1657  fHistEventsPerSpill = fHelperRandom->Poisson(fXSecMassPOT*fTotalHistFlux);
1658  return true;
1659  }
1660 
1661  //--------------------------------------------------
1662  bool GENIEHelper::Sample(simb::MCTruth &truth, simb::MCFlux &flux, simb::GTruth &gtruth)
1663  {
1664  // set the top volume for the geometry
1665  fGeoManager->SetTopVolume(fGeoManager->FindVolumeFast(fTopVolume.c_str()));
1666 
1667  if ( fGenieEventRecord ) delete fGenieEventRecord;
1668 
1669  // ART Framework plays games with gRandom, undo that if requested
1670  TRandom* old_gRandom = gRandom;
1671  if (fUseHelperRndGen4GENIE) gRandom = fHelperRandom;
1672 
1673  fGenieEventRecord = fDriver->GenerateEvent();
1674 
1675  if (fForceApplyFlxWgt) {
1676  //const issue: double flxweight = fDriver->FluxDriver().Weight();
1677  //genie::GFluxI& flx = const_cast<genie::GFluxI&>(fDriver->FluxDriver());
1678  //double flxweight = flx.Weight();
1679  // use the flux driver handle we already have, rather than fetch const one
1680  // from the GMCJDriver fDriver
1681  double flxweight = fFluxD->Weight();
1682  double curweight = fGenieEventRecord->Weight();
1683  fGenieEventRecord->SetWeight(flxweight*curweight);
1684  }
1685 
1686  if (fUseHelperRndGen4GENIE) gRandom = old_gRandom;
1687 
1688  // now check if we produced a viable event record
1689  bool viableInteraction = true;
1690  if ( ! fGenieEventRecord ) viableInteraction = false;
1691 
1692  // update the spill total information, then check to see
1693  // if we got an event record that was valid
1694 
1695  genie::flux::GFluxExposureI* fexposure =
1696  dynamic_cast<genie::flux::GFluxExposureI*>(fFluxD);
1697  if ( fexposure ) {
1698  fSpillExposure =
1699  (fexposure->GetTotalExposure()/fDriver->GlobProbScale()) - fTotalExposure;
1700  }
1701  // use GENIE2ART code to fill MCFlux
1702  evgb::FillMCFlux(fFluxD,flux);
1703 
1704  // if no interaction generated return false
1705  if ( ! viableInteraction ) return false;
1706 
1707  // fill the MCTruth & GTruth information as we have a good interaction
1708  // these two objects are enough to reconstruct the GENIE record
1709  // use the new external functions in GENIE2ART
1710 
1711  // choose a time within the spill (ns) to shift the vertex times by:
1712  double timeoffset = 0;
1713  if ( ! fTimeShifter ) {
1714  timeoffset = fHelperRandom->Uniform()*fRandomTimeOffset;
1715  } else {
1716  timeoffset = fTimeShifter->TimeOffset();
1717  }
1718  // mf::LogInfo("GENIEHelper") << "TimeShifter adding " << timeoffset;
1719  double spilltime = fGlobalTimeOffset + timeoffset;
1720 
1721  std::unordered_map<std::string, std::string> genConfig;
1722 
1723  if(fFluxType.find("PowerSpectrum") != std::string::npos){
1724  genConfig.emplace("SpectralIndex", std::to_string(fAtmoSpectralIndex));
1725  }
1726 
1727  evgb::FillMCTruth(fGenieEventRecord, spilltime, truth,
1728  __GENIE_RELEASE__, fTuneName, fAddGenieVtxTime );
1729  evgb::FillGTruth(fGenieEventRecord, gtruth);
1730 
1731  // check to see if we are using flux ntuples but want to
1732  // make n events per spill
1733  if ( ( fEventsPerSpill > 0 ) &&
1734  ( fFluxType.find("tree_") == 0 ) ) ++fSpillEvents;
1735 
1736  // now check if using either histogram or mono fluxes, using
1737  // either n events per spill or basing events on POT per spill for the
1738  // histogram case
1739  if ( fFluxType.find("histogram") == 0 ) {
1740  // set the flag in the parent object that says the
1741  // fluxes came from histograms and fill related values
1742  flux.fFluxType = simb::kHistPlusFocus;
1743 
1744  // save the fluxes - fluxes were added to the vector in the same
1745  // order that the flavors appear in fGenFlavors
1746  int ctr = 0;
1747  int bin = fFluxHistograms[0]->FindBin(truth.GetNeutrino().Nu().E());
1748  std::vector<double> fluxes(6, 0.);
1749  for(std::vector<int>::iterator i = fGenFlavors.begin(); i != fGenFlavors.end(); i++){
1750  if(*i == 12) fluxes[kNue] = fFluxHistograms[ctr]->GetBinContent(bin);
1751  if(*i == -12) fluxes[kNueBar] = fFluxHistograms[ctr]->GetBinContent(bin);
1752  if(*i == 14) fluxes[kNuMu] = fFluxHistograms[ctr]->GetBinContent(bin);
1753  if(*i == -14) fluxes[kNuMuBar] = fFluxHistograms[ctr]->GetBinContent(bin);
1754  if(*i == 16) fluxes[kNuTau] = fFluxHistograms[ctr]->GetBinContent(bin);
1755  if(*i == -16) fluxes[kNuTauBar] = fFluxHistograms[ctr]->GetBinContent(bin);
1756  ++ctr;
1757  }
1758 
1759  // get the flux for each neutrino flavor of this energy
1760  flux.SetFluxGen(fluxes[kNue], fluxes[kNueBar],
1761  fluxes[kNuMu], fluxes[kNuMuBar],
1762  fluxes[kNuTau], fluxes[kNuTauBar]);
1763 
1764  ++fSpillEvents;
1765  }
1766  else if ( fFluxType.find("mono") == 0 ||
1767  fFluxType.find("function") == 0 ){
1768  ++fSpillEvents;
1769  }
1770  else if ( fFluxType.find("atmo_FLUKA") == 0 ||
1771  fFluxType.find("atmo_BARTOL") == 0 ||
1772  fFluxType.find("atmo_PowerSpectrum") == 0 ||
1773  fFluxType.find("atmo_BGLRS") == 0 ||
1774  fFluxType.find("atmo_HAKKM") == 0 ||
1775  fFluxType.find("atmo_HONDA") == 0 ) {
1776  if ( fEventsPerSpill > 0 ) ++fSpillEvents;
1777  flux.fFluxType = simb::kHistPlusFocus;
1778  }
1779 
1780 
1781  // fill these after the Pack[NuMI|Simple]Flux because those
1782  // will Reset() the values at the start
1783  TLorentzVector *vertex = fGenieEventRecord->Vertex();
1784  TLorentzVector nuray_pos = fFluxD->Position();
1785  TVector3 ray2vtx = nuray_pos.Vect() - vertex->Vect();
1786  flux.fgenx = nuray_pos.X();
1787  flux.fgeny = nuray_pos.Y();
1788  flux.fgenz = nuray_pos.Z();
1789  flux.fgen2vtx = ray2vtx.Mag();
1790 
1791  genie::flux::GFluxBlender* blender =
1792  dynamic_cast<genie::flux::GFluxBlender*>(fFluxD2GMCJD);
1793  if ( blender ) {
1794  if ( fUseBlenderDist ) flux.fdk2gen = blender->TravelDist();
1795  // / if mixing flavors print the state of the blender
1796  if ( fDebugFlags & 0x02 ) blender->PrintState();
1797  }
1798 
1799  if ( fDebugFlags & 0x04 ) {
1800  mf::LogInfo("GENIEHelper") << "vertex loc " << vertex->X() << ","
1801  << vertex->Y() << "," << vertex->Z() << std::endl
1802  << " flux ray start " << nuray_pos.X() << ","
1803  << nuray_pos.Y() << "," << nuray_pos.Z() << std::endl
1804  << " ray2vtx = " << flux.fgen2vtx
1805  << " dk2ray = " << flux.fdk2gen;
1806  }
1807  if ( fGHepPrintLevel >= 0 ) {
1808  std::cout << *fGenieEventRecord;
1809  }
1810 
1811  // set the top volume of the geometry back to the world volume
1812  fGeoManager->SetTopVolume(fGeoManager->FindVolumeFast(fWorldVolume.c_str()));
1813 
1814  return true;
1815  }
1816 
1817  //---------------------------------------------------------
1818  void GENIEHelper::BuildFluxRotation()
1819  {
1820  // construct fFluxRotation matrix
1821  // from fFluxRotCfg + fFluxRotValues
1822  if ( fFluxRotCfg == "" ||
1823  ( fFluxRotCfg.find("none") != std::string::npos ) ) return;
1824 
1825  size_t nval = fFluxRotValues.size();
1826 
1827  bool verbose = ( fFluxRotCfg.find("verbose") != std::string::npos );
1828  if (verbose) {
1829  std::ostringstream indata;
1830  indata << "BuildFluxRotation: Cfg \"" << fFluxRotCfg << "\"\n"
1831  << " " << nval << " values\n";
1832  for (size_t i=0; i<nval; ++i) {
1833  indata << " [" << std::setw(2) << i << "] " << fFluxRotValues[i] << "\n";
1834  }
1835  mf::LogInfo("GENIEHelper") << indata.str();
1836  }
1837 
1838  // interpret as a full 3x3 array
1839  if ( fFluxRotCfg.find("newxyz") != std::string::npos ||
1840  fFluxRotCfg.find("3x3") != std::string::npos ) {
1841  if ( nval == 9 ) {
1842  TRotation fTempRot;
1843  TVector3 newX = TVector3(fFluxRotValues[0],
1844  fFluxRotValues[1],
1845  fFluxRotValues[2]);
1846  TVector3 newY = TVector3(fFluxRotValues[3],
1847  fFluxRotValues[4],
1848  fFluxRotValues[5]);
1849  TVector3 newZ = TVector3(fFluxRotValues[6],
1850  fFluxRotValues[7],
1851  fFluxRotValues[8]);
1852  fTempRot.RotateAxes(newX,newY,newZ);
1853  // weirdly necessary; frame vs. obj rotation
1854  fFluxRotation = new TRotation(fTempRot.Inverse());
1855  return;
1856  } else {
1857  throw cet::exception("BadFluxRotation")
1858  << "specified: " << fFluxRotCfg << "\n"
1859  << " but nval=" << nval << ", need 9";
1860  }
1861  }
1862 
1863  // another possibility ... series of rotations around particular axes
1864  if ( fFluxRotCfg.find("series") != std::string::npos ) {
1865  TRotation fTempRot;
1866  // if series then cfg should also have series of rot{X|Y|Z}{deg|rad}
1867  std::vector<std::string> strs = genie::utils::str::Split(fFluxRotCfg," ,;(){}[]");
1868  size_t nrot = -1;
1869  for (size_t j=0; j<strs.size(); ++j) {
1870  std::string what = strs[j];
1871  if ( what == "" ) continue;
1872  // lower case for convenience
1873  std::transform(what.begin(),what.end(),what.begin(),::tolower);
1874  if ( what == "series" ) continue;
1875  if ( what == "verbose" ) continue;
1876  if ( what.find("rot") != 0 ) {
1877  mf::LogWarning("GENIEHelper")
1878  << "processing series rotation saw keyword \"" << what << "\" -- ignoring";
1879  continue;
1880  }
1881  char axis = what[3];
1882  // check that axis is sensibly x, y or z
1883  if ( axis != 'x' && axis != 'y' && axis != 'z' ) {
1884  throw cet::exception("BadFluxRotation")
1885  << "specified: " << fFluxRotCfg << "\n"
1886  << " keyword '" << what << "': bad axis '" << axis << "'";
1887  }
1888  std::string units = what.substr(4);
1889  // don't worry if written fully as "radians" or "degrees"
1890  if (units.size() > 3) units.erase(3);
1891  if ( units != "" && units != "rad" && units != "deg" ) {
1892  throw cet::exception("BadFluxRotation")
1893  << "specified: " << fFluxRotCfg << "\n"
1894  << " keyword '" << what << "': bad units '" << units << "'";
1895  }
1896  // no units? assume degrees
1897  double scale = (( units == "rad" ) ? 1.0 : TMath::DegToRad() );
1898 
1899  ++nrot;
1900  if ( nrot >= nval ) {
1901  // not enough values
1902  throw cet::exception("BadFluxRotation")
1903  << "specified: " << fFluxRotCfg << "\n"
1904  << " asking for rotation [" << nrot << "] "
1905  << what << " but nval=" << nval;
1906  }
1907  double rot = scale * fFluxRotValues[nrot];
1908  if ( axis == 'x' ) fTempRot.RotateX(rot);
1909  if ( axis == 'y' ) fTempRot.RotateY(rot);
1910  if ( axis == 'z' ) fTempRot.RotateZ(rot);
1911 
1912  } // loop over tokens in cfg string
1913 
1914  // weirdly necessary; frame vs. obj rotation
1915  fFluxRotation = new TRotation(fTempRot.Inverse());
1916 
1917  if ( nrot+1 != nval ) {
1918  // call out possible user mistake
1919  mf::LogWarning("GENIEHelper")
1920  << "BuildFluxRotation only used " << nrot+1 << " of "
1921  << nval << " FluxRotValues";
1922  }
1923  return;
1924  }
1925 
1926  // could put other interpretations here ...
1927 
1928  throw cet::exception("BadFluxRotation")
1929  << "specified: " << fFluxRotCfg << "\n"
1930  << " nval=" << nval << ", but don't know how to interpret that";
1931 
1932  }
1933  //---------------------------------------------------------
1934  void GENIEHelper::ExpandFluxPaths()
1935  {
1936  // expand any wildcards in the paths variable
1937  // if unset and using the old DIRECT method allow it to fall back
1938  // to using FW_SEARCH_PATH ... but not for the new ifdhc approach
1939 
1940  std::string initial = fFluxSearchPaths;
1941 
1942  if ( fFluxCopyMethod == "DIRECT" && fFluxSearchPaths == "" ) {
1943  fFluxSearchPaths = cet::getenv("FW_SEARCH_PATH");
1944  }
1945  fFluxSearchPaths = gSystem->ExpandPathName(fFluxSearchPaths.c_str());
1946 
1947  mf::LogInfo("GENIEHelper")
1948  << "ExpandFluxPaths initially: \"" << initial << "\"\n"
1949  << " final result: \"" << fFluxSearchPaths << "\"\n"
1950  << " using: \"" << fFluxCopyMethod << "\" method";
1951  }
1952  //---------------------------------------------------------
1953  void GENIEHelper::ExpandFluxFilePatternsDirect()
1954  {
1955  // Using the the fFluxSearchPaths list of directories, apply the
1956  // user supplied pattern as a suffix to find the flux files.
1957  // The userpattern might include simple wildcard globs (in contrast
1958  // to proper regexp patterns).
1959 
1960  // After expanding the list to individual files, randomize them
1961  // and start selecting until a size limit is about to be
1962  // exceeded (though a minimum there needs to be one file, not
1963  // matter the limit).
1964 
1965  // older GENIE (pre r3669) can't handle vectors/sets of file names
1966  // but an only accept a pattern that will resolve to files. This
1967  // collection could exceed the desired size threshold, but for
1968  // pick the collection that expands to the largest list
1969 
1970  bool randomizeFiles = false;
1971  if ( fFluxType.find("tree_") == 0 ) randomizeFiles = true;
1972 
1973  std::vector<std::string> dirs;
1974  cet::split_path(fFluxSearchPaths,dirs);
1975  if ( dirs.empty() ) dirs.push_back(std::string()); // at least null string
1976 
1977  glob_t g;
1978  int flags = GLOB_TILDE; // expand ~ home directories
1979 
1980  std::ostringstream patterntext; // for info/error messages
1981  std::ostringstream dirstext; // for info/error messages
1982 
1983  std::vector<std::string>::const_iterator uitr = fFluxFilePatterns.begin();
1984  int ipatt = 0;
1985 
1986  for ( ; uitr != fFluxFilePatterns.end(); ++uitr, ++ipatt ) {
1987  std::string userpattern = *uitr;
1988  patterntext << "\n\t" << userpattern;
1989 
1990  std::vector<std::string>::const_iterator ditr = dirs.begin();
1991  for ( ; ditr != dirs.end(); ++ditr ) {
1992  std::string dalt = *ditr;
1993  // if non-null, does it end with a "/"? if not add one
1994  size_t len = dalt.size();
1995  if ( len > 0 && dalt.rfind('/') != len-1 ) dalt.append("/");
1996  if ( uitr == fFluxFilePatterns.begin() ) dirstext << "\n\t" << dalt;
1997 
1998  std::string filepatt = dalt + userpattern;
1999 
2000  glob(filepatt.c_str(),flags,NULL,&g);
2001  if ( g.gl_pathc > 0 ) flags |= GLOB_APPEND; // next glob() will append to list
2002 
2003  } // loop over FluxSearchPaths dirs
2004  } // loop over user patterns
2005 
2006  std::ostringstream paretext;
2007  std::ostringstream flisttext;
2008 
2009  int nfiles = g.gl_pathc;
2010 
2011  if ( nfiles == 0 ) {
2012  paretext << "\n expansion resulted in a null list for flux files";
2013 
2014  } else if ( ! randomizeFiles ) {
2015  // some sets of files should be left in order
2016  // and no size limitations imposed ... just copy the list
2017 
2018  paretext << "\n list of files will be processed in order";
2019  for (int i=0; i<nfiles; ++i) {
2020  std::string afile(g.gl_pathv[i]);
2021  fSelectedFluxFiles.push_back(afile);
2022 
2023  flisttext << "[" << setw(3) << i << "] "
2024  << afile << "\n";
2025  }
2026  } else {
2027 
2028  // now pull from the list randomly
2029  // do this by assigning a random number to each;
2030  // ordering that list; and pulling in that order
2031 
2032  paretext << "list of " << nfiles
2033  << " will be randomized and pared down to "
2034  << fMaxFluxFileMB << " MB or "
2035  << fMaxFluxFileNumber << " files";
2036 
2037  double* order = new double[nfiles];
2038  int* indices = new int[nfiles];
2039  fHelperRandom->RndmArray(nfiles,order);
2040  // assign random # for their relative order
2041 
2042  TMath::Sort((int)nfiles,order,indices,false);
2043 
2044  long long int sumBytes = 0; // accumulated size in bytes
2045  long long int maxBytes = (long long int)fMaxFluxFileMB * 1024ll * 1024ll;
2046 
2047  FileStat_t fstat;
2048  for (int i=0; i < TMath::Min(nfiles,fMaxFluxFileNumber); ++i) {
2049  int indx = indices[i];
2050  std::string afile(g.gl_pathv[indx]);
2051  bool keep = true;
2052 
2053  gSystem->GetPathInfo(afile.c_str(),fstat);
2054  sumBytes += fstat.fSize;
2055  // skip those that would push sum above total
2056  // but always accept at least one (the first)
2057  if ( sumBytes > maxBytes && i != 0 ) keep = false;
2058 
2059  flisttext << "[" << setw(3) << i << "] "
2060  << "=> g[" << setw(3) << indx << "] "
2061  << ((keep)?"keep":"skip") << " "
2062  << setw(6) << (sumBytes/(1024ll*1024ll)) << " "
2063  << afile << "\n";
2064 
2065  if ( keep ) fSelectedFluxFiles.push_back(afile);
2066  else break; // <voice name=Scotty> Captain, she can't take any more</voice>
2067 
2068  }
2069  delete [] order;
2070  delete [] indices;
2071 
2072  }
2073 
2074  mf::LogInfo("GENIEHelper")
2075  << "ExpandFluxFilePatternsDirect initially found " << nfiles
2076  << " files for user patterns:"
2077  << patterntext.str() << "\n using FluxSearchPaths of: "
2078  << dirstext.str() << "\n" << paretext.str();
2079  //<< "\"" << cet::getenv("FW_SEARCH_PATH") << "\"";
2080 
2081  mf::LogDebug("GENIEHelper") << "\n" << flisttext.str();
2082 
2083  // done with glob list
2084  globfree(&g);
2085 
2086  // no null path allowed for at least these
2087  if ( fFluxType.find("tree_") == 0 ) {
2088  size_t nfiles = fSelectedFluxFiles.size();
2089  if ( nfiles == 0 ) {
2090  mf::LogError("GENIEHelper")
2091  << "For \"" << fFluxType <<"\" "
2092  << "(e.g. \"dk2nu\', \"ntuple\" (\"numi\") or \"simple\") "
2093  << " specification must resolve to at least one file"
2094  << "\n none were found. DIRECT user pattern(s): "
2095  << patterntext.str()
2096  << "\n using FluxSearchPaths of: "
2097  << dirstext.str();
2098 
2099  throw cet::exception("NoFluxFiles")
2100  << "no flux files found for: " << patterntext.str() << "\n"
2101  << " in: " << dirstext.str();
2102 
2103  }
2104  }
2105 
2106  } // ExpandFluxFilePatternsDirect
2107 
2108  //---------------------------------------------------------
2109  void GENIEHelper::ExpandFluxFilePatternsIFDH()
2110  {
2111  // Using the the FluxSearchPaths list of directories, apply the
2112  // user supplied pattern as a suffix to find the flux files.
2113  // The userpattern might include simple wildcard globs (in contrast
2114  // to proper regexp patterns).
2115 
2116  // After expanding the list to individual files, randomize them
2117  // and start selecting until a size limit is about to be
2118  // exceeded (though a minimum there needs to be one file, not
2119  // matter the limit).
2120 
2121  // Use the IFDH interface to get the list of files and sizes;
2122  // after sorting/selecting use IFDH to make a local copy
2123 
2124 #ifdef NO_IFDH_LIB
2125  std::ostringstream fmesg;
2126  std::string marker = "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n";
2127  fmesg << marker
2128  << __FILE__ << ":" << __LINE__
2129  << "\nno IFDH implemented on this platform\n"
2130  << marker;
2131  // make sure the message goes everywhere
2132  std::cout << fmesg.str() << std::flush;
2133  std::cerr << fmesg.str();
2134  throw cet::exception("Attempt to use ifdh class") << fmesg.str();
2135  assert(0);
2136 #else
2137  // if "method" just an identifier and not a scheme then clear it
2138  if ( fFluxCopyMethod.find("IFDH") == 0 ) fFluxCopyMethod = "";
2139 
2140  bool randomizeFiles = false;
2141  if ( fFluxType.find("tree_") == 0 ) randomizeFiles = true;
2142 
2143  #ifdef USE_IFDH_SERVICE
2144  art::ServiceHandle<IFDH> ifdhp;
2145  #else
2146  if ( ! fIFDH ) fIFDH = new ifdh_ns::ifdh;
2147  #endif
2148 
2149  std::string spaths = fFluxSearchPaths;
2150 
2151  const char* ifdh_debug_env = std::getenv("IFDH_DEBUG_LEVEL");
2152  if ( ifdh_debug_env ) {
2153  mf::LogInfo("GENIEHelper") << "IFDH_DEBUG_LEVEL: " << ifdh_debug_env;
2154  fIFDH->set_debug(ifdh_debug_env);
2155  }
2156 
2157  // filenames + size
2158  std::vector<std::pair<std::string,long>>
2159  partiallist, fulllist, selectedlist, locallist;
2160 
2161  std::ostringstream patterntext; // stringification of vector of patterns
2162  std::ostringstream fulltext; // for info on full list of matches
2163  std::ostringstream selectedtext; // for info on selected files
2164  std::ostringstream localtext; // for info on local files
2165  fulltext << "search paths: " << spaths;
2166 
2167  //std::vector<std::string>::const_iterator uitr = fFluxFilePatterns.begin();
2168 
2169  // loop over possible patterns
2170  // IFDH handles various stems but done a list of globs
2171  size_t ipatt=0;
2172  auto uitr = fFluxFilePatterns.begin();
2173  for ( ; uitr != fFluxFilePatterns.end(); ++uitr, ++ipatt ) {
2174  std::string userpattern = *uitr;
2175  patterntext << "\npattern [" << std::setw(3) << ipatt << "] " << userpattern;
2176  fulltext << "\npattern [" << std::setw(3) << ipatt << "] " << userpattern;
2177 
2178  #ifdef USE_IFDH_SERVICE
2179  partiallist = ifdhp->findMatchingFiles(spaths,userpattern);
2180  #else
2181  partiallist = fIFDH->findMatchingFiles(spaths,userpattern);
2182  #endif
2183  fulllist.insert(fulllist.end(),partiallist.begin(),partiallist.end());
2184 
2185  // make a complete list ...
2186  fulltext << " found " << partiallist.size() << " files";
2187  for (auto pitr = partiallist.begin(); pitr != partiallist.end(); ++pitr) {
2188  fulltext << "\n " << std::setw(10) << pitr->second << " " << pitr->first;
2189  }
2190 
2191  partiallist.clear();
2192  } // loop over user patterns
2193 
2194  size_t nfiles = fulllist.size();
2195 
2196  mf::LogInfo("GENIEHelper")
2197  << "ExpandFluxFilePatternsIFDH initially found " << nfiles << " files";
2198  mf::LogDebug("GENIEHelper")
2199  << fulltext.str();
2200 
2201  if ( nfiles == 0 ) {
2202  selectedtext << "\n expansion resulted in a null list for flux files";
2203  } else if ( ! randomizeFiles ) {
2204  // some sets of files should be left in order
2205  // and no size limitations imposed ... just copy the list
2206 
2207  selectedtext << "\n list of files will be processed in order";
2208  selectedlist.insert(selectedlist.end(),fulllist.begin(),fulllist.end());
2209 
2210  } else {
2211 
2212  // for list needing size based trimming and randomization ...
2213  // pull from the list randomly until a cummulative limit is reached
2214  // do this by assigning a random number to each;
2215  // ordering that list; and pulling in that order
2216 
2217  selectedtext << "list of " << nfiles
2218  << " will be randomized and pared down to "
2219  << fMaxFluxFileMB << " MB or "
2220  << fMaxFluxFileNumber << " files";
2221 
2222  double* order = new double[nfiles];
2223  int* indices = new int[nfiles];
2224  fHelperRandom->RndmArray(nfiles,order);
2225  // assign random # for their relative order
2226 
2227  TMath::Sort((int)nfiles,order,indices,false);
2228 
2229  long long int sumBytes = 0; // accumulated size in bytes
2230  long long int maxBytes = (long long int)fMaxFluxFileMB * 1024ll * 1024ll;
2231 
2232  for (size_t i=0; i < TMath::Min(nfiles,size_t(fMaxFluxFileNumber)); ++i) {
2233  int indx = indices[i];
2234  bool keep = true;
2235 
2236  auto p = fulllist[indx]; // this pair <name,size>
2237  sumBytes += p.second;
2238  // skip those that would push sum above total
2239  // but always accept at least one (the first)
2240  if ( sumBytes > maxBytes && i != 0 ) keep = false;
2241 
2242  selectedtext << "\n[" << setw(3) << i << "] "
2243  << "=> [" << setw(3) << indx << "] "
2244  << ((keep)?"keep":"SKIP") << " "
2245  << std::setw(6) << (sumBytes/(1024ll*1024ll)) << " MB "
2246  << p.first;
2247 
2248  if ( keep ) selectedlist.push_back(p);
2249  else break; // <voice name=Scotty> Captain, she can't take any more</voice>
2250 
2251  }
2252  delete [] order;
2253  delete [] indices;
2254  }
2255 
2256  mf::LogInfo("GENIEHelper")
2257  << selectedtext.str();
2258 
2259  // have a selected list of remote files
2260  // get paths to local copies
2261  #ifdef USE_IFDH_SERVICE
2262  locallist = ifdhp->fetchSharedFiles(selectedlist,fFluxCopyMethod);
2263  #else
2264  locallist = fIFDH->fetchSharedFiles(selectedlist,fFluxCopyMethod);
2265  #endif
2266 
2267  localtext << "final list of files:";
2268  size_t i=0;
2269  for (auto litr = locallist.begin(); litr != locallist.end(); ++litr, ++i) {
2270  fSelectedFluxFiles.push_back(litr->first);
2271  localtext << "\n\t[" << std::setw(3) << i << "]\t" << litr->first;
2272  }
2273 
2274  mf::LogInfo("GENIEHelper")
2275  << localtext.str();
2276 
2277  // no null path allowed for at least these
2278  if ( fFluxType.find("tree_") == 0 ) {
2279  size_t nfiles = fSelectedFluxFiles.size();
2280  if ( nfiles == 0 ) {
2281  mf::LogError("GENIEHelper")
2282  << "For \"" << fFluxType <<"\" "
2283  << "(e.g. \"dk2nu\', \"ntuple\" (\"numi\") or \"simple\") "
2284  << "specification must resolve to at least one file"
2285  << "\n none were found. IFDH user pattern(s): "
2286  << patterntext.str()
2287  << "\n using FluxSearchPaths of: "
2288  << spaths;
2289 
2290  throw cet::exception("NoFluxFiles")
2291  << "no flux files found for: " << patterntext.str() << "\n"
2292  << " in " << spaths;
2293 
2294  }
2295  }
2296 #endif // 'else' code only if NO_IFDH_LIB not defined
2297  } // ExpandFluxFilePatternsIFDH
2298 
2299  //---------------------------------------------------------
2300  void GENIEHelper::SetGXMLPATH()
2301  {
2304 
2305  // priority order
2306  // (fcl file paths):(existing user environment):(FW_SEARCH_PATH)
2307  //
2308  // fcl parameters might be the explicit GXMLPATH value
2309  // or a pair in the fEnvironment vector
2310  // but the final "official result" should the fGXMLPATH
2311 
2312  // start with fGXMLPATH set from pset "GXMLPATH" value
2313 
2314  // find it in the vector, if it exists
2315  int indxGXMLPATH = -1;
2316  for (size_t i = 0; i < fEnvironment.size(); i += 2) {
2317  if ( fEnvironment[i].find("GXMLPATH") == 0 ) {
2318  if ( fGXMLPATH != "" ) fGXMLPATH += ":";
2319  fGXMLPATH += fEnvironment[i+1];
2320  indxGXMLPATH = i;
2321  /*
2322  throw cet::exception("UsingGXMLPATH")
2323  << "using Environment fcl parameter GXMLPATH: "
2324  << fEnvironment[indxGXMLPATH+1]
2325  << ", use fcl parameter GXMLPATH instead.";
2326  */
2327  break;
2328  }
2329  }
2330 
2331  const char* gxmlpath_env = std::getenv("GXMLPATH");
2332  if ( gxmlpath_env ) {
2333  if ( fGXMLPATH != "" ) fGXMLPATH += ":";
2334  fGXMLPATH += gxmlpath_env;
2335  }
2336  const char* fwpath_env = std::getenv("FW_SEARCH_PATH");
2337  if ( fwpath_env ) {
2338  if ( fGXMLPATH != "" ) fGXMLPATH += ":";
2339  fGXMLPATH += fwpath_env;
2340  }
2341 
2342  // refresh fEnvironment (in case fEnvironment is still being used)
2343  if ( indxGXMLPATH < 0 ) {
2344  // nothing in fcl parameters
2345  indxGXMLPATH=fEnvironment.size();
2346  fEnvironment.push_back("GXMLPATH");
2347  fEnvironment.push_back(fGXMLPATH);
2348  } else {
2349  // update value
2350  fEnvironment[indxGXMLPATH+1] = fGXMLPATH;
2351  }
2352 
2353  // now set it externally for use by GENIE
2354  gSystem->Setenv("GXMLPATH",fGXMLPATH.c_str());
2355 
2356  }
2357 
2358  //---------------------------------------------------------
2359  void GENIEHelper::SetGMSGLAYOUT()
2360  {
2364 
2365  // start with GMSGLAYOUT set from pset "GMSGLAYOUT" value
2366 
2367  // find it in the vector, if it exists
2368  // this will override the top level fcl value
2369  for (size_t i = 0; i < fEnvironment.size(); i += 2) {
2370  if ( fEnvironment[i].find("GMSGLAYOUT") == 0 ) {
2371  fGMSGLAYOUT = fEnvironment[i+1];
2372  break;
2373  }
2374  }
2375 
2376  // now set it externally, if we have a value
2377  if ( fGMSGLAYOUT != "" ) {
2378  gSystem->Setenv("GMSGLAYOUT",fGMSGLAYOUT.c_str());
2379  }
2380  }
2381 
2382  //---------------------------------------------------------
2383  void GENIEHelper::StartGENIEMessenger(std::string prodmodestr)
2384  {
2389 
2392 
2393  int indxGPRODMODE = -1;
2394  int indxGMSGCONF = -1;
2395 
2396  for (size_t i = 0; i < fEnvironment.size(); i += 2) {
2397  if ( fEnvironment[i].find("GPRODMODE") == 0 ) {
2398  indxGPRODMODE = i;
2399  continue;
2400  }
2401  if ( fEnvironment[i].find("GMSGCONF") == 0 ) {
2402  indxGMSGCONF = i;
2403  continue;
2404  }
2405  }
2406  // GMSGCONF set in fEnvironment tack it on to current setting
2407  if ( indxGMSGCONF >= 0 ) {
2408  if ( fGENIEMsgThresholds != "" ) fGENIEMsgThresholds += ":";
2409  fGENIEMsgThresholds += fEnvironment[indxGMSGCONF+1];
2410  } else {
2411  // nothing in fcl parameters, make a placeholder
2412  indxGMSGCONF = fEnvironment.size();
2413  fEnvironment.push_back("GMSGCONF");
2414  fEnvironment.push_back("");
2415  }
2416 
2417  bool prodmode = StringToBool(prodmodestr);
2418  if ( indxGPRODMODE >= 0 ) {
2419  // user tried to set GPRODMODE via Environment fcl values
2420  prodmode |= StringToBool(fEnvironment[indxGPRODMODE+1]);
2421  }
2422 
2423  if ( prodmode ) {
2424  // okay they really meant it
2425  // PREpend "Messenger_whisper.xml" to existing value
2426  // don't forget the colon if necessary
2427  std::string newval = "Messenger_whisper.xml";
2428  if ( fGENIEMsgThresholds != "" ) {
2429  newval += ":";
2430  newval += fGENIEMsgThresholds;
2431  }
2432  fGENIEMsgThresholds = newval;
2433  }
2434 
2435  // update fEnvironment value
2436  if ( indxGMSGCONF >= 0 ) {
2437  fEnvironment[indxGMSGCONF+1] = fGENIEMsgThresholds;
2438  }
2439 
2440  mf::LogInfo("GENIEHelper")
2441  << "StartGENIEMessenger ProdMode=" << ((prodmode)?"yes":"no")
2442  << " read from: " << fGENIEMsgThresholds;
2443 
2444  genie::utils::app_init::MesgThresholds(fGENIEMsgThresholds);
2445 
2446  }
2447 
2448  //---------------------------------------------------------
2449  void GENIEHelper::ReadXSecTable()
2450  {
2453 
2454  // priority order:
2455  // fcl fEnvironment GSPLOAD
2456  // fcl XSecTable
2457  // $GSPLOAD in environment
2458  // default 'gxspl-FNALsmall.xml'
2459 
2460  if ( fXSecTable == "" ) {
2461  // stand-alone value is not set
2462  const char* gspload_alt = std::getenv("GSPLOAD");
2463  if ( ! gspload_alt ) {
2464  const char* gspload_dflt = "gxspl-FNALsmall.xml"; // fall back
2465  gspload_alt = gspload_dflt;
2466  } else {
2467  throw cet::exception("$GSPLOAD")
2468  << "using env variable $GSPLOAD: " << gspload_alt
2469  << ", use fcl parameter 'XSecTable' instead.";
2470  }
2471  fXSecTable = std::string(gspload_alt);
2472  }
2473 
2474  // find GSPLOAD in the vector, if it exists
2475  int indxGSPLOAD = -1;
2476  for (size_t i = 0; i < fEnvironment.size(); i += 2) {
2477  if ( fEnvironment[i].find("GSPLOAD") == 0 ) {
2478  indxGSPLOAD = i;
2479  throw cet::exception("UsingGSPLOAD")
2480  << "using Environment fcl parameter GSPLOAD: "
2481  << fEnvironment[indxGSPLOAD+1]
2482  << ", use fcl parameter 'XSecTable' instead. "
2483  << __FILE__ << ":" << __LINE__ << "\n";
2484  continue;
2485  }
2486  }
2487 
2488  if ( indxGSPLOAD < 0 ) {
2489  // nothing in fcl parameters
2490  indxGSPLOAD=fEnvironment.size();
2491  fEnvironment.push_back("GSPLOAD");
2492  fEnvironment.push_back(fXSecTable);
2493  } else {
2494  fXSecTable = fEnvironment[indxGSPLOAD+1]; // get two in sync
2495  }
2496 
2497  // currently GENIE doesn't internally use GXMLPATH when looking for
2498  // spline files, but instead wants a fully expanded path.
2499  // Do the expansion here using the extended GXMLPATH list
2500  // of locations (which included $FW_SEARCH_PATH)
2501  mf::LogDebug("GENIEHelper") << "GSPLOAD as originally: " << fXSecTable;
2502 
2503  // RWH cet::search_path returns "" if the input string is actually
2504  // the full path to the file .. this is not really what one wants,
2505  // one just wan't the full path to the file; seems to work if
2506  // "/" is made to be another possible PATH
2507  cet::search_path spGXML( "/:" + fGXMLPATH );
2508  std::string fullpath;
2509  spGXML.find_file(fXSecTable, fullpath);
2510 
2511  if ( fullpath == "" ) {
2512  mf::LogError("GENIEHelper")
2513  << "could not resolve full path for spline file XSecTable/GSPLOAD "
2514  << "\"" << fXSecTable << "\" using: "
2515  << fGXMLPATH;
2516  throw cet::exception("UnresolvedGSPLOAD")
2517  << "can't find XSecTable/GSPLOAD file: " << fXSecTable;
2518  }
2519  fXSecTable = fullpath;
2520  fEnvironment[indxGSPLOAD+1] = fXSecTable; // get two in sync
2521 
2522  mf::LogInfo("GENIEHelper")
2523  << "XSecTable/GSPLOAD full path \"" << fXSecTable << "\"";
2524 
2525  TStopwatch xtime;
2526  xtime.Start();
2527 
2528  // can't use gSystem->Unsetenv() as it is really gSystem->Setenv(name,"")
2529  unsetenv("GSPLOAD"); // MUST!!! ensure that it isn't set externally
2530  genie::utils::app_init::XSecTable(fXSecTable,true);
2531 
2532  xtime.Stop();
2533  mf::LogInfo("GENIEHelper")
2534  << "Time to read GENIE XSecTable: "
2535  << " Real " << xtime.RealTime() << " s,"
2536  << " CPU " << xtime.CpuTime() << " s"
2537  << " from " << fXSecTable;
2538 
2539  }
2540 
2541  //---------------------------------------------------------
2542  bool GENIEHelper::StringToBool(std::string v)
2543  {
2544  if (v == "true") return true; // C++ style
2545  if (v == "false") return false;
2546  if (v == "kTRUE") return true; // ROOT style
2547  if (v == "kFALSE") return false;
2548  if (v == "TRUE") return true; // Some other reasonable variations...
2549  if (v == "FALSE") return false;
2550  if (v == "True") return true;
2551  if (v == "False") return false;
2552  if (v == "on") return true;
2553  if (v == "off") return false;
2554  if (v == "On") return true;
2555  if (v == "Off") return false;
2556  if (v == "ON") return true;
2557  if (v == "OFF") return false;
2558  if (v == "YES") return true;
2559  if (v == "NO") return false;
2560  if (v == "Yes") return true;
2561  if (v == "No") return false;
2562  if (v == "yes") return true;
2563  if (v == "no") return false;
2564  if (v == "1") return true;
2565  if (v == "0") return false;
2566 
2567  return false; // by default invalid strings are false
2568  }
2569 
2570 } // namespace evgb
2571 
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:234
simb::MCFlux::fgenz
double fgenz
Definition: MCFlux.h:102
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:27
evgb::GENIEHelper::fFluxFilePatterns
std::vector< std::string > fFluxFilePatterns
wildcard patterns files containing histograms or ntuples, or txt
Definition: GENIEHelper.h:143
evgb::kNuTauBar
static const int kNuTauBar
Definition: GENIEHelper.cxx:211
evgb::GENIEHelper::fUseHelperRndGen4GENIE
bool fUseHelperRndGen4GENIE
use fHelperRandom for gRandom during Sample()
Definition: GENIEHelper.h:136
evgb::GENIEHelper::fEventsPerSpill
double fEventsPerSpill
Definition: GENIEHelper.h:160
simb::MCTruth::GetNeutrino
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
evgb::GENIEHelper::fDetLocation
std::string fDetLocation
name of flux window location
Definition: GENIEHelper.h:156
evgb::GENIEHelper::fAtmoEmax
double fAtmoEmax
atmo: Maximum energy of neutrinos in GeV
Definition: GENIEHelper.h:188
evgb::GENIEHelper::fMaxFluxFileNumber
int fMaxFluxFileNumber
maximum # of flux files
Definition: GENIEHelper.h:146
evgb::GENIEHelper::fMixerBaseline
double fMixerBaseline
baseline distance if genie flux can&#39;t calculate it
Definition: GENIEHelper.h:203
evgb::GENIEHelper::fTuneName
std::string fTuneName
GENIE R-3 Tune name (defines model configuration)
Definition: GENIEHelper.h:196
evgb::FillMCTruth
void FillMCTruth(const genie::EventRecord *grec, double spillTime, simb::MCTruth &mctruth, const std::string &genieVersion="unknown", const std::string &genieTune="unknown", bool addGenieVtxTime=false, const std::unordered_map< std::string, std::string > genConfig={})
Definition: GENIE2ART.cxx:182
simb::MCFlux::fgeny
double fgeny
Definition: MCFlux.h:101
evgb::GENIEHelper::fTotalHistFlux
double fTotalHistFlux
total flux of neutrinos from flux histograms for used flavors
Definition: GENIEHelper.h:173
simb::MCFlux::fgen2vtx
double fgen2vtx
distance from ray origin to event vtx
Definition: MCFlux.h:104
mf::MaybeLogger_
Definition: MessageLogger.h:113
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:209
evgb::GENIEHelper::fIFDH
ifdh_ns::ifdh * fIFDH
(optional) flux file handling
Definition: GENIEHelper.h:133
evgb::FillMCFlux
void FillMCFlux(genie::GFluxI *fdriver, simb::MCFlux &mcflux)
Definition: GENIE2ART.cxx:834
evgb::GENIEHelper::fGXMLPATH
std::string fGXMLPATH
locations for GENIE XML files
Definition: GENIEHelper.h:198
evgb::GENIEHelper::fFiducialCut
std::string fFiducialCut
configuration for geometry selector
Definition: GENIEHelper.h:205
simb::MCNeutrino::Nu
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
simb::MCFlux::fgenx
double fgenx
origin of ray from flux generator
Definition: MCFlux.h:100
evgb::GENIEHelper::Sample
bool Sample(simb::MCTruth &truth, simb::MCFlux &flux, simb::GTruth &gtruth)
Definition: GENIEHelper.cxx:1662
simb::MCFlux::fFluxType
simb::flux_code_ fFluxType
Definition: MCFlux.h:98
evgb::kNuTau
static const int kNuTau
Definition: GENIEHelper.cxx:210
evgb::GetRandomNumberSeed
unsigned int GetRandomNumberSeed()
Definition: evgenbase.h:22
evgb::GENIEHelper::fRandomTimeOffset
double fRandomTimeOffset
additional random time shift (ns) added to every particle time
Definition: GENIEHelper.h:181
evgb::GENIEHelper::fBeamCenter
TVector3 fBeamCenter
center of beam for histogram fluxes - must be in meters
Definition: GENIEHelper.h:175
evgb::GENIEHelper::fGlobalTimeOffset
double fGlobalTimeOffset
overall time shift (ns) added to every particle time
Definition: GENIEHelper.h:180
evgb::GENIEHelper::fSpillTimeConfig
std::string fSpillTimeConfig
alternative to flat spill distribution
Definition: GENIEHelper.h:182
evgb::GENIEHelper::fEventGeneratorList
std::string fEventGeneratorList
control over event topologies, was $GEVGL [Default]
Definition: GENIEHelper.h:197
evgb::GENIEHelper::fGeoManager
TGeoManager * fGeoManager
pointer to ROOT TGeoManager
Definition: GENIEHelper.h:123
evgb::kNue
static const int kNue
Definition: GENIEHelper.cxx:206
util::abs
constexpr auto abs(T v)
Returns the absolute value of the argument.
Definition: constexpr_math.h:40
evgb::GENIEHelper::TotalHistFlux
double TotalHistFlux()
Definition: GENIEHelper.cxx:591
std
STL namespace.
GENIE2ART.h
Functions for transforming GENIE objects into ART objects (and back)
evgb::GENIEHelper::fGENIEMsgThresholds
std::string fGENIEMsgThresholds
additional XML file setting Messager level thresholds (":" separated)
Definition: GENIEHelper.h:200
genie::flux::GPowerSpectrumAtmoFlux
A driver for a power spectrum atmospheric neutrino flux. Elements extensively reused from GAtmoFlux...
Definition: GPowerSpectrumAtmoFlux.h:28
evgb::GENIEHelper::fGenFlavors
std::vector< int > fGenFlavors
pdg codes for flavors to generate
Definition: GENIEHelper.h:186
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
mf::LogError
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Definition: MessageLogger.h:211
evgb::GENIEHelper::fTimeShifter
evgb::EvtTimeShiftI * fTimeShifter
generator for time offset within a spill
Definition: GENIEHelper.h:137
simb::MCFlux::fdk2gen
double fdk2gen
distance from decay to ray origin
Definition: MCFlux.h:103
tf
Definition: tf_graph.h:25
MCParticle.h
Particle class.
evgb::GENIEHelper::fDriver
genie::GMCJDriver * fDriver
Definition: GENIEHelper.h:130
evgb::GENIEHelper::fXSecMassPOT
double fXSecMassPOT
product of cross section, mass and POT/spill for histogram fluxes
Definition: GENIEHelper.h:172
evgb::GENIEHelper::ExpandFluxFilePatternsDirect
void ExpandFluxFilePatternsDirect()
Definition: GENIEHelper.cxx:1953
evgb::GENIEHelper::SqueezeFilePatterns
void SqueezeFilePatterns()
Definition: GENIEHelper.cxx:761
evgb::GENIEHelper::fSelectedFluxFiles
std::vector< std::string > fSelectedFluxFiles
flux files selected after wildcard expansion and subset selection
Definition: GENIEHelper.h:144
MCFlux.h
object containing MC flux information
evgb::GENIEHelper::fGenieEventRecord
genie::EventRecord * fGenieEventRecord
last generated event
Definition: GENIEHelper.h:126
evgb::GENIEHelper::fXSecTable
std::string fXSecTable
cross section file (was $GSPLOAD)
Definition: GENIEHelper.h:195
evgb::GENIEHelper::fFluxUpstreamZ
double fFluxUpstreamZ
z where flux starts from (if non-default, simple/ntuple only)
Definition: GENIEHelper.h:159
evgb::GENIEHelper::fFluxHistograms
std::vector< TH1D * > fFluxHistograms
histograms for each nu species
Definition: GENIEHelper.h:157
simb::kHistPlusFocus
Flux for positive horn focus.
Definition: MCFlux.h:18
genie::flux::GPowerSpectrumAtmoFlux::SetFlavors
void SetFlavors(std::vector< int > flavors)
Definition: GPowerSpectrumAtmoFlux.cxx:329
genie::flux::GPowerSpectrumAtmoFlux::AddFluxFile
void AddFluxFile(int neutrino_pdg, string filename)
Definition: GPowerSpectrumAtmoFlux.cxx:402
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
evgb::GENIEHelper::fSurroundingMass
double fSurroundingMass
Definition: GENIEHelper.h:178
scale
Double_t scale
Definition: plot.C:24
evgb::GENIEHelper::fGHepPrintLevel
int fGHepPrintLevel
GHepRecord::SetPrintLevel(), -1=no-print.
Definition: GENIEHelper.h:201
evgb::GENIEHelper::fFluxRotCfg
std::string fFluxRotCfg
how to interpret fFluxRotValues
Definition: GENIEHelper.h:150
evgb::GENIEHelper::fBeamDirection
TVector3 fBeamDirection
direction of the beam for histogram fluxes
Definition: GENIEHelper.h:174
evgb::GENIEHelper::fBeamName
std::string fBeamName
name of the beam we are simulating
Definition: GENIEHelper.h:149
evgb::GENIEHelper::fWorldVolume
std::string fWorldVolume
name of the world volume in the ROOT geometry
Definition: GENIEHelper.h:155
util::to_string
decltype(auto) constexpr to_string(T &&obj)
ADL-aware version of std::to_string.
evgb::GENIEHelper::fForceApplyFlxWgt
bool fForceApplyFlxWgt
apply GFluxI::Weight() before returning event
Definition: GENIEHelper.h:184
evgb::GENIEHelper::fTotalExposure
double fTotalExposure
pot used from flux ntuple
Definition: GENIEHelper.h:166
evgb::GENIEHelper::fHelperRandom
TRandom3 * fHelperRandom
random # generator for GENIEHelper
Definition: GENIEHelper.h:135
evgb::GENIEHelper::fTopVolume
std::string fTopVolume
top volume in the ROOT geometry in which to generate events
Definition: GENIEHelper.h:154
geo::vect::indices
constexpr std::array< std::size_t, geo::vect::dimension< Vector >)> indices()
Returns a sequence of indices valid for a vector of the specified type.
Definition: geo_vectors_utils.h:2113
evgb::GENIEHelper::fUseBlenderDist
bool fUseBlenderDist
get neutrino&#39;s travel distance from blender (default: true)
Definition: GENIEHelper.h:204
evgb::EvtTimeShiftI::PrintConfig
virtual void PrintConfig(bool verbose=true)=0
provide a means of printing the configuration
evgb::EvtTimeShiftI::IsRandomGeneratorSeeded
bool IsRandomGeneratorSeeded() const
Definition: EvtTimeShiftI.h:63
evgb::kNuMuBar
static const int kNuMuBar
Definition: GENIEHelper.cxx:209
evgb::GENIEHelper::fGeomD
genie::GeomAnalyzerI * fGeomD
Definition: GENIEHelper.h:127
evgb::GENIEHelper::fFluxSearchPaths
std::string fFluxSearchPaths
colon separated set of path stems
Definition: GENIEHelper.h:142
GENIEHelper.h
Wrapper for generating neutrino interactions with GENIE.
evgb::GENIEHelper::fGMSGLAYOUT
std::string fGMSGLAYOUT
format for GENIE log message [BASIC]|SIMPLE (SIMPLE=no timestamps)
Definition: GENIEHelper.h:199
bin
float bin[41]
Definition: plottest35.C:14
evgb::GENIEHelper::fBeamRadius
double fBeamRadius
radius of cylindar for histogram fluxes - must be in meters
Definition: GENIEHelper.h:176
evgb::EvtTimeShiftI::GetRandomGenerator
TRandom * GetRandomGenerator() const
Definition: EvtTimeShiftI.h:61
evgb::GENIEHelper::fMixerConfig
std::string fMixerConfig
configuration string for genie GFlavorMixerI
Definition: GENIEHelper.h:202
evgb::GENIEHelper::StringToBool
bool StringToBool(std::string v)
Definition: GENIEHelper.cxx:2542
evgb::GENIEHelper::fFluxD
genie::GFluxI * fFluxD
real flux driver
Definition: GENIEHelper.h:128
evgb::GENIEHelper::StartGENIEMessenger
void StartGENIEMessenger(std::string prodmode)
Definition: GENIEHelper.cxx:2383
simb::MCFlux::SetFluxGen
void SetFluxGen(double nue, double nuebar, double numu, double numubar, double nutau, double nutaubar)
Definition: MCFlux.cxx:214
evgb::GENIEHelper::fAddGenieVtxTime
bool fAddGenieVtxTime
incorporate time from flux window to interaction point and (possibily) proton-on-target to flux windo...
Definition: GENIEHelper.h:183
art::left
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
weight
double weight
Definition: plottest35.C:25
evgb::GENIEHelper::fMonoEnergy
double fMonoEnergy
energy of monoenergetic neutrinos
Definition: GENIEHelper.h:167
evgb::GENIEHelper::fDetectorMass
double fDetectorMass
mass of the detector in kg
Definition: GENIEHelper.h:177
evgb::GENIEHelper::fFluxRotValues
std::vector< double > fFluxRotValues
parameters for rotation
Definition: GENIEHelper.h:151
util::get
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
evgb::GENIEHelper::fAtmoSpectralIndex
double fAtmoSpectralIndex
atmo: Spectral index for power spectrum generation
Definition: GENIEHelper.h:192
evgb::GENIEHelper::InitializeFiducialSelection
void InitializeFiducialSelection()
Definition: GENIEHelper.cxx:930
evgb::SetEventGeneratorListAndTune
void SetEventGeneratorListAndTune(const std::string &evtlistname="", const std::string &tunename="${GENIE_XSEC_TUNE}")
Definition: GENIE2ART.cxx:128
evgb::GENIEHelper::fEnvironment
std::vector< std::string > fEnvironment
environmental variables and settings used by genie
Definition: GENIEHelper.h:194
evgenbase.h
mf::LogDebug
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Definition: MessageLogger.h:208
evgb::GENIEHelper::fMaxPathOutInfo
std::string fMaxPathOutInfo
output info if writing PathLengthList from GeomScan
Definition: GENIEHelper.h:207
evgb::GENIEHelper::fSpillExposure
double fSpillExposure
total exposure (i.e. pot) for this spill
Definition: GENIEHelper.h:165
EvtTimeShiftFactory.h
A class for generating concrete EvtTimeShiftI derived classes based on the factory pattern...
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:210
evgb::GENIEHelper::fAtmoEmin
double fAtmoEmin
atmo: Minimum energy of neutrinos in GeV
Definition: GENIEHelper.h:187
evgb::EvtTimeShiftFactory::Instance
static EvtTimeShiftFactory & Instance()
Definition: EvtTimeShiftFactory.cxx:35
evgb::GENIEHelper::fGeoFile
std::string fGeoFile
name of file containing the Geometry description
Definition: GENIEHelper.h:124
evgb::GENIEHelper::InitializeRockBoxSelection
void InitializeRockBoxSelection()
Definition: GENIEHelper.cxx:1085
evgb::GENIEHelper::fAtmoRl
double fAtmoRl
atmo: radius of the sphere on where the neutrinos are generated
Definition: GENIEHelper.h:189
genie::flux::GPowerSpectrumAtmoFlux::SetRadii
void SetRadii(double Rlongitudinal, double Rtransverse)
Definition: GPowerSpectrumAtmoFlux.cxx:316
evgb::GENIEHelper::ExpandFluxFilePatternsIFDH
void ExpandFluxFilePatternsIFDH()
Definition: GENIEHelper.cxx:2109
evgb::GENIEHelper::fDebugFlags
unsigned int fDebugFlags
set bits to enable debug info
Definition: GENIEHelper.h:208
evgb::EvtTimeShiftI::TimeOffset
virtual double TimeOffset()=0
evgb
Physics generators for neutrinos, cosmic rays, and others.
Definition: CRYHelper.cxx:33
evgb::GENIEHelper::fFluxCopyMethod
std::string fFluxCopyMethod
"DIRECT" = old direct access method, otherwise = ifdh approach schema ("" okay)
Definition: GENIEHelper.h:147
ff
TFile ff[ntarg]
Definition: Style.C:26
evgb::FillGTruth
void FillGTruth(const genie::EventRecord *grec, simb::GTruth &gtruth)
Definition: GENIE2ART.cxx:391
evgb::GENIEHelper::fFluxCleanup
std::string fFluxCleanup
"ALWAYS", "/var/tmp", "NEVER"
Definition: GENIEHelper.h:148
simb::MCTruth
Event generator information.
Definition: MCTruth.h:32
evgb::GENIEHelper::fHistEventsPerSpill
double fHistEventsPerSpill
number of events per spill for histogram fluxes - changes each spill
Definition: GENIEHelper.h:163
evgb::GENIEHelper::fFluxRotation
TRotation * fFluxRotation
rotation for atmos / astro flux coord systems
Definition: GENIEHelper.h:152
e
Float_t e
Definition: plot.C:35
evgb::GENIEHelper::fGeomScan
std::string fGeomScan
configuration for geometry scan to determine max pathlengths
Definition: GENIEHelper.h:206
evgb::GENIEHelper::fSpillEvents
int fSpillEvents
total events for this spill
Definition: GENIEHelper.h:164
w
Float_t w
Definition: plot.C:20
evgb::GENIEHelper::fFluxD2GMCJD
genie::GFluxI * fFluxD2GMCJD
flux driver passed to genie GMCJDriver, might be GFluxBlender
Definition: GENIEHelper.h:129
util::quantities::centimeter
centimeter_as<> centimeter
Type of space stored in centimeters, in double precision.
Definition: spacetime.h:441
evgb::kNueBar
static const int kNueBar
Definition: GENIEHelper.cxx:207
evgb::GENIEHelper::fFunctionalFlux
std::string fFunctionalFlux
Definition: GENIEHelper.h:168
evgb::GENIEHelper::fMaxFluxFileMB
int fMaxFluxFileMB
maximum size of flux files (MB)
Definition: GENIEHelper.h:145
evgb::kNuMu
static const int kNuMu
Definition: GENIEHelper.cxx:208
evgb::GENIEHelper::fPOTPerSpill
double fPOTPerSpill
number of pot per spill
Definition: GENIEHelper.h:162
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
evgb::GENIEHelper::fFluxType
std::string fFluxType
Definition: GENIEHelper.h:139
evgb::GENIEHelper::GENIEHelper
GENIEHelper(fhicl::ParameterSet const &pset, TGeoManager *rootGeom, std::string const &rootFile, double const &detectorMass)
Definition: GENIEHelper.cxx:214
evgb::EvtTimeShiftFactory::GetEvtTimeShift
evgb::EvtTimeShiftI * GetEvtTimeShift(const std::string &name, const std::string &config="") const
Definition: EvtTimeShiftFactory.cxx:50
vertex
vertex reconstruction
Definition: AggregateVertexAna_module.cc:32