LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
TestGENIEHelper_module.cc
Go to the documentation of this file.
1 
7 #include <cassert>
8 #include <cstdlib>
9 #include <string>
10 #include <sstream>
11 #include <vector>
12 #include <map>
13 #include <unistd.h>
14 #include <string>
15 
16 // ROOT includes
17 #include "TStopwatch.h"
18 #include "TGeoManager.h"
19 
20 // Framework includes
25 #include "fhiclcpp/ParameterSet.h"
29 
30 
33 
35 
41 
42 //--- Dk2Nu additions
43 //--- BEGIN
44 #include "dk2nu/tree/dk2nu.h"
45 #include "dk2nu/tree/NuChoice.h"
46 #include "dk2nu/genie/GDk2NuFlux.h"
47 //--- END
48 
49 // GENIE includes
50 #ifdef GENIE_PRE_R3
51  #include "GENIE/EVGCore/EventRecord.h"
52 #else
53  #include "GENIE/Framework/EventGen/EventRecord.h"
54 #endif
55 
56 //#undef PUT_DK2NU_ASSN
57 #define PUT_DK2NU_ASSN 1
58 
60 namespace evgen {
61 
64 
65  public:
66 
67  explicit TestGENIEHelper(fhicl::ParameterSet const &pset);
68  virtual ~TestGENIEHelper();
69 
70  void produce(art::Event& evt);
71  void beginJob();
72  void beginRun(art::Run &run);
73  void endSubRun(art::SubRun &sr);
74 
75  private:
76 
78  TStopwatch fStopwatch;
80  unsigned int fDebugFlags;
81 
82  };
83 }
84 
85 namespace evgen {
86 
87  //___________________________________________________________________________
89  : EDProducer{pset}
90  , fGENIEHelp (0)
91  //, fPassEmptySpills (pset.get< bool >("PassEmptySpills"))
92  //, fSpillCounter (0)
93  //, fPOTPerSpill (pset.get< double >("POTPerSpill", 5.0e13))
94  , fEventsPerSpill (pset.get< double >("EventsPerSpill", 1))
95  //, fTotalExposure (0)
96  //, fTotalPOTLimit (pset.get< double >("TotalPOTLimit"))
97  , fDebugFlags (pset.get< unsigned int >("DebugFlags", 0))
98  {
99  fStopwatch.Start();
100 
101  produces< std::vector<simb::MCTruth> >();
102  produces< std::vector<simb::MCFlux> >();
103  produces< std::vector<simb::GTruth> >();
104  //produces< sumdata::SpillData >();
105  //produces< sumdata::POTSum, art::InSubRun >();
106  //produces< sumdata::RunData, art::InRun >();
107  // Associate every truth with the flux it came from
108  produces< art::Assns<simb::MCTruth, simb::MCFlux> >();
109  produces< art::Assns<simb::MCTruth, simb::GTruth> >();
110 
111  //--- Dk2Nu additions
112  //--- BEGIN
113  produces< std::vector<bsim::Dk2Nu> >();
114  produces< std::vector<bsim::NuChoice> >();
115 #ifdef PUT_DK2NU_ASSN
116  produces< art::Assns<simb::MCTruth, bsim::Dk2Nu> >();
117  produces< art::Assns<simb::MCTruth, bsim::NuChoice> >();
118 #endif
119  //--- END
120 
121  //art::ServiceHandle<geo::Geometry> geo;
122  std::string geomFileName =
123  pset.get<std::string>("GeomFileName");
124  mf::LogInfo("TestGENIEHelper") << "using GeomFileName '"
125  << geomFileName << "'";
126  TGeoManager::Import(geomFileName.c_str());
127  double detectorMass = 1; // no generally used (except for histogram)
128 
129  fGENIEHelp = new evgb::GENIEHelper(pset,
130  gGeoManager,
131  geomFileName,
132  detectorMass);
133  /*
134  geo->ROOTGeoManager(),
135  geo->ROOTFile(),
136  geo->TotalMass(pset.get< std::string>("TopVolume").c_str()));
137  */
138  }
139 
140  //___________________________________________________________________________
142  {
143  fStopwatch.Stop();
144  mf::LogInfo("TestGENIEHelper") << "real time to produce file: "
145  << fStopwatch.RealTime();
146  delete fGENIEHelp; // clean up, and let dtor do its thing
147  }
148 
149  //___________________________________________________________________________
151  {
152  }
153 
154  //___________________________________________________________________________
156  {
157  // grab the geometry object to see what geometry we are using
158  // art::ServiceHandle<geo::Geometry> geo;
159 
160  /*
161  std::unique_ptr<sumdata::RunData>
162  runcol(new sumdata::RunData(geo->DetId(),
163  geo->FileBaseName(),
164  geo->ExtractGDML()));
165 
166  run.put(std::move(runcol));
167  */
168 
169  std::cerr << " *** TestGENIEHelper::beginRun() begin "
170  << std::endl << std::flush;
171 
172 
173  // initialize the GENIEHelper here rather than in beginJob to
174  // avoid problems with the Geometry reloading at a run boundary.
175  // If we ever make more than one run in a single job we will have
176  // to re-evaluate
178  //fTotalExposure = 0.0;
179 
180  std::cerr << " *** TestGENIEHelper::beginRun() done "
181  << std::endl << std::flush;
182 
183  return;
184  }
185 
186  //___________________________________________________________________________
188  {
189  /*
190  std::unique_ptr< sumdata::POTSum > p(new sumdata::POTSum);
191 
192  // p->totpot = fGENIEHelp->TotalExposure();
193  // p->totgoodpot = fGENIEHelp->TotalExposure();
194  p->totpot = fTotalExposure;
195  p->totgoodpot = fTotalExposure;
196  p->totspills = fSpillCounter;
197  p->goodspills = fSpillCounter;
198  p->Print(std::cout);
199 
200  sr.put(std::move(p));
201  */
202 
203  mf::LogInfo("TestGENIEHelper") << "Total Exposure was "
205 
206  }
207 
208  //___________________________________________________________________________
210  {
211  // A temporary value is needed to store the spill exposure that GENIEHelper uses. TestGENIEHelper
212  // needs to remember the number after GENIEHelper has reset it to zero for the purposes of
213  // updating fTotalExposure.
214  //double SpillExpTemp = 0.0;
215 
216  std::unique_ptr< std::vector<simb::MCTruth> > truthcol(new std::vector<simb::MCTruth>);
217  std::unique_ptr< std::vector<simb::MCFlux> > fluxcol (new std::vector<simb::MCFlux >);
218  std::unique_ptr< std::vector<simb::GTruth> > gtruthcol (new std::vector<simb::GTruth >);
219  std::unique_ptr< art::Assns<simb::MCTruth, simb::GTruth> > tgtassn(new art::Assns<simb::MCTruth, simb::GTruth>);
220  std::unique_ptr< art::Assns<simb::MCTruth, simb::MCFlux> > assns(new art::Assns<simb::MCTruth, simb::MCFlux>);
221 
222  std::cerr << " ******************************* TestGENIEHelper::produce() " << std::endl << std::flush;
223  std::cout << " stopwatch at produce() ";
224  fStopwatch.Print("um"); fStopwatch.Continue();
225  std::cout << std::flush;
226 
227  //--- Dk2Nu additions
228  //--- BEGIN
229  std::unique_ptr< std::vector<bsim::Dk2Nu> >
230  dk2nucol(new std::vector<bsim::Dk2Nu>);
231  std::unique_ptr< std::vector<bsim::NuChoice> >
232  nuchoicecol(new std::vector<bsim::NuChoice>);
233 
234  std::unique_ptr< art::Assns<simb::MCTruth, bsim::Dk2Nu> >
236  std::unique_ptr< art::Assns<simb::MCTruth, bsim::NuChoice> >
238  //--- END
239 
240  while ( ! fGENIEHelp->Stop() ) {
241 
242  simb::MCTruth truth;
243  simb::MCFlux flux;
244  simb::GTruth gTruth;
245 
246  std::cerr << " *** TestGENIEHelper::produce() about to sample "
247  << truthcol->size()
248  << std::endl << std::flush;
249  std::cout << " stopwatch before Sample() ";
250  fStopwatch.Print("um"); fStopwatch.Continue();
251  std::cout << std::flush;
252 
253 
254  // GENIEHelper returns a false in the sample method if
255  // either no neutrino was generated, or the interaction
256  // occurred beyond the detector's z extent - ie something we
257  // would never see anyway.
258  if ( fGENIEHelp->Sample(truth, flux, gTruth ) ) {
259 
260  std::cout << " stopwatch after Sample() ";
261  fStopwatch.Print("um"); fStopwatch.Continue();
262  std::cout << std::flush;
263 
264  truthcol ->push_back(truth);
265  gtruthcol->push_back(gTruth);
266  fluxcol ->push_back(flux);
267 
268  evgb::util::CreateAssn(*this, evt, *truthcol, *fluxcol, *assns,
269  fluxcol->size()-1, fluxcol->size());
270 
271  evgb::util::CreateAssn(*this, evt, *truthcol, *gtruthcol, *tgtassn,
272  gtruthcol->size()-1, gtruthcol->size());
273 
274  //--- Dk2Nu additions
275  //--- BEGIN
276  genie::GFluxI* fdriver = fGENIEHelp->GetFluxDriver(true);
277  genie::flux::GDk2NuFlux* dk2nuDriver =
278  dynamic_cast<genie::flux::GDk2NuFlux*>(fdriver);
279  if ( dk2nuDriver ) {
280  const bsim::Dk2Nu& dk2nuObj = dk2nuDriver->GetDk2Nu();
281  dk2nucol ->push_back(dk2nuObj);
282  const bsim::NuChoice& nuchoiceObj = dk2nuDriver->GetNuChoice();
283  nuchoicecol->push_back(nuchoiceObj);
284 
285  if ( (fDebugFlags & 0x10 ) != 0 ) {
286  std::cout << "---------==== creation dump" << std::endl;
287  genie::EventRecord* gevtrec = fGENIEHelp->GetGenieEventRecord();
288  if ( gevtrec ) std::cout << *gevtrec << std::endl;
289  std::cout << dk2nuObj << std::endl;
290  const bsim::Decay& decay = dk2nuObj.decay;
291  std::cout << " necm " << decay.necm
292  << " muparp4 " << decay.muparpx << " "
293  << decay.muparpy << " " << decay.muparpz << " "
294  << decay.mupare << std::endl;
295  std::cout << nuchoiceObj << std::endl;
296  }
297 
298 #ifdef PUT_DK2NU_ASSN
299  evgb::util::CreateAssn(*this, evt, *truthcol, *dk2nucol, *dk2nuassn,
300  dk2nucol->size()-1, dk2nucol->size());
301  evgb::util::CreateAssn(*this, evt, *truthcol, *nuchoicecol, *nuchoiceassn,
302  nuchoicecol->size()-1, nuchoicecol->size());
303 #endif
304  }
305  //--- END
306 
307  std::cerr << " *** TestGENIEHelper::produce() sample success size "
308  << truthcol->size()
309  << std::endl << std::flush;
310  //RWH//if ( dk2nuDriver ) { std::cout << dk2nuDriver->GetDk2Nu() << std::endl; }
311  std::cout << " stopwatch after push_back + CreateAssn ";
312  fStopwatch.Print("um"); fStopwatch.Continue();
313  std::cout << std::flush;
314 
315 
316  } // end if genie was able to make an event
317 
318  } // end event generation loop
319 
320  // put the collections in the event
321  evt.put(std::move(truthcol));
322  evt.put(std::move(fluxcol));
323  evt.put(std::move(gtruthcol));
324  evt.put(std::move(assns));
325  evt.put(std::move(tgtassn));
326 
327  std::cerr << " *** TestGENIEHelper::produce() done "
328  << " event " << evt.event()
329  << std::endl << std::flush;
330 
331  //--- Dk2Nu additions
332  //--- BEGIN
333  // in the constructor we said these were produced ...
334  // ... so we have to put them in the record, even if just empty
335  evt.put(std::move(dk2nucol));
336  evt.put(std::move(nuchoicecol));
337 
338 #ifdef PUT_DK2NU_ASSN
339  std::cerr << " *** TestGENIEHelper::produce()"
340  << " put dk2nuAssn + nuchoiceAssn ** "
341  << " event " << evt.event()
342  << std::endl << std::flush;
343 
344  evt.put(std::move(dk2nuassn));
345  evt.put(std::move(nuchoiceassn));
346 
347  std::cerr << " *** TestGENIEHelper::produce() finished put "
348  << " event " << evt.event()
349  << std::endl << std::flush;
350 
351 #endif
352  //--- END
353 
354  } // end of produce method
355 }// namespace
356 
genie::GFluxI * GetFluxDriver(bool base=true)
Definition: GENIEHelper.h:91
void produce(art::Event &evt)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
genie::EventRecord * GetGenieEventRecord()
Definition: GENIEHelper.h:84
bool Sample(simb::MCTruth &truth, simb::MCFlux &flux, simb::GTruth &gtruth)
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
Particle class.
object containing MC flux information
Definition: Run.h:37
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
A module to check the results from the Monte Carlo generator.
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
int fEventsPerSpill
negative for Poisson()
EventNumber_t event() const
Definition: Event.cc:41
Wrapper for generating neutrino interactions with GENIE.
void endSubRun(art::SubRun &sr)
TestGENIEHelper(fhicl::ParameterSet const &pset)
double TotalExposure() const
Definition: GENIEHelper.h:71
evgb::GENIEHelper * fGENIEHelp
GENIEHelper object.
TCEvent evt
Definition: DataStructs.cxx:8
Event generator information.
Definition: MCTruth.h:32
Event Generation using GENIE, cosmics or single particles.