LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
GENIEGen_module.cc
Go to the documentation of this file.
1 //
3 //
4 // GENIE neutrino event generator
5 //
6 // brebel@fnal.gov
7 //
9 
10 #include <cstdlib>
11 #include <iomanip>
12 #include <iostream>
13 #include <memory>
14 #include <string>
15 
16 // ROOT includes
17 #include "TDatabasePDG.h"
18 #include "TH1.h"
19 #include "TH2.h"
20 #include "TStopwatch.h"
21 
22 // Framework includes
29 #include "art_root_io/TFileService.h"
31 #include "fhiclcpp/ParameterSet.h"
33 
34 // art extensions
36 
37 // LArSoft includes
42 #include "lardataalg/MCDumpers/MCDumpers.h" // sim::dump namespace
48 
49 // dk2nu extensions
50 #include "dk2nu/genie/GDk2NuFlux.h"
51 #include "dk2nu/tree/NuChoice.h"
52 #include "dk2nu/tree/dk2nu.h"
53 
54 #include "GENIE/Framework/EventGen/EventRecord.h"
57 
59 namespace evgen {
91  class GENIEGen : public art::EDProducer {
92  public:
93  explicit GENIEGen(fhicl::ParameterSet const& pset);
94  virtual ~GENIEGen();
95 
96  void produce(art::Event& evt);
97  void beginJob();
98  void beginRun(art::Run& run);
99  void beginSubRun(art::SubRun& sr);
100  void endSubRun(art::SubRun& sr);
101 
102  private:
103  std::string ParticleStatus(int StatusCode);
104  std::string ReactionChannel(int ccnc, int mode);
105 
106  void FillHistograms(simb::MCTruth mc);
107 
109  bool
111  std::vector<double> fVtxPosHistRange;
112 
114  TStopwatch fStopwatch;
115 
119 
120  double fPrevTotPOT;
122 
123  TH1F* fGenerated[6];
124 
125  TH1F* fVertexX;
126  TH1F* fVertexY;
127  TH1F* fVertexZ;
128 
129  TH2F* fVertexXY;
130  TH2F* fVertexXZ;
131  TH2F* fVertexYZ;
132 
133  TH1F* fDCosX;
134  TH1F* fDCosY;
135  TH1F* fDCosZ;
136 
137  TH1F* fMuMomentum;
138  TH1F* fMuDCosX;
139  TH1F* fMuDCosY;
140  TH1F* fMuDCosZ;
141 
142  TH1F* fEMomentum;
143  TH1F* fEDCosX;
144  TH1F* fEDCosY;
145  TH1F* fEDCosZ;
146 
147  TH1F* fCCMode;
148  TH1F* fNCMode;
149 
150  TH1F* fDeltaE;
151  TH1F* fECons;
152  };
153 }
154 
155 namespace evgen {
156 
157  //____________________________________________________________________________
159  : EDProducer{pset}
160  , fGENIEHelp(0)
161  , fDefinedVtxHistRange(pset.get<bool>("DefinedVtxHistRange"))
162  , fVtxPosHistRange(pset.get<std::vector<double>>("VtxPosHistRange"))
163  , fPassEmptySpills(pset.get<bool>("PassEmptySpills"))
164  , fGlobalTimeOffset(pset.get<double>("GlobalTimeOffset", 0))
165  , fRandomTimeOffset(pset.get<double>("RandomTimeOffset", 1600.)) // BNB default value
166  , fBeamType(::sim::kBNB)
167  {
168  fStopwatch.Start();
169 
170  produces<std::vector<simb::MCTruth>>();
171  produces<std::vector<simb::MCFlux>>();
172  produces<std::vector<simb::GTruth>>();
173  produces<sumdata::RunData, art::InRun>();
174  produces<sumdata::POTSummary, art::InSubRun>();
175  produces<art::Assns<simb::MCTruth, simb::MCFlux>>();
176  produces<art::Assns<simb::MCTruth, simb::GTruth>>();
177  produces<std::vector<sim::BeamGateInfo>>();
178 
179  // dk2nu additions
180  if (pset.get<std::string>("FluxType").find("dk2nu") != std::string::npos) {
181  produces<std::vector<bsim::Dk2Nu>>();
182  produces<std::vector<bsim::NuChoice>>();
183  produces<art::Assns<simb::MCTruth, bsim::Dk2Nu>>();
184  produces<art::Assns<simb::MCTruth, bsim::NuChoice>>();
185  }
186 
187  std::string beam_type_name = pset.get<std::string>("BeamName");
188 
189  if (beam_type_name == "numi")
190 
192 
193  else if (beam_type_name == "booster")
194 
196 
197  else
198 
200 
202 
203  signed int temp_seed; // the seed read by GENIEHelper is a signed integer...
204  fhicl::ParameterSet GENIEconfig(pset);
205  if (!GENIEconfig.get_if_present("RandomSeed",
206  temp_seed)) { // TODO use has_key() when it becomes available
207  // no RandomSeed specified; check for the LArSoft-style "Seed" instead:
208  // obtain the random seed from a service,
209  // unless overridden in configuration with key "Seed"
210  unsigned int seed;
211  if (!GENIEconfig.get_if_present("Seed", seed))
212  seed = art::ServiceHandle<rndm::NuRandomService>()->getSeed();
213 
214  // The seed is not passed to RandomNumberGenerator,
215  // since GENIE uses a TRandom generator that is owned by the GENIEHelper.
216  // Instead, we explicitly configure the random seed for GENIEHelper:
217  GENIEconfig.put("RandomSeed", seed);
218  } // if no RandomSeed present
219 
220  double detectorMass = 0;
221  // detectorMass is _only_ needed by GENIEHelper in the case of
222  // histogram flux and non-fixed # of event/spill (ie. POTPerSpill non-zero)
223  if (pset.get<std::string>("FluxType").find("histogram") == 0 &&
224  pset.get<double>("EventsPerSpill") == 0.0 && pset.get<double>("POTPerSpill") > 0.0) {
225  TStopwatch timer;
226  timer.Start();
227  detectorMass = geo->TotalMass(pset.get<std::string>("TopVolume").c_str());
228  timer.Stop();
229  mf::LogInfo("GENIEProductionTime")
230  << "real time to calculate TotalMass: " << timer.RealTime() << " sec";
231  }
232 
233  fGENIEHelp =
234  new evgb::GENIEHelper(GENIEconfig, geo->ROOTGeoManager(), geo->ROOTFile(), detectorMass);
235  }
236 
237  //____________________________________________________________________________
239  {
240  if (fGENIEHelp) delete fGENIEHelp;
241  fStopwatch.Stop();
242  mf::LogInfo("GENIEProductionTime") << "real time to produce file: " << fStopwatch.RealTime();
243  }
244 
245  //____________________________________________________________________________
247  {
249 
250  fPrevTotPOT = 0.;
251  fPrevTotGoodPOT = 0.;
252 
253  // Get access to the TFile service.
255 
256  fGenerated[0] = tfs->make<TH1F>("fGenerated_necc", "", 100, 0.0, 20.0);
257  fGenerated[1] = tfs->make<TH1F>("fGenerated_nebcc", "", 100, 0.0, 20.0);
258  fGenerated[2] = tfs->make<TH1F>("fGenerated_nmcc", "", 100, 0.0, 20.0);
259  fGenerated[3] = tfs->make<TH1F>("fGenerated_nmbcc", "", 100, 0.0, 20.0);
260  fGenerated[4] = tfs->make<TH1F>("fGenerated_nnc", "", 100, 0.0, 20.0);
261  fGenerated[5] = tfs->make<TH1F>("fGenerated_nbnc", "", 100, 0.0, 20.0);
262 
263  fDCosX = tfs->make<TH1F>("fDCosX", ";dx/ds", 200, -1., 1.);
264  fDCosY = tfs->make<TH1F>("fDCosY", ";dy/ds", 200, -1., 1.);
265  fDCosZ = tfs->make<TH1F>("fDCosZ", ";dz/ds", 200, -1., 1.);
266 
267  fMuMomentum = tfs->make<TH1F>("fMuMomentum", ";p_{#mu} (GeV/c)", 500, 0., 50.);
268  fMuDCosX = tfs->make<TH1F>("fMuDCosX", ";dx/ds;", 200, -1., 1.);
269  fMuDCosY = tfs->make<TH1F>("fMuDCosY", ";dy/ds;", 200, -1., 1.);
270  fMuDCosZ = tfs->make<TH1F>("fMuDCosZ", ";dz/ds;", 200, -1., 1.);
271 
272  fEMomentum = tfs->make<TH1F>("fEMomentum", ";p_{e} (GeV/c)", 500, 0., 50.);
273  fEDCosX = tfs->make<TH1F>("fEDCosX", ";dx/ds;", 200, -1., 1.);
274  fEDCosY = tfs->make<TH1F>("fEDCosY", ";dy/ds;", 200, -1., 1.);
275  fEDCosZ = tfs->make<TH1F>("fEDCosZ", ";dz/ds;", 200, -1., 1.);
276 
277  fCCMode = tfs->make<TH1F>("fCCMode", ";CC Interaction Mode;", 4, 0., 4.);
278  fCCMode->GetXaxis()->SetBinLabel(1, "QE");
279  fCCMode->GetXaxis()->SetBinLabel(2, "Res");
280  fCCMode->GetXaxis()->SetBinLabel(3, "DIS");
281  fCCMode->GetXaxis()->SetBinLabel(4, "Coh");
282  fCCMode->GetXaxis()->CenterLabels();
283 
284  fNCMode = tfs->make<TH1F>("fNCMode", ";NC Interaction Mode;", 4, 0., 4.);
285  fNCMode->GetXaxis()->SetBinLabel(1, "QE");
286  fNCMode->GetXaxis()->SetBinLabel(2, "Res");
287  fNCMode->GetXaxis()->SetBinLabel(3, "DIS");
288  fNCMode->GetXaxis()->SetBinLabel(4, "Coh");
289  fNCMode->GetXaxis()->CenterLabels();
290 
291  fDeltaE = tfs->make<TH1F>("fDeltaE", ";#Delta E_{#nu} (GeV);", 200, -1., 1.);
292  fECons = tfs->make<TH1F>("fECons", ";#Delta E(#nu,lepton);", 500, -5., 5.);
293 
294  if (fDefinedVtxHistRange == false) {
296  double x = 2.1 * geo->DetHalfWidth();
297  double y = 2.1 * geo->DetHalfHeight();
298  double z = 2. * geo->DetLength();
299  int xdiv = TMath::Nint(2 * x / 5.);
300  int ydiv = TMath::Nint(2 * y / 5.);
301  int zdiv = TMath::Nint(2 * z / 5.);
302 
303  fVertexX = tfs->make<TH1F>("fVertexX", ";x (cm)", xdiv, -0.1 * x, x);
304  fVertexY = tfs->make<TH1F>("fVertexY", ";y (cm)", ydiv, -y, y);
305  fVertexZ = tfs->make<TH1F>("fVertexZ", ";z (cm)", zdiv, -0.1 * z, z);
306 
307  fVertexXY = tfs->make<TH2F>("fVertexXY", ";x (cm);y (cm)", xdiv, -0.1 * x, x, ydiv, -y, y);
308  fVertexXZ =
309  tfs->make<TH2F>("fVertexXZ", ";z (cm);x (cm)", zdiv, -0.2 * z, z, xdiv, -0.1 * x, x);
310  fVertexYZ = tfs->make<TH2F>("fVertexYZ", ";z (cm);y (cm)", zdiv, -0.2 * z, z, ydiv, -y, y);
311  }
312  else {
313  int xdiv = TMath::Nint((fVtxPosHistRange[1] - fVtxPosHistRange[0]) / 5.);
314  int ydiv = TMath::Nint((fVtxPosHistRange[3] - fVtxPosHistRange[2]) / 5.);
315  int zdiv = TMath::Nint((fVtxPosHistRange[5] - fVtxPosHistRange[4]) / 5.);
316 
317  fVertexX =
318  tfs->make<TH1F>("fVertexX", ";x (cm)", xdiv, fVtxPosHistRange[0], fVtxPosHistRange[1]);
319  fVertexY =
320  tfs->make<TH1F>("fVertexY", ";y (cm)", ydiv, fVtxPosHistRange[2], fVtxPosHistRange[3]);
321  fVertexZ =
322  tfs->make<TH1F>("fVertexZ", ";z (cm)", zdiv, fVtxPosHistRange[4], fVtxPosHistRange[5]);
323 
324  fVertexXY = tfs->make<TH2F>("fVertexXY",
325  ";x (cm);y (cm)",
326  xdiv,
327  fVtxPosHistRange[0],
328  fVtxPosHistRange[1],
329  ydiv,
330  fVtxPosHistRange[2],
331  fVtxPosHistRange[3]);
332  fVertexXZ = tfs->make<TH2F>("fVertexXZ",
333  ";z (cm);x (cm)",
334  zdiv,
335  fVtxPosHistRange[4],
336  fVtxPosHistRange[5],
337  xdiv,
338  fVtxPosHistRange[0],
339  fVtxPosHistRange[1]);
340  fVertexYZ = tfs->make<TH2F>("fVertexYZ",
341  ";z (cm);y (cm)",
342  zdiv,
343  fVtxPosHistRange[4],
344  fVtxPosHistRange[5],
345  ydiv,
346  fVtxPosHistRange[2],
347  fVtxPosHistRange[3]);
348  }
349  }
350 
351  //____________________________________________________________________________
353  {
355  run.put(std::make_unique<sumdata::RunData>(geo->DetectorName()), art::fullRun());
356  }
357 
358  //____________________________________________________________________________
360  {
361 
364 
365  return;
366  }
367 
368  //____________________________________________________________________________
370  {
371 
372  auto p = std::make_unique<sumdata::POTSummary>();
373 
374  p->totpot = fGENIEHelp->TotalExposure() - fPrevTotPOT;
375  p->totgoodpot = fGENIEHelp->TotalExposure() - fPrevTotGoodPOT;
376 
377  sr.put(std::move(p), art::subRunFragment());
378 
379  return;
380  }
381 
382  //____________________________________________________________________________
384  {
385  std::unique_ptr<std::vector<simb::MCTruth>> truthcol(new std::vector<simb::MCTruth>);
386  std::unique_ptr<std::vector<simb::MCFlux>> fluxcol(new std::vector<simb::MCFlux>);
387  std::unique_ptr<std::vector<simb::GTruth>> gtruthcol(new std::vector<simb::GTruth>);
388  std::unique_ptr<art::Assns<simb::MCTruth, simb::MCFlux>> tfassn(
390  std::unique_ptr<art::Assns<simb::MCTruth, simb::GTruth>> tgtassn(
392  std::unique_ptr<std::vector<sim::BeamGateInfo>> gateCollection(
393  new std::vector<sim::BeamGateInfo>);
394 
395  std::unique_ptr<std::vector<bsim::Dk2Nu>> dk2nucol(new std::vector<bsim::Dk2Nu>);
396  std::unique_ptr<std::vector<bsim::NuChoice>> nuchoicecol(new std::vector<bsim::NuChoice>);
397  std::unique_ptr<art::Assns<simb::MCTruth, bsim::Dk2Nu>> dk2nuassn(
399  std::unique_ptr<art::Assns<simb::MCTruth, bsim::NuChoice>> nuchoiceassn(
401 
402  genie::flux::GDk2NuFlux* dk2nuDriver =
403  dynamic_cast<genie::flux::GDk2NuFlux*>(fGENIEHelp->GetFluxDriver(true));
404  while (truthcol->size() < 1) {
405  while (!fGENIEHelp->Stop()) {
406 
407  simb::MCTruth truth;
408  simb::MCFlux flux;
409  simb::GTruth gTruth;
410 
411  // GENIEHelper returns a false in the sample method if
412  // either no neutrino was generated, or the interaction
413  // occurred beyond the detector's z extent - ie something we
414  // would never see anyway.
415  if (fGENIEHelp->Sample(truth, flux, gTruth)) {
416 
417  truthcol->push_back(truth);
418  fluxcol->push_back(flux);
419  gtruthcol->push_back(gTruth);
420  auto const truthPtr = art::PtrMaker<simb::MCTruth>{evt}(truthcol->size() - 1);
421  tfassn->addSingle(truthPtr, art::PtrMaker<simb::MCFlux>{evt}(fluxcol->size() - 1));
422  tgtassn->addSingle(truthPtr, art::PtrMaker<simb::GTruth>{evt}(gtruthcol->size() - 1));
423 
424  FillHistograms(truth);
425 
426  if (dk2nuDriver) {
427  const bsim::Dk2Nu& dk2nuObj = dk2nuDriver->GetDk2Nu();
428  dk2nucol->push_back(dk2nuObj);
429  const bsim::NuChoice& nuchoiceObj = dk2nuDriver->GetNuChoice();
430  nuchoicecol->push_back(nuchoiceObj);
432  evt, *truthcol, *dk2nucol, *dk2nuassn, dk2nucol->size() - 1, dk2nucol->size());
433  util::CreateAssn(evt,
434  *truthcol,
435  *nuchoicecol,
436  *nuchoiceassn,
437  nuchoicecol->size() - 1,
438  nuchoicecol->size());
439  }
440 
441  // check that the process code is not unsupported by GENIE
442  // (see issue #18025 for reference);
443  // if it is, print all the information we can about this truth record
444  if (truth.NeutrinoSet() &&
446  mf::LogWarning log("GENIEmissingProcessMapping");
447  log << "Found an interaction that is not represented by the interaction type code in "
448  "GENIE:"
449  "\nMCTruth record:"
450  "\n";
451  sim::dump::DumpMCTruth(log, truth, 2U); // 2 trajectory points per line
452  log << "\nGENIE truth record:"
453  "\n";
454  sim::dump::DumpGTruth(log, gTruth);
455  } // if
456 
457  } // end if genie was able to make an event
458 
459  } // end event generation loop
460 
461  // check to see if we are to pass empty spills
462  if (truthcol->size() < 1 && fPassEmptySpills) {
463  MF_LOG_DEBUG("GENIEGen") << "no events made for this spill but "
464  << "passing it on and ending the event anyway";
465  break;
466  }
467 
468  } // end loop while no interactions are made
469 
470  // Create a simulated "beam gate" for these neutrino events.
471  // We're creating a vector of these because, in a
472  // distant-but-possible future, we may be generating more than one
473  // beam gate within a simulated time window.
474  gateCollection->push_back(sim::BeamGateInfo(fGlobalTimeOffset, fRandomTimeOffset, fBeamType));
475 
476  // put the collections in the event
477  evt.put(std::move(truthcol));
478  evt.put(std::move(fluxcol));
479  evt.put(std::move(gtruthcol));
480  evt.put(std::move(tfassn));
481  evt.put(std::move(tgtassn));
482  evt.put(std::move(gateCollection));
483 
484  // dk2nu additions
485  if (dk2nuDriver) {
486  evt.put(std::move(dk2nucol));
487  evt.put(std::move(nuchoicecol));
488  evt.put(std::move(dk2nuassn));
489  evt.put(std::move(nuchoiceassn));
490  }
491 
492  return;
493  }
494 
495  //......................................................................
496  std::string GENIEGen::ParticleStatus(int StatusCode)
497  {
498  int code = StatusCode;
499  std::string ParticleStatusName;
500 
501  switch (code) {
502  case -1: ParticleStatusName = "kIStUndefined"; break;
503  case 0: ParticleStatusName = "kIStInitialState"; break;
504  case 1: ParticleStatusName = "kIStStableFinalState"; break;
505  case 2: ParticleStatusName = "kIStIntermediateState"; break;
506  case 3: ParticleStatusName = "kIStDecayedState"; break;
507  case 11: ParticleStatusName = "kIStNucleonTarget"; break;
508  case 12: ParticleStatusName = "kIStDISPreFragmHadronicState"; break;
509  case 13: ParticleStatusName = "kIStPreDecayResonantState"; break;
510  case 14: ParticleStatusName = "kIStHadronInTheNucleus"; break;
511  case 15: ParticleStatusName = "kIStFinalStateNuclearRemnant"; break;
512  case 16: ParticleStatusName = "kIStNucleonClusterTarget"; break;
513  default: ParticleStatusName = "Status Unknown";
514  }
515  return ParticleStatusName;
516  }
517 
518  //......................................................................
519  std::string GENIEGen::ReactionChannel(int ccnc, int mode)
520  {
521  std::string ReactionChannelName = " ";
522 
523  if (ccnc == 0)
524  ReactionChannelName = "kCC";
525  else if (ccnc == 1)
526  ReactionChannelName = "kNC";
527  else
528  std::cout << "Current mode unknown!! " << std::endl;
529 
530  if (mode == 0)
531  ReactionChannelName += "_kQE";
532  else if (mode == 1)
533  ReactionChannelName += "_kRes";
534  else if (mode == 2)
535  ReactionChannelName += "_kDIS";
536  else if (mode == 3)
537  ReactionChannelName += "_kCoh";
538  else
539  std::cout << "interaction mode unknown!! " << std::endl;
540 
541  return ReactionChannelName;
542  }
543 
544  //......................................................................
546  {
547  // Decide which histograms to put the spectrum in
548  int id = -1;
549  if (mc.GetNeutrino().CCNC() == simb::kCC) {
550  fCCMode->Fill(mc.GetNeutrino().Mode());
551  if (mc.GetNeutrino().Nu().PdgCode() == 12)
552  id = 0;
553  else if (mc.GetNeutrino().Nu().PdgCode() == -12)
554  id = 1;
555  else if (mc.GetNeutrino().Nu().PdgCode() == 14)
556  id = 2;
557  else if (mc.GetNeutrino().Nu().PdgCode() == -14)
558  id = 3;
559  else
560  return;
561  }
562  else {
563  fNCMode->Fill(mc.GetNeutrino().Mode());
564  if (mc.GetNeutrino().Nu().PdgCode() > 0)
565  id = 4;
566  else
567  id = 5;
568  }
569  if (id == -1) abort();
570 
571  // Fill the specta histograms
572  fGenerated[id]->Fill(mc.GetNeutrino().Nu().E());
573 
576  simb::MCNeutrino mcnu = mc.GetNeutrino();
577  const simb::MCParticle nu = mcnu.Nu();
578 
579  fVertexX->Fill(nu.Vx());
580  fVertexY->Fill(nu.Vy());
581  fVertexZ->Fill(nu.Vz());
582 
583  fVertexXY->Fill(nu.Vx(), nu.Vy());
584  fVertexXZ->Fill(nu.Vz(), nu.Vx());
585  fVertexYZ->Fill(nu.Vz(), nu.Vy());
586 
587  double mom = nu.P();
588  if (std::abs(mom) > 0.) {
589  fDCosX->Fill(nu.Px() / mom);
590  fDCosY->Fill(nu.Py() / mom);
591  fDCosZ->Fill(nu.Pz() / mom);
592  }
593 
594  MF_LOG_DEBUG("GENIEInteractionInformation")
595  << std::endl
596  << "REACTION: " << ReactionChannel(mc.GetNeutrino().CCNC(), mc.GetNeutrino().Mode())
597  << std::endl
598  << "-----------> Particles in the Stack = " << mc.NParticles() << std::endl
599  << std::setiosflags(std::ios::left) << std::setw(20) << "PARTICLE"
600  << std::setiosflags(std::ios::left) << std::setw(32) << "STATUS" << std::setw(18) << "E (GeV)"
601  << std::setw(18) << "m (GeV/c2)" << std::setw(18) << "Ek (GeV)" << std::endl
602  << std::endl;
603 
604  const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
605 
606  // Loop over the particle stack for this event
607  for (int i = 0; i < mc.NParticles(); ++i) {
609  std::string name = databasePDG->GetParticle(part.PdgCode())->GetName();
610  int code = part.StatusCode();
611  std::string status = ParticleStatus(code);
612  double mass = part.Mass();
613  double energy = part.E();
614  double Ek = (energy - mass); // Kinetic Energy (GeV)
615  if (status == "kIStStableFinalState" || status == "kIStHadronInTheNucleus")
616  MF_LOG_DEBUG("GENIEFinalState")
617  << std::setiosflags(std::ios::left) << std::setw(20) << name
618  << std::setiosflags(std::ios::left) << std::setw(32) << status << std::setw(18) << energy
619  << std::setw(18) << mass << std::setw(18) << Ek << std::endl;
620  else
621  MF_LOG_DEBUG("GENIEFinalState")
622  << std::setiosflags(std::ios::left) << std::setw(20) << name
623  << std::setiosflags(std::ios::left) << std::setw(32) << status << std::setw(18) << energy
624  << std::setw(18) << mass << std::endl;
625  }
626 
627  if (mc.GetNeutrino().CCNC() == simb::kCC) {
628 
631  for (int i = 0; i < mc.NParticles(); ++i) {
633  if (abs(part.PdgCode()) == 11) {
634  fEMomentum->Fill(part.P());
635  fEDCosX->Fill(part.Px() / part.P());
636  fEDCosY->Fill(part.Py() / part.P());
637  fEDCosZ->Fill(part.Pz() / part.P());
638  fECons->Fill(nu.E() - part.E());
639  break;
640  }
641  else if (abs(part.PdgCode()) == 13) {
642  fMuMomentum->Fill(part.P());
643  fMuDCosX->Fill(part.Px() / part.P());
644  fMuDCosY->Fill(part.Py() / part.P());
645  fMuDCosZ->Fill(part.Pz() / part.P());
646  fECons->Fill(nu.E() - part.E());
647  break;
648  }
649  } // end loop over particles
650  } //end if CC interaction
651 
652  return;
653  }
654 
655 }
656 
657 namespace evgen {
658 
660 
661 }
TH1F * fVertexX
vertex location of generated events in x
std::string ParticleStatus(int StatusCode)
Float_t x
Definition: compare.C:6
double E(const int i=0) const
Definition: MCParticle.h:234
genie::GFluxI * GetFluxDriver(bool base=true)
Definition: GENIEHelper.h:91
TH1F * fMuDCosZ
direction cosine of outgoing mu in z
TH2F * fVertexXZ
vertex location in xz
int PdgCode() const
Definition: MCParticle.h:213
int CCNC() const
Definition: MCNeutrino.h:148
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:77
double Py(const int i=0) const
Definition: MCParticle.h:232
void DumpGTruth(Stream &&out, simb::GTruth const &truth, std::string indent, std::string firstIndent)
Dumps the content of the GENIE truth in the output stream.
Definition: MCDumpers.h:377
constexpr auto fullRun()
TH1F * fDeltaE
difference in neutrino energy from MCTruth::Enu() vs TParticle
TH2F * fVertexXY
vertex location in xy
void beginSubRun(art::SubRun &sr)
std::string ReactionChannel(int ccnc, int mode)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void beginRun(art::Run &run)
Unknown view.
Definition: geo_types.h:142
TH1F * fMuDCosY
direction cosine of outgoing mu in y
Float_t y
Definition: compare.C:6
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
bool Sample(simb::MCTruth &truth, simb::MCFlux &flux, simb::GTruth &gtruth)
EDProducer(fhicl::ParameterSet const &pset)
Definition: EDProducer.cc:6
Length_t DetHalfWidth(TPCID const &tpcid=tpc_zero) const
Returns the half width of the active volume of the specified TPC.
Double_t z
Definition: plot.C:276
double Px(const int i=0) const
Definition: MCParticle.h:231
constexpr auto abs(T v)
Returns the absolute value of the argument.
std::string const & ROOTFile() const
Returns the full directory path to the geometry file source.
Definition: GeometryCore.h:184
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Run.h:121
offset to account for adding in Nuance codes to this enum
Definition: MCNeutrino.h:95
void endSubRun(art::SubRun &sr)
int NParticles() const
Definition: MCTruth.h:75
TGeoManager * ROOTGeoManager() const
Access to the ROOT geometry description manager.
Length_t DetLength(TPCID const &tpcid=tpc_zero) const
Returns the length of the active volume of the specified TPC.
TH1F * fNCMode
CC interaction mode.
object containing MC flux information
Definition: Run.h:37
::sim::BeamType_t fBeamType
The width of a simulated "beam gate".
NuMI.
Definition: BeamTypes.h:12
double TotalMass() const
Returns the total mass [kg] of the specified volume (default: world).
Definition: GeometryCore.h:319
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: Event.h:77
TH1F * fDCosZ
direction cosine in z
int InteractionType() const
Definition: MCNeutrino.h:150
TString part[npart]
Definition: Style.C:32
void produce(art::Event &evt)
TH1F * fEDCosY
direction cosine of outgoing e in y
Utility functions to print MC truth information.
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
GENIEGen(fhicl::ParameterSet const &pset)
long seed
Definition: chem4.cc:67
int fPassEmptySpills
whether or not to kill evnets with no interactions
TH1F * fCCMode
CC interaction mode.
double fPrevTotGoodPOT
Total good POT from subruns previous to current subrun.
evgb::GENIEHelper * fGENIEHelp
GENIEHelper object.
double P(const int i=0) const
Definition: MCParticle.h:235
double fGlobalTimeOffset
keep track of how long it takes to run the job
TH1F * fEDCosX
direction cosine of outgoing e in x
double energy
Definition: plottest35.C:25
TH1F * fMuDCosX
direction cosine of outgoing mu in x
void DumpMCTruth(Stream &&out, simb::MCTruth const &truth, unsigned int pointsPerLine, std::string indent, std::string firstIndent)
Dumps the content of the specified MC truth in the output stream.
Definition: MCDumpers.h:338
TH1F * fGenerated[6]
Spectra as generated.
Wrapper for generating neutrino interactions with GENIE.
std::vector< double > fVtxPosHistRange
use defined hist range; it is useful to have for asymmetric ranges like in DP FD. ...
double fPrevTotPOT
The type of beam.
An art service to assist in the distribution of guaranteed unique seeds to all engines within an art ...
PutHandle< PROD > put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: SubRun.h:126
TH1F * fDCosY
direction cosine in y
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
constexpr auto subRunFragment()
double Vx(const int i=0) const
Definition: MCParticle.h:222
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
TH1F * fEMomentum
momentum of outgoing electrons
BNB.
Definition: BeamTypes.h:11
TH2F * fVertexYZ
vertex location in yz
Utility object to perform functions of association.
double TotalExposure() const
Definition: GENIEHelper.h:71
TH1F * fECons
histogram to determine if energy is conserved in the event
std::optional< T > get_if_present(std::string const &key) const
Definition: ParameterSet.h:267
#define MF_LOG_DEBUG(id)
double Pz(const int i=0) const
Definition: MCParticle.h:233
double Vz(const int i=0) const
Definition: MCParticle.h:224
TH1F * fEDCosZ
direction cosine of outgoing e in z
TH1F * fDCosX
direction cosine in x
Length_t DetHalfHeight(TPCID const &tpcid=tpc_zero) const
Returns the half height of the active volume of the specified TPC.
std::string const & DetectorName() const
Returns a string with the name of the detector, as configured.
Definition: GeometryCore.h:203
bool NeutrinoSet() const
Definition: MCTruth.h:78
TCEvent evt
Definition: DataStructs.cxx:8
Event generator information.
Definition: MCTruth.h:32
TStopwatch fStopwatch
Event generator information.
Definition: MCNeutrino.h:18
Namespace collecting geometry-related classes utilities.
TH1F * fMuMomentum
momentum of outgoing muons
int Mode() const
Definition: MCNeutrino.h:149
Event Generation using GENIE, cosmics or single particles.
std::string ParticleStatusName(int code)
Describes the status of a particle (simb::MCParticle::StatusCode()).
TH1F * fVertexY
vertex location of generated events in y
BeamType_t
Defines category of beams to be stored in sim::BeamGateInfo.
Definition: BeamTypes.h:9
void FillHistograms(simb::MCTruth mc)
TH1F * fVertexZ
vertex location of generated events in z
double Vy(const int i=0) const
Definition: MCParticle.h:223
void put(std::string const &key)
art framework interface to geometry description
A module to check the results from the Monte Carlo generator.
double fRandomTimeOffset
The start of a simulated "beam gate".