CAFMaker_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief This module creates Common Analysis Files.
3 /// \author rocco@physics.umn.edu
4 /// \date July 2012
5 ////////////////////////////////////////////////////////////////////////
6 
9 #include "CAFMaker/FillPIDs.h"
11 #include "CAFMaker/FillReco.h"
12 #include "CAFMaker/FillTruth.h"
13 #include "CAFMaker/Utils.h"
14 
15 // C/C++ includes
16 #include <fenv.h>
17 #include <time.h>
18 #include <algorithm>
19 #include <cmath>
20 #include <iomanip>
21 #include <iostream>
22 #include <iterator>
23 #include <map>
24 #include <sstream>
25 #include <string>
26 #include <vector>
27 
28 #ifdef DARWINBUILD
29 #include <libgen.h>
30 #endif
31 
32 // ROOT includes
33 #include "TFile.h"
34 #include "TH1D.h"
35 #include "TTree.h"
36 
37 // Framework includes
52 #include "cetlib_except/exception.h"
53 #include "fhiclcpp/ParameterSet.h"
55 
56 #include <IFDH_service.h>
57 
58 // NOvASoft includes
59 // Geo/Cal/Util
60 #include "Calibrator/Calibrator.h"
62 #include "Geometry/Geometry.h"
63 #include "Geometry/LiveGeometry.h"
65 #include "GeometryObjects/Geo.h"
68 #include "NovaDAQConventions/DAQConventions.h"
69 #include "RawData/RawTrigger.h"
71 #include "Utilities/func/MathUtil.h"
72 
73 // DataQuality
75 
76 // Sim/BackTracker
77 #include "MCCheater/BackTracker.h"
78 #include "MCReweight/FluxWeights.h"
85 
86 // RecoBase
87 #include "RecoBase/CellHit.h"
88 #include "RecoBase/Cluster.h"
89 #include "RecoBase/Energy.h"
90 #include "RecoBase/FilterList.h"
91 #include "RecoBase/FitSum.h"
92 #include "RecoBase/Prong.h"
93 #include "RecoBase/RecoHit.h"
94 #include "RecoBase/Shower.h"
95 #include "RecoBase/Track.h"
96 #include "RecoBase/Vertex.h"
97 #include "RecoBase/WeightedProng.h"
98 
99 // Other Reco
100 #include "CVN/func/CVNFeatures.h"
101 #include "CVN/func/PixelMap.h"
103 #include "CVN/func/Result.h"
104 #include "CVN/func/RegResult.h"
105 #include "CVN/func/TrainingData.h"
106 #include "CosRej/CosRejObj.h"
107 #include "CosRej/NueCosRej.h"
108 #include "CosRej/TrkCntObj.h"
109 #include "LEM/PIDDetails.h"
111 #include "MEFinder/MEClusters.h"
113 #include "NCID/NCCosRej.h"
114 #include "NumuEnergy/NumuE.h"
115 #include "Preselection/PreselObj.h"
116 #include "Preselection/Veto.h"
117 #include "ReMId/ReMId.h"
118 #include "RecVarPID/RVP.h"
119 #include "RecoJMShower/EID.h"
120 #include "RecoJMShower/JMShower.h"
122 #include "ShowerLID/EventLID.h"
123 #include "ShowerLID/ShowerLID.h"
124 #include "ShowerLID/ShowerPID.h"
125 #include "ShowerLID/SliceLID.h"
128 #include "XnuePID/Xnue.h"
130 #include "TrackInfo/TrackInfoObj.h"
131 
132 // SummData
135 #include "SummaryData/POTSum.h"
136 #include "SummaryData/SpillData.h"
137 
138 // OscLib - for woscdumb
139 #include "OscLib/OscCalcDumb.h"
140 
141 // StandardRecord
144 // CAFMaker support
145 #include "CAFMaker/Blinding.h"
146 
147 #include "Utilities/AssociationUtil.h"
148 
149 // Note it is convention to NOT declare ANY "using namespace xxx" here
150 // especially "using namespace std"
151 namespace caf {
152 /// Module to create Common Analysis Files from ART files
153 class CAFMaker : public art::EDProducer {
154  public:
155  // Allows 'nova --print-description' to work
157 
158  explicit CAFMaker(const Parameters& params);
159  virtual ~CAFMaker();
160 
161  void produce(art::Event& evt) noexcept;
162 
164 
165  void beginJob();
166  void endJob();
167  virtual void beginRun(art::Run& r);
168  virtual void beginSubRun(art::SubRun& sr);
169  virtual void endSubRun(art::SubRun& sr);
170 
171  protected:
173 
175 
176  bool fIsRealData; // use this instead of evt.isRealData(), see init in
177  // produce(evt)
178  double fTotalPOT;
180  double fTotalEvents;
184  int fCycle;
185  int fBatch;
186  int fMask;
187 
188  TFile* fFile;
189  TTree* fRecTree;
190  TH1D* hPOT;
191  TH1D* hSinglePOT;
194  TH1D* hEvents;
195  TH1D* hLivetime;
201 
203  fDetID; ///< Detector ID in nova daq convention namespace typedef
204  Det_t fDet; ///< Detector ID in caf namespace typedef
205 
206  // Creating separate friend tree and branches for spill info.
207  // This tree will store one entry per event, even empty ones
208  // so that POT counting can be done right.
209  TTree* fSpillTree;
210 
211  //Seperate tree for CVN cosmic rejeciton. This operates on slices of time
212  //and not physics slices. So it must be stored seperately.
214 
215  TTree* fNuTree;
216 
217  void InitializeOutfile();
218 
219  void AddMetadataToFile(std::map<std::string, std::string> metadata);
220 
221  void AddMCTruthToVec(const cheat::NeutrinoEffPur* effPur,
222  std::vector<cheat::TrackIDE>& allTracks,
223  std::vector<cheat::TrackIDE>& sliceTracks,
224  std::vector<cheat::TrackIDE>& allTracksBirks,
225  std::vector<cheat::TrackIDE>& sliceTracksBirks,
226  std::vector<SRNeutrino>* vec, const art::Event& evt,
227  const SRSpill& spill);
228 
229  void drawSliceTruthTable(
230  const std::vector<std::vector<cheat::NeutrinoEffPur>>& sliceTruthTable,
231  const std::vector<int>& faveIds);
232  void drawSliceTruthTable(
233  const std::vector<std::vector<cheat::NeutrinoEffPur>>& sliceTruthTable,
234  const std::vector<int>& faveIds, const std::vector<int>& otherFaveIds);
235 
236  /// Equivalent of FindManyP except a return that is !isValid() prints a
237  /// messsage and aborts if StrictMode is true.
238  template <class T, class U>
239  art::FindManyP<T> FindManyPStrict(const U& from, const art::Event& evt,
240  const art::InputTag& label) const;
241 
242  /// \brief Retrieve an object from an association, with error handling
243  ///
244  /// This can go wrong in two ways: either the FindManyP itself is
245  /// invalid, or the result for the requested index is empty. In most
246  /// cases these have the same response, so conflating them here
247  /// saves redundancy elsewhere.
248  ///
249  /// \param fm The FindManyP object describing the association
250  /// \param idx Which element of the FindManyP to look it
251  /// \param[out] ret The product retrieved
252  /// \return Whether \a ret was filled
253  template <class T>
254  bool GetAssociatedProduct(const art::FindManyP<T>& fm, int idx, T& ret) const;
255 
256  /// Equivalent of evt.getByLabel(label, handle) except failedToGet
257  /// prints a message and aborts if StrictMode is true.
258  template <class T>
259  void GetByLabelStrict(const art::Event& evt, const std::string& label,
260  art::Handle<T>& handle) const;
261 
262  /// Equivalent of evt.getByLabel(label, handle) except failedToGet
263  /// prints a message.
264  template <class T>
265  void GetByLabelIfExists(const art::Event& evt, const std::string& label,
266  art::Handle<T>& handle) const;
267 
268  /// \param pset The parameter set
269  /// \param name Pass "foo.bar.baz" as {"foo", "bar", "baz"}
270  /// \param[out] ret Value of the key, not set if we return false
271  /// \return Whether the key was found
272  template <class T>
273  bool GetPsetParameter(const fhicl::ParameterSet& pset,
274  const std::vector<std::string>& name, T& ret) const;
275 
276  static bool EssentiallyEqual(double a, double b, double precision = 0.0001) {
277  return a <= (b + precision) && a >= (b - precision);
278  }
279 
280  static bool sortTrackLength(const SRTrack& a, const SRTrack& b) {
281  return a.len > b.len;
282  }
283 
285  const art::Ptr<rb::Track>& b) {
286  return a->TotalLength() > b->TotalLength();
287  }
288 
290  const art::Ptr<rb::Prong>& b) {
291  return a->TotalLength() > b->TotalLength();
292  }
293 
294  static bool sortShowerEnergy(const SRFuzzyKProng& a, const SRFuzzyKProng& b) {
295  return a.shwlid.shwE > b.shwlid.shwE;
296  }
297 
298  static bool sortProngCalE(const SRFuzzyKProng& a, const SRFuzzyKProng& b) {
299  return a.calE > b.calE;
300  }
301 
302  static bool sortElasticProng(const SRElastic& a, const SRElastic& b) {
303  return a.fuzzyk.png.size() > b.fuzzyk.png.size();
304  }
305 
306  static bool sortHoughProng(const SRHoughVertex& a, const SRHoughVertex& b) {
307  return a.fuzzyk.png.size() > b.fuzzyk.png.size();
308  }
309 
310  static bool sortMuIdLength(const SRMuId& a, const SRMuId& b) {
311  return a.len > b.len;
312  }
313 
314  static bool sortRemid(const SRKalmanTrack& a, const SRKalmanTrack& b) {
315  return a.rempid > b.rempid;
316  }
317 
319  const cheat::NeutrinoWithIndex& b) {
320  return a.truthColIndex < b.truthColIndex;
321  }
322 
323  std::pair<int, int> calcFirstLastLivePlane(int plane, std::bitset<14> binary,
325 
326  double SimpleOscProb(const simb::MCFlux& flux,
327  const simb::MCNeutrino& nu) const;
328 
332  art::Handle<std::vector<rawdata::RawTrigger>> trigs,
333  SRSpill& spill, art::Event& evt) const;
334 
335  enum class MixingType {
336  eRockOverlay,
340  eNone,
341  ePileup,
342  eSingles,
343  eOverlay,
346  };
347 
348  std::map<std::string, MixingType> mapStringMixingType = {
349  {"rock_overlay", MixingType::eRockOverlay},
350  {"cosmic_overlay", MixingType::eCosmicOverlay},
351  {"rock_cosmic_overlay", MixingType::eRockCosmicOverlay},
352  {"rock_singles_overlay", MixingType::eRockSinglesOverlay},
353  {"mcspill_singles_overlay", MixingType::eMCSpillSinglesOverlay},
354  {"dataspill_singles_overlay", MixingType::eDataSpillSinglesOverlay},
355  {"none", MixingType::eNone},
356  {"pileup", MixingType::ePileup},
357  {"singles", MixingType::eSingles},
358  {"overlay", MixingType::eOverlay}};
359 
360  std::map<MixingType, std::string> mapMixingTypeString = {
361  {MixingType::eRockOverlay, "rock_overlay"},
362  {MixingType::eCosmicOverlay, "cosmic_overlay"},
363  {MixingType::eRockCosmicOverlay, "rock_cosmic_overlay"},
364  {MixingType::eRockSinglesOverlay, "rock_singles_overlay"},
365  {MixingType::eMCSpillSinglesOverlay, "mcspill_singles_overlay"},
366  {MixingType::eDataSpillSinglesOverlay, "dataspill_singles_overlay"},
367  {MixingType::eNone, "none"},
368  {MixingType::ePileup, "pileup"},
369  {MixingType::eSingles, "singles"},
370  {MixingType::eOverlay, "overlay"}};
371 };
372 
373 //.......................................................................
374 
375 //.......................................................................
377  : fParams(params()), fIsRealData(false), fFile(0) {
378  // We update this, so have to pull it out of the config
380 
381  // Normally CAFMaker is run without an output ART stream, so these go
382  // nowhere, but can be occasionally useful for filtering as part of an
383  // ART job.
384  produces<std::vector<caf::StandardRecord>>();
385  produces<art::Assns<caf::StandardRecord, rb::Cluster>>();
386  produces<caf::SRSpill>();
387  if (fParams.FillNuTree()) produces<std::vector<caf::SRSpillTruthBranch>>();
388 }
389 
390 //......................................................................
392 
393 //......................................................................
395  if (!fFile) {
396  // Filename wasn't set in the FCL, and this is the
397  // first file we've seen
398  char temp[fb.fileName().size() + 1];
399  std::strcpy(temp, fb.fileName().c_str());
400  fCafFilename = basename(temp);
401  const size_t dotpos = fCafFilename.find('.');
402  assert(dotpos != std::string::npos); // Must have a dot, surely?
403  fCafFilename.resize(dotpos);
404  fCafFilename += fParams.FileExtension();
405 
407  }
408 }
409 
410 //......................................................................
412  if (!fCafFilename.empty()) InitializeOutfile();
413 }
414 
415 //......................................................................
417  fDetID = geom->DetId();
418  fDet = (Det_t)fDetID;
419 
420  if (!rh->GetDiBlockMaskFromCondb() && !fParams.ApplyingFilter()) { // get for each run; here if subrun
421  // masking not used and there is only
422  // one stream with calhit run
424  r.getByLabel(fParams.HitsLabel(), mask); // why doesn't fParams.HitInput() work here?
425  if (!mask.isValid()) {
426  mf::LogError("CAFMaker") << "Could not find a diBlock mask from " << fParams.HitsLabel() << ". Bailing.";
427  abort();
428  }
429 
430  fMask = *mask;
431  }
432 
433 
434  if (!rh->GetDiBlockMaskFromCondb() && fParams.ApplyingFilter()) { // get for each run; here if subrun
435  // masking not used and there are two
436  // streams one of which does not have
437  // calhits to save space and time
438  std::unique_ptr<int> mask(new int(rh->GoodDiBlockMask(0,true)));
439 
440  fMask = *mask;
441  }
442 }
443 
444 //......................................................................
446  if (rh->GetDiBlockMaskFromCondb() && !fParams.ApplyingFilter()) { // get if we are using subrun-level
447  // masks and there is only
448  // one stream with calhit run
450  sr.getByLabel(fParams.HitsLabel(), mask); // why doesn't fParams.HitInput() work here?
451  if (!mask.isValid()) {
452  mf::LogError("CAFMaker") << "Could not find a diBlock mask from " << fParams.HitsLabel() << ". Bailing.";
453  abort();
454  }
455 
456  fMask = *mask;
457  }
458 
459  if (rh->GetDiBlockMaskFromCondb() && fParams.ApplyingFilter()) { // get if we are using subrun-level
460  // masks and there are two
461  // streams one of which does not have
462  // calhits to save space and time
463  std::unique_ptr<int> mask(new int(rh->GoodDiBlockMask(sr.subRun(),false)));
464 
465  fMask = *mask;
466  }
467 }
468 
469 //......................................................................
471  assert(!fFile);
472  assert(!fCafFilename.empty());
473 
474  mf::LogInfo("CAFMaker") << "Output filename is " << fCafFilename;
475 
476  fFile = new TFile(fCafFilename.c_str(), "RECREATE");
477 
478  hPOT = new TH1D("TotalPOT", "TotalPOT;; POT", 1, 0, 1);
479  hSinglePOT =
480  new TH1D("TotalSinglePOT", "TotalSinglePOT;; Single POT", 1, 0, 1);
482  new TH1D("TotalTrueNonswapNue",
483  "TotalTrueNonswap;; Total True Nonswap Nue", 1, 0, 1);
484  hTotalTrueSingleNue = new TH1D(
485  "TotalTrueSingleNue", "TotalTrueSingle;; Total True Single Nue", 1, 0, 1);
486  hEvents = new TH1D("TotalEvents", "TotalEvents;; Events", 1, 0, 1);
487  hLivetime = new TH1D("TotalLivetime", "TotalLiveTime;;Seconds", 1, 0, 1);
488 
489  fRecTree = new TTree("recTree", "records");
490 
491  // Tell the tree it's expecting StandardRecord objects
492  StandardRecord* rec = 0;
493  fRecTree->Branch("rec", "caf::StandardRecord", &rec);
494 
495  // Create tree to store spill info for POT counting
496  fSpillTree = new TTree("spillTree", "spills");
497 
498  caf::SRSpill* dummy = 0;
499  fSpillTree->Branch("spill", "caf::SRSpill", &dummy);
500 
501  // Make the spillTree a friend of the fRecTree
502  fRecTree->AddFriend("spillTree");
503 
504 
505  if (fParams.FillNuTree()) {
506  fNuTree = new TTree("nuTree", "spillTruth");
507  caf::SRSpillTruthBranch* dummy2 = 0;
508  fNuTree->Branch("mc", "caf::SRSpillTruthBranch", &dummy2);
509  } else {
510  fNuTree = 0;
511  }
512 
513  fTotalPOT = 0;
514  fTotalSinglePOT = 0;
517  fTotalEvents = 0;
518  fTotalLivetime = 0;
519  fCycle = -5;
520  fBatch = -5;
521 
522  // Global information about the processing details:
523 
524  std::map<TString, TString> envmap;
525  std::string envstr;
526  // Environ comes from unistd.h
527  // environ is not present on OSX for some reason, so just use getenv to
528  // grab the variables we care about.
529 #ifdef DARWINBUILD
530  std::set<TString> variables;
531  variables.insert("USER");
532  variables.insert("HOSTNAME");
533  variables.insert("PWD");
534  variables.insert("SRT_PUBLIC_CONTEXT");
535  variables.insert("SRT_PRIVATE_CONTEXT");
536  for (auto var : variables) envmap[var] = getenv(var);
537 #else
538  for (char** penv = environ; *penv; ++penv) {
539  const std::string pair = *penv;
540  envstr += pair;
541  envstr += "\n";
542  const size_t split = pair.find("=");
543  if (split == std::string::npos) continue; // Huh?
544  const std::string key = pair.substr(0, split);
545  const std::string value = pair.substr(split + 1);
546  envmap[key] = value;
547  }
548 #endif
549 
550  // Get the command-line we were invoked with. What I'd really like is
551  // just the fcl script and list of input filenames in a less hacky
552  // fashion. I'm not sure that's actually possible in ART.
553  // TODO: ask the artists.
554  FILE* cmdline = fopen("/proc/self/cmdline", "rb");
555  char* arg = 0;
556  size_t size = 0;
558  while (getdelim(&arg, &size, 0, cmdline) != -1) {
559  cmd += arg;
560  cmd += " ";
561  }
562  free(arg);
563  fclose(cmdline);
564 
565  fFile->mkdir("env")->cd();
566 
567  TObjString(envmap["USER"]).Write("user");
568  TObjString(envmap["HOSTNAME"]).Write("hostname");
569  TObjString(envmap["PWD"]).Write("pwd");
570  TObjString(envmap["SRT_PUBLIC_CONTEXT"]).Write("publiccontext");
571  TObjString(envmap["SRT_PRIVATE_CONTEXT"]).Write("privatecontext");
572  // Default constructor is "now"
573  TObjString(TTimeStamp().AsString()).Write("date");
574  TObjString(cmd.c_str()).Write("cmd");
575  TObjString(fCafFilename.c_str()).Write("output");
576  TObjString(envstr.c_str()).Write("env");
577 }
578 
579 void CAFMaker::AddMetadataToFile(std::map<std::string, std::string> metadata) {
580  assert(fFile && "CAFMaker: Trying to add metadata to an uninitialized file");
581 
582  fFile->mkdir("metadata")->cd();
583 
584  for (auto const& pair : metadata) {
585  TObjString(pair.second.c_str()).Write(pair.first.c_str());
586  }
587 }
588 
589 //......................................................................
591  std::vector<cheat::TrackIDE>& allTracks,
592  std::vector<cheat::TrackIDE>& sliceTracks,
593  std::vector<cheat::TrackIDE>& allTracksBirks,
594  std::vector<cheat::TrackIDE>& sliceTracksBirks,
595  std::vector<SRNeutrino>* vec,
596  const art::Event& evt, const SRSpill& spill) {
597  // Get information about truNu
598  vec->push_back(SRNeutrino());
599  SRNeutrino& srTruth = vec->back();
600 
601  srTruth.run = evt.run();
602  srTruth.subrun = evt.subRun();
603  srTruth.evt = evt.id().event();
604  srTruth.cycle = fCycle;
605 
606  srTruth.isFHC = spill.isFHC;
607  srTruth.is0HC = spill.is0HC;
608  srTruth.isRHC = spill.isRHC;
609 
610  srTruth.det = fDet;
611 
612  const art::Ptr<simb::MCTruth>& tru = effPur->neutrinoInt;
613 
615 
616  // this is very rudimentary caching
617  static std::string oldVersionStr = std::string{};
618  static std::vector<unsigned int> oldVersionVec = std::vector<unsigned int>{};
619  static std::unordered_map<std::string, std::string> oldGenConfig = std::unordered_map<std::string, std::string>{};
620  static std::string oldGenConfigStr = std::string{};
621  if (tru->GeneratorInfo().generatorVersion != oldVersionStr || tru->GeneratorInfo().generatorConfig != oldGenConfig)
622  {
623  oldVersionStr = tru->GeneratorInfo().generatorVersion;
624  oldVersionVec = DecodeGeneratorVersion(tru->GeneratorInfo().generatorVersion);
625  oldGenConfig = tru->GeneratorInfo().generatorConfig;
626  oldGenConfigStr = "";
627  for (const auto & configPair: oldGenConfig)
628  {
629  if (configPair.first == "tune")
630  oldGenConfigStr = configPair.second;
631  else
632  mf::LogWarning("CAFMaker") << "Don't know how to store generator config parameter named '" << configPair.first << "'";
633  }
634  }
635  srTruth.genVersion = oldVersionVec;
636  if (!oldGenConfigStr.empty())
637  srTruth.genConfigString = oldGenConfigStr;
638 
639 
640  srTruth.pdg = tru->GetNeutrino().Nu().PdgCode();
641  srTruth.hitnuc = tru->GetNeutrino().HitNuc();
642 
643  srTruth.mode = tru->GetNeutrino().Mode(); // to avoid depending on simb
644  srTruth.iscc = (tru->GetNeutrino().CCNC() == simb::kCC);
645 
646  if (srTruth.mode == kElectronScattering) {
647  // For a numu or nutau scattering off an electron, the exchanged particle
648  // is a Z. But for a nue, there are two diagrams, one with a Z and one
649  // with a W. Both the final states have the electron ejected, so the
650  // interaction is a quantum-mechanical superposition between the
651  // two. It's neither really CC or NC. We do need to set this boolean
652  // though. Treating nu-on-e scattering as NC (ie not oscillating it) is
653  // wrong, the cross-section for nue is different to numu and nutau so
654  // oscillations do have an effect on the final spectrum. If we mark it CC
655  // the different oscillated fluxes will be weighted together in the fit
656  // appropriately, so this is safe, even if not perfectly physically
657  // accurate.
658  srTruth.iscc = true;
659  }
660  srTruth.inttype = tru->GetNeutrino().InteractionType();
661 
662  srTruth.q2 = tru->GetNeutrino().QSqr();
663  srTruth.x = tru->GetNeutrino().X();
664  srTruth.y = tru->GetNeutrino().Y();
665  srTruth.W2 = tru->GetNeutrino().W() * tru->GetNeutrino().W();
666  srTruth.E = tru->GetNeutrino().Nu().Momentum().E();
667  srTruth.visE = effPur->energyTotal;
668  srTruth.visEinslc = effPur->energySlice;
669  srTruth.eff = effPur->efficiency;
670  srTruth.pur = effPur->purity;
671  srTruth.nhitslc = effPur->nSliceHits;
672  srTruth.nhittot = effPur->nTotalHits;
673  srTruth.time = tru->GetNeutrino().Nu().Position().T();
674  srTruth.vtx = tru->GetNeutrino().Nu().Position().Vect();
675  srTruth.p = tru->GetNeutrino().Nu().Momentum();
676 
677  float visEBirks = 0;
678  for (std::vector<cheat::TrackIDE>::iterator it = allTracksBirks.begin();
679  it != allTracksBirks.end(); ++it) {
680  if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
681  visEBirks += it->energy;
682  }
683  float visEinslcBirks = 0;
684  for (std::vector<cheat::TrackIDE>::iterator it = sliceTracksBirks.begin();
685  it != sliceTracksBirks.end(); ++it) {
686  if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
687  visEinslcBirks += it->energy;
688  }
689 
690  srTruth.visEBirks = visEBirks;
691  srTruth.visEinslcBirks = visEinslcBirks;
692 
693  // Get tgtZ / tgtA from the MCTruth - which stores pdg code of nucleus
694  int tarpdg = tru->GetNeutrino().Target();
695  // pdgcode of nucleus = 10LZZZAAAI
696  srTruth.tgtZ = (tarpdg / 10000) % 1000;
697  srTruth.tgtA = (tarpdg / 10) % 1000;
698 
699  srTruth.genweight = tru->GetNeutrino().Nu().Weight();
700 
701  TVector3 nuVtx = tru->GetNeutrino().Nu().Position().Vect();
702  bool inDet = fabs(nuVtx.X()) <= geom->DetHalfWidth() &&
703  fabs(nuVtx.Y()) <= geom->DetHalfHeight() && nuVtx.Z() >= 0 &&
704  nuVtx.Z() <= geom->DetLength();
705  srTruth.isvtxcont = inDet;
706 
707  // rescue the parts of GTruth that we always want, put them in SRNeutrino
709  pv.push_back(tru);
711  FindManyPStrict<simb::GTruth>(pv, evt, fParams.GeneratorLabel());
712 
713  if (fmGtruth.isValid()) {
714  std::vector<art::Ptr<simb::GTruth>> gtruths = fmGtruth.at(0);
715 
716  if (!gtruths.empty()) {
717  const simb::GTruth& gtruth = *gtruths[0];
718 
719  // Useful for identifying pre-FSI generated hadron content
720  srTruth.npiplus = gtruth.fNumPiPlus;
721  srTruth.npiminus = gtruth.fNumPiMinus;
722  srTruth.npizero = gtruth.fNumPi0;
723  srTruth.nproton = gtruth.fNumProton;
724  srTruth.nneutron = gtruth.fNumNeutron;
725 
726  // Descriptive variables relating to the thrown cross sections
727  srTruth.ischarm = gtruth.fIsCharm;
728  srTruth.isseaquark = gtruth.fIsSeaQuark;
729  srTruth.resnum = gtruth.fResNum;
730  srTruth.xsec = gtruth.fXsec;
731  }
732  }
733 
734  // now stuff from the flux record
735  // art::FindManyP<simb::MCFlux> fmFlux(pv, evt, fParams.GeneratorLabel());
737  FindManyPStrict<simb::MCFlux>(pv, evt, fParams.GeneratorLabel());
738  if (fmFlux.isValid()) {
739  std::vector<art::Ptr<simb::MCFlux>> fluxes = fmFlux.at(0);
740 
741  if (!fluxes.empty()) {
742  const simb::MCFlux& flux = *fluxes[0];
743 
744  srTruth.pdgorig = flux.fntype;
745  srTruth.woscdumb = SimpleOscProb(flux, tru->GetNeutrino());
746 
747  SRBeam& beam = srTruth.beam;
748 
749  beam.tv = TVector3(flux.ftvx, flux.ftvy, flux.ftvz);
750  beam.tp = TVector3(flux.ftpx, flux.ftpy, flux.ftpz);
751  beam.runjob = flux.frun;
752  beam.potnum = flux.fevtno;
753  beam.tptype = flux.ftptype;
754  beam.nimpwt = flux.fnimpwt;
755  beam.ndecay = flux.fndecay;
756  beam.v = TVector3(flux.fvx, flux.fvy, flux.fvz);
757  beam.pdp = TVector3(flux.fpdpx, flux.fpdpy, flux.fpdpz);
758  beam.ppdxdz = flux.fppdxdz;
759  beam.ppdydz = flux.fppdydz;
760  beam.pppz = flux.fpppz;
761  beam.ppenergy = flux.fppenergy;
762  beam.ppmedium = flux.fppmedium;
763  beam.ptype = flux.fptype;
764  beam.ppv = TVector3(flux.fppvx, flux.fppvy, flux.fppvz);
765  beam.muparp = TVector3(flux.fmuparpx, flux.fmuparpy, flux.fmuparpz);
766  beam.mupare = flux.fmupare;
767  beam.necm = flux.fnecm;
768  beam.tgen = flux.ftgen;
769  beam.tgptype = flux.ftgptype;
770  beam.dk2gen = flux.fdk2gen;
771  beam.gen2vtx = flux.fgen2vtx;
772  beam.dk2vtx = flux.fdk2gen + flux.fgen2vtx;
773 
774  } else {
775  // this does not work for data-singlenu-mc overlay, but we don't want the
776  // singlenu flux anyway
777  if (!fParams.DataSpillHasMC())
778  throw cet::exception(
779  "Missing Association Error",
780  "Failed to find MCFlux object assoicated with MCTruth. This "
781  "probably means something went wrong with overlays.");
782  }
783  }
784 
785  // Must be after the setting of 'vtx' and 'beam' info
786  srTruth.L = TrueNeutrinoDistance(fDetID, srTruth);
787 
788  // Get the neutrino's daughters (and other progeny)
789  std::vector<const sim::Particle*> particles = bt->MCTruthToParticles(tru);
790 
791  // Get the sim::TrueEnergy's out of the event.
794  std::vector<sim::TrueEnergy> TrueEnergies;
795  if (!tes.failedToGet()) TrueEnergies = *tes;
796 
797  // We're going to loop over daughters to add other primaries to truNu.
798  // First we're gonna get the lepton and add it as the zeroth entry,
799  // so let's get its information.
800  int lepPdg = tru->GetNeutrino().Lepton().PdgCode();
801  double lepEnergy = tru->GetNeutrino().Lepton().E();
802  int lepCount = 0; // count the found leptons, for debugging.
803  int lepTrackId = 0; // We'll set this later
804  // Now we'll loop over the lepton and add it
805  float visENeutron = 0;
806  for (std::vector<cheat::TrackIDE>::iterator it = allTracks.begin();
807  it != allTracks.end(); ++it) {
808  if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
809  int mId = bt->TrackIDToParticle(it->trackID)->Mother();
810  while (mId != 0) {
811  if (bt->TrackIDToParticle(mId)->PdgCode() == 2112) {
812  visENeutron += it->energy;
813  break;
814  }
815  mId = bt->TrackIDToParticle(mId)->Mother();
816  }
817  }
818  float visEinslcNeutron = 0;
819  for (std::vector<cheat::TrackIDE>::iterator it = sliceTracks.begin();
820  it != sliceTracks.end(); ++it) {
821  if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
822  int mId = bt->TrackIDToParticle(it->trackID)->Mother();
823  while (mId != 0) {
824  if (bt->TrackIDToParticle(mId)->PdgCode() == 2112) {
825  visEinslcNeutron += it->energy;
826  break;
827  }
828  mId = bt->TrackIDToParticle(mId)->Mother();
829  }
830  }
831  float visENeutronBirks = 0;
832  for (std::vector<cheat::TrackIDE>::iterator it = allTracksBirks.begin();
833  it != allTracksBirks.end(); ++it) {
834  if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
835  int mId = bt->TrackIDToParticle(it->trackID)->Mother();
836  while (mId != 0) {
837  if (bt->TrackIDToParticle(mId)->PdgCode() == 2112) {
838  visENeutronBirks += it->energy;
839  break;
840  }
841  mId = bt->TrackIDToParticle(mId)->Mother();
842  }
843  }
844  float visEinslcNeutronBirks = 0;
845  for (std::vector<cheat::TrackIDE>::iterator it = sliceTracksBirks.begin();
846  it != sliceTracksBirks.end(); ++it) {
847  if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
848  int mId = bt->TrackIDToParticle(it->trackID)->Mother();
849  while (mId != 0) {
850  if (bt->TrackIDToParticle(mId)->PdgCode() == 2112) {
851  visEinslcNeutronBirks += it->energy;
852  break;
853  }
854  mId = bt->TrackIDToParticle(mId)->Mother();
855  }
856  }
857  srTruth.visENeutron = visENeutron;
858  srTruth.visEinslcNeutron = visEinslcNeutron;
859  srTruth.visENeutronBirks = visENeutronBirks;
860  srTruth.visEinslcNeutronBirks = visEinslcNeutronBirks;
861 
862  for (const sim::Particle* part : particles) {
863  // Add particle if it's a primary and it has the correct PDG code
864  // and the correct energy.
865  if (part->Mother() == 0 && part->PdgCode() == lepPdg &&
866  EssentiallyEqual(part->E(), lepEnergy, .01)) {
867  lepTrackId = part->TrackId();
868  ++lepCount;
869  AddParticleToVec(*part, allTracks, sliceTracks, allTracksBirks,
870  sliceTracksBirks, &srTruth.prim, TrueEnergies);
871  }
872  }
873 
874  if (lepCount != 1) {
875  switch (lepCount) {
876  case 0:
877  if (!srTruth.iscc) break;
878 
879  mf::LogError("CAFMaker")
880  << "*****************************************"
881  << "\n"
882  << "Uh, oh, CAFMaker did not find primary lepton. "
883  << "\n"
884  << "We're going to need to find a way to deal with this."
885  << "\n"
886  << "Info about the particle: "
887  << "\n"
888  << " pdg: " << srTruth.pdg << "\n"
889  << " pdgorig: " << srTruth.pdgorig << "\n"
890  << " E: " << srTruth.E << "\n"
891  << " visE: " << srTruth.visE << "\n"
892  << " nhittot: " << srTruth.nhittot << "\n"
893  << " X: " << srTruth.vtx.X() << "\n"
894  << " Y: " << srTruth.vtx.Y() << "\n"
895  << " Z: " << srTruth.vtx.Z() << "\n"
896  << " iscc: " << srTruth.iscc << "\n"
897  << " infid(nd): "
898  << (srTruth.vtx.X() < 200 && srTruth.vtx.X() > -200 &&
899  srTruth.vtx.Y() < 200 && srTruth.vtx.Y() > -200 &&
900  srTruth.vtx.Z() > 0 && srTruth.vtx.X() < 1600)
901  << "\n"
902 
903  << "*****************************************"
904  << "\n";
905  break;
906  default:
907  mf::LogError("CAFMaker")
908  << "*****************************************"
909  << "\n"
910  << "Uh, oh, CAFMaker found too many primary leptons in event. "
911  << "\n"
912  << "We're going to need to find a way to deal with this."
913  << "\n"
914  << "*****************************************";
915  break;
916  }
917  // assert(lepCount == 1);
918  }
919 
920  for (const sim::Particle* part : particles) {
921  // Add particle if it's a primary and it's not the lepton
922  if (part->TrackId() != lepTrackId && part->Mother() == 0)
923  AddParticleToVec(*part, allTracks, sliceTracks, allTracksBirks,
924  sliceTracksBirks, &srTruth.prim, TrueEnergies);
925  }
926 
927  FindAndAddMichels(particles, allTracks, &srTruth.michel);
928 }
929 
930 //......................................................................
931 /// This function is here for debugging purposes.
932 /// It will draw a table of with slices as rows and neutrinos as columns.
933 /// The favorites will be written with a star next to them.
934 /// If two faveIds vectors are provided, the second will have
935 /// an @ drawn next to them
937  const std::vector<std::vector<cheat::NeutrinoEffPur>>& sliceTruthTable,
938  const std::vector<int>& faveIds) {
939  std::vector<int> otherFaveIds(faveIds.size(), -1);
940  drawSliceTruthTable(sliceTruthTable, faveIds, otherFaveIds);
941 }
942 
943 //......................................................................
944 
945 /// This function is here for debugging purposes.
946 /// It will draw a table of with slices as rows and neutrinos as columns.
947 /// The favorites will be written with a star next to them.
948 /// If two faveIds vectors are provided, the second will have
949 /// an @ drawn next to them
951  const std::vector<std::vector<cheat::NeutrinoEffPur>>& sliceTruthTable,
952  const std::vector<int>& faveIds, const std::vector<int>& otherFaveIds) {
953  std::cout.precision(2);
954  std::cout.width(2);
955  std::cout.setf(std::ios::fixed, std::ios::floatfield);
956  int iSlc = 0;
957  for (auto slc : sliceTruthTable) {
958  int iNu = 0;
959  for (auto nu = slc.begin(); nu != slc.end(); ++nu, ++iNu) {
960  std::cout << nu->efficiency;
961  if (faveIds[iSlc] == iNu)
962  std::cout << "*";
963  else
964  std::cout << " ";
965  if (otherFaveIds[iSlc] == iNu)
966  std::cout << "@";
967  else
968  std::cout << " ";
969  std::cout << " ";
970  }
971  std::cout << std::endl;
972  }
973 }
974 
975 //......................................................................
976 template <class T, class U>
978  const art::Event& evt,
979  const art::InputTag& tag) const {
980  art::FindManyP<T> ret(from, evt, tag);
981 
982  if (!tag.label().empty() && !ret.isValid() && fParams.StrictMode()) {
983  std::cout << "CAFMaker: No Assn from '"
984  << abi::__cxa_demangle(typeid(from).name(), 0, 0, 0) << "' to '"
985  << abi::__cxa_demangle(typeid(T).name(), 0, 0, 0)
986  << "' found under label '" << tag << "'. "
987  << "Set 'StrictMode: false' to continue anyway." << std::endl;
988  abort();
989  }
990 
991  return ret;
992 }
993 
994 //......................................................................
995 template <class T>
997  T& ret) const {
998  if (!fm.isValid()) return false;
999 
1000  const std::vector<art::Ptr<T>> prods = fm.at(idx);
1001 
1002  if (prods.empty()) return false;
1003 
1004  ret = *prods[0];
1005 
1006  return true;
1007 }
1008 
1009 //......................................................................
1010 template <class T>
1012  art::Handle<T>& handle) const {
1013  evt.getByLabel(label, handle);
1014  if (!label.empty() && handle.failedToGet() && fParams.StrictMode()) {
1015  std::cout << "CAFMaker: No product of type '"
1016  << abi::__cxa_demangle(typeid(*handle).name(), 0, 0, 0)
1017  << "' found under label '" << label << "'. "
1018  << "Set 'StrictMode: false' to continue anyway." << std::endl;
1019  abort();
1020  }
1021 }
1022 
1023 //......................................................................
1024 template <class T>
1026  const std::string& label,
1027  art::Handle<T>& handle) const {
1028  evt.getByLabel(label, handle);
1029  if (!label.empty() && handle.failedToGet() && fParams.StrictMode()) {
1030  std::cout << "CAFMaker: No product of type '"
1031  << abi::__cxa_demangle(typeid(*handle).name(), 0, 0, 0)
1032  << "' found under label '" << label << "'. "
1033  << "Continuing without it." << std::endl;
1034  }
1035 }
1036 
1037 //......................................................................
1038 template <class T>
1040  const std::vector<std::string>& name,
1041  T& ret) const {
1042  fhicl::ParameterSet p = pset;
1043  for (unsigned int i = 0; i < name.size() - 1; ++i) {
1044  if (!p.has_key(name[i])) return false;
1045  p = p.get<fhicl::ParameterSet>(name[i]);
1046  }
1047  if (!p.has_key(name.back())) return false;
1048  ret = p.get<T>(name.back());
1049  return true;
1050 }
1051 
1052 //......................................................................
1053 std::pair<int, int> CAFMaker::calcFirstLastLivePlane(int plane,
1054  std::bitset<14> binary,
1055  caf::Det_t det) {
1056  // if not the FD, return full ND size
1057  if (det != caf::kFARDET) return std::pair<int, int>(0, 213);
1058  int testDB = (plane / 64);
1059  int minblock = testDB;
1060  int maxblock = testDB;
1061  // find block boundaries;
1062  for (int i = testDB - 1; i >= 0; --i) {
1063  if (binary[i])
1064  minblock = i;
1065  else
1066  break;
1067  }
1068  for (int i = testDB + 1; i < 14; ++i) {
1069  if (binary[i])
1070  maxblock = i;
1071  else
1072  break;
1073  }
1074  return std::pair<int, int>(64 * (minblock), 64 * (maxblock + 1) - 1);
1075 }
1076 
1077 //......................................................................
1079  // Normally CAFMaker is run without an output ART stream, so these go
1080  // nowhere, but can be occasionally useful for filtering as part of
1081  // an ART job.
1082  std::unique_ptr<std::vector<caf::StandardRecord>> srcol(
1083  new std::vector<caf::StandardRecord>);
1084  std::unique_ptr<art::Assns<caf::StandardRecord, rb::Cluster>> srAssn(
1086  std::unique_ptr<caf::SRSpill> spillcol(new SRSpill);
1087  std::unique_ptr<std::vector<caf::SRSpillTruthBranch>> spilltruthcol(
1088  new std::vector<SRSpillTruthBranch>);
1089 
1090  fTotalEvents += 1;
1091  // Adjust for art misfeature isRealData confusing data-mc overlay
1092  fIsRealData = evt.isRealData() || fParams.DataSpillHasMC();
1093  // -- Case -- : Data? == isRealData || DataSpillHasMC :
1094  // bt->HaveTruthInfo
1095  // Normal Data : true true false : false
1096  // Data Overlaid on MC : true false true : true
1097  // Just MC : false false false : true
1098 
1099  // Get a handle on the slices
1100 
1103 
1104  /* // Define a vector holding slices
1105  art::PtrVector<rb::Cluster> subevts;
1106  for(unsigned int i = 0; i < slices->size(); ++i)
1107  subevts.push_back(art::Ptr<rb::Cluster>(slices, i));
1108  art::Handle<std::vector<rb::Cluster> > hparentSlices;
1109  art::PtrVector<rb::Cluster> parentSubevts;*/
1110 
1111  if (fParams.IsSinglesOverlay()) {
1113  art::Handle<std::vector<simb::MCTruth>> singlemctruthcol;
1114  GetByLabelStrict(evt, fParams.GeneratorLabel(), mctruthcol);
1115  GetByLabelStrict(evt, fParams.SingleMixerLabel(), singlemctruthcol);
1116 
1117  if (mctruthcol->empty()) {
1118  mf::LogWarning("CAFMaker") << "Error retrieving MCTruth list";
1119  return;
1120  }
1121  if (singlemctruthcol->empty()) {
1122  mf::LogWarning("CAFMaker")
1123  << "Error retrieving MCTruth list from singles";
1124  return;
1125  }
1126 
1127  for (size_t i_intx = 0; i_intx < mctruthcol->size(); ++i_intx) {
1128  LOG_DEBUG("outnNueCCOverlay") << "start loop";
1129  simb::MCTruth const& truth = mctruthcol->at(i_intx);
1130  simb::MCNeutrino const& mc_neutrino = truth.GetNeutrino();
1131  simb::MCParticle const& mc_particle = mc_neutrino.Nu();
1132 
1133  if (mc_neutrino.CCNC() == 0 && fabs(mc_particle.PdgCode()) == 12) {
1135  }
1136  LOG_DEBUG("CAFMaker") << "end loop over mctruth";
1137  }
1138 
1139  for (size_t i_intx = 0; i_intx < singlemctruthcol->size(); ++i_intx) {
1140  LOG_DEBUG("outnNueCCOverlay") << "start loop";
1141  simb::MCTruth const& truth = singlemctruthcol->at(i_intx);
1142  simb::MCNeutrino const& mc_neutrino = truth.GetNeutrino();
1143  simb::MCParticle const& mc_particle = mc_neutrino.Nu();
1144 
1145  if (mc_neutrino.CCNC() == 0 && fabs(mc_particle.PdgCode()) == 12) {
1147  }
1148  LOG_DEBUG("CAFMaker") << "end loop over mctruth";
1149  }
1150  }
1151 
1153 
1154  // in case we are running over MRE, make a vector of
1155  // parent slices, so we can pass it to backtracker
1156  // for filling the truth branch
1157 
1158  if (fParams.MRMode()) {
1159  GetByLabelStrict(evt, fParams.MRParentSliceLabel(), parentSlices);
1160  }
1161 
1162  // Get all of the hits from calhits
1163 
1165  GetByLabelStrict(evt, fParams.HitsLabel(), allHits);
1166 
1167  // We're gonna get a vector of NeutrinoEffPurs for allHits
1168  std::vector<cheat::NeutrinoWithIndex> allNus = bt->allMCTruth();
1169  std::vector<std::vector<cheat::NeutrinoEffPur>> sliceTruthTable;
1170 
1171 
1172  // Try to get POT spill info if this is real data
1174 
1175  // treat as data spill in MC-dataspill overlay
1176  if (!fIsRealData) {
1178  } else {
1179  GetByLabelStrict(evt, fParams.NuMIBeamLabel(), spillPot);
1180  }
1181 
1182  // Get the EventQuality information
1184  GetByLabelStrict(evt, fParams.SpillQualLabel(), spillQual);
1185 
1186  // Get DAQ Header information
1189 
1192 
1194  FillSpillVars(spillPot, spillQual, daq, trigs, spill, evt);
1195 
1196  // Fill the spill tree once per event so that POT counting can be
1197  // done properly.
1198  // This tree will have an entry even for empty events with no physics
1199  // slices.
1200  caf::SRSpill* pspill = &spill;
1201  fSpillTree->SetBranchAddress("spill", &pspill);
1202 
1204  if (fParams.FillCosmicCVN()){
1205  GetByLabelStrict(evt, fParams.CosmicCVNLabel(), cosmiccvn);
1206  std::vector<cvntf::CVNCosmicFiltList> cosmicList = *cosmiccvn;
1207  for (unsigned int iC = 0; iC < cosmicList[0].ListSize(); ++iC){
1208  cvntf::CVNCosmicFilt cvn = cosmicList[0].GetTimeSlice(iC);
1209  caf::SRCosmicCVN tcosmic;
1210  tcosmic.nHits = cvn.nHits;
1211  tcosmic.timeWinMax = cvn.timeWinMax;
1212  tcosmic.timeWinMin = cvn.timeWinMin;
1213  tcosmic.numuVal = cvn.numuVal;
1214  tcosmic.nueVal = cvn.nueVal;
1215  tcosmic.nutauVal = cvn.nutauVal;
1216  tcosmic.ncVal = cvn.ncVal;
1217  tcosmic.cosmicVal = cvn.cosmicVal;
1218  tcosmic.passSel = cvn.passSel;
1219  spill.cosmiccvn.push_back(tcosmic);
1220  }
1221  spill.ncosmiccvn = spill.cosmiccvn.size();
1222  }
1223  fSpillTree->Fill();
1224 
1225  // Pull the MC cycle number from metadata
1226  // before filling nuTree so cycle is filled
1227  // Assumes metadata module is being run...
1228  std::map<std::string, std::string> metadata =
1230  // Check it exists in the metadata, otherwise it defaults to -5
1231  if (metadata.count("simulated.cycle"))
1232  fCycle = std::stoi(metadata["simulated.cycle"]);
1233  if (metadata.count("simulated.batch"))
1234  fBatch = std::stoi(metadata["simulated.batch"]);
1235 
1236  // Fill the nuTree one per event, in sync with the spillTree
1237  if (fParams.FillNuTree()) {
1238  caf::SRSpillTruthBranch spillTruth;
1239  caf::SRSpillTruthBranch* pspillTruth = &spillTruth;
1240  fNuTree->SetBranchAddress("mc", &pspillTruth);
1241 
1242  for (const cheat::NeutrinoWithIndex& nuIdx :
1243  bt->allMCTruth()) { // NeutrinoInteractions()){
1244  // AddMCTruthToVec wants one of these, but most of this info is only
1245  // relevant if the truth is associated with a slice.
1246  const cheat::NeutrinoEffPur ep = {nuIdx.neutrinoInt, 0, 0, 0, 0, 0, 0, 0};
1247  std::vector<cheat::TrackIDE> dummy;
1248  std::vector<SRNeutrino> nu;
1249  AddMCTruthToVec(&ep, dummy, dummy, dummy, dummy, &nu, evt, spill);
1250  spillTruth.nu = nu[0];
1251 
1252  art::Ptr<simb::MCTruth> mctruth = nuIdx.neutrinoInt;
1253  const std::vector<art::Ptr<simb::MCTruth>> mctv(1, mctruth);
1254 
1256  FindManyPStrict<fxwgt::FluxWeights>(mctv, evt,
1258  fxwgt::FluxWeights flxrw;
1259  if (GetAssociatedProduct(fmFLXRW, 0, flxrw))
1260  spillTruth.nu.rwgt.ppfx = FluxReweights(flxrw);
1261  else
1262  spillTruth.nu.rwgt.ppfx.setDefault();
1263 
1265  FindManyPStrict<rwgt::GENIEReweightTable>(mctv, evt,
1268  if (GetAssociatedProduct(fmRW, 0, rw))
1269  spillTruth.nu.rwgt.genie = GenieReweightTable(rw);
1270 
1271  fNuTree->Fill();
1272  spilltruthcol->push_back(spillTruth);
1273  } // end for nuIdx
1274  } // end if FillNuTree
1275 
1276  // Applying an event filter but we need the unfiltered events for POT/LT accounting.
1277  // Thus we need the same version of cafmaker to run over the filtered events and !filtered.
1278  // However, the point is to save time so we don't want slicer or calhit on the !filtered stream.
1279  if(fParams.ApplyingFilter() && !slices.isValid()) {
1280  *spillcol = spill;
1281 
1282  evt.put(std::move(srcol));
1283  evt.put(std::move(spillcol));
1284  if (fParams.FillNuTree()) evt.put(std::move(spilltruthcol));
1285  evt.put(std::move(srAssn));
1286 
1287  return;
1288  }
1289 
1291  FindManyPStrict<murem::MRCCParent>(slices, evt,
1293  if (fParams.MRMode()) {
1294  std::vector<std::vector<cheat::NeutrinoEffPur>> temptable // I know :)
1295  = bt->SlicesToMCTruthsTable(*parentSlices);
1296  // some mre slices may not have matched parents
1297  // so sort the table so that nth entry in table
1298  // is the truth for the parent of the nth mre
1299  // slice. If nth mre slice has no parent,
1300  // make an empty effpur entry.
1301  for (unsigned int isli = 0; isli < slices->size(); isli++) {
1303  if (GetAssociatedProduct(fmMrccPar, isli, parent)) {
1304  int idx = parent.ParentIndex();
1305  sliceTruthTable.push_back(temptable[idx]);
1306  } else {
1308  cheat::NeutrinoEffPur effPur =
1309  cheat::NeutrinoEffPur{null, -1, -1, 0, -1, -1, -1, -1};
1310  std::vector<cheat::NeutrinoEffPur> emptyEffPur(allNus.size(), effPur);
1311  sliceTruthTable.push_back(emptyEffPur);
1312  }
1313  } // end loop over mr slices
1314  } // end of mr mode
1315  else
1316  sliceTruthTable = bt->SlicesToMCTruthsTable(*slices);
1317 
1318  // Find favorite neutrino IDs using several metrics. Total energy
1319  // is the one we will use, for the favorite neutrino, but the others
1320  // will be stored.
1321  std::vector<int> faveIdsEff = bt->SliceToOrderedNuIds(
1322  allNus, sliceTruthTable, cheat::EffMetric, cheat::EffMetric);
1323  std::vector<int> faveIdsEnergy = bt->SliceToOrderedNuIds(
1324  allNus, sliceTruthTable, cheat::EnergyMetric, cheat::EnergyMetric);
1325  std::vector<int> faveIdsPur = bt->SliceToOrderedNuIds(
1326  allNus, sliceTruthTable, cheat::PurMetric, cheat::PurMetric);
1327  std::vector<int> faveIdsEffPur = bt->SliceToOrderedNuIds(
1328  allNus, sliceTruthTable, cheat::EffPurMetric, cheat::EffPurMetric);
1329  std::vector<int> faveIdsEffThenPur = bt->SliceToOrderedNuIds(
1330  allNus, sliceTruthTable, cheat::EffMetric, cheat::PurMetric);
1331 
1332  // Get all truth track hits in the event
1333  std::vector<cheat::TrackIDE> allTracks;
1334  std::vector<cheat::TrackIDE> sliceTracks;
1335  std::vector<cheat::TrackIDE> allTracksBirks;
1336  std::vector<cheat::TrackIDE> sliceTracksBirks;
1337  if (bt->HaveTruthInfo()) allTracks = bt->HitsToTrackIDE(*allHits);
1338  if (bt->HaveTruthInfo()) allTracksBirks = bt->HitsToTrackIDE(*allHits, true);
1339 
1340  // Get all association finder objects
1341  auto fmDiscrete2d =
1342  FindManyPStrict<rb::Track>(slices, evt, fParams.DiscreteTrack2dLabel());
1343  auto fmDiscrete =
1344  FindManyPStrict<rb::Track>(slices, evt, fParams.DiscreteTrackLabel());
1345  auto fmKalman =
1346  FindManyPStrict<rb::Track>(slices, evt, fParams.KalmanTrackLabel());
1347  auto fmCosmic =
1348  FindManyPStrict<rb::Track>(slices, evt, fParams.CosmicTrackLabel());
1349  auto fmWindow =
1350  FindManyPStrict<rb::Track>(slices, evt, fParams.WindowTrackLabel());
1351  auto fmElastic =
1352  FindManyPStrict<rb::Vertex>(slices, evt, fParams.ElasticArmsLabel());
1353  auto fmHoughVertex =
1354  FindManyPStrict<rb::Vertex>(slices, evt, fParams.HoughVertexLabel());
1355  auto fmVertexDT =
1356  FindManyPStrict<rb::Vertex>(slices, evt, fParams.VertexDTLabel());
1357  auto fmSlcME = FindManyPStrict<me::SlcME>(slices, evt, fParams.MELabel());
1358  auto fmJmShower =
1359  FindManyPStrict<jmshower::JMShower>(slices, evt, fParams.JMShowerLabel());
1360  auto fmNueSel =
1361  FindManyPStrict<jmshower::EID>(slices, evt, fParams.NueSelLabel());
1363  slices, evt, fParams.LIDLabel()); // one lid object per slice
1365  slices, evt, fParams.SlcLIDLabel()); // one lid object per slice
1366  art::FindMany<slid::ShowerLID> fmShwLID(slices, evt, fParams.LIDLabel());
1367  art::FindMany<slid::ShowerPID> fmShwPID(slices, evt, fParams.SPIDLabel());
1368  auto fmSliceLID = FindManyPStrict<SliceLID::Prediction>(
1369  slices, evt, fParams.SliceLIDLabel());
1370  auto fmNumuEnergy =
1371  FindManyPStrict<numue::NumuE>(slices, evt, fParams.NumuEnergyLabel());
1372  auto fmNumuLSTMEnergy = FindManyPStrict<LSTME::LSTMEnergy>(
1373  slices, evt, fParams.NumuLSTMEnergyLabel());
1374  auto fmNuePresel = FindManyPStrict<presel::PreselObj>(
1375  slices, evt, fParams.NuePreselLabel());
1376  auto fmRockPresel = FindManyPStrict<presel::PreselObj>(
1377  slices, evt, fParams.RockPreselLabel());
1378  auto fmVeto = FindManyPStrict<presel::Veto>(slices, evt, fParams.VetoLabel());
1379  auto fmNueVeto =
1380  FindManyPStrict<presel::Veto>(slices, evt, fParams.NueVetoLabel());
1381  auto fmLem =
1382  FindManyPStrict<lem::PIDDetails>(slices, evt, fParams.LemLabel());
1383  auto fmRvp = FindManyPStrict<rvp::RVP>(slices, evt, fParams.RvpLabel());
1384  auto fmXnue = FindManyPStrict<xnue::Xnue>(slices, evt, fParams.XnueLabel());
1385  auto fmCosRej =
1386  FindManyPStrict<cosrej::CosRejObj>(slices, evt, fParams.CosRejLabel());
1387  auto fmNueCosRej =
1388  FindManyPStrict<cosrej::NueCosRej>(slices, evt, fParams.NueCosRejLabel());
1389  auto fmNCCosRej = FindManyPStrict<ncid::NCCosRej>(
1390  slices, evt, fParams.NCCosRejLabel()); /// NC Cosmic Rejection
1391  auto fmNCPi0BkgRej = FindManyPStrict<ncpi0::NCPi0BkgRej>(
1392  slices, evt, fParams.NCPi0BkgRejLabel()); /// NC Pi0 Bkg Rejection
1393 
1394  auto fmNuonE =
1395  FindManyPStrict<cvn::Result>(slices, evt, fParams.NuonELabel());
1396 
1397  auto fmCVN_LoosePreselPtP =
1398  FindManyPStrict<cvn::Result>(slices, evt, fParams.CVNLabelLoosePreselPtP());
1399  auto fmCVN_OldPresel =
1400  FindManyPStrict<cvn::Result>(slices, evt, fParams.CVNLabelOldPresel());
1401  auto fmCVN_NoCosmics =
1402  FindManyPStrict<cvn::Result>(slices, evt, fParams.CVNLabelNoCosmics());
1403 
1404  auto fmWS =
1405  FindManyPStrict<cvn::Result>(slices, evt, fParams.WSCNNLabel());
1406 
1407  auto fmCVNFeatures =
1408  FindManyPStrict<cvn::Features>(slices, evt, fParams.CVNLabelLoosePreselPtP());
1409  auto fmCVNPixelMaps =
1410  FindManyPStrict<cvn::PixelMap>(slices, evt, fParams.CVNPixelMapLabel());
1411  auto fmCVNTrainingData =
1412  FindManyPStrict<cvn::TrainingData> (slices, evt, fParams.CVNTrainingDataLabel());
1413  auto fmEvtRegCVN =
1414  FindManyPStrict<cvn::RegNuResult>(slices, evt, fParams.RegCVNLabel());
1415  auto fmHadRegCVN =
1416  FindManyPStrict<cvn::RegHadronResult>(slices, evt, fParams.RegCVNLabel());
1417 
1418 
1420  if (fParams.IsSinglesOverlay()) {
1421  GetByLabelStrict(evt, fParams.SingleMixerLabel(), singlepot);
1422  if (!singlepot.failedToGet()) {
1423  float totsinglepot = singlepot->spillpot;
1424  fTotalSinglePOT += totsinglepot;
1425  }
1426  }
1427 
1428  int nNonNoise = 0;
1429  for (const rb::Cluster& slice : *slices) {
1430  if (slice.IsNoise() || slice.NCell() == 0) continue;
1431  ++nNonNoise;
1432  }
1433  for (unsigned int sliceId = 0; sliceId < slices->size(); ++sliceId) {
1434  const rb::Cluster& slice = (*slices)[sliceId];
1435 
1436  if (slice.IsNoise() || slice.NCell() == 0) continue;
1437  // Because we don't care about the noise slice and slices with no hits.
1439  StandardRecord* prec = &rec; // TTree wants a pointer-to-pointer
1440  fRecTree->SetBranchAddress("rec", &prec);
1441 
1442  if (bt->HaveTruthInfo()) sliceTracks = bt->HitsToTrackIDE(slice.AllCells());
1443  if (bt->HaveTruthInfo())
1444  sliceTracksBirks = bt->HitsToTrackIDE(slice.AllCells(), true);
1445  //#######################################################
1446  // Fill header.
1447  //#######################################################
1448  // Get metadata information for header
1449  unsigned int run = evt.run();
1450  unsigned int subrun = evt.subRun();
1451  unsigned int spillNum = evt.id().event();
1452 
1453  art::Timestamp ts = evt.time();
1454  const time_t timeSec = ts.timeHigh();
1455 
1456  struct tm* timeStruct = localtime(&timeSec);
1457 
1458  unsigned int maskstatus = 0;
1459  unsigned int dibmask = fMask;
1460 
1461  if (rh->GoodDiBlockMask() == ((1 << 16) - 1))
1462  maskstatus = 0;
1463  else if (rh->GoodDiBlockMask() == ((1 << 17) - 1))
1464  maskstatus = 2;
1465  else
1466  maskstatus = 1;
1467 
1468  unsigned int dibfirst = 0;
1469  unsigned int diblast = 0;
1470  if (dibmask) {
1471  int iD;
1472  iD = 0;
1473  while (!((dibmask >> iD) & 1)) iD++;
1474  dibfirst = iD + 1;
1475  iD = 0;
1476  while (dibmask >> iD) iD++;
1477  diblast = iD;
1478  }
1479 
1480  unsigned int nbadchan = bc->NBadInSubRun(subrun);
1481  unsigned int ntotchan = rh->NAnalysisChannels();
1482 
1483  rec.hdr = SRHeader();
1484 
1485  rec.hdr.run = run;
1486  rec.hdr.subrun = subrun;
1487  rec.hdr.cycle = fCycle;
1488  rec.hdr.batch = fBatch;
1489  rec.hdr.evt = spillNum;
1490  rec.hdr.subevt = sliceId;
1491  rec.hdr.ismc = !fIsRealData;
1492  rec.hdr.det = fDet;
1493  rec.hdr.blind = 0;
1494  rec.hdr.filt = rb::IsFiltered(evt, slices, sliceId);
1495 
1496  rec.hdr.subevtstarttime = slice.MinTNS();
1497  rec.hdr.subevtendtime = slice.MaxTNS();
1498  rec.hdr.subevtmeantime = slice.MeanTNS();
1499  rec.hdr.unixtime = timeSec;
1500 
1501  rec.hdr.dibfirst = dibfirst;
1502  rec.hdr.diblast = diblast;
1503  rec.hdr.dibmask = dibmask;
1504  rec.hdr.maskstatus = maskstatus;
1505 
1506  rec.hdr.year = timeStruct->tm_year + 1900;
1507  rec.hdr.month = timeStruct->tm_mon + 1;
1508  rec.hdr.day = timeStruct->tm_mday;
1509  rec.hdr.doy = timeStruct->tm_yday + 1;
1510  rec.hdr.hour = timeStruct->tm_hour + 1;
1511  rec.hdr.minute = timeStruct->tm_min;
1512  rec.hdr.second = timeStruct->tm_sec;
1513 
1514  rec.hdr.nbadchan = nbadchan;
1515  rec.hdr.ntotchan = ntotchan;
1516 
1517  rec.hdr.gain = rh->DetGainSetting();
1519 
1520  //#######################################################
1521  // Fill Spill info.
1522  //#######################################################
1523  // Store info on NuMI beam quality and if the pot as good or bad.
1524  // This is meant to be used to refine a good spill, not for
1525  // POT counting since events have multiple slices and empty
1526  // spill are missing from this tree. Use the friend tree for
1527  // POT counting.
1528  rec.spill = spill;
1529 
1530  //trim list of Cosmic CVN scores down to only the time window the
1531  //physics slice matches
1532  for (int iC = 0; iC < (int)rec.spill.cosmiccvn.size(); ++iC){
1533 
1534  if ((slice.MeanTNS() > rec.spill.cosmiccvn[iC].timeWinMin) &&
1535  (slice.MeanTNS() < rec.spill.cosmiccvn[iC].timeWinMax)) continue;
1536  rec.spill.cosmiccvn.erase(rec.spill.cosmiccvn.begin()+iC);
1537  --iC;
1538 
1539  }
1540  //if more then one time window matches the slice, keep the one with the
1541  //lowest cosmic CVN score
1542  if (rec.spill.cosmiccvn.size() > 1){
1543  SRCosmicCVN tcos = rec.spill.cosmiccvn[0];
1544  for (unsigned int iC = 1; iC < rec.spill.cosmiccvn.size(); ++iC){
1545  if (rec.spill.cosmiccvn[iC].cosmicVal < tcos.cosmicVal)
1546  tcos = rec.spill.cosmiccvn[iC];
1547  }
1548  rec.spill.cosmiccvn.clear();
1549  rec.spill.cosmiccvn.push_back(tcos);
1550  }
1551  rec.spill.ncosmiccvn = rec.spill.cosmiccvn.size();
1552 
1553 
1554  //#######################################################
1555  // Add slice info.
1556  //#######################################################
1557  FillSliceVars(slice, rec.slc);
1558  rec.slc.nnonnoise = nNonNoise;
1559  // find the slice with minimum mean time distance from this slice
1560  float mintimediff = 1e7; // initialize to a big number
1561  signed int mintimejsli = -1;
1562  for (unsigned int jsli = 0; jsli < slices->size(); jsli++) {
1563  const rb::Cluster& slicej = (*slices)[jsli];
1564  if (jsli == sliceId || slicej.IsNoise()) continue;
1565  float diff = fabs(slicej.MeanTNS() - slice.MeanTNS());
1566  if (diff < mintimediff) {
1567  mintimejsli = jsli;
1568  }
1569  mintimediff = std::min(mintimediff, diff);
1570  }
1571  if (mintimejsli != -1) {
1572  const rb::Cluster& slicejmin = (*slices)[mintimejsli];
1573  rec.slc.closestslicetime = slicejmin.MeanTNS() - slice.MeanTNS();
1574  rec.slc.closestslicenhit = slicejmin.NCell();
1575  rec.slc.closestslicecalE = slicejmin.CalorimetricEnergy();
1576  TVector3 minp, maxp;
1577  minp = slicejmin.MinXYZ();
1578  maxp = slicejmin.MaxXYZ();
1585  float mincelldistX = 55555;
1586  float mincelldistY = 55555;
1587  float mincelldist = 55555;
1588  float cellD = 2. * geom->Plane(0)->Cell(0)->HalfD();
1589  float cellW = 2. * geom->Plane(0)->Cell(0)->HalfW();
1590  for (unsigned int hitIdxX = 0; hitIdxX < slice.NXCell(); ++hitIdxX) {
1591  const art::Ptr<rb::CellHit>& chitX = slice.XCell(hitIdxX);
1592  for (unsigned int minhitIdxX = 0; minhitIdxX < slicejmin.NXCell();
1593  ++minhitIdxX) {
1594  const art::Ptr<rb::CellHit>& minchitX = slicejmin.XCell(minhitIdxX);
1595  float dPlaneX = (float)chitX->Plane() - (float)minchitX->Plane();
1596  float dCellX = (float)chitX->Cell() - (float)minchitX->Cell();
1597  float distX = sqrt((dPlaneX * cellD) * (dPlaneX * cellD) +
1598  (dCellX * cellW) * (dCellX * cellW));
1599  if (distX < mincelldistX) mincelldistX = distX;
1600  }
1601  }
1602  for (unsigned int hitIdxY = 0; hitIdxY < slice.NYCell(); ++hitIdxY) {
1603  const art::Ptr<rb::CellHit>& chitY = slice.YCell(hitIdxY);
1604  for (unsigned int minhitIdxY = 0; minhitIdxY < slicejmin.NYCell();
1605  ++minhitIdxY) {
1606  const art::Ptr<rb::CellHit>& minchitY = slicejmin.YCell(minhitIdxY);
1607  float dPlaneY = (float)chitY->Plane() - (float)minchitY->Plane();
1608  float dCellY = (float)chitY->Cell() - (float)minchitY->Cell();
1609  float distY = sqrt((dPlaneY * cellD) * (dPlaneY * cellD) +
1610  (dCellY * cellW) * (dCellY * cellW));
1611  if (distY < mincelldistY) mincelldistY = distY;
1612  }
1613  }
1614  if (mincelldistX < mincelldistY) {
1615  mincelldist = mincelldistX;
1616  } else {
1617  mincelldist = mincelldistY;
1618  }
1619  rec.slc.closestslicemindist = mincelldist;
1620  } else { // protection when mintimejsli == -1
1621  rec.slc.closestslicetime = -5;
1622  rec.slc.closestslicenhit = 0;
1623  rec.slc.closestslicecalE = -5;
1624  rec.slc.closestsliceminfromtop = -5;
1625  rec.slc.closestsliceminfrombottom = -5;
1626  rec.slc.closestsliceminfromfront = -5;
1627  rec.slc.closestsliceminfromback = -5;
1628  rec.slc.closestsliceminfromeast = -5;
1629  rec.slc.closestsliceminfromwest = -5;
1630  rec.slc.closestslicemindist = 55555;
1631  }
1632  if (fDet == kFARDET) {
1633  rec.sel.contain.nplanestofront = rec.slc.firstplane - (dibfirst - 1) * 64;
1634  rec.sel.contain.nplanestoback = ((diblast)*64) - 1 - rec.slc.lastplane;
1635  }
1636 
1637  if (fDet == kNEARDET) {
1639  rec.sel.contain.nplanestoback = 214 - rec.slc.lastplane;
1640  }
1641 
1642  //#######################################################
1643  // Adding reconstructed objects.
1644  //#######################################################
1645 
1646  // Get all the tracks associated with this slice from discrete
1647  // tracker. This depends on the findMany object created above.
1648  if (fmDiscrete.isValid()) {
1649  std::vector<art::Ptr<rb::Track>> discTracks = fmDiscrete.at(sliceId);
1650  // Sort the tracks in descending length order
1651  std::sort(discTracks.begin(), discTracks.end(), sortRBTrackLength);
1652  rec.trk.discrete.tracks.resize(discTracks.size());
1653  for (size_t trackId = 0; trackId < discTracks.size(); ++trackId) {
1654  FillTrackVars((*discTracks[trackId]), rec.trk.discrete.tracks[trackId],
1655  *slices, allTracks, sliceTracks, allTracksBirks,
1656  sliceTracksBirks, sliceId);
1657  }
1658  }
1659 
1660  // if(fmDiscrete2d.isValid())
1661  // {
1662  // std::vector< art::Ptr<rb::Track> > discTracks =
1663  // fmDiscrete2d.at(sliceId);
1664  // Sort the tracks in descending length order
1665  // std::sort(discTracks.begin(), discTracks.end(),
1666  // sortRBTrackLength);
1667  // rec.trk.discrete.tracks2d.resize(discTracks.size());
1668  // for(size_t trackId = 0; trackId < discTracks.size(); ++trackId){
1669  // FillTrackVars(*discTracks[trackId],
1670  // rec.trk.discrete.tracks2d[trackId], *slices, sliceId);
1671  // }
1672  // }
1673 
1674  if (fmCosmic.isValid()) {
1675  std::vector<art::Ptr<rb::Track>> cosTracks = fmCosmic.at(sliceId);
1676  art::FindManyP<me::TrkME> fmCosmicTrkME(cosTracks, evt, fParams.MELabel());
1677 
1678  // Sort the tracks in descending length order
1679  std::sort(cosTracks.begin(), cosTracks.end(), sortRBTrackLength);
1680 
1681  rec.trk.cosmic.tracks.reserve(cosTracks.size());
1682  rec.trk.cosmic.tracks.resize(cosTracks.size());
1683  for (size_t trackId = 0; trackId < cosTracks.size(); ++trackId) {
1684  FillTrackVars((*cosTracks[trackId]), rec.trk.cosmic.tracks[trackId],
1685  *slices, allTracks, sliceTracks, allTracksBirks,
1686  sliceTracksBirks, sliceId);
1687 
1688  // Start of MR Shower
1689  // Fill the MR Shower informations
1690  art::FindManyP<rb::Cluster> fmDiF(cosTracks, evt,
1692  if (fmDiF.isValid()) {
1693  std::vector<art::Ptr<rb::Cluster>> difs = fmDiF.at(trackId);
1694  for (unsigned int difId = 0; difId < difs.size(); ++difId) {
1695  SRMRProperties srMR; // we're going to accumulate info into this
1696  art::Ptr<rb::Cluster> dif = difs[difId];
1697  // Filling DiF Cluster variables
1698  FillDiFVars(*dif, *slices, allTracks, sliceTracks, allTracksBirks,
1699  sliceTracksBirks, sliceId, srMR);
1700 
1701  // -- CVN - Loose Presel PtP Cut
1702  art::FindManyP<cvn::Result> fmDiFCVNResult_LoosePreselPtP(
1704  cvn::Result cvnDiFResult_LoosePreselPtP;
1705  if (GetAssociatedProduct(fmDiFCVNResult_LoosePreselPtP, difId,
1706  cvnDiFResult_LoosePreselPtP)) {
1707  FillCVNResultVars(cvnDiFResult_LoosePreselPtP, srMR.cvnloosepreselptp);
1708  } else {
1710  }
1711  // -- CVN - Old Presel
1712  art::FindManyP<cvn::Result> fmDiFCVNResult_OldPresel(
1714  cvn::Result cvnDiFResult_OldPresel;
1715  if (GetAssociatedProduct(fmDiFCVNResult_OldPresel, difId,
1716  cvnDiFResult_OldPresel)) {
1717  FillCVNResultVars(cvnDiFResult_OldPresel, srMR.cvnoldpresel);
1718  } else {
1719  srMR.cvnoldpresel.setDefault();
1720  }
1721  // -- CVN - No Cosmics used in training
1722  art::FindManyP<cvn::Result> fmDiFCVNResult_NoCosmics(
1724  cvn::Result cvnDiFResult_NoCosmics;
1725  if (GetAssociatedProduct(fmDiFCVNResult_NoCosmics, difId,
1726  cvnDiFResult_NoCosmics)) {
1727  FillCVNResultVars(cvnDiFResult_NoCosmics, srMR.cvnnocosmics, true);
1728  } else {
1729  srMR.cvnnocosmics.setDefault();
1730  }
1731 
1732  // Filling DiF LID
1733  art::FindManyP<rb::Shower> fmDiFShw({dif}, evt,
1735  if (fmDiFShw.isValid()) {
1736  std::vector<art::Ptr<rb::Shower>> difshws = fmDiFShw.at(difId);
1737  for (art::Ptr<rb::Shower> difshw : difshws) {
1738  if (difshw) {
1739  FillDiFShowerVars(*difshw, *slices, allTracks, sliceTracks,
1740  allTracksBirks, sliceTracksBirks, sliceId,
1741  srMR);
1742  } // end for if difshw
1743  } // end for loop difshw
1744  } // end if Shw assn DiF valid
1745 
1747  {dif}, evt, fParams.DiFShowerLIDLabel());
1748  if (fmDiFLID.isValid()) {
1749  std::vector<art::Ptr<slid::ShowerLID>> diflids =
1750  fmDiFLID.at(difId);
1751  for (art::Ptr<slid::ShowerLID> diflid : diflids) {
1752  if (diflid) {
1753  FillSlidVars(*diflid, srMR.lid);
1754  } // end for if diflid
1755  } // end for loop diflid
1756  } // end if LID assn DiF valid
1757 
1758  // Now the DiF cluster is ready
1759  rec.trk.cosmic.tracks[trackId].mrdif.push_back(srMR);
1760  } // end for dif cluster loop
1761  } // end if Track assn with DiF is valid
1762  // End of MR Shower
1763 
1764  if (fmCosmicTrkME.isValid()) {
1765  // Get MichelE from MEFinder
1766  std::vector<art::Ptr<me::TrkME>> michels = fmCosmicTrkME.at(trackId);
1767  for (art::Ptr<me::TrkME> michel : michels) {
1768  AddTrkMEToVec(*michel, &rec.me.trkcosmic,
1769  rec.trk.cosmic.tracks.back(), *slices, allTracks,
1770  sliceTracks, allTracksBirks, sliceTracksBirks,
1771  sliceId);
1772  }
1773  } // end
1774 
1775  } // end for trackId
1776  }
1777 
1778  // Do the same for windowtrack tracks
1779  if (fmWindow.isValid()) {
1780  std::vector<art::Ptr<rb::Track>> winTracks = fmWindow.at(sliceId);
1781  // Sort the tracks in descending length order
1782  std::sort(winTracks.begin(), winTracks.end(), sortRBTrackLength);
1783  rec.trk.window.tracks.reserve(winTracks.size());
1784  rec.trk.window.tracks.resize(winTracks.size());
1785  for (size_t trackId = 0; trackId < winTracks.size(); ++trackId) {
1786  FillTrackVars((*winTracks[trackId]), rec.trk.window.tracks[trackId],
1787  *slices, allTracks, sliceTracks, allTracksBirks,
1788  sliceTracksBirks, sliceId);
1789  } // end for trackId
1790  }
1791 
1792  // Now do the same for kalmanTrack
1793  if (fmKalman.isValid()) {
1794  std::vector<art::Ptr<rb::Track>> kalTracks = fmKalman.at(sliceId);
1795  // Sort the tracks in descending length order to start
1796  std::sort(kalTracks.begin(), kalTracks.end(), sortRBTrackLength);
1797  rec.trk.kalman.idxlongest = 0;
1798 
1799  art::FindManyP<remid::ReMId> fmRemid(kalTracks, evt,
1800  fParams.RemidLabel());
1801 
1803  FindManyPStrict<cosrej::TrkCntObj>(kalTracks, evt,
1804  fParams.CosRejLabel());
1805  art::FindManyP<me::TrkME> fmKalmanME(kalTracks, evt, fParams.MELabel());
1806  art::FindManyP<rb::Energy> fmKalOverlapE(
1807  kalTracks, evt,fParams.OverlapEKalLabel());
1808 
1810  FindManyPStrict<trackinfo::TrackInfoObj>(kalTracks, evt,
1812 
1813  // Give SRRemid the pid values of track with largest PID
1814  unsigned int bestPidIdx =
1816 
1817  unsigned int bestPidIdx3d = 999;
1818 
1819 
1820 
1821  for (size_t trackId = 0; trackId < kalTracks.size(); ++trackId) {
1822  if (kalTracks[trackId]->View() != geo::kXorY) {
1823  rec.trk.kalman.tracks2d.push_back(SRTrack());
1824 
1825  // Fill the Track Vars with overlapping E if that info is present,
1826  // otherwise use the default FillTrackVars function.
1827  if (fmKalOverlapE.isValid()) {
1828  std::vector<art::Ptr<rb::Energy>> overlapEs =
1829  fmKalOverlapE.at(trackId);
1830  if (!overlapEs.empty()) {
1831  FillTrackVarsWithOverlapE((*kalTracks[trackId]), (*overlapEs[0]),
1832  rec.trk.kalman.tracks2d.back(), *slices,
1833  allTracks, sliceTracks, allTracksBirks,
1834  sliceTracksBirks, sliceId);
1835  } else {
1836  FillTrackVars((*kalTracks[trackId]),
1837  rec.trk.kalman.tracks2d.back(), *slices, allTracks,
1838  sliceTracks, allTracksBirks, sliceTracksBirks,
1839  sliceId);
1840  }
1841  } else {
1842  FillTrackVars((*kalTracks[trackId]), rec.trk.kalman.tracks2d.back(),
1843  *slices, allTracks, sliceTracks, allTracksBirks,
1844  sliceTracksBirks, sliceId);
1845  }
1846  if(fmTrkInfo.isValid()){
1847  std::vector<art::Ptr<trackinfo::TrackInfoObj> > trkinfos = fmTrkInfo.at(trackId);
1848  //if(!trkinfos.empty())
1849  // FillTrackInfoVars(*trkinfos[0], rec.trk.kalman.tracks.back());
1850  }
1851  } else {
1852  rec.trk.kalman.tracks.push_back(SRKalmanTrack());
1853 
1854  // Fill the Track Vars with overlapping E if that info is present,
1855  // otherwise use the default FillTrackVars function.
1856  if (fmKalOverlapE.isValid()) {
1857  std::vector<art::Ptr<rb::Energy>> overlapEs =
1858  fmKalOverlapE.at(trackId);
1859  if (!overlapEs.empty()) {
1860  FillTrackVarsWithOverlapE((*kalTracks[trackId]), (*overlapEs[0]),
1861  rec.trk.kalman.tracks.back(), *slices,
1862  allTracks, sliceTracks, allTracksBirks,
1863  sliceTracksBirks, sliceId);
1864  } else {
1865  FillTrackVars((*kalTracks[trackId]), rec.trk.kalman.tracks.back(),
1866  *slices, allTracks, sliceTracks, allTracksBirks,
1867  sliceTracksBirks, sliceId);
1868  }
1869  } else {
1870  FillTrackVars((*kalTracks[trackId]), rec.trk.kalman.tracks.back(),
1871  *slices, allTracks, sliceTracks, allTracksBirks,
1872  sliceTracksBirks, sliceId);
1873  }
1874 
1875  if(fmTrkInfo.isValid()){
1876  std::vector<art::Ptr<trackinfo::TrackInfoObj> > trkinfos = fmTrkInfo.at(trackId);
1877  if(!trkinfos.empty())
1878  FillTrackInfoVars(*trkinfos[0], rec.trk.kalman.tracks.back());
1879  }
1880 
1881  if (trackId == bestPidIdx)
1882  bestPidIdx3d = rec.trk.kalman.tracks.size() - 1;
1883 
1884  if (fmRemid.isValid()) {
1885  // Get MuID from ReMId
1886  std::vector<art::Ptr<remid::ReMId>> remids = fmRemid.at(trackId);
1887 
1888  if (!remids.empty()) {
1889  FillMuIdVars(*remids[0], rec.trk.kalman.tracks.back());
1890  } else {
1891  std::cout << "Missing Association Error: Failed to find a ReMId "
1892  "object with this track."
1893  << std::endl;
1894  abort();
1895  }
1896  }
1897 
1898  if (fmTrkCnt.isValid()) {
1899  // Get track containment info from TrkCntObj
1900  std::vector<art::Ptr<cosrej::TrkCntObj>> trkcnts =
1901  fmTrkCnt.at(trackId);
1902 
1903  if (!trkcnts.empty())
1904  FillTrackContainmentVars(*trkcnts[0],
1905  rec.trk.kalman.tracks.back());
1906  }
1907 
1908  if (fmKalmanME.isValid()) {
1909  // Get MichelE from MEFinder
1910  std::vector<art::Ptr<me::TrkME>> michels = fmKalmanME.at(trackId);
1911 
1912  for (art::Ptr<me::TrkME> michel : michels) {
1913  AddTrkMEToVec(*michel, &rec.me.trkkalman,
1914  rec.trk.kalman.tracks.back(), *slices, allTracks,
1915  sliceTracks, allTracksBirks, sliceTracksBirks,
1916  sliceId);
1917  }
1918  } // end
1919  }
1920  } // end for trackId
1921 
1922  if (!rec.trk.kalman.tracks.empty() && bestPidIdx3d != 999) {
1923  CopyRemidVars(rec.trk.kalman.tracks[bestPidIdx3d], rec.sel.remid);
1924  std::vector<art::Ptr<remid::ReMId>> remids = fmRemid.at(bestPidIdx3d);
1925  rec.trk.kalman.idxremid = bestPidIdx3d;
1926  } else {
1927  rec.sel.remid.setDefault();
1928  }
1929 
1930  } else {
1931  rec.sel.remid.setDefault();
1932  }
1933 
1934  std::sort(rec.trk.kalman.tracks.begin(), rec.trk.kalman.tracks.end(),
1935  sortRemid);
1936  rec.trk.kalman.idxremid = 0;
1937  unsigned int tempidx = 0, templen = 0;
1938  for (size_t trackId = 0; trackId < rec.trk.kalman.tracks.size();
1939  ++trackId) {
1940  if (rec.trk.kalman.tracks[trackId].len > templen) {
1941  tempidx = trackId;
1942  templen = rec.trk.kalman.tracks[trackId].len;
1943  }
1944  }
1945 
1946  rec.trk.kalman.idxlongest = tempidx;
1947  rec.trk.kalman.fillSizes();
1948  rec.trk.discrete.fillSizes();
1949  rec.trk.window.fillSizes();
1950  rec.trk.cosmic.fillSizes();
1951 
1952  // pre-fill nue energy variables with default
1953  rec.energy.nue.setDefault();
1954  rec.energy.nue.lid.setDefault();
1955 
1956  // get vertices from ElasticArms
1957  if (fmElastic.isValid()) {
1958 
1959 
1960  std::vector<art::Ptr<rb::Vertex>> elastics = fmElastic.at(sliceId);
1961  if(elastics.size() > 0){ //this should always be size 1 or 0
1962 
1963  rec.vtx.elastic = SRElastic();
1964  int vtxId = 0; //always check the 0th entry since there should only be 1.
1965  //will remove this and replace all vtxId with 0 once we are certain the
1966  //code all works.
1967  SRElastic& srVtx = rec.vtx.elastic;
1968  srVtx.time = elastics[vtxId]->GetT();
1969  srVtx.vtx = elastics[vtxId]->GetXYZ();
1970  srVtx.IsValid = true;
1971  // for each vertex (currently only one elastic arms ever), look for
1972  // prongs (currently only fuzzyk)
1973  art::FindManyP<rb::Prong> fmFuzzyProng2D = FindManyPStrict<rb::Prong>(
1974  elastics, evt, fParams.FuzzyK2DLabel()); // 2D prongs first
1975  if (fmFuzzyProng2D.isValid()) {
1976  std::vector<art::Ptr<rb::Prong>> prongs = fmFuzzyProng2D.at(vtxId);
1977  std::sort(prongs.begin(), prongs.end(),
1978  sortRBProngLength); // 2D prongs are sorted by length
1980  prongs, evt,fParams.FuzzyKWeight2DLabel());
1981 
1982  for (unsigned int iPng = 0; iPng < prongs.size(); ++iPng) {
1983  art::Ptr<rb::Prong> prong = prongs[iPng];
1984  srVtx.fuzzyk.png2d.push_back(SRProng());
1985  FillProngVars(*prong, srVtx.fuzzyk.png2d.back(), *slices, allTracks,
1986  sliceTracks, allTracksBirks, sliceTracksBirks,
1987  sliceId);
1988  art::FindManyP<rb::PID> fmCVNParticleResult(
1989  prongs, evt, fParams.CVNParticle2DLabel());
1990  if (fmCVNParticleResult.isValid()) {
1991  std::vector<art::Ptr<rb::PID>> cvnparts =
1992  fmCVNParticleResult.at(iPng);
1993  if(!cvnparts.empty()){
1994  FillCVNParticleResultVars(cvnparts, srVtx.fuzzyk.png2d.back().cvnpart);
1995  }
1996  }
1997  art::FindManyP<cvn::PixelMap> fmCVNProngPixelMaps(
1998  prongs, evt, fParams.CVNPixelMapLabel());
1999  art::FindManyP<cvn::ProngTrainingData> fmCVNProngTrainingData(
2000  prongs, evt, fParams.CVNTrainingDataLabel());
2001 
2002  if (fmCVNProngPixelMaps.isValid() && fParams.FillPixelMaps()) {
2003  std::vector<art::Ptr<cvn::PixelMap>> cvnProngPixelMaps =
2004  fmCVNProngPixelMaps.at(iPng);
2005  for (unsigned int ip = 0; ip < cvnProngPixelMaps.size(); ++ip) {
2006  srVtx.fuzzyk.png2d.back().cvnmaps.push_back(SRPixelMap());
2007  SRPixelMap& pmap = srVtx.fuzzyk.png2d.back().cvnmaps.back();
2008  FillCVNPixelMaps(*cvnProngPixelMaps[ip], pmap,
2010  }
2011  }
2012 
2013  cvn::ProngTrainingData cvnProngTrainingData;
2014  if (GetAssociatedProduct(fmCVNProngTrainingData, iPng,
2015  cvnProngTrainingData) &&
2017  srVtx.fuzzyk.png2d.back().prongtrainingdata.push_back(
2019  SRProngTrainingData& srProngTrainingData =
2020  srVtx.fuzzyk.png2d.back().prongtrainingdata.back();
2021  FillCVNProngTrainingData(cvnProngTrainingData,
2022  srProngTrainingData);
2023  }
2024 
2025  if (fmPngPngWt.isValid()) {
2026  std::vector<art::Ptr<rb::WeightedProng>> prongWt =
2027  fmPngPngWt.at(iPng);
2028  srVtx.fuzzyk.png2d.back().weightedCalE =
2029  prongWt[0]->CalorimetricEnergy();
2030  }
2031 
2032  } // end loop over 2D prongs
2033  } // end if 2D prongs are valid for this vertex
2034 
2035  // 3D prongs next
2036  art::FindManyP<rb::Prong> fmFuzzyProng3D = FindManyPStrict<rb::Prong>(
2037  elastics, evt,fParams.FuzzyK3DLabel());
2038  unsigned int nshwlids = 0;
2039  if (fmFuzzyProng3D.isValid()) {
2040  std::vector<art::Ptr<rb::Prong>> prongs = fmFuzzyProng3D.at(vtxId);
2041  if (!prongs.empty()) {
2042  std::sort(prongs.begin(), prongs.end(),
2043  sortRBProngLength); // start sorting by length
2044  unsigned int pngId = 0;
2045  srVtx.fuzzyk.longestidx = 0;
2046  art::FindOneP<rb::Shower> fmShowerLID(prongs, evt,
2048  art::FindManyP<rb::Track> fmBPF(prongs, evt,
2050  art::FindManyP<rb::PID> fmCVNParticleResult(
2051  prongs, evt, fParams.CVNParticleLabel());
2052  art::FindManyP<cvn::PixelMap> fmCVNProngPixelMaps(
2053  prongs, evt, fParams.CVNPixelMapLabel());
2054  art::FindManyP<cvn::ProngTrainingData> fmCVNProngTrainingData(
2055  prongs, evt, fParams.CVNTrainingDataLabel());
2056  art::FindManyP<cvn::RegProngResult> fmRegCVNResult(prongs, evt,
2057  fParams.RegCVNLabel());
2058 
2060  prongs, evt,fParams.FuzzyKWeight3DLabel());
2061  float highE = 0;
2062  for (unsigned int iPng = 0; iPng < prongs.size(); ++iPng) {
2063  art::Ptr<rb::Prong> prong = prongs[iPng];
2064  srVtx.fuzzyk.png.push_back(SRFuzzyKProng());
2065  srVtx.fuzzyk.png.back().setDefault();
2066  FillProngVars(*prong, srVtx.fuzzyk.png.back(), *slices, allTracks,
2067  sliceTracks, allTracksBirks, sliceTracksBirks,
2068  sliceId);
2069  if (fmPngPngWt.isValid()) {
2070  std::vector<art::Ptr<rb::WeightedProng>> prongWt =
2071  fmPngPngWt.at(iPng);
2072  srVtx.fuzzyk.png.back().weightedCalE =
2073  prongWt[0]->CalorimetricEnergy();
2074  }
2075 
2076  if (fmCVNProngPixelMaps.isValid() && fParams.FillPixelMaps()) {
2077  std::vector<art::Ptr<cvn::PixelMap>> cvnProngPixelMaps =
2078  fmCVNProngPixelMaps.at(iPng);
2079  for (unsigned int ip = 0; ip < cvnProngPixelMaps.size(); ++ip) {
2080  srVtx.fuzzyk.png.back().cvnmaps.push_back(SRPixelMap());
2081  SRPixelMap& pmap = srVtx.fuzzyk.png.back().cvnmaps.back();
2082  FillCVNPixelMaps(*cvnProngPixelMaps[ip], pmap,
2084  }
2085  }
2086 
2087  cvn::ProngTrainingData cvnProngTrainingData;
2088  if (GetAssociatedProduct(fmCVNProngTrainingData, iPng,
2089  cvnProngTrainingData) &&
2091  srVtx.fuzzyk.png.back().prongtrainingdata.push_back(
2093  SRProngTrainingData& srProngTrainingData =
2094  srVtx.fuzzyk.png.back().prongtrainingdata.back();
2095  FillCVNProngTrainingData(cvnProngTrainingData,
2096  srProngTrainingData);
2097  }
2098 
2099  if (fmCVNParticleResult.isValid()) {
2100  std::vector<art::Ptr<rb::PID>> cvnparts =
2101  fmCVNParticleResult.at(iPng);
2102  FillCVNParticleResultVars(cvnparts, srVtx.fuzzyk.png.back().cvnpart);
2103  }
2104  if (fmRegCVNResult.isValid()) {
2105  std::vector<art::Ptr<cvn::RegProngResult>> regcvn =
2106  fmRegCVNResult.at(iPng);
2107  if(regcvn.size()>0)
2108  srVtx.fuzzyk.png.back().regcvn.prongE = regcvn[0]->Output();
2109  else
2110  srVtx.fuzzyk.png.back().regcvn.prongE = -5.;
2111  }
2112  if (fmShowerLID.isValid()) {
2113  art::Ptr<rb::Shower> slid = fmShowerLID.at(pngId);
2114  // look for associated LID reconstructions
2115  art::FindOneP<slid::ShowerLID> fmLID({slid}, evt,
2116  fParams.LIDLabel());
2117  if (slid) {
2118  // fill fuzzyk.png.shwlid; assume 1 shower per prong
2119  FillShowerVars(*slid, srVtx.fuzzyk.png.back(), *slices,
2120  allTracks, sliceTracks, allTracksBirks,
2121  sliceTracksBirks, sliceId);
2122  nshwlids++;
2123  if (fmLID.isValid()) {
2124  art::Ptr<slid::ShowerLID> lid = fmLID.at(0);
2125  // fill fuzzyk.png.shwlid.(lid/lidE); assume 1 LID info
2126  // object per shower
2127  if (lid) {
2128  FillSlidVars(*lid, srVtx.fuzzyk.png.back().shwlid);
2129  FillLIDEnergyVars(*lid, srVtx.fuzzyk.png.back());
2130  srVtx.fuzzyk.png.back().shwlid.shwE =
2131  srVtx.fuzzyk.png.back().shwlid.lidE.shwE;
2132  if (srVtx.fuzzyk.png.back().shwlid.shwE > highE) {
2133  highE = srVtx.fuzzyk.png.back().shwlid.shwE;
2134  // FillNueEnergyVars(slice, *slid,
2135  //rec.energy.nue);
2136  }
2137  }
2138  } // end if LID assn valid
2139  }
2140  } // end if showerlid is valid
2141 
2142  if (fmBPF.isValid()) {
2143 
2144  // BPF attempts (and rarely fails) to fit each prong under 3 different
2145  // particle assumptions: muon, pion, proton. So for each prong, there
2146  // will be at most 3 tracks. Each track _should_ have one and only one
2147  // FitSum, OverlapE, and TrkCntObj associated with it.
2148 
2149  // So... ...get all of these objects etc.
2150  std::vector<art::Ptr<rb::Track>> bpfTracks = fmBPF.at(pngId);
2151  art::FindManyP<rb::FitSum> fmBPFFitSums(
2152  bpfTracks, evt, fParams.BreakPointTrackLabel());
2153  art::FindManyP<rb::Energy> fmBPFOverlapE(
2154  bpfTracks, evt,fParams.OverlapEBPFLabel());
2155  art::FindManyP<cosrej::TrkCntObj> fmTrkCntBPF =
2156  FindManyPStrict<cosrej::TrkCntObj>(
2157  bpfTracks, evt, fParams.CosRejBPFLabel());
2159  FindManyPStrict<trackinfo::TrackInfoObj>(
2160  bpfTracks, evt, fParams.TrackInfoBPF());
2161 
2162  for (size_t trackId = 0; trackId < bpfTracks.size();
2163  ++trackId) {
2164  if (fmBPFFitSums.isValid()) {
2165  std::vector<art::Ptr<rb::FitSum>> bpfFitSums =
2166  fmBPFFitSums.at(trackId);
2167  if (bpfFitSums.size() < 1) {
2168  std::cerr << "bpfFitSums.size() < 1, this should be 1, "
2169  "since each track has a fitsum!"
2170  << std::endl;
2171  continue;
2172  }
2173 
2174  //
2175  // Here we should fill the track and fitsum info
2176  //
2177 
2178  // For easy reference, pull out this specific prong.
2179  SRFuzzyKProng& srPng = srVtx.fuzzyk.png.back();
2180 
2181 
2182 
2183  // Sorry for the cut and paste, but it is only ~20 lines of actual code.
2184  // Fill the track if it is for the muon, then repeat for pion and proton.
2185 
2186  // If this is the muon, fill the muon track info...
2187  if(abs(bpfFitSums[0]->PDG()) == 13) {
2188 
2189  srPng.bpf.muon.IsValid = true;
2190 
2191  // Fill the Track Vars with overlapping E if that info is
2192  // present, otherwise use the default FillTrackVars
2193  // function.
2194  if (fmBPFOverlapE.isValid()) {
2195  std::vector<art::Ptr<rb::Energy>> overlapEs =
2196  fmBPFOverlapE.at(trackId);
2197  if (!overlapEs.empty()) {
2199  (*bpfTracks[trackId]), (*overlapEs[0]), srPng.bpf.muon,
2200  *slices, allTracks, sliceTracks, allTracksBirks,
2201  sliceTracksBirks, sliceId);
2202  } else {
2203  FillTrackVars((*bpfTracks[trackId]), srPng.bpf.muon,
2204  *slices, allTracks, sliceTracks,
2205  allTracksBirks, sliceTracksBirks,
2206  sliceId);
2207  }
2208  } else {
2209  FillTrackVars((*bpfTracks[trackId]), srPng.bpf.muon, *slices,
2210  allTracks, sliceTracks, allTracksBirks,
2211  sliceTracksBirks, sliceId);
2212  }
2213 
2214  // Fill the FitSum info
2215  FillTrackVarsBpfFitSum((*bpfFitSums[0]), srPng.bpf.muon);
2216  if (fmTrkCntBPF.isValid()) {
2217  // Get track containment info from TrkCntObj
2218  std::vector<art::Ptr<cosrej::TrkCntObj>> trkcnts =
2219  fmTrkCntBPF.at(trackId);
2220 
2221  if (!trkcnts.empty())
2222  FillTrackContainmentVars(*trkcnts[0], srPng.bpf.muon);
2223  }
2224 
2225  // Fill the TrackInfo info (yes I said info twice...)
2226  if(fmTrkInfoBPF.isValid()){
2227  std::vector<art::Ptr<trackinfo::TrackInfoObj> > trkinfos = fmTrkInfoBPF.at(trackId);
2228  if(!trkinfos.empty())
2229  FillTrackInfoVars(*trkinfos[0], srPng.bpf.muon);
2230  }
2231 
2232  } // end if(PDG == 13)
2233 
2234  // If this is the pion, fill the pion track info...
2235  else if(abs(bpfFitSums[0]->PDG()) == 211) {
2236 
2237  srPng.bpf.pion.IsValid = true;
2238 
2239  // Fill the Track Vars with overlapping E if that info is
2240  // present, otherwise use the default FillTrackVars
2241  // function.
2242  if (fmBPFOverlapE.isValid()) {
2243  std::vector<art::Ptr<rb::Energy>> overlapEs =
2244  fmBPFOverlapE.at(trackId);
2245  if (!overlapEs.empty()) {
2247  (*bpfTracks[trackId]), (*overlapEs[0]), srPng.bpf.pion,
2248  *slices, allTracks, sliceTracks, allTracksBirks,
2249  sliceTracksBirks, sliceId);
2250  } else {
2251  FillTrackVars((*bpfTracks[trackId]), srPng.bpf.pion,
2252  *slices, allTracks, sliceTracks,
2253  allTracksBirks, sliceTracksBirks,
2254  sliceId);
2255  }
2256  } else {
2257  FillTrackVars((*bpfTracks[trackId]), srPng.bpf.pion, *slices,
2258  allTracks, sliceTracks, allTracksBirks,
2259  sliceTracksBirks, sliceId);
2260  }
2261 
2262  // Fill the FitSum info
2263  FillTrackVarsBpfFitSum((*bpfFitSums[0]), srPng.bpf.pion);
2264  if (fmTrkCntBPF.isValid()) {
2265  // Get track containment info from TrkCntObj
2266  std::vector<art::Ptr<cosrej::TrkCntObj>> trkcnts =
2267  fmTrkCntBPF.at(trackId);
2268 
2269  if (!trkcnts.empty())
2270  FillTrackContainmentVars(*trkcnts[0], srPng.bpf.pion);
2271  }
2272 
2273  // Fill the TrackInfo info (yes I said info twice...)
2274  if(fmTrkInfoBPF.isValid()){
2275  std::vector<art::Ptr<trackinfo::TrackInfoObj> > trkinfos = fmTrkInfoBPF.at(trackId);
2276  if(!trkinfos.empty())
2277  FillTrackInfoVars(*trkinfos[0], srPng.bpf.pion);
2278  }
2279 
2280  } // end if(PDG == 211)
2281 
2282  // If this is the proton, fill the proton track info...
2283  else if(abs(bpfFitSums[0]->PDG()) == 2212) {
2284 
2285  srPng.bpf.proton.IsValid = true;
2286 
2287  // Fill the Track Vars with overlapping E if that info is
2288  // present, otherwise use the default FillTrackVars
2289  // function.
2290  if (fmBPFOverlapE.isValid()) {
2291  std::vector<art::Ptr<rb::Energy>> overlapEs =
2292  fmBPFOverlapE.at(trackId);
2293  if (!overlapEs.empty()) {
2295  (*bpfTracks[trackId]), (*overlapEs[0]), srPng.bpf.proton,
2296  *slices, allTracks, sliceTracks, allTracksBirks,
2297  sliceTracksBirks, sliceId);
2298  } else {
2299  FillTrackVars((*bpfTracks[trackId]), srPng.bpf.proton,
2300  *slices, allTracks, sliceTracks,
2301  allTracksBirks, sliceTracksBirks,
2302  sliceId);
2303  }
2304  } else {
2305  FillTrackVars((*bpfTracks[trackId]), srPng.bpf.proton, *slices,
2306  allTracks, sliceTracks, allTracksBirks,
2307  sliceTracksBirks, sliceId);
2308  }
2309 
2310  // Fill the FitSum info
2311  FillTrackVarsBpfFitSum((*bpfFitSums[0]), srPng.bpf.proton);
2312  if (fmTrkCntBPF.isValid()) {
2313  // Get track containment info from TrkCntObj
2314  std::vector<art::Ptr<cosrej::TrkCntObj>> trkcnts =
2315  fmTrkCntBPF.at(trackId);
2316 
2317  if (!trkcnts.empty())
2318  FillTrackContainmentVars(*trkcnts[0], srPng.bpf.proton);
2319  }
2320 
2321  // Fill the TrackInfo info (yes I said info twice...)
2322  if(fmTrkInfoBPF.isValid()){
2323  std::vector<art::Ptr<trackinfo::TrackInfoObj> > trkinfos = fmTrkInfoBPF.at(trackId);
2324  if(!trkinfos.empty())
2325  FillTrackInfoVars(*trkinfos[0], srPng.bpf.proton);
2326  }
2327 
2328  } // end if(PDG == 2212)
2329 
2330  // otherwise, the PDG code was not one of the three track assumptions
2331  else {
2332  // throw a warning but continue...
2333  std::cout << "\n\nWARNING: BPF FitSum object contained non-standard PDG code "
2334  << bpfFitSums[0]->PDG()
2335  << ". Continuing but not filling info for this track.\n\n";
2336  } // end else
2337 
2338  } // fmBPFFitSums.isValid()
2339  else std::cerr << "fmBPFFitSums is not valid" << std::endl;
2340 
2341  } // loop trackId
2342 
2343  } // fmBPF.isValid()
2344 
2345  pngId++;
2346 
2347  } // end loop over 3D prongs
2348  } // there are non-zero 3D prongs
2349  rec.vtx.elastic.fuzzyk.nshwlid = nshwlids;
2350  rec.vtx.elastic.fuzzyk.fillSizes();
2351 
2352  // energy from prong orphaned hits
2353  art::FindManyP<rb::Prong> fmProngOrph(
2354  elastics, evt,fParams.FuzzyKOrphLabel());
2355  if (fmProngOrph.isValid()) {
2356  std::vector<art::Ptr<rb::Prong>> prongsOrph = fmProngOrph.at(vtxId);
2357  if (prongsOrph.size() > 0) {
2358  rec.vtx.elastic.fuzzyk.orphCalE =
2359  prongsOrph[0]->CalorimetricEnergy();
2360  }
2361  }
2362 
2363  std::sort(srVtx.fuzzyk.png.begin(), srVtx.fuzzyk.png.end(), sortProngCalE); // Use prong calE to sort
2364 
2365  float longest = 0;
2366  for (unsigned int pngIdx = 0; pngIdx < srVtx.fuzzyk.png.size();
2367  pngIdx++) {
2368  if (srVtx.fuzzyk.png[pngIdx].len > longest) {
2369  longest = srVtx.fuzzyk.png[pngIdx].len;
2370  srVtx.fuzzyk.longestidx = pngIdx; // note the longest prong index
2371  }
2372  }
2373 
2374  if (nshwlids >
2375  0) { // fill nue energy for primary shower if there were showers
2376  rec.energy.nue.lid.E =
2377  rec.vtx.elastic.fuzzyk.png[0].shwlid.lidE.E;
2378  rec.energy.nue.lid.depE =
2379  rec.vtx.elastic.fuzzyk.png[0].shwlid.lidE.depE;
2380  rec.energy.nue.lid.shwE =
2381  rec.vtx.elastic.fuzzyk.png[0].shwlid.lidE.shwE;
2382  rec.energy.nue.lid.hadE =
2383  rec.vtx.elastic.fuzzyk.png[0].shwlid.lidE.hadE;
2384  }
2385 
2386  } // end if 3D prongs are valid for this vertex
2387  } // end loop over vertices
2388 
2389  } // end ElasticArms
2390 
2391  FillNueEnergyVars(slice, rec.vtx, rec.energy.nue);
2392 
2393  // HoughVertex
2394 
2395  if (fmHoughVertex.isValid()) {
2396  std::vector<art::Ptr<rb::Vertex>> houghs = fmHoughVertex.at(sliceId);
2397 
2398  for (unsigned int vtxId = 0; vtxId < houghs.size(); ++vtxId) {
2399  rec.vtx.hough.push_back(SRHoughVertex());
2400  SRHoughVertex& srVtx = rec.vtx.hough.back();
2401 
2402  srVtx.time = houghs[vtxId]->GetT();
2403  srVtx.vtx = houghs[vtxId]->GetXYZ();
2404 
2405  /* if(fmFuzzyProng3D.isValid()){
2406  std::vector< art::Ptr<rb::Prong> > prongs = fmFuzzyProng3D.at(vtxId);
2407  // Sort prongs in descending length order
2408  std::sort(prongs.begin(), prongs.end(),sortRBProngLength);
2409  for(art::Ptr<rb::Prong> prong: prongs){
2410  if(prong->NXCell() == 0 || prong->NYCell() == 0){
2411  srVtx.fuzzyk.png2d.push_back(SRProng());
2412  FillProngVars(*prong, srVtx.fuzzyk.png2d.back(), slicess,
2413  sliceId);
2414  }// end 2d prong case
2415  else{
2416  srVtx.fuzzyk.png.push_back(SRProng());
2417  FillProngVars(*prong, srVtx.fuzzyk.png.back(), *slices, sliceId);
2418  }// end 3d prong case
2419  }// end loop over prongs
2420  }// end if prongs are valid for this vertex
2421  // Check if there is a valid instance label for 2D prongs
2422  if(fmFuzzyProng2D.isValid() && fFuzzyK2dAssnLabel != ""){
2423  std::vector< art::Ptr<rb::Prong> > prongs = fmFuzzyProng2D.at(vtxId);
2424  // Sort prongs in descending length order
2425  std::sort(prongs.begin(), prongs.end(),sortRBProngLength);
2426  for(art::Ptr<rb::Prong> prong: prongs){
2427  if(prong->NXCell() == 0 || prong->NYCell() == 0){
2428  srVtx.fuzzyk.png2d.push_back(SRProng());
2429  FillProngVars(*prong, srVtx.fuzzyk.png2d.back(), *slices,
2430  sliceId);
2431  }// end 2d prong case
2432  else{
2433  srVtx.fuzzyk.png.push_back(SRProng());
2434  FillProngVars(*prong, srVtx.fuzzyk.png.back(), *slices, sliceId);
2435  }// end 3d prong case
2436  }// end loop over prongs
2437  }// end if prongs are valid for this vertex
2438  */
2439  } // end for vtxId
2440  // Sort the vector
2441  // std::sort(rec.vtx.hough.begin(), rec.vtx.hough.end(),
2442  // sortHoughProng);
2443  }
2444 
2445  // Get all the VertexDT vertices associated with the slice
2446  if (fmVertexDT.isValid()) {
2447  std::vector<art::Ptr<rb::Vertex>> vtxs = fmVertexDT.at(sliceId);
2448  for (unsigned int vtxId = 0; vtxId < vtxs.size(); ++vtxId) {
2449  rec.vtx.vdt.push_back(SRVertexDT());
2450  SRVertexDT& srVtx = rec.vtx.vdt.back();
2451  srVtx.time = vtxs[vtxId]->GetT();
2452  srVtx.vtx = vtxs[vtxId]->GetXYZ();
2453  }
2454  }
2455  rec.vtx.fillSizes();
2456 
2457  // Get all MEFs aassociated with slice
2458  if (fmSlcME.isValid()) {
2459  std::vector<art::Ptr<me::SlcME>> michels = fmSlcME.at(sliceId);
2460  for (art::Ptr<me::SlcME> michel : michels) {
2461  AddSlcMEToVec(*michel, &rec.me.slc, *slices, allTracks, sliceTracks,
2462  allTracksBirks, sliceTracksBirks, sliceId);
2463  }
2464  }
2465  rec.me.fillSizes();
2466 
2467  // Get slice level PID from LIDBuilder
2468  if (foLID.isValid()) {
2469  cet::maybe_ref<art::Ptr<slid::EventLID> const> rlid(foLID.at(sliceId));
2470  art::Ptr<slid::EventLID> elid = rlid.ref();
2471  // If the pointer is valid
2472  if (elid) {
2473  rec.sel.lid.ann = elid->Value();
2474  rec.sel.lid.anne = elid->ValueWithE();
2475  rec.sel.lid.annepi0 = elid->ValueEPi0();
2476  rec.sel.lid.annecos = elid->ValueECos();
2477  } else {
2478  rec.sel.lid.setDefault();
2479  }
2480  }
2481 
2482  if (foSlcLID.isValid()) {
2484  foSlcLID.at(sliceId));
2485  art::Ptr<slid::SliceLID> slclid = rslclid.ref();
2486  if (slclid) {
2487  rec.sel.lid.rnn = slclid->Value();
2488  } else {
2489  rec.sel.lid.setDefault();
2490  }
2491  }
2492 
2493  /*
2494  if(fmShwLID.isValid()){
2495  std::vector<const slid::ShowerLID* > shwlids = fmShwLID.at(sliceId);
2496  // Sort here since the vector may not be
2497  sort(shwlids.begin(),shwlids.end(), slid::CompareByEnergy);
2498  art::FindOneP<rb::Shower> fos(shwlids, evt, fLIDLabel());
2499  for(unsigned int shwId = 0; shwId < shwlids.size(); ++shwId){
2500  cet::maybe_ref<art::Ptr<rb::Shower> const> rshw(fos.at(shwId));
2501  art::Ptr<rb::Shower> shw = rshw.ref();
2502  if (shw){
2503  rec.shw.shwlid.push_back(SRShower());
2504  FillShowerVars(*shw, rec.shw.shwlid.back(), *slices, sliceId);
2505  }
2506 
2507  rec.sel.elecid.shwlid.push_back(SRSLid());
2508  FillSlidVars(*shwlids[shwId], rec.sel.elecid.shwlid.back().shwlid);
2509 
2510  // Fill energy information for leading energy shower
2511  // which also serves as a slice pid
2512  if (shwId == 0){
2513  FillLIDEnergyVars(*shwlids[shwId], rec.energy.elecid);
2514  }
2515 
2516  // Fill energy vector with the other shower energies
2517  rec.energy.shwlid.push_back(SRSLidEnergy());
2518  FillLIDEnergyVars(*shwlids[shwId], rec.energy.shwlid.back());
2519 
2520  if (shw && shwId == 0){
2521  FillNueEnergyVars(slice, *shw, rec.energy.nue);
2522  }
2523  } // end for loop on showerlids
2524  } // end if showerlids are valid
2525 
2526  if(fmShwPID.isValid()){
2527  std::vector<const slid::ShowerPID* > shwpids = fmShwPID.at(sliceId);
2528  // Sort here since the vector may not be
2529  sort(shwpids.begin(),shwpids.end(), slid::CompareByEnergy);
2530  art::FindOneP<rb::Shower> fos(shwpids, evt, fSPIDLabel());
2531  for(unsigned int shwId = 0; shwId < shwpids.size(); ++shwId){
2532  rec.sel.elecid.shwpid.push_back(SRSPid());
2533  FillSpidVars(*shwpids[shwId], rec.sel.elecid.shwpid.back());
2534  } // end for loop on showerpids
2535  } // end if showerpids are valid
2536 
2537  } // end if event id is valid
2538  else{
2539  rec.sel.lid .setDefault();
2540  rec.energy.nue.lid .setDefault();
2541  }
2542 } // end if slice lids are valid
2543 else{
2544  rec.sel.lid .setDefault();
2545  rec.energy.lid .setDefault();
2546 }
2547 
2548 rec.shw.fillSizes();
2549 rec.energy.fillSizes();
2550 rec.sel.elecid.fillSizes();
2551  */
2553  bool isLemmed = GetAssociatedProduct(fmLem, sliceId, lem);
2554  std::vector<std::string> vec = {fParams.LEMNuePresel()};
2555 
2556  if (isLemmed && !(fParams.ApplyLEMNuePresel() &&
2557  rb::IsFiltered(evt, slices, sliceId, vec))) {
2558  FillLEMVars(lem, rec.sel.lem);
2559  } else
2560  rec.sel.lem.setDefault();
2561 
2562  rvp::RVP rvp;
2563  if (GetAssociatedProduct(fmRvp, sliceId, rvp))
2564  FillRVPVars(rvp, rec.sel.rvp);
2565  else
2566  rec.sel.rvp.setDefault();
2567 
2568  xnue::Xnue xnue;
2569  if (GetAssociatedProduct(fmXnue, sliceId, xnue)) {
2570  FillXnueVars(xnue, rec.sel.xnuepid);
2571  } else
2572  rec.sel.xnuepid.setDefault();
2573 
2574  numue::NumuE numuE;
2575  if (GetAssociatedProduct(fmNumuEnergy, sliceId, numuE)) {
2576  FillNumuEnergyVars(numuE, rec.energy.numu);
2577  } else {
2578  rec.energy.numu.setDefault();
2579  }
2580 
2581 
2582  cvn::RegNuResult evtregcvn;
2583  if (GetAssociatedProduct(fmEvtRegCVN, sliceId, evtregcvn)) {
2584  rec.energy.nue.regcvnEvtE = evtregcvn.Output();
2585  } else {
2586  rec.energy.nue.regcvnEvtE = -5.;
2587  }
2588 
2589 
2590 
2591  LSTME::LSTMEnergy numulstmE;
2592  if (GetAssociatedProduct(fmNumuLSTMEnergy, sliceId, numulstmE)) {
2593  FillNumuLSTMEnergyVars(numulstmE, rec.energy.numu);
2594  } else {
2595  rec.energy.numu.setLSTMDefault();
2596  }
2597 
2598  cvn::RegHadronResult hadregcvn;
2599  if (GetAssociatedProduct(fmHadRegCVN, sliceId, hadregcvn)) {
2600  rec.energy.numu.regcvnhadE = hadregcvn.Output();
2601  } else {
2602  rec.energy.numu.regcvnhadE = -5.;
2603  }
2604 
2605  SliceLID::Prediction sliceLID;
2606  if (GetAssociatedProduct(fmSliceLID, sliceId, sliceLID)) {
2607  FillSliceLID(sliceLID, rec.sel.slicelid);
2608  } else {
2609  rec.sel.slicelid.setDefault();
2610  }
2611 
2613  if (GetAssociatedProduct(fmCosRej, sliceId, cosrej)) {
2614  FillCosRejVars(cosrej, rec.sel.cosrej, rec.sel.contain);
2615 
2616  float kDistAllTop = rec.sel.nuecosrej.distallpngtop;
2617  if (std::isnan(kDistAllTop)) kDistAllTop = -1000.0;
2618 
2620  if (std::isnan(kDistAllBottom)) kDistAllBottom = -1000.0;
2621 
2623  if (std::isnan(kDistAllEast)) kDistAllEast = -1000.0;
2624 
2626  if (std::isnan(kDistAllWest)) kDistAllWest = -1000.0;
2627 
2629  if (std::isnan(kDistAllBack)) kDistAllBack = -1000.0;
2630 
2632  if (std::isnan(kDistAllFront)) kDistAllFront = -1000.0;
2633 
2634  if (fDet ==
2635  kFARDET) { // adopted from CAFAna/Cuts/NumuCuts2017.h on 8.16.17
2636  std::pair<int, int> planes = calcFirstLastLivePlane(
2637  rec.slc.firstplane, std::bitset<14>(rec.hdr.dibmask));
2638  int planestofront = rec.slc.firstplane - planes.first;
2639  int planestoback = planes.second - rec.slc.lastplane;
2640  bool kcut1 =
2641  (rec.sel.contain.kalfwdcell > 6 && rec.sel.contain.kalbakcell > 6 &&
2642  rec.sel.contain.cosfwdcell > 0 && rec.sel.contain.cosbakcell > 7 &&
2643  planestofront > 1 && planestoback > 1);
2644  bool kcut2 = kDistAllTop > 60 && kDistAllBottom > 12 &&
2645  kDistAllEast > 16 && kDistAllWest > 12 &&
2646  kDistAllFront > 18 && kDistAllBack > 18;
2647  rec.sel.contain.numucontain = (kcut1 && kcut2);
2648 
2649  if (cosrej.CosFwdCell() > 0 && cosrej.CosBakCell() > 0 &&
2650  cosrej.KalFwdCell() > 10 && cosrej.KalBakCell() > 10 &&
2651  rec.slc.ncellsfromedge > 1 && rec.sel.contain.nplanestofront > 1 &&
2652  rec.sel.contain.nplanestoback > 1) {
2653  rec.sel.contain.numucontainSA = true;
2654  } else
2655  rec.sel.contain.numucontainSA = false;
2656  }
2657 
2658  else if (fDet == kNEARDET) { // adopted from CAFAna/Cuts/NumuCuts2017.h
2659  // on 8.16.17
2660  if (rec.vtx.elastic.IsValid < 1)
2661  rec.sel.contain.numucontain = false;
2662  else {
2663  for (unsigned int i = 0; i < rec.vtx.elastic.fuzzyk.nshwlid; ++i) {
2664  TVector3 start = rec.vtx.elastic.fuzzyk.png[i].shwlid.start;
2665  TVector3 stop = rec.vtx.elastic.fuzzyk.png[i].shwlid.stop;
2666  if (std::min(start.X(), stop.X()) < -180.0)
2667  rec.sel.contain.numucontain = false;
2668  if (std::max(start.X(), stop.X()) > 180.0)
2669  rec.sel.contain.numucontain = false;
2670  if (std::min(start.Y(), stop.Y()) < -180.0)
2671  rec.sel.contain.numucontain = false;
2672  if (std::max(start.Y(), stop.Y()) > 180.0)
2673  rec.sel.contain.numucontain = false;
2674  if (std::min(start.Z(), stop.Z()) < 20.0)
2675  rec.sel.contain.numucontain = false;
2676  if (std::max(start.Z(), stop.Z()) > 1525.0)
2677  rec.sel.contain.numucontain = false;
2678  }
2679  if (rec.trk.kalman.ntracks < 1) rec.sel.contain.numucontain = false;
2680  for (unsigned int i = 0; i < rec.trk.kalman.ntracks; ++i) {
2681  if (i == rec.trk.kalman.idxremid)
2682  continue;
2683  else if (rec.trk.kalman.tracks[i].start.Z() > 1275 ||
2684  rec.trk.kalman.tracks[i].stop.Z() > 1275)
2685  rec.sel.contain.numucontain = false;
2686  }
2687  }
2688 
2689  rec.sel.contain.numucontain =
2690  (rec.trk.kalman.ntracks > rec.trk.kalman.idxremid &&
2691  rec.slc.firstplane > 1 && rec.slc.lastplane < 212 &&
2692  rec.trk.kalman.tracks[0].start.Z() < 1100 &&
2693  (rec.trk.kalman.tracks[0].stop.Z() < 1275 ||
2694  rec.sel.contain.kalyposattrans < 55) &&
2695  rec.sel.contain.kalfwdcellnd > 5 &&
2696  rec.sel.contain.kalbakcellnd > 10);
2697 
2698  if (cosrej.KalFwdCellND() > 4 && cosrej.KalBakCellND() > 8 &&
2699  rec.slc.ncellsfromedge > 1 && rec.slc.firstplane > 1 &&
2700  rec.slc.lastplane < 212 &&
2701  rec.trk.kalman.tracks[0].start.Z() < 1150 &&
2702  (rec.trk.kalman.tracks[0].stop.Z() < 1275 ||
2703  rec.sel.contain.kalyposattrans < 55) &&
2705  0.03) {
2706  rec.sel.contain.numucontainSA = true;
2707  } else
2708  rec.sel.contain.numucontainSA = false;
2709  }
2710  }
2711 
2712  cosrej::NueCosRej nuecosrej;
2713  if (GetAssociatedProduct(fmNueCosRej, sliceId, nuecosrej)) {
2714  FillNueCosRejVars(nuecosrej, rec.sel.nuecosrej);
2715  /*
2716  // nue containment boolean - nue group please feel free to update.
2717  // If variables move from sandbox this needs to be changed too
2718  if(fDet == kFARDET){
2719  if(nuecosrej.ProngMaxX() < 750 && nuecosrej.ProngMinX() > -750 &&
2720  nuecosrej.ProngMaxY() < 750 && nuecosrej.ProngMinY() > -750 &&
2721  std::min(nuecosrej.StartDistToFront(),
2722  nuecosrej.StopDistToFront()) > 50 &&
2723  std::min(nuecosrej.StartDistToBack(),
2724  nuecosrej.StopDistToBack()) > 50){
2725  rec.sel.contain.nuecontain = true;
2726  }
2727  }
2728 
2729  if(fDet == kNEARDET){
2730  if(nuecosrej.ProngMaxX() < 185 && nuecosrej.ProngMinX() > -185 &&
2731  nuecosrej.ProngMaxY() < 185 && nuecosrej.ProngMinY() > -185 &&
2732  std::min(nuecosrej.StartDistToFront(),
2733  nuecosrej.StopDistToFront()) > 50 &&
2734  std::min(nuecosrej.StartDistToBack(),
2735  nuecosrej.StopDistToBack()) > 50){
2736  rec.sel.contain.nuecontain = true;
2737  }
2738  }
2739  */
2740  } else {
2741  rec.sel.nuecosrej.setDefault();
2742  }
2743 
2744  ncid::NCCosRej nccosrej; ///////////// NC Cosmic Rejection
2745  if (GetAssociatedProduct(fmNCCosRej, sliceId, nccosrej)) {
2746  FillNCCosRejVars(nccosrej, rec.sel.nccosrej);
2747  } else {
2748  rec.sel.nccosrej.setDefault();
2749  }
2750 
2751  ncpi0::NCPi0BkgRej ncpi0bkgrej;
2752  if (GetAssociatedProduct(fmNCPi0BkgRej, sliceId, ncpi0bkgrej)) {
2753  FillNCPi0BkgRejVars(ncpi0bkgrej, rec.sel.ncpi0bkgrej);
2754  } else {
2755  rec.sel.ncpi0bkgrej.setDefault();
2756  }
2757 
2758  presel::PreselObj nuepre;
2759  if (GetAssociatedProduct(fmNuePresel, sliceId, nuepre)) {
2760  FillNuePreselVars(nuepre, rec.sel.nuepre);
2761  } else {
2762  rec.sel.nuepre.setDefault();
2763  }
2764 
2765  presel::PreselObj rockpre;
2766  if (GetAssociatedProduct(fmRockPresel, sliceId, rockpre)) {
2767  FillRockPreselVars(rockpre, rec.sel.rockpre);
2768  } else {
2769  rec.sel.rockpre.setDefault();
2770  }
2771 
2772  presel::Veto veto, nueveto;
2773  if (GetAssociatedProduct(fmVeto, sliceId, veto) ||
2774  GetAssociatedProduct(fmNueVeto, sliceId, nueveto)) {
2775  FillVetoVars(veto, nueveto, rec.sel.veto);
2776  } else {
2777  rec.sel.veto.setDefault();
2778  }
2779 
2780  // --- Horrible hack to appease CAFAna.
2781  rec.sel.cvn.setDefault();
2782 
2783  cvn::Result nuoneResult;
2784  if (GetAssociatedProduct(fmNuonE, sliceId, nuoneResult)) {
2785  FillNuonEResultVars(nuoneResult, rec.sel.nuone);
2786  } else {
2787  rec.sel.nuone.setDefault();
2788  }
2789 
2790  // --- Loose Presel PtP Cut
2791  cvn::Result cvnResult_LoosePreselPtP;
2792  if (GetAssociatedProduct(fmCVN_LoosePreselPtP, sliceId, cvnResult_LoosePreselPtP)) {
2793  FillCVNResultVars(cvnResult_LoosePreselPtP, rec.sel.cvnloosepreselptp);
2794  } else {
2796  }
2797  // --- Old Presel
2798  cvn::Result cvnResult_OldPresel;
2799  if (GetAssociatedProduct(fmCVN_OldPresel, sliceId, cvnResult_OldPresel)) {
2800  FillCVNResultVars(cvnResult_OldPresel, rec.sel.cvnoldpresel);
2801  } else {
2802  rec.sel.cvnoldpresel.setDefault();
2803  }
2804  // --- No Cosmics used in training
2805  cvn::Result cvnResult_NoCosmics;
2806  if (GetAssociatedProduct(fmCVN_NoCosmics, sliceId, cvnResult_NoCosmics)) {
2807  FillCVNResultVars(cvnResult_NoCosmics, rec.sel.cvnnocosmics);
2808  } else {
2809  rec.sel.cvnnocosmics.setDefault();
2810  }
2811  // --- Wrong Sign
2812  cvn::Result cvnResult_ws;
2813  if (GetAssociatedProduct(fmWS, sliceId, cvnResult_ws)) {
2814  rec.sel.wsid = cvnResult_ws.fOutput[0];
2815  }
2816 
2817  cvn::Features cvnFeatures;
2818  if (GetAssociatedProduct(fmCVNFeatures, sliceId, cvnFeatures) &&
2820  FillCVNFeaturesVars(cvnFeatures, rec.training.cvnfeatures,
2822  }
2823 
2824  std::cout<<"------------> wwj: FillCVNPixelMaps"<<std::endl;
2825  std::cout<<fmCVNPixelMaps.size()<<std::endl;
2826  if (fmCVNPixelMaps.isValid() && fParams.FillPixelMaps()) {
2827  std::vector<art::Ptr<cvn::PixelMap>> cvnPixelMaps =
2828  fmCVNPixelMaps.at(sliceId);
2829  std::cout<<cvnPixelMaps.size()<<std::endl;
2830  for (unsigned int ip = 0; ip < cvnPixelMaps.size(); ++ip) {
2831  rec.training.cvnmaps.push_back(SRPixelObjMap());
2832  SRPixelObjMap& pmap = rec.training.cvnmaps.back();
2833  FillCVNPixelMaps(*cvnPixelMaps[ip], pmap, fParams.UseGeVPixelMaps());
2834  }
2835  }
2836 
2837  std::cout<<"------------> wwj: FillCVNTrainingData"<<std::endl;
2838  cvn::TrainingData cvnTrainingData;
2839  if (GetAssociatedProduct(fmCVNTrainingData, sliceId, cvnTrainingData) &&
2841  rec.training.trainingdata.push_back(SRTrainingData());
2842  SRTrainingData& srTrainingData = rec.training.trainingdata.back();
2843  FillCVNTrainingData(cvnTrainingData, srTrainingData);
2844  }
2845 
2846  // also store every hit in the slice. Doing this in a fixed size array the
2847  // size of the detector since the zeros compress better
2848  if (fParams.FillPixelMaps()) {
2849  rec.training.slicemaps.push_back(SRSliceMap());
2850  SRSliceMap& smap = rec.training.slicemaps.back();
2851  float corr = 15123.7; // Rough PE to GeV? Rev 17083
2852  float peCorrChunk;
2853  // This value is 2.0 / (PE/GeV)
2854  if (fParams.UseGeVPixelMaps())
2855  peCorrChunk = 1.322e-4 / 255.0;
2856  else
2857  peCorrChunk = 2.0 / 255.0;
2858  for (unsigned int ih = 0; ih < slice.NCell(); ih++) {
2859  rb::RecoHit rh = slice.RecoHit(ih);
2860  art::Ptr<rb::CellHit> ch = slice.Cell(ih);
2861  // convert planes to 0 to 447 in each view
2862  unsigned int pl = (ch->Plane() / 2);
2863  unsigned int chan = ch->Cell() + 384 * (pl + 448 * ch->View());
2864  // convert cell energy to an unsigned char. Scales same as in CVN
2865  float hitPE = ch->PE();
2866  float gev = (rh.IsCalibrated()) ? rh.GeV() : hitPE / corr;
2867  float pe = (rh.IsCalibrated()) ? rh.PECorr() : hitPE;
2868  float tv = (fParams.UseGeVPixelMaps()) ? gev : pe;
2869  tv /= 2100.0; // scale factor from CVN Pixel maps
2870  float truncateCorr = ceil(tv / peCorrChunk);
2871  unsigned char value;
2872  if (!fParams.UseEnergy()) {
2873  if (truncateCorr > 0)
2874  truncateCorr = 1;
2875  else
2876  truncateCorr = 0;
2877  }
2878  if (truncateCorr > 255)
2879  value = 255;
2880  else
2881  value = (unsigned char)truncateCorr;
2882  smap.slicemap[chan] = value;
2883  }
2884  }
2885 
2887  if (GetAssociatedProduct(fmMrccPar, sliceId, parent))
2888  FillMRCCParentInfo(parent, rec.parent.mrccpar, evt);
2889  else
2890  rec.parent.mrccpar.setDefault();
2891 
2892  //#######################################################
2893  // Fill truth information
2894  //#######################################################
2895 
2896  // Set mc branch values to default
2897  rec.mc.setDefault();
2898  if (bt->HaveTruthInfo()) {
2899  // Get all truth track hits associated with slice
2900 
2901  int currentNuIdx = 0;
2902  unsigned int nNuGlobal =
2903  0; // count of neutrino interactions in this slice
2904 
2905  for (unsigned int iEffPur = 0; iEffPur < sliceTruthTable[sliceId].size();
2906  ++iEffPur) {
2907  // continue if truth left no hits in slice
2908  if (!(sliceTruthTable[sliceId][iEffPur].nSliceHits > 0)) continue;
2909 
2910  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
2911  faveIdsEnergy[sliceId]) {
2912  cheat::NeutrinoEffPur* faveNu = &sliceTruthTable[sliceId][iEffPur];
2913 
2914  // If NOT Neutrinos
2915  if (!faveNu->neutrinoInt->NeutrinoSet()) {
2916  // Add faveNu to allcosmics
2917  if (fParams.FullTruthInfo())
2918  AddCosmicTruthToVec(faveNu, allTracks, sliceTracks,
2919  allTracksBirks, sliceTracksBirks,
2920  &rec.mc.allcosmics);
2921 
2922  // Add faveNu to cosmic
2923  AddCosmicTruthToVec(faveNu, allTracks, sliceTracks, allTracksBirks,
2924  sliceTracksBirks, &rec.mc.cosmic);
2925  }
2926  // Otherwise we are a neutrino
2927  else {
2928  // Add faveNu to allnus
2929  if (fParams.FullTruthInfo())
2930  AddMCTruthToVec(faveNu, allTracks, sliceTracks, allTracksBirks,
2931  sliceTracksBirks, &rec.mc.allnus, evt, spill);
2932 
2933  // Add faveNu to nu
2934  AddMCTruthToVec(faveNu, allTracks, sliceTracks, allTracksBirks,
2935  sliceTracksBirks, &rec.mc.nu, evt, spill);
2936  rec.mc.faveidxenergy = currentNuIdx;
2937 
2939  pv.push_back(faveNu->neutrinoInt);
2940 
2941  if (rec.mc.nu.size() > 0) {
2942  // if a neutrino is matched to the slice, fill the mc reweighting
2943  // object.
2944 
2946  pv, evt, fParams.ReweightLabel());
2948  if (GetAssociatedProduct(fmRW, 0, rw))
2949  rec.mc.nu[0].rwgt.genie = GenieReweightTable(rw);
2950 
2952  pv, evt, fParams.FluxReweightLabel());
2953  fxwgt::FluxWeights flxrw;
2954  if (GetAssociatedProduct(fmFLXRW, 0, flxrw))
2955  rec.mc.nu[0].rwgt.ppfx = FluxReweights(flxrw);
2956  else
2957  rec.mc.nu[0].rwgt.ppfx.setDefault();
2958  }
2959 
2960  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
2961  faveIdsEff[sliceId])
2962  rec.mc.faveidxeff = currentNuIdx;
2963  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
2964  faveIdsPur[sliceId])
2965  rec.mc.faveidxpur = currentNuIdx;
2966  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
2967  faveIdsEffPur[sliceId])
2968  rec.mc.faveidxeffpur = currentNuIdx;
2969  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
2970  faveIdsEffThenPur[sliceId])
2971  rec.mc.faveidxeffthenpur = currentNuIdx;
2972 
2973  currentNuIdx++;
2974  }
2975  }
2976  }
2977 
2978  for (unsigned int iEffPur = 0; iEffPur < sliceTruthTable[sliceId].size();
2979  ++iEffPur) {
2980  if (sliceTruthTable[sliceId][iEffPur].neutrinoInt &&
2981  sliceTruthTable[sliceId][iEffPur].neutrinoInt->NeutrinoSet())
2982  nNuGlobal++;
2983  // continue if truth left no hits in slice
2984  if (!(sliceTruthTable[sliceId][iEffPur].nSliceHits > 0)) continue;
2985 
2986  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex !=
2987  faveIdsEnergy[sliceId]) {
2988  // Cosmics
2989  if (!sliceTruthTable[sliceId][iEffPur].neutrinoInt->NeutrinoSet()) {
2990  if (fParams.FullTruthInfo())
2991  AddCosmicTruthToVec(&sliceTruthTable[sliceId][iEffPur], allTracks,
2992  sliceTracks, allTracksBirks, sliceTracksBirks,
2993  &rec.mc.allcosmics);
2994  }
2995  // Neutrino
2996  else {
2997  if (fParams.FullTruthInfo())
2998  AddMCTruthToVec(&sliceTruthTable[sliceId][iEffPur], allTracks,
2999  sliceTracks, allTracksBirks, sliceTracksBirks,
3000  &rec.mc.allnus, evt, spill);
3001 
3002  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
3003  faveIdsEff[sliceId])
3004  rec.mc.faveidxeff = currentNuIdx;
3005  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
3006  faveIdsPur[sliceId])
3007  rec.mc.faveidxpur = currentNuIdx;
3008  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
3009  faveIdsEffPur[sliceId])
3010  rec.mc.faveidxeffpur = currentNuIdx;
3011  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
3012  faveIdsEffThenPur[sliceId])
3013  rec.mc.faveidxeffthenpur = currentNuIdx;
3014 
3015  currentNuIdx++;
3016  }
3017  }
3018  } // end loop to fill all* vectors
3019 
3020  // fill size fields for mc branch
3021  rec.mc.nnu = rec.mc.nu.size();
3022  rec.mc.ncosmic = rec.mc.cosmic.size();
3023  if (fParams.FullTruthInfo()) {
3024  rec.mc.nallnus = rec.mc.allnus.size();
3025  rec.mc.nallcosmics = rec.mc.allcosmics.size();
3026  }
3027  rec.mc.global.nnu = nNuGlobal;
3028  } // end MC block
3029 
3031 
3032  fRecTree->Fill();
3033  srcol->push_back(rec);
3034  util::CreateAssn(*this, evt, *srcol, art::Ptr<rb::Cluster>(slices, sliceId),
3035  *srAssn);
3036  } // end loop over slices
3037 
3038  *spillcol = spill;
3039 
3040  evt.put(std::move(srcol));
3041  evt.put(std::move(spillcol));
3042  if (fParams.FillNuTree()) evt.put(std::move(spilltruthcol));
3043  evt.put(std::move(srAssn));
3044 }
3045 
3048 
3049  // If this is MC get the spill info from the generator
3050  if (!fIsRealData) {
3052  }
3053  // If we have real data try to get the spill info from the NuMI beam
3054  else {
3056  }
3057 
3058  if (!pot.failedToGet()) {
3059  // NOTE: The POT that comes out of IFDBSpillInfo has an implied
3060  // factor of 1E12 that is being left off.
3061  // For now the data POT is being multiplied by this factor so that
3062  // the number is consistent with what we get in MC.
3063  // TODO: Is this what we want to do in the future?
3064  float subrunpot = pot->totgoodpot;
3065  if (fIsRealData) subrunpot *= 1e12;
3066  fTotalPOT += subrunpot;
3067 
3068  mf::LogInfo("CAFMaker")
3069  << "POT in subrun " << sr.subRun() << " = " << pot->totgoodpot;
3070  }
3071 
3073  sr.getByLabel(fParams.ExposureLabel(), expo);
3074 
3075  if (!expo.failedToGet()) {
3076  fTotalLivetime += expo->totlivetime;
3077  mf::LogInfo("CAFMaker")
3078  << "Exposure in subrun " << sr.subRun() << " = " << expo->totlivetime;
3079  }
3080 }
3081 
3082 //......................................................................
3084  bool EmptyFile = false;
3085  if (fTotalEvents == 0) {
3086  // Among other things this makes us unable to determine the file_type
3087  // metdata field.
3088  if (fParams.AbortOnEmpty()) {
3089  mf::LogError("CAFMaker") << "No events processed in this file. Aborting rather than "
3090  << "produce an empty CAF."
3091  << std::endl;
3092  abort();
3093  }
3094  else {
3095  mf::LogWarning("CAFMaker") << "Making an empty CAF file." << std::endl;
3096  EmptyFile = true;
3097  }
3098  }
3099 
3100  // Make sure the recTree is safely in the file before starting on the POT
3101  // etc histograms. This way, absence of the histograms should be a clear
3102  // sign something went wrong with the tree.
3103  fFile->Write();
3104 
3105  fFile->cd();
3106  hPOT->Fill(.5, fTotalPOT);
3107  hEvents->Fill(.5, fTotalEvents);
3108  hLivetime->Fill(.5, fTotalLivetime);
3109 
3110  hPOT->Write();
3111  hEvents->Write();
3112  hLivetime->Write();
3113  fFile->Write();
3114 
3115  if (fParams.IsSinglesOverlay()) {
3118  hSinglePOT->Fill(0.5, fTotalSinglePOT);
3119 
3120  hSinglePOT->Write();
3121  hTotalTrueNonswapNue->Write();
3122  hTotalTrueSingleNue->Write();
3123  }
3124 
3125  std::map<std::string, std::string> metadata =
3127  metadata["data_tier"] = fParams.DataTier();
3128  if (EmptyFile)
3129  metadata["file_type"] = "empty";
3130  else if (fIsRealData)
3131  metadata["file_type"] = "importedDetector";
3132  else
3133  metadata["file_type"] = "importedSimulated";
3134 
3135  // Gymnastics to mark "restrictedbox" in nova.special for FD numi data
3136  if (metadata.count("online.stream")) {
3137  // Then this is a real data file and the stream exists in the metadata
3138  if (!metadata["online.stream"].compare("0") &&
3140  // Then this is a numi trigger file from the far detector
3141  if (!fParams.EnableBlindness()) {
3142  // Then we are not blinding it, need to mark it as such in the metadata
3143  std::string tag("restrictedbox");
3144 
3145  if (metadata.count("nova.special") &&
3146  std::string("").compare(metadata["nova.special"]) &&
3147  std::string("none").compare(
3148  metadata["nova.special"])) { // Then the metadata is not the
3149  // same as none, we need to append
3150  // to it
3151  // Use a stringstream to jam the strings together and append the tag
3152 
3153  std::stringstream ss;
3154  ss.str();
3155  ss << metadata["nova.special"] << "-" << tag;
3156  metadata["nova.special"] = ss.str();
3157  } else
3158  metadata["nova.special"] = tag;
3159  }
3160  }
3161  }
3162 
3163  AddMetadataToFile(metadata);
3164 }
3165 
3166 //......................................................................
3168  const simb::MCNeutrino& nu) const {
3169  if (nu.CCNC() == 1 && nu.Nu().PdgCode() == flux.fntype) {
3170  // NC
3171  return 1;
3172  }
3173  if (nu.CCNC() == 1 && nu.Nu().PdgCode() != flux.fntype) {
3174  // Swapped NC. Doesn't exist. If people want to play clever
3175  // reweighting games to reclaim these, they can figure that out,
3176  // but this allows the naive implementation to work.
3177  return 0;
3178  }
3179 
3181  return calc.P(flux.fntype, nu.Nu().PdgCode(), nu.Nu().E());
3182 }
3183 
3184 //......................................................................
3189  art::Handle<std::vector<rawdata::RawTrigger>> trigs, SRSpill& spill,
3190  art::Event& evt) const {
3191  // Add in run, subrun, evt information so that the spill tree can be
3192  // mapped to the record tree.
3193  spill.run = evt.run();
3194  spill.subrun = evt.subRun();
3195  spill.evt = evt.id().event();
3196 
3197  // If we fail to identify a spill, we want it set to bad
3198 
3199  // sumdata::SpillData variables
3200  spill.det = fDet;
3201  spill.ismc = !fIsRealData;
3202  spill.isgoodspill = false;
3203  spill.spilltimesec = 0;
3204  spill.spilltimensec = 0;
3205  spill.gpsspilltimesec = 0;
3206  spill.gpsspilltimensec = 0;
3207  spill.deltaspilltimensec = 0;
3208  spill.spillpot = 0.0;
3209  spill.livetime = 0.0;
3210  spill.hornI = 0.0;
3211  spill.posx = 0.0;
3212  spill.posy = 0.0;
3213  spill.widthx = 0.0;
3214  spill.widthy = 0.0;
3215 
3216  // Calculate diblock masking. Also in header, but needed here for
3217  // mass accounting.
3218  unsigned int maskstatus = 0;
3219  unsigned int dibmask = fMask;
3220 
3221  if (rh->GoodDiBlockMask() == ((1 << 16) - 1))
3222  maskstatus = 0;
3223  else if (rh->GoodDiBlockMask() == ((1 << 17) - 1))
3224  maskstatus = 2;
3225  else
3226  maskstatus = 1;
3227 
3228  unsigned int dibfirst = 0;
3229  unsigned int diblast = 0;
3230  if (dibmask) {
3231  int iD;
3232  iD = 0;
3233  while (!((dibmask >> iD) & 1)) iD++;
3234  dibfirst = iD + 1;
3235  iD = 0;
3236  while (dibmask >> iD) iD++;
3237  diblast = iD;
3238  }
3239 
3240  spill.dibfirst = dibfirst;
3241  spill.diblast = diblast;
3242  spill.dibmask = dibmask;
3243  spill.maskstatus = maskstatus;
3244 
3245  // sumdata::EventQuality vars, filled by DQSpillFlags
3246  spill.nmissingdcms = 0;
3247  spill.fracdcm3hits = 0.0;
3248  spill.nouttimehits = 0;
3249  spill.nnoisyapds = 0;
3250 
3251  spill.nmicroslices = 0;
3252  spill.ndcms = 0;
3253 
3254  spill.nmissingdcmslg = 0; // LiveGeometry check
3255  spill.nbaddcmslg = 0; // LiveGeometry check
3256  spill.dcmedgematchfrac = 0;
3257 
3258  // Need to flag each spill if the pot was good or bad
3259  // each slice will get the same tag
3260 
3261  // If there was no spill info for this slice mark it bad. log message this
3262  // should not happen
3263  if (spillPot.failedToGet()) {
3264  mf::LogError("CAFMaker") << "Spill Data not found for real data event";
3265  spill.isgoodspill = false;
3266  }
3267  // If we have a SpillData object, get all the variables
3268  else {
3269  spill.isgoodspill = spillPot->goodbeam;
3270  spill.spilltimesec = spillPot->spilltimesec;
3271  spill.spilltimensec = spillPot->spilltimensec;
3272  spill.gpsspilltimesec = spillPot->gpsspilltimensec;
3273  spill.gpsspilltimensec = spillPot->gpsspilltimensec;
3274  spill.deltaspilltimensec = spillPot->deltaspilltimensec;
3275  spill.spillpot = spillPot->spillpot;
3276  spill.hornI = spillPot->hornI;
3277  spill.isRHC = spillPot->isRHC;
3278  spill.is0HC = spillPot->is0HC;
3279  spill.isFHC = !spill.isRHC && !spill.is0HC;
3280  spill.widthx = spillPot->widthx;
3281  spill.widthy = spillPot->widthy;
3282  spill.posx = spillPot->posx;
3283  spill.posy = spillPot->posy;
3284 
3285  spill.intx.resize(spillPot->intx.size());
3286  spill.inty.resize(spillPot->inty.size());
3287  spill.bposx.resize(spillPot->bposx.size());
3288  spill.bposy.resize(spillPot->bposy.size());
3289  for (size_t i = 0; i < spillPot->intx.size(); ++i)
3290  spill.intx[i] = spillPot->intx[i];
3291  for (size_t i = 0; i < spillPot->inty.size(); ++i)
3292  spill.inty[i] = spillPot->inty[i];
3293  for (size_t i = 0; i < spillPot->bposx.size(); ++i)
3294  spill.bposx[i] = spillPot->bposx[i];
3295  for (size_t i = 0; i < spillPot->bposy.size(); ++i)
3296  spill.bposy[i] = spillPot->bposy[i];
3297 
3298  if (!fIsRealData) {
3299  // MC will flag each spill as good since it doesn't make bad spills
3300  spill.isgoodspill = true;
3301  } else {
3302  // NOTE: there is an implied factor of 1e12 when getting
3303  // IFDBSpillInfo POT. Adding this back in to avoid confusion
3304  // between data and MC.
3305  spill.spillpot *= 1e12;
3306  }
3307  }
3308 
3309  // If there was no spill info for this slice mark it bad.
3310  // Log message this should not happen
3311  if (spillQual.failedToGet()) {
3312  mf::LogError("CAFMaker")
3313  << "Spill EventQuality not found for real data event";
3314  } else {
3315  spill.nmissingdcmslg = spillQual->nmissingdcmslg;
3316  spill.nmissingdcms = spillQual->nmissingdcms;
3317  spill.fracdcm3hits = spillQual->fracdcm3hits;
3318  spill.nouttimehits = spillQual->nouttimehits;
3319  spill.nnoisyapds = spillQual->nnoisyapds;
3320  spill.dcmedgematchfrac = spillQual->dcmedgematchfrac;
3321  spill.nmicroslices = spillQual->nmicroslices;
3322  }
3323 
3324  if (trigs.failedToGet()) {
3325  mf::LogError("CAFMaker") << "No trigger info found for this event";
3326  } else {
3327  bool once = true;
3328  for (const rawdata::RawTrigger& trig : *trigs) {
3329  // Units of trigger field are 500ns
3330  spill.livetime += trig.fTriggerRange_TriggerLength * 500e-9;
3331  if (once) {
3332  once = false;
3333  spill.trigger = trig.fTriggerMask_TriggerType;
3334  } else {
3335  if (spill.trigger != trig.fTriggerMask_TriggerType) {
3336  mf::LogError("CAFMaker") << "Multiple trigger types found in one "
3337  "Event. Storing the first value seen.";
3338  }
3339  }
3340  }
3341  }
3342 
3343  // dq::DAQEventSummary variables
3344  if (daq.failedToGet()) {
3345  mf::LogError("CAFMaker") << "No DAQ summary info found for this event.";
3346  } else {
3347  spill.ndiblocks = daq->fNDiblocks;
3348  spill.eventincomplete = daq->fEventIncomplete;
3349  spill.emptydatablock = daq->fNemptyDataBlock;
3350  spill.nmicroblocks = daq->fNmicroBlocks;
3351  spill.nemptymicroslice = daq->fNemptyMicroSlice;
3355  spill.nnanoslices = daq->fNtotalNanoSlices;
3364  }
3365 
3366  spill.nbaddcmslg = livegeom->NBadDCMBC();
3367  spill.ndcms = rh->NDCMs();
3368 }
3369 
3371 
3372 } // end namespace caf
3373 ////////////////////////////////////////////////////////////////////////
std::vector< float > intx
Definition: SRSpill.h:47
art::ServiceHandle< geo::LiveGeometry > livegeom
int cosfwdcell
cosmictrack projected # cells from end point forwards to det edge
Definition: SRContain.h:30
float x
Bjorken x = (k-k&#39;)^2 / (2*p.q), [Dimensionless].
Definition: SRNeutrino.h:75
double E(const int i=0) const
Definition: MCParticle.h:232
Information about the neutrino production. Docs from http://www.hep.utexas.edu/~zarko/wwwgnumi/v19/v1...
Definition: SRBeam.h:14
Atom< string > TrueEnergyLabel
float visEinslcNeutron
Sum of FLS hits that made CellHits from this neutrino in this subevent [GeV] that were daughters of n...
Definition: SRNeutrino.h:30
caf::generator_ CAFGeneratorEnum(simb::Generator_t simbGeneratorEnum)
Definition: FillTruth.cxx:786
float time
Time [ns].
Definition: SRElastic.h:25
Det_t det
Detector, ND = 1, FD = 2, NDOS = 3.
Definition: SRHeader.h:28
Near Detector underground.
Definition: SREnums.h:10
Atom< string > NumuEnergyLabel
Atom< string > DiscreteTrackLabel
SRShowerLID shwlid
Shower information.
Definition: SRFuzzyKProng.h:18
float closestslicemindist
minimum distance to the closest slice in time domain
Definition: SRSlice.h:56
T max(const caf::Proxy< T > &a, T b)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
void split(double tt, double *fr)
float visEinslcNeutronBirks
Sum of FLS hits that made CellHits from this neutrino in this subevent [GeV] that were daughters of n...
Definition: SRNeutrino.h:34
Det_t
Which NOvA detector?
Definition: SREnums.h:7
virtual void setDefault()
Definition: SRNueEnergy.cxx:21
int frun
Definition: MCFlux.h:35
const XML_Char * name
Definition: expat.h:151
unsigned long int spilltimesec
Spill time in seconds [s].
Definition: SRSpill.h:33
unsigned int nmissingdcmslg
# of DCMS with 63 or more bad FEBs (LiveGeometry, subset of baddcmslg)
Definition: SRSpill.h:68
bool isRHC
is the beam in antineutrino mode, aka RHC
Definition: SpillData.h:28
SRSliceLID slicelid
Output of SliceLID classifier.
Definition: SRIDBranch.h:56
SRBpfTrack muon
The track reconstructed under the muon assumption.
Definition: SRBpf.h:14
double TrueNeutrinoDistance(novadaq::cnv::DetId det, const SRNeutrino &nu)
Definition: FillTruth.cxx:724
Atom< string > FileExtension
SRVector3D vtx
Vertex position in detector coordinates [cm].
Definition: SRNeutrino.h:45
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
float visE
Sum of FLS hits that made CellHits from this neutrino [GeV].
Definition: SRNeutrino.h:27
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
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.
static MetadataManager & getInstance()
art::ServiceHandle< cheat::BackTracker > bt
Atom< InputTag > FuzzyK2DLabel
bool is0HC
horn off
Definition: SpillData.h:29
SRPresel rockpre
Official rock preselection information.
Definition: SRIDBranch.h:43
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
std::vector< float > bposx
Definition: SRSpill.h:49
Atom< bool > StrictMode
int cycle
MC repition index on a given run,subrun.
Definition: SRNeutrino.h:94
SRBpf bpf
Container class for BreakPointFitter tracks.
Definition: SRFuzzyKProng.h:19
short ncosmic
Number of cosmics in cosmic vector (0 or 1)
Definition: SRTruthBranch.h:38
float annecos
e/cosmic ann output with energy for the slice, currently the same as the most energetic shower ...
Definition: SRELid.h:28
double QSqr() const
Definition: MCNeutrino.h:157
float shwE
Energy of shower [GeV].
Definition: SRSLidEnergy.h:21
double HalfD() const
Definition: CellGeo.cxx:205
Atom< string > ClusterLabel
Variables describing Michel E&#39;s found around the end of a track.
Definition: SRPixelMap.h:13
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
Far Detector at Ash River.
Definition: SREnums.h:11
float visENeutronBirks
Sum of FLS hits that made CellHits from this neutrino [GeV] that were daughters of neutrons with birk...
Definition: SRNeutrino.h:33
void AddMCTruthToVec(const cheat::NeutrinoEffPur *effPur, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, std::vector< SRNeutrino > *vec, const art::Event &evt, const SRSpill &spill)
int nouttimehits
of out-beam-window hits
Definition: EventQuality.h:25
Atom< InputTag > OverlapEKalLabel
Atom< string > DiFShowerCVNLabelOldPresel
double fpdpx
Definition: MCFlux.h:55
float unixtime
unix time of spill
Definition: SRHeader.h:45
unsigned short maskstatus
0 no mask found in DB, 1 mask used ok, 2 masking turned off. If 0 or 2 dibmask is instead the configu...
Definition: SRHeader.h:35
SRCVNResult cvnnocosmics
Output from CVN - No cosmics ued in training (many-class PID)
Definition: SRIDBranch.h:54
static bool sortTrackLength(const SRTrack &a, const SRTrack &b)
std::vector< double > bposy
Definition: SpillData.h:33
bool blind
if true, record has been corrupted for blindness
Definition: SRHeader.h:29
float visEinslc
Sum of FLS hits that made CellHits from this neutrino in this subevent [GeV].
Definition: SRNeutrino.h:28
SRHeader hdr
Header branch: run, subrun, etc.
const Var kDistAllBottom([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngbottom)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngbottom);})
Distance of all showers in slice from the bottom edge of detector.
Definition: NueVars.h:33
float orphCalE
calorimetric energy of hits that don&#39;t appear in any FuzzyK prongs
Definition: SRFuzzyK.h:30
Atom< string > NueSelLabel
art::ServiceHandle< geo::Geometry > geom
short faveidxeffthenpur
Index of favorite in allnus when neutrinos are sorted by efficiency and slices break ties by purity...
Definition: SRTruthBranch.h:35
set< int >::iterator it
unsigned int subrun
subrun number
Definition: SRHeader.h:22
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
bool is0HC
Definition: SRSpill.h:44
double ftpx
Definition: MCFlux.h:79
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
unsigned int npizero
Number of &#39;s after neutrino reaction, before FSI.
Definition: SRNeutrino.h:62
void InitializeOutfile()
SRBeam beam
Information about neutrino production.
Definition: SRNeutrino.h:89
double Weight() const
Definition: MCParticle.h:253
float distallpngbottom
Definition: SRNueCosRej.h:95
SubRunNumber_t subRun() const
Definition: SubRun.h:44
This is a helper class for ParticleIDAlg that provides a tidy structure in which to hold the dE/dx hi...
int NAnalysisChannels(int sr=0)
Definition: RunHistory.cxx:285
Neutrino energy output from RegCVN.
Definition: RegResult.h:29
Variables describing Michel E&#39;s found around the end of a track.
Definition: SRSliceMap.h:13
Atom< string > RockPreselLabel
void AddParticleToVec(const sim::Particle &part, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, std::vector< SRTrueParticle > *vec, const std::vector< sim::TrueEnergy > &TrueEnergies)
Definition: FillTruth.cxx:255
const Var kDistAllWest([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngwest)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngwest);})
Distance of all showers in slice from the west edge of detector.
Definition: NueVars.h:36
std::map< MixingType, std::string > mapMixingTypeString
void FillVetoVars(const presel::Veto &veto, const presel::Veto &nueveto, caf::SRVeto &srveto)
Definition: FillPIDs.cxx:520
std::vector< SRProng > png2d
Vector of 2D prong objects.
Definition: SRFuzzyK.h:20
unsigned int lastplane
last plane
Definition: SRSlice.h:27
SRCosRej cosrej
Output from CosRej (Cosmic Rejection)
Definition: SRIDBranch.h:45
int fppmedium
Definition: MCFlux.h:62
unsigned int firstplane
first plane
Definition: SRSlice.h:26
Atom< string > CosmicCVNLabel
Atom< bool > UseGeVPixelMaps
void FillNumuLSTMEnergyVars(const LSTME::LSTMEnergy &artE, caf::SRNumuEnergy &cafE)
size_t ntracks
Definition: SRKalman.h:23
double fgen2vtx
distance from ray origin to event vtx
Definition: MCFlux.h:104
float q2
Squared momentum transfer [GeV^2].
Definition: SRNeutrino.h:74
Atom< string > CVNParticleLabel
float closestsliceminfromback
minimum distance to edge of detector in the closest slice
Definition: SRSlice.h:53
A vertex found by the VertexDT algorithm.
Definition: SRVertexDT.h:13
Atom< string > NueVetoLabel
X or Y views.
Definition: PlaneGeo.h:30
const Var kDistAllTop([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngtop)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngtop);})
Distance of all showers in slice from the top edge of detector.
Definition: NueVars.h:30
std::vector< SRVertexDT > vdt
Vector of vertices found by VertexDT.
virtual ~CAFMaker()
unsigned int idxremid
index number of the best ReMId track
Definition: SRKalman.h:32
float E
Energy [GeV].
Definition: SRSLidEnergy.h:18
Atom< string > VetoLabel
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
int nanoslicebufferempty
# of nanoslices reporting BufferEmpty
Definition: SRSpill.h:91
unsigned short Plane() const
Definition: CellHit.h:39
Atom< string > KalmanTrackLabel
double SimpleOscProb(const simb::MCFlux &flux, const simb::MCNeutrino &nu) const
Atom< bool > FullTruthInfo
float ncVal
Cosmic CVN nc score for each time slice.
int inttype
Interaction type enum int_type::[...].
Definition: SRNeutrino.h:58
unsigned short day
day of spill within month
Definition: SRHeader.h:39
float woscdumb
Simplest possible oscillation weight.
Definition: SRNeutrino.h:54
std::vector< SRCosmicCVN > cosmiccvn
Contain cosmic CVN scores for all time windows in event.
Definition: SRSpill.h:100
Atom< string > BreakPointTrackLabel
int Mother() const
Definition: MCParticle.h:212
void fillSizes()
Definition: SRKalman.cxx:14
std::vector< SRTrkME > trkcosmic
Definition: SRMichelE.h:24
void GetByLabelStrict(const art::Event &evt, const std::string &label, art::Handle< T > &handle) const
std::vector< SRFuzzyKProng > png
Vector of 3D prong objects.
Definition: SRFuzzyK.h:19
float regcvnhadE
Hadronic Energy predicted by Regression CNN [GeV].
Definition: SRNumuEnergy.h:58
geo::View_t View() const
Definition: CellHit.h:41
Module to create Common Analysis Files from ART files.
int ftptype
Definition: MCFlux.h:82
SRMichelE me
Michel electron branch.
Atom< string > DiFShowerLIDLabel
std::map< std::string, MixingType > mapStringMixingType
A reconstructed shower from the MRProperties module.
void FillTrackVars(const rb::Track &trk, SRTrack &srTrk, const std::vector< rb::Cluster > &sliceList, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, int sliceIdx)
Definition: FillReco.cxx:278
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const Var kDistAllBack([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngback)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngback);})
Distance of all showers in slice from the back edge of detector.
Definition: NueVars.h:42
bool DetFineTimingSetting() const
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
Atom< string > SlcLIDLabel
std::vector< NeutrinoWithIndex > allMCTruth() const
void FillTrackContainmentVars(const cosrej::TrkCntObj &trkcnt, SRTrack &srTrk)
Function to fill SRTrack containment information.
Definition: FillReco.cxx:552
void GetByLabelIfExists(const art::Event &evt, const std::string &label, art::Handle< T > &handle) const
std::vector< SRCosmic > allcosmics
vector holding all Cosmics
Definition: SRTruthBranch.h:29
T sqrt(T number)
Definition: d0nt_math.hpp:156
float len
Track length of identified track.
Definition: SRMuId.h:23
std::vector< SRTrueParticle > prim
Primary daughters, lepton comes first in vector.
Definition: SRNeutrino.h:81
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:35
float Y() const
Definition: SRVector3D.h:33
void FillCVNFeaturesVars(const cvn::Features &features, caf::SRCVNFeatures &srcvnfeatures, int maxcomponents)
Definition: FillPIDs.cxx:663
Atom< bool > ApplyLEMNuePresel
float Output() const
Definition: RegResult.h:22
bool isFHC
Flags for horn direction.
Definition: SRSpill.h:43
unsigned int npiplus
Number of &#39;s after neutrino reaction, before FSI.
Definition: SRNeutrino.h:60
int fNdataBlockMissingData
of occurances of isMissingData
float cosmicVal
Cosmic CVN cosmic score for each time slice.
int nanosliceadcerror
# of nanoslices reporting ADCError
Definition: SRSpill.h:96
std::vector< SRTrainingData > trainingdata
Collection of labels associated with the Pixel maps.
SRTrainingBranch training
Extra training information for prototyping PIDs etc.
SRMCReweight rwgt
Definition: SRNeutrino.h:87
std::vector< float > inty
Definition: SRSpill.h:48
std::vector< unsigned int > genVersion
Version of the generator that created this neutrino interaction.
Definition: SRNeutrino.h:84
int fNumNeutron
number of neutrons after reaction, before FSI
Definition: GTruth.h:76
double fvx
Definition: MCFlux.h:52
float kalyposattrans
Y position of Kalman track and transition (ND only, use to check if went through air gap) ...
Definition: SRContain.h:50
int fNmicroBlocks
How many microblocks?
std::string const & fileName() const
Definition: FileBlock.h:38
short faveidxpur
Index of favorite in allnus when sorted by purity.
Definition: SRTruthBranch.h:33
int ndiblocks
# of diblocks reporting in event
Definition: SRSpill.h:81
Results for Regression CVN.
Float_t ss
Definition: plot.C:24
int HitNuc() const
Definition: MCNeutrino.h:152
short faveidxeff
Index of favorite in allnus when sorted by slicer efficiency.
Definition: SRTruthBranch.h:31
void produces(std::string const &instanceName={}, Persistable const persistable=Persistable::Yes)
virtual void endSubRun(art::SubRun &sr)
double HalfW() const
Definition: CellGeo.cxx:191
Truth info for all neutrinos in the spill.
A collection of associated CellHits.
Definition: Cluster.h:47
Atom< string > CVNTrainingDataLabel
double DetLength() const
Atom< string > NuMIBeamLabel
float time
Time [ns].
Definition: SRVertexDT.h:22
static bool sortProngCalE(const SRFuzzyKProng &a, const SRFuzzyKProng &b)
virtual void setDefault()
Definition: SRCVNResult.cxx:28
std::pair< int, int > calcFirstLastLivePlane(int plane, std::bitset< 14 > binary, caf::Det_t det=caf::kFARDET)
Atom< InputTag > FuzzyKWeight3DLabel
SRVector3D vtx
Vertex position in detector coordinates. [cm].
Definition: SRVertexDT.h:23
std::vector< float > bposy
Definition: SRSpill.h:50
Atom< string > ReweightLabel
Atom< bool > IsSinglesOverlay
bool finetiming
Is fine timing enabled in this run?
Definition: SRHeader.h:55
OStream cerr
Definition: OStream.cxx:7
int fNanoSliceADCError
of nanoslices reporting ADCError
unsigned int run
run number
Definition: SRHeader.h:21
SRTrackBase discrete
3D tracks produced by DiscreteTrack
Definition: SRTrackBranch.h:25
double fXsec
cross section of interaction
Definition: GTruth.h:30
int fNanoSliceCommError
of nanoslices reporting CommError
int nTotalHits
Total number of hits the neutrino left in the detector.
Definition: BackTracker.h:54
Atom< bool > FillNuTree
Atom< bool > FillPixelMaps
std::map< std::string, std::string > & GetMetadata()
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
int nemptymicroslice
# of empty micro slices
Definition: SRSpill.h:84
caf::StandardRecord * rec
Definition: tutCAFMacro.C:20
float visENeutron
Sum of FLS hits that made CellHits from this neutrino [GeV] that were daughters of neutrons...
Definition: SRNeutrino.h:29
Atom< string > SliceLIDLabel
void abs(TH1 *hist)
int fNumPiPlus
number of pi pluses after reaction, before FSI
Definition: GTruth.h:78
unsigned int longestidx
index of longest prong
Definition: SRFuzzyK.h:23
unsigned int nouttimehits
# of out-of-time hits
Definition: SRSpill.h:66
double energyTotal
Sum of FLS hits from the neutrino contributing to any hit in the event.
Definition: BackTracker.h:52
bool numucontainSA
is this contained by second analysis Numu Standards?
Definition: SRContain.h:19
int KalBakCellND() const
Definition: CosRejObj.cxx:670
Atom< string > DiscreteTrack2dLabel
Atom< bool > FullTrainingInfo
float hadE
Hadronic energy [GeV].
Definition: SRSLidEnergy.h:22
SRNueEnergy nue
Nue energy variables.
Atom< string > DiFShowerCVNLabelNoCosmics
Atom< string > SpillQualLabel
int fNumPiMinus
number of pi minuses after reaction, before FSI
Definition: GTruth.h:79
double fppdxdz
Definition: MCFlux.h:58
Atom< string > GeneratorLabel
TString ip
Definition: loadincs.C:5
double corr(TGraph *g, double thresh)
Atom< bool > UseEnergy
Atom< string > CosRejLabel
unsigned int run
run number
Definition: SRSpill.h:25
SRCVNResult cvnoldpresel
Output from CVN - Preselection used in Prod3/4 (many-class PID)
Definition: SRIDBranch.h:53
double DistToFront(TVector3 vertex)
SRCVNResult cvnloosepreselptp
Output from CVN - Loose Presel plus PtP cut (many-class PID)
Definition: SRIDBranch.h:52
int nanoslicedatanotpresent
# of nanoslices reporting !DataPresent
Definition: SRSpill.h:89
std::vector< SRHoughVertex > hough
Vector of vertices found by HoughVertex.
int tgtA
A of target nucleus.
Definition: SRNeutrino.h:72
int ftgen
Definition: MCFlux.h:83
float closestslicecalE
Calorimetric energy of the closest-in-time slice (GeV)
Definition: SRSlice.h:45
const PlaneGeo * Plane(unsigned int i) const
double DistToBack(TVector3 vertex)
Atom< string > CVNLabelLoosePreselPtP
bool ismc
data or MC? True if MC
Definition: SRHeader.h:27
int emptydatablock
# of empty data blocks
Definition: SRSpill.h:82
int cosbakcell
cosmictrack projected # cells from start point backwards to det edge
Definition: SRContain.h:35
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
float nutauVal
Cosmic CVN nutau score for each time slice.
int nanoslicenolinkstatus
# of nanoslices reporting !LinkPresent
Definition: SRSpill.h:90
Atom< string > ExposureLabel
void setDefault()
Definition: SRELid.cxx:27
double fdk2gen
distance from decay to ray origin
Definition: MCFlux.h:103
DEFINE_ART_MODULE(TestTMapFile)
std::vector< SRTrkME > trkkalman
Definition: SRMichelE.h:22
virtual double P(int flavBefore, int flavAfter, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
Definition: OscCalcDumb.h:21
Det_t det
Detector, ND = 1, FD = 2, NDOS = 3.
Definition: SRSpill.h:29
Atom< string > MRParentSliceLabel
Features, basic output of CVN neural net.
Definition: CVNFeatures.h:15
float subevtstarttime
time of beginning of slice within spill [ns]
Definition: SRHeader.h:47
std::string fCafFilename
PixelMap for CVN.
unsigned int evt
ART event number, indexes trigger windows.
Definition: SRNeutrino.h:93
SRRvp rvp
Output from RecoVariablePID (RVP)
Definition: SRIDBranch.h:44
Defines an enumeration for prong classification.
Particle class.
unsigned int nproton
Number of protons after neutrino reaction, before FSI.
Definition: SRNeutrino.h:63
osc::OscCalcDumb calc
virtual void setDefault()
Definition: SRRemid.cxx:29
unsigned int nmissingdcms
# of missing DCMs
Definition: SRSpill.h:64
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
float closestsliceminfromeast
minimum distance to edge of detector in the closest slice
Definition: SRSlice.h:54
std::string genConfigString
String associated with generator configuration. (For GENIE 3, this is the "Comprehensive Model Config...
Definition: SRNeutrino.h:85
int KalFwdCell() const
Definition: CosRejObj.cxx:652
SRFluxWeights FluxReweights(const fxwgt::FluxWeights &flxwgts)
Definition: FillTruth.cxx:776
int nmissingdcms
of missing DCMs
Definition: EventQuality.h:23
float posy
y position on target
Definition: SRSpill.h:52
float dcmedgematchfrac
How many hits at the DCM edge are matched in the adjacent DCM?
Definition: SRSpill.h:74
int kalbakcellnd
Kalmantrack projected # cells from start point backwards to det edge, including muon catcher...
Definition: SRContain.h:46
Atom< bool > ApplyingFilter
unsigned short minute
minute of spill
Definition: SRHeader.h:42
double efficiency
Efficiency (based on FLS energy) of neutrino interaction relative to slice.
Definition: BackTracker.h:48
Loaders::FluxType flux
std::vector< float > fOutput
Vector of outputs from neural net.
Definition: Result.h:30
Det_t fDet
Detector ID in caf namespace typedef.
float genweight
Weight, if any, assigned by the generator.
Definition: SRNeutrino.h:79
art::FindManyP< T > FindManyPStrict(const U &from, const art::Event &evt, const art::InputTag &label) const
unsigned int subrun
subrun number
Definition: SRNeutrino.h:92
void FillProngVars(const rb::Prong &prng, T &srPrng, const std::vector< rb::Cluster > &sliceList, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, int sliceIdx)
Definition: FillReco.cxx:176
void FillSlidVars(const slid::ShowerLID &slid, SRShowerLID &shwlid)
Definition: FillPIDs.cxx:168
Atom< InputTag > FuzzyKWeight2DLabel
object containing MC flux information
novadaq::cnv::DetId fDetID
Detector ID in nova daq convention namespace typedef.
bool isseaquark
Did neutrino scatter off a sea quark.
Definition: SRNeutrino.h:67
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
Definition: Run.h:31
void FillSliceLID(const SliceLID::Prediction &artSliceLID, caf::SRSliceLID &cafSliceLID)
Definition: FillPIDs.cxx:686
static bool sortNuWithIdx(const cheat::NeutrinoWithIndex &a, const cheat::NeutrinoWithIndex &b)
int ndatablockmissingdata
# of occurances of isMissingData
Definition: SRSpill.h:86
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
std::vector< int > SliceToOrderedNuIds(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable, std::function< double(const cheat::NeutrinoEffPur &)> slMetric, std::function< double(const cheat::NeutrinoEffPur &)> nuMetric) const
Given a vector of indexed neutrino interaction and a vector of slice truth tables, returns a vector of neutrino interaction indices ordered by best match to the corresponding slice. Here, best match is determined according to the given cheat::NeutrinoEffPur functions for the slice and the nu.
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
The SRNeutrino is a representation of neutrino interaction information.
Definition: SRNeutrino.h:19
Atom< string > SPIDLabel
SRCVNResult cvnnocosmics
Output from CVN - No cosmics ued in training (many-class PID)
unsigned int evt
ART event number, indexes trigger windows.
Definition: SRHeader.h:25
TVector3 MaxXYZ() const
Definition: Cluster.cxx:492
const char * label
float pur
Slicer purity for this truth interaction.
Definition: SRNeutrino.h:39
unsigned int closestslicenhit
Number of hits in the closest-in-time slice.
Definition: SRSlice.h:44
float fracdcm3hits
fraction of DCM3 hits in horizontal modules
Definition: SRSpill.h:65
int nmicroslices
of micro slices
Definition: EventQuality.h:27
TH1D * hTotalTrueSingleNue
Atom< string > DataTier
unsigned short dibfirst
first diblock in detector configuration (1-14)
Definition: SRSpill.h:56
Atom< string > RegCVNLabel
int nnanoslices
# of nano slices in the event
Definition: SRSpill.h:88
Atom< string > CVNPixelMapLabel
const Var kDistAllEast([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngeast)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngeast);})
Distance of all showers in slice from the east edge of detector.
Definition: NueVars.h:39
Atom< string > SingleMixerLabel
float nueVal
Cosmic CVN nue score for each time slice.
float hornI
Horn current.
Definition: SRSpill.h:42
void CopyRemidVars(const SRKalmanTrack &srTrk, caf::SRRemid &remid)
Definition: FillPIDs.cxx:676
PID
Definition: FillPIDs.h:13
std::vector< double > bposx
Definition: SpillData.h:32
unsigned short gain
Global gain setting of the detector.
Definition: SRHeader.h:54
bool IsValid
This defaults to false, and only gets set to true in CAFMaker if there is a valid BPF track...
Definition: SRBpfTrack.h:15
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
void FillNCPi0BkgRejVars(const ncpi0::NCPi0BkgRej &ncpi0bkgrej, caf::SRNCPi0BkgRej &srncpi0bkgrej)
Definition: FillPIDs.cxx:496
TTree * fCosmicCVNTree
fclose(fg1)
float X() const
Definition: SRVector3D.h:32
int ftgptype
Definition: MCFlux.h:84
bool isValid() const
Definition: Handle.h:189
unsigned short Cell() const
Definition: CellHit.h:40
TH1D * hTotalTrueNonswapNue
A class detailing the cuts made on a particular slice.
Definition: Veto.h:12
virtual void setDefault()
int fResNum
resonance number
Definition: GTruth.h:80
unsigned int nbaddcmslg
# of DCMS with too many bad channels (LiveGeometry)
Definition: SRSpill.h:69
int fNDiblocks
of diblocks reporting in event
SRLorentzVector p
True momentum [GeV].
Definition: SRNeutrino.h:44
short nallcosmics
Number of cosmics in allcosmics vector.
Definition: SRTruthBranch.h:40
int fNumProton
number of protons after reaction, before FSI
Definition: GTruth.h:75
Atom< string > DiFShowerShwLabel
double MinTNS() const
Definition: Cluster.cxx:482
void FillCVNProngTrainingData(const cvn::ProngTrainingData &cvnpdata, caf::SRProngTrainingData &srpdata)
Definition: FillPIDs.cxx:641
double DistToTop(TVector3 vertex)
int InteractionType() const
Definition: MCNeutrino.h:150
unsigned int npiminus
Number of &#39;s after neutrino reaction, before FSI.
Definition: SRNeutrino.h:61
Atom< string > JMShowerLabel
int nanoslicepacketerror
# of nanoslices reporting PacketError
Definition: SRSpill.h:94
unsigned int subrun
subrun number
Definition: SRSpill.h:26
TString part[npart]
Definition: Style.C:32
int fNemptyMicroSlice
How many empty micro slices?
Atom< bool > EnableBlindness
virtual double TotalLength() const
Distance along prong to reach last cell hit.
Definition: Prong.cxx:186
SRKalman kalman
Tracks produced by KalmanTrack.
Definition: SRTrackBranch.h:24
Atom< string > NCCosRejLabel
Far Detector at Ash River, MN.
int resnum
Straight from GENIE, resonance number.
Definition: SRNeutrino.h:68
string cmd
Definition: run_hadd.py:52
bool iscc
true if charged-current interaction, false if not.
Definition: SRNeutrino.h:57
Atom< string > DAQEventSummaryLabel
short faveidxeffpur
Index of favorite in allnus when sorted by product of efficiency and purity.
Definition: SRTruthBranch.h:34
int kalfwdcellnd
Kalmantrack projected # cells from end point forwards to det edge, including muon catcher...
Definition: SRContain.h:41
std::vector< double > intx
Definition: SpillData.h:30
unsigned int nshwlid
number of shwlid showers - either 0 or number of 3d prongs
Definition: SRFuzzyK.h:24
void FillLEMVars(const lem::PIDDetails &lem, caf::SRLem &srlem)
Definition: FillPIDs.cxx:282
SRGlobalTruth global
Definition: SRTruthBranch.h:42
void FillNuePreselVars(const presel::PreselObj &nuepre, caf::SRPresel &srnuepre)
Definition: FillPIDs.cxx:505
void FillCosRejVars(const cosrej::CosRejObj &cosrej, caf::SRCosRej &srcosrej, caf::SRContain &srcontain)
Definition: FillPIDs.cxx:360
bool GetPsetParameter(const fhicl::ParameterSet &pset, const std::vector< std::string > &name, T &ret) const
int fNanoSliceNoLinkStatus
of nanoslices reporting !LinkPresent
void FillXnueVars(const xnue::Xnue &xnue, caf::SRXnue &srxnue)
Definition: FillPIDs.cxx:345
float L
True distance from hadron/muon decay to neutrino interaction [km].
Definition: SRNeutrino.h:35
double DistToBottom(TVector3 vertex)
Result for CVN.
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
float ndhadcalcatE
Near detector – hadronic calorimetric energy NOT on the muon track in muon catcher [GeV]...
Definition: SRNumuEnergy.h:44
SRFuzzyK fuzzyk
Primary 3D prong object.
Definition: SRHoughVertex.h:29
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Definition: Cluster.cxx:439
Atom< string > LemLabel
Atom< std::string > CAFFilename
double W() const
Definition: MCNeutrino.h:154
int fptype
Definition: MCFlux.h:63
if(dump)
float distallpngfront
Definition: SRNueCosRej.h:99
double ftvx
Definition: MCFlux.h:76
Atom< string > XnueLabel
Atom< InputTag > FuzzyKOrphLabel
const XML_Char int const XML_Char * value
Definition: expat.h:331
void FindAndAddMichels(std::vector< const sim::Particle * > particles, std::vector< cheat::TrackIDE > &allTracks, std::vector< SRTrueMichelE > *michelVec)
Definition: FillTruth.cxx:661
Atom< InputTag > OverlapEBPFLabel
signed long long int deltaspilltimensec
Definition: SpillData.h:24
float anne
ann output with energy for the slice, currently the same as the most energetic shower ...
Definition: SRELid.h:26
double fracdcm3hits
fraction of DCM3 hits in horizontal modules
Definition: EventQuality.h:24
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
SRCVNResult cvn
Horrible hack to appease CAFAna.
Definition: SRIDBranch.h:51
void FillDiFVars(const rb::Cluster &slice, const std::vector< rb::Cluster > &sliceList, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, int sliceIdx, SRMRProperties &srMR)
Definition: FillReco.cxx:94
float W2
Invariant mass of final state squared. [GeV^2].
Definition: SRNeutrino.h:77
double Y() const
Definition: MCNeutrino.h:156
Atom< string > RawDataLabel
void FillNueCosRejVars(const cosrej::NueCosRej &nuecosrej, caf::SRNueCosRej &srnuecosrej)
Definition: FillPIDs.cxx:424
Atom< string > NuonELabel
Store flux weigths for neutrino correction.
Definition: FluxWeights.h:15
float closestsliceminfrombottom
minimum distance to edge of detector in the closest slice
Definition: SRSlice.h:51
Atom< bool > DataSpillHasMC
int fndecay
Definition: MCFlux.h:50
float dcmedgematchfrac
Low values mean out-of-sync detector.
Definition: EventQuality.h:29
unsigned int idxlongest
Definition: SRKalman.h:33
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
bool numucontain
is this contained by Numu Standards? (can use as general containment)
Definition: SRContain.h:18
void FillRVPVars(const rvp::RVP &rvp, caf::SRRvp &srrvp)
Definition: FillPIDs.cxx:323
virtual void beginRun(art::Run &r)
unsigned short diblast
last diblock in detector configuration (1-14)
Definition: SRSpill.h:57
void FillMRCCParentInfo(const murem::MRCCParent &parent, caf::SRMRCCParent &srparent, const art::Event &evt)
const double a
Atom< string > ReclusShowerLabel
double fpdpz
Definition: MCFlux.h:57
std::string getenv(std::string const &name)
Atom< string > CVNLabelOldPresel
int fNumPi0
number of pi0 after reaction, before FSI
Definition: GTruth.h:77
float time
Time [ns].
Definition: SRHoughVertex.h:25
float eff
Slicer efficiency for this truth interaction.
Definition: SRNeutrino.h:38
SRNueCosRej nuecosrej
Output from NueCosRej (Nue Cosmic Rejection)
Definition: SRIDBranch.h:46
T get(std::string const &key) const
Definition: ParameterSet.h:231
bool isgoodspill
Was the pot for a spill good? (only applicable to data, default true)
Definition: SRSpill.h:32
bool filt
if true, record has ben filtered
Definition: SRHeader.h:30
SRRemid remid
Output from RecoMuonID (ReMId) package.
Definition: SRIDBranch.h:38
SRNCCosRej nccosrej
Output from NCCosRej (NC Cosmic Rejection)
Definition: SRIDBranch.h:47
int fNanoSliceDataNotPresent
of nanoslices reporting !DataPresent
simb::Generator_t generator
event generator that generated this event
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
bool fEventIncomplete
Is the event incomplete?
double widthx
mm
Definition: SpillData.h:36
unsigned long int gpsspilltimensec
Spill time from GPS [ns].
Definition: SRSpill.h:36
unsigned int nmicroslices
# of micro slices
Definition: SRSpill.h:75
int KalBakCell() const
Definition: CosRejObj.cxx:664
Atom< string > CVNLabelNoCosmics
Atom< string > WSCNNLabel
Atom< bool > AbortOnEmpty
short ncosmiccvn
Definition: SRSpill.h:101
#define pot
float PE() const
Definition: CellHit.h:42
unsigned int evt
ART event number, indexes trigger windows.
Definition: SRSpill.h:27
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
A potential interaction point found by the HoughVertex algorithm.
Definition: SRHoughVertex.h:19
virtual void setDefault()
Definition: SRSliceLID.cxx:20
float rnn
ann output for the slice, currently the same as most energetic shower
Definition: SRELid.h:25
bool fIsCharm
did the interaction produce a charmed hadron?
Definition: GTruth.h:71
double fmuparpz
Definition: MCFlux.h:69
double EnergyMetric(const cheat::NeutrinoEffPur &ep)
Function for NeutrinoEffPur&#39;s nu interaction to slice energy.
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
Collect Geo headers and supply basic geometry functions.
bool isRHC
Definition: SRSpill.h:45
void fillSizes()
Definition: SRFuzzyK.cxx:24
Atom< string > HitsLabel
unsigned int run
run number
Definition: SRNeutrino.h:91
void setDefault()
Definition: SRXnue.cxx:28
static bool sortRBProngLength(const art::Ptr< rb::Prong > &a, const art::Ptr< rb::Prong > &b)
unsigned int nbadchan
Number of bad channels in a subrun. Ignores channels in diblocks that are masked off for analysis...
Definition: SRHeader.h:51
Atom< string > HoughVertexLabel
unsigned long int spilltimensec
Definition: SpillData.h:21
float depE
Total energy deposited in shower [GeV].
Definition: SRSLidEnergy.h:20
float subevtendtime
Slice end time [ns].
Definition: SRHeader.h:48
unsigned short doy
day of spill within year
Definition: SRHeader.h:40
float widthy
Spill width in y dimension.
Definition: SRSpill.h:54
Definition: RVP.h:16
Atom< int > NumCVNPrincipalComponents
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
int CosBakCell() const
Definition: CosRejObj.cxx:640
void FillNumuEnergyVars(const numue::NumuE &e, caf::SRNumuEnergy &sre)
void respondToOpenInputFile(const art::FileBlock &fb)
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
void FillCVNParticleResultVars(const std::vector< art::Ptr< rb::PID > > &cvnparts, caf::SRCVNParticleResult &cvnpart)
Definition: FillPIDs.cxx:89
int CosFwdCell() const
Definition: CosRejObj.cxx:628
double X() const
Definition: MCNeutrino.h:155
void setDefault()
Definition: SRPresel.cxx:22
std::vector< SRSlcME > slc
Definition: SRMichelE.h:26
SRCVNResult cvnoldpresel
Output from CVN - Preselection used in Prod3/4 (many-class PID)
int fevtno
Definition: MCFlux.h:36
int nanosliceoverflowerror
# of nanoslices reporting OverflowError
Definition: SRSpill.h:95
double hornI
kA
Definition: SpillData.h:27
unsigned int nnu
WARNING: this variable is duplicated across slices.
Definition: SRGlobalTruth.h:23
SRNuonEResult nuone
Ouput of nuone classifier.
Definition: SRIDBranch.h:57
Eigen::VectorXd vec
int fNanoSliceOverflowError
of nanoslices reporting OverflowError
unsigned short maskstatus
Definition: SRSpill.h:62
static bool sortHoughProng(const SRHoughVertex &a, const SRHoughVertex &b)
SRVector3D vtx
Vertex position in detector coordinates. [cm].
Definition: SRElastic.h:26
virtual void setDefault()
double fppdydz
Definition: MCFlux.h:59
int fNanoSliceBufferFull
of nanoslices reporting BufferFull
static const double fm
Definition: Units.h:83
double ftpz
Definition: MCFlux.h:81
double ftvy
Definition: MCFlux.h:77
double DetHalfHeight() const
unsigned long int gpsspilltimesec
Spill time from GPS [s].
Definition: SRSpill.h:35
float subevtmeantime
Slice mean time [ns].
Definition: SRHeader.h:49
std::vector< unsigned int > DecodeGeneratorVersion(const std::string &versionString)
Definition: FillTruth.cxx:803
int nplanestoback
number of planes between the back of the detector (configuration dependent) and hit with the largest ...
Definition: SRContain.h:24
static bool sortRBTrackLength(const art::Ptr< rb::Track > &a, const art::Ptr< rb::Track > &b)
unsigned short second
second of spill
Definition: SRHeader.h:43
float widthx
Spill width in x dimension.
Definition: SRSpill.h:53
Atom< string > LIDLabel
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
SRSpill spill
Beam spill branch: pot, beam current, etc.
float ann
ann output for the slice, currently the same as most energetic shower
Definition: SRELid.h:24
Vertex location in position and time.
float Z() const
Definition: SRVector3D.h:34
double posy
mm
Definition: SpillData.h:35
Atom< bool > FillTrainingData
Definition: View.py:1
Atom< string > LEMNuePresel
std::vector< SRNeutrino > allnus
vector holding all Neutrinos
Definition: SRTruthBranch.h:28
int fNemptyDataBlock
How many empty data blocks?
void FillDiFShowerVars(const rb::Shower &shw, const std::vector< rb::Cluster > &sliceList, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, int sliceIdx, SRMRProperties &srMR)
Definition: FillReco.cxx:501
void FillTrackInfoVars(const trackinfo::TrackInfoObj &trkinfo, SRTrack &srTrk)
Function to fill additional SRTrack information.
Definition: FillReco.cxx:573
unsigned int nhitslc
Number of hits recorded in this slice by this truth interaction.
Definition: SRNeutrino.h:40
short pdgorig
Unoscillated (unswapped) pdg code.
Definition: SRNeutrino.h:48
void BlindThisRecord(StandardRecord *rec)
Definition: Blinding.cxx:162
int mode
interaction mode from enum mode_type::[QE, RES, COH, ...]
Definition: SRNeutrino.h:56
double MaxTNS() const
Definition: Cluster.cxx:528
Definition: run.py:1
void FillSliceVars(const rb::Cluster &slice, caf::SRSlice &srslice, bool allowEmpty)
Definition: FillReco.cxx:344
Atom< InputTag > ElasticArmsLabel
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
The TrainingData objects contains a PixelMap and the output class type, and any other bit that goes i...
void FillNCCosRejVars(const ncid::NCCosRej &nccosrej, caf::SRNCCosRej &srnccosrej)
Definition: FillPIDs.cxx:487
double fpdpy
Definition: MCFlux.h:56
OStream cout
Definition: OStream.cxx:6
float regcvnEvtE
Regression CVN neutrino energy estimate [GeV].
Definition: SRNueEnergy.h:25
Atom< string > MCSpillDataLabel
void FillShowerVars(const rb::Shower &shw, caf::SRFuzzyKProng &srshw, const std::vector< rb::Cluster > &sliceList, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, const int &sliceIdx)
Definition: FillPIDs.cxx:135
Atom< InputTag > FuzzyK3DLabel
float livetime
Length of readout [s].
Definition: SRSpill.h:40
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
void FillTrackVarsBpfFitSum(const rb::FitSum &fitsum, SRBpfTrack &srBpfTrk)
Definition: FillReco.cxx:313
SRShowerLID lid
LID PID information for the MR Shower.
int tgtZ
Z of target nucleus.
Definition: SRNeutrino.h:71
SRBpfTrack pion
The track reconstructed under the pion assumption.
Definition: SRBpf.h:15
float shwE
reconstructed shower energy [GeV]
Definition: SRShowerLID.h:28
short pdg
pdg code
Definition: SRNeutrino.h:23
void FillMuIdVars(const remid::ReMId &remid, SRKalmanTrack &srTrk)
Definition: FillPIDs.cxx:36
int fNmicroSliceDataNotPresent
of microslices with !DataPresent
EventNumber_t event() const
Definition: EventID.h:116
unsigned short month
month of spill
Definition: SRHeader.h:38
std::string generatorVersion
event generator version
virtual void setDefault()
Definition: SRVeto.cxx:35
unsigned int ncellsfromedge
minimum number of cells to edge of detector
Definition: SRSlice.h:30
CAFMaker(const Parameters &params)
bool isvtxcont
Checks if neutrino true vertex is within detector.
Definition: SRNeutrino.h:36
Atom< string > CVNParticle2DLabel
double fvy
Definition: MCFlux.h:53
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
double fmupare
Definition: MCFlux.h:70
unsigned int nnonnoise
Definition: SRSlice.h:31
void FillLIDEnergyVars(const slid::ShowerLID &slid, caf::SRFuzzyKProng &png)
Definition: FillPIDs.cxx:272
int nmicroslicedatanotpresent
# of microslices with !DataPresent
Definition: SRSpill.h:87
std::vector< TrackIDE > HitsToTrackIDE(const std::vector< const rb::CellHit * > &hits, bool useBirksE=false) const
Returns vector of TrackIDE structs contributing to the given collection of hits.
unsigned short dibmask
Definition: SRSpill.h:58
void fillSizes()
Definition: SRMichelE.cxx:28
Atom< string > NCPi0BkgRejLabel
Cosmic Rejection PIDs for Numu analysis.
Definition: FillParentInfo.h:9
Atom< string > RvpLabel
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
double EffPurMetric(const cheat::NeutrinoEffPur &ep)
Function for NeutrinoEffPur&#39;s nu interaction to slice efficiency * purity.
The StandardRecord is the primary top-level object in the Common Analysis File trees.
virtual void setDefault()
int ndcms
# of DCMs in partition; may not = # of LIVE DCMs = (hdr.diblast-hdr.dibfirst+1)*12 ...
Definition: SRSpill.h:76
signed long long int deltaspilltimensec
Delta time [ns].
Definition: SRSpill.h:37
static bool sortMuIdLength(const SRMuId &a, const SRMuId &b)
float len
track length [cm]
Definition: SRTrack.h:41
int nanoslicebufferfull
# of nanoslices reporting BufferFull
Definition: SRSpill.h:92
float ndhadcaltranE
Near detector – hadronic calorimetric energy NOT on the muon track in transition plane [GeV]...
Definition: SRNumuEnergy.h:43
Atom< string > CosmicTrackLabel
Atom< string > NumuLSTMEnergyLabel
unsigned int nnoisyapds
# of noisy APDs
Definition: SRSpill.h:67
Result, basic output of CVN neural net.
Definition: Result.h:15
void FillTrackVarsWithOverlapE(const rb::Track &trk, const rb::Energy &energy, SRTrack &srTrk, const std::vector< rb::Cluster > &sliceList, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, int sliceIdx)
Function to fill SRTrack branch if overlapping energy object is present.
Definition: FillReco.cxx:298
const Var kDistAllFront([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngfront)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngfront);})
Distance of all showers in slice from the front edge of detector.
Definition: NueVars.h:45
float GeV() const
Definition: RecoHit.cxx:69
int kalfwdcell
Kalmantrack projected # cells from end point forwards to det edge.
Definition: SRContain.h:40
int Target() const
Definition: MCNeutrino.h:151
int nSliceHits
Number of hits from this neutrino in this slice.
Definition: BackTracker.h:53
SRNumuEnergy numu
Numu energy estimator.
double DetHalfWidth() const
bool GetAssociatedProduct(const art::FindManyP< T > &fm, int idx, T &ret) const
Retrieve an object from an association, with error handling.
const art::Ptr< simb::MCTruth > & TrackIDToMCTruth(int const &id) const
double DistToEdgeXMinus(TVector3 vertex)
float visEBirks
Sum of FLS hits that made CellHits from this neutrino [GeV] with birks suppression.
Definition: SRNeutrino.h:31
int nanoslicecommerror
# of nanoslices reporting CommError
Definition: SRSpill.h:93
double fppenergy
Definition: MCFlux.h:61
Atom< string > WindowTrackLabel
Atom< string > TrackInfoBPF