CAFMaker_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   This module creates Common Analysis Files.
/// \author  rocco@physics.umn.edu
/// \date    July 2012
////////////////////////////////////////////////////////////////////////

#include "CAFMaker/CAFMakerParams.h"
#include "CAFMaker/FillEnergies.h"
#include "CAFMaker/FillPIDs.h"
#include "CAFMaker/FillParentInfo.h"
#include "CAFMaker/FillReco.h"
#include "CAFMaker/FillTruth.h"
#include "CAFMaker/Utils.h"

// C/C++ includes
#include <fenv.h>
#include <time.h>
#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <map>
#include <sstream>
#include <string>
#include <vector>

#ifdef DARWINBUILD
#include <libgen.h>
#endif

// ROOT includes
#include "TFile.h"
#include "TH1D.h"
#include "TTree.h"

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/FileBlock.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "canvas/Persistency/Common/FindMany.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/FindOne.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "cetlib_except/exception.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include <IFDH_service.h>

// NOvASoft includes
// Geo/Cal/Util
#include "Calibrator/Calibrator.h"
#include "ChannelInfo/BadChanList.h"
#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/PlaneGeo.h"
#include "Metadata/MetadataManager.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "RawData/RawTrigger.h"
#include "RunHistory/service/RunHistoryService.h"
#include "Utilities/func/MathUtil.h"

// DataQuality
#include "DataQuality/DAQHeader/DAQEventSummary.h"

// Sim/BackTracker
#include "MCCheater/BackTracker.h"
#include "MCReweight/FluxWeights.h"
#include "Simulation/GENIEReweightTable.h"
#include "nusimdata/SimulationBase/GTruth.h"
#include "nusimdata/SimulationBase/MCFlux.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"
#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"

// RecoBase
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Energy.h"
#include "RecoBase/FilterList.h"
#include "RecoBase/FitSum.h"
#include "RecoBase/Prong.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Shower.h"
#include "RecoBase/Track.h"
#include "RecoBase/Vertex.h"
#include "RecoBase/WeightedProng.h"

// Other Reco
#include "CVN/func/CVNFeatures.h"
#include "CVN/func/PixelMap.h"
#include "CVN/func/ProngTrainingData.h"
#include "CVN/func/Result.h"
#include "CVN/func/RegResult.h"
#include "CVN/func/TrainingData.h"
#include "CosRej/CosRejObj.h"
#include "CosRej/NueCosRej.h"
#include "CosRej/TrkCntObj.h"
#include "LEM/PIDDetails.h"
#include "LSTME/prediction/LSTMEnergy.h"
#include "MEFinder/MEClusters.h"
#include "MuonRemove/art/MRCCParent.h"
#include "NCID/NCCosRej.h"
#include "NumuEnergy/NumuE.h"
#include "Preselection/PreselObj.h"
#include "Preselection/Veto.h"
#include "ReMId/ReMId.h"
#include "RecVarPID/RVP.h"
#include "RecoJMShower/EID.h"
#include "RecoJMShower/JMShower.h"
#include "ShowerLID/DedxDistribution.h"
#include "ShowerLID/EventLID.h"
#include "ShowerLID/ShowerLID.h"
#include "ShowerLID/ShowerPID.h"
#include "ShowerLID/SliceLID.h"
#include "SliceLID/prediction/Prediction.h"
#include "XSecReco/art/NCPi0BkgRej.h"
#include "XnuePID/Xnue.h"
#include "TensorFlowProducts/CVNCosmicFiltList.h"
#include "TrackInfo/TrackInfoObj.h"

// SummData
#include "SummaryData/CosmicExposure.h"
#include "SummaryData/EventQuality.h"
#include "SummaryData/POTSum.h"
#include "SummaryData/SpillData.h"

// OscLib - for woscdumb
#include "OscLib/OscCalcDumb.h"

// StandardRecord
#include "StandardRecord/SRSpillTruthBranch.h"
#include "StandardRecord/StandardRecord.h"
// CAFMaker support
#include "CAFMaker/Blinding.h"

#include "Utilities/AssociationUtil.h"

// Note it is convention to NOT declare ANY "using namespace xxx" here
// especially "using namespace std"
namespace caf {
/// Module to create Common Analysis Files from ART files
class CAFMaker : public art::EDProducer {
 public:
  // Allows 'nova --print-description' to work
  using Parameters = art::EDProducer::Table<CAFMakerParams>;

  explicit CAFMaker(const Parameters& params);
  virtual ~CAFMaker();

  void produce(art::Event& evt) noexcept;

  void respondToOpenInputFile(const art::FileBlock& fb);

  void beginJob();
  void endJob();
  virtual void beginRun(art::Run& r);
  virtual void beginSubRun(art::SubRun& sr);
  virtual void endSubRun(art::SubRun& sr);

 protected:
  CAFMakerParams fParams;

  std::string fCafFilename;

  bool fIsRealData;  // use this instead of evt.isRealData(), see init in
                     // produce(evt)
  double fTotalPOT;
  double fTotalSinglePOT;
  double fTotalEvents;
  double fTotalLivetime;
  int fTotalTrueNonswapNue;
  int fTotalTrueSingleNue;
  int fCycle;
  int fBatch;
  int fMask;

  TFile* fFile;
  TTree* fRecTree;
  TH1D* hPOT;
  TH1D* hSinglePOT;
  TH1D* hTotalTrueNonswapNue;
  TH1D* hTotalTrueSingleNue;
  TH1D* hEvents;
  TH1D* hLivetime;
  art::ServiceHandle<cheat::BackTracker> bt;
  art::ServiceHandle<geo::Geometry> geom;
  art::ServiceHandle<geo::LiveGeometry> livegeom;
  art::ServiceHandle<nova::dbi::RunHistoryService> rh;
  art::ServiceHandle<chaninfo::BadChanList> bc;

  novadaq::cnv::DetId
      fDetID;  ///< Detector ID in nova daq convention namespace typedef
  Det_t fDet;  ///< Detector ID in caf namespace typedef

  // Creating separate friend tree and branches for spill info.
  // This tree will store one entry per event, even empty ones
  // so that POT counting can be done right.
  TTree* fSpillTree;

  //Seperate tree for CVN cosmic rejeciton.  This operates on slices of time
  //and not physics slices.  So it must be stored seperately.
  TTree* fCosmicCVNTree;

  TTree* fNuTree;

  void InitializeOutfile();

  void AddMetadataToFile(std::map<std::string, std::string> metadata);

  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);

  void drawSliceTruthTable(
      const std::vector<std::vector<cheat::NeutrinoEffPur>>& sliceTruthTable,
      const std::vector<int>& faveIds);
  void drawSliceTruthTable(
      const std::vector<std::vector<cheat::NeutrinoEffPur>>& sliceTruthTable,
      const std::vector<int>& faveIds, const std::vector<int>& otherFaveIds);

  /// Equivalent of FindManyP except a return that is !isValid() prints a
  /// messsage and aborts if StrictMode is true.
  template <class T, class U>
  art::FindManyP<T> FindManyPStrict(const U& from, const art::Event& evt,
                                    const art::InputTag& label) const;

  /// \brief Retrieve an object from an association, with error handling
  ///
  /// This can go wrong in two ways: either the FindManyP itself is
  /// invalid, or the result for the requested index is empty. In most
  /// cases these have the same response, so conflating them here
  /// saves redundancy elsewhere.
  ///
  /// \param      fm  The FindManyP object describing the association
  /// \param      idx Which element of the FindManyP to look it
  /// \param[out] ret The product retrieved
  /// \return          Whether \a ret was filled
  template <class T>
  bool GetAssociatedProduct(const art::FindManyP<T>& fm, int idx, T& ret) const;

  /// Equivalent of evt.getByLabel(label, handle) except failedToGet
  /// prints a message and aborts if StrictMode is true.
  template <class T>
  void GetByLabelStrict(const art::Event& evt, const std::string& label,
                        art::Handle<T>& handle) const;

  /// Equivalent of evt.getByLabel(label, handle) except failedToGet
  /// prints a message.
  template <class T>
  void GetByLabelIfExists(const art::Event& evt, const std::string& label,
                          art::Handle<T>& handle) const;

  /// \param      pset The parameter set
  /// \param      name Pass "foo.bar.baz" as {"foo", "bar", "baz"}
  /// \param[out] ret  Value of the key, not set if we return false
  /// \return          Whether the key was found
  template <class T>
  bool GetPsetParameter(const fhicl::ParameterSet& pset,
                        const std::vector<std::string>& name, T& ret) const;

  static bool EssentiallyEqual(double a, double b, double precision = 0.0001) {
    return a <= (b + precision) && a >= (b - precision);
  }

  static bool sortTrackLength(const SRTrack& a, const SRTrack& b) {
    return a.len > b.len;
  }

  static bool sortRBTrackLength(const art::Ptr<rb::Track>& a,
                                const art::Ptr<rb::Track>& b) {
    return a->TotalLength() > b->TotalLength();
  }

  static bool sortRBProngLength(const art::Ptr<rb::Prong>& a,
                                const art::Ptr<rb::Prong>& b) {
    return a->TotalLength() > b->TotalLength();
  }

  static bool sortShowerEnergy(const SRFuzzyKProng& a, const SRFuzzyKProng& b) {
    return a.shwlid.shwE > b.shwlid.shwE;
  }

  static bool sortProngCalE(const SRFuzzyKProng& a, const SRFuzzyKProng& b) {
    return a.calE > b.calE;
  }

  static bool sortElasticProng(const SRElastic& a, const SRElastic& b) {
    return a.fuzzyk.png.size() > b.fuzzyk.png.size();
  }

  static bool sortHoughProng(const SRHoughVertex& a, const SRHoughVertex& b) {
    return a.fuzzyk.png.size() > b.fuzzyk.png.size();
  }

  static bool sortMuIdLength(const SRMuId& a, const SRMuId& b) {
    return a.len > b.len;
  }

  static bool sortRemid(const SRKalmanTrack& a, const SRKalmanTrack& b) {
    return a.rempid > b.rempid;
  }

  static bool sortNuWithIdx(const cheat::NeutrinoWithIndex& a,
                            const cheat::NeutrinoWithIndex& b) {
    return a.truthColIndex < b.truthColIndex;
  }

  std::pair<int, int> calcFirstLastLivePlane(int plane, std::bitset<14> binary,
                                             caf::Det_t det = caf::kFARDET);

  double SimpleOscProb(const simb::MCFlux& flux,
                       const simb::MCNeutrino& nu) const;

  void FillSpillVars(art::Handle<sumdata::SpillData> spillPot,
                     art::Handle<sumdata::EventQuality> spillQual,
                     art::Handle<dq::DAQEventSummary> daq,
                     art::Handle<std::vector<rawdata::RawTrigger>> trigs,
                     SRSpill& spill, art::Event& evt) const;

  enum class MixingType {
    eRockOverlay,
    eCosmicOverlay,
    eRockCosmicOverlay,
    eRockSinglesOverlay,
    eNone,
    ePileup,
    eSingles,
    eOverlay,
    eMCSpillSinglesOverlay,
    eDataSpillSinglesOverlay
  };

  std::map<std::string, MixingType> mapStringMixingType = {
      {"rock_overlay", MixingType::eRockOverlay},
      {"cosmic_overlay", MixingType::eCosmicOverlay},
      {"rock_cosmic_overlay", MixingType::eRockCosmicOverlay},
      {"rock_singles_overlay", MixingType::eRockSinglesOverlay},
      {"mcspill_singles_overlay", MixingType::eMCSpillSinglesOverlay},
      {"dataspill_singles_overlay", MixingType::eDataSpillSinglesOverlay},
      {"none", MixingType::eNone},
      {"pileup", MixingType::ePileup},
      {"singles", MixingType::eSingles},
      {"overlay", MixingType::eOverlay}};

  std::map<MixingType, std::string> mapMixingTypeString = {
      {MixingType::eRockOverlay, "rock_overlay"},
      {MixingType::eCosmicOverlay, "cosmic_overlay"},
      {MixingType::eRockCosmicOverlay, "rock_cosmic_overlay"},
      {MixingType::eRockSinglesOverlay, "rock_singles_overlay"},
      {MixingType::eMCSpillSinglesOverlay, "mcspill_singles_overlay"},
      {MixingType::eDataSpillSinglesOverlay, "dataspill_singles_overlay"},
      {MixingType::eNone, "none"},
      {MixingType::ePileup, "pileup"},
      {MixingType::eSingles, "singles"},
      {MixingType::eOverlay, "overlay"}};
};

//.......................................................................

//.......................................................................
CAFMaker::CAFMaker(const Parameters& params)
    : fParams(params()), fIsRealData(false), fFile(0) {
  // We update this, so have to pull it out of the config
  fCafFilename = fParams.CAFFilename();

  // Normally CAFMaker is run without an output ART stream, so these go
  // nowhere, but can be occasionally useful for filtering as part of an
  // ART job.
  produces<std::vector<caf::StandardRecord>>();
  produces<art::Assns<caf::StandardRecord, rb::Cluster>>();
  produces<caf::SRSpill>();
  if (fParams.FillNuTree()) produces<std::vector<caf::SRSpillTruthBranch>>();
}

//......................................................................
CAFMaker::~CAFMaker() {}

//......................................................................
void CAFMaker::respondToOpenInputFile(const art::FileBlock& fb) {
  if (!fFile) {
    // Filename wasn't set in the FCL, and this is the
    // first file we've seen
    char temp[fb.fileName().size() + 1];
    std::strcpy(temp, fb.fileName().c_str());
    fCafFilename = basename(temp);
    const size_t dotpos = fCafFilename.find('.');
    assert(dotpos != std::string::npos);  // Must have a dot, surely?
    fCafFilename.resize(dotpos);
    fCafFilename += fParams.FileExtension();

    InitializeOutfile();
  }
}

//......................................................................
void CAFMaker::beginJob() {
  if (!fCafFilename.empty()) InitializeOutfile();
}

//......................................................................
void CAFMaker::beginRun(art::Run& r) {
  fDetID = geom->DetId();
  fDet = (Det_t)fDetID;

  if (!rh->GetDiBlockMaskFromCondb() && !fParams.ApplyingFilter()) {  // get for each run; here if subrun
                                                                      // masking not used and there is only 
                                                                      // one stream with calhit run
    art::Handle<int> mask;
    r.getByLabel(fParams.HitsLabel(), mask);  // why doesn't fParams.HitInput() work here?
    if (!mask.isValid()) {
      mf::LogError("CAFMaker") << "Could not find a diBlock mask from " << fParams.HitsLabel() << ". Bailing.";
      abort();
    }
      
    fMask = *mask;
  }


  if (!rh->GetDiBlockMaskFromCondb() && fParams.ApplyingFilter()) {  // get for each run; here if subrun
                                                                     // masking not used and there are two 
                                                                     // streams one of which does not have 
                                                                     // calhits to save space and time
    std::unique_ptr<int> mask(new int(rh->GoodDiBlockMask(0,true)));
          
    fMask = *mask;
  }
}

//......................................................................
void CAFMaker::beginSubRun(art::SubRun& sr) {
  if (rh->GetDiBlockMaskFromCondb() && !fParams.ApplyingFilter()) {  // get if we are using subrun-level
                                                                     // masks and there is only 
                                                                     // one stream with calhit run
    art::Handle<int> mask;
    sr.getByLabel(fParams.HitsLabel(), mask);  // why doesn't fParams.HitInput() work here?
    if (!mask.isValid()) {
      mf::LogError("CAFMaker") << "Could not find a diBlock mask from " << fParams.HitsLabel() << ". Bailing.";
      abort();
    }

    fMask = *mask;
  }

  if (rh->GetDiBlockMaskFromCondb() && fParams.ApplyingFilter()) {  // get if we are using subrun-level
                                                                    // masks and there are two 
                                                                    // streams one of which does not have 
                                                                    // calhits to save space and time
    std::unique_ptr<int> mask(new int(rh->GoodDiBlockMask(sr.subRun(),false)));
    
    fMask = *mask;
  }
}

//......................................................................
void CAFMaker::InitializeOutfile() {
  assert(!fFile);
  assert(!fCafFilename.empty());

  mf::LogInfo("CAFMaker") << "Output filename is " << fCafFilename;

  fFile = new TFile(fCafFilename.c_str(), "RECREATE");

  hPOT = new TH1D("TotalPOT", "TotalPOT;; POT", 1, 0, 1);
  hSinglePOT =
      new TH1D("TotalSinglePOT", "TotalSinglePOT;; Single POT", 1, 0, 1);
  hTotalTrueNonswapNue =
      new TH1D("TotalTrueNonswapNue",
               "TotalTrueNonswap;; Total True Nonswap Nue", 1, 0, 1);
  hTotalTrueSingleNue = new TH1D(
      "TotalTrueSingleNue", "TotalTrueSingle;; Total True Single Nue", 1, 0, 1);
  hEvents = new TH1D("TotalEvents", "TotalEvents;; Events", 1, 0, 1);
  hLivetime = new TH1D("TotalLivetime", "TotalLiveTime;;Seconds", 1, 0, 1);

  fRecTree = new TTree("recTree", "records");

  // Tell the tree it's expecting StandardRecord objects
  StandardRecord* rec = 0;
  fRecTree->Branch("rec", "caf::StandardRecord", &rec);

  // Create tree to store spill info for POT counting
  fSpillTree = new TTree("spillTree", "spills");

  caf::SRSpill* dummy = 0;
  fSpillTree->Branch("spill", "caf::SRSpill", &dummy);

  // Make the spillTree a friend of the fRecTree
  fRecTree->AddFriend("spillTree");


  if (fParams.FillNuTree()) {
    fNuTree = new TTree("nuTree", "spillTruth");
    caf::SRSpillTruthBranch* dummy2 = 0;
    fNuTree->Branch("mc", "caf::SRSpillTruthBranch", &dummy2);
  } else {
    fNuTree = 0;
  }

  fTotalPOT = 0;
  fTotalSinglePOT = 0;
  fTotalTrueNonswapNue = 0;
  fTotalTrueSingleNue = 0;
  fTotalEvents = 0;
  fTotalLivetime = 0;
  fCycle = -5;
  fBatch = -5;

  // Global information about the processing details:

  std::map<TString, TString> envmap;
  std::string envstr;
  // Environ comes from unistd.h
  // environ is not present on OSX for some reason, so just use getenv to
  // grab the variables we care about.
#ifdef DARWINBUILD
  std::set<TString> variables;
  variables.insert("USER");
  variables.insert("HOSTNAME");
  variables.insert("PWD");
  variables.insert("SRT_PUBLIC_CONTEXT");
  variables.insert("SRT_PRIVATE_CONTEXT");
  for (auto var : variables) envmap[var] = getenv(var);
#else
  for (char** penv = environ; *penv; ++penv) {
    const std::string pair = *penv;
    envstr += pair;
    envstr += "\n";
    const size_t split = pair.find("=");
    if (split == std::string::npos) continue;  // Huh?
    const std::string key = pair.substr(0, split);
    const std::string value = pair.substr(split + 1);
    envmap[key] = value;
  }
#endif

  // Get the command-line we were invoked with. What I'd really like is
  // just the fcl script and list of input filenames in a less hacky
  // fashion. I'm not sure that's actually possible in ART.
  // TODO: ask the artists.
  FILE* cmdline = fopen("/proc/self/cmdline", "rb");
  char* arg = 0;
  size_t size = 0;
  std::string cmd;
  while (getdelim(&arg, &size, 0, cmdline) != -1) {
    cmd += arg;
    cmd += " ";
  }
  free(arg);
  fclose(cmdline);

  fFile->mkdir("env")->cd();

  TObjString(envmap["USER"]).Write("user");
  TObjString(envmap["HOSTNAME"]).Write("hostname");
  TObjString(envmap["PWD"]).Write("pwd");
  TObjString(envmap["SRT_PUBLIC_CONTEXT"]).Write("publiccontext");
  TObjString(envmap["SRT_PRIVATE_CONTEXT"]).Write("privatecontext");
  // Default constructor is "now"
  TObjString(TTimeStamp().AsString()).Write("date");
  TObjString(cmd.c_str()).Write("cmd");
  TObjString(fCafFilename.c_str()).Write("output");
  TObjString(envstr.c_str()).Write("env");
}

void CAFMaker::AddMetadataToFile(std::map<std::string, std::string> metadata) {
  assert(fFile && "CAFMaker: Trying to add metadata to an uninitialized file");

  fFile->mkdir("metadata")->cd();

  for (auto const& pair : metadata) {
    TObjString(pair.second.c_str()).Write(pair.first.c_str());
  }
}

//......................................................................
void CAFMaker::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) {
  // Get information about truNu
  vec->push_back(SRNeutrino());
  SRNeutrino& srTruth = vec->back();

  srTruth.run = evt.run();
  srTruth.subrun = evt.subRun();
  srTruth.evt = evt.id().event();
  srTruth.cycle = fCycle;

  srTruth.isFHC = spill.isFHC;
  srTruth.is0HC = spill.is0HC;
  srTruth.isRHC = spill.isRHC;

  srTruth.det = fDet;

  const art::Ptr<simb::MCTruth>& tru = effPur->neutrinoInt;

  srTruth.generator = CAFGeneratorEnum(tru->GeneratorInfo().generator);

  // this is very rudimentary caching
  static std::string oldVersionStr = std::string{};
  static std::vector<unsigned int> oldVersionVec = std::vector<unsigned int>{};
  static std::unordered_map<std::string, std::string> oldGenConfig = std::unordered_map<std::string, std::string>{};
  static std::string oldGenConfigStr = std::string{};
  if (tru->GeneratorInfo().generatorVersion != oldVersionStr || tru->GeneratorInfo().generatorConfig != oldGenConfig)
  {
    oldVersionStr = tru->GeneratorInfo().generatorVersion;
    oldVersionVec = DecodeGeneratorVersion(tru->GeneratorInfo().generatorVersion);
    oldGenConfig = tru->GeneratorInfo().generatorConfig;
    oldGenConfigStr = "";
    for (const auto & configPair: oldGenConfig)
    {
      if (configPair.first == "tune")
        oldGenConfigStr = configPair.second;
      else
        mf::LogWarning("CAFMaker") << "Don't know how to store generator config parameter named '" << configPair.first << "'";
    }
  }
  srTruth.genVersion = oldVersionVec;
  if (!oldGenConfigStr.empty())
    srTruth.genConfigString = oldGenConfigStr;


  srTruth.pdg = tru->GetNeutrino().Nu().PdgCode();
  srTruth.hitnuc = tru->GetNeutrino().HitNuc();

  srTruth.mode = tru->GetNeutrino().Mode();  // to avoid depending on simb
  srTruth.iscc = (tru->GetNeutrino().CCNC() == simb::kCC);

  if (srTruth.mode == kElectronScattering) {
    // For a numu or nutau scattering off an electron, the exchanged particle
    // is a Z. But for a nue, there are two diagrams, one with a Z and one
    // with a W. Both the final states have the electron ejected, so the
    // interaction is a quantum-mechanical superposition between the
    // two. It's neither really CC or NC. We do need to set this boolean
    // though. Treating nu-on-e scattering as NC (ie not oscillating it) is
    // wrong, the cross-section for nue is different to numu and nutau so
    // oscillations do have an effect on the final spectrum. If we mark it CC
    // the different oscillated fluxes will be weighted together in the fit
    // appropriately, so this is safe, even if not perfectly physically
    // accurate.
    srTruth.iscc = true;
  }
  srTruth.inttype = tru->GetNeutrino().InteractionType();

  srTruth.q2 = tru->GetNeutrino().QSqr();
  srTruth.x = tru->GetNeutrino().X();
  srTruth.y = tru->GetNeutrino().Y();
  srTruth.W2 = tru->GetNeutrino().W() * tru->GetNeutrino().W();
  srTruth.E = tru->GetNeutrino().Nu().Momentum().E();
  srTruth.visE = effPur->energyTotal;
  srTruth.visEinslc = effPur->energySlice;
  srTruth.eff = effPur->efficiency;
  srTruth.pur = effPur->purity;
  srTruth.nhitslc = effPur->nSliceHits;
  srTruth.nhittot = effPur->nTotalHits;
  srTruth.time = tru->GetNeutrino().Nu().Position().T();
  srTruth.vtx = tru->GetNeutrino().Nu().Position().Vect();
  srTruth.p = tru->GetNeutrino().Nu().Momentum();

  float visEBirks = 0;
  for (std::vector<cheat::TrackIDE>::iterator it = allTracksBirks.begin();
       it != allTracksBirks.end(); ++it) {
    if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
    visEBirks += it->energy;
  }
  float visEinslcBirks = 0;
  for (std::vector<cheat::TrackIDE>::iterator it = sliceTracksBirks.begin();
       it != sliceTracksBirks.end(); ++it) {
    if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
    visEinslcBirks += it->energy;
  }

  srTruth.visEBirks = visEBirks;
  srTruth.visEinslcBirks = visEinslcBirks;

  // Get tgtZ / tgtA from the MCTruth - which stores pdg code of nucleus
  int tarpdg = tru->GetNeutrino().Target();
  // pdgcode of nucleus = 10LZZZAAAI
  srTruth.tgtZ = (tarpdg / 10000) % 1000;
  srTruth.tgtA = (tarpdg / 10) % 1000;

  srTruth.genweight = tru->GetNeutrino().Nu().Weight();

  TVector3 nuVtx = tru->GetNeutrino().Nu().Position().Vect();
  bool inDet = fabs(nuVtx.X()) <= geom->DetHalfWidth() &&
               fabs(nuVtx.Y()) <= geom->DetHalfHeight() && nuVtx.Z() >= 0 &&
               nuVtx.Z() <= geom->DetLength();
  srTruth.isvtxcont = inDet;

  // rescue the parts of GTruth that we always want, put them in SRNeutrino
  art::PtrVector<simb::MCTruth> pv;
  pv.push_back(tru);
  art::FindManyP<simb::GTruth> fmGtruth =
      FindManyPStrict<simb::GTruth>(pv, evt, fParams.GeneratorLabel());

  if (fmGtruth.isValid()) {
    std::vector<art::Ptr<simb::GTruth>> gtruths = fmGtruth.at(0);

    if (!gtruths.empty()) {
      const simb::GTruth& gtruth = *gtruths[0];

      // Useful for identifying pre-FSI generated hadron content
      srTruth.npiplus = gtruth.fNumPiPlus;
      srTruth.npiminus = gtruth.fNumPiMinus;
      srTruth.npizero = gtruth.fNumPi0;
      srTruth.nproton = gtruth.fNumProton;
      srTruth.nneutron = gtruth.fNumNeutron;

      // Descriptive variables relating to the thrown cross sections
      srTruth.ischarm = gtruth.fIsCharm;
      srTruth.isseaquark = gtruth.fIsSeaQuark;
      srTruth.resnum = gtruth.fResNum;
      srTruth.xsec = gtruth.fXsec;
    }
  }

  // now stuff from the flux record
  //    art::FindManyP<simb::MCFlux> fmFlux(pv, evt, fParams.GeneratorLabel());
  art::FindManyP<simb::MCFlux> fmFlux =
      FindManyPStrict<simb::MCFlux>(pv, evt, fParams.GeneratorLabel());
  if (fmFlux.isValid()) {
    std::vector<art::Ptr<simb::MCFlux>> fluxes = fmFlux.at(0);

    if (!fluxes.empty()) {
      const simb::MCFlux& flux = *fluxes[0];

      srTruth.pdgorig = flux.fntype;
      srTruth.woscdumb = SimpleOscProb(flux, tru->GetNeutrino());

      SRBeam& beam = srTruth.beam;

      beam.tv = TVector3(flux.ftvx, flux.ftvy, flux.ftvz);
      beam.tp = TVector3(flux.ftpx, flux.ftpy, flux.ftpz);
      beam.runjob = flux.frun;
      beam.potnum = flux.fevtno;
      beam.tptype = flux.ftptype;
      beam.nimpwt = flux.fnimpwt;
      beam.ndecay = flux.fndecay;
      beam.v = TVector3(flux.fvx, flux.fvy, flux.fvz);
      beam.pdp = TVector3(flux.fpdpx, flux.fpdpy, flux.fpdpz);
      beam.ppdxdz = flux.fppdxdz;
      beam.ppdydz = flux.fppdydz;
      beam.pppz = flux.fpppz;
      beam.ppenergy = flux.fppenergy;
      beam.ppmedium = flux.fppmedium;
      beam.ptype = flux.fptype;
      beam.ppv = TVector3(flux.fppvx, flux.fppvy, flux.fppvz);
      beam.muparp = TVector3(flux.fmuparpx, flux.fmuparpy, flux.fmuparpz);
      beam.mupare = flux.fmupare;
      beam.necm = flux.fnecm;
      beam.tgen = flux.ftgen;
      beam.tgptype = flux.ftgptype;
      beam.dk2gen = flux.fdk2gen;
      beam.gen2vtx = flux.fgen2vtx;
      beam.dk2vtx = flux.fdk2gen + flux.fgen2vtx;

    } else {
      // this does not work for data-singlenu-mc overlay, but we don't want the
      // singlenu flux anyway
      if (!fParams.DataSpillHasMC())
        throw cet::exception(
            "Missing Association Error",
            "Failed to find MCFlux object assoicated with MCTruth.  This "
            "probably means something went wrong with overlays.");
    }
  }

  // Must be after the setting of 'vtx' and 'beam' info
  srTruth.L = TrueNeutrinoDistance(fDetID, srTruth);

  // Get the neutrino's daughters (and other progeny)
  std::vector<const sim::Particle*> particles = bt->MCTruthToParticles(tru);

  // Get the sim::TrueEnergy's out of the event.
  art::Handle<std::vector<sim::TrueEnergy>> tes;
  GetByLabelStrict(evt, fParams.TrueEnergyLabel(), tes);
  std::vector<sim::TrueEnergy> TrueEnergies;
  if (!tes.failedToGet()) TrueEnergies = *tes;

  // We're going to loop over daughters to add other primaries to truNu.
  // First we're gonna get the lepton and add it as the zeroth entry,
  // so let's get its information.
  int lepPdg = tru->GetNeutrino().Lepton().PdgCode();
  double lepEnergy = tru->GetNeutrino().Lepton().E();
  int lepCount = 0;    // count the found leptons, for debugging.
  int lepTrackId = 0;  // We'll set this later
  // Now we'll loop over the lepton and add it
  float visENeutron = 0;
  for (std::vector<cheat::TrackIDE>::iterator it = allTracks.begin();
       it != allTracks.end(); ++it) {
    if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
    int mId = bt->TrackIDToParticle(it->trackID)->Mother();
    while (mId != 0) {
      if (bt->TrackIDToParticle(mId)->PdgCode() == 2112) {
        visENeutron += it->energy;
        break;
      }
      mId = bt->TrackIDToParticle(mId)->Mother();
    }
  }
  float visEinslcNeutron = 0;
  for (std::vector<cheat::TrackIDE>::iterator it = sliceTracks.begin();
       it != sliceTracks.end(); ++it) {
    if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
    int mId = bt->TrackIDToParticle(it->trackID)->Mother();
    while (mId != 0) {
      if (bt->TrackIDToParticle(mId)->PdgCode() == 2112) {
        visEinslcNeutron += it->energy;
        break;
      }
      mId = bt->TrackIDToParticle(mId)->Mother();
    }
  }
  float visENeutronBirks = 0;
  for (std::vector<cheat::TrackIDE>::iterator it = allTracksBirks.begin();
       it != allTracksBirks.end(); ++it) {
    if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
    int mId = bt->TrackIDToParticle(it->trackID)->Mother();
    while (mId != 0) {
      if (bt->TrackIDToParticle(mId)->PdgCode() == 2112) {
        visENeutronBirks += it->energy;
        break;
      }
      mId = bt->TrackIDToParticle(mId)->Mother();
    }
  }
  float visEinslcNeutronBirks = 0;
  for (std::vector<cheat::TrackIDE>::iterator it = sliceTracksBirks.begin();
       it != sliceTracksBirks.end(); ++it) {
    if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
    int mId = bt->TrackIDToParticle(it->trackID)->Mother();
    while (mId != 0) {
      if (bt->TrackIDToParticle(mId)->PdgCode() == 2112) {
        visEinslcNeutronBirks += it->energy;
        break;
      }
      mId = bt->TrackIDToParticle(mId)->Mother();
    }
  }
  srTruth.visENeutron = visENeutron;
  srTruth.visEinslcNeutron = visEinslcNeutron;
  srTruth.visENeutronBirks = visENeutronBirks;
  srTruth.visEinslcNeutronBirks = visEinslcNeutronBirks;

  for (const sim::Particle* part : particles) {
    // Add particle if it's a primary and it has the correct PDG code
    // and the correct energy.
    if (part->Mother() == 0 && part->PdgCode() == lepPdg &&
        EssentiallyEqual(part->E(), lepEnergy, .01)) {
      lepTrackId = part->TrackId();
      ++lepCount;
      AddParticleToVec(*part, allTracks, sliceTracks, allTracksBirks,
                       sliceTracksBirks, &srTruth.prim, TrueEnergies);
    }
  }

  if (lepCount != 1) {
    switch (lepCount) {
      case 0:
        if (!srTruth.iscc) break;

        mf::LogError("CAFMaker")
            << "*****************************************"
            << "\n"
            << "Uh, oh, CAFMaker did not find primary lepton. "
            << "\n"
            << "We're going to need to find a way to deal with this."
            << "\n"
            << "Info about the particle: "
            << "\n"
            << "  pdg: " << srTruth.pdg << "\n"
            << "  pdgorig: " << srTruth.pdgorig << "\n"
            << "  E: " << srTruth.E << "\n"
            << "  visE: " << srTruth.visE << "\n"
            << "  nhittot: " << srTruth.nhittot << "\n"
            << "  X: " << srTruth.vtx.X() << "\n"
            << "  Y: " << srTruth.vtx.Y() << "\n"
            << "  Z: " << srTruth.vtx.Z() << "\n"
            << "  iscc: " << srTruth.iscc << "\n"
            << "  infid(nd): "
            << (srTruth.vtx.X() < 200 && srTruth.vtx.X() > -200 &&
                srTruth.vtx.Y() < 200 && srTruth.vtx.Y() > -200 &&
                srTruth.vtx.Z() > 0 && srTruth.vtx.X() < 1600)
            << "\n"

            << "*****************************************"
            << "\n";
        break;
      default:
        mf::LogError("CAFMaker")
            << "*****************************************"
            << "\n"
            << "Uh, oh, CAFMaker found too many primary leptons in event. "
            << "\n"
            << "We're going to need to find a way to deal with this."
            << "\n"
            << "*****************************************";
        break;
    }
    // assert(lepCount == 1);
  }

  for (const sim::Particle* part : particles) {
    // Add particle if it's a primary and it's not the lepton
    if (part->TrackId() != lepTrackId && part->Mother() == 0)
      AddParticleToVec(*part, allTracks, sliceTracks, allTracksBirks,
                       sliceTracksBirks, &srTruth.prim, TrueEnergies);
  }

  FindAndAddMichels(particles, allTracks, &srTruth.michel);
}

//......................................................................
/// This function is here for debugging purposes.
/// It will draw a table of with slices as rows and neutrinos as columns.
/// The favorites will be written with a star next to them.
/// If two faveIds vectors are provided, the second will have
/// an @ drawn next to them
void CAFMaker::drawSliceTruthTable(
    const std::vector<std::vector<cheat::NeutrinoEffPur>>& sliceTruthTable,
    const std::vector<int>& faveIds) {
  std::vector<int> otherFaveIds(faveIds.size(), -1);
  drawSliceTruthTable(sliceTruthTable, faveIds, otherFaveIds);
}

//......................................................................

/// This function is here for debugging purposes.
/// It will draw a table of with slices as rows and neutrinos as columns.
/// The favorites will be written with a star next to them.
/// If two faveIds vectors are provided, the second will have
/// an @ drawn next to them
void CAFMaker::drawSliceTruthTable(
    const std::vector<std::vector<cheat::NeutrinoEffPur>>& sliceTruthTable,
    const std::vector<int>& faveIds, const std::vector<int>& otherFaveIds) {
  std::cout.precision(2);
  std::cout.width(2);
  std::cout.setf(std::ios::fixed, std::ios::floatfield);
  int iSlc = 0;
  for (auto slc : sliceTruthTable) {
    int iNu = 0;
    for (auto nu = slc.begin(); nu != slc.end(); ++nu, ++iNu) {
      std::cout << nu->efficiency;
      if (faveIds[iSlc] == iNu)
        std::cout << "*";
      else
        std::cout << " ";
      if (otherFaveIds[iSlc] == iNu)
        std::cout << "@";
      else
        std::cout << " ";
      std::cout << "         ";
    }
    std::cout << std::endl;
  }
}

//......................................................................
template <class T, class U>
art::FindManyP<T> CAFMaker::FindManyPStrict(const U& from,
                                            const art::Event& evt,
                                            const art::InputTag& tag) const {
  art::FindManyP<T> ret(from, evt, tag);

  if (!tag.label().empty() && !ret.isValid() && fParams.StrictMode()) {
    std::cout << "CAFMaker: No Assn from '"
              << abi::__cxa_demangle(typeid(from).name(), 0, 0, 0) << "' to '"
              << abi::__cxa_demangle(typeid(T).name(), 0, 0, 0)
              << "' found under label '" << tag << "'. "
              << "Set 'StrictMode: false' to continue anyway." << std::endl;
    abort();
  }

  return ret;
}

//......................................................................
template <class T>
bool CAFMaker::GetAssociatedProduct(const art::FindManyP<T>& fm, int idx,
                                    T& ret) const {
  if (!fm.isValid()) return false;

  const std::vector<art::Ptr<T>> prods = fm.at(idx);

  if (prods.empty()) return false;

  ret = *prods[0];

  return true;
}

//......................................................................
template <class T>
void CAFMaker::GetByLabelStrict(const art::Event& evt, const std::string& label,
                                art::Handle<T>& handle) const {
  evt.getByLabel(label, handle);
  if (!label.empty() && handle.failedToGet() && fParams.StrictMode()) {
    std::cout << "CAFMaker: No product of type '"
              << abi::__cxa_demangle(typeid(*handle).name(), 0, 0, 0)
              << "' found under label '" << label << "'. "
              << "Set 'StrictMode: false' to continue anyway." << std::endl;
    abort();
  }
}

//......................................................................
template <class T>
void CAFMaker::GetByLabelIfExists(const art::Event& evt,
                                  const std::string& label,
                                  art::Handle<T>& handle) const {
  evt.getByLabel(label, handle);
  if (!label.empty() && handle.failedToGet() && fParams.StrictMode()) {
    std::cout << "CAFMaker: No product of type '"
              << abi::__cxa_demangle(typeid(*handle).name(), 0, 0, 0)
              << "' found under label '" << label << "'. "
              << "Continuing without it." << std::endl;
  }
}

//......................................................................
template <class T>
bool CAFMaker::GetPsetParameter(const fhicl::ParameterSet& pset,
                                const std::vector<std::string>& name,
                                T& ret) const {
  fhicl::ParameterSet p = pset;
  for (unsigned int i = 0; i < name.size() - 1; ++i) {
    if (!p.has_key(name[i])) return false;
    p = p.get<fhicl::ParameterSet>(name[i]);
  }
  if (!p.has_key(name.back())) return false;
  ret = p.get<T>(name.back());
  return true;
}

//......................................................................
std::pair<int, int> CAFMaker::calcFirstLastLivePlane(int plane,
                                                     std::bitset<14> binary,
                                                     caf::Det_t det) {
  // if not the FD, return full ND size
  if (det != caf::kFARDET) return std::pair<int, int>(0, 213);
  int testDB = (plane / 64);
  int minblock = testDB;
  int maxblock = testDB;
  // find block boundaries;
  for (int i = testDB - 1; i >= 0; --i) {
    if (binary[i])
      minblock = i;
    else
      break;
  }
  for (int i = testDB + 1; i < 14; ++i) {
    if (binary[i])
      maxblock = i;
    else
      break;
  }
  return std::pair<int, int>(64 * (minblock), 64 * (maxblock + 1) - 1);
}

//......................................................................
void CAFMaker::produce(art::Event& evt) noexcept {
  // Normally CAFMaker is run without an output ART stream, so these go
  // nowhere, but can be occasionally useful for filtering as part of
  // an ART job.
  std::unique_ptr<std::vector<caf::StandardRecord>> srcol(
      new std::vector<caf::StandardRecord>);
  std::unique_ptr<art::Assns<caf::StandardRecord, rb::Cluster>> srAssn(
      new art::Assns<caf::StandardRecord, rb::Cluster>);
  std::unique_ptr<caf::SRSpill> spillcol(new SRSpill);
  std::unique_ptr<std::vector<caf::SRSpillTruthBranch>> spilltruthcol(
      new std::vector<SRSpillTruthBranch>);

  fTotalEvents += 1;
  // Adjust for art misfeature isRealData confusing data-mc overlay
  fIsRealData = evt.isRealData() || fParams.DataSpillHasMC();
  //   -- Case --        :  Data?  ==  isRealData || DataSpillHasMC :
  //   bt->HaveTruthInfo
  // Normal Data         :  true          true           false      :      false
  // Data Overlaid on MC :  true          false          true       :      true
  // Just MC             :  false         false          false      :      true

  // Get a handle on the slices

  art::Handle<std::vector<rb::Cluster>> slices;
  GetByLabelStrict(evt, fParams.ClusterLabel(), slices);

  /*    // Define a vector holding slices
  art::PtrVector<rb::Cluster> subevts;
  for(unsigned int i = 0; i < slices->size(); ++i)
    subevts.push_back(art::Ptr<rb::Cluster>(slices, i));
  art::Handle<std::vector<rb::Cluster> > hparentSlices;
  art::PtrVector<rb::Cluster> parentSubevts;*/

  if (fParams.IsSinglesOverlay()) {
    art::Handle<std::vector<simb::MCTruth>> mctruthcol;
    art::Handle<std::vector<simb::MCTruth>> singlemctruthcol;
    GetByLabelStrict(evt, fParams.GeneratorLabel(), mctruthcol);
    GetByLabelStrict(evt, fParams.SingleMixerLabel(), singlemctruthcol);

    if (mctruthcol->empty()) {
      mf::LogWarning("CAFMaker") << "Error retrieving MCTruth list";
      return;
    }
    if (singlemctruthcol->empty()) {
      mf::LogWarning("CAFMaker")
          << "Error retrieving MCTruth list from singles";
      return;
    }

    for (size_t i_intx = 0; i_intx < mctruthcol->size(); ++i_intx) {
      LOG_DEBUG("outnNueCCOverlay") << "start loop";
      simb::MCTruth const& truth = mctruthcol->at(i_intx);
      simb::MCNeutrino const& mc_neutrino = truth.GetNeutrino();
      simb::MCParticle const& mc_particle = mc_neutrino.Nu();

      if (mc_neutrino.CCNC() == 0 && fabs(mc_particle.PdgCode()) == 12) {
        ++fTotalTrueNonswapNue;
      }
      LOG_DEBUG("CAFMaker") << "end loop over mctruth";
    }

    for (size_t i_intx = 0; i_intx < singlemctruthcol->size(); ++i_intx) {
      LOG_DEBUG("outnNueCCOverlay") << "start loop";
      simb::MCTruth const& truth = singlemctruthcol->at(i_intx);
      simb::MCNeutrino const& mc_neutrino = truth.GetNeutrino();
      simb::MCParticle const& mc_particle = mc_neutrino.Nu();

      if (mc_neutrino.CCNC() == 0 && fabs(mc_particle.PdgCode()) == 12) {
        ++fTotalTrueSingleNue;
      }
      LOG_DEBUG("CAFMaker") << "end loop over mctruth";
    }
  }

  art::Handle<std::vector<rb::Cluster>> parentSlices;

  // in case we are running over MRE, make a vector of
  // parent slices, so we can pass it to backtracker
  // for filling the truth branch

  if (fParams.MRMode()) {
    GetByLabelStrict(evt, fParams.MRParentSliceLabel(), parentSlices);
  }

  // Get all of the hits from calhits

  art::Handle<std::vector<rb::CellHit>> allHits;
  GetByLabelStrict(evt, fParams.HitsLabel(), allHits);

  // We're gonna get a vector of NeutrinoEffPurs for allHits
  std::vector<cheat::NeutrinoWithIndex> allNus = bt->allMCTruth();
  std::vector<std::vector<cheat::NeutrinoEffPur>> sliceTruthTable;


  // Try to get POT spill info if this is real data
  art::Handle<sumdata::SpillData> spillPot;

  // treat as data spill in MC-dataspill overlay
  if (!fIsRealData) {
    GetByLabelStrict(evt, fParams.MCSpillDataLabel(), spillPot);
  } else {
    GetByLabelStrict(evt, fParams.NuMIBeamLabel(), spillPot);
  }

  // Get the EventQuality information
  art::Handle<sumdata::EventQuality> spillQual;
  GetByLabelStrict(evt, fParams.SpillQualLabel(), spillQual);

  // Get DAQ Header information
  art::Handle<dq::DAQEventSummary> daq;
  GetByLabelStrict(evt, fParams.DAQEventSummaryLabel(), daq);

  art::Handle<std::vector<rawdata::RawTrigger>> trigs;
  GetByLabelStrict(evt, fParams.RawDataLabel(), trigs);

  caf::SRSpill spill;
  FillSpillVars(spillPot, spillQual, daq, trigs, spill, evt);

  // Fill the spill tree once per event so that POT counting can be
  // done properly.
  // This tree will have an entry even for empty events with no physics
  // slices.
  caf::SRSpill* pspill = &spill;
  fSpillTree->SetBranchAddress("spill", &pspill);

  art::Handle<std::vector<cvntf::CVNCosmicFiltList> > cosmiccvn;
  if (fParams.FillCosmicCVN()){
    GetByLabelStrict(evt, fParams.CosmicCVNLabel(), cosmiccvn);
    std::vector<cvntf::CVNCosmicFiltList> cosmicList = *cosmiccvn;
    for (unsigned int iC = 0; iC < cosmicList[0].ListSize(); ++iC){
      cvntf::CVNCosmicFilt cvn = cosmicList[0].GetTimeSlice(iC);
      caf::SRCosmicCVN tcosmic;
      tcosmic.nHits      = cvn.nHits;
      tcosmic.timeWinMax = cvn.timeWinMax;
      tcosmic.timeWinMin = cvn.timeWinMin;
      tcosmic.numuVal    = cvn.numuVal;
      tcosmic.nueVal     = cvn.nueVal;
      tcosmic.nutauVal   = cvn.nutauVal;
      tcosmic.ncVal      = cvn.ncVal;
      tcosmic.cosmicVal  = cvn.cosmicVal;
      tcosmic.passSel    = cvn.passSel;
      spill.cosmiccvn.push_back(tcosmic);
    }
    spill.ncosmiccvn = spill.cosmiccvn.size();
  }
  fSpillTree->Fill();

  // Pull the MC cycle number from metadata
  //   before filling nuTree so cycle is filled
  // Assumes metadata module is being run...
  std::map<std::string, std::string> metadata =
      meta::MetadataManager::getInstance().GetMetadata();
  // Check it exists in the metadata, otherwise it defaults to -5
  if (metadata.count("simulated.cycle"))
    fCycle = std::stoi(metadata["simulated.cycle"]);
  if (metadata.count("simulated.batch"))
    fBatch = std::stoi(metadata["simulated.batch"]);

  // Fill the nuTree one per event, in sync with the spillTree
  if (fParams.FillNuTree()) {
    caf::SRSpillTruthBranch spillTruth;
    caf::SRSpillTruthBranch* pspillTruth = &spillTruth;
    fNuTree->SetBranchAddress("mc", &pspillTruth);

    for (const cheat::NeutrinoWithIndex& nuIdx :
         bt->allMCTruth()) {  // NeutrinoInteractions()){
      // AddMCTruthToVec wants one of these, but most of this info is only
      // relevant if the truth is associated with a slice.
      const cheat::NeutrinoEffPur ep = {nuIdx.neutrinoInt, 0, 0, 0, 0, 0, 0, 0};
      std::vector<cheat::TrackIDE> dummy;
      std::vector<SRNeutrino> nu;
      AddMCTruthToVec(&ep, dummy, dummy, dummy, dummy, &nu, evt, spill);
      spillTruth.nu = nu[0];

      art::Ptr<simb::MCTruth> mctruth = nuIdx.neutrinoInt;
      const std::vector<art::Ptr<simb::MCTruth>> mctv(1, mctruth);

      art::FindManyP<fxwgt::FluxWeights> fmFLXRW =
          FindManyPStrict<fxwgt::FluxWeights>(mctv, evt,
                                              fParams.FluxReweightLabel());
      fxwgt::FluxWeights flxrw;
      if (GetAssociatedProduct(fmFLXRW, 0, flxrw))
        spillTruth.nu.rwgt.ppfx = FluxReweights(flxrw);
      else
        spillTruth.nu.rwgt.ppfx.setDefault();

      art::FindManyP<rwgt::GENIEReweightTable> fmRW =
          FindManyPStrict<rwgt::GENIEReweightTable>(mctv, evt,
                                                    fParams.ReweightLabel());
      rwgt::GENIEReweightTable rw;
      if (GetAssociatedProduct(fmRW, 0, rw))
        spillTruth.nu.rwgt.genie = GenieReweightTable(rw);

      fNuTree->Fill();
      spilltruthcol->push_back(spillTruth);
    }  // end for nuIdx
  }    // end if FillNuTree

  // Applying an event filter but we need the unfiltered events for POT/LT accounting.
  // Thus we need the same version of cafmaker to run over the filtered events and !filtered.
  // However, the point is to save time so we don't want slicer or calhit on the !filtered stream.
  if(fParams.ApplyingFilter() && !slices.isValid()) {
    *spillcol = spill;

    evt.put(std::move(srcol));
    evt.put(std::move(spillcol));
    if (fParams.FillNuTree()) evt.put(std::move(spilltruthcol));
    evt.put(std::move(srAssn));
    
    return;
  }

  art::FindManyP<murem::MRCCParent> fmMrccPar =
      FindManyPStrict<murem::MRCCParent>(slices, evt,
                                         fParams.MRCCParentLabel());
  if (fParams.MRMode()) {
    std::vector<std::vector<cheat::NeutrinoEffPur>> temptable  // I know :)
        = bt->SlicesToMCTruthsTable(*parentSlices);
    // some mre slices may not have matched parents
    // so sort the table so that nth entry in table
    // is the truth for the parent of the nth mre
    // slice. If nth mre slice has no parent,
    // make an empty effpur entry.
    for (unsigned int isli = 0; isli < slices->size(); isli++) {
      murem::MRCCParent parent;
      if (GetAssociatedProduct(fmMrccPar, isli, parent)) {
        int idx = parent.ParentIndex();
        sliceTruthTable.push_back(temptable[idx]);
      } else {
        art::Ptr<simb::MCTruth> null;
        cheat::NeutrinoEffPur effPur =
            cheat::NeutrinoEffPur{null, -1, -1, 0, -1, -1, -1, -1};
        std::vector<cheat::NeutrinoEffPur> emptyEffPur(allNus.size(), effPur);
        sliceTruthTable.push_back(emptyEffPur);
      }
    }  // end loop over mr slices
  }    // end of mr mode
  else
    sliceTruthTable = bt->SlicesToMCTruthsTable(*slices);

  // Find favorite neutrino IDs using several metrics.  Total energy
  // is the one we will use, for the favorite neutrino, but the others
  // will be stored.
  std::vector<int> faveIdsEff = bt->SliceToOrderedNuIds(
      allNus, sliceTruthTable, cheat::EffMetric, cheat::EffMetric);
  std::vector<int> faveIdsEnergy = bt->SliceToOrderedNuIds(
      allNus, sliceTruthTable, cheat::EnergyMetric, cheat::EnergyMetric);
  std::vector<int> faveIdsPur = bt->SliceToOrderedNuIds(
      allNus, sliceTruthTable, cheat::PurMetric, cheat::PurMetric);
  std::vector<int> faveIdsEffPur = bt->SliceToOrderedNuIds(
      allNus, sliceTruthTable, cheat::EffPurMetric, cheat::EffPurMetric);
  std::vector<int> faveIdsEffThenPur = bt->SliceToOrderedNuIds(
      allNus, sliceTruthTable, cheat::EffMetric, cheat::PurMetric);

  // Get all truth track hits in the event
  std::vector<cheat::TrackIDE> allTracks;
  std::vector<cheat::TrackIDE> sliceTracks;
  std::vector<cheat::TrackIDE> allTracksBirks;
  std::vector<cheat::TrackIDE> sliceTracksBirks;
  if (bt->HaveTruthInfo()) allTracks = bt->HitsToTrackIDE(*allHits);
  if (bt->HaveTruthInfo()) allTracksBirks = bt->HitsToTrackIDE(*allHits, true);

  // Get all association finder objects
  auto fmDiscrete2d =
      FindManyPStrict<rb::Track>(slices, evt, fParams.DiscreteTrack2dLabel());
  auto fmDiscrete =
      FindManyPStrict<rb::Track>(slices, evt, fParams.DiscreteTrackLabel());
  auto fmKalman =
      FindManyPStrict<rb::Track>(slices, evt, fParams.KalmanTrackLabel());
  auto fmCosmic =
      FindManyPStrict<rb::Track>(slices, evt, fParams.CosmicTrackLabel());
  auto fmWindow =
      FindManyPStrict<rb::Track>(slices, evt, fParams.WindowTrackLabel());
  auto fmElastic =
      FindManyPStrict<rb::Vertex>(slices, evt, fParams.ElasticArmsLabel());
  auto fmHoughVertex =
      FindManyPStrict<rb::Vertex>(slices, evt, fParams.HoughVertexLabel());
  auto fmVertexDT =
      FindManyPStrict<rb::Vertex>(slices, evt, fParams.VertexDTLabel());
  auto fmSlcME = FindManyPStrict<me::SlcME>(slices, evt, fParams.MELabel());
  auto fmJmShower =
      FindManyPStrict<jmshower::JMShower>(slices, evt, fParams.JMShowerLabel());
  auto fmNueSel =
      FindManyPStrict<jmshower::EID>(slices, evt, fParams.NueSelLabel());
  art::FindOneP<slid::EventLID> foLID(
      slices, evt, fParams.LIDLabel());  // one lid object per slice
  art::FindOneP<slid::SliceLID> foSlcLID(
      slices, evt, fParams.SlcLIDLabel());  // one lid object per slice
  art::FindMany<slid::ShowerLID> fmShwLID(slices, evt, fParams.LIDLabel());
  art::FindMany<slid::ShowerPID> fmShwPID(slices, evt, fParams.SPIDLabel());
  auto fmSliceLID = FindManyPStrict<SliceLID::Prediction>(
      slices, evt, fParams.SliceLIDLabel());
  auto fmNumuEnergy =
      FindManyPStrict<numue::NumuE>(slices, evt, fParams.NumuEnergyLabel());
  auto fmNumuLSTMEnergy = FindManyPStrict<LSTME::LSTMEnergy>(
      slices, evt, fParams.NumuLSTMEnergyLabel());
  auto fmNuePresel = FindManyPStrict<presel::PreselObj>(
      slices, evt, fParams.NuePreselLabel());
  auto fmRockPresel = FindManyPStrict<presel::PreselObj>(
      slices, evt, fParams.RockPreselLabel());
  auto fmVeto = FindManyPStrict<presel::Veto>(slices, evt, fParams.VetoLabel());
  auto fmNueVeto =
      FindManyPStrict<presel::Veto>(slices, evt, fParams.NueVetoLabel());
  auto fmLem =
      FindManyPStrict<lem::PIDDetails>(slices, evt, fParams.LemLabel());
  auto fmRvp = FindManyPStrict<rvp::RVP>(slices, evt, fParams.RvpLabel());
  auto fmXnue = FindManyPStrict<xnue::Xnue>(slices, evt, fParams.XnueLabel());
  auto fmCosRej =
      FindManyPStrict<cosrej::CosRejObj>(slices, evt, fParams.CosRejLabel());
  auto fmNueCosRej =
      FindManyPStrict<cosrej::NueCosRej>(slices, evt, fParams.NueCosRejLabel());
  auto fmNCCosRej = FindManyPStrict<ncid::NCCosRej>(
      slices, evt, fParams.NCCosRejLabel());  /// NC Cosmic Rejection
  auto fmNCPi0BkgRej = FindManyPStrict<ncpi0::NCPi0BkgRej>(
      slices, evt, fParams.NCPi0BkgRejLabel());  /// NC Pi0 Bkg Rejection

  auto fmNuonE = 
      FindManyPStrict<cvn::Result>(slices, evt, fParams.NuonELabel());

  auto fmCVN_LoosePreselPtP =
      FindManyPStrict<cvn::Result>(slices, evt, fParams.CVNLabelLoosePreselPtP());
  auto fmCVN_OldPresel =
      FindManyPStrict<cvn::Result>(slices, evt, fParams.CVNLabelOldPresel());
  auto fmCVN_NoCosmics =
      FindManyPStrict<cvn::Result>(slices, evt, fParams.CVNLabelNoCosmics());

  auto fmWS =
    FindManyPStrict<cvn::Result>(slices, evt, fParams.WSCNNLabel());

  auto fmCVNFeatures =
    FindManyPStrict<cvn::Features>(slices, evt, fParams.CVNLabelLoosePreselPtP());
  auto fmCVNPixelMaps =
    FindManyPStrict<cvn::PixelMap>(slices, evt, fParams.CVNPixelMapLabel());
  auto fmCVNTrainingData = 
    FindManyPStrict<cvn::TrainingData> (slices, evt, fParams.CVNTrainingDataLabel());
  auto fmEvtRegCVN =
      FindManyPStrict<cvn::RegNuResult>(slices, evt, fParams.RegCVNLabel());
  auto fmHadRegCVN =
      FindManyPStrict<cvn::RegHadronResult>(slices, evt, fParams.RegCVNLabel());


  art::Handle<sumdata::SpillData> singlepot;
  if (fParams.IsSinglesOverlay()) {
    GetByLabelStrict(evt, fParams.SingleMixerLabel(), singlepot);
    if (!singlepot.failedToGet()) {
      float totsinglepot = singlepot->spillpot;
      fTotalSinglePOT += totsinglepot;
    }
  }

  int nNonNoise = 0;
  for (const rb::Cluster& slice : *slices) {
    if (slice.IsNoise() || slice.NCell() == 0) continue;
    ++nNonNoise;
  }
  for (unsigned int sliceId = 0; sliceId < slices->size(); ++sliceId) {
    const rb::Cluster& slice = (*slices)[sliceId];

    if (slice.IsNoise() || slice.NCell() == 0) continue;
    // Because we don't care about the noise slice and slices with no hits.
    StandardRecord rec;
    StandardRecord* prec = &rec;  // TTree wants a pointer-to-pointer
    fRecTree->SetBranchAddress("rec", &prec);

    if (bt->HaveTruthInfo()) sliceTracks = bt->HitsToTrackIDE(slice.AllCells());
    if (bt->HaveTruthInfo())
      sliceTracksBirks = bt->HitsToTrackIDE(slice.AllCells(), true);
    //#######################################################
    // Fill header.
    //#######################################################
    // Get metadata information for header
    unsigned int run = evt.run();
    unsigned int subrun = evt.subRun();
    unsigned int spillNum = evt.id().event();

    art::Timestamp ts = evt.time();
    const time_t timeSec = ts.timeHigh();

    struct tm* timeStruct = localtime(&timeSec);

    unsigned int maskstatus = 0;
    unsigned int dibmask = fMask;

    if (rh->GoodDiBlockMask() == ((1 << 16) - 1))
      maskstatus = 0;
    else if (rh->GoodDiBlockMask() == ((1 << 17) - 1))
      maskstatus = 2;
    else
      maskstatus = 1;

    unsigned int dibfirst = 0;
    unsigned int diblast = 0;
    if (dibmask) {
      int iD;
      iD = 0;
      while (!((dibmask >> iD) & 1)) iD++;
      dibfirst = iD + 1;
      iD = 0;
      while (dibmask >> iD) iD++;
      diblast = iD;
    }

    unsigned int nbadchan = bc->NBadInSubRun(subrun);
    unsigned int ntotchan = rh->NAnalysisChannels();

    rec.hdr = SRHeader();

    rec.hdr.run = run;
    rec.hdr.subrun = subrun;
    rec.hdr.cycle = fCycle;
    rec.hdr.batch = fBatch;
    rec.hdr.evt = spillNum;
    rec.hdr.subevt = sliceId;
    rec.hdr.ismc = !fIsRealData;
    rec.hdr.det = fDet;
    rec.hdr.blind = 0;
    rec.hdr.filt = rb::IsFiltered(evt, slices, sliceId);

    rec.hdr.subevtstarttime = slice.MinTNS();
    rec.hdr.subevtendtime = slice.MaxTNS();
    rec.hdr.subevtmeantime = slice.MeanTNS();
    rec.hdr.unixtime = timeSec;

    rec.hdr.dibfirst = dibfirst;
    rec.hdr.diblast = diblast;
    rec.hdr.dibmask = dibmask;
    rec.hdr.maskstatus = maskstatus;

    rec.hdr.year = timeStruct->tm_year + 1900;
    rec.hdr.month = timeStruct->tm_mon + 1;
    rec.hdr.day = timeStruct->tm_mday;
    rec.hdr.doy = timeStruct->tm_yday + 1;
    rec.hdr.hour = timeStruct->tm_hour + 1;
    rec.hdr.minute = timeStruct->tm_min;
    rec.hdr.second = timeStruct->tm_sec;

    rec.hdr.nbadchan = nbadchan;
    rec.hdr.ntotchan = ntotchan;

    rec.hdr.gain = rh->DetGainSetting();
    rec.hdr.finetiming = rh->DetFineTimingSetting();

    //#######################################################
    // Fill Spill info.
    //#######################################################
    // Store info on NuMI beam quality and if the pot as good or bad.
    // This is meant to be used to refine a good spill, not for
    // POT counting since events have multiple slices and empty
    // spill are missing from this tree.  Use the friend tree for
    // POT counting.
    rec.spill = spill;

    //trim list of Cosmic CVN scores down to only the time window the
    //physics slice matches
    for (int iC = 0; iC < (int)rec.spill.cosmiccvn.size(); ++iC){

      if ((slice.MeanTNS() > rec.spill.cosmiccvn[iC].timeWinMin) &&
          (slice.MeanTNS() < rec.spill.cosmiccvn[iC].timeWinMax)) continue;
      rec.spill.cosmiccvn.erase(rec.spill.cosmiccvn.begin()+iC);
      --iC;

    }
    //if more then one time window matches the slice, keep the one with the 
    //lowest cosmic CVN score
    if (rec.spill.cosmiccvn.size() > 1){
      SRCosmicCVN tcos = rec.spill.cosmiccvn[0];
      for (unsigned int iC = 1; iC < rec.spill.cosmiccvn.size(); ++iC){
        if (rec.spill.cosmiccvn[iC].cosmicVal < tcos.cosmicVal) 
          tcos = rec.spill.cosmiccvn[iC];
      }
      rec.spill.cosmiccvn.clear();
      rec.spill.cosmiccvn.push_back(tcos);
    }
    rec.spill.ncosmiccvn = rec.spill.cosmiccvn.size();


    //#######################################################
    // Add slice info.
    //#######################################################
    FillSliceVars(slice, rec.slc);
    rec.slc.nnonnoise = nNonNoise;
    // find the slice with minimum mean time distance from this slice
    float mintimediff = 1e7;  // initialize to a big number
    signed int mintimejsli = -1;
    for (unsigned int jsli = 0; jsli < slices->size(); jsli++) {
      const rb::Cluster& slicej = (*slices)[jsli];
      if (jsli == sliceId || slicej.IsNoise()) continue;
      float diff = fabs(slicej.MeanTNS() - slice.MeanTNS());
      if (diff < mintimediff) {
        mintimejsli = jsli;
      }
      mintimediff = std::min(mintimediff, diff);
    }
    if (mintimejsli != -1) {
      const rb::Cluster& slicejmin = (*slices)[mintimejsli];
      rec.slc.closestslicetime = slicejmin.MeanTNS() - slice.MeanTNS();
      rec.slc.closestslicenhit = slicejmin.NCell();
      rec.slc.closestslicecalE = slicejmin.CalorimetricEnergy();
      TVector3 minp, maxp;
      minp = slicejmin.MinXYZ();
      maxp = slicejmin.MaxXYZ();
      rec.slc.closestsliceminfromtop = livegeom->DistToTop(maxp);
      rec.slc.closestsliceminfrombottom = livegeom->DistToBottom(minp);
      rec.slc.closestsliceminfromfront = livegeom->DistToFront(minp);
      rec.slc.closestsliceminfromback = livegeom->DistToBack(maxp);
      rec.slc.closestsliceminfromeast = livegeom->DistToEdgeXMinus(minp);
      rec.slc.closestsliceminfromwest = livegeom->DistToEdgeXPlus(maxp);
      float mincelldistX = 55555;
      float mincelldistY = 55555;
      float mincelldist = 55555;
      float cellD = 2. * geom->Plane(0)->Cell(0)->HalfD();
      float cellW = 2. * geom->Plane(0)->Cell(0)->HalfW();
      for (unsigned int hitIdxX = 0; hitIdxX < slice.NXCell(); ++hitIdxX) {
        const art::Ptr<rb::CellHit>& chitX = slice.XCell(hitIdxX);
        for (unsigned int minhitIdxX = 0; minhitIdxX < slicejmin.NXCell();
             ++minhitIdxX) {
          const art::Ptr<rb::CellHit>& minchitX = slicejmin.XCell(minhitIdxX);
          float dPlaneX = (float)chitX->Plane() - (float)minchitX->Plane();
          float dCellX = (float)chitX->Cell() - (float)minchitX->Cell();
          float distX = sqrt((dPlaneX * cellD) * (dPlaneX * cellD) +
                             (dCellX * cellW) * (dCellX * cellW));
          if (distX < mincelldistX) mincelldistX = distX;
        }
      }
      for (unsigned int hitIdxY = 0; hitIdxY < slice.NYCell(); ++hitIdxY) {
        const art::Ptr<rb::CellHit>& chitY = slice.YCell(hitIdxY);
        for (unsigned int minhitIdxY = 0; minhitIdxY < slicejmin.NYCell();
             ++minhitIdxY) {
          const art::Ptr<rb::CellHit>& minchitY = slicejmin.YCell(minhitIdxY);
          float dPlaneY = (float)chitY->Plane() - (float)minchitY->Plane();
          float dCellY = (float)chitY->Cell() - (float)minchitY->Cell();
          float distY = sqrt((dPlaneY * cellD) * (dPlaneY * cellD) +
                             (dCellY * cellW) * (dCellY * cellW));
          if (distY < mincelldistY) mincelldistY = distY;
        }
      }
      if (mincelldistX < mincelldistY) {
        mincelldist = mincelldistX;
      } else {
        mincelldist = mincelldistY;
      }
      rec.slc.closestslicemindist = mincelldist;
    } else {  // protection when mintimejsli == -1
      rec.slc.closestslicetime = -5;
      rec.slc.closestslicenhit = 0;
      rec.slc.closestslicecalE = -5;
      rec.slc.closestsliceminfromtop = -5;
      rec.slc.closestsliceminfrombottom = -5;
      rec.slc.closestsliceminfromfront = -5;
      rec.slc.closestsliceminfromback = -5;
      rec.slc.closestsliceminfromeast = -5;
      rec.slc.closestsliceminfromwest = -5;
      rec.slc.closestslicemindist = 55555;
    }
    if (fDet == kFARDET) {
      rec.sel.contain.nplanestofront = rec.slc.firstplane - (dibfirst - 1) * 64;
      rec.sel.contain.nplanestoback = ((diblast)*64) - 1 - rec.slc.lastplane;
    }

    if (fDet == kNEARDET) {
      rec.sel.contain.nplanestofront = rec.slc.firstplane;
      rec.sel.contain.nplanestoback = 214 - rec.slc.lastplane;
    }

    //#######################################################
    // Adding reconstructed objects.
    //#######################################################

    // Get all the tracks associated with this slice from discrete
    // tracker. This depends on the findMany object created above.
    if (fmDiscrete.isValid()) {
      std::vector<art::Ptr<rb::Track>> discTracks = fmDiscrete.at(sliceId);
      // Sort the tracks in descending length order
      std::sort(discTracks.begin(), discTracks.end(), sortRBTrackLength);
      rec.trk.discrete.tracks.resize(discTracks.size());
      for (size_t trackId = 0; trackId < discTracks.size(); ++trackId) {
        FillTrackVars((*discTracks[trackId]), rec.trk.discrete.tracks[trackId],
                      *slices, allTracks, sliceTracks, allTracksBirks,
                      sliceTracksBirks, sliceId);
      }
    }

    //      if(fmDiscrete2d.isValid())
    //      {
    //        std::vector< art::Ptr<rb::Track> > discTracks =
    //        fmDiscrete2d.at(sliceId);
    // Sort the tracks in descending length order
    //        std::sort(discTracks.begin(), discTracks.end(),
    //        sortRBTrackLength);
    //        rec.trk.discrete.tracks2d.resize(discTracks.size());
    //        for(size_t trackId = 0; trackId < discTracks.size(); ++trackId){
    //          FillTrackVars(*discTracks[trackId],
    //          rec.trk.discrete.tracks2d[trackId], *slices, sliceId);
    //        }
    //      }

    if (fmCosmic.isValid()) {
      std::vector<art::Ptr<rb::Track>> cosTracks = fmCosmic.at(sliceId);
      art::FindManyP<me::TrkME> fmCosmicTrkME(cosTracks, evt, fParams.MELabel());

      // Sort the tracks in descending length order
      std::sort(cosTracks.begin(), cosTracks.end(), sortRBTrackLength);

      rec.trk.cosmic.tracks.reserve(cosTracks.size());
      rec.trk.cosmic.tracks.resize(cosTracks.size());
      for (size_t trackId = 0; trackId < cosTracks.size(); ++trackId) {
        FillTrackVars((*cosTracks[trackId]), rec.trk.cosmic.tracks[trackId],
                      *slices, allTracks, sliceTracks, allTracksBirks,
                      sliceTracksBirks, sliceId);

        // Start of MR Shower
        // Fill the MR Shower informations
        art::FindManyP<rb::Cluster> fmDiF(cosTracks, evt,
                                          fParams.DiFShowerClusterLabel());
        if (fmDiF.isValid()) {
          std::vector<art::Ptr<rb::Cluster>> difs = fmDiF.at(trackId);
          for (unsigned int difId = 0; difId < difs.size(); ++difId) {
            SRMRProperties srMR;  // we're going to accumulate info into this
            art::Ptr<rb::Cluster> dif = difs[difId];
            // Filling DiF Cluster variables
            FillDiFVars(*dif, *slices, allTracks, sliceTracks, allTracksBirks,
                        sliceTracksBirks, sliceId, srMR);

            // -- CVN - Loose Presel PtP Cut
            art::FindManyP<cvn::Result> fmDiFCVNResult_LoosePreselPtP(
                {dif}, evt, fParams.DiFShowerCVNLabelLoosePreselPtP());
            cvn::Result cvnDiFResult_LoosePreselPtP;
            if (GetAssociatedProduct(fmDiFCVNResult_LoosePreselPtP, difId,
                                     cvnDiFResult_LoosePreselPtP)) {
              FillCVNResultVars(cvnDiFResult_LoosePreselPtP, srMR.cvnloosepreselptp);
            } else {
              srMR.cvnloosepreselptp.setDefault();
            }
            // -- CVN - Old Presel
            art::FindManyP<cvn::Result> fmDiFCVNResult_OldPresel(
                {dif}, evt, fParams.DiFShowerCVNLabelOldPresel());
            cvn::Result cvnDiFResult_OldPresel;
            if (GetAssociatedProduct(fmDiFCVNResult_OldPresel, difId,
                                     cvnDiFResult_OldPresel)) {
              FillCVNResultVars(cvnDiFResult_OldPresel, srMR.cvnoldpresel);
            } else {
              srMR.cvnoldpresel.setDefault();
            }
            // -- CVN - No Cosmics used in training
            art::FindManyP<cvn::Result> fmDiFCVNResult_NoCosmics(
                {dif}, evt, fParams.DiFShowerCVNLabelNoCosmics());
            cvn::Result cvnDiFResult_NoCosmics;
            if (GetAssociatedProduct(fmDiFCVNResult_NoCosmics, difId,
                                     cvnDiFResult_NoCosmics)) {
              FillCVNResultVars(cvnDiFResult_NoCosmics, srMR.cvnnocosmics, true);
            } else {
              srMR.cvnnocosmics.setDefault();
            }

            // Filling DiF LID
            art::FindManyP<rb::Shower> fmDiFShw({dif}, evt,
                                                fParams.DiFShowerShwLabel());
            if (fmDiFShw.isValid()) {
              std::vector<art::Ptr<rb::Shower>> difshws = fmDiFShw.at(difId);
              for (art::Ptr<rb::Shower> difshw : difshws) {
                if (difshw) {
                  FillDiFShowerVars(*difshw, *slices, allTracks, sliceTracks,
                                    allTracksBirks, sliceTracksBirks, sliceId,
                                    srMR);
                }  // end for if difshw
              }    // end for loop difshw
            }      // end if Shw assn DiF valid

            art::FindManyP<slid::ShowerLID> fmDiFLID(
                {dif}, evt, fParams.DiFShowerLIDLabel());
            if (fmDiFLID.isValid()) {
              std::vector<art::Ptr<slid::ShowerLID>> diflids =
                  fmDiFLID.at(difId);
              for (art::Ptr<slid::ShowerLID> diflid : diflids) {
                if (diflid) {
                  FillSlidVars(*diflid, srMR.lid);
                }  // end for if diflid
              }    // end for loop diflid
            }      // end if LID assn DiF valid

            // Now the DiF cluster is ready
            rec.trk.cosmic.tracks[trackId].mrdif.push_back(srMR);
          }  // end for dif cluster loop
        }    // end if Track assn with DiF is valid
        // End of MR Shower

        if (fmCosmicTrkME.isValid()) {
          // Get MichelE from MEFinder
          std::vector<art::Ptr<me::TrkME>> michels = fmCosmicTrkME.at(trackId);
          for (art::Ptr<me::TrkME> michel : michels) {
            AddTrkMEToVec(*michel, &rec.me.trkcosmic,
                          rec.trk.cosmic.tracks.back(), *slices, allTracks,
                          sliceTracks, allTracksBirks, sliceTracksBirks,
                          sliceId);
          }
        }  // end

      }  // end for trackId
    }

    // Do the same for windowtrack tracks
    if (fmWindow.isValid()) {
      std::vector<art::Ptr<rb::Track>> winTracks = fmWindow.at(sliceId);
      // Sort the tracks in descending length order
      std::sort(winTracks.begin(), winTracks.end(), sortRBTrackLength);
      rec.trk.window.tracks.reserve(winTracks.size());
      rec.trk.window.tracks.resize(winTracks.size());
      for (size_t trackId = 0; trackId < winTracks.size(); ++trackId) {
        FillTrackVars((*winTracks[trackId]), rec.trk.window.tracks[trackId],
                      *slices, allTracks, sliceTracks, allTracksBirks,
                      sliceTracksBirks, sliceId);
      }  // end for trackId
    }

    // Now do the same for kalmanTrack
    if (fmKalman.isValid()) {
      std::vector<art::Ptr<rb::Track>> kalTracks = fmKalman.at(sliceId);
      // Sort the tracks in descending length order to start
      std::sort(kalTracks.begin(), kalTracks.end(), sortRBTrackLength);
      rec.trk.kalman.idxlongest = 0;

      art::FindManyP<remid::ReMId> fmRemid(kalTracks, evt,
                                           fParams.RemidLabel());

      art::FindManyP<cosrej::TrkCntObj> fmTrkCnt =
          FindManyPStrict<cosrej::TrkCntObj>(kalTracks, evt,
                                             fParams.CosRejLabel());
      art::FindManyP<me::TrkME> fmKalmanME(kalTracks, evt, fParams.MELabel());
      art::FindManyP<rb::Energy> fmKalOverlapE(
          kalTracks, evt,fParams.OverlapEKalLabel());

      art::FindManyP<trackinfo::TrackInfoObj> fmTrkInfo = 
        FindManyPStrict<trackinfo::TrackInfoObj>(kalTracks, evt, 
                                                 fParams.TrackInfoKalman());

      // Give SRRemid the pid values of track with largest PID
      unsigned int bestPidIdx =
          remid::HighestPIDTrack(kalTracks, fParams.RemidLabel(), evt);

      unsigned int bestPidIdx3d = 999;



      for (size_t trackId = 0; trackId < kalTracks.size(); ++trackId) {
        if (kalTracks[trackId]->View() != geo::kXorY) {
          rec.trk.kalman.tracks2d.push_back(SRTrack());

          // Fill the Track Vars with overlapping E if that info is present,
          // otherwise use the default FillTrackVars function.
          if (fmKalOverlapE.isValid()) {
            std::vector<art::Ptr<rb::Energy>> overlapEs =
                fmKalOverlapE.at(trackId);
            if (!overlapEs.empty()) {
              FillTrackVarsWithOverlapE((*kalTracks[trackId]), (*overlapEs[0]),
                                        rec.trk.kalman.tracks2d.back(), *slices,
                                        allTracks, sliceTracks, allTracksBirks,
                                        sliceTracksBirks, sliceId);
            } else {
              FillTrackVars((*kalTracks[trackId]),
                            rec.trk.kalman.tracks2d.back(), *slices, allTracks,
                            sliceTracks, allTracksBirks, sliceTracksBirks,
                            sliceId);
            }
          } else {
            FillTrackVars((*kalTracks[trackId]), rec.trk.kalman.tracks2d.back(),
                          *slices, allTracks, sliceTracks, allTracksBirks,
                          sliceTracksBirks, sliceId);
          }
          if(fmTrkInfo.isValid()){
            std::vector<art::Ptr<trackinfo::TrackInfoObj> >  trkinfos = fmTrkInfo.at(trackId);
            //if(!trkinfos.empty())
            //  FillTrackInfoVars(*trkinfos[0], rec.trk.kalman.tracks.back());
          }
        } else {
          rec.trk.kalman.tracks.push_back(SRKalmanTrack());

          // Fill the Track Vars with overlapping E if that info is present,
          // otherwise use the default FillTrackVars function.
          if (fmKalOverlapE.isValid()) {
            std::vector<art::Ptr<rb::Energy>> overlapEs =
                fmKalOverlapE.at(trackId);
            if (!overlapEs.empty()) {
              FillTrackVarsWithOverlapE((*kalTracks[trackId]), (*overlapEs[0]),
                                        rec.trk.kalman.tracks.back(), *slices,
                                        allTracks, sliceTracks, allTracksBirks,
                                        sliceTracksBirks, sliceId);
            } else {
              FillTrackVars((*kalTracks[trackId]), rec.trk.kalman.tracks.back(),
                            *slices, allTracks, sliceTracks, allTracksBirks,
                            sliceTracksBirks, sliceId);
            }
          } else {
            FillTrackVars((*kalTracks[trackId]), rec.trk.kalman.tracks.back(),
                          *slices, allTracks, sliceTracks, allTracksBirks,
                          sliceTracksBirks, sliceId);
          }

          if(fmTrkInfo.isValid()){
            std::vector<art::Ptr<trackinfo::TrackInfoObj> >  trkinfos = fmTrkInfo.at(trackId);
            if(!trkinfos.empty())
              FillTrackInfoVars(*trkinfos[0], rec.trk.kalman.tracks.back());
          }

          if (trackId == bestPidIdx)
            bestPidIdx3d = rec.trk.kalman.tracks.size() - 1;

          if (fmRemid.isValid()) {
            // Get MuID from ReMId
            std::vector<art::Ptr<remid::ReMId>> remids = fmRemid.at(trackId);

            if (!remids.empty()) {
              FillMuIdVars(*remids[0], rec.trk.kalman.tracks.back());
            } else {
              std::cout << "Missing Association Error: Failed to find a ReMId "
                           "object with this track."
                        << std::endl;
              abort();
            }
          }

          if (fmTrkCnt.isValid()) {
            // Get track containment info from TrkCntObj
            std::vector<art::Ptr<cosrej::TrkCntObj>> trkcnts =
                fmTrkCnt.at(trackId);

            if (!trkcnts.empty())
              FillTrackContainmentVars(*trkcnts[0],
                                       rec.trk.kalman.tracks.back());
          }

          if (fmKalmanME.isValid()) {
            // Get MichelE from MEFinder
            std::vector<art::Ptr<me::TrkME>> michels = fmKalmanME.at(trackId);

            for (art::Ptr<me::TrkME> michel : michels) {
              AddTrkMEToVec(*michel, &rec.me.trkkalman,
                            rec.trk.kalman.tracks.back(), *slices, allTracks,
                            sliceTracks, allTracksBirks, sliceTracksBirks,
                            sliceId);
            }
          }  // end
        }
      }  // end for trackId

      if (!rec.trk.kalman.tracks.empty() && bestPidIdx3d != 999) {
        CopyRemidVars(rec.trk.kalman.tracks[bestPidIdx3d], rec.sel.remid);
        std::vector<art::Ptr<remid::ReMId>> remids = fmRemid.at(bestPidIdx3d);
        rec.trk.kalman.idxremid = bestPidIdx3d;
      } else {
        rec.sel.remid.setDefault();
      }

    } else {
      rec.sel.remid.setDefault();
    }

    std::sort(rec.trk.kalman.tracks.begin(), rec.trk.kalman.tracks.end(),
              sortRemid);
    rec.trk.kalman.idxremid = 0;
    unsigned int tempidx = 0, templen = 0;
    for (size_t trackId = 0; trackId < rec.trk.kalman.tracks.size();
         ++trackId) {
      if (rec.trk.kalman.tracks[trackId].len > templen) {
        tempidx = trackId;
        templen = rec.trk.kalman.tracks[trackId].len;
      }
    }

    rec.trk.kalman.idxlongest = tempidx;
    rec.trk.kalman.fillSizes();
    rec.trk.discrete.fillSizes();
    rec.trk.window.fillSizes();
    rec.trk.cosmic.fillSizes();

    // pre-fill nue energy variables with default
    rec.energy.nue.setDefault();
    rec.energy.nue.lid.setDefault();

    // get vertices from ElasticArms
    if (fmElastic.isValid()) {


      std::vector<art::Ptr<rb::Vertex>> elastics = fmElastic.at(sliceId);
      if(elastics.size() > 0){ //this should always be size 1 or 0
        
        rec.vtx.elastic = SRElastic();
        int vtxId = 0; //always check the 0th entry since there should only be 1.
        //will remove this and replace all vtxId with 0 once we are certain the 
        //code all works. 
        SRElastic& srVtx = rec.vtx.elastic;
        srVtx.time = elastics[vtxId]->GetT();
        srVtx.vtx = elastics[vtxId]->GetXYZ();
        srVtx.IsValid = true;
        // for each vertex (currently only one elastic arms ever), look for
        // prongs (currently only fuzzyk)
        art::FindManyP<rb::Prong> fmFuzzyProng2D = FindManyPStrict<rb::Prong>(
            elastics, evt, fParams.FuzzyK2DLabel());  // 2D prongs first
        if (fmFuzzyProng2D.isValid()) {
          std::vector<art::Ptr<rb::Prong>> prongs = fmFuzzyProng2D.at(vtxId);
          std::sort(prongs.begin(), prongs.end(),
                    sortRBProngLength);  // 2D prongs are sorted by length
          art::FindManyP<rb::WeightedProng> fmPngPngWt(
              prongs, evt,fParams.FuzzyKWeight2DLabel());

          for (unsigned int iPng = 0; iPng < prongs.size(); ++iPng) {
            art::Ptr<rb::Prong> prong = prongs[iPng];
            srVtx.fuzzyk.png2d.push_back(SRProng());
            FillProngVars(*prong, srVtx.fuzzyk.png2d.back(), *slices, allTracks,
                          sliceTracks, allTracksBirks, sliceTracksBirks,
                          sliceId);
            art::FindManyP<rb::PID> fmCVNParticleResult(
                prongs, evt, fParams.CVNParticle2DLabel());
            if (fmCVNParticleResult.isValid()) {
              std::vector<art::Ptr<rb::PID>> cvnparts =
                 fmCVNParticleResult.at(iPng);
              if(!cvnparts.empty()){
                FillCVNParticleResultVars(cvnparts, srVtx.fuzzyk.png2d.back().cvnpart);
                }
            }
            art::FindManyP<cvn::PixelMap> fmCVNProngPixelMaps(
                prongs, evt, fParams.CVNPixelMapLabel());
            art::FindManyP<cvn::ProngTrainingData> fmCVNProngTrainingData(
                prongs, evt, fParams.CVNTrainingDataLabel());

            if (fmCVNProngPixelMaps.isValid() && fParams.FillPixelMaps()) {
              std::vector<art::Ptr<cvn::PixelMap>> cvnProngPixelMaps =
                  fmCVNProngPixelMaps.at(iPng);
              for (unsigned int ip = 0; ip < cvnProngPixelMaps.size(); ++ip) {
                srVtx.fuzzyk.png2d.back().cvnmaps.push_back(SRPixelMap());
                SRPixelMap& pmap = srVtx.fuzzyk.png2d.back().cvnmaps.back();
                FillCVNPixelMaps(*cvnProngPixelMaps[ip], pmap,
                                 fParams.UseGeVPixelMaps());
              }
            }

            cvn::ProngTrainingData cvnProngTrainingData;
            if (GetAssociatedProduct(fmCVNProngTrainingData, iPng,
                                     cvnProngTrainingData) &&
                fParams.FillTrainingData()) {
              srVtx.fuzzyk.png2d.back().prongtrainingdata.push_back(
                  SRProngTrainingData());
              SRProngTrainingData& srProngTrainingData =
                  srVtx.fuzzyk.png2d.back().prongtrainingdata.back();
              FillCVNProngTrainingData(cvnProngTrainingData,
                                       srProngTrainingData);
            }

            if (fmPngPngWt.isValid()) {
              std::vector<art::Ptr<rb::WeightedProng>> prongWt =
                  fmPngPngWt.at(iPng);
              srVtx.fuzzyk.png2d.back().weightedCalE =
                  prongWt[0]->CalorimetricEnergy();
            }

          }  // end loop over 2D prongs
        }    // end if 2D prongs are valid for this vertex

        // 3D prongs next
        art::FindManyP<rb::Prong> fmFuzzyProng3D = FindManyPStrict<rb::Prong>(
            elastics, evt,fParams.FuzzyK3DLabel());
        unsigned int nshwlids = 0;
        if (fmFuzzyProng3D.isValid()) {
          std::vector<art::Ptr<rb::Prong>> prongs = fmFuzzyProng3D.at(vtxId);
          if (!prongs.empty()) {
            std::sort(prongs.begin(), prongs.end(),
                      sortRBProngLength);  // start sorting by length
            unsigned int pngId = 0;
            srVtx.fuzzyk.longestidx = 0;
            art::FindOneP<rb::Shower> fmShowerLID(prongs, evt,
                                                  fParams.ReclusShowerLabel());
            art::FindManyP<rb::Track> fmBPF(prongs, evt,
                                            fParams.BreakPointTrackLabel());
            art::FindManyP<rb::PID> fmCVNParticleResult(
                prongs, evt, fParams.CVNParticleLabel());
            art::FindManyP<cvn::PixelMap> fmCVNProngPixelMaps(
                prongs, evt, fParams.CVNPixelMapLabel());
            art::FindManyP<cvn::ProngTrainingData> fmCVNProngTrainingData(
                prongs, evt, fParams.CVNTrainingDataLabel());
            art::FindManyP<cvn::RegProngResult> fmRegCVNResult(prongs, evt,
                                                        fParams.RegCVNLabel());

            art::FindManyP<rb::WeightedProng> fmPngPngWt(
                prongs, evt,fParams.FuzzyKWeight3DLabel());
            float highE = 0;
            for (unsigned int iPng = 0; iPng < prongs.size(); ++iPng) {
              art::Ptr<rb::Prong> prong = prongs[iPng];
              srVtx.fuzzyk.png.push_back(SRFuzzyKProng());
              srVtx.fuzzyk.png.back().setDefault();
              FillProngVars(*prong, srVtx.fuzzyk.png.back(), *slices, allTracks,
                            sliceTracks, allTracksBirks, sliceTracksBirks,
                            sliceId);
              if (fmPngPngWt.isValid()) {
                std::vector<art::Ptr<rb::WeightedProng>> prongWt =
                    fmPngPngWt.at(iPng);
                srVtx.fuzzyk.png.back().weightedCalE =
                    prongWt[0]->CalorimetricEnergy();
              }

              if (fmCVNProngPixelMaps.isValid() && fParams.FillPixelMaps()) {
                std::vector<art::Ptr<cvn::PixelMap>> cvnProngPixelMaps =
                    fmCVNProngPixelMaps.at(iPng);
                for (unsigned int ip = 0; ip < cvnProngPixelMaps.size(); ++ip) {
                  srVtx.fuzzyk.png.back().cvnmaps.push_back(SRPixelMap());
                  SRPixelMap& pmap = srVtx.fuzzyk.png.back().cvnmaps.back();
                  FillCVNPixelMaps(*cvnProngPixelMaps[ip], pmap,
                                   fParams.UseGeVPixelMaps());
                }
              }

              cvn::ProngTrainingData cvnProngTrainingData;
              if (GetAssociatedProduct(fmCVNProngTrainingData, iPng,
                                       cvnProngTrainingData) &&
                  fParams.FillTrainingData()) {
                srVtx.fuzzyk.png.back().prongtrainingdata.push_back(
                    SRProngTrainingData());
                SRProngTrainingData& srProngTrainingData =
                    srVtx.fuzzyk.png.back().prongtrainingdata.back();
                FillCVNProngTrainingData(cvnProngTrainingData,
                                         srProngTrainingData);
              }

              if (fmCVNParticleResult.isValid()) {
                std::vector<art::Ptr<rb::PID>> cvnparts =
                    fmCVNParticleResult.at(iPng);
                FillCVNParticleResultVars(cvnparts, srVtx.fuzzyk.png.back().cvnpart);
              }
              if (fmRegCVNResult.isValid()) {
                std::vector<art::Ptr<cvn::RegProngResult>> regcvn =
                    fmRegCVNResult.at(iPng);
                if(regcvn.size()>0)
                  srVtx.fuzzyk.png.back().regcvn.prongE = regcvn[0]->Output();
                else
                  srVtx.fuzzyk.png.back().regcvn.prongE = -5.;  
              }
              if (fmShowerLID.isValid()) {
                art::Ptr<rb::Shower> slid = fmShowerLID.at(pngId);
                // look for associated LID reconstructions
                art::FindOneP<slid::ShowerLID> fmLID({slid}, evt,
                                                     fParams.LIDLabel());
                if (slid) {
                  // fill fuzzyk.png.shwlid; assume 1 shower per prong
                  FillShowerVars(*slid, srVtx.fuzzyk.png.back(), *slices,
                                 allTracks, sliceTracks, allTracksBirks,
                                 sliceTracksBirks, sliceId);
                  nshwlids++;
                  if (fmLID.isValid()) {
                    art::Ptr<slid::ShowerLID> lid = fmLID.at(0);
                    // fill fuzzyk.png.shwlid.(lid/lidE); assume 1 LID info
                    // object per shower
                    if (lid) {
                      FillSlidVars(*lid, srVtx.fuzzyk.png.back().shwlid);
                      FillLIDEnergyVars(*lid, srVtx.fuzzyk.png.back());
                      srVtx.fuzzyk.png.back().shwlid.shwE =
                          srVtx.fuzzyk.png.back().shwlid.lidE.shwE;
                      if (srVtx.fuzzyk.png.back().shwlid.shwE > highE) {
                        highE = srVtx.fuzzyk.png.back().shwlid.shwE;
                        //                        FillNueEnergyVars(slice, *slid,
                        //rec.energy.nue);
                      }
                    }
                  }  // end if LID assn valid
                }
              }  // end if showerlid is valid

              if (fmBPF.isValid()) {

                // BPF attempts (and rarely fails) to fit each prong under 3 different
                // particle assumptions:  muon, pion, proton. So for each prong, there
                // will be at most 3 tracks. Each track _should_ have one and only one
                // FitSum, OverlapE, and TrkCntObj associated with it.

                // So...   ...get all of these objects etc.
                std::vector<art::Ptr<rb::Track>> bpfTracks = fmBPF.at(pngId);
                art::FindManyP<rb::FitSum> fmBPFFitSums(
                    bpfTracks, evt, fParams.BreakPointTrackLabel());
                art::FindManyP<rb::Energy> fmBPFOverlapE(
                    bpfTracks, evt,fParams.OverlapEBPFLabel());
                art::FindManyP<cosrej::TrkCntObj> fmTrkCntBPF =
                    FindManyPStrict<cosrej::TrkCntObj>(
                        bpfTracks, evt, fParams.CosRejBPFLabel());
                art::FindManyP<trackinfo::TrackInfoObj> fmTrkInfoBPF = 
                  FindManyPStrict<trackinfo::TrackInfoObj>(
                    bpfTracks, evt, fParams.TrackInfoBPF());

                for (size_t trackId = 0; trackId < bpfTracks.size();
                     ++trackId) {
                  if (fmBPFFitSums.isValid()) {
                    std::vector<art::Ptr<rb::FitSum>> bpfFitSums =
                        fmBPFFitSums.at(trackId);
                    if (bpfFitSums.size() < 1) {
                      std::cerr << "bpfFitSums.size() < 1, this should be 1, "
                                   "since each track has a fitsum!"
                                << std::endl;
                      continue;
                    }

                    //
                    // Here we should fill the track and fitsum info
                    //

                    // For easy reference, pull out this specific prong.
                    SRFuzzyKProng& srPng = srVtx.fuzzyk.png.back();



                    // Sorry for the cut and paste, but it is only ~20 lines of actual code.
                    // Fill the track if it is for the muon, then repeat for pion and proton.

                    // If this is the muon, fill the muon track info...
                    if(abs(bpfFitSums[0]->PDG()) == 13) {

                      srPng.bpf.muon.IsValid = true;

                      // Fill the Track Vars with overlapping E if that info is
                      // present, otherwise use the default FillTrackVars
                      // function.
                      if (fmBPFOverlapE.isValid()) {
                        std::vector<art::Ptr<rb::Energy>> overlapEs =
                          fmBPFOverlapE.at(trackId);
                        if (!overlapEs.empty()) {
                          FillTrackVarsWithOverlapE(
                                                    (*bpfTracks[trackId]), (*overlapEs[0]), srPng.bpf.muon,
                                                    *slices, allTracks, sliceTracks, allTracksBirks,
                                                    sliceTracksBirks, sliceId);
                        } else {
                          FillTrackVars((*bpfTracks[trackId]), srPng.bpf.muon,
                                        *slices, allTracks, sliceTracks,
                                        allTracksBirks, sliceTracksBirks,
                                        sliceId);
                        }
                      } else {
                        FillTrackVars((*bpfTracks[trackId]), srPng.bpf.muon, *slices,
                                      allTracks, sliceTracks, allTracksBirks,
                                      sliceTracksBirks, sliceId);
                      }
                      
                      // Fill the FitSum info
                      FillTrackVarsBpfFitSum((*bpfFitSums[0]), srPng.bpf.muon);
                      if (fmTrkCntBPF.isValid()) {
                        // Get track containment info from TrkCntObj
                        std::vector<art::Ptr<cosrej::TrkCntObj>> trkcnts =
                          fmTrkCntBPF.at(trackId);
                        
                        if (!trkcnts.empty())
                          FillTrackContainmentVars(*trkcnts[0], srPng.bpf.muon);
                      }

                      // Fill the TrackInfo info (yes I said info twice...)
                      if(fmTrkInfoBPF.isValid()){
                        std::vector<art::Ptr<trackinfo::TrackInfoObj> >  trkinfos = fmTrkInfoBPF.at(trackId);
                        if(!trkinfos.empty())
                          FillTrackInfoVars(*trkinfos[0], srPng.bpf.muon);
                      }

                    } // end if(PDG == 13)

                    // If this is the pion, fill the pion track info...
                    else if(abs(bpfFitSums[0]->PDG()) == 211) {

                      srPng.bpf.pion.IsValid = true;

                      // Fill the Track Vars with overlapping E if that info is
                      // present, otherwise use the default FillTrackVars
                      // function.
                      if (fmBPFOverlapE.isValid()) {
                        std::vector<art::Ptr<rb::Energy>> overlapEs =
                          fmBPFOverlapE.at(trackId);
                        if (!overlapEs.empty()) {
                          FillTrackVarsWithOverlapE(
                                                    (*bpfTracks[trackId]), (*overlapEs[0]), srPng.bpf.pion,
                                                    *slices, allTracks, sliceTracks, allTracksBirks,
                                                    sliceTracksBirks, sliceId);
                        } else {
                          FillTrackVars((*bpfTracks[trackId]), srPng.bpf.pion,
                                        *slices, allTracks, sliceTracks,
                                        allTracksBirks, sliceTracksBirks,
                                        sliceId);
                        }
                      } else {
                        FillTrackVars((*bpfTracks[trackId]), srPng.bpf.pion, *slices,
                                      allTracks, sliceTracks, allTracksBirks,
                                      sliceTracksBirks, sliceId);
                      }
                      
                      // Fill the FitSum info
                      FillTrackVarsBpfFitSum((*bpfFitSums[0]), srPng.bpf.pion);
                      if (fmTrkCntBPF.isValid()) {
                        // Get track containment info from TrkCntObj
                        std::vector<art::Ptr<cosrej::TrkCntObj>> trkcnts =
                          fmTrkCntBPF.at(trackId);
                        
                        if (!trkcnts.empty())
                          FillTrackContainmentVars(*trkcnts[0], srPng.bpf.pion);
                      }

                      // Fill the TrackInfo info (yes I said info twice...)
                      if(fmTrkInfoBPF.isValid()){
                        std::vector<art::Ptr<trackinfo::TrackInfoObj> >  trkinfos = fmTrkInfoBPF.at(trackId);
                        if(!trkinfos.empty())
                          FillTrackInfoVars(*trkinfos[0], srPng.bpf.pion);
                      }

                    } // end if(PDG == 211)

                    // If this is the proton, fill the proton track info...
                    else if(abs(bpfFitSums[0]->PDG()) == 2212) {

                      srPng.bpf.proton.IsValid = true;

                      // Fill the Track Vars with overlapping E if that info is
                      // present, otherwise use the default FillTrackVars
                      // function.
                      if (fmBPFOverlapE.isValid()) {
                        std::vector<art::Ptr<rb::Energy>> overlapEs =
                          fmBPFOverlapE.at(trackId);
                        if (!overlapEs.empty()) {
                          FillTrackVarsWithOverlapE(
                                                    (*bpfTracks[trackId]), (*overlapEs[0]), srPng.bpf.proton,
                                                    *slices, allTracks, sliceTracks, allTracksBirks,
                                                    sliceTracksBirks, sliceId);
                        } else {
                          FillTrackVars((*bpfTracks[trackId]), srPng.bpf.proton,
                                        *slices, allTracks, sliceTracks,
                                        allTracksBirks, sliceTracksBirks,
                                        sliceId);
                        }
                      } else {
                        FillTrackVars((*bpfTracks[trackId]), srPng.bpf.proton, *slices,
                                      allTracks, sliceTracks, allTracksBirks,
                                      sliceTracksBirks, sliceId);
                      }
                      
                      // Fill the FitSum info
                      FillTrackVarsBpfFitSum((*bpfFitSums[0]), srPng.bpf.proton);
                      if (fmTrkCntBPF.isValid()) {
                        // Get track containment info from TrkCntObj
                        std::vector<art::Ptr<cosrej::TrkCntObj>> trkcnts =
                          fmTrkCntBPF.at(trackId);
                        
                        if (!trkcnts.empty())
                          FillTrackContainmentVars(*trkcnts[0], srPng.bpf.proton);
                      }

                      // Fill the TrackInfo info (yes I said info twice...)
                      if(fmTrkInfoBPF.isValid()){
                        std::vector<art::Ptr<trackinfo::TrackInfoObj> >  trkinfos = fmTrkInfoBPF.at(trackId);
                        if(!trkinfos.empty())
                          FillTrackInfoVars(*trkinfos[0], srPng.bpf.proton);
                      }

                    } // end if(PDG == 2212)

                    // otherwise, the PDG code was not one of the three track assumptions
                    else {
                      // throw a warning but continue...
                      std::cout << "\n\nWARNING:   BPF FitSum object contained non-standard PDG code "
                                << bpfFitSums[0]->PDG()
                                << ". Continuing but not filling info for this track.\n\n";
                    } // end else

                  } // fmBPFFitSums.isValid()
                  else std::cerr << "fmBPFFitSums is not valid" << std::endl;

                }  // loop trackId

              }    // fmBPF.isValid()

              pngId++;

            }  // end loop over 3D prongs
          }    // there are non-zero 3D prongs
          rec.vtx.elastic.fuzzyk.nshwlid = nshwlids;
          rec.vtx.elastic.fuzzyk.fillSizes();

          // energy from prong orphaned hits
          art::FindManyP<rb::Prong> fmProngOrph(
              elastics, evt,fParams.FuzzyKOrphLabel());
          if (fmProngOrph.isValid()) {
            std::vector<art::Ptr<rb::Prong>> prongsOrph = fmProngOrph.at(vtxId);
            if (prongsOrph.size() > 0) {
              rec.vtx.elastic.fuzzyk.orphCalE =
                  prongsOrph[0]->CalorimetricEnergy();
            }
          }

          std::sort(srVtx.fuzzyk.png.begin(), srVtx.fuzzyk.png.end(), sortProngCalE);  // Use prong calE to sort

          float longest = 0;
          for (unsigned int pngIdx = 0; pngIdx < srVtx.fuzzyk.png.size();
               pngIdx++) {
            if (srVtx.fuzzyk.png[pngIdx].len > longest) {
              longest = srVtx.fuzzyk.png[pngIdx].len;
              srVtx.fuzzyk.longestidx = pngIdx;  // note the longest prong index
            }
          }

          if (nshwlids >
              0) {  // fill nue energy for primary shower if there were showers
            rec.energy.nue.lid.E =
                rec.vtx.elastic.fuzzyk.png[0].shwlid.lidE.E;
            rec.energy.nue.lid.depE =
                rec.vtx.elastic.fuzzyk.png[0].shwlid.lidE.depE;
            rec.energy.nue.lid.shwE =
                rec.vtx.elastic.fuzzyk.png[0].shwlid.lidE.shwE;
            rec.energy.nue.lid.hadE =
                rec.vtx.elastic.fuzzyk.png[0].shwlid.lidE.hadE;
          }

        }  // end if 3D prongs are valid for this vertex
      }    // end loop over vertices

    }  // end ElasticArms

    FillNueEnergyVars(slice, rec.vtx, rec.energy.nue);

    // HoughVertex

    if (fmHoughVertex.isValid()) {
      std::vector<art::Ptr<rb::Vertex>> houghs = fmHoughVertex.at(sliceId);

      for (unsigned int vtxId = 0; vtxId < houghs.size(); ++vtxId) {
        rec.vtx.hough.push_back(SRHoughVertex());
        SRHoughVertex& srVtx = rec.vtx.hough.back();

        srVtx.time = houghs[vtxId]->GetT();
        srVtx.vtx = houghs[vtxId]->GetXYZ();

        /*          if(fmFuzzyProng3D.isValid()){
          std::vector< art::Ptr<rb::Prong> > prongs = fmFuzzyProng3D.at(vtxId);
          // Sort prongs in descending length order
          std::sort(prongs.begin(), prongs.end(),sortRBProngLength);
          for(art::Ptr<rb::Prong> prong: prongs){
            if(prong->NXCell() == 0 || prong->NYCell() == 0){
              srVtx.fuzzyk.png2d.push_back(SRProng());
              FillProngVars(*prong, srVtx.fuzzyk.png2d.back(), slicess,
        sliceId);
            }// end 2d prong case
            else{
              srVtx.fuzzyk.png.push_back(SRProng());
              FillProngVars(*prong, srVtx.fuzzyk.png.back(), *slices, sliceId);
            }// end 3d prong case
          }// end loop over prongs
        }// end if prongs are valid for this vertex
        // Check if there is a valid instance label for 2D prongs
        if(fmFuzzyProng2D.isValid() && fFuzzyK2dAssnLabel != ""){
          std::vector< art::Ptr<rb::Prong> > prongs = fmFuzzyProng2D.at(vtxId);
          // Sort prongs in descending length order
          std::sort(prongs.begin(), prongs.end(),sortRBProngLength);
          for(art::Ptr<rb::Prong> prong: prongs){
            if(prong->NXCell() == 0 || prong->NYCell() == 0){
              srVtx.fuzzyk.png2d.push_back(SRProng());
              FillProngVars(*prong, srVtx.fuzzyk.png2d.back(), *slices,
        sliceId);
            }// end 2d prong case
            else{
              srVtx.fuzzyk.png.push_back(SRProng());
              FillProngVars(*prong, srVtx.fuzzyk.png.back(), *slices, sliceId);
            }// end 3d prong case
          }// end loop over prongs
        }// end if prongs are valid for this vertex
        */
      }  // end for vtxId
      // Sort the vector
      //        std::sort(rec.vtx.hough.begin(), rec.vtx.hough.end(),
      //        sortHoughProng);
    }

    // Get all the VertexDT vertices associated with the slice
    if (fmVertexDT.isValid()) {
      std::vector<art::Ptr<rb::Vertex>> vtxs = fmVertexDT.at(sliceId);
      for (unsigned int vtxId = 0; vtxId < vtxs.size(); ++vtxId) {
        rec.vtx.vdt.push_back(SRVertexDT());
        SRVertexDT& srVtx = rec.vtx.vdt.back();
        srVtx.time = vtxs[vtxId]->GetT();
        srVtx.vtx = vtxs[vtxId]->GetXYZ();
      }
    }
    rec.vtx.fillSizes();

    // Get all MEFs aassociated with slice
    if (fmSlcME.isValid()) {
      std::vector<art::Ptr<me::SlcME>> michels = fmSlcME.at(sliceId);
      for (art::Ptr<me::SlcME> michel : michels) {
        AddSlcMEToVec(*michel, &rec.me.slc, *slices, allTracks, sliceTracks,
                      allTracksBirks, sliceTracksBirks, sliceId);
      }
    }
    rec.me.fillSizes();

    // Get slice level PID from LIDBuilder
    if (foLID.isValid()) {
      cet::maybe_ref<art::Ptr<slid::EventLID> const> rlid(foLID.at(sliceId));
      art::Ptr<slid::EventLID> elid = rlid.ref();
      // If the pointer is valid
      if (elid) {
        rec.sel.lid.ann = elid->Value();
        rec.sel.lid.anne = elid->ValueWithE();
        rec.sel.lid.annepi0 = elid->ValueEPi0();
        rec.sel.lid.annecos = elid->ValueECos();
      } else {
        rec.sel.lid.setDefault();
      }
    }

    if (foSlcLID.isValid()) {
      cet::maybe_ref<art::Ptr<slid::SliceLID> const> rslclid(
          foSlcLID.at(sliceId));
      art::Ptr<slid::SliceLID> slclid = rslclid.ref();
      if (slclid) {
        rec.sel.lid.rnn = slclid->Value();
      } else {
        rec.sel.lid.setDefault();
      }
    }

    /*
    if(fmShwLID.isValid()){
      std::vector<const  slid::ShowerLID* > shwlids = fmShwLID.at(sliceId);
        // Sort here since the vector may not be
      sort(shwlids.begin(),shwlids.end(), slid::CompareByEnergy);
      art::FindOneP<rb::Shower> fos(shwlids, evt, fLIDLabel());
      for(unsigned int shwId = 0; shwId < shwlids.size(); ++shwId){
        cet::maybe_ref<art::Ptr<rb::Shower> const> rshw(fos.at(shwId));
        art::Ptr<rb::Shower> shw = rshw.ref();
        if (shw){
          rec.shw.shwlid.push_back(SRShower());
          FillShowerVars(*shw, rec.shw.shwlid.back(), *slices, sliceId);
        }

        rec.sel.elecid.shwlid.push_back(SRSLid());
        FillSlidVars(*shwlids[shwId], rec.sel.elecid.shwlid.back().shwlid);

        // Fill energy information for leading energy shower
              // which also serves as a slice pid
        if (shwId == 0){
          FillLIDEnergyVars(*shwlids[shwId], rec.energy.elecid);
        }

        // Fill energy vector with the other shower energies
        rec.energy.shwlid.push_back(SRSLidEnergy());
        FillLIDEnergyVars(*shwlids[shwId], rec.energy.shwlid.back());

        if (shw && shwId == 0){
          FillNueEnergyVars(slice, *shw, rec.energy.nue);
        }
      } // end for loop on showerlids
    } // end if showerlids are valid

    if(fmShwPID.isValid()){
      std::vector<const  slid::ShowerPID* > shwpids = fmShwPID.at(sliceId);
      // Sort here since the vector may not be
      sort(shwpids.begin(),shwpids.end(), slid::CompareByEnergy);
      art::FindOneP<rb::Shower> fos(shwpids, evt, fSPIDLabel());
      for(unsigned int shwId = 0; shwId < shwpids.size(); ++shwId){
        rec.sel.elecid.shwpid.push_back(SRSPid());
        FillSpidVars(*shwpids[shwId], rec.sel.elecid.shwpid.back());
      } // end for loop on showerpids
    } // end if showerpids are valid

  } // end if event id is valid
  else{
    rec.sel.lid        .setDefault();
    rec.energy.nue.lid .setDefault();
  }
} // end if slice lids are valid
else{
  rec.sel.lid    .setDefault();
  rec.energy.lid .setDefault();
}

rec.shw.fillSizes();
rec.energy.fillSizes();
rec.sel.elecid.fillSizes();
    */
    lem::PIDDetails lem;
    bool isLemmed = GetAssociatedProduct(fmLem, sliceId, lem);
    std::vector<std::string> vec = {fParams.LEMNuePresel()};

    if (isLemmed && !(fParams.ApplyLEMNuePresel() &&
                      rb::IsFiltered(evt, slices, sliceId, vec))) {
      FillLEMVars(lem, rec.sel.lem);
    } else
      rec.sel.lem.setDefault();

    rvp::RVP rvp;
    if (GetAssociatedProduct(fmRvp, sliceId, rvp))
      FillRVPVars(rvp, rec.sel.rvp);
    else
      rec.sel.rvp.setDefault();

    xnue::Xnue xnue;
    if (GetAssociatedProduct(fmXnue, sliceId, xnue)) {
      FillXnueVars(xnue, rec.sel.xnuepid);
    } else
      rec.sel.xnuepid.setDefault();

    numue::NumuE numuE;
    if (GetAssociatedProduct(fmNumuEnergy, sliceId, numuE)) {
      FillNumuEnergyVars(numuE, rec.energy.numu);
    } else {
      rec.energy.numu.setDefault();
    }


    cvn::RegNuResult evtregcvn;
    if (GetAssociatedProduct(fmEvtRegCVN, sliceId, evtregcvn)) {
      rec.energy.nue.regcvnEvtE  = evtregcvn.Output();
    } else {
      rec.energy.nue.regcvnEvtE = -5.;
    }



    LSTME::LSTMEnergy numulstmE;
    if (GetAssociatedProduct(fmNumuLSTMEnergy, sliceId, numulstmE)) {
      FillNumuLSTMEnergyVars(numulstmE, rec.energy.numu);
    } else {
      rec.energy.numu.setLSTMDefault();
    }
    
    cvn::RegHadronResult hadregcvn;
    if (GetAssociatedProduct(fmHadRegCVN, sliceId, hadregcvn)) {
      rec.energy.numu.regcvnhadE  = hadregcvn.Output();
    } else {
      rec.energy.numu.regcvnhadE = -5.;
    }

    SliceLID::Prediction sliceLID;
    if (GetAssociatedProduct(fmSliceLID, sliceId, sliceLID)) {
      FillSliceLID(sliceLID, rec.sel.slicelid);
    } else {
      rec.sel.slicelid.setDefault();
    }

    cosrej::CosRejObj cosrej;
    if (GetAssociatedProduct(fmCosRej, sliceId, cosrej)) {
      FillCosRejVars(cosrej, rec.sel.cosrej, rec.sel.contain);

      float kDistAllTop = rec.sel.nuecosrej.distallpngtop;
      if (std::isnan(kDistAllTop)) kDistAllTop = -1000.0;

      float kDistAllBottom = rec.sel.nuecosrej.distallpngbottom;
      if (std::isnan(kDistAllBottom)) kDistAllBottom = -1000.0;

      float kDistAllEast = rec.sel.nuecosrej.distallpngeast;
      if (std::isnan(kDistAllEast)) kDistAllEast = -1000.0;

      float kDistAllWest = rec.sel.nuecosrej.distallpngwest;
      if (std::isnan(kDistAllWest)) kDistAllWest = -1000.0;

      float kDistAllBack = rec.sel.nuecosrej.distallpngback;
      if (std::isnan(kDistAllBack)) kDistAllBack = -1000.0;

      float kDistAllFront = rec.sel.nuecosrej.distallpngfront;
      if (std::isnan(kDistAllFront)) kDistAllFront = -1000.0;

      if (fDet ==
          kFARDET) {  // adopted from CAFAna/Cuts/NumuCuts2017.h on 8.16.17
        std::pair<int, int> planes = calcFirstLastLivePlane(
            rec.slc.firstplane, std::bitset<14>(rec.hdr.dibmask));
        int planestofront = rec.slc.firstplane - planes.first;
        int planestoback = planes.second - rec.slc.lastplane;
        bool kcut1 =
            (rec.sel.contain.kalfwdcell > 6 && rec.sel.contain.kalbakcell > 6 &&
             rec.sel.contain.cosfwdcell > 0 && rec.sel.contain.cosbakcell > 7 &&
             planestofront > 1 && planestoback > 1);
        bool kcut2 = kDistAllTop > 60 && kDistAllBottom > 12 &&
                     kDistAllEast > 16 && kDistAllWest > 12 &&
                     kDistAllFront > 18 && kDistAllBack > 18;
        rec.sel.contain.numucontain = (kcut1 && kcut2);

        if (cosrej.CosFwdCell() > 0 && cosrej.CosBakCell() > 0 &&
            cosrej.KalFwdCell() > 10 && cosrej.KalBakCell() > 10 &&
            rec.slc.ncellsfromedge > 1 && rec.sel.contain.nplanestofront > 1 &&
            rec.sel.contain.nplanestoback > 1) {
          rec.sel.contain.numucontainSA = true;
        } else
          rec.sel.contain.numucontainSA = false;
      }

      else if (fDet == kNEARDET) {  // adopted from CAFAna/Cuts/NumuCuts2017.h
                                    // on 8.16.17
        if (rec.vtx.elastic.IsValid < 1)
          rec.sel.contain.numucontain = false;
        else {
          for (unsigned int i = 0; i < rec.vtx.elastic.fuzzyk.nshwlid; ++i) {
            TVector3 start = rec.vtx.elastic.fuzzyk.png[i].shwlid.start;
            TVector3 stop = rec.vtx.elastic.fuzzyk.png[i].shwlid.stop;
            if (std::min(start.X(), stop.X()) < -180.0)
              rec.sel.contain.numucontain = false;
            if (std::max(start.X(), stop.X()) > 180.0)
              rec.sel.contain.numucontain = false;
            if (std::min(start.Y(), stop.Y()) < -180.0)
              rec.sel.contain.numucontain = false;
            if (std::max(start.Y(), stop.Y()) > 180.0)
              rec.sel.contain.numucontain = false;
            if (std::min(start.Z(), stop.Z()) < 20.0)
              rec.sel.contain.numucontain = false;
            if (std::max(start.Z(), stop.Z()) > 1525.0)
              rec.sel.contain.numucontain = false;
          }
          if (rec.trk.kalman.ntracks < 1) rec.sel.contain.numucontain = false;
          for (unsigned int i = 0; i < rec.trk.kalman.ntracks; ++i) {
            if (i == rec.trk.kalman.idxremid)
              continue;
            else if (rec.trk.kalman.tracks[i].start.Z() > 1275 ||
                     rec.trk.kalman.tracks[i].stop.Z() > 1275)
              rec.sel.contain.numucontain = false;
          }
        }

        rec.sel.contain.numucontain =
            (rec.trk.kalman.ntracks > rec.trk.kalman.idxremid &&
             rec.slc.firstplane > 1 && rec.slc.lastplane < 212 &&
             rec.trk.kalman.tracks[0].start.Z() < 1100 &&
             (rec.trk.kalman.tracks[0].stop.Z() < 1275 ||
              rec.sel.contain.kalyposattrans < 55) &&
             rec.sel.contain.kalfwdcellnd > 5 &&
             rec.sel.contain.kalbakcellnd > 10);

        if (cosrej.KalFwdCellND() > 4 && cosrej.KalBakCellND() > 8 &&
            rec.slc.ncellsfromedge > 1 && rec.slc.firstplane > 1 &&
            rec.slc.lastplane < 212 &&
            rec.trk.kalman.tracks[0].start.Z() < 1150 &&
            (rec.trk.kalman.tracks[0].stop.Z() < 1275 ||
             rec.sel.contain.kalyposattrans < 55) &&
            (rec.energy.numu.ndhadcaltranE + rec.energy.numu.ndhadcalcatE) <
                0.03) {
          rec.sel.contain.numucontainSA = true;
        } else
          rec.sel.contain.numucontainSA = false;
      }
    }

    cosrej::NueCosRej nuecosrej;
    if (GetAssociatedProduct(fmNueCosRej, sliceId, nuecosrej)) {
      FillNueCosRejVars(nuecosrej, rec.sel.nuecosrej);
      /*
      // nue containment boolean - nue group please feel free to update.
            // If variables move from sandbox this needs to be changed too
      if(fDet == kFARDET){
        if(nuecosrej.ProngMaxX() < 750 && nuecosrej.ProngMinX() > -750 &&
           nuecosrej.ProngMaxY() < 750 && nuecosrej.ProngMinY() > -750 &&
           std::min(nuecosrej.StartDistToFront(),
                    nuecosrej.StopDistToFront()) > 50 &&
           std::min(nuecosrej.StartDistToBack(),
                    nuecosrej.StopDistToBack()) > 50){
          rec.sel.contain.nuecontain = true;
        }
      }

      if(fDet == kNEARDET){
        if(nuecosrej.ProngMaxX() < 185 && nuecosrej.ProngMinX() > -185 &&
           nuecosrej.ProngMaxY() < 185 && nuecosrej.ProngMinY() > -185 &&
           std::min(nuecosrej.StartDistToFront(),
                    nuecosrej.StopDistToFront()) > 50 &&
           std::min(nuecosrej.StartDistToBack(),
                    nuecosrej.StopDistToBack()) > 50){
          rec.sel.contain.nuecontain = true;
        }
      }
      */
    } else {
      rec.sel.nuecosrej.setDefault();
    }

    ncid::NCCosRej nccosrej;  /////////////  NC Cosmic Rejection
    if (GetAssociatedProduct(fmNCCosRej, sliceId, nccosrej)) {
      FillNCCosRejVars(nccosrej, rec.sel.nccosrej);
    } else {
      rec.sel.nccosrej.setDefault();
    }

    ncpi0::NCPi0BkgRej ncpi0bkgrej;
    if (GetAssociatedProduct(fmNCPi0BkgRej, sliceId, ncpi0bkgrej)) {
      FillNCPi0BkgRejVars(ncpi0bkgrej, rec.sel.ncpi0bkgrej);
    } else {
      rec.sel.ncpi0bkgrej.setDefault();
    }

    presel::PreselObj nuepre;
    if (GetAssociatedProduct(fmNuePresel, sliceId, nuepre)) {
      FillNuePreselVars(nuepre, rec.sel.nuepre);
    } else {
      rec.sel.nuepre.setDefault();
    }

    presel::PreselObj rockpre;
    if (GetAssociatedProduct(fmRockPresel, sliceId, rockpre)) {
      FillRockPreselVars(rockpre, rec.sel.rockpre);
    } else {
      rec.sel.rockpre.setDefault();
    }

    presel::Veto veto, nueveto;
    if (GetAssociatedProduct(fmVeto, sliceId, veto) ||
        GetAssociatedProduct(fmNueVeto, sliceId, nueveto)) {
      FillVetoVars(veto, nueveto, rec.sel.veto);
    } else {
      rec.sel.veto.setDefault();
    }
    
    // --- Horrible hack to appease CAFAna.
    rec.sel.cvn.setDefault();

    cvn::Result nuoneResult;
    if (GetAssociatedProduct(fmNuonE, sliceId, nuoneResult)) {
      FillNuonEResultVars(nuoneResult, rec.sel.nuone);
    } else {
      rec.sel.nuone.setDefault();
    }

    // --- Loose Presel PtP Cut
    cvn::Result cvnResult_LoosePreselPtP;
    if (GetAssociatedProduct(fmCVN_LoosePreselPtP, sliceId, cvnResult_LoosePreselPtP)) {
      FillCVNResultVars(cvnResult_LoosePreselPtP, rec.sel.cvnloosepreselptp);
    } else {
      rec.sel.cvnloosepreselptp.setDefault();
    }
    // --- Old Presel
    cvn::Result cvnResult_OldPresel;
    if (GetAssociatedProduct(fmCVN_OldPresel, sliceId, cvnResult_OldPresel)) {
      FillCVNResultVars(cvnResult_OldPresel, rec.sel.cvnoldpresel);
    } else {
      rec.sel.cvnoldpresel.setDefault();
    }
    // --- No Cosmics used in training
    cvn::Result cvnResult_NoCosmics;
    if (GetAssociatedProduct(fmCVN_NoCosmics, sliceId, cvnResult_NoCosmics)) {
      FillCVNResultVars(cvnResult_NoCosmics, rec.sel.cvnnocosmics);
    } else {
      rec.sel.cvnnocosmics.setDefault();
    }
    // --- Wrong Sign
    cvn::Result cvnResult_ws;
    if (GetAssociatedProduct(fmWS, sliceId, cvnResult_ws)) {
      rec.sel.wsid = cvnResult_ws.fOutput[0];
    }

    cvn::Features cvnFeatures;
    if (GetAssociatedProduct(fmCVNFeatures, sliceId, cvnFeatures) &&
        fParams.FullTrainingInfo()) {
      FillCVNFeaturesVars(cvnFeatures, rec.training.cvnfeatures,
                          fParams.NumCVNPrincipalComponents());
    }
    
    std::cout<<"------------> wwj: FillCVNPixelMaps"<<std::endl;
    std::cout<<fmCVNPixelMaps.size()<<std::endl;
    if (fmCVNPixelMaps.isValid() && fParams.FillPixelMaps()) {
      std::vector<art::Ptr<cvn::PixelMap>> cvnPixelMaps =
          fmCVNPixelMaps.at(sliceId);
      std::cout<<cvnPixelMaps.size()<<std::endl;
      for (unsigned int ip = 0; ip < cvnPixelMaps.size(); ++ip) {
        rec.training.cvnmaps.push_back(SRPixelObjMap());
        SRPixelObjMap& pmap = rec.training.cvnmaps.back();
        FillCVNPixelMaps(*cvnPixelMaps[ip], pmap, fParams.UseGeVPixelMaps());
      }
    }

    std::cout<<"------------> wwj: FillCVNTrainingData"<<std::endl;
    cvn::TrainingData cvnTrainingData;
    if (GetAssociatedProduct(fmCVNTrainingData, sliceId, cvnTrainingData) &&
        fParams.FillTrainingData()) {
      rec.training.trainingdata.push_back(SRTrainingData());
      SRTrainingData& srTrainingData = rec.training.trainingdata.back();
      FillCVNTrainingData(cvnTrainingData, srTrainingData);
    }

    // also store every hit in the slice.  Doing this in a fixed size array the
    // size of the detector since the zeros compress better
    if (fParams.FillPixelMaps()) {
      rec.training.slicemaps.push_back(SRSliceMap());
      SRSliceMap& smap = rec.training.slicemaps.back();
      float corr = 15123.7;  // Rough PE to GeV? Rev 17083
      float peCorrChunk;
      // This value is 2.0 / (PE/GeV)
      if (fParams.UseGeVPixelMaps())
        peCorrChunk = 1.322e-4 / 255.0;
      else
        peCorrChunk = 2.0 / 255.0;
      for (unsigned int ih = 0; ih < slice.NCell(); ih++) {
        rb::RecoHit rh = slice.RecoHit(ih);
        art::Ptr<rb::CellHit> ch = slice.Cell(ih);
        // convert planes to 0 to 447 in each view
        unsigned int pl = (ch->Plane() / 2);
        unsigned int chan = ch->Cell() + 384 * (pl + 448 * ch->View());
        // convert cell energy to an unsigned char.  Scales same as in CVN
        float hitPE = ch->PE();
        float gev = (rh.IsCalibrated()) ? rh.GeV() : hitPE / corr;
        float pe = (rh.IsCalibrated()) ? rh.PECorr() : hitPE;
        float tv = (fParams.UseGeVPixelMaps()) ? gev : pe;
        tv /= 2100.0;  // scale factor from CVN Pixel maps
        float truncateCorr = ceil(tv / peCorrChunk);
        unsigned char value;
        if (!fParams.UseEnergy()) {
          if (truncateCorr > 0)
            truncateCorr = 1;
          else
            truncateCorr = 0;
        }
        if (truncateCorr > 255)
          value = 255;
        else
          value = (unsigned char)truncateCorr;
        smap.slicemap[chan] = value;
      }
    }

    murem::MRCCParent parent;
    if (GetAssociatedProduct(fmMrccPar, sliceId, parent))
      FillMRCCParentInfo(parent, rec.parent.mrccpar, evt);
    else
      rec.parent.mrccpar.setDefault();

    //#######################################################
    // Fill truth information
    //#######################################################

    // Set mc branch values to default
    rec.mc.setDefault();
    if (bt->HaveTruthInfo()) {
      // Get all truth track hits associated with slice

      int currentNuIdx = 0;
      unsigned int nNuGlobal =
          0;  // count of neutrino interactions in this slice

      for (unsigned int iEffPur = 0; iEffPur < sliceTruthTable[sliceId].size();
           ++iEffPur) {
        // continue if truth left no hits in slice
        if (!(sliceTruthTable[sliceId][iEffPur].nSliceHits > 0)) continue;

        if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
            faveIdsEnergy[sliceId]) {
          cheat::NeutrinoEffPur* faveNu = &sliceTruthTable[sliceId][iEffPur];

          // If NOT Neutrinos
          if (!faveNu->neutrinoInt->NeutrinoSet()) {
            // Add faveNu to allcosmics
            if (fParams.FullTruthInfo())
              AddCosmicTruthToVec(faveNu, allTracks, sliceTracks,
                                  allTracksBirks, sliceTracksBirks,
                                  &rec.mc.allcosmics);

            // Add faveNu to cosmic
            AddCosmicTruthToVec(faveNu, allTracks, sliceTracks, allTracksBirks,
                                sliceTracksBirks, &rec.mc.cosmic);
          }
          // Otherwise we are a neutrino
          else {
            // Add faveNu to allnus
            if (fParams.FullTruthInfo())
              AddMCTruthToVec(faveNu, allTracks, sliceTracks, allTracksBirks,
                              sliceTracksBirks, &rec.mc.allnus, evt, spill);

            // Add faveNu to nu
            AddMCTruthToVec(faveNu, allTracks, sliceTracks, allTracksBirks,
                            sliceTracksBirks, &rec.mc.nu, evt, spill);
            rec.mc.faveidxenergy = currentNuIdx;

            art::PtrVector<simb::MCTruth> pv;
            pv.push_back(faveNu->neutrinoInt);

            if (rec.mc.nu.size() > 0) {
              // if a neutrino is matched to the slice, fill the mc reweighting
              // object.

              art::FindManyP<rwgt::GENIEReweightTable> fmRW(
                  pv, evt, fParams.ReweightLabel());
              rwgt::GENIEReweightTable rw;
              if (GetAssociatedProduct(fmRW, 0, rw))
                rec.mc.nu[0].rwgt.genie = GenieReweightTable(rw);

              art::FindManyP<fxwgt::FluxWeights> fmFLXRW(
                  pv, evt, fParams.FluxReweightLabel());
              fxwgt::FluxWeights flxrw;
              if (GetAssociatedProduct(fmFLXRW, 0, flxrw))
                rec.mc.nu[0].rwgt.ppfx = FluxReweights(flxrw);
              else
                rec.mc.nu[0].rwgt.ppfx.setDefault();
            }

            if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
                faveIdsEff[sliceId])
              rec.mc.faveidxeff = currentNuIdx;
            if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
                faveIdsPur[sliceId])
              rec.mc.faveidxpur = currentNuIdx;
            if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
                faveIdsEffPur[sliceId])
              rec.mc.faveidxeffpur = currentNuIdx;
            if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
                faveIdsEffThenPur[sliceId])
              rec.mc.faveidxeffthenpur = currentNuIdx;

            currentNuIdx++;
          }
        }
      }

      for (unsigned int iEffPur = 0; iEffPur < sliceTruthTable[sliceId].size();
           ++iEffPur) {
        if (sliceTruthTable[sliceId][iEffPur].neutrinoInt &&
            sliceTruthTable[sliceId][iEffPur].neutrinoInt->NeutrinoSet())
          nNuGlobal++;
        // continue if truth left no hits in slice
        if (!(sliceTruthTable[sliceId][iEffPur].nSliceHits > 0)) continue;

        if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex !=
            faveIdsEnergy[sliceId]) {
          // Cosmics
          if (!sliceTruthTable[sliceId][iEffPur].neutrinoInt->NeutrinoSet()) {
            if (fParams.FullTruthInfo())
              AddCosmicTruthToVec(&sliceTruthTable[sliceId][iEffPur], allTracks,
                                  sliceTracks, allTracksBirks, sliceTracksBirks,
                                  &rec.mc.allcosmics);
          }
          // Neutrino
          else {
            if (fParams.FullTruthInfo())
              AddMCTruthToVec(&sliceTruthTable[sliceId][iEffPur], allTracks,
                              sliceTracks, allTracksBirks, sliceTracksBirks,
                              &rec.mc.allnus, evt, spill);

            if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
                faveIdsEff[sliceId])
              rec.mc.faveidxeff = currentNuIdx;
            if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
                faveIdsPur[sliceId])
              rec.mc.faveidxpur = currentNuIdx;
            if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
                faveIdsEffPur[sliceId])
              rec.mc.faveidxeffpur = currentNuIdx;
            if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
                faveIdsEffThenPur[sliceId])
              rec.mc.faveidxeffthenpur = currentNuIdx;

            currentNuIdx++;
          }
        }
      }  // end loop to fill all* vectors

      // fill size fields for mc branch
      rec.mc.nnu = rec.mc.nu.size();
      rec.mc.ncosmic = rec.mc.cosmic.size();
      if (fParams.FullTruthInfo()) {
        rec.mc.nallnus = rec.mc.allnus.size();
        rec.mc.nallcosmics = rec.mc.allcosmics.size();
      }
      rec.mc.global.nnu = nNuGlobal;
    }  // end MC block

    if (fParams.EnableBlindness()) BlindThisRecord(&rec);

    fRecTree->Fill();
    srcol->push_back(rec);
    util::CreateAssn(*this, evt, *srcol, art::Ptr<rb::Cluster>(slices, sliceId),
                     *srAssn);
  }  // end loop over slices

  *spillcol = spill;

  evt.put(std::move(srcol));
  evt.put(std::move(spillcol));
  if (fParams.FillNuTree()) evt.put(std::move(spilltruthcol));
  evt.put(std::move(srAssn));
}

void CAFMaker::endSubRun(art::SubRun& sr) {
  art::Handle<sumdata::POTSum> pot;

  // If this is MC get the spill info from the generator
  if (!fIsRealData) {
    sr.getByLabel(fParams.GeneratorLabel(), pot);
  }
  // If we have real data try to get the spill info from the NuMI beam
  else {
    sr.getByLabel(fParams.NuMIBeamLabel(), pot);
  }

  if (!pot.failedToGet()) {
    // NOTE: The POT that comes out of IFDBSpillInfo has an implied
    // factor of 1E12 that is being left off.
    // For now the data POT is being multiplied by this factor so that
    // the number is consistent with what we get in MC.
    // TODO: Is this what we want to do in the future?
    float subrunpot = pot->totgoodpot;
    if (fIsRealData) subrunpot *= 1e12;
    fTotalPOT += subrunpot;

    mf::LogInfo("CAFMaker")
        << "POT in subrun " << sr.subRun() << " = " << pot->totgoodpot;
  }

  art::Handle<sumdata::CosmicExposure> expo;
  sr.getByLabel(fParams.ExposureLabel(), expo);

  if (!expo.failedToGet()) {
    fTotalLivetime += expo->totlivetime;
    mf::LogInfo("CAFMaker")
        << "Exposure in subrun " << sr.subRun() << " = " << expo->totlivetime;
  }
}

//......................................................................
void CAFMaker::endJob() {
  bool EmptyFile = false;
  if (fTotalEvents == 0) {
    // Among other things this makes us unable to determine the file_type
    // metdata field.
    if (fParams.AbortOnEmpty()) {
      mf::LogError("CAFMaker")  << "No events processed in this file. Aborting rather than "
                                << "produce an empty CAF."
                                << std::endl;
      abort();
    }
    else {
      mf::LogWarning("CAFMaker") << "Making an empty CAF file." << std::endl;
      EmptyFile = true;
    }
  }

  // Make sure the recTree is safely in the file before starting on the POT
  // etc histograms. This way, absence of the histograms should be a clear
  // sign something went wrong with the tree.
  fFile->Write();

  fFile->cd();
  hPOT->Fill(.5, fTotalPOT);
  hEvents->Fill(.5, fTotalEvents);
  hLivetime->Fill(.5, fTotalLivetime);

  hPOT->Write();
  hEvents->Write();
  hLivetime->Write();
  fFile->Write();

  if (fParams.IsSinglesOverlay()) {
    hTotalTrueNonswapNue->Fill(0.5, fTotalTrueNonswapNue - fTotalTrueSingleNue);
    hTotalTrueSingleNue->Fill(0.5, fTotalTrueSingleNue);
    hSinglePOT->Fill(0.5, fTotalSinglePOT);

    hSinglePOT->Write();
    hTotalTrueNonswapNue->Write();
    hTotalTrueSingleNue->Write();
  }

  std::map<std::string, std::string> metadata =
      meta::MetadataManager::getInstance().GetMetadata();
  metadata["data_tier"] = fParams.DataTier();
  if (EmptyFile)
    metadata["file_type"] = "empty";
  else if (fIsRealData)
    metadata["file_type"] = "importedDetector";
  else
    metadata["file_type"] = "importedSimulated";

  // Gymnastics to mark "restrictedbox" in nova.special for FD numi data
  if (metadata.count("online.stream")) {
    // Then this is a real data file and the stream exists in the metadata
    if (!metadata["online.stream"].compare("0") &&
        fDetID == novadaq::cnv::kFARDET) {
      // Then this is a numi trigger file from the far detector
      if (!fParams.EnableBlindness()) {
        // Then we are not blinding it, need to mark it as such in the metadata
        std::string tag("restrictedbox");

        if (metadata.count("nova.special") &&
            std::string("").compare(metadata["nova.special"]) &&
            std::string("none").compare(
                metadata["nova.special"])) {  // Then the metadata is not the
                                              // same as none, we need to append
                                              // to it
          // Use a stringstream to jam the strings together and append the tag

          std::stringstream ss;
          ss.str();
          ss << metadata["nova.special"] << "-" << tag;
          metadata["nova.special"] = ss.str();
        } else
          metadata["nova.special"] = tag;
      }
    }
  }

  AddMetadataToFile(metadata);
}

//......................................................................
double CAFMaker::SimpleOscProb(const simb::MCFlux& flux,
                               const simb::MCNeutrino& nu) const {
  if (nu.CCNC() == 1 && nu.Nu().PdgCode() == flux.fntype) {
    // NC
    return 1;
  }
  if (nu.CCNC() == 1 && nu.Nu().PdgCode() != flux.fntype) {
    // Swapped NC. Doesn't exist. If people want to play clever
    // reweighting games to reclaim these, they can figure that out,
    // but this allows the naive implementation to work.
    return 0;
  }

  osc::OscCalcDumb calc;
  return calc.P(flux.fntype, nu.Nu().PdgCode(), nu.Nu().E());
}

//......................................................................
void CAFMaker::FillSpillVars(
    art::Handle<sumdata::SpillData> spillPot,
    art::Handle<sumdata::EventQuality> spillQual,
    art::Handle<dq::DAQEventSummary> daq,
    art::Handle<std::vector<rawdata::RawTrigger>> trigs, SRSpill& spill,
    art::Event& evt) const {
  // Add in run, subrun, evt information so that the spill tree can be
  // mapped to the record tree.
  spill.run = evt.run();
  spill.subrun = evt.subRun();
  spill.evt = evt.id().event();

  // If we fail to identify a spill, we want it set to bad

  // sumdata::SpillData variables
  spill.det = fDet;
  spill.ismc = !fIsRealData;
  spill.isgoodspill = false;
  spill.spilltimesec = 0;
  spill.spilltimensec = 0;
  spill.gpsspilltimesec = 0;
  spill.gpsspilltimensec = 0;
  spill.deltaspilltimensec = 0;
  spill.spillpot = 0.0;
  spill.livetime = 0.0;
  spill.hornI = 0.0;
  spill.posx = 0.0;
  spill.posy = 0.0;
  spill.widthx = 0.0;
  spill.widthy = 0.0;

  // Calculate diblock masking.  Also in header, but needed here for
  // mass accounting.
  unsigned int maskstatus = 0;
  unsigned int dibmask = fMask;

  if (rh->GoodDiBlockMask() == ((1 << 16) - 1))
    maskstatus = 0;
  else if (rh->GoodDiBlockMask() == ((1 << 17) - 1))
    maskstatus = 2;
  else
    maskstatus = 1;

  unsigned int dibfirst = 0;
  unsigned int diblast = 0;
  if (dibmask) {
    int iD;
    iD = 0;
    while (!((dibmask >> iD) & 1)) iD++;
    dibfirst = iD + 1;
    iD = 0;
    while (dibmask >> iD) iD++;
    diblast = iD;
  }

  spill.dibfirst = dibfirst;
  spill.diblast = diblast;
  spill.dibmask = dibmask;
  spill.maskstatus = maskstatus;

  // sumdata::EventQuality vars, filled by DQSpillFlags
  spill.nmissingdcms = 0;
  spill.fracdcm3hits = 0.0;
  spill.nouttimehits = 0;
  spill.nnoisyapds = 0;

  spill.nmicroslices = 0;
  spill.ndcms = 0;

  spill.nmissingdcmslg = 0;  // LiveGeometry check
  spill.nbaddcmslg = 0;      // LiveGeometry check
  spill.dcmedgematchfrac = 0;

  // Need to flag each spill if the pot was good or bad
  // each slice will get the same tag

  // If there was no spill info for this slice mark it bad. log message this
  // should not happen
  if (spillPot.failedToGet()) {
    mf::LogError("CAFMaker") << "Spill Data not found for real data event";
    spill.isgoodspill = false;
  }
  // If we have a SpillData object, get all the variables
  else {
    spill.isgoodspill = spillPot->goodbeam;
    spill.spilltimesec = spillPot->spilltimesec;
    spill.spilltimensec = spillPot->spilltimensec;
    spill.gpsspilltimesec = spillPot->gpsspilltimensec;
    spill.gpsspilltimensec = spillPot->gpsspilltimensec;
    spill.deltaspilltimensec = spillPot->deltaspilltimensec;
    spill.spillpot = spillPot->spillpot;
    spill.hornI = spillPot->hornI;
    spill.isRHC = spillPot->isRHC;
    spill.is0HC = spillPot->is0HC;
    spill.isFHC = !spill.isRHC && !spill.is0HC;
    spill.widthx = spillPot->widthx;
    spill.widthy = spillPot->widthy;
    spill.posx = spillPot->posx;
    spill.posy = spillPot->posy;

    spill.intx.resize(spillPot->intx.size());
    spill.inty.resize(spillPot->inty.size());
    spill.bposx.resize(spillPot->bposx.size());
    spill.bposy.resize(spillPot->bposy.size());
    for (size_t i = 0; i < spillPot->intx.size(); ++i)
      spill.intx[i] = spillPot->intx[i];
    for (size_t i = 0; i < spillPot->inty.size(); ++i)
      spill.inty[i] = spillPot->inty[i];
    for (size_t i = 0; i < spillPot->bposx.size(); ++i)
      spill.bposx[i] = spillPot->bposx[i];
    for (size_t i = 0; i < spillPot->bposy.size(); ++i)
      spill.bposy[i] = spillPot->bposy[i];

    if (!fIsRealData) {
      // MC will flag each spill as good since it doesn't make bad spills
      spill.isgoodspill = true;
    } else {
      // NOTE: there is an implied factor of 1e12 when getting
      // IFDBSpillInfo POT.  Adding this back in to avoid confusion
      // between data and MC.
      spill.spillpot *= 1e12;
    }
  }

  // If there was no spill info for this slice mark it bad.
  // Log message this should not happen
  if (spillQual.failedToGet()) {
    mf::LogError("CAFMaker")
        << "Spill EventQuality not found for real data event";
  } else {
    spill.nmissingdcmslg = spillQual->nmissingdcmslg;
    spill.nmissingdcms = spillQual->nmissingdcms;
    spill.fracdcm3hits = spillQual->fracdcm3hits;
    spill.nouttimehits = spillQual->nouttimehits;
    spill.nnoisyapds = spillQual->nnoisyapds;
    spill.dcmedgematchfrac = spillQual->dcmedgematchfrac;
    spill.nmicroslices = spillQual->nmicroslices;
  }

  if (trigs.failedToGet()) {
    mf::LogError("CAFMaker") << "No trigger info found for this event";
  } else {
    bool once = true;
    for (const rawdata::RawTrigger& trig : *trigs) {
      // Units of trigger field are 500ns
      spill.livetime += trig.fTriggerRange_TriggerLength * 500e-9;
      if (once) {
        once = false;
        spill.trigger = trig.fTriggerMask_TriggerType;
      } else {
        if (spill.trigger != trig.fTriggerMask_TriggerType) {
          mf::LogError("CAFMaker") << "Multiple trigger types found in one "
                                      "Event. Storing the first value seen.";
        }
      }
    }
  }

  // dq::DAQEventSummary variables
  if (daq.failedToGet()) {
    mf::LogError("CAFMaker") << "No DAQ summary info found for this event.";
  } else {
    spill.ndiblocks = daq->fNDiblocks;
    spill.eventincomplete = daq->fEventIncomplete;
    spill.emptydatablock = daq->fNemptyDataBlock;
    spill.nmicroblocks = daq->fNmicroBlocks;
    spill.nemptymicroslice = daq->fNemptyMicroSlice;
    spill.ndroppedmicroblocks = daq->fNDroppedMicroBlocks;
    spill.ndatablockmissingdata = daq->fNdataBlockMissingData;
    spill.nmicroslicedatanotpresent = daq->fNmicroSliceDataNotPresent;
    spill.nnanoslices = daq->fNtotalNanoSlices;
    spill.nanoslicedatanotpresent = daq->fNanoSliceDataNotPresent;
    spill.nanoslicenolinkstatus = daq->fNanoSliceNoLinkStatus;
    spill.nanoslicebufferempty = daq->fNanoSliceBufferEmpty;
    spill.nanoslicebufferfull = daq->fNanoSliceBufferFull;
    spill.nanoslicecommerror = daq->fNanoSliceCommError;
    spill.nanoslicepacketerror = daq->fNanoSlicePacketError;
    spill.nanosliceoverflowerror = daq->fNanoSliceOverflowError;
    spill.nanosliceadcerror = daq->fNanoSliceADCError;
  }

  spill.nbaddcmslg = livegeom->NBadDCMBC();
  spill.ndcms = rh->NDCMs();
}

DEFINE_ART_MODULE(CAFMaker)

}  // end namespace caf
////////////////////////////////////////////////////////////////////////