MakeGENIEReweightTable_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Use NuReweight to compute +/-1,2sigma shifts for all systematics
3 /// \author Christopher Backhouse - bckhouse@caltech.edu
4 ////////////////////////////////////////////////////////////////////////
5 
6 // GENIE: Macros in here conflict with art's MessageFacility. So include it
7 // first and then unset them again
8 #include "GENIE/Framework/Messenger/Messenger.h"
9 #undef LOG_INFO
10 #undef LOG_DEBUG
11 #undef LOG_ERROR
12 
13 #include <cassert>
14 #include <ctime>
15 #include <iomanip>
16 #include <memory>
17 #include <string>
18 #include <unordered_set>
19 
20 #include "TH1D.h"
21 
28 #include "cetlib_except/exception.h"
29 #include "fhiclcpp/ParameterSet.h"
31 
32 #include "Geometry/Geometry.h"
34 #include "nugen/EventGeneratorBase/GENIE/GENIE2ART.h"
35 #include "nugen/NuReweight/art/NuReweight.h"
40 #include "Utilities/AssociationUtil.h"
41 
42 #include "GENIE/RwFramework/GReWeight.h"
43 
44 namespace rwgt
45 {
46  // It's a shame such a thing isn't officially defined in any of the headers
47  const int kMaxReweightIdx = rwgt::fReweightAxFFCCQEshape;
48 
49  /// \brief Use NuReweight to compute +/-1,2sigma shifts for all systematics
51  {
52  public:
53  explicit MakeGENIEReweightTable(const fhicl::ParameterSet& pset);
54  virtual ~MakeGENIEReweightTable();
55 
56  void produce(art::Event& evt) override;
57  void beginJob() override;
58  void endJob() override;
59 
60  private:
63 
64  bool fConsiderNDRock; ///< Include rock interaction in the ND?
65  bool fConsiderFDRock; ///< Include rock interaction in the FD?
66 
67  bool fTimeKnobs; ///< Keep tabs on the reweight knobs & measure the time they take?
68 
70 
71  std::vector<std::string> fKnobNamesToSkip;
72  // organized by tune name
73  std::unordered_map<std::string, std::unordered_set<rwgt::ReweightLabel_t, std::hash<int>>> fKnobsToSkip;
74 
75  std::vector<rwgt::NuReweight*> fRWsMinus2;
76  std::vector<rwgt::NuReweight*> fRWsMinus1;
77  std::vector<rwgt::NuReweight*> fRWsPlus1;
78  std::vector<rwgt::NuReweight*> fRWsPlus2;
79 
80  std::map<rwgt::ReweightLabel_t, TH1D*> fKnobTimes;
82  };
83 
84  //___________________________________________________________________________
86  : fGeneratorToken(consumes<std::vector<simb::MCTruth>>(pset.get<std::string>("GeneratorLabel"))),
87  fGeneratorLabel(pset.get<std::string>("GeneratorLabel")),
88  fConsiderNDRock(pset.get<bool>("ConsiderNDRock")),
89  fConsiderFDRock(pset.get<bool>("ConsiderFDRock")),
90  fTimeKnobs(pset.get<bool>("TimeKnobs")),
91  fSkipCCQEWgtsForHyperons(pset.get<bool>("SkipCCQEWgtsForHyperons"))
92 // fKnobNamesToSkip(pset.get<std::unordered_map<std::string, std::vector<std::string>>>("KnobsToSkip"))
93  {
94  produces<std::vector<rwgt::GENIEReweightTable>>();
95 
96  produces<art::Assns<rwgt::GENIEReweightTable, simb::MCTruth>>();
97 
98  mf::LogInfo("MakeGENIEReweightTable") << "Skipping the following reweight knobs:";
99  auto knobsToSkipLists = pset.get<fhicl::ParameterSet>("KnobsToSkip");
100  if (!knobsToSkipLists.has_key("default"))
101  {
102  mf::LogWarning("MakeGENIEReweightTable") << "No 'default' key in KnobsToSkip table. Assuming empty list...";
103  fKnobsToSkip["default"]; // default-initializes...
104  }
105  for (const auto & knobListKey : knobsToSkipLists.get_names())
106  {
107  mf::LogInfo("MakeGENIEReweightTable") << " for tune '" << knobListKey << "'" << std::endl;
108  for (const auto &knobName : knobsToSkipLists.get<std::vector<std::string>>(knobListKey))
109  {
111  throw cet::exception("MakeGENIEReweightTable",
112  "Unknown GENIE reweight knob name: " + knobName);
113 
114  this->fKnobsToSkip[knobListKey].insert(rwgt::REWEIGHT_KNOB_INDICES.at(knobName));
115  mf::LogInfo("MakeGENIEReweightTable")
116  << " " << knobName
117  << " (index: " << rwgt::REWEIGHT_KNOB_INDICES.at(knobName) << ")";
118  }
119  }
120  }
121 
122  //___________________________________________________________________________
124  {
125  for(rwgt::NuReweight* rw: fRWsMinus2) delete rw;
126  for(rwgt::NuReweight* rw: fRWsMinus1) delete rw;
127  for(rwgt::NuReweight* rw: fRWsPlus1) delete rw;
128  for(rwgt::NuReweight* rw: fRWsPlus2) delete rw;
129  }
130 
131  //___________________________________________________________________________
133  {
134  mf::LogInfo("MakeGENIEReweightTable") << "Initializing GENIE reweight calculators...";
135 
136  // for GENIE 3, need to set the tune before you do anything.
137  // (this function is a no-op for GENIE 2 so it's safe anyway.)
138  // note that this sets the tune to the one that the genie_xsec UPS product
139  // is configured to load. 99.9999% of the time that's the one you want,
140  // but if you're trying to create reweights for MC that was generated
141  // in a different release than the one you're using,
142  // you may need to add some FCL handles to pass options to this function.
143  evgb::SetEventGeneratorListAndTune();
144 
145  // Don't need to see every tiny thing these modules are thinking...
147  log4cpp::Priority::WARN);
149  log4cpp::Priority::WARN);
151  log4cpp::Priority::WARN);
152 
153  // The way NuReweight constructs the EventRecord it inevitably trips this
156 
157  genie::Messenger::Instance()->SetPriorityLevel("GENIEReweight",
158  log4cpp::Priority::WARN);
159 
160  fRWsMinus2.resize(kMaxReweightIdx+1, 0);
161  fRWsMinus1.resize(kMaxReweightIdx+1, 0);
162  fRWsPlus1. resize(kMaxReweightIdx+1, 0);
163  fRWsPlus2. resize(kMaxReweightIdx+1, 0);
164 
165  // Start at 1. Zero is kReweightNull. Not worth dealing with shifting all
166  // the indices by one everywhere, just leave it empty.
167  for(int i = 1; i <= kMaxReweightIdx; ++i){
168  for(int sigma: {-2, -1, +1, +2}){
169  rwgt::NuReweight* rw = new rwgt::NuReweight;
170 
171  switch(sigma){
172  case -2: fRWsMinus2[i] = rw; break;
173  case -1: fRWsMinus1[i] = rw; break;
174  case +1: fRWsPlus1 [i] = rw; break;
175  case +2: fRWsPlus2 [i] = rw; break;
176  }
177  double shift = sigma;
178 
179  //Catch the cases where shape parameters are called instead of rate parameters.
180  switch(i) {
181  case rwgt::fReweightNormCCQE:
182  case rwgt::fReweightNormCCQEenu:
183  case rwgt::fReweightMaCCQEshape:
184  rw->MaQEshape();
185  break;
186 
187  case rwgt::fReweightNormCCRES:
188  case rwgt::fReweightMaCCRESshape:
189  case rwgt::fReweightMvCCRESshape:
190  rw->CCRESshape();
191  break;
192 
193  case rwgt::fReweightNormNCRES:
194  case rwgt::fReweightMaNCRESshape:
195  case rwgt::fReweightMvNCRESshape:
196  rw->NCRESshape();
197  break;
198 
199  case rwgt::fReweightAhtBYshape:
200  case rwgt::fReweightBhtBYshape:
201  case rwgt::fReweightCV1uBYshape:
202  case rwgt::fReweightCV2uBYshape:
203  rw->DIS_BYshape();
204  break;
205 
206  case rwgt::fReweightCCQEMomDistroFGtoSF:
207  case rwgt::fReweightTheta_Delta2Npi:
208  // These are continuous knobs.
209  // Encode -2sigma = .25, -1sigma = .5, +1sigma = .75, +2sigma = 1
210  shift = (shift < 0) ? (shift+3.0)/4.0 : (shift + 2)/4.0;
211  break;
212 
213  case rwgt::fReweightRnubarnuCC:
214 // case rwgt::fReweightDISNuclMod: // disabled in GENIE v3. may return later?
215  // These are binary switches. Make the negative sigmas be "off" and
216  // the positive be "on".
217  shift = (shift < 0) ? 0 : 1;
218  break;
219  }
220 
221  rw->AddReweightValue(rwgt::ReweightLabel_t(i), shift);
222  rw->Configure();
223  }
224  }
225 
226  mf::LogInfo("MakeGENIEReweightTable") << "Done initializing GENIE reweight calculators...";
227 
228  }
229  //___________________________________________________________________________
231  {
232  std::cout << "Statistics regarding time spent in each reweight knob:" << std::endl;
233  std::cout << "knob # | total time (s) | mean time (s) | median time (s)" << std::endl;
234  std::cout << "---------------------------------------------------------" << std::endl;
235 
236  for (auto & knobTimePair : this->fKnobTimes)
237  {
238  std::cout << std::setw(7) << knobTimePair.first << "|";
239  TH1D * knobHist = knobTimePair.second;
240 
241  double sum = 0;
242  for (int bin = 1; bin < knobHist->GetNbinsX()+1; bin++)
243  sum += pow(10, knobHist->GetBinCenter(bin)) * knobHist->GetBinContent(bin);
244  std::cout << std::setw(16) << sum << "|";
245 
246  std::cout << std::setw(15);
247  if (knobHist->GetEntries() > 0)
248  std::cout << pow(10, knobHist->GetMean());
249  else
250  std::cout << "-";
251  std::cout << "|";
252 
253  std::cout << std::setw(16);
254  int n = knobHist->GetXaxis()->GetNbins();
255  std::vector<double> x(n);
256  knobHist->GetXaxis()->GetCenter( &x[0] );
257  const double * y = knobHist->GetArray();
258  std::cout << pow(10, TMath::Median(n, &x[0], &y[1]));
259 
260  std::cout << std::endl;
261  }
262  }
263 
264  //___________________________________________________________________________
266  {
268 
269  std::unique_ptr<std::vector<rwgt::GENIEReweightTable>> tablecol(new std::vector<rwgt::GENIEReweightTable>);
270 
271  std::unique_ptr<art::Assns<rwgt::GENIEReweightTable, simb::MCTruth>> assns(new art::Assns<rwgt::GENIEReweightTable, simb::MCTruth>);
272 
273 
275  evt.getByToken(fGeneratorToken, mctruths);
276  assert(!mctruths.failedToGet());
277 
278  art::FindManyP<simb::GTruth> fmp(mctruths, evt, fGeneratorLabel);
279  assert(fmp.isValid());
280 
281  for(unsigned int truthIdx = 0; truthIdx < mctruths->size(); ++truthIdx){
282  const std::vector<art::Ptr<simb::GTruth>> gtruths = fmp.at(truthIdx);
283  assert(gtruths.size() == 1);
284 
285  const simb::MCTruth& mctruth = (*mctruths)[truthIdx];
286 
287  // obviously we won't have a reweight table if this isn't a GENIE event
288  if (mctruth.GeneratorInfo().generator != simb::Generator_t::kGENIE)
289  continue;
290 
291  const simb::GTruth& gtruth = *gtruths[0];
292 
293  mf::LogDebug("MakeGENIEReweightTable") << "Reweighting event:\n"
294  << "Interaction type: " << (mctruth.GetNeutrino().CCNC() == simb::kCC ? "CC" : "NC") << " " << ", mode = " << mctruth.GetNeutrino().Mode();
295 
296  // If we want to skip rock events in this detector
297  if((geom->DetId() == novadaq::cnv::kNEARDET && !fConsiderNDRock) ||
298  (geom->DetId() == novadaq::cnv::kFARDET && !fConsiderFDRock)){
299  // Where's the true vertex?
300  const TVector3 vtx = mctruth.GetNeutrino().Nu().Position().Vect();
301 
302  // If the vertex is outside the detector, skip this event
303  if(!geom->isInsideDetectorBigBox(vtx.X(), vtx.Y(), vtx.Z())) continue;
304  }
305 
306  // Signal of a nu/e interaction, which can cause trouble
307  if(gtruth.fgQ2 < 0) continue;
308 
310 
311  std::string tuneName = "default";
312  const auto * knobsToSkip = &this->fKnobsToSkip.at("default");
313  if (!mctruth.GeneratorInfo().generatorVersion.empty()
314  && (mctruth.GeneratorInfo().generatorVersion != "unknown")
315  && std::atoi(&mctruth.GeneratorInfo().generatorVersion[0]) >= 3
316  && mctruth.GeneratorInfo().generatorConfig.find("tune") != mctruth.GeneratorInfo().generatorConfig.end())
317  {
318  tuneName = mctruth.GeneratorInfo().generatorConfig.at("tune");
319  if (this->fKnobsToSkip.find(tuneName) == this->fKnobsToSkip.end())
320  {
321  mf::LogWarning("MakeGENIEReweightTable") << "No tune-specific knob list for tune '" << tuneName
322  << "'. Using default list.";
323  tuneName = "default";
324  }
325  else
326  knobsToSkip = &this->fKnobsToSkip.at(tuneName);
327  }
328  mf::LogDebug("MakeGENIEReweightTable") << "Using knob list for tune: '" << tuneName << "'";
329  bool hyperonCCQE = ((gtruth.fIsStrange || gtruth.fIsCharm)
330  && mctruth.GetNeutrino().Mode() == simb::kQE
331  && mctruth.GetNeutrino().CCNC() == simb::kCC);
332  for(int systIdx = 1; systIdx <= kMaxReweightIdx; ++systIdx){
333  if (knobsToSkip->find(ReweightLabel_t(systIdx)) != knobsToSkip->end())
334  continue;
335 
336  // CCQE knobs sometimes barf on CCQE-strange antinu events. assume it probably does the same on CCQE-charm nu events too
337  if (hyperonCCQE && fSkipCCQEWgtsForHyperons)
338  {
339  genie::rew::GReWeight * wgtCalc = fRWsPlus1[systIdx]->WeightCalculator();
340  bool skipThisKnob = false;
341  for (const auto & wgtrName : wgtCalc->WghtCalcNames())
342  {
343  const genie::rew::GReWeightI * wgtr = wgtCalc->WghtCalc(wgtrName);
344  if (wgtr->AppliesTo(genie::kScQuasiElastic, true))
345  {
346  skipThisKnob = true;
347  break;
348  }
349  }
350  if (skipThisKnob)
351  {
352  mf::LogWarning("MakeGENIEReweightTable") << "Skipping knob "
353  << genie::rew::GSyst::AsString(static_cast<genie::rew::GSyst_t>(systIdx))
354  << " for QES event:";
355  mf::LogWarning("MakeGENIEReweightTable") << gtruth;
356  continue;
357  }
358  }
359 
360  std::clock_t start = std::clock();
361 
362  const double wm2 = fRWsMinus2[systIdx]->CalcWeight(mctruth, gtruth);
363  const double wm1 = fRWsMinus1[systIdx]->CalcWeight(mctruth, gtruth);
364  const double wp1 = fRWsPlus1 [systIdx]->CalcWeight(mctruth, gtruth);
365  const double wp2 = fRWsPlus2 [systIdx]->CalcWeight(mctruth, gtruth);
366 
367  double duration = ( std::clock() - start ) / (double) CLOCKS_PER_SEC;
368 
369  rwgt::EReweightLabel systLabel = rwgt::EReweightLabel(systIdx);
370  if (this->fTimeKnobs)
371  {
372  if (this->fKnobTimes.find(systLabel) == this->fKnobTimes.end())
373  this->fKnobTimes[systLabel] = this->fFileService->make<TH1D>( Form("knob_%d", systIdx), Form(";Log_{10}(Time reweight knob %d used in event) (s)", systIdx), 100, -7, 3);
374 
375  if (duration <= 0)
376  this->fKnobTimes[systLabel]->Fill(-7);
377  else
378  this->fKnobTimes[systLabel]->Fill(log10(duration));
379  }
380 
381  table.SetWeights(systIdx, wm2, wm1, wp1, wp2);
382  } // end for systIdx
383 
384  tablecol->push_back(table);
385 
386  util::CreateAssn(*this, evt, *tablecol,
387  art::Ptr<simb::MCTruth>(mctruths, truthIdx),
388  *assns);
389  } // end for truthIdx
390 
391 
392  evt.put(std::move(tablecol));
393  evt.put(std::move(assns));
394  }
395 
397 
398 } // namespace
399 
400 
401 // If the GENIE enum values are ever out of sync between generation-time and
402 // analysis time that can cause all sorts of confusion. Hopefully GENIE never
403 // change their enum values, but if it does happen then at least this will warn
404 // us.
405 
406 // Message has to be a string literal so I can't do "const std::string" here.
407 #define m "GENIE reweight enumeration values have changed."
408 
409 static_assert(rwgt::kReweightNull == 0, m);
410 static_assert(rwgt::fReweightMaNCEL == 1, m);
411 static_assert(rwgt::fReweightEtaNCEL == 2, m);
412 static_assert(rwgt::fReweightNormCCQE == 3, m);
413 static_assert(rwgt::fReweightNormCCQEenu == 4, m);
414 static_assert(rwgt::fReweightMaCCQEshape == 5, m);
415 static_assert(rwgt::fReweightMaCCQE == 6, m);
416 static_assert(rwgt::fReweightVecCCQEshape == 7, m);
417 static_assert(rwgt::fReweightNormCCRES == 8, m);
418 static_assert(rwgt::fReweightMaCCRESshape == 9, m);
419 static_assert(rwgt::fReweightMvCCRESshape == 10, m);
420 static_assert(rwgt::fReweightMaCCRES == 11, m);
421 static_assert(rwgt::fReweightMvCCRES == 12, m);
422 static_assert(rwgt::fReweightNormNCRES == 13, m);
423 static_assert(rwgt::fReweightMaNCRESshape == 14, m);
424 static_assert(rwgt::fReweightMvNCRESshape == 15, m);
425 static_assert(rwgt::fReweightMaNCRES == 16, m);
426 static_assert(rwgt::fReweightMvNCRES == 17, m);
427 static_assert(rwgt::fReweightMaCOHpi == 18, m);
428 static_assert(rwgt::fReweightR0COHpi == 19, m);
429 static_assert(rwgt::fReweightRvpCC1pi == 20, m);
430 static_assert(rwgt::fReweightRvpCC2pi == 21, m);
431 static_assert(rwgt::fReweightRvpNC1pi == 22, m);
432 static_assert(rwgt::fReweightRvpNC2pi == 23, m);
433 static_assert(rwgt::fReweightRvnCC1pi == 24, m);
434 static_assert(rwgt::fReweightRvnCC2pi == 25, m);
435 static_assert(rwgt::fReweightRvnNC1pi == 26, m);
436 static_assert(rwgt::fReweightRvnNC2pi == 27, m);
437 static_assert(rwgt::fReweightRvbarpCC1pi == 28, m);
438 static_assert(rwgt::fReweightRvbarpCC2pi == 29, m);
439 static_assert(rwgt::fReweightRvbarpNC1pi == 30, m);
440 static_assert(rwgt::fReweightRvbarpNC2pi == 31, m);
441 static_assert(rwgt::fReweightRvbarnCC1pi == 32, m);
442 static_assert(rwgt::fReweightRvbarnCC2pi == 33, m);
443 static_assert(rwgt::fReweightRvbarnNC1pi == 34, m);
444 static_assert(rwgt::fReweightRvbarnNC2pi == 35, m);
445 static_assert(rwgt::fReweightAhtBY == 36, m);
446 static_assert(rwgt::fReweightBhtBY == 37, m);
447 static_assert(rwgt::fReweightCV1uBY == 38, m);
448 static_assert(rwgt::fReweightCV2uBY == 39, m);
449 static_assert(rwgt::fReweightAhtBYshape == 40, m);
450 static_assert(rwgt::fReweightBhtBYshape == 41, m);
451 static_assert(rwgt::fReweightCV1uBYshape == 42, m);
452 static_assert(rwgt::fReweightCV2uBYshape == 43, m);
453 static_assert(rwgt::fReweightNormDISCC == 44, m);
454 static_assert(rwgt::fReweightRnubarnuCC == 45, m);
455 static_assert(rwgt::fReweightDISNuclMod == 46, m);
456 static_assert(rwgt::fReweightNC == 47, m);
457 static_assert(rwgt::fReweightAGKY_xF1pi == 48, m);
458 static_assert(rwgt::fReweightAGKY_pT1pi == 49, m);
459 static_assert(rwgt::fReweightFormZone == 50, m);
460 static_assert(rwgt::fReweightMFP_pi == 51, m);
461 static_assert(rwgt::fReweightMFP_N == 52, m);
462 static_assert(rwgt::fReweightFrCEx_pi == 53, m);
463 //static_assert(rwgt::fReweightFrElas_pi == 54, m);
464 static_assert(rwgt::fReweightFrInel_pi == 55, m);
465 static_assert(rwgt::fReweightFrAbs_pi == 56, m);
466 static_assert(rwgt::fReweightFrPiProd_pi == 57, m);
467 static_assert(rwgt::fReweightFrCEx_N == 58, m);
468 //static_assert(rwgt::fReweightFrElas_N == 59, m);
469 static_assert(rwgt::fReweightFrInel_N == 60, m);
470 static_assert(rwgt::fReweightFrAbs_N == 61, m);
471 static_assert(rwgt::fReweightFrPiProd_N == 62, m);
472 static_assert(rwgt::fReweightCCQEPauliSupViaKF == 63, m);
473 static_assert(rwgt::fReweightCCQEMomDistroFGtoSF == 64, m);
474 static_assert(rwgt::fReweightBR1gamma == 65, m);
475 static_assert(rwgt::fReweightBR1eta == 66, m);
476 static_assert(rwgt::fReweightTheta_Delta2Npi == 67, m);
477 static_assert(rwgt::fReweightZNormCCQE == 68, m);
478 static_assert(rwgt::fReweightZExpA1CCQE == 69, m);
479 static_assert(rwgt::fReweightZExpA2CCQE == 70, m);
480 static_assert(rwgt::fReweightZExpA3CCQE == 71, m);
481 static_assert(rwgt::fReweightZExpA4CCQE == 72, m);
482 static_assert(rwgt::fReweightAxFFCCQEshape == 73, m);
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
bool fTimeKnobs
Keep tabs on the reweight knobs & measure the time they take?
bool fConsiderNDRock
Include rock interaction in the ND?
void SetWeights(unsigned int i, double m2, double m1, double p1, double p2)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
const std::unordered_map< std::string, ReweightLabel_t > REWEIGHT_KNOB_INDICES
std::vector< rwgt::NuReweight * > fRWsPlus1
constexpr T pow(T x)
Definition: pow.h:75
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
Use NuReweight to compute +/-1,2sigma shifts for all systematics.
DEFINE_ART_MODULE(TestTMapFile)
std::vector< rwgt::NuReweight * > fRWsMinus2
bool fIsStrange
strange production // added version 13
Definition: GTruth.h:73
std::vector< rwgt::NuReweight * > fRWsPlus2
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
const art::ProductToken< std::vector< simb::MCTruth > > fGeneratorToken
::xsd::cxx::tree::duration< char, simple_type > duration
Definition: Database.h:188
Far Detector at Ash River, MN.
std::vector< rwgt::NuReweight * > fRWsMinus1
Use NuReweight to compute +/-1,2sigma shifts for all systematics.
Definition: FillTruth.h:18
void resize(T &x, std::vector< int > dims)
Definition: resize.hpp:41
static Messenger * Instance(void)
Definition: Messenger.cxx:71
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
#define ERROR
Near Detector in the NuMI cavern.
bool fIsCharm
did the interaction produce a charmed hadron?
Definition: GTruth.h:71
bool isInsideDetectorBigBox(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Big Box?
std::vector< std::string > fKnobNamesToSkip
float bin[41]
Definition: plottest35.C:14
This class describes a particle created in the detector Monte Carlo simulation.
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
T * make(ARGS...args) const
std::map< rwgt::ReweightLabel_t, TH1D * > fKnobTimes
T log10(T number)
Definition: d0nt_math.hpp:120
void SetPriorityLevel(const char *stream, log4cpp::Priority::Value p)
Definition: Messenger.cxx:110
std::unordered_map< std::string, std::unordered_set< rwgt::ReweightLabel_t, std::hash< int > > > fKnobsToSkip
MakeGENIEReweightTable(const fhicl::ParameterSet &pset)
double fgQ2
Definition: GTruth.h:61
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
Store +/-1,2sigma shifts for all GENIE reweighting systematics.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
assert(nhit_max >=nhit_nbins)
const char * AsString(Resonance_t res)
resonance id -> string
Event generator information.
Definition: MCTruth.h:32
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
Double_t sum
Definition: plot.C:31
ProductToken< T > consumes(InputTag const &)
bool fConsiderFDRock
Include rock interaction in the FD?
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196
void produce(art::Event &evt) override
art::ServiceHandle< art::TFileService > fFileService
enum BeamMode string