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