Public Types | Public Member Functions | Static Public Member Functions | Protected Types | Protected Member Functions | Static Protected Member Functions | Protected Attributes | List of all members
caf::CAFMaker Class Reference

Module to create Common Analysis Files from ART files. More...

Inheritance diagram for caf::CAFMaker:
art::EDProducer art::ProducerBase art::Consumer art::EngineCreator art::ProductRegistryHelper

Public Types

using Parameters = art::EDProducer::Table< CAFMakerParams >
 
using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
template<typename UserConfig , typename KeysToIgnore = void>
using Table = ProducerBase::Table< UserConfig, KeysToIgnore >
 

Public Member Functions

 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)
 
template<typename PROD , BranchType B = InEvent>
ProductID getProductID (std::string const &instanceName={}) const
 
template<typename PROD , BranchType B>
ProductID getProductID (ModuleDescription const &moduleDescription, std::string const &instanceName) const
 
bool modifiesEvent () const
 
template<typename T , BranchType = InEvent>
ProductToken< Tconsumes (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< Tconsumes (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TconsumesView (InputTag const &it)
 
template<typename T , BranchType = InEvent>
ProductToken< TmayConsume (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ProductToken< TmayConsume (InputTag const &it)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , art::BranchType BT>
art::ViewToken< TmayConsumeView (InputTag const &it)
 
base_engine_tcreateEngine (seed_t seed)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make)
 
base_engine_tcreateEngine (seed_t seed, std::string const &kind_of_engine_to_make, label_t const &engine_label)
 
seed_t get_seed_value (fhicl::ParameterSet const &pset, char const key[]="seed", seed_t const implicit_seed=-1)
 

Static Public Member Functions

static cet::exempt_ptr< Consumernon_module_context ()
 

Protected Types

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

Protected Member Functions

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)
 
template<class T , class U >
art::FindManyP< TFindManyPStrict (const U &from, const art::Event &evt, const art::InputTag &label) const
 
template<class T >
bool GetAssociatedProduct (const art::FindManyP< T > &fm, int idx, T &ret) const
 Retrieve an object from an association, with error handling. More...
 
template<class T >
void GetByLabelStrict (const art::Event &evt, const std::string &label, art::Handle< T > &handle) const
 
template<class T >
void GetByLabelIfExists (const art::Event &evt, const std::string &label, art::Handle< T > &handle) const
 
template<class T >
bool GetPsetParameter (const fhicl::ParameterSet &pset, const std::vector< std::string > &name, T &ret) const
 
std::pair< int, intcalcFirstLastLivePlane (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
 
CurrentProcessingContext const * currentContext () const
 
void validateConsumedProduct (BranchType const bt, ProductInfo const &pi)
 
void prepareForJob (fhicl::ParameterSet const &pset)
 
void showMissingConsumes () const
 

Static Protected Member Functions

static bool EssentiallyEqual (double a, double b, double precision=0.0001)
 
static bool sortTrackLength (const SRTrack &a, const SRTrack &b)
 
static bool sortRBTrackLength (const art::Ptr< rb::Track > &a, const art::Ptr< rb::Track > &b)
 
static bool sortRBProngLength (const art::Ptr< rb::Prong > &a, const art::Ptr< rb::Prong > &b)
 
static bool sortShowerEnergy (const SRFuzzyKProng &a, const SRFuzzyKProng &b)
 
static bool sortProngCalE (const SRFuzzyKProng &a, const SRFuzzyKProng &b)
 
static bool sortElasticProng (const SRElastic &a, const SRElastic &b)
 
static bool sortHoughProng (const SRHoughVertex &a, const SRHoughVertex &b)
 
static bool sortMuIdLength (const SRMuId &a, const SRMuId &b)
 
static bool sortRemid (const SRKalmanTrack &a, const SRKalmanTrack &b)
 
static bool compMuonID (const SRKalmanTrack &a, const SRKalmanTrack &b)
 
static bool sortNuWithIdx (const cheat::NeutrinoWithIndex &a, const cheat::NeutrinoWithIndex &b)
 

Protected Attributes

CAFMakerParams fParams
 
std::string fCafFilename
 
bool fIsRealData
 
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::BackTrackerbt
 
art::ServiceHandle< geo::Geometrygeom
 
art::ServiceHandle< geo::LiveGeometrylivegeom
 
art::ServiceHandle< nova::dbi::RunHistoryServicerh
 
art::ServiceHandle< chaninfo::BadChanListbc
 
novadaq::cnv::DetId fDetID
 Detector ID in nova daq convention namespace typedef. More...
 
Det_t fDet
 Detector ID in caf namespace typedef. More...
 
TTree * fSpillTree
 
TTree * fCosmicCVNTree
 
TTree * fNuTree
 
std::map< std::string, MixingTypemapStringMixingType
 
std::map< MixingType, std::stringmapMixingTypeString
 

Detailed Description

Module to create Common Analysis Files from ART files.

Definition at line 162 of file CAFMaker_module.cc.

Member Typedef Documentation

using art::EDProducer::ModuleType = EDProducer
inherited

Definition at line 34 of file EDProducer.h.

Definition at line 165 of file CAFMaker_module.cc.

template<typename UserConfig , typename KeysToIgnore = void>
using art::EDProducer::Table = ProducerBase::Table<UserConfig, KeysToIgnore>
inherited

Definition at line 43 of file EDProducer.h.

using art::EDProducer::WorkerType = WorkerT<EDProducer>
inherited

Definition at line 35 of file EDProducer.h.

Member Enumeration Documentation

enum caf::CAFMaker::MixingType
strongprotected
Enumerator
eRockOverlay 
eCosmicOverlay 
eRockCosmicOverlay 
eRockSinglesOverlay 
eNone 
ePileup 
eSingles 
eOverlay 
eMCSpillSinglesOverlay 
eDataSpillSinglesOverlay 

Definition at line 348 of file CAFMaker_module.cc.

348  {
349  eRockOverlay,
350  eCosmicOverlay,
351  eRockCosmicOverlay,
352  eRockSinglesOverlay,
353  eNone,
354  ePileup,
355  eSingles,
356  eOverlay,
357  eMCSpillSinglesOverlay,
358  eDataSpillSinglesOverlay
359  };

Constructor & Destructor Documentation

caf::CAFMaker::CAFMaker ( const Parameters params)
explicit

Definition at line 389 of file CAFMaker_module.cc.

References caf::CAFMakerParams::CAFFilename, fCafFilename, caf::CAFMakerParams::FillNuTree, fParams, and art::ProductRegistryHelper::produces().

390  : fParams(params()), fIsRealData(false), fFile(0) {
391  // We update this, so have to pull it out of the config
393 
394  // Normally CAFMaker is run without an output ART stream, so these go
395  // nowhere, but can be occasionally useful for filtering as part of an
396  // ART job.
397  produces<std::vector<caf::StandardRecord>>();
398  produces<art::Assns<caf::StandardRecord, rb::Cluster>>();
399  produces<caf::SRSpill>();
400  if (fParams.FillNuTree()) produces<std::vector<caf::SRSpillTruthBranch>>();
401 }
void produces(std::string const &instanceName={}, Persistable const persistable=Persistable::Yes)
Atom< bool > FillNuTree
std::string fCafFilename
Atom< std::string > CAFFilename
CAFMakerParams fParams
caf::CAFMaker::~CAFMaker ( )
virtual

Definition at line 404 of file CAFMaker_module.cc.

404 {}

Member Function Documentation

void caf::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 
)
protected

Definition at line 605 of file CAFMaker_module.cc.

References caf::AddParticleToVec(), caf::AddPreFSI(), POTSpillRate::beam, caf::SRNeutrino::beam, bt, caf::CAFGeneratorEnum(), simb::MCNeutrino::CCNC(), caf::SRNeutrino::cycle, caf::CAFMakerParams::DataSpillHasMC, caf::DecodeGeneratorVersion(), caf::SRNeutrino::det, geo::GeometryBase::DetHalfHeight(), geo::GeometryBase::DetHalfWidth(), geo::GeometryBase::DetLength(), caf::SRNeutrino::E, simb::MCParticle::E(), caf::SRNeutrino::eff, cheat::NeutrinoEffPur::efficiency, cheat::NeutrinoEffPur::energySlice, cheat::NeutrinoEffPur::energyTotal, EssentiallyEqual(), art::EventID::event(), evt, caf::SRNeutrino::evt, stan::math::fabs(), art::Handle< T >::failedToGet(), fCycle, fDet, fDetID, simb::MCFlux::fdk2gen, simb::MCFlux::fevtno, simb::MCFlux::fgen2vtx, caf::FindAndAddMichels(), simb::GTruth::fIsCharm, simb::GTruth::fIsSeaQuark, flux, simb::MCFlux::fmupare, simb::MCFlux::fmuparpx, simb::MCFlux::fmuparpy, simb::MCFlux::fmuparpz, simb::MCFlux::fndecay, simb::MCFlux::fnecm, simb::MCFlux::fnimpwt, simb::MCFlux::fntype, simb::GTruth::fNumNeutron, simb::GTruth::fNumPi0, simb::GTruth::fNumPiMinus, simb::GTruth::fNumPiPlus, simb::GTruth::fNumProton, fParams, simb::MCFlux::fpdpx, simb::MCFlux::fpdpy, simb::MCFlux::fpdpz, simb::MCFlux::fppdxdz, simb::MCFlux::fppdydz, simb::MCFlux::fppenergy, simb::MCFlux::fppmedium, simb::MCFlux::fpppz, simb::MCFlux::fppvx, simb::MCFlux::fppvy, simb::MCFlux::fppvz, simb::MCFlux::fptype, simb::GTruth::fResNum, simb::MCFlux::frun, simb::MCFlux::ftgen, simb::MCFlux::ftgptype, simb::MCFlux::ftptype, simb::MCFlux::ftpx, simb::MCFlux::ftpy, simb::MCFlux::ftpz, simb::MCFlux::ftvx, simb::MCFlux::ftvy, simb::MCFlux::ftvz, simb::MCFlux::fvx, simb::MCFlux::fvy, simb::MCFlux::fvz, simb::GTruth::fXsec, caf::SRNeutrino::genConfigString, simb::MCGeneratorInfo::generator, caf::SRNeutrino::generator, simb::MCGeneratorInfo::generatorConfig, simb::MCTruth::GeneratorInfo(), caf::CAFMakerParams::GeneratorLabel, simb::MCGeneratorInfo::generatorVersion, caf::SRNeutrino::genVersion, caf::SRNeutrino::genweight, geom, GetByLabelStrict(), simb::MCTruth::GetNeutrino(), simb::MCTruth::GetParticle(), caf::SRNeutrino::hitnuc, simb::MCNeutrino::HitNuc(), caf::SRNeutrino::hitnucp, art::Event::id(), simb::MCNeutrino::InteractionType(), caf::SRNeutrino::inttype, ipart, caf::SRSpill::is0HC, caf::SRNeutrino::is0HC, caf::SRNeutrino::iscc, caf::SRNeutrino::ischarm, caf::SRSpill::isFHC, caf::SRNeutrino::isFHC, caf::SRSpill::isRHC, caf::SRNeutrino::isRHC, caf::SRNeutrino::isseaquark, caf::SRNeutrino::isvtxcont, it, simb::kCC, caf::kElectronScattering, caf::SRNeutrino::L, simb::MCNeutrino::Lepton(), cheat::BackTracker::MCTruthToParticles(), caf::SRNeutrino::michel, simb::MCNeutrino::Mode(), caf::SRNeutrino::mode, simb::MCParticle::Momentum(), simb::MCParticle::Mother(), cheat::NeutrinoEffPur::neutrinoInt, caf::SRNeutrino::nhitslc, caf::SRNeutrino::nhittot, caf::SRNeutrino::nneutron, simb::MCTruth::NParticles(), caf::SRNeutrino::npiminus, caf::SRNeutrino::npiplus, caf::SRNeutrino::npizero, caf::SRNeutrino::nproton, cheat::NeutrinoEffPur::nSliceHits, cheat::NeutrinoEffPur::nTotalHits, simb::MCNeutrino::Nu(), caf::SRNeutrino::p, part, caf::SRNeutrino::pdg, simb::MCParticle::PdgCode(), caf::SRNeutrino::pdgorig, simb::MCParticle::Position(), caf::SRNeutrino::prim, caf::SRNeutrino::pur, cheat::NeutrinoEffPur::purity, art::PtrVector< T >::push_back(), caf::SRNeutrino::q2, simb::MCNeutrino::QSqr(), caf::SRNeutrino::resnum, art::Event::run(), caf::SRNeutrino::run, SimpleOscProb(), simb::MCParticle::StatusCode(), string, art::Event::subRun(), caf::SRNeutrino::subrun, simb::MCNeutrino::Target(), caf::SRNeutrino::tgtA, caf::SRNeutrino::tgtZ, caf::SRNeutrino::time, cheat::BackTracker::TrackIDToMCTruth(), cheat::BackTracker::TrackIDToParticle(), caf::CAFMakerParams::TrueEnergyLabel, caf::TrueNeutrinoDistance(), caf::SRBeam::tv, caf::SRNeutrino::visE, caf::SRNeutrino::visEBirks, caf::SRNeutrino::visEinslc, caf::SRNeutrino::visEinslcBirks, caf::SRNeutrino::visEinslcNeutron, caf::SRNeutrino::visEinslcNeutronBirks, caf::SRNeutrino::visENeutron, caf::SRNeutrino::visENeutronBirks, caf::SRNeutrino::vtx, simb::MCNeutrino::W(), caf::SRNeutrino::W2, simb::MCParticle::Weight(), caf::SRNeutrino::woscdumb, caf::SRVector3D::X(), simb::MCNeutrino::X(), caf::SRNeutrino::x, caf::SRNeutrino::xsec, caf::SRVector3D::Y(), simb::MCNeutrino::Y(), caf::SRNeutrino::y, and caf::SRVector3D::Z().

Referenced by produce().

611  {
612  // Get information about truNu
613  vec->push_back(SRNeutrino());
614  SRNeutrino& srTruth = vec->back();
615 
616  srTruth.run = evt.run();
617  srTruth.subrun = evt.subRun();
618  srTruth.evt = evt.id().event();
619  srTruth.cycle = fCycle;
620 
621  srTruth.isFHC = spill.isFHC;
622  srTruth.is0HC = spill.is0HC;
623  srTruth.isRHC = spill.isRHC;
624 
625  srTruth.det = fDet;
626 
627  const art::Ptr<simb::MCTruth>& tru = effPur->neutrinoInt;
628 
629  srTruth.generator = CAFGeneratorEnum(tru->GeneratorInfo().generator);
630 
631  // this is very rudimentary caching
632  static std::string oldVersionStr = std::string{};
633  static std::vector<unsigned int> oldVersionVec = std::vector<unsigned int>{};
634  static std::unordered_map<std::string, std::string> oldGenConfig = std::unordered_map<std::string, std::string>{};
635  static std::string oldGenConfigStr = std::string{};
636  if (tru->GeneratorInfo().generatorVersion != oldVersionStr || tru->GeneratorInfo().generatorConfig != oldGenConfig)
637  {
638  oldVersionStr = tru->GeneratorInfo().generatorVersion;
639  oldVersionVec = DecodeGeneratorVersion(tru->GeneratorInfo().generatorVersion);
640  oldGenConfig = tru->GeneratorInfo().generatorConfig;
641  oldGenConfigStr = "";
642  for (const auto & configPair: oldGenConfig)
643  {
644  if (configPair.first == "tune")
645  oldGenConfigStr = configPair.second;
646  else
647  mf::LogWarning("CAFMaker") << "Don't know how to store generator config parameter named '" << configPair.first << "'";
648  }
649  }
650  srTruth.genVersion = oldVersionVec;
651  if (!oldGenConfigStr.empty())
652  srTruth.genConfigString = oldGenConfigStr;
653 
654 
655  srTruth.pdg = tru->GetNeutrino().Nu().PdgCode();
656  srTruth.hitnuc = tru->GetNeutrino().HitNuc();
657  // Get momentum of struck nucleon within the nucleus
658  for ( int ipart = 0; ipart < tru->NParticles(); ++ipart )
659  {
660  if ( tru->GetParticle( ipart ).StatusCode() == 11 )
661  {
662  srTruth.hitnucp = SRLorentzVector( tru->GetParticle( ipart ).Momentum() );
663  }
664  }
665 
666  srTruth.mode = tru->GetNeutrino().Mode(); // to avoid depending on simb
667  srTruth.iscc = (tru->GetNeutrino().CCNC() == simb::kCC);
668 
669  if (srTruth.mode == kElectronScattering) {
670  // For a numu or nutau scattering off an electron, the exchanged particle
671  // is a Z. But for a nue, there are two diagrams, one with a Z and one
672  // with a W. Both the final states have the electron ejected, so the
673  // interaction is a quantum-mechanical superposition between the
674  // two. It's neither really CC or NC. We do need to set this boolean
675  // though. Treating nu-on-e scattering as NC (ie not oscillating it) is
676  // wrong, the cross-section for nue is different to numu and nutau so
677  // oscillations do have an effect on the final spectrum. If we mark it CC
678  // the different oscillated fluxes will be weighted together in the fit
679  // appropriately, so this is safe, even if not perfectly physically
680  // accurate.
681  srTruth.iscc = true;
682  }
683  srTruth.inttype = tru->GetNeutrino().InteractionType();
684 
685  srTruth.q2 = tru->GetNeutrino().QSqr();
686  srTruth.x = tru->GetNeutrino().X();
687  srTruth.y = tru->GetNeutrino().Y();
688  srTruth.W2 = tru->GetNeutrino().W() * tru->GetNeutrino().W();
689  srTruth.E = tru->GetNeutrino().Nu().Momentum().E();
690  srTruth.visE = effPur->energyTotal;
691  srTruth.visEinslc = effPur->energySlice;
692  srTruth.eff = effPur->efficiency;
693  srTruth.pur = effPur->purity;
694  srTruth.nhitslc = effPur->nSliceHits;
695  srTruth.nhittot = effPur->nTotalHits;
696  srTruth.time = tru->GetNeutrino().Nu().Position().T();
697  srTruth.vtx = tru->GetNeutrino().Nu().Position().Vect();
698  srTruth.p = tru->GetNeutrino().Nu().Momentum();
699 
700  float visEBirks = 0;
701  for (std::vector<cheat::TrackIDE>::iterator it = allTracksBirks.begin();
702  it != allTracksBirks.end(); ++it) {
703  if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
704  visEBirks += it->energy;
705  }
706  float visEinslcBirks = 0;
707  for (std::vector<cheat::TrackIDE>::iterator it = sliceTracksBirks.begin();
708  it != sliceTracksBirks.end(); ++it) {
709  if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
710  visEinslcBirks += it->energy;
711  }
712 
713  srTruth.visEBirks = visEBirks;
714  srTruth.visEinslcBirks = visEinslcBirks;
715 
716  // Get tgtZ / tgtA from the MCTruth - which stores pdg code of nucleus
717  int tarpdg = tru->GetNeutrino().Target();
718  // pdgcode of nucleus = 10LZZZAAAI
719  srTruth.tgtZ = (tarpdg / 10000) % 1000;
720  srTruth.tgtA = (tarpdg / 10) % 1000;
721 
722  srTruth.genweight = tru->GetNeutrino().Nu().Weight();
723 
724  TVector3 nuVtx = tru->GetNeutrino().Nu().Position().Vect();
725  bool inDet = fabs(nuVtx.X()) <= geom->DetHalfWidth() &&
726  fabs(nuVtx.Y()) <= geom->DetHalfHeight() && nuVtx.Z() >= 0 &&
727  nuVtx.Z() <= geom->DetLength();
728  srTruth.isvtxcont = inDet;
729 
730  // rescue the parts of GTruth that we always want, put them in SRNeutrino
732  pv.push_back(tru);
734  FindManyPStrict<simb::GTruth>(pv, evt, fParams.GeneratorLabel());
735 
736  if (fmGtruth.isValid()) {
737  std::vector<art::Ptr<simb::GTruth>> gtruths = fmGtruth.at(0);
738 
739  if (!gtruths.empty()) {
740  const simb::GTruth& gtruth = *gtruths[0];
741 
742  // Useful for identifying pre-FSI generated hadron content
743  srTruth.npiplus = gtruth.fNumPiPlus;
744  srTruth.npiminus = gtruth.fNumPiMinus;
745  srTruth.npizero = gtruth.fNumPi0;
746  srTruth.nproton = gtruth.fNumProton;
747  srTruth.nneutron = gtruth.fNumNeutron;
748 
749  // Descriptive variables relating to the thrown cross sections
750  srTruth.ischarm = gtruth.fIsCharm;
751  srTruth.isseaquark = gtruth.fIsSeaQuark;
752  srTruth.resnum = gtruth.fResNum;
753  srTruth.xsec = gtruth.fXsec;
754  }
755  }
756 
757  // now stuff from the flux record
758  // art::FindManyP<simb::MCFlux> fmFlux(pv, evt, fParams.GeneratorLabel());
760  FindManyPStrict<simb::MCFlux>(pv, evt, fParams.GeneratorLabel());
761  if (fmFlux.isValid()) {
762  std::vector<art::Ptr<simb::MCFlux>> fluxes = fmFlux.at(0);
763 
764  if (!fluxes.empty()) {
765  const simb::MCFlux& flux = *fluxes[0];
766 
767  srTruth.pdgorig = flux.fntype;
768  srTruth.woscdumb = SimpleOscProb(flux, tru->GetNeutrino());
769 
770  SRBeam& beam = srTruth.beam;
771 
772  beam.tv = TVector3(flux.ftvx, flux.ftvy, flux.ftvz);
773  beam.tp = TVector3(flux.ftpx, flux.ftpy, flux.ftpz);
774  beam.runjob = flux.frun;
775  beam.potnum = flux.fevtno;
776  beam.tptype = flux.ftptype;
777  beam.nimpwt = flux.fnimpwt;
778  beam.ndecay = flux.fndecay;
779  beam.v = TVector3(flux.fvx, flux.fvy, flux.fvz);
780  beam.pdp = TVector3(flux.fpdpx, flux.fpdpy, flux.fpdpz);
781  beam.ppdxdz = flux.fppdxdz;
782  beam.ppdydz = flux.fppdydz;
783  beam.pppz = flux.fpppz;
784  beam.ppenergy = flux.fppenergy;
785  beam.ppmedium = flux.fppmedium;
786  beam.ptype = flux.fptype;
787  beam.ppv = TVector3(flux.fppvx, flux.fppvy, flux.fppvz);
788  beam.muparp = TVector3(flux.fmuparpx, flux.fmuparpy, flux.fmuparpz);
789  beam.mupare = flux.fmupare;
790  beam.necm = flux.fnecm;
791  beam.tgen = flux.ftgen;
792  beam.tgptype = flux.ftgptype;
793  beam.dk2gen = flux.fdk2gen;
794  beam.gen2vtx = flux.fgen2vtx;
795  beam.dk2vtx = flux.fdk2gen + flux.fgen2vtx;
796 
797  } else {
798  // this does not work for data-singlenu-mc overlay, but we don't want the
799  // singlenu flux anyway
800  if (!fParams.DataSpillHasMC())
801  throw cet::exception(
802  "Missing Association Error",
803  "Failed to find MCFlux object assoicated with MCTruth. This "
804  "probably means something went wrong with overlays.");
805  }
806  }
807 
808  // Must be after the setting of 'vtx' and 'beam' info
809  srTruth.L = TrueNeutrinoDistance(fDetID, srTruth);
810 
811  // Get the neutrino's daughters (and other progeny)
812  std::vector<const sim::Particle*> particles = bt->MCTruthToParticles(tru);
813 
814  // Get the sim::TrueEnergy's out of the event.
817  std::vector<sim::TrueEnergy> TrueEnergies;
818  if (!tes.failedToGet()) TrueEnergies = *tes;
819 
820  // We're going to loop over daughters to add other primaries to truNu.
821  // First we're gonna get the lepton and add it as the zeroth entry,
822  // so let's get its information.
823  int lepPdg = tru->GetNeutrino().Lepton().PdgCode();
824  double lepEnergy = tru->GetNeutrino().Lepton().E();
825  int lepCount = 0; // count the found leptons, for debugging.
826  int lepTrackId = 0; // We'll set this later
827  // Now we'll loop over the lepton and add it
828  float visENeutron = 0;
829  for (std::vector<cheat::TrackIDE>::iterator it = allTracks.begin();
830  it != allTracks.end(); ++it) {
831  if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
832  int mId = bt->TrackIDToParticle(it->trackID)->Mother();
833  while (mId != 0) {
834  if (bt->TrackIDToParticle(mId)->PdgCode() == 2112) {
835  visENeutron += it->energy;
836  break;
837  }
838  mId = bt->TrackIDToParticle(mId)->Mother();
839  }
840  }
841  float visEinslcNeutron = 0;
842  for (std::vector<cheat::TrackIDE>::iterator it = sliceTracks.begin();
843  it != sliceTracks.end(); ++it) {
844  if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
845  int mId = bt->TrackIDToParticle(it->trackID)->Mother();
846  while (mId != 0) {
847  if (bt->TrackIDToParticle(mId)->PdgCode() == 2112) {
848  visEinslcNeutron += it->energy;
849  break;
850  }
851  mId = bt->TrackIDToParticle(mId)->Mother();
852  }
853  }
854  float visENeutronBirks = 0;
855  for (std::vector<cheat::TrackIDE>::iterator it = allTracksBirks.begin();
856  it != allTracksBirks.end(); ++it) {
857  if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
858  int mId = bt->TrackIDToParticle(it->trackID)->Mother();
859  while (mId != 0) {
860  if (bt->TrackIDToParticle(mId)->PdgCode() == 2112) {
861  visENeutronBirks += it->energy;
862  break;
863  }
864  mId = bt->TrackIDToParticle(mId)->Mother();
865  }
866  }
867  float visEinslcNeutronBirks = 0;
868  for (std::vector<cheat::TrackIDE>::iterator it = sliceTracksBirks.begin();
869  it != sliceTracksBirks.end(); ++it) {
870  if (bt->TrackIDToMCTruth(it->trackID) != tru) continue;
871  int mId = bt->TrackIDToParticle(it->trackID)->Mother();
872  while (mId != 0) {
873  if (bt->TrackIDToParticle(mId)->PdgCode() == 2112) {
874  visEinslcNeutronBirks += it->energy;
875  break;
876  }
877  mId = bt->TrackIDToParticle(mId)->Mother();
878  }
879  }
880  srTruth.visENeutron = visENeutron;
881  srTruth.visEinslcNeutron = visEinslcNeutron;
882  srTruth.visENeutronBirks = visENeutronBirks;
883  srTruth.visEinslcNeutronBirks = visEinslcNeutronBirks;
884 
885  for (const sim::Particle* part : particles) {
886  // Add particle if it's a primary and it has the correct PDG code
887  // and the correct energy.
888  if (part->Mother() == 0 && part->PdgCode() == lepPdg &&
889  EssentiallyEqual(part->E(), lepEnergy, .01)) {
890  lepTrackId = part->TrackId();
891  ++lepCount;
892  AddParticleToVec(*part, allTracks, sliceTracks, allTracksBirks,
893  sliceTracksBirks, &srTruth.prim, TrueEnergies);
894  }
895  }
896 
897  if (lepCount != 1) {
898  switch (lepCount) {
899  case 0:
900  if (!srTruth.iscc) break;
901 
902  mf::LogError("CAFMaker")
903  << "*****************************************"
904  << "\n"
905  << "Uh, oh, CAFMaker did not find primary lepton. "
906  << "\n"
907  << "We're going to need to find a way to deal with this."
908  << "\n"
909  << "Info about the particle: "
910  << "\n"
911  << " pdg: " << srTruth.pdg << "\n"
912  << " pdgorig: " << srTruth.pdgorig << "\n"
913  << " E: " << srTruth.E << "\n"
914  << " visE: " << srTruth.visE << "\n"
915  << " nhittot: " << srTruth.nhittot << "\n"
916  << " X: " << srTruth.vtx.X() << "\n"
917  << " Y: " << srTruth.vtx.Y() << "\n"
918  << " Z: " << srTruth.vtx.Z() << "\n"
919  << " iscc: " << srTruth.iscc << "\n"
920  << " infid(nd): "
921  << (srTruth.vtx.X() < 200 && srTruth.vtx.X() > -200 &&
922  srTruth.vtx.Y() < 200 && srTruth.vtx.Y() > -200 &&
923  srTruth.vtx.Z() > 0 && srTruth.vtx.X() < 1600)
924  << "\n"
925 
926  << "*****************************************"
927  << "\n";
928  break;
929  default:
930  mf::LogError("CAFMaker")
931  << "*****************************************"
932  << "\n"
933  << "Uh, oh, CAFMaker found too many primary leptons in event. "
934  << "\n"
935  << "We're going to need to find a way to deal with this."
936  << "\n"
937  << "*****************************************";
938  break;
939  }
940  // assert(lepCount == 1);
941  }
942 
943  for (const sim::Particle* part : particles) {
944  // Add particle if it's a primary and it's not the lepton
945  if (part->TrackId() != lepTrackId && part->Mother() == 0)
946  AddParticleToVec(*part, allTracks, sliceTracks, allTracksBirks,
947  sliceTracksBirks, &srTruth.prim, TrueEnergies);
948  }
949 
950  FindAndAddMichels(particles, allTracks, &srTruth.michel);
951 
952  AddPreFSI( tru, srTruth );
953 }
double E(const int i=0) const
Definition: MCParticle.h:232
Atom< string > TrueEnergyLabel
caf::generator_ CAFGeneratorEnum(simb::Generator_t simbGeneratorEnum)
Definition: FillTruth.cxx:905
int frun
Definition: MCFlux.h:35
double TrueNeutrinoDistance(novadaq::cnv::DetId det, const SRNeutrino &nu)
Definition: FillTruth.cxx:817
SubRunNumber_t subRun() const
Definition: Event.h:72
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
art::ServiceHandle< cheat::BackTracker > bt
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
double QSqr() const
Definition: MCNeutrino.h:157
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
double fpdpx
Definition: MCFlux.h:55
art::ServiceHandle< geo::Geometry > geom
set< int >::iterator it
double ftpx
Definition: MCFlux.h:79
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Int_t ipart
Definition: Style.C:10
double Weight() const
Definition: MCParticle.h:253
void AddPreFSI(const art::Ptr< simb::MCTruth > &truth, SRNeutrino &nu)
Definition: FillTruth.cxx:726
void AddParticleToVec(const sim::Particle &part, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, std::vector< SRTrueParticle > *vec, const std::vector< sim::TrueEnergy > &TrueEnergies)
Definition: FillTruth.cxx:258
int fppmedium
Definition: MCFlux.h:62
double fgen2vtx
distance from ray origin to event vtx
Definition: MCFlux.h:104
double SimpleOscProb(const simb::MCFlux &flux, const simb::MCNeutrino &nu) const
int Mother() const
Definition: MCParticle.h:212
void GetByLabelStrict(const art::Event &evt, const std::string &label, art::Handle< T > &handle) const
int ftptype
Definition: MCFlux.h:82
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
int fNumNeutron
number of neutrons after reaction, before FSI
Definition: GTruth.h:76
double fvx
Definition: MCFlux.h:52
int HitNuc() const
Definition: MCNeutrino.h:152
double DetLength() const
double fXsec
cross section of interaction
Definition: GTruth.h:30
int nTotalHits
Total number of hits the neutrino left in the detector.
Definition: BackTracker.h:54
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
int fNumPiPlus
number of pi pluses after reaction, before FSI
Definition: GTruth.h:78
double energyTotal
Sum of FLS hits from the neutrino contributing to any hit in the event.
Definition: BackTracker.h:52
int StatusCode() const
Definition: MCParticle.h:210
int fNumPiMinus
number of pi minuses after reaction, before FSI
Definition: GTruth.h:79
double fppdxdz
Definition: MCFlux.h:58
Atom< string > GeneratorLabel
int NParticles() const
Definition: MCTruth.h:74
int ftgen
Definition: MCFlux.h:83
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
double fdk2gen
distance from decay to ray origin
Definition: MCFlux.h:103
double efficiency
Efficiency (based on FLS energy) of neutrino interaction relative to slice.
Definition: BackTracker.h:48
Loaders::FluxType flux
Det_t fDet
Detector ID in caf namespace typedef.
novadaq::cnv::DetId fDetID
Detector ID in nova daq convention namespace typedef.
int ftgptype
Definition: MCFlux.h:84
int fResNum
resonance number
Definition: GTruth.h:80
int fNumProton
number of protons after reaction, before FSI
Definition: GTruth.h:75
int InteractionType() const
Definition: MCNeutrino.h:150
TString part[npart]
Definition: Style.C:32
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
double W() const
Definition: MCNeutrino.h:154
int fptype
Definition: MCFlux.h:63
double ftvx
Definition: MCFlux.h:76
void FindAndAddMichels(std::vector< const sim::Particle * > particles, std::vector< cheat::TrackIDE > &allTracks, std::vector< SRTrueMichelE > *michelVec)
Definition: FillTruth.cxx:664
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
double Y() const
Definition: MCNeutrino.h:156
Atom< bool > DataSpillHasMC
int fndecay
Definition: MCFlux.h:50
double fpdpz
Definition: MCFlux.h:57
int fNumPi0
number of pi0 after reaction, before FSI
Definition: GTruth.h:77
simb::Generator_t generator
event generator that generated this event
int evt
bool fIsCharm
did the interaction produce a charmed hadron?
Definition: GTruth.h:71
double fmuparpz
Definition: MCFlux.h:69
double X() const
Definition: MCNeutrino.h:155
int fevtno
Definition: MCFlux.h:36
Eigen::VectorXd vec
double fppdydz
Definition: MCFlux.h:59
double ftpz
Definition: MCFlux.h:81
double ftvy
Definition: MCFlux.h:77
double DetHalfHeight() const
std::vector< unsigned int > DecodeGeneratorVersion(const std::string &versionString)
Definition: FillTruth.cxx:922
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
double fpdpy
Definition: MCFlux.h:56
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
EventNumber_t event() const
Definition: EventID.h:116
std::string generatorVersion
event generator version
double fvy
Definition: MCFlux.h:53
double fmupare
Definition: MCFlux.h:70
int Target() const
Definition: MCNeutrino.h:151
int nSliceHits
Number of hits from this neutrino in this slice.
Definition: BackTracker.h:53
double DetHalfWidth() const
const art::Ptr< simb::MCTruth > & TrackIDToMCTruth(int const &id) const
double fppenergy
Definition: MCFlux.h:61
CAFMakerParams fParams
int fntype
Definition: MCFlux.h:51
double fnecm
Definition: MCFlux.h:71
double energySlice
Sum of FLS hits from the neutrino contributing to hits included in the slice.
Definition: BackTracker.h:51
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double fmuparpy
Definition: MCFlux.h:68
double fppvx
Definition: MCFlux.h:64
double fpppz
Definition: MCFlux.h:60
double fppvy
Definition: MCFlux.h:65
double fppvz
Definition: MCFlux.h:66
double ftpy
Definition: MCFlux.h:80
double fvz
Definition: MCFlux.h:54
RunNumber_t run() const
Definition: Event.h:77
int Mode() const
Definition: MCNeutrino.h:149
bool fIsSeaQuark
Definition: GTruth.h:50
std::unordered_map< std::string, std::string > generatorConfig
free-form field that can be used to keep track of generator configuration (e.g. GENIE tune) ...
static bool EssentiallyEqual(double a, double b, double precision=0.0001)
EventID id() const
Definition: Event.h:56
double fmuparpx
Definition: MCFlux.h:67
bool failedToGet() const
Definition: Handle.h:196
art::Ptr< simb::MCTruth > neutrinoInt
Truth information about neutrino interaction.
Definition: BackTracker.h:47
double ftvz
Definition: MCFlux.h:78
double fnimpwt
Definition: MCFlux.h:72
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const
double purity
Purity (based on FLS energy) of neutrino interaction relative to slice.
Definition: BackTracker.h:49
const simb::MCGeneratorInfo & GeneratorInfo() const
Definition: MCTruth.h:72
enum BeamMode string
void caf::CAFMaker::AddMetadataToFile ( std::map< std::string, std::string metadata)
protected

Definition at line 594 of file CAFMaker_module.cc.

References ana::assert(), fFile, and Write().

Referenced by endJob().

594  {
595  assert(fFile && "CAFMaker: Trying to add metadata to an uninitialized file");
596 
597  fFile->mkdir("metadata")->cd();
598 
599  for (auto const& pair : metadata) {
600  TObjString(pair.second.c_str()).Write(pair.first.c_str());
601  }
602 }
assert(nhit_max >=nhit_nbins)
gm Write()
void caf::CAFMaker::beginJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 424 of file CAFMaker_module.cc.

References fCafFilename, and InitializeOutfile().

424  {
425  if (!fCafFilename.empty()) InitializeOutfile();
426 }
void InitializeOutfile()
std::string fCafFilename
void caf::CAFMaker::beginRun ( art::Run r)
virtual

Reimplemented from art::EDProducer.

Definition at line 429 of file CAFMaker_module.cc.

References caf::CAFMakerParams::ApplyingFilter, geo::GeometryBase::DetId(), fDet, fDetID, fMask, fParams, geom, art::DataViewImpl::getByLabel(), nova::dbi::RunHistory::GetDiBlockMaskFromCondb(), nova::dbi::RunHistory::GoodDiBlockMask(), caf::CAFMakerParams::HitsLabel, art::Handle< T >::isValid(), lem_server::mask, and rh.

429  {
430  fDetID = geom->DetId();
431  fDet = (Det_t)fDetID;
432 
433  if (!rh->GetDiBlockMaskFromCondb() && !fParams.ApplyingFilter()) { // get for each run; here if subrun
434  // masking not used and there is only
435  // one stream with calhit run
437  r.getByLabel(fParams.HitsLabel(), mask); // why doesn't fParams.HitInput() work here?
438  if (!mask.isValid()) {
439 
440  mf::LogError("CAFMaker") << "Could not find a diBlock mask from " << fParams.HitsLabel() << ". Bailing.";
441  abort();
442  }
443 
444  fMask = *mask;
445  }
446 
447 
448  if (!rh->GetDiBlockMaskFromCondb() && fParams.ApplyingFilter()) { // get for each run; here if subrun
449  // masking not used and there are two
450  // streams one of which does not have
451  // calhits to save space and time
452  std::unique_ptr<int> mask(new int(rh->GoodDiBlockMask(0,true)));
453 
454  fMask = *mask;
455  }
456 }
Det_t
Which NOvA detector?
Definition: SREnums.h:7
art::ServiceHandle< geo::Geometry > geom
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Atom< bool > ApplyingFilter
Det_t fDet
Detector ID in caf namespace typedef.
novadaq::cnv::DetId fDetID
Detector ID in nova daq convention namespace typedef.
bool isValid() const
Definition: Handle.h:189
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
Atom< string > HitsLabel
CAFMakerParams fParams
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
bool GetDiBlockMaskFromCondb() const
Definition: RunHistory.h:368
int GoodDiBlockMask(int subrun=-1, bool reload=false)
art::ServiceHandle< nova::dbi::RunHistoryService > rh
void caf::CAFMaker::beginSubRun ( art::SubRun sr)
virtual

Reimplemented from art::EDProducer.

Definition at line 459 of file CAFMaker_module.cc.

References caf::CAFMakerParams::ApplyingFilter, fMask, fParams, art::DataViewImpl::getByLabel(), nova::dbi::RunHistory::GetDiBlockMaskFromCondb(), nova::dbi::RunHistory::GoodDiBlockMask(), caf::CAFMakerParams::HitsLabel, art::Handle< T >::isValid(), lem_server::mask, rh, and art::SubRun::subRun().

459  {
460  if (rh->GetDiBlockMaskFromCondb() && !fParams.ApplyingFilter()) { // get if we are using subrun-level
461  // masks and there is only
462  // one stream with calhit run
464  sr.getByLabel(fParams.HitsLabel(), mask); // why doesn't fParams.HitInput() work here?
465  if (!mask.isValid()) {
466 
467  mf::LogError("CAFMaker") << "Could not find a diBlock mask from " << fParams.HitsLabel() << ". Bailing.";
468  abort();
469  }
470 
471  fMask = *mask;
472  }
473 
474  if (rh->GetDiBlockMaskFromCondb() && fParams.ApplyingFilter()) { // get if we are using subrun-level
475  // masks and there are two
476  // streams one of which does not have
477  // calhits to save space and time
478  std::unique_ptr<int> mask(new int(rh->GoodDiBlockMask(sr.subRun(),false)));
479 
480  fMask = *mask;
481  }
482 }
SubRunNumber_t subRun() const
Definition: SubRun.h:44
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
Atom< bool > ApplyingFilter
bool isValid() const
Definition: Handle.h:189
Atom< string > HitsLabel
CAFMakerParams fParams
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
bool GetDiBlockMaskFromCondb() const
Definition: RunHistory.h:368
int GoodDiBlockMask(int subrun=-1, bool reload=false)
art::ServiceHandle< nova::dbi::RunHistoryService > rh
std::pair< int, int > caf::CAFMaker::calcFirstLastLivePlane ( int  plane,
std::bitset< 14 >  binary,
caf::Det_t  det = caf::kFARDET 
)
protected

Definition at line 1078 of file CAFMaker_module.cc.

References MECModelEnuComparisons::i, and caf::kFARDET.

Referenced by produce(), and sortNuWithIdx().

1080  {
1081  // if not the FD, return full ND size
1082  if (det != caf::kFARDET) return std::pair<int, int>(0, 213);
1083  int testDB = (plane / 64);
1084  int minblock = testDB;
1085  int maxblock = testDB;
1086  // find block boundaries;
1087  for (int i = testDB - 1; i >= 0; --i) {
1088  if (binary[i])
1089  minblock = i;
1090  else
1091  break;
1092  }
1093  for (int i = testDB + 1; i < 14; ++i) {
1094  if (binary[i])
1095  maxblock = i;
1096  else
1097  break;
1098  }
1099  return std::pair<int, int>(64 * (minblock), 64 * (maxblock + 1) - 1);
1100 }
Far Detector at Ash River.
Definition: SREnums.h:11
static bool caf::CAFMaker::compMuonID ( const SRKalmanTrack a,
const SRKalmanTrack b 
)
inlinestaticprotected

Definition at line 327 of file CAFMaker_module.cc.

References caf::SRKalmanTrack::muonid.

Referenced by produce().

327  {
328  return a.muonid < b.muonid; // Use standard comparitor syntax for std::max_element
329  }
const double a
const hit & b
Definition: hits.cxx:21
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::consumes ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::consumes ( InputTag const &  it)
inherited

Definition at line 146 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

147 {
148  if (!moduleContext_)
149  return ProductToken<T>::invalid();
150 
151  consumables_[BT].emplace_back(ConsumableType::Product,
152  TypeID{typeid(T)},
153  it.label(),
154  it.instance(),
155  it.process());
156  return ProductToken<T>{it};
157 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::consumesMany ( )
inherited

Definition at line 161 of file Consumer.h.

References T.

162 {
163  if (!moduleContext_)
164  return;
165 
166  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
167 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::consumesView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::consumesView ( InputTag const &  it)
inherited

Definition at line 171 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

172 {
173  if (!moduleContext_)
174  return ViewToken<T>::invalid();
175 
176  consumables_[BT].emplace_back(ConsumableType::ViewElement,
177  TypeID{typeid(T)},
178  it.label(),
179  it.instance(),
180  it.process());
181  return ViewToken<T>{it};
182 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make 
)
inherited
base_engine_t& art::EngineCreator::createEngine ( seed_t  seed,
std::string const &  kind_of_engine_to_make,
label_t const &  engine_label 
)
inherited
CurrentProcessingContext const* art::EDProducer::currentContext ( ) const
protectedinherited
void caf::CAFMaker::drawSliceTruthTable ( const std::vector< std::vector< cheat::NeutrinoEffPur >> &  sliceTruthTable,
const std::vector< int > &  faveIds 
)
protected

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

Definition at line 961 of file CAFMaker_module.cc.

963  {
964  std::vector<int> otherFaveIds(faveIds.size(), -1);
965  drawSliceTruthTable(sliceTruthTable, faveIds, otherFaveIds);
966 }
void drawSliceTruthTable(const std::vector< std::vector< cheat::NeutrinoEffPur >> &sliceTruthTable, const std::vector< int > &faveIds)
void caf::CAFMaker::drawSliceTruthTable ( const std::vector< std::vector< cheat::NeutrinoEffPur >> &  sliceTruthTable,
const std::vector< int > &  faveIds,
const std::vector< int > &  otherFaveIds 
)
protected

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

Definition at line 975 of file CAFMaker_module.cc.

References om::cout, allTimeWatchdog::endl, and PandAna.Demos.tute_pid_validation::slc.

977  {
978  std::cout.precision(2);
979  std::cout.width(2);
980  std::cout.setf(std::ios::fixed, std::ios::floatfield);
981  int iSlc = 0;
982  for (auto slc : sliceTruthTable) {
983  int iNu = 0;
984  for (auto nu = slc.begin(); nu != slc.end(); ++nu, ++iNu) {
985  std::cout << nu->efficiency;
986  if (faveIds[iSlc] == iNu)
987  std::cout << "*";
988  else
989  std::cout << " ";
990  if (otherFaveIds[iSlc] == iNu)
991  std::cout << "@";
992  else
993  std::cout << " ";
994  std::cout << " ";
995  }
996  std::cout << std::endl;
997  }
998 }
OStream cout
Definition: OStream.cxx:6
void caf::CAFMaker::endJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 3372 of file CAFMaker_module.cc.

References caf::CAFMakerParams::AbortOnEmpty, AddMetadataToFile(), compare(), caf::CAFMakerParams::DataTier, caf::CAFMakerParams::EnableBlindness, allTimeWatchdog::endl, fDetID, fFile, fIsRealData, fParams, fTotalEvents, fTotalLivetime, fTotalPOT, fTotalSinglePOT, fTotalTrueNonswapNue, fTotalTrueSingleNue, meta::MetadataManager::getInstance(), meta::MetadataManager::GetMetadata(), hEvents, hLivetime, hPOT, hSinglePOT, hTotalTrueNonswapNue, hTotalTrueSingleNue, caf::CAFMakerParams::IsSinglesOverlay, novadaq::cnv::kFARDET, runNovaSAM::metadata, ss, string, and getGoodRuns4SAM::tag.

3372  {
3373  bool EmptyFile = false;
3374  if (fTotalEvents == 0) {
3375  // Among other things this makes us unable to determine the file_type
3376  // metdata field.
3377  if (fParams.AbortOnEmpty()) {
3378  mf::LogError("CAFMaker") << "No events processed in this file. Aborting rather than "
3379  << "produce an empty CAF."
3380  << std::endl;
3381  abort();
3382  }
3383  else {
3384  mf::LogWarning("CAFMaker") << "Making an empty CAF file." << std::endl;
3385  EmptyFile = true;
3386  }
3387  }
3388 
3389  // Make sure the recTree is safely in the file before starting on the POT
3390  // etc histograms. This way, absence of the histograms should be a clear
3391  // sign something went wrong with the tree.
3392  fFile->Write();
3393 
3394  fFile->cd();
3395  hPOT->Fill(.5, fTotalPOT);
3396  hEvents->Fill(.5, fTotalEvents);
3397  hLivetime->Fill(.5, fTotalLivetime);
3398 
3399  hPOT->Write();
3400  hEvents->Write();
3401  hLivetime->Write();
3402  fFile->Write();
3403 
3404  if (fParams.IsSinglesOverlay()) {
3407  hSinglePOT->Fill(0.5, fTotalSinglePOT);
3408 
3409  hSinglePOT->Write();
3410  hTotalTrueNonswapNue->Write();
3411  hTotalTrueSingleNue->Write();
3412  }
3413 
3414  std::map<std::string, std::string> metadata =
3416  metadata["data_tier"] = fParams.DataTier();
3417  if (EmptyFile)
3418  metadata["file_type"] = "empty";
3419  else if (fIsRealData)
3420  metadata["file_type"] = "importedDetector";
3421  else
3422  metadata["file_type"] = "importedSimulated";
3423 
3424  // Gymnastics to mark "restrictedbox" in nova.special for FD numi data
3425  if (metadata.count("online.stream")) {
3426  // Then this is a real data file and the stream exists in the metadata
3427  if (!metadata["online.stream"].compare("0") &&
3429  // Then this is a numi trigger file from the far detector
3430  if (!fParams.EnableBlindness()) {
3431  // Then we are not blinding it, need to mark it as such in the metadata
3432  std::string tag("restrictedbox");
3433 
3434  if (metadata.count("nova.special") &&
3435  std::string("").compare(metadata["nova.special"]) &&
3436  std::string("none").compare(
3437  metadata["nova.special"])) { // Then the metadata is not the
3438  // same as none, we need to append
3439  // to it
3440  // Use a stringstream to jam the strings together and append the tag
3441 
3442  std::stringstream ss;
3443  ss.str();
3444  ss << metadata["nova.special"] << "-" << tag;
3445  metadata["nova.special"] = ss.str();
3446  } else
3447  metadata["nova.special"] = tag;
3448  }
3449  }
3450  }
3451 
3452  AddMetadataToFile(metadata);
3453 }
static MetadataManager & getInstance()
Float_t ss
Definition: plot.C:24
Atom< bool > IsSinglesOverlay
std::map< std::string, std::string > & GetMetadata()
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
novadaq::cnv::DetId fDetID
Detector ID in nova daq convention namespace typedef.
TH1D * hTotalTrueSingleNue
Atom< string > DataTier
TH1D * hTotalTrueNonswapNue
Atom< bool > EnableBlindness
Far Detector at Ash River, MN.
Atom< bool > AbortOnEmpty
CAFMakerParams fParams
void AddMetadataToFile(std::map< std::string, std::string > metadata)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
bool compare(const GFluxGenerator &g1, const GFluxGenerator &g2)
enum BeamMode string
void caf::CAFMaker::endSubRun ( art::SubRun sr)
virtual

Reimplemented from art::EDProducer.

Definition at line 3335 of file CAFMaker_module.cc.

References caf::CAFMakerParams::ExposureLabel, art::Handle< T >::failedToGet(), fIsRealData, fParams, fTotalLivetime, fTotalPOT, caf::CAFMakerParams::GeneratorLabel, art::DataViewImpl::getByLabel(), caf::CAFMakerParams::NuMIBeamLabel, pot, art::SubRun::subRun(), sumdata::POTSum::totgoodpot, and sumdata::CosmicExposure::totlivetime.

3335  {
3337 
3338  // If this is MC get the spill info from the generator
3339  if (!fIsRealData) {
3341  }
3342  // If we have real data try to get the spill info from the NuMI beam
3343  else {
3345  }
3346 
3347  if (!pot.failedToGet()) {
3348  // NOTE: The POT that comes out of IFDBSpillInfo has an implied
3349  // factor of 1E12 that is being left off.
3350  // For now the data POT is being multiplied by this factor so that
3351  // the number is consistent with what we get in MC.
3352  // TODO: Is this what we want to do in the future?
3353  float subrunpot = pot->totgoodpot;
3354  if (fIsRealData) subrunpot *= 1e12;
3355  fTotalPOT += subrunpot;
3356 
3357  mf::LogInfo("CAFMaker")
3358  << "POT in subrun " << sr.subRun() << " = " << pot->totgoodpot;
3359  }
3360 
3362  sr.getByLabel(fParams.ExposureLabel(), expo);
3363 
3364  if (!expo.failedToGet()) {
3365  fTotalLivetime += expo->totlivetime;
3366  mf::LogInfo("CAFMaker")
3367  << "Exposure in subrun " << sr.subRun() << " = " << expo->totlivetime;
3368  }
3369 }
SubRunNumber_t subRun() const
Definition: SubRun.h:44
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Atom< string > NuMIBeamLabel
Atom< string > GeneratorLabel
Atom< string > ExposureLabel
#define pot
CAFMakerParams fParams
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double totgoodpot
normalized by 10^12 POT
Definition: POTSum.h:28
bool failedToGet() const
Definition: Handle.h:196
static bool caf::CAFMaker::EssentiallyEqual ( double  a,
double  b,
double  precision = 0.0001 
)
inlinestaticprotected

Definition at line 285 of file CAFMaker_module.cc.

Referenced by AddMCTruthToVec().

285  {
286  return a <= (b + precision) && a >= (b - precision);
287  }
const hit & b
Definition: hits.cxx:21
void caf::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
protected

Definition at line 3474 of file CAFMaker_module.cc.

References sumdata::SpillData::bposx, caf::SRSpill::bposx, sumdata::SpillData::bposy, caf::SRSpill::bposy, sumdata::EventQuality::dcmedgematchfrac, caf::SRSpill::dcmedgematchfrac, DEFINE_ART_MODULE(), sumdata::SpillData::deltaspilltimensec, caf::SRSpill::deltaspilltimensec, caf::SRSpill::det, caf::SRSpill::dibfirst, caf::SRSpill::diblast, caf::SRSpill::dibmask, e, caf::SRSpill::emptydatablock, art::EventID::event(), caf::SRSpill::eventincomplete, caf::SRSpill::evt, art::Handle< T >::failedToGet(), fDet, dq::DAQEventSummary::fEventIncomplete, fIsRealData, fMask, dq::DAQEventSummary::fNanoSliceADCError, dq::DAQEventSummary::fNanoSliceBufferEmpty, dq::DAQEventSummary::fNanoSliceBufferFull, dq::DAQEventSummary::fNanoSliceCommError, dq::DAQEventSummary::fNanoSliceDataNotPresent, dq::DAQEventSummary::fNanoSliceNoLinkStatus, dq::DAQEventSummary::fNanoSliceOverflowError, dq::DAQEventSummary::fNanoSlicePacketError, dq::DAQEventSummary::fNdataBlockMissingData, dq::DAQEventSummary::fNDiblocks, dq::DAQEventSummary::fNDroppedMicroBlocks, dq::DAQEventSummary::fNemptyDataBlock, dq::DAQEventSummary::fNemptyMicroSlice, dq::DAQEventSummary::fNmicroBlocks, dq::DAQEventSummary::fNmicroSliceDataNotPresent, dq::DAQEventSummary::fNtotalNanoSlices, for, sumdata::EventQuality::fracdcm3hits, caf::SRSpill::fracdcm3hits, sumdata::SpillData::goodbeam, nova::dbi::RunHistory::GoodDiBlockMask(), sumdata::SpillData::gpsspilltimensec, caf::SRSpill::gpsspilltimensec, caf::SRSpill::gpsspilltimesec, sumdata::SpillData::hornI, caf::SRSpill::hornI, MECModelEnuComparisons::i, art::Event::id(), if(), sumdata::SpillData::intx, caf::SRSpill::intx, sumdata::SpillData::inty, caf::SRSpill::inty, sumdata::SpillData::is0HC, caf::SRSpill::is0HC, caf::SRSpill::isFHC, caf::SRSpill::isgoodspill, caf::SRSpill::ismc, sumdata::SpillData::isRHC, caf::SRSpill::isRHC, livegeom, caf::SRSpill::livetime, caf::SRSpill::maskstatus, caf::SRSpill::nanosliceadcerror, caf::SRSpill::nanoslicebufferempty, caf::SRSpill::nanoslicebufferfull, caf::SRSpill::nanoslicecommerror, caf::SRSpill::nanoslicedatanotpresent, caf::SRSpill::nanoslicenolinkstatus, caf::SRSpill::nanosliceoverflowerror, caf::SRSpill::nanoslicepacketerror, geo::LiveGeometry::NBadDCMBC(), caf::SRSpill::nbaddcmslg, caf::SRSpill::ndatablockmissingdata, caf::SRSpill::ndcms, nova::dbi::RunHistory::NDCMs(), caf::SRSpill::ndiblocks, caf::SRSpill::ndroppedmicroblocks, caf::SRSpill::nemptymicroslice, caf::SRSpill::nmicroblocks, caf::SRSpill::nmicroslicedatanotpresent, sumdata::EventQuality::nmicroslices, caf::SRSpill::nmicroslices, sumdata::EventQuality::nmissingdcms, caf::SRSpill::nmissingdcms, sumdata::EventQuality::nmissingdcmslg, caf::SRSpill::nmissingdcmslg, caf::SRSpill::nnanoslices, sumdata::EventQuality::nnoisyapds, caf::SRSpill::nnoisyapds, sumdata::EventQuality::nouttimehits, caf::SRSpill::nouttimehits, sumdata::SpillData::posx, caf::SRSpill::posx, sumdata::SpillData::posy, caf::SRSpill::posy, rh, caf::SRSpill::run, art::Event::run(), sumdata::SpillData::spillpot, caf::SRSpill::spillpot, sumdata::SpillData::spilltimensec, caf::SRSpill::spilltimensec, sumdata::SpillData::spilltimesec, caf::SRSpill::spilltimesec, caf::SRSpill::subrun, art::Event::subRun(), caf::SRSpill::trigger, sumdata::SpillData::widthx, caf::SRSpill::widthx, sumdata::SpillData::widthy, and caf::SRSpill::widthy.

Referenced by produce(), and sortNuWithIdx().

3479  {
3480  // Add in run, subrun, evt information so that the spill tree can be
3481  // mapped to the record tree.
3482  spill.run = evt.run();
3483  spill.subrun = evt.subRun();
3484  spill.evt = evt.id().event();
3485 
3486  // If we fail to identify a spill, we want it set to bad
3487 
3488  // sumdata::SpillData variables
3489  spill.det = fDet;
3490  spill.ismc = !fIsRealData;
3491  spill.isgoodspill = false;
3492  spill.spilltimesec = 0;
3493  spill.spilltimensec = 0;
3494  spill.gpsspilltimesec = 0;
3495  spill.gpsspilltimensec = 0;
3496  spill.deltaspilltimensec = 0;
3497  spill.spillpot = 0.0;
3498  spill.livetime = 0.0;
3499  spill.hornI = 0.0;
3500  spill.posx = 0.0;
3501  spill.posy = 0.0;
3502  spill.widthx = 0.0;
3503  spill.widthy = 0.0;
3504 
3505  // Calculate diblock masking. Also in header, but needed here for
3506  // mass accounting.
3507  unsigned int maskstatus = 0;
3508  unsigned int dibmask = fMask;
3509 
3510  if (rh->GoodDiBlockMask() == ((1 << 16) - 1))
3511  maskstatus = 0;
3512  else if (rh->GoodDiBlockMask() == ((1 << 17) - 1))
3513  maskstatus = 2;
3514  else
3515  maskstatus = 1;
3516 
3517  unsigned int dibfirst = 0;
3518  unsigned int diblast = 0;
3519  if (dibmask) {
3520  int iD;
3521  iD = 0;
3522  while (!((dibmask >> iD) & 1)) iD++;
3523  dibfirst = iD + 1;
3524  iD = 0;
3525  while (dibmask >> iD) iD++;
3526  diblast = iD;
3527  }
3528 
3529  spill.dibfirst = dibfirst;
3530  spill.diblast = diblast;
3531  spill.dibmask = dibmask;
3532  spill.maskstatus = maskstatus;
3533 
3534  // sumdata::EventQuality vars, filled by DQSpillFlags
3535  spill.nmissingdcms = 0;
3536  spill.fracdcm3hits = 0.0;
3537  spill.nouttimehits = 0;
3538  spill.nnoisyapds = 0;
3539 
3540  spill.nmicroslices = 0;
3541  spill.ndcms = 0;
3542 
3543  spill.nmissingdcmslg = 0; // LiveGeometry check
3544  spill.nbaddcmslg = 0; // LiveGeometry check
3545  spill.dcmedgematchfrac = 0;
3546 
3547  // Need to flag each spill if the pot was good or bad
3548  // each slice will get the same tag
3549 
3550  // If there was no spill info for this slice mark it bad. log message this
3551  // should not happen
3552  if (spillPot.failedToGet()) {
3553  mf::LogError("CAFMaker") << "Spill Data not found for real data event";
3554  spill.isgoodspill = false;
3555  }
3556  // If we have a SpillData object, get all the variables
3557  else {
3558  spill.isgoodspill = spillPot->goodbeam;
3559  spill.spilltimesec = spillPot->spilltimesec;
3560  spill.spilltimensec = spillPot->spilltimensec;
3561  spill.gpsspilltimesec = spillPot->gpsspilltimensec;
3562  spill.gpsspilltimensec = spillPot->gpsspilltimensec;
3563  spill.deltaspilltimensec = spillPot->deltaspilltimensec;
3564  spill.spillpot = spillPot->spillpot;
3565  spill.hornI = spillPot->hornI;
3566  spill.isRHC = spillPot->isRHC;
3567  spill.is0HC = spillPot->is0HC;
3568  spill.isFHC = !spill.isRHC && !spill.is0HC;
3569  spill.widthx = spillPot->widthx;
3570  spill.widthy = spillPot->widthy;
3571  spill.posx = spillPot->posx;
3572  spill.posy = spillPot->posy;
3573 
3574  spill.intx.resize(spillPot->intx.size());
3575  spill.inty.resize(spillPot->inty.size());
3576  spill.bposx.resize(spillPot->bposx.size());
3577  spill.bposy.resize(spillPot->bposy.size());
3578  for (size_t i = 0; i < spillPot->intx.size(); ++i)
3579  spill.intx[i] = spillPot->intx[i];
3580  for (size_t i = 0; i < spillPot->inty.size(); ++i)
3581  spill.inty[i] = spillPot->inty[i];
3582  for (size_t i = 0; i < spillPot->bposx.size(); ++i)
3583  spill.bposx[i] = spillPot->bposx[i];
3584  for (size_t i = 0; i < spillPot->bposy.size(); ++i)
3585  spill.bposy[i] = spillPot->bposy[i];
3586 
3587  if (!fIsRealData) {
3588  // MC will flag each spill as good since it doesn't make bad spills
3589  spill.isgoodspill = true;
3590  } else {
3591  // NOTE: there is an implied factor of 1e12 when getting
3592  // IFDBSpillInfo POT. Adding this back in to avoid confusion
3593  // between data and MC.
3594  spill.spillpot *= 1e12;
3595  }
3596  }
3597 
3598  // If there was no spill info for this slice mark it bad.
3599  // Log message this should not happen
3600  if (spillQual.failedToGet()) {
3601  mf::LogError("CAFMaker")
3602  << "Spill EventQuality not found for real data event";
3603  } else {
3604  spill.nmissingdcmslg = spillQual->nmissingdcmslg;
3605  spill.nmissingdcms = spillQual->nmissingdcms;
3606  spill.fracdcm3hits = spillQual->fracdcm3hits;
3607  spill.nouttimehits = spillQual->nouttimehits;
3608  spill.nnoisyapds = spillQual->nnoisyapds;
3609  spill.dcmedgematchfrac = spillQual->dcmedgematchfrac;
3610  spill.nmicroslices = spillQual->nmicroslices;
3611  }
3612 
3613  if (trigs.failedToGet()) {
3614  mf::LogError("CAFMaker") << "No trigger info found for this event";
3615  } else {
3616  bool once = true;
3617  for (const rawdata::RawTrigger& trig : *trigs) {
3618  // Units of trigger field are 500ns
3619  spill.livetime += trig.fTriggerRange_TriggerLength * 500e-9;
3620  if (once) {
3621  once = false;
3622  spill.trigger = trig.fTriggerMask_TriggerType;
3623  } else {
3624  if (spill.trigger != trig.fTriggerMask_TriggerType) {
3625  mf::LogError("CAFMaker") << "Multiple trigger types found in one "
3626  "Event. Storing the first value seen.";
3627  }
3628  }
3629  }
3630  }
3631 
3632  // dq::DAQEventSummary variables
3633  if (daq.failedToGet()) {
3634  mf::LogError("CAFMaker") << "No DAQ summary info found for this event.";
3635  } else {
3636  spill.ndiblocks = daq->fNDiblocks;
3637  spill.eventincomplete = daq->fEventIncomplete;
3638  spill.emptydatablock = daq->fNemptyDataBlock;
3639  spill.nmicroblocks = daq->fNmicroBlocks;
3640  spill.nemptymicroslice = daq->fNemptyMicroSlice;
3641  spill.ndroppedmicroblocks = daq->fNDroppedMicroBlocks;
3642  spill.ndatablockmissingdata = daq->fNdataBlockMissingData;
3643  spill.nmicroslicedatanotpresent = daq->fNmicroSliceDataNotPresent;
3644  spill.nnanoslices = daq->fNtotalNanoSlices;
3645  spill.nanoslicedatanotpresent = daq->fNanoSliceDataNotPresent;
3646  spill.nanoslicenolinkstatus = daq->fNanoSliceNoLinkStatus;
3647  spill.nanoslicebufferempty = daq->fNanoSliceBufferEmpty;
3648  spill.nanoslicebufferfull = daq->fNanoSliceBufferFull;
3649  spill.nanoslicecommerror = daq->fNanoSliceCommError;
3650  spill.nanoslicepacketerror = daq->fNanoSlicePacketError;
3651  spill.nanosliceoverflowerror = daq->fNanoSliceOverflowError;
3652  spill.nanosliceadcerror = daq->fNanoSliceADCError;
3653  }
3654 
3655  spill.nbaddcmslg = livegeom->NBadDCMBC();
3656  spill.ndcms = rh->NDCMs();
3657 }
art::ServiceHandle< geo::LiveGeometry > livegeom
bool isRHC
is the beam in antineutrino mode, aka RHC
Definition: SpillData.h:28
SubRunNumber_t subRun() const
Definition: Event.h:72
bool is0HC
horn off
Definition: SpillData.h:29
int nouttimehits
of out-beam-window hits
Definition: EventQuality.h:25
std::vector< double > bposy
Definition: SpillData.h:33
int fNdataBlockMissingData
of occurances of isMissingData
int fNmicroBlocks
How many microblocks?
int fNanoSliceADCError
of nanoslices reporting ADCError
int fNanoSliceCommError
of nanoslices reporting CommError
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
int nmissingdcms
of missing DCMs
Definition: EventQuality.h:23
Det_t fDet
Detector ID in caf namespace typedef.
int nmicroslices
of micro slices
Definition: EventQuality.h:27
std::vector< double > bposx
Definition: SpillData.h:32
int fNDiblocks
of diblocks reporting in event
int fNemptyMicroSlice
How many empty micro slices?
std::vector< double > intx
Definition: SpillData.h:30
int fNanoSliceNoLinkStatus
of nanoslices reporting !LinkPresent
if(dump)
signed long long int deltaspilltimensec
Definition: SpillData.h:24
double fracdcm3hits
fraction of DCM3 hits in horizontal modules
Definition: EventQuality.h:24
float dcmedgematchfrac
Low values mean out-of-sync detector.
Definition: EventQuality.h:29
int fNanoSliceDataNotPresent
of nanoslices reporting !DataPresent
bool fEventIncomplete
Is the event incomplete?
double widthx
mm
Definition: SpillData.h:36
unsigned long int spilltimensec
Definition: SpillData.h:21
double hornI
kA
Definition: SpillData.h:27
int fNanoSliceOverflowError
of nanoslices reporting OverflowError
int fNanoSliceBufferFull
of nanoslices reporting BufferFull
double posy
mm
Definition: SpillData.h:35
int fNemptyDataBlock
How many empty data blocks?
int fNmicroSliceDataNotPresent
of microslices with !DataPresent
EventNumber_t event() const
Definition: EventID.h:116
int fNanoSlicePacketError
of nanoslices reporting PacketError
unsigned long int spilltimesec
Definition: SpillData.h:20
int fNtotalNanoSlices
of nano slices in the event
int fNanoSliceBufferEmpty
of nanoslices reporting BufferEmpty
int GoodDiBlockMask(int subrun=-1, bool reload=false)
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
#define for
Definition: msvc_pragmas.h:3
art::ServiceHandle< nova::dbi::RunHistoryService > rh
double posx
mm
Definition: SpillData.h:34
int nnoisyapds
of noisy APDs
Definition: EventQuality.h:26
unsigned long int gpsspilltimensec
Definition: SpillData.h:23
int fNDroppedMicroBlocks
How many dropped micro blocks?
double widthy
mm
Definition: SpillData.h:37
double spillpot
POT for spill normalized by 10^12.
Definition: SpillData.h:26
EventID id() const
Definition: Event.h:56
bool failedToGet() const
Definition: Handle.h:196
std::vector< double > inty
Definition: SpillData.h:31
template<class T , class U >
art::FindManyP< T > caf::CAFMaker::FindManyPStrict ( const U &  from,
const art::Event evt,
const art::InputTag label 
) const
protected

Equivalent of FindManyP except a return that is !isValid() prints a messsage and aborts if StrictMode is true.

Definition at line 1002 of file CAFMaker_module.cc.

References om::cout, allTimeWatchdog::endl, fParams, art::InputTag::label(), runNovaSAM::ret, caf::CAFMakerParams::StrictMode, and T.

1004  {
1005  art::FindManyP<T> ret(from, evt, tag);
1006 
1007  if (!tag.label().empty() && !ret.isValid() && fParams.StrictMode()) {
1008  std::cout << "CAFMaker: No Assn from '"
1009  << abi::__cxa_demangle(typeid(from).name(), 0, 0, 0) << "' to '"
1010  << abi::__cxa_demangle(typeid(T).name(), 0, 0, 0)
1011  << "' found under label '" << tag << "'. "
1012  << "Set 'StrictMode: false' to continue anyway." << std::endl;
1013  abort();
1014  }
1015 
1016  return ret;
1017 }
const XML_Char * name
Definition: expat.h:151
Atom< bool > StrictMode
OStream cout
Definition: OStream.cxx:6
CAFMakerParams fParams
double T
Definition: Xdiff_gwt.C:5
seed_t art::EngineCreator::get_seed_value ( fhicl::ParameterSet const &  pset,
char const  key[] = "seed",
seed_t const  implicit_seed = -1 
)
inherited
template<class T >
bool caf::CAFMaker::GetAssociatedProduct ( const art::FindManyP< T > &  fm,
int  idx,
T ret 
) const
protected

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.

Parameters
fmThe FindManyP object describing the association
idxWhich element of the FindManyP to look it
[out]retThe product retrieved
Returns
Whether ret was filled

Definition at line 1021 of file CAFMaker_module.cc.

Referenced by produce().

1022  {
1023  if (!fm.isValid()) return false;
1024 
1025  const std::vector<art::Ptr<T>> prods = fm.at(idx);
1026 
1027  if (prods.empty()) return false;
1028 
1029  ret = *prods[0];
1030 
1031  return true;
1032 }
template<class T >
void caf::CAFMaker::GetByLabelIfExists ( const art::Event evt,
const std::string label,
art::Handle< T > &  handle 
) const
protected

Equivalent of evt.getByLabel(label, handle) except failedToGet prints a message.

Definition at line 1050 of file CAFMaker_module.cc.

References om::cout, allTimeWatchdog::endl, art::Handle< T >::failedToGet(), fParams, art::DataViewImpl::getByLabel(), and caf::CAFMakerParams::StrictMode.

1052  {
1053  evt.getByLabel(label, handle);
1054  if (!label.empty() && handle.failedToGet() && fParams.StrictMode()) {
1055  std::cout << "CAFMaker: No product of type '"
1056  << abi::__cxa_demangle(typeid(*handle).name(), 0, 0, 0)
1057  << "' found under label '" << label << "'. "
1058  << "Continuing without it." << std::endl;
1059  }
1060 }
Atom< bool > StrictMode
const char * label
OStream cout
Definition: OStream.cxx:6
CAFMakerParams fParams
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
bool failedToGet() const
Definition: Handle.h:196
template<class T >
void caf::CAFMaker::GetByLabelStrict ( const art::Event evt,
const std::string label,
art::Handle< T > &  handle 
) const
protected

Equivalent of evt.getByLabel(label, handle) except failedToGet prints a message and aborts if StrictMode is true.

Definition at line 1036 of file CAFMaker_module.cc.

References om::cout, allTimeWatchdog::endl, art::Handle< T >::failedToGet(), fParams, art::DataViewImpl::getByLabel(), and caf::CAFMakerParams::StrictMode.

Referenced by AddMCTruthToVec(), and produce().

1037  {
1038  evt.getByLabel(label, handle);
1039  if (!label.empty() && handle.failedToGet() && fParams.StrictMode()) {
1040  std::cout << "CAFMaker: No product of type '"
1041  << abi::__cxa_demangle(typeid(*handle).name(), 0, 0, 0)
1042  << "' found under label '" << label << "'. "
1043  << "Set 'StrictMode: false' to continue anyway." << std::endl;
1044  abort();
1045  }
1046 }
Atom< bool > StrictMode
const char * label
OStream cout
Definition: OStream.cxx:6
CAFMakerParams fParams
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
bool failedToGet() const
Definition: Handle.h:196
template<typename PROD , BranchType B>
ProductID art::EDProducer::getProductID ( std::string const &  instanceName = {}) const
inlineinherited

Definition at line 123 of file EDProducer.h.

References art::EDProducer::moduleDescription_.

Referenced by skim::NueSkimmer::CopyMichelSlice(), and skim::NueSkimmer::CopyMichelTrack().

124  {
125  return ProducerBase::getProductID<PROD, B>(moduleDescription_,
126  instanceName);
127  }
ModuleDescription moduleDescription_
Definition: EDProducer.h:115
template<typename PROD , BranchType B>
ProductID art::ProducerBase::getProductID ( ModuleDescription const &  moduleDescription,
std::string const &  instanceName 
) const
inherited

Definition at line 56 of file ProducerBase.h.

References art::ModuleDescription::moduleLabel().

Referenced by art::ProducerBase::modifiesEvent().

58  {
59  auto const& pd =
60  get_ProductDescription<PROD>(B, md.moduleLabel(), instanceName);
61  return pd.productID();
62  }
template<class T >
bool caf::CAFMaker::GetPsetParameter ( const fhicl::ParameterSet pset,
const std::vector< std::string > &  name,
T ret 
) const
protected
Parameters
psetThe parameter set
namePass "foo.bar.baz" as {"foo", "bar", "baz"}
[out]retValue of the key, not set if we return false
Returns
Whether the key was found

Definition at line 1064 of file CAFMaker_module.cc.

References fhicl::ParameterSet::get(), fhicl::ParameterSet::has_key(), MECModelEnuComparisons::i, and T.

1066  {
1067  fhicl::ParameterSet p = pset;
1068  for (unsigned int i = 0; i < name.size() - 1; ++i) {
1069  if (!p.has_key(name[i])) return false;
1070  p = p.get<fhicl::ParameterSet>(name[i]);
1071  }
1072  if (!p.has_key(name.back())) return false;
1073  ret = p.get<T>(name.back());
1074  return true;
1075 }
const XML_Char * name
Definition: expat.h:151
const char * p
Definition: xmltok.h:285
T get(std::string const &key) const
Definition: ParameterSet.h:231
bool has_key(std::string const &key) const
double T
Definition: Xdiff_gwt.C:5
void caf::CAFMaker::InitializeOutfile ( )
protected

Definition at line 485 of file CAFMaker_module.cc.

References ana::assert(), genie::utils::res::AsString(), run_hadd::cmd, APDHVSetting::dummy, datagram_client::envmap, fBatch, fCafFilename, fclose(), fCycle, fFile, caf::CAFMakerParams::FillNuTree, fNuTree, fParams, fRecTree, fSpillTree, fTotalEvents, fTotalLivetime, fTotalPOT, fTotalSinglePOT, fTotalTrueNonswapNue, fTotalTrueSingleNue, cet::getenv(), hEvents, hLivetime, hPOT, hSinglePOT, hTotalTrueNonswapNue, hTotalTrueSingleNue, findDuplicateFiles::key, rec, gen_flatrecord::size, split(), string, PandAna.Demos.tute_pid_validation::var, and Write().

Referenced by beginJob(), and respondToOpenInputFile().

485  {
486  assert(!fFile);
487  assert(!fCafFilename.empty());
488 
489  mf::LogInfo("CAFMaker") << "Output filename is " << fCafFilename;
490 
491  fFile = new TFile(fCafFilename.c_str(), "RECREATE");
492 
493  hPOT = new TH1D("TotalPOT", "TotalPOT;; POT", 1, 0, 1);
494  hSinglePOT =
495  new TH1D("TotalSinglePOT", "TotalSinglePOT;; Single POT", 1, 0, 1);
497  new TH1D("TotalTrueNonswapNue",
498  "TotalTrueNonswap;; Total True Nonswap Nue", 1, 0, 1);
499  hTotalTrueSingleNue = new TH1D(
500  "TotalTrueSingleNue", "TotalTrueSingle;; Total True Single Nue", 1, 0, 1);
501  hEvents = new TH1D("TotalEvents", "TotalEvents;; Events", 1, 0, 1);
502  hLivetime = new TH1D("TotalLivetime", "TotalLiveTime;;Seconds", 1, 0, 1);
503 
504  fRecTree = new TTree("recTree", "records");
505 
506  // Tell the tree it's expecting StandardRecord objects
507  StandardRecord* rec = 0;
508  fRecTree->Branch("rec", "caf::StandardRecord", &rec);
509 
510  // Create tree to store spill info for POT counting
511  fSpillTree = new TTree("spillTree", "spills");
512 
513  caf::SRSpill* dummy = 0;
514  fSpillTree->Branch("spill", "caf::SRSpill", &dummy);
515 
516  // Make the spillTree a friend of the fRecTree
517  fRecTree->AddFriend("spillTree");
518 
519 
520  if (fParams.FillNuTree()) {
521  fNuTree = new TTree("nuTree", "spillTruth");
522  caf::SRSpillTruthBranch* dummy2 = 0;
523  fNuTree->Branch("mc", "caf::SRSpillTruthBranch", &dummy2);
524  } else {
525  fNuTree = 0;
526  }
527 
528  fTotalPOT = 0;
529  fTotalSinglePOT = 0;
532  fTotalEvents = 0;
533  fTotalLivetime = 0;
534  fCycle = -5;
535  fBatch = -5;
536 
537  // Global information about the processing details:
538 
539  std::map<TString, TString> envmap;
540  std::string envstr;
541  // Environ comes from unistd.h
542  // environ is not present on OSX for some reason, so just use getenv to
543  // grab the variables we care about.
544 #ifdef DARWINBUILD
545  std::set<TString> variables;
546  variables.insert("USER");
547  variables.insert("HOSTNAME");
548  variables.insert("PWD");
549  variables.insert("SRT_PUBLIC_CONTEXT");
550  variables.insert("SRT_PRIVATE_CONTEXT");
551  for (auto var : variables) envmap[var] = getenv(var);
552 #else
553  for (char** penv = environ; *penv; ++penv) {
554  const std::string pair = *penv;
555  envstr += pair;
556  envstr += "\n";
557  const size_t split = pair.find("=");
558  if (split == std::string::npos) continue; // Huh?
559  const std::string key = pair.substr(0, split);
560  const std::string value = pair.substr(split + 1);
561  envmap[key] = value;
562  }
563 #endif
564 
565  // Get the command-line we were invoked with. What I'd really like is
566  // just the fcl script and list of input filenames in a less hacky
567  // fashion. I'm not sure that's actually possible in ART.
568  // TODO: ask the artists.
569  FILE* cmdline = fopen("/proc/self/cmdline", "rb");
570  char* arg = 0;
571  size_t size = 0;
573  while (getdelim(&arg, &size, 0, cmdline) != -1) {
574  cmd += arg;
575  cmd += " ";
576  }
577  free(arg);
578  fclose(cmdline);
579 
580  fFile->mkdir("env")->cd();
581 
582  TObjString(envmap["USER"]).Write("user");
583  TObjString(envmap["HOSTNAME"]).Write("hostname");
584  TObjString(envmap["PWD"]).Write("pwd");
585  TObjString(envmap["SRT_PUBLIC_CONTEXT"]).Write("publiccontext");
586  TObjString(envmap["SRT_PRIVATE_CONTEXT"]).Write("privatecontext");
587  // Default constructor is "now"
588  TObjString(TTimeStamp().AsString()).Write("date");
589  TObjString(cmd.c_str()).Write("cmd");
590  TObjString(fCafFilename.c_str()).Write("output");
591  TObjString(envstr.c_str()).Write("env");
592 }
void split(double tt, double *fr)
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Truth info for all neutrinos in the spill.
Atom< bool > FillNuTree
caf::StandardRecord * rec
Definition: tutCAFMacro.C:20
std::string fCafFilename
TH1D * hTotalTrueSingleNue
fclose(fg1)
TH1D * hTotalTrueNonswapNue
string cmd
Definition: run_hadd.py:52
const XML_Char int const XML_Char * value
Definition: expat.h:331
std::string getenv(std::string const &name)
CAFMakerParams fParams
assert(nhit_max >=nhit_nbins)
const char * AsString(Resonance_t res)
resonance id -> string
gm Write()
enum BeamMode string
template<typename T , BranchType = InEvent>
ProductToken<T> art::Consumer::mayConsume ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ProductToken<T> art::Consumer::mayConsume ( InputTag const &  it)
inherited

Definition at line 189 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

190 {
191  if (!moduleContext_)
192  return ProductToken<T>::invalid();
193 
194  consumables_[BT].emplace_back(ConsumableType::Product,
195  TypeID{typeid(T)},
196  it.label(),
197  it.instance(),
198  it.process());
199  return ProductToken<T>{it};
200 }
set< int >::iterator it
static ProductToken< T > invalid()
Definition: ProductToken.h:47
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename T , art::BranchType BT>
void art::Consumer::mayConsumeMany ( )
inherited

Definition at line 204 of file Consumer.h.

References T.

205 {
206  if (!moduleContext_)
207  return;
208 
209  consumables_[BT].emplace_back(ConsumableType::Many, TypeID{typeid(T)});
210 }
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
template<typename Element , BranchType = InEvent>
ViewToken<Element> art::Consumer::mayConsumeView ( InputTag const &  )
inherited
template<typename T , art::BranchType BT>
art::ViewToken<T> art::Consumer::mayConsumeView ( InputTag const &  it)
inherited

Definition at line 214 of file Consumer.h.

References art::InputTag::instance(), PandAna.reco_validation.prod5_pid_validation::invalid, art::InputTag::label(), art::InputTag::process(), and T.

215 {
216  if (!moduleContext_)
217  return ViewToken<T>::invalid();
218 
219  consumables_[BT].emplace_back(ConsumableType::ViewElement,
220  TypeID{typeid(T)},
221  it.label(),
222  it.instance(),
223  it.process());
224  return ViewToken<T>{it};
225 }
set< int >::iterator it
static ViewToken< Element > invalid()
Definition: ProductToken.h:75
ConsumableProducts consumables_
Definition: Consumer.h:137
double T
Definition: Xdiff_gwt.C:5
bool moduleContext_
Definition: Consumer.h:135
bool art::ProducerBase::modifiesEvent ( ) const
inlineinherited

Definition at line 40 of file ProducerBase.h.

References art::ProducerBase::getProductID(), and string.

41  {
42  return true;
43  }
static cet::exempt_ptr<Consumer> art::Consumer::non_module_context ( )
staticinherited
void art::Consumer::prepareForJob ( fhicl::ParameterSet const &  pset)
protectedinherited
void caf::CAFMaker::produce ( art::Event evt)
virtualnoexcept

NC Cosmic Rejection

NC Pi0 Bkg Rejection

Implements art::EDProducer.

Definition at line 1103 of file CAFMaker_module.cc.

References abs(), caf::AddCosmicTruthToVec(), AddMCTruthToVec(), caf::AddSlcMEToVec(), caf::AddTrkMEToVec(), rb::Cluster::AllCells(), caf::SRTruthBranch::allcosmics, cheat::BackTracker::allMCTruth(), caf::SRTruthBranch::allnus, caf::SRELid::ann, caf::SRELid::anne, caf::SRELid::annecos, caf::SRELid::annepi0, caf::CAFMakerParams::ApplyingFilter, caf::CAFMakerParams::ApplyLEMNuePresel, caf::SRHeader::batch, bc, caf::SRHeader::blind, caf::BlindThisRecord(), caf::SRFuzzyKProng::bpf, caf::CAFMakerParams::BreakPointTrackLabel, bt, calcFirstLastLivePlane(), rb::Cluster::CalorimetricEnergy(), simb::MCNeutrino::CCNC(), stan::math::ceil(), rb::CellHit::Cell(), geo::PlaneGeo::Cell(), rb::Cluster::Cell(), om::cerr, febshutoff_auto::chan, caf::SRSlice::closestslicecalE, caf::SRSlice::closestslicemindist, caf::SRSlice::closestsliceminfromback, caf::SRSlice::closestsliceminfrombottom, caf::SRSlice::closestsliceminfromeast, caf::SRSlice::closestsliceminfromfront, caf::SRSlice::closestsliceminfromtop, caf::SRSlice::closestsliceminfromwest, caf::SRSlice::closestslicenhit, caf::SRSlice::closestslicetime, caf::CAFMakerParams::ClusterLabel, compMuonID(), caf::SRIDBranch::contain, caf::CopyMuonIDVars(), caf::CopyRemidVars(), calib::corr(), caf::SRContain::cosbakcell, cosrej::CosRejObj::CosBakCell(), caf::SRContain::cosfwdcell, cosrej::CosRejObj::CosFwdCell(), caf::SRTruthBranch::cosmic, caf::SRTrackBranch::cosmic, caf::SRSpill::cosmiccvn, caf::CAFMakerParams::CosmicCVNLabel, caf::CAFMakerParams::CosmicTrackLabel, cvntf::CVNCosmicFilt::cosmicVal, caf::SRCosmicCVN::cosmicVal, caf::SRIDBranch::cosrej, caf::CAFMakerParams::CosRejBPFLabel, caf::CAFMakerParams::CosRejLabel, om::cout, util::CreateAssn(), caf::SRIDBranch::cvn, caf::CAFMakerParams::CVNLabelLoosePreselPtP, caf::CAFMakerParams::CVNLabelLoosePreselPtPOppHorn, caf::CAFMakerParams::CVNLabelNeutronProng, caf::CAFMakerParams::CVNLabelNoCosmics, caf::CAFMakerParams::CVNLabelNoCosmicsOppHorn, caf::CAFMakerParams::CVNLabelOldPresel, caf::CAFMakerParams::CVNLabelOldPreselOppHorn, caf::SRMRProperties::cvnloosepreselptp, caf::SRIDBranch::cvnloosepreselptp, caf::SRIDBranch::cvnloosepreselptp_opphorn, caf::SRTrainingBranch::cvnmaps, caf::SRTrainingBranch::cvnmapsfilled, caf::SRMRProperties::cvnnocosmics, caf::SRIDBranch::cvnnocosmics, caf::SRIDBranch::cvnnocosmics_opphorn, caf::SRMRProperties::cvnoldpresel, caf::SRIDBranch::cvnoldpresel, caf::SRIDBranch::cvnoldpresel_opphorn, caf::CAFMakerParams::CVNParticle2DLabel, caf::CAFMakerParams::CVNParticleLabel, caf::CAFMakerParams::CVNParticleOppHornLabel, caf::CAFMakerParams::CVNPixelMapLabel, caf::CAFMakerParams::CVNTrainingDataLabel, caf::SRHeader::cycle, caf::CAFMakerParams::DAQEventSummaryLabel, caf::CAFMakerParams::DataSpillHasMC, caf::SRHeader::day, caf::SRSLidEnergy::depE, caf::SRHeader::det, nova::dbi::RunHistory::DetFineTimingSetting(), nova::dbi::RunHistory::DetGainSetting(), caf::SRHeader::dibfirst, caf::SRHeader::diblast, caf::SRHeader::dibmask, release_diff::diff, caf::CAFMakerParams::DiFShowerClusterLabel, caf::CAFMakerParams::DiFShowerCVNLabelLoosePreselPtP, caf::CAFMakerParams::DiFShowerCVNLabelNoCosmics, caf::CAFMakerParams::DiFShowerCVNLabelOldPresel, caf::CAFMakerParams::DiFShowerLIDLabel, caf::CAFMakerParams::DiFShowerShwLabel, caf::SRTrackBranch::discrete, caf::CAFMakerParams::DiscreteTrack2dLabel, caf::CAFMakerParams::DiscreteTrackLabel, caf::SRNueCosRej::distallpngback, caf::SRNueCosRej::distallpngbottom, caf::SRNueCosRej::distallpngeast, caf::SRNueCosRej::distallpngfront, caf::SRNueCosRej::distallpngtop, caf::SRNueCosRej::distallpngwest, geo::LiveGeometry::DistToBack(), geo::LiveGeometry::DistToBottom(), geo::LiveGeometry::DistToEdgeXMinus(), geo::LiveGeometry::DistToEdgeXPlus(), geo::LiveGeometry::DistToFront(), geo::LiveGeometry::DistToTop(), caf::SRHeader::doy, APDHVSetting::dummy, caf::SRSLidEnergy::E, e, cheat::EffMetric(), cheat::EffPurMetric(), caf::SRVertexBranch::elastic, caf::CAFMakerParams::ElasticArmsLabel, caf::CAFMakerParams::EnableBlindness, allTimeWatchdog::endl, caf::StandardRecord::energy, cheat::EnergyMetric(), lem_server::ep, evt, caf::SRHeader::evt, stan::math::fabs(), caf::SRTruthBranch::faveidxeff, caf::SRTruthBranch::faveidxeffpur, caf::SRTruthBranch::faveidxeffthenpur, caf::SRTruthBranch::faveidxenergy, caf::SRTruthBranch::faveidxpur, fBatch, fCycle, fDet, caf::CAFMakerParams::FillCosmicCVN, caf::FillCosRejVars(), caf::FillCVNNeutronDaughterResultVars(), caf::FillCVNParticleResultVars(), caf::FillCVNPixelMaps(), caf::FillCVNProngTrainingData(), caf::FillCVNResultVars(), caf::FillCVNTrainingData(), caf::FillDiFShowerVars(), caf::FillDiFVars(), caf::FillLEMVars(), caf::FillLIDEnergyVars(), caf::FillMRCCParentInfo(), caf::FillMuIdVars(), caf::FillMuonIDVars(), caf::FillNCCosRejVars(), caf::FillNCPi0BkgRejVars(), caf::FillNDRecoBPFVars(), caf::FillNDRecoPngVars(), caf::FillNDRecoTrkVars(), caf::FillNueCosRejVars(), caf::FillNueEnergyVars(), caf::FillNuePreselVars(), caf::FillNumuEnergyVars(), caf::FillNumuLSTMEnergyVars(), caf::FillNuonEResultVars(), caf::CAFMakerParams::FillNuTree, caf::CAFMakerParams::FillPixelMaps, caf::FillProngVars(), caf::FillRockPreselVars(), caf::FillRVPVars(), caf::FillShowerVars(), caf::SRTrackBase::fillSizes(), caf::SRFuzzyK::fillSizes(), caf::SRMichelE::fillSizes(), caf::SRVertexBranch::fillSizes(), caf::SRKalman::fillSizes(), caf::FillSliceLID(), caf::FillSliceVars(), caf::FillSlidVars(), FillSpillVars(), caf::FillTrackContainmentVars(), caf::FillTrackInfoVars(), caf::FillTrackVars(), caf::FillTrackVarsBpfFitSum(), caf::FillTrackVarsWithOverlapE(), caf::CAFMakerParams::FillTrainingData, caf::FillVetoVars(), caf::FillXnueVars(), caf::SRHeader::filt, caf::SRHeader::finetiming, caf::SRSlice::firstplane, fIsRealData, check_time_usage::float, caf::CAFMakerParams::FluxReweightLabel, caf::FluxReweights(), fMask, fNuTree, cvn::Result::fOutput, fParams, fRecTree, fSpillTree, fTotalEvents, fTotalSinglePOT, fTotalTrueNonswapNue, fTotalTrueSingleNue, caf::CAFMakerParams::FullTruthInfo, caf::SRElastic::fuzzyk, caf::CAFMakerParams::FuzzyK2DLabel, caf::CAFMakerParams::FuzzyK3DLabel, caf::CAFMakerParams::FuzzyKOrphLabel, caf::CAFMakerParams::FuzzyKWeight2DLabel, caf::CAFMakerParams::FuzzyKWeight3DLabel, caf::CAFMakerParams::G4ReweightLabel, caf::SRHeader::gain, caf::SRMCReweight::geant4, caf::Geant4Reweights(), caf::CAFMakerParams::GeneratorLabel, caf::SRMCReweight::genie, caf::GenieReweightTable(), geom, GetAssociatedProduct(), GetByLabelStrict(), meta::MetadataManager::getInstance(), meta::MetadataManager::GetMetadata(), simb::MCTruth::GetNeutrino(), rb::RecoHit::GeV(), caf::SRTruthBranch::global, nova::dbi::RunHistory::GoodDiBlockMask(), caf::SRSLidEnergy::hadE, geo::CellGeo::HalfD(), geo::CellGeo::HalfW(), cheat::BackTracker::HaveTruthInfo(), caf::StandardRecord::hdr, muonid::HighestPIDTrack(), remid::HighestPIDTrack(), caf::CAFMakerParams::HitsLabel, cheat::BackTracker::HitsToTrackIDE(), caf::SRVertexBranch::hough, caf::CAFMakerParams::HoughVertexLabel, caf::SRHeader::hour, MECModelEnuComparisons::i, compare_h5_caf::idx, caf::SRKalman::idxlongest, caf::SRKalman::idxmuonid, caf::SRKalman::idxremid, if(), makeTrainCVSamples::int, ip, rb::RecoHit::IsCalibrated(), rb::IsFiltered(), caf::SRHeader::ismc, std::isnan(), rb::Cluster::IsNoise(), caf::CAFMakerParams::IsSinglesOverlay, caf::SRBpfTrack::IsValid, caf::SRElastic::IsValid, art::Handle< T >::isValid(), caf::CAFMakerParams::JMShowerLabel, caf::SRContain::kalbakcell, cosrej::CosRejObj::KalBakCell(), caf::SRContain::kalbakcellnd, cosrej::CosRejObj::KalBakCellND(), caf::SRContain::kalfwdcell, cosrej::CosRejObj::KalFwdCell(), caf::SRContain::kalfwdcellnd, cosrej::CosRejObj::KalFwdCellND(), caf::SRTrackBranch::kalman, caf::CAFMakerParams::KalmanTrackLabel, caf::SRContain::kalyposattrans, ana::kDistAllBack, ana::kDistAllBottom, ana::kDistAllEast, ana::kDistAllFront, ana::kDistAllTop, ana::kDistAllWest, caf::kFARDET, caf::kNEARDET, geo::kXorY, caf::SRSlice::lastplane, caf::SRIDBranch::lem, caf::CAFMakerParams::LemLabel, caf::CAFMakerParams::LEMNuePresel, caf::SRNueEnergy::lid, caf::SRIDBranch::lid, caf::SRMRProperties::lid, caf::CAFMakerParams::LIDLabel, livegeom, LOG_DEBUG, caf::SRFuzzyK::longestidx, caf::SRNumuEnergy::lstmmuon, caf::SRNumuEnergy::lstmmuon_opphorn, caf::SRNumuEnergy::lstmnu, caf::SRNumuEnergy::lstmnu_opphorn, caf::SRHeader::maskstatus, std::max(), rb::Cluster::MaxTNS(), rb::Cluster::MaxXYZ(), caf::StandardRecord::mc, caf::CAFMakerParams::MCSpillDataLabel, caf::StandardRecord::me, rb::Cluster::MeanTNS(), caf::CAFMakerParams::MELabel, runNovaSAM::metadata, std::min(), rb::Cluster::MinTNS(), caf::SRHeader::minute, rb::Cluster::MinXYZ(), caf::SRHeader::month, caf::SRParentBranch::mrccpar, caf::CAFMakerParams::MRCCParentLabel, caf::CAFMakerParams::MRMode, caf::CAFMakerParams::MRParentSliceLabel, caf::SRBpf::muon, caf::SRIDBranch::muonid, caf::CAFMakerParams::MuonIDLabel, caf::SRTruthBranch::nallcosmics, caf::SRTruthBranch::nallnus, nova::dbi::RunHistory::NAnalysisChannels(), caf::SRHeader::nbadchan, chaninfo::BadChanList::NBadInSubRun(), caf::SRIDBranch::nccosrej, caf::CAFMakerParams::NCCosRejLabel, rb::Cluster::NCell(), caf::SRSlice::ncellsfromedge, caf::SRTruthBranch::ncosmic, caf::SRSpill::ncosmiccvn, caf::SRIDBranch::ncpi0bkgrej, caf::CAFMakerParams::NCPi0BkgRejLabel, cvntf::CVNCosmicFilt::ncVal, caf::SRCosmicCVN::ncVal, caf::SRNumuEnergy::ndhadcalcatE, caf::SRNumuEnergy::ndhadcaltranE, caf::CAFMakerParams::NDRecoLabel, cheat::NeutrinoEffPur::neutrinoInt, simb::MCTruth::NeutrinoSet(), cvntf::CVNCosmicFilt::nHits, caf::SRCosmicCVN::nHits, caf::SRSlice::nnonnoise, caf::SRGlobalTruth::nnu, caf::SRTruthBranch::nnu, caf::SRContain::nplanestoback, caf::SRContain::nplanestofront, caf::SRFuzzyK::nshwlid, caf::SRHeader::ntotchan, caf::SRKalman::ntracks, caf::SRSpillTruthBranch::nu, caf::SRTruthBranch::nu, simb::MCNeutrino::Nu(), caf::SREnergyBranch::nue, caf::SRIDBranch::nuecosrej, caf::CAFMakerParams::NueCosRejLabel, caf::SRIDBranch::nuepre, caf::CAFMakerParams::NuePreselLabel, caf::CAFMakerParams::NueSelLabel, cvntf::CVNCosmicFilt::nueVal, caf::SRCosmicCVN::nueVal, caf::CAFMakerParams::NueVetoLabel, caf::CAFMakerParams::NuMIBeamLabel, caf::SREnergyBranch::numu, caf::SRContain::numucontain, caf::SRContain::numucontainSA, caf::CAFMakerParams::NumuEnergyLabel, caf::CAFMakerParams::NumuLSTMEnergyLabel, caf::CAFMakerParams::NumuLSTMEnergyOppHornLabel, cvntf::CVNCosmicFilt::numuVal, caf::SRCosmicCVN::numuVal, caf::SRIDBranch::nuone, caf::SRIDBranch::nuone_opphorn, caf::CAFMakerParams::NuonELabel, caf::CAFMakerParams::NuonEOppHornLabel, cvntf::CVNCosmicFilt::nutauVal, caf::SRCosmicCVN::nutauVal, rb::Cluster::NXCell(), rb::Cluster::NYCell(), caf::SRFuzzyK::orphCalE, cvn::RegResult::Output(), caf::CAFMakerParams::OverlapEBPFLabel, caf::CAFMakerParams::OverlapEKalLabel, caf::StandardRecord::parent, file_size_ana::parent, murem::MRCCParent::ParentIndex(), cvntf::CVNCosmicFilt::passSel, caf::SRCosmicCVN::passSel, simb::MCParticle::PdgCode(), rb::CellHit::PE(), rb::RecoHit::PECorr(), caf::SRBpf::pion, std_candles::pl, rb::CellHit::Plane(), geo::GeometryBase::Plane(), caf::SRFuzzyK::png, caf::SRFuzzyK::png2d, caf::SRMCReweight::ppfx, caf::SRElastic::prong2dvertexenergyvolume10, caf::SRElastic::prong2dvertexenergyvolume20, caf::SRElastic::prong2dvertexenergyvolume30, caf::SRElastic::prong2dvertexenergyvolume40, caf::SRElastic::prong3dvertexenergyvolume10, caf::SRElastic::prong3dvertexenergyvolume20, caf::SRElastic::prong3dvertexenergyvolume30, caf::SRElastic::prong3dvertexenergyvolume40, caf::SRBpf::proton, cheat::PurMetric(), art::PtrVector< T >::push_back(), caf::CAFMakerParams::RawDataLabel, rec, caf::CAFMakerParams::ReclusShowerLabel, rb::Cluster::RecoHit(), caf::SRNueEnergy::regcvnEvtE, caf::SRNueEnergy::regcvnEvtE_opphorn, caf::SRNumuEnergy::regcvnhadE, caf::SRNumuEnergy::regcvnhadE_opphorn, caf::CAFMakerParams::RegCVNLabel, caf::CAFMakerParams::RegCVNOppHornLabel, caf::SRIDBranch::remid, caf::CAFMakerParams::RemidLabel, caf::CAFMakerParams::ReweightLabel, rh, caf::SRELid::rnn, caf::SRIDBranch::rockpre, caf::CAFMakerParams::RockPreselLabel, caf::SRHeader::run, updateRunHistoryTables::run, caf::SRIDBranch::rvp, caf::CAFMakerParams::RvpLabel, caf::SRNeutrino::rwgt, caf::SRHeader::second, caf::StandardRecord::sel, caf::SRNCPi0BkgRej::setDefault(), caf::SRNCCosRej::setDefault(), caf::SRNueCosRej::setDefault(), caf::SRSliceLID::setDefault(), caf::SRPresel::setDefault(), caf::SRFluxWeights::setDefault(), caf::SRMuonID::setDefault(), caf::SRSLidEnergy::setDefault(), caf::SRGeant4Weights::setDefault(), caf::SRXnue::setDefault(), caf::SRNueEnergy::setDefault(), caf::SRNuonEResult::setDefault(), caf::SRCVNResult::setDefault(), caf::SRELid::setDefault(), caf::SRRemid::setDefault(), caf::SRRvp::setDefault(), caf::SRVeto::setDefault(), caf::SRTruthBranch::setDefault(), caf::SRNumuEnergy::setDefault(), caf::SRMRCCParent::setDefault(), caf::SRLem::setDefault(), caf::SRNumuEnergy::setLSTMDefault(), caf::SRSLidEnergy::shwE, caf::CAFMakerParams::SingleMixerLabel, caf::SRMichelE::slc, caf::StandardRecord::slc, caf::CAFMakerParams::SlcLIDLabel, caf::SRIDBranch::slicelid, caf::SRIDBranch::slicelid_opphorn, caf::CAFMakerParams::SliceLIDLabel, caf::CAFMakerParams::SliceLIDOppHornLabel, caf::SRSliceMap::slicemap, caf::SRTrainingBranch::slicemaps, cheat::BackTracker::SlicesToMCTruthsTable(), cheat::BackTracker::SliceToOrderedNuIds(), caf::SRElastic::slicevertexenergyvolume10, caf::SRElastic::slicevertexenergyvolume20, caf::SRElastic::slicevertexenergyvolume30, caf::SRElastic::slicevertexenergyvolume40, sortProngCalE(), sortRBProngLength(), sortRBTrackLength(), sortRemid(), caf::CAFMakerParams::SPIDLabel, caf::StandardRecord::spill, caf_analysis::spill, caf::CAFMakerParams::SpillQualLabel, caf::CAFMakerParams::SPProngCVN5labelLabel, caf::CAFMakerParams::SPProngCVNNumuCCEMIDLabel, std::sqrt(), febshutoff_auto::start, caf::SRHeader::subevt, caf::SRHeader::subevtendtime, caf::SRHeader::subevtmeantime, caf::SRHeader::subevtstarttime, caf::SRHeader::subrun, getGoodRuns4SAM::subrun, caf::SRVertexDT::time, caf::SRElastic::time, caf::SRHoughVertex::time, art::Timestamp::timeHigh(), cvntf::CVNCosmicFilt::timeWinMax, caf::SRCosmicCVN::timeWinMax, cvntf::CVNCosmicFilt::timeWinMin, caf::SRCosmicCVN::timeWinMin, caf::CAFMakerParams::TrackInfoBPF, caf::CAFMakerParams::TrackInfoKalman, caf::SRTrackBase::tracks, caf::SRKalman::tracks, caf::SRKalman::tracks2d, caf::StandardRecord::training, caf::SRTrainingBranch::trainingdata, caf::StandardRecord::trk, caf::SRMichelE::trkcosmic, caf::SRMichelE::trkkalman, caf::SRHeader::unixtime, caf::CAFMakerParams::UseEnergy, caf::CAFMakerParams::UseGeVPixelMaps, caf::SRVertexBranch::vdt, caf::CAFMakerParams::VertexDTLabel, caf::SRIDBranch::veto, caf::CAFMakerParams::VetoLabel, rb::CellHit::View(), caf::SRVertexDT::vtx, caf::SRHoughVertex::vtx, caf::SRElastic::vtx, caf::StandardRecord::vtx, caf::SRTrackBranch::window, caf::CAFMakerParams::WindowTrackLabel, caf::CAFMakerParams::WSCNNLabel, caf::SRIDBranch::wsid, rb::Cluster::XCell(), caf::CAFMakerParams::XnueLabel, caf::SRIDBranch::xnuepid, rb::Cluster::YCell(), and caf::SRHeader::year.

1103  {
1104  // Normally CAFMaker is run without an output ART stream, so these go
1105  // nowhere, but can be occasionally useful for filtering as part of
1106  // an ART job.
1107  std::unique_ptr<std::vector<caf::StandardRecord>> srcol(
1108  new std::vector<caf::StandardRecord>);
1109  std::unique_ptr<art::Assns<caf::StandardRecord, rb::Cluster>> srAssn(
1111  std::unique_ptr<caf::SRSpill> spillcol(new SRSpill);
1112  std::unique_ptr<std::vector<caf::SRSpillTruthBranch>> spilltruthcol(
1113  new std::vector<SRSpillTruthBranch>);
1114 
1115  fTotalEvents += 1;
1116  // Adjust for art misfeature isRealData confusing data-mc overlay
1118  // -- Case -- : Data? == isRealData || DataSpillHasMC :
1119  // bt->HaveTruthInfo
1120  // Normal Data : true true false : false
1121  // Data Overlaid on MC : true false true : true
1122  // Just MC : false false false : true
1123 
1124  // Get a handle on the slices
1125 
1127  GetByLabelStrict(evt, fParams.ClusterLabel(), slices);
1128 
1129  /* // Define a vector holding slices
1130  art::PtrVector<rb::Cluster> subevts;
1131  for(unsigned int i = 0; i < slices->size(); ++i)
1132  subevts.push_back(art::Ptr<rb::Cluster>(slices, i));
1133  art::Handle<std::vector<rb::Cluster> > hparentSlices;
1134  art::PtrVector<rb::Cluster> parentSubevts;*/
1135 
1136  if (fParams.IsSinglesOverlay()) {
1138  art::Handle<std::vector<simb::MCTruth>> singlemctruthcol;
1139  GetByLabelStrict(evt, fParams.GeneratorLabel(), mctruthcol);
1140  GetByLabelStrict(evt, fParams.SingleMixerLabel(), singlemctruthcol);
1141 
1142  if (mctruthcol->empty()) {
1143  mf::LogWarning("CAFMaker") << "Error retrieving MCTruth list";
1144  return;
1145  }
1146  if (singlemctruthcol->empty()) {
1147  mf::LogWarning("CAFMaker")
1148  << "Error retrieving MCTruth list from singles";
1149  return;
1150  }
1151 
1152  for (size_t i_intx = 0; i_intx < mctruthcol->size(); ++i_intx) {
1153  LOG_DEBUG("outnNueCCOverlay") << "start loop";
1154  simb::MCTruth const& truth = mctruthcol->at(i_intx);
1155  simb::MCNeutrino const& mc_neutrino = truth.GetNeutrino();
1156  simb::MCParticle const& mc_particle = mc_neutrino.Nu();
1157 
1158  if (mc_neutrino.CCNC() == 0 && fabs(mc_particle.PdgCode()) == 12) {
1160  }
1161  LOG_DEBUG("CAFMaker") << "end loop over mctruth";
1162  }
1163 
1164  for (size_t i_intx = 0; i_intx < singlemctruthcol->size(); ++i_intx) {
1165  LOG_DEBUG("outnNueCCOverlay") << "start loop";
1166  simb::MCTruth const& truth = singlemctruthcol->at(i_intx);
1167  simb::MCNeutrino const& mc_neutrino = truth.GetNeutrino();
1168  simb::MCParticle const& mc_particle = mc_neutrino.Nu();
1169 
1170  if (mc_neutrino.CCNC() == 0 && fabs(mc_particle.PdgCode()) == 12) {
1172  }
1173  LOG_DEBUG("CAFMaker") << "end loop over mctruth";
1174  }
1175  }
1176 
1178 
1179  // in case we are running over MRE, make a vector of
1180  // parent slices, so we can pass it to backtracker
1181  // for filling the truth branch
1182 
1183  if (fParams.MRMode()) {
1184  GetByLabelStrict(evt, fParams.MRParentSliceLabel(), parentSlices);
1185  }
1186 
1187  // Get all of the hits from calhits
1188 
1190  GetByLabelStrict(evt, fParams.HitsLabel(), allHits);
1191 
1192  // We're gonna get a vector of NeutrinoEffPurs for allHits
1193  std::vector<cheat::NeutrinoWithIndex> allNus = bt->allMCTruth();
1194  std::vector<std::vector<cheat::NeutrinoEffPur>> sliceTruthTable;
1195 
1196 
1197  // Try to get POT spill info if this is real data
1199 
1200  // treat as data spill in MC-dataspill overlay
1201  if (!fIsRealData) {
1202  GetByLabelStrict(evt, fParams.MCSpillDataLabel(), spillPot);
1203  } else {
1204  GetByLabelStrict(evt, fParams.NuMIBeamLabel(), spillPot);
1205  }
1206 
1207  // Get the EventQuality information
1209  GetByLabelStrict(evt, fParams.SpillQualLabel(), spillQual);
1210 
1211  // Get DAQ Header information
1214 
1216  GetByLabelStrict(evt, fParams.RawDataLabel(), trigs);
1217 
1219  FillSpillVars(spillPot, spillQual, daq, trigs, spill, evt);
1220 
1221  // Fill the spill tree once per event so that POT counting can be
1222  // done properly.
1223  // This tree will have an entry even for empty events with no physics
1224  // slices.
1225  caf::SRSpill* pspill = &spill;
1226  fSpillTree->SetBranchAddress("spill", &pspill);
1227 
1229  if (fParams.FillCosmicCVN()){
1230  GetByLabelStrict(evt, fParams.CosmicCVNLabel(), cosmiccvn);
1231  std::vector<cvntf::CVNCosmicFiltList> cosmicList = *cosmiccvn;
1232  for (unsigned int iC = 0; iC < cosmicList[0].ListSize(); ++iC){
1233  cvntf::CVNCosmicFilt cvn = cosmicList[0].GetTimeSlice(iC);
1234  caf::SRCosmicCVN tcosmic;
1235  tcosmic.nHits = cvn.nHits;
1236  tcosmic.timeWinMax = cvn.timeWinMax;
1237  tcosmic.timeWinMin = cvn.timeWinMin;
1238  tcosmic.numuVal = cvn.numuVal;
1239  tcosmic.nueVal = cvn.nueVal;
1240  tcosmic.nutauVal = cvn.nutauVal;
1241  tcosmic.ncVal = cvn.ncVal;
1242  tcosmic.cosmicVal = cvn.cosmicVal;
1243  tcosmic.passSel = cvn.passSel;
1244  spill.cosmiccvn.push_back(tcosmic);
1245  }
1246  spill.ncosmiccvn = spill.cosmiccvn.size();
1247  }
1248  fSpillTree->Fill();
1249 
1250  // Pull the MC cycle number from metadata
1251  // before filling nuTree so cycle is filled
1252  // Assumes metadata module is being run...
1253  std::map<std::string, std::string> metadata =
1255  // Check it exists in the metadata, otherwise it defaults to -5
1256  if (metadata.count("simulated.cycle"))
1257  fCycle = std::stoi(metadata["simulated.cycle"]);
1258  if (metadata.count("simulated.batch"))
1259  fBatch = std::stoi(metadata["simulated.batch"]);
1260 
1261  // Fill the nuTree one per event, in sync with the spillTree
1262  if (fParams.FillNuTree()) {
1263  caf::SRSpillTruthBranch spillTruth;
1264  caf::SRSpillTruthBranch* pspillTruth = &spillTruth;
1265  fNuTree->SetBranchAddress("mc", &pspillTruth);
1266 
1267  for (const cheat::NeutrinoWithIndex& nuIdx :
1268  bt->allMCTruth()) { // NeutrinoInteractions()){
1269  // AddMCTruthToVec wants one of these, but most of this info is only
1270  // relevant if the truth is associated with a slice.
1271  const cheat::NeutrinoEffPur ep = {nuIdx.neutrinoInt, 0, 0, 0, 0, 0, 0, 0};
1272  std::vector<cheat::TrackIDE> dummy;
1273  std::vector<SRNeutrino> nu;
1274  AddMCTruthToVec(&ep, dummy, dummy, dummy, dummy, &nu, evt, spill);
1275  spillTruth.nu = nu[0];
1276 
1277  art::Ptr<simb::MCTruth> mctruth = nuIdx.neutrinoInt;
1278  const std::vector<art::Ptr<simb::MCTruth>> mctv(1, mctruth);
1279 
1281  FindManyPStrict<fxwgt::FluxWeights>(mctv, evt,
1283  fxwgt::FluxWeights flxrw;
1284  if (GetAssociatedProduct(fmFLXRW, 0, flxrw))
1285  spillTruth.nu.rwgt.ppfx = FluxReweights(flxrw);
1286  else
1287  spillTruth.nu.rwgt.ppfx.setDefault();
1288 
1290  FindManyPStrict<rwgt::GENIEReweightTable>(mctv, evt,
1293  if (GetAssociatedProduct(fmRW, 0, rw)){
1294  spillTruth.nu.rwgt.genie = GenieReweightTable(rw);
1295  }
1296 
1298  FindManyPStrict<g4rwgt::G4WeightTable>(mctv, evt,
1301  if (GetAssociatedProduct(fmG4RW, 0, g4w))
1302  spillTruth.nu.rwgt.geant4 = Geant4Reweights(g4w);
1303  else
1304  spillTruth.nu.rwgt.geant4.setDefault();
1305 
1306 
1307  fNuTree->Fill();
1308  spilltruthcol->push_back(spillTruth);
1309  } // end for nuIdx
1310  } // end if FillNuTree
1311 
1312  // Applying an event filter but we need the unfiltered events for POT/LT accounting.
1313  // Thus we need the same version of cafmaker to run over the filtered events and !filtered.
1314  // However, the point is to save time so we don't want slicer or calhit on the !filtered stream.
1315  if(fParams.ApplyingFilter() && !slices.isValid()) {
1316  *spillcol = spill;
1317 
1318  evt.put(std::move(srcol));
1319  evt.put(std::move(spillcol));
1320  if (fParams.FillNuTree()) evt.put(std::move(spilltruthcol));
1321  evt.put(std::move(srAssn));
1322 
1323  return;
1324  }
1325 
1327  FindManyPStrict<murem::MRCCParent>(slices, evt,
1329  if (fParams.MRMode()) {
1330  std::vector<std::vector<cheat::NeutrinoEffPur>> temptable // I know :)
1331  = bt->SlicesToMCTruthsTable(*parentSlices);
1332  // some mre slices may not have matched parents
1333  // so sort the table so that nth entry in table
1334  // is the truth for the parent of the nth mre
1335  // slice. If nth mre slice has no parent,
1336  // make an empty effpur entry.
1337  for (unsigned int isli = 0; isli < slices->size(); isli++) {
1339  if (GetAssociatedProduct(fmMrccPar, isli, parent)) {
1340  int idx = parent.ParentIndex();
1341  sliceTruthTable.push_back(temptable[idx]);
1342  } else {
1344  cheat::NeutrinoEffPur effPur =
1345  cheat::NeutrinoEffPur{null, -1, -1, 0, -1, -1, -1, -1};
1346  std::vector<cheat::NeutrinoEffPur> emptyEffPur(allNus.size(), effPur);
1347  sliceTruthTable.push_back(emptyEffPur);
1348  }
1349  } // end loop over mr slices
1350  } // end of mr mode
1351  else
1352  sliceTruthTable = bt->SlicesToMCTruthsTable(*slices);
1353 
1354  // Find favorite neutrino IDs using several metrics. Total energy
1355  // is the one we will use, for the favorite neutrino, but the others
1356  // will be stored.
1357  std::vector<int> faveIdsEff = bt->SliceToOrderedNuIds(
1358  allNus, sliceTruthTable, cheat::EffMetric, cheat::EffMetric);
1359  std::vector<int> faveIdsEnergy = bt->SliceToOrderedNuIds(
1360  allNus, sliceTruthTable, cheat::EnergyMetric, cheat::EnergyMetric);
1361  std::vector<int> faveIdsPur = bt->SliceToOrderedNuIds(
1362  allNus, sliceTruthTable, cheat::PurMetric, cheat::PurMetric);
1363  std::vector<int> faveIdsEffPur = bt->SliceToOrderedNuIds(
1364  allNus, sliceTruthTable, cheat::EffPurMetric, cheat::EffPurMetric);
1365  std::vector<int> faveIdsEffThenPur = bt->SliceToOrderedNuIds(
1366  allNus, sliceTruthTable, cheat::EffMetric, cheat::PurMetric);
1367 
1368  // Get all truth track hits in the event
1369  std::vector<cheat::TrackIDE> allTracks;
1370  std::vector<cheat::TrackIDE> sliceTracks;
1371  std::vector<cheat::TrackIDE> allTracksBirks;
1372  std::vector<cheat::TrackIDE> sliceTracksBirks;
1373  if (bt->HaveTruthInfo()) allTracks = bt->HitsToTrackIDE(*allHits);
1374  if (bt->HaveTruthInfo()) allTracksBirks = bt->HitsToTrackIDE(*allHits, true);
1375 
1376  // Get all association finder objects
1377  auto fmDiscrete2d =
1378  FindManyPStrict<rb::Track>(slices, evt, fParams.DiscreteTrack2dLabel());
1379  auto fmDiscrete =
1380  FindManyPStrict<rb::Track>(slices, evt, fParams.DiscreteTrackLabel());
1381  auto fmKalman =
1382  FindManyPStrict<rb::Track>(slices, evt, fParams.KalmanTrackLabel());
1383  auto fmCosmic =
1384  FindManyPStrict<rb::Track>(slices, evt, fParams.CosmicTrackLabel());
1385  auto fmWindow =
1386  FindManyPStrict<rb::Track>(slices, evt, fParams.WindowTrackLabel());
1387  auto fmElastic =
1388  FindManyPStrict<rb::Vertex>(slices, evt, fParams.ElasticArmsLabel());
1389  auto fmHoughVertex =
1390  FindManyPStrict<rb::Vertex>(slices, evt, fParams.HoughVertexLabel());
1391  auto fmVertexDT =
1392  FindManyPStrict<rb::Vertex>(slices, evt, fParams.VertexDTLabel());
1393  auto fmSlcME = FindManyPStrict<me::SlcME>(slices, evt, fParams.MELabel());
1394  auto fmJmShower =
1395  FindManyPStrict<jmshower::JMShower>(slices, evt, fParams.JMShowerLabel());
1396  auto fmNueSel =
1397  FindManyPStrict<jmshower::EID>(slices, evt, fParams.NueSelLabel());
1399  slices, evt, fParams.LIDLabel()); // one lid object per slice
1401  slices, evt, fParams.SlcLIDLabel()); // one lid object per slice
1402  art::FindMany<slid::ShowerLID> fmShwLID(slices, evt, fParams.LIDLabel());
1403  art::FindMany<slid::ShowerPID> fmShwPID(slices, evt, fParams.SPIDLabel());
1404  auto fmSliceLID = FindManyPStrict<SliceLID::Prediction>(
1405  slices, evt, fParams.SliceLIDLabel());
1406  auto fmSliceLID_opphorn = FindManyPStrict<SliceLID::Prediction>(
1407  slices, evt, fParams.SliceLIDOppHornLabel());
1408  auto fmNumuEnergy =
1409  FindManyPStrict<numue::NumuE>(slices, evt, fParams.NumuEnergyLabel());
1410  auto fmNumuLSTMEnergy = FindManyPStrict<LSTME::LSTMEnergy>(
1411  slices, evt, fParams.NumuLSTMEnergyLabel());
1412  auto fmNumuLSTMEnergy_opphorn = FindManyPStrict<LSTME::LSTMEnergy>(
1414  auto fmNuePresel = FindManyPStrict<presel::PreselObj>(
1415  slices, evt, fParams.NuePreselLabel());
1416  auto fmRockPresel = FindManyPStrict<presel::PreselObj>(
1417  slices, evt, fParams.RockPreselLabel());
1418  auto fmVeto = FindManyPStrict<presel::Veto>(slices, evt, fParams.VetoLabel());
1419  auto fmNueVeto =
1420  FindManyPStrict<presel::Veto>(slices, evt, fParams.NueVetoLabel());
1421  auto fmLem =
1422  FindManyPStrict<lem::PIDDetails>(slices, evt, fParams.LemLabel());
1423  auto fmRvp = FindManyPStrict<rvp::RVP>(slices, evt, fParams.RvpLabel());
1424  auto fmXnue = FindManyPStrict<xnue::Xnue>(slices, evt, fParams.XnueLabel());
1425  auto fmCosRej =
1426  FindManyPStrict<cosrej::CosRejObj>(slices, evt, fParams.CosRejLabel());
1427  auto fmNueCosRej =
1428  FindManyPStrict<cosrej::NueCosRej>(slices, evt, fParams.NueCosRejLabel());
1429  auto fmNCCosRej = FindManyPStrict<ncid::NCCosRej>(
1430  slices, evt, fParams.NCCosRejLabel()); /// NC Cosmic Rejection
1431  auto fmNCPi0BkgRej = FindManyPStrict<ncpi0::NCPi0BkgRej>(
1432  slices, evt, fParams.NCPi0BkgRejLabel()); /// NC Pi0 Bkg Rejection
1433 
1434  auto fmNuonE =
1435  FindManyPStrict<cvn::Result>(slices, evt, fParams.NuonELabel());
1436 
1437  auto fmNuonE_opphorn =
1438  FindManyPStrict<cvn::Result>(slices, evt, fParams.NuonEOppHornLabel());
1439 
1440  auto fmCVN_LoosePreselPtP =
1441  FindManyPStrict<cvn::Result>(slices, evt, fParams.CVNLabelLoosePreselPtP());
1442 
1443  auto fmCVN_LoosePreselPtP_opphorn =
1444  FindManyPStrict<cvn::Result>(slices, evt, fParams.CVNLabelLoosePreselPtPOppHorn());
1445 
1446 
1447 
1448  auto fmCVN_OldPresel =
1449  FindManyPStrict<cvn::Result>(slices, evt, fParams.CVNLabelOldPresel());
1450 
1451  auto fmCVN_OldPresel_opphorn =
1452  FindManyPStrict<cvn::Result>(slices, evt, fParams.CVNLabelOldPreselOppHorn());
1453 
1454  auto fmCVN_NoCosmics =
1455  FindManyPStrict<cvn::Result>(slices, evt, fParams.CVNLabelNoCosmics());
1456 
1457  auto fmCVN_NoCosmics_opphorn =
1458  FindManyPStrict<cvn::Result>(slices, evt, fParams.CVNLabelNoCosmicsOppHorn());
1459 
1460  auto fmWS =
1461  FindManyPStrict<cvn::Result>(slices, evt, fParams.WSCNNLabel());
1462 
1463  //auto fmCVNFeatures =
1464  // FindManyPStrict<cvn::Features>(slices, evt, fParams.CVNLabelLoosePreselPtP());
1465  auto fmCVNPixelMaps =
1466  FindManyPStrict<cvn::PixelMap>(slices, evt, fParams.CVNPixelMapLabel());
1467  auto fmCVNTrainingData =
1468  FindManyPStrict<cvn::TrainingData> (slices, evt, fParams.CVNTrainingDataLabel());
1469  auto fmEvtRegCVN =
1470  FindManyPStrict<cvn::RegNuResult>(slices, evt, fParams.RegCVNLabel());
1471  auto fmEvtRegCVN_opphorn =
1472  FindManyPStrict<cvn::RegNuResult>(slices, evt, fParams.RegCVNOppHornLabel());
1473  auto fmHadRegCVN =
1474  FindManyPStrict<cvn::RegHadronResult>(slices, evt, fParams.RegCVNLabel());
1475  auto fmHadRegCVN_opphorn =
1476  FindManyPStrict<cvn::RegHadronResult>(slices, evt, fParams.RegCVNOppHornLabel());
1477 
1478  //--------------------- Filling Vertex Energy variables------------------------------------
1479  auto fmNDVertexReco =
1480  FindManyPStrict<ndreco::NDRecoVertexObj>(slices, evt, fParams.NDRecoLabel());
1481 
1482 
1484  if (fParams.IsSinglesOverlay()) {
1485  GetByLabelStrict(evt, fParams.SingleMixerLabel(), singlepot);
1486  if (!singlepot.failedToGet()) {
1487  float totsinglepot = singlepot->spillpot;
1488  fTotalSinglePOT += totsinglepot;
1489  }
1490  }
1491 
1492  int nNonNoise = 0;
1493  for (const rb::Cluster& slice : *slices) {
1494  if (slice.IsNoise() || slice.NCell() == 0) continue;
1495  ++nNonNoise;
1496  }
1497  for (unsigned int sliceId = 0; sliceId < slices->size(); ++sliceId) {
1498  const rb::Cluster& slice = (*slices)[sliceId];
1499 
1500  if (slice.IsNoise() || slice.NCell() == 0) continue;
1501  // Because we don't care about the noise slice and slices with no hits.
1502  StandardRecord rec;
1503  StandardRecord* prec = &rec; // TTree wants a pointer-to-pointer
1504  fRecTree->SetBranchAddress("rec", &prec);
1505 
1506  if (bt->HaveTruthInfo()) sliceTracks = bt->HitsToTrackIDE(slice.AllCells());
1507  if (bt->HaveTruthInfo())
1508  sliceTracksBirks = bt->HitsToTrackIDE(slice.AllCells(), true);
1509  //#######################################################
1510  // Fill header.
1511  //#######################################################
1512  // Get metadata information for header
1513  unsigned int run = evt.run();
1514  unsigned int subrun = evt.subRun();
1515  unsigned int spillNum = evt.id().event();
1516 
1517  art::Timestamp ts = evt.time();
1518  const time_t timeSec = ts.timeHigh();
1519 
1520  struct tm* timeStruct = localtime(&timeSec);
1521 
1522  unsigned int maskstatus = 0;
1523  unsigned int dibmask = fMask;
1524 
1525  if (rh->GoodDiBlockMask() == ((1 << 16) - 1))
1526  maskstatus = 0;
1527  else if (rh->GoodDiBlockMask() == ((1 << 17) - 1))
1528  maskstatus = 2;
1529  else
1530  maskstatus = 1;
1531 
1532  unsigned int dibfirst = 0;
1533  unsigned int diblast = 0;
1534  if (dibmask) {
1535  int iD;
1536  iD = 0;
1537  while (!((dibmask >> iD) & 1)) iD++;
1538  dibfirst = iD + 1;
1539  iD = 0;
1540  while (dibmask >> iD) iD++;
1541  diblast = iD;
1542  }
1543 
1544  unsigned int nbadchan = bc->NBadInSubRun(subrun);
1545  unsigned int ntotchan = rh->NAnalysisChannels();
1546 
1547  rec.hdr = SRHeader();
1548 
1549  rec.hdr.run = run;
1550  rec.hdr.subrun = subrun;
1551  rec.hdr.cycle = fCycle;
1552  rec.hdr.batch = fBatch;
1553  rec.hdr.evt = spillNum;
1554  rec.hdr.subevt = sliceId;
1555  rec.hdr.ismc = !fIsRealData;
1556  rec.hdr.det = fDet;
1557  rec.hdr.blind = 0;
1558  rec.hdr.filt = rb::IsFiltered(evt, slices, sliceId);
1559 
1560  rec.hdr.subevtstarttime = slice.MinTNS();
1561  rec.hdr.subevtendtime = slice.MaxTNS();
1562  rec.hdr.subevtmeantime = slice.MeanTNS();
1563  rec.hdr.unixtime = timeSec;
1564 
1565  rec.hdr.dibfirst = dibfirst;
1566  rec.hdr.diblast = diblast;
1567  rec.hdr.dibmask = dibmask;
1568  rec.hdr.maskstatus = maskstatus;
1569 
1570  rec.hdr.year = timeStruct->tm_year + 1900;
1571  rec.hdr.month = timeStruct->tm_mon + 1;
1572  rec.hdr.day = timeStruct->tm_mday;
1573  rec.hdr.doy = timeStruct->tm_yday + 1;
1574  rec.hdr.hour = timeStruct->tm_hour + 1;
1575  rec.hdr.minute = timeStruct->tm_min;
1576  rec.hdr.second = timeStruct->tm_sec;
1577 
1578  rec.hdr.nbadchan = nbadchan;
1579  rec.hdr.ntotchan = ntotchan;
1580 
1581  rec.hdr.gain = rh->DetGainSetting();
1582  rec.hdr.finetiming = rh->DetFineTimingSetting();
1583 
1584  //#######################################################
1585  // Fill Spill info.
1586  //#######################################################
1587  // Store info on NuMI beam quality and if the pot as good or bad.
1588  // This is meant to be used to refine a good spill, not for
1589  // POT counting since events have multiple slices and empty
1590  // spill are missing from this tree. Use the friend tree for
1591  // POT counting.
1592  rec.spill = spill;
1593 
1594  //trim list of Cosmic CVN scores down to only the time window the
1595  //physics slice matches
1596  for (int iC = 0; iC < (int)rec.spill.cosmiccvn.size(); ++iC){
1597 
1598  if ((slice.MeanTNS() > rec.spill.cosmiccvn[iC].timeWinMin) &&
1599  (slice.MeanTNS() < rec.spill.cosmiccvn[iC].timeWinMax)) continue;
1600  rec.spill.cosmiccvn.erase(rec.spill.cosmiccvn.begin()+iC);
1601  --iC;
1602 
1603  }
1604  //if more then one time window matches the slice, keep the one with the
1605  //lowest cosmic CVN score
1606  if (rec.spill.cosmiccvn.size() > 1){
1607  SRCosmicCVN tcos = rec.spill.cosmiccvn[0];
1608  for (unsigned int iC = 1; iC < rec.spill.cosmiccvn.size(); ++iC){
1609  if (rec.spill.cosmiccvn[iC].cosmicVal < tcos.cosmicVal)
1610  tcos = rec.spill.cosmiccvn[iC];
1611  }
1612  rec.spill.cosmiccvn.clear();
1613  rec.spill.cosmiccvn.push_back(tcos);
1614  }
1615  rec.spill.ncosmiccvn = rec.spill.cosmiccvn.size();
1616 
1617 
1618  //#######################################################
1619  // Add slice info.
1620  //#######################################################
1621 
1622  FillSliceVars(slice, rec.slc);
1623  rec.slc.nnonnoise = nNonNoise;
1624  // find the slice with minimum mean time distance from this slice
1625  float mintimediff = 1e7; // initialize to a big number
1626  signed int mintimejsli = -1;
1627  for (unsigned int jsli = 0; jsli < slices->size(); jsli++) {
1628  const rb::Cluster& slicej = (*slices)[jsli];
1629  if (jsli == sliceId || slicej.IsNoise()) continue;
1630  float diff = fabs(slicej.MeanTNS() - slice.MeanTNS());
1631  if (diff < mintimediff) {
1632  mintimejsli = jsli;
1633  }
1634  mintimediff = std::min(mintimediff, diff);
1635  }
1636  if (mintimejsli != -1) {
1637  const rb::Cluster& slicejmin = (*slices)[mintimejsli];
1638  rec.slc.closestslicetime = slicejmin.MeanTNS() - slice.MeanTNS();
1639  rec.slc.closestslicenhit = slicejmin.NCell();
1640  rec.slc.closestslicecalE = slicejmin.CalorimetricEnergy();
1641  TVector3 minp, maxp;
1642  minp = slicejmin.MinXYZ();
1643  maxp = slicejmin.MaxXYZ();
1644  rec.slc.closestsliceminfromtop = livegeom->DistToTop(maxp);
1645  rec.slc.closestsliceminfrombottom = livegeom->DistToBottom(minp);
1646  rec.slc.closestsliceminfromfront = livegeom->DistToFront(minp);
1647  rec.slc.closestsliceminfromback = livegeom->DistToBack(maxp);
1648  rec.slc.closestsliceminfromeast = livegeom->DistToEdgeXMinus(minp);
1649  rec.slc.closestsliceminfromwest = livegeom->DistToEdgeXPlus(maxp);
1650  float mincelldistX = 55555;
1651  float mincelldistY = 55555;
1652  float mincelldist = 55555;
1653  float cellD = 2. * geom->Plane(0)->Cell(0)->HalfD();
1654  float cellW = 2. * geom->Plane(0)->Cell(0)->HalfW();
1655  for (unsigned int hitIdxX = 0; hitIdxX < slice.NXCell(); ++hitIdxX) {
1656  const art::Ptr<rb::CellHit>& chitX = slice.XCell(hitIdxX);
1657  for (unsigned int minhitIdxX = 0; minhitIdxX < slicejmin.NXCell();
1658  ++minhitIdxX) {
1659  const art::Ptr<rb::CellHit>& minchitX = slicejmin.XCell(minhitIdxX);
1660  float dPlaneX = (float)chitX->Plane() - (float)minchitX->Plane();
1661  float dCellX = (float)chitX->Cell() - (float)minchitX->Cell();
1662  float distX = sqrt((dPlaneX * cellD) * (dPlaneX * cellD) +
1663  (dCellX * cellW) * (dCellX * cellW));
1664  if (distX < mincelldistX) mincelldistX = distX;
1665  }
1666  }
1667  for (unsigned int hitIdxY = 0; hitIdxY < slice.NYCell(); ++hitIdxY) {
1668  const art::Ptr<rb::CellHit>& chitY = slice.YCell(hitIdxY);
1669  for (unsigned int minhitIdxY = 0; minhitIdxY < slicejmin.NYCell();
1670  ++minhitIdxY) {
1671  const art::Ptr<rb::CellHit>& minchitY = slicejmin.YCell(minhitIdxY);
1672  float dPlaneY = (float)chitY->Plane() - (float)minchitY->Plane();
1673  float dCellY = (float)chitY->Cell() - (float)minchitY->Cell();
1674  float distY = sqrt((dPlaneY * cellD) * (dPlaneY * cellD) +
1675  (dCellY * cellW) * (dCellY * cellW));
1676  if (distY < mincelldistY) mincelldistY = distY;
1677  }
1678  }
1679  if (mincelldistX < mincelldistY) {
1680  mincelldist = mincelldistX;
1681  } else {
1682  mincelldist = mincelldistY;
1683  }
1684  rec.slc.closestslicemindist = mincelldist;
1685  } else { // protection when mintimejsli == -1
1686  rec.slc.closestslicetime = -5;
1687  rec.slc.closestslicenhit = 0;
1688  rec.slc.closestslicecalE = -5;
1689  rec.slc.closestsliceminfromtop = -5;
1690  rec.slc.closestsliceminfrombottom = -5;
1691  rec.slc.closestsliceminfromfront = -5;
1692  rec.slc.closestsliceminfromback = -5;
1693  rec.slc.closestsliceminfromeast = -5;
1694  rec.slc.closestsliceminfromwest = -5;
1695  rec.slc.closestslicemindist = 55555;
1696  }
1697  if (fDet == kFARDET) {
1698  rec.sel.contain.nplanestofront = rec.slc.firstplane - (dibfirst - 1) * 64;
1699  rec.sel.contain.nplanestoback = ((diblast)*64) - 1 - rec.slc.lastplane;
1700  }
1701 
1702  if (fDet == kNEARDET) {
1703  rec.sel.contain.nplanestofront = rec.slc.firstplane;
1704  rec.sel.contain.nplanestoback = 214 - rec.slc.lastplane;
1705  }
1706 
1707  //#######################################################
1708  // Adding reconstructed objects.
1709  //#######################################################
1710 
1711  // Get all the tracks associated with this slice from discrete
1712  // tracker. This depends on the findMany object created above.
1713  if (fmDiscrete.isValid()) {
1714  std::vector<art::Ptr<rb::Track>> discTracks = fmDiscrete.at(sliceId);
1715  // Sort the tracks in descending length order
1716  std::sort(discTracks.begin(), discTracks.end(), sortRBTrackLength);
1717  rec.trk.discrete.tracks.resize(discTracks.size());
1718  for (size_t trackId = 0; trackId < discTracks.size(); ++trackId) {
1719  FillTrackVars((*discTracks[trackId]), rec.trk.discrete.tracks[trackId],
1720  *slices, allTracks, sliceTracks, allTracksBirks,
1721  sliceTracksBirks, sliceId);
1722  }
1723  }
1724 
1725  // if(fmDiscrete2d.isValid())
1726  // {
1727  // std::vector< art::Ptr<rb::Track> > discTracks =
1728  // fmDiscrete2d.at(sliceId);
1729  // Sort the tracks in descending length order
1730  // std::sort(discTracks.begin(), discTracks.end(),
1731  // sortRBTrackLength);
1732  // rec.trk.discrete.tracks2d.resize(discTracks.size());
1733  // for(size_t trackId = 0; trackId < discTracks.size(); ++trackId){
1734  // FillTrackVars(*discTracks[trackId],
1735  // rec.trk.discrete.tracks2d[trackId], *slices, sliceId);
1736  // }
1737  // }
1738 
1739  if (fmCosmic.isValid()) {
1740  std::vector<art::Ptr<rb::Track>> cosTracks = fmCosmic.at(sliceId);
1741  art::FindManyP<me::TrkME> fmCosmicTrkME(cosTracks, evt, fParams.MELabel());
1742 
1743  // Sort the tracks in descending length order
1744  std::sort(cosTracks.begin(), cosTracks.end(), sortRBTrackLength);
1745 
1746  rec.trk.cosmic.tracks.reserve(cosTracks.size());
1747  rec.trk.cosmic.tracks.resize(cosTracks.size());
1748  for (size_t trackId = 0; trackId < cosTracks.size(); ++trackId) {
1749  FillTrackVars((*cosTracks[trackId]), rec.trk.cosmic.tracks[trackId],
1750  *slices, allTracks, sliceTracks, allTracksBirks,
1751  sliceTracksBirks, sliceId);
1752 
1753  // Start of MR Shower
1754  // Fill the MR Shower informations
1755  art::FindManyP<rb::Cluster> fmDiF(cosTracks, evt,
1757  if (fmDiF.isValid()) {
1758  std::vector<art::Ptr<rb::Cluster>> difs = fmDiF.at(trackId);
1759  for (unsigned int difId = 0; difId < difs.size(); ++difId) {
1760  SRMRProperties srMR; // we're going to accumulate info into this
1761  art::Ptr<rb::Cluster> dif = difs[difId];
1762  // Filling DiF Cluster variables
1763  FillDiFVars(*dif, *slices, allTracks, sliceTracks, allTracksBirks,
1764  sliceTracksBirks, sliceId, srMR);
1765 
1766  // -- CVN - Loose Presel PtP Cut
1767  art::FindManyP<cvn::Result> fmDiFCVNResult_LoosePreselPtP(
1769  cvn::Result cvnDiFResult_LoosePreselPtP;
1770  if (GetAssociatedProduct(fmDiFCVNResult_LoosePreselPtP, difId,
1771  cvnDiFResult_LoosePreselPtP)) {
1772  FillCVNResultVars(cvnDiFResult_LoosePreselPtP, srMR.cvnloosepreselptp);
1773  } else {
1774  srMR.cvnloosepreselptp.setDefault();
1775  }
1776  // -- CVN - Old Presel
1777  art::FindManyP<cvn::Result> fmDiFCVNResult_OldPresel(
1779  cvn::Result cvnDiFResult_OldPresel;
1780  if (GetAssociatedProduct(fmDiFCVNResult_OldPresel, difId,
1781  cvnDiFResult_OldPresel)) {
1782  FillCVNResultVars(cvnDiFResult_OldPresel, srMR.cvnoldpresel);
1783  } else {
1784  srMR.cvnoldpresel.setDefault();
1785  }
1786  // -- CVN - No Cosmics used in training
1787  art::FindManyP<cvn::Result> fmDiFCVNResult_NoCosmics(
1789  cvn::Result cvnDiFResult_NoCosmics;
1790  if (GetAssociatedProduct(fmDiFCVNResult_NoCosmics, difId,
1791  cvnDiFResult_NoCosmics)) {
1792  FillCVNResultVars(cvnDiFResult_NoCosmics, srMR.cvnnocosmics, true);
1793  } else {
1794  srMR.cvnnocosmics.setDefault();
1795  }
1796 
1797  // Filling DiF LID
1798  art::FindManyP<rb::Shower> fmDiFShw({dif}, evt,
1800  if (fmDiFShw.isValid()) {
1801  std::vector<art::Ptr<rb::Shower>> difshws = fmDiFShw.at(difId);
1802  for (art::Ptr<rb::Shower> difshw : difshws) {
1803  if (difshw) {
1804  FillDiFShowerVars(*difshw, *slices, allTracks, sliceTracks,
1805  allTracksBirks, sliceTracksBirks, sliceId,
1806  srMR);
1807  } // end for if difshw
1808  } // end for loop difshw
1809  } // end if Shw assn DiF valid
1810 
1812  {dif}, evt, fParams.DiFShowerLIDLabel());
1813  if (fmDiFLID.isValid()) {
1814  std::vector<art::Ptr<slid::ShowerLID>> diflids =
1815  fmDiFLID.at(difId);
1816  for (art::Ptr<slid::ShowerLID> diflid : diflids) {
1817  if (diflid) {
1818  FillSlidVars(*diflid, srMR.lid);
1819  } // end for if diflid
1820  } // end for loop diflid
1821  } // end if LID assn DiF valid
1822 
1823  // Now the DiF cluster is ready
1824  rec.trk.cosmic.tracks[trackId].mrdif.push_back(srMR);
1825  } // end for dif cluster loop
1826  } // end if Track assn with DiF is valid
1827  // End of MR Shower
1828 
1829  if (fmCosmicTrkME.isValid()) {
1830  // Get MichelE from MEFinder
1831  std::vector<art::Ptr<me::TrkME>> michels = fmCosmicTrkME.at(trackId);
1832  for (art::Ptr<me::TrkME> michel : michels) {
1833  AddTrkMEToVec(*michel, &rec.me.trkcosmic,
1834  rec.trk.cosmic.tracks.back(), *slices, allTracks,
1835  sliceTracks, allTracksBirks, sliceTracksBirks,
1836  sliceId);
1837  }
1838  } // end
1839 
1840  } // end for trackId
1841  }
1842 
1843  // Do the same for windowtrack tracks
1844  if (fmWindow.isValid()) {
1845  std::vector<art::Ptr<rb::Track>> winTracks = fmWindow.at(sliceId);
1846  // Sort the tracks in descending length order
1847  std::sort(winTracks.begin(), winTracks.end(), sortRBTrackLength);
1848  rec.trk.window.tracks.reserve(winTracks.size());
1849  rec.trk.window.tracks.resize(winTracks.size());
1850  for (size_t trackId = 0; trackId < winTracks.size(); ++trackId) {
1851  FillTrackVars((*winTracks[trackId]), rec.trk.window.tracks[trackId],
1852  *slices, allTracks, sliceTracks, allTracksBirks,
1853  sliceTracksBirks, sliceId);
1854  } // end for trackId
1855  }
1856 
1857  // Now do the same for kalmanTrack
1858  if (fmKalman.isValid()) {
1859  std::vector<art::Ptr<rb::Track>> kalTracks = fmKalman.at(sliceId);
1860  // Sort the tracks in descending length order to start
1861  std::sort(kalTracks.begin(), kalTracks.end(), sortRBTrackLength);
1862  rec.trk.kalman.idxlongest = 0;
1863 
1864  art::FindManyP<remid::ReMId> fmRemid(kalTracks, evt,
1865  fParams.RemidLabel());
1866 
1867  art::FindManyP<muonid::MuonID> fmMuonid(kalTracks, evt,
1868  fParams.MuonIDLabel());
1869 
1871  FindManyPStrict<cosrej::TrkCntObj>(kalTracks, evt,
1872  fParams.CosRejLabel());
1874  FindManyPStrict<ndreco::NDRecoTrkObj>(kalTracks, evt,
1875  fParams.NDRecoLabel());
1876  art::FindManyP<me::TrkME> fmKalmanME(kalTracks, evt, fParams.MELabel());
1877  art::FindManyP<rb::Energy> fmKalOverlapE(
1878  kalTracks, evt,fParams.OverlapEKalLabel());
1879 
1881  FindManyPStrict<trackinfo::TrackInfoObj>(kalTracks, evt,
1883 
1884  // Give SRRemid the pid values of track with largest PID
1885  unsigned int bestPidIdx =
1887 
1888  unsigned int bestPidIdx3d = 999;
1889 
1890  // Give SRMuonID the pid values of track with largest PID
1891  unsigned int bestMuonidIdx =
1893  unsigned int bestMuonidIdx3d = 999;
1894 
1895 
1896 
1897  for (size_t trackId = 0; trackId < kalTracks.size(); ++trackId) {
1898  if (kalTracks[trackId]->View() != geo::kXorY) {
1899  rec.trk.kalman.tracks2d.push_back(SRTrack());
1900 
1901  // Fill the Track Vars with overlapping E if that info is present,
1902  // otherwise use the default FillTrackVars function.
1903  if (fmKalOverlapE.isValid()) {
1904  std::vector<art::Ptr<rb::Energy>> overlapEs =
1905  fmKalOverlapE.at(trackId);
1906  if (!overlapEs.empty()) {
1907  FillTrackVarsWithOverlapE((*kalTracks[trackId]), (*overlapEs[0]),
1908  rec.trk.kalman.tracks2d.back(), *slices,
1909  allTracks, sliceTracks, allTracksBirks,
1910  sliceTracksBirks, sliceId);
1911  } else {
1912  FillTrackVars((*kalTracks[trackId]),
1913  rec.trk.kalman.tracks2d.back(), *slices, allTracks,
1914  sliceTracks, allTracksBirks, sliceTracksBirks,
1915  sliceId);
1916  }
1917  } else {
1918  FillTrackVars((*kalTracks[trackId]), rec.trk.kalman.tracks2d.back(),
1919  *slices, allTracks, sliceTracks, allTracksBirks,
1920  sliceTracksBirks, sliceId);
1921  }
1922  if(fmTrkInfo.isValid()){
1923  std::vector<art::Ptr<trackinfo::TrackInfoObj> > trkinfos = fmTrkInfo.at(trackId);
1924  //if(!trkinfos.empty())
1925  // FillTrackInfoVars(*trkinfos[0], rec.trk.kalman.tracks.back());
1926  }
1927  } else {
1928  rec.trk.kalman.tracks.push_back(SRKalmanTrack());
1929 
1930  // Fill the Track Vars with overlapping E if that info is present,
1931  // otherwise use the default FillTrackVars function.
1932  if (fmKalOverlapE.isValid()) {
1933  std::vector<art::Ptr<rb::Energy>> overlapEs =
1934  fmKalOverlapE.at(trackId);
1935  if (!overlapEs.empty()) {
1936  FillTrackVarsWithOverlapE((*kalTracks[trackId]), (*overlapEs[0]),
1937  rec.trk.kalman.tracks.back(), *slices,
1938  allTracks, sliceTracks, allTracksBirks,
1939  sliceTracksBirks, sliceId);
1940  } else {
1941  FillTrackVars((*kalTracks[trackId]), rec.trk.kalman.tracks.back(),
1942  *slices, allTracks, sliceTracks, allTracksBirks,
1943  sliceTracksBirks, sliceId);
1944  }
1945  } else {
1946  FillTrackVars((*kalTracks[trackId]), rec.trk.kalman.tracks.back(),
1947  *slices, allTracks, sliceTracks, allTracksBirks,
1948  sliceTracksBirks, sliceId);
1949  }
1950 
1951  if(fmTrkInfo.isValid()){
1952  std::vector<art::Ptr<trackinfo::TrackInfoObj> > trkinfos = fmTrkInfo.at(trackId);
1953  if(!trkinfos.empty())
1954  FillTrackInfoVars(*trkinfos[0], rec.trk.kalman.tracks.back());
1955  }
1956 
1957  if (trackId == bestPidIdx)
1958  bestPidIdx3d = rec.trk.kalman.tracks.size() - 1;
1959 
1960  if (trackId == bestMuonidIdx)
1961  bestMuonidIdx3d = rec.trk.kalman.tracks.size() - 1;
1962 
1963  if (fmRemid.isValid()) {
1964  // Get MuID from ReMId
1965  std::vector<art::Ptr<remid::ReMId>> remids = fmRemid.at(trackId);
1966 
1967  if (!remids.empty()) {
1968  FillMuIdVars(*remids[0], rec.trk.kalman.tracks.back());
1969  } else {
1970  std::cout << "Missing Association Error: Failed to find a ReMId "
1971  "object with this track."
1972  << std::endl;
1973  abort();
1974  }
1975  }
1976 
1977  if (fmMuonid.isValid()) {
1978  std::vector<art::Ptr<muonid::MuonID>> muonids = fmMuonid.at(trackId);
1979 
1980  if (!muonids.empty())
1981  FillMuonIDVars(*muonids[0], rec.trk.kalman.tracks.back());
1982  else {
1983  std::cout << "Missing Association Error: Failed to find a MuonID "
1984  "object with this track."
1985  << std::endl;
1986  abort();
1987  }
1988  }
1989 
1990  if (fmTrkCnt.isValid()) {
1991  // Get track containment info from TrkCntObj
1992  std::vector<art::Ptr<cosrej::TrkCntObj>> trkcnts =
1993  fmTrkCnt.at(trackId);
1994 
1995  if (!trkcnts.empty())
1996  FillTrackContainmentVars(*trkcnts[0],
1997  rec.trk.kalman.tracks.back());
1998  }
1999 
2000  if (fmNDTrkReco.isValid()) {
2001  // Get pion energy reconstruction variables info from NDRecoObj
2002  std::vector<art::Ptr<ndreco::NDRecoTrkObj>> ndtrkrecos =
2003  fmNDTrkReco.at(trackId);
2004 
2005  if (!ndtrkrecos.empty())
2006  FillNDRecoTrkVars(*ndtrkrecos[0], rec.trk.kalman.tracks.back());
2007  }
2008 
2009  if (fmKalmanME.isValid()) {
2010  // Get MichelE from MEFinder
2011  std::vector<art::Ptr<me::TrkME>> michels = fmKalmanME.at(trackId);
2012 
2013  for (art::Ptr<me::TrkME> michel : michels) {
2014  AddTrkMEToVec(*michel, &rec.me.trkkalman,
2015  rec.trk.kalman.tracks.back(), *slices, allTracks,
2016  sliceTracks, allTracksBirks, sliceTracksBirks,
2017  sliceId);
2018  }
2019  } // end
2020  }
2021  } // end for trackId
2022 
2023  if (!rec.trk.kalman.tracks.empty() && bestPidIdx3d != 999) {
2024  CopyRemidVars(rec.trk.kalman.tracks[bestPidIdx3d], rec.sel.remid);
2025  std::vector<art::Ptr<remid::ReMId>> remids = fmRemid.at(bestPidIdx3d);
2026  rec.trk.kalman.idxremid = bestPidIdx3d;
2027  } else {
2028  rec.sel.remid.setDefault();
2029  }
2030  if (!rec.trk.kalman.tracks.empty() && bestMuonidIdx3d != 999) {
2031  CopyMuonIDVars(rec.trk.kalman.tracks[bestMuonidIdx3d], rec.sel.muonid);
2032  std::vector<art::Ptr<muonid::MuonID>> muonids = fmMuonid.at(bestMuonidIdx3d);
2033  rec.trk.kalman.idxmuonid = bestMuonidIdx3d;
2034  }
2035  else
2036  rec.sel.muonid.setDefault();
2037 
2038  } else {
2039  rec.sel.remid.setDefault();
2040  rec.sel.muonid.setDefault();
2041  }
2042 
2043  std::sort(rec.trk.kalman.tracks.begin(), rec.trk.kalman.tracks.end(),
2044  sortRemid);
2045  rec.trk.kalman.idxremid = 0;
2046 
2047  // Now that everything is unsorted, need to re-find best muonid track
2048  rec.trk.kalman.idxmuonid = std::max_element(rec.trk.kalman.tracks.begin(),
2049  rec.trk.kalman.tracks.end(),
2050  compMuonID)
2051  - rec.trk.kalman.tracks.begin();
2052 
2053  unsigned int tempidx = 0, templen = 0;
2054  for (size_t trackId = 0; trackId < rec.trk.kalman.tracks.size();
2055  ++trackId) {
2056  if (rec.trk.kalman.tracks[trackId].len > templen) {
2057  tempidx = trackId;
2058  templen = rec.trk.kalman.tracks[trackId].len;
2059  }
2060  }
2061 
2062  rec.trk.kalman.idxlongest = tempidx;
2063  rec.trk.kalman.fillSizes();
2064  rec.trk.discrete.fillSizes();
2065  rec.trk.window.fillSizes();
2066  rec.trk.cosmic.fillSizes();
2067 
2068  // pre-fill nue energy variables with default
2069  rec.energy.nue.setDefault();
2070  rec.energy.nue.lid.setDefault();
2071 
2072  // get vertices from ElasticArms
2073  if (fmElastic.isValid()) {
2074 
2075 
2076  std::vector<art::Ptr<rb::Vertex>> elastics = fmElastic.at(sliceId);
2077  if(elastics.size() > 0){ //this should always be size 1 or 0
2078 
2079  rec.vtx.elastic = SRElastic();
2080  int vtxId = 0; //always check the 0th entry since there should only be 1.
2081  //will remove this and replace all vtxId with 0 once we are certain the
2082  //code all works.
2083  SRElastic& srVtx = rec.vtx.elastic;
2084  srVtx.time = elastics[vtxId]->GetT();
2085  srVtx.vtx = elastics[vtxId]->GetXYZ();
2086 
2087  if (fmNDVertexReco.isValid())
2088  {
2089  std::vector<art::Ptr<ndreco::NDRecoVertexObj>> ndvertexrecos =
2090  fmNDVertexReco.at(sliceId);
2091  if (ndvertexrecos.size() > 0)
2092  {
2093  srVtx.slicevertexenergyvolume10 = ndvertexrecos[0]->SliceVertexE10();
2094  srVtx.slicevertexenergyvolume20 = ndvertexrecos[0]->SliceVertexE20();
2095  srVtx.slicevertexenergyvolume30 = ndvertexrecos[0]->SliceVertexE30();
2096  srVtx.slicevertexenergyvolume40 = ndvertexrecos[0]->SliceVertexE40();
2097  srVtx.prong3dvertexenergyvolume10 = ndvertexrecos[0]->Prong3DVertexE10();
2098  srVtx.prong3dvertexenergyvolume20 = ndvertexrecos[0]->Prong3DVertexE20();
2099  srVtx.prong3dvertexenergyvolume30 = ndvertexrecos[0]->Prong3DVertexE30();
2100  srVtx.prong3dvertexenergyvolume40 = ndvertexrecos[0]->Prong3DVertexE40();
2101  srVtx.prong2dvertexenergyvolume10 = ndvertexrecos[0]->Prong2DVertexE10();
2102  srVtx.prong2dvertexenergyvolume20 = ndvertexrecos[0]->Prong2DVertexE20();
2103  srVtx.prong2dvertexenergyvolume30 = ndvertexrecos[0]->Prong2DVertexE30();
2104  srVtx.prong2dvertexenergyvolume40 = ndvertexrecos[0]->Prong2DVertexE40();
2105  }
2106  }
2107 
2108  srVtx.IsValid = true;
2109  // for each vertex (currently only one elastic arms ever), look for
2110  // prongs (currently only fuzzyk)
2111  art::FindManyP<rb::Prong> fmFuzzyProng2D = FindManyPStrict<rb::Prong>(
2112  elastics, evt, fParams.FuzzyK2DLabel()); // 2D prongs first
2113  if (fmFuzzyProng2D.isValid()) {
2114  std::vector<art::Ptr<rb::Prong>> prongs = fmFuzzyProng2D.at(vtxId);
2115  std::sort(prongs.begin(), prongs.end(),
2116  sortRBProngLength); // 2D prongs are sorted by length
2118  prongs, evt,fParams.FuzzyKWeight2DLabel());
2119 
2120  for (unsigned int iPng = 0; iPng < prongs.size(); ++iPng) {
2121  art::Ptr<rb::Prong> prong = prongs[iPng];
2122  srVtx.fuzzyk.png2d.push_back(SRProng());
2123  FillProngVars(*prong, srVtx.fuzzyk.png2d.back(), *slices, allTracks,
2124  sliceTracks, allTracksBirks, sliceTracksBirks,
2125  sliceId);
2126 
2127  art::FindManyP<rb::PID> fmCVNNeutronDaughterResult(
2128  prongs, evt, fParams.CVNLabelNeutronProng());
2129  if (fmCVNNeutronDaughterResult.isValid()){
2130  std::vector<art::Ptr<rb::PID>> cvnneutrons =
2131  fmCVNNeutronDaughterResult.at(iPng);
2132  if(!cvnneutrons.empty()){
2133  FillCVNNeutronDaughterResultVars(cvnneutrons, srVtx.fuzzyk.png2d.back().cvnneutron);
2134  }
2135  }
2136 
2137  art::FindManyP<rb::PID> fmCVNParticleResult(
2138  prongs, evt, fParams.CVNParticle2DLabel());
2139  if (fmCVNParticleResult.isValid()) {
2140  std::vector<art::Ptr<rb::PID>> cvnparts =
2141  fmCVNParticleResult.at(iPng);
2142  if(!cvnparts.empty()){
2143  FillCVNParticleResultVars(cvnparts, srVtx.fuzzyk.png2d.back().cvnpart);
2144  }
2145  }
2146  art::FindManyP<cvn::PixelMap> fmCVNProngPixelMaps(
2147  prongs, evt, fParams.CVNPixelMapLabel());
2148  art::FindManyP<cvn::ProngTrainingData> fmCVNProngTrainingData(
2149  prongs, evt, fParams.CVNTrainingDataLabel());
2150 
2151  if (fmCVNProngPixelMaps.isValid() && fParams.FillPixelMaps()) {
2152  std::vector<art::Ptr<cvn::PixelMap>> cvnProngPixelMaps =
2153  fmCVNProngPixelMaps.at(iPng);
2154  for (unsigned int ip = 0; ip < cvnProngPixelMaps.size(); ++ip) {
2155  srVtx.fuzzyk.png2d.back().cvnmaps.push_back(SRPixelMap());
2156  SRPixelMap& pmap = srVtx.fuzzyk.png2d.back().cvnmaps.back();
2157  FillCVNPixelMaps(*cvnProngPixelMaps[ip], pmap,
2159  }
2160  }
2161 
2162  cvn::ProngTrainingData cvnProngTrainingData;
2163  if (GetAssociatedProduct(fmCVNProngTrainingData, iPng,
2164  cvnProngTrainingData) &&
2166  srVtx.fuzzyk.png2d.back().prongtrainingdata.push_back(
2167  SRProngTrainingData());
2168  SRProngTrainingData& srProngTrainingData =
2169  srVtx.fuzzyk.png2d.back().prongtrainingdata.back();
2170  FillCVNProngTrainingData(cvnProngTrainingData,
2171  srProngTrainingData);
2172  }
2173 
2174  if (fmPngPngWt.isValid()) {
2175  std::vector<art::Ptr<rb::WeightedProng>> prongWt =
2176  fmPngPngWt.at(iPng);
2177  srVtx.fuzzyk.png2d.back().weightedCalE =
2178  prongWt[0]->CalorimetricEnergy();
2179  }
2180 
2181  } // end loop over 2D prongs
2182  } // end if 2D prongs are valid for this vertex
2183 
2184  // 3D prongs next
2185  art::FindManyP<rb::Prong> fmFuzzyProng3D = FindManyPStrict<rb::Prong>(
2186  elastics, evt,fParams.FuzzyK3DLabel());
2187  unsigned int nshwlids = 0;
2188  if (fmFuzzyProng3D.isValid()) {
2189  std::vector<art::Ptr<rb::Prong>> prongs = fmFuzzyProng3D.at(vtxId);
2190  if (!prongs.empty()) {
2191  std::sort(prongs.begin(), prongs.end(),
2192  sortRBProngLength); // start sorting by length
2193  unsigned int pngId = 0;
2194  srVtx.fuzzyk.longestidx = 0;
2195  art::FindOneP<rb::Shower> fmShowerLID(prongs, evt,
2197  art::FindManyP<rb::Track> fmBPF(prongs, evt,
2199  art::FindManyP<rb::PID> fmCVNNeutronDaughterResult(
2200  prongs, evt, fParams.CVNLabelNeutronProng() );
2201  art::FindManyP<rb::PID> fmCVNParticleResult(
2202  prongs, evt, fParams.CVNParticleLabel());
2203  art::FindManyP<rb::PID> fmCVNParticleResult_opphorn(
2204  prongs, evt, fParams.CVNParticleOppHornLabel());
2205  art::FindManyP<rb::PID> fmSPProngCVN5labelParticleResult(
2206  prongs, evt, fParams.SPProngCVN5labelLabel());
2207  art::FindManyP<rb::PID> fmSPProngCVNNumuCCEMIDParticleResult(prongs,
2208  evt,
2210  art::FindManyP<cvn::PixelMap> fmCVNProngPixelMaps(
2211  prongs, evt, fParams.CVNPixelMapLabel());
2212  art::FindManyP<cvn::ProngTrainingData> fmCVNProngTrainingData(
2213  prongs, evt, fParams.CVNTrainingDataLabel());
2214  art::FindManyP<cvn::RegProngResult> fmRegCVNResult(prongs, evt,
2215  fParams.RegCVNLabel());
2216 
2218  prongs, evt,fParams.FuzzyKWeight3DLabel());
2220  FindManyPStrict<ndreco::NDRecoPngObj>(prongs, evt,
2221  fParams.NDRecoLabel());
2222  float highE = 0;
2223  for (unsigned int iPng = 0; iPng < prongs.size(); ++iPng) {
2224  art::Ptr<rb::Prong> prong = prongs[iPng];
2225  srVtx.fuzzyk.png.push_back(SRFuzzyKProng());
2226  srVtx.fuzzyk.png.back().setDefault();
2227  FillProngVars(*prong, srVtx.fuzzyk.png.back(), *slices, allTracks,
2228  sliceTracks, allTracksBirks, sliceTracksBirks,
2229  sliceId);
2230  if (fmPngPngWt.isValid()) {
2231  std::vector<art::Ptr<rb::WeightedProng>> prongWt =
2232  fmPngPngWt.at(iPng);
2233  srVtx.fuzzyk.png.back().weightedCalE =
2234  prongWt[0]->CalorimetricEnergy();
2235  }
2236 
2237  if (fmCVNProngPixelMaps.isValid() && fParams.FillPixelMaps()) {
2238  std::vector<art::Ptr<cvn::PixelMap>> cvnProngPixelMaps =
2239  fmCVNProngPixelMaps.at(iPng);
2240  for (unsigned int ip = 0; ip < cvnProngPixelMaps.size(); ++ip) {
2241  srVtx.fuzzyk.png.back().cvnmaps.push_back(SRPixelMap());
2242  SRPixelMap& pmap = srVtx.fuzzyk.png.back().cvnmaps.back();
2243  FillCVNPixelMaps(*cvnProngPixelMaps[ip], pmap,
2245  }
2246  }
2247 
2248  cvn::ProngTrainingData cvnProngTrainingData;
2249  if (GetAssociatedProduct(fmCVNProngTrainingData, iPng,
2250  cvnProngTrainingData) &&
2252  srVtx.fuzzyk.png.back().prongtrainingdata.push_back(
2253  SRProngTrainingData());
2254  SRProngTrainingData& srProngTrainingData =
2255  srVtx.fuzzyk.png.back().prongtrainingdata.back();
2256 
2257  FillCVNProngTrainingData(cvnProngTrainingData,
2258  srProngTrainingData);
2259  srVtx.fuzzyk.png.back().cvnneutron.setDefault();
2260 
2261  }
2262 
2263  if (fmCVNNeutronDaughterResult.isValid()){
2264  std::vector<art::Ptr<rb::PID>> cvnneutrons =
2265  fmCVNNeutronDaughterResult.at(iPng);
2266  FillCVNNeutronDaughterResultVars(cvnneutrons, srVtx.fuzzyk.png.back().cvnneutron);
2267  }
2268 
2269  if (fmCVNParticleResult.isValid()) {
2270  std::vector<art::Ptr<rb::PID>> cvnparts =
2271  fmCVNParticleResult.at(iPng);
2272  FillCVNParticleResultVars(cvnparts, srVtx.fuzzyk.png.back().cvnpart);
2273  }
2274  if (fmCVNParticleResult_opphorn.isValid()) {
2275  std::vector<art::Ptr<rb::PID>> cvnparts =
2276  fmCVNParticleResult_opphorn.at(iPng);
2277  FillCVNParticleResultVars(cvnparts, srVtx.fuzzyk.png.back().cvnpart_opphorn);
2278  }
2279  if (fmSPProngCVN5labelParticleResult.isValid()) {
2280  std::vector<art::Ptr<rb::PID>> cvnparts = fmSPProngCVN5labelParticleResult.at(iPng);
2281  FillCVNParticleResultVars(cvnparts, srVtx.fuzzyk.png.back().spprongcvnpart5label);
2282  }
2283  if (fmSPProngCVNNumuCCEMIDParticleResult.isValid()) {
2284  std::vector<art::Ptr<rb::PID>> cvnparts = fmSPProngCVNNumuCCEMIDParticleResult.at(iPng);
2285  FillCVNParticleResultVars(cvnparts, srVtx.fuzzyk.png.back().spprongcvnpartnumuccemid);
2286  }
2287  if (fmRegCVNResult.isValid()) {
2288  std::vector<art::Ptr<cvn::RegProngResult>> regcvn =
2289  fmRegCVNResult.at(iPng);
2290  if(regcvn.size()>0)
2291  srVtx.fuzzyk.png.back().regcvn.prongE = regcvn[0]->Output();
2292  else
2293  srVtx.fuzzyk.png.back().regcvn.prongE = -5.;
2294  }
2295  if (fmNDPngReco.isValid()) {
2296  std::vector<art::Ptr<ndreco::NDRecoPngObj>> ndpngrecos =
2297  fmNDPngReco.at(iPng);
2298 
2299  if (!ndpngrecos.empty())
2300  FillNDRecoPngVars(*ndpngrecos[0], srVtx.fuzzyk.png.back());
2301  }
2302 
2303  if (fmShowerLID.isValid()) {
2304  art::Ptr<rb::Shower> slid = fmShowerLID.at(pngId);
2305  // look for associated LID reconstructions
2306  art::FindOneP<slid::ShowerLID> fmLID({slid}, evt,
2307  fParams.LIDLabel());
2308  if (slid) {
2309  // fill fuzzyk.png.shwlid; assume 1 shower per prong
2310  FillShowerVars(*slid, srVtx.fuzzyk.png.back(), *slices,
2311  allTracks, sliceTracks, allTracksBirks,
2312  sliceTracksBirks, sliceId);
2313  nshwlids++;
2314  if (fmLID.isValid()) {
2315  art::Ptr<slid::ShowerLID> lid = fmLID.at(0);
2316  // fill fuzzyk.png.shwlid.(lid/lidE); assume 1 LID info
2317  // object per shower
2318  if (lid) {
2319  FillSlidVars(*lid, srVtx.fuzzyk.png.back().shwlid);
2320  FillLIDEnergyVars(*lid, srVtx.fuzzyk.png.back());
2321  srVtx.fuzzyk.png.back().shwlid.shwE =
2322  srVtx.fuzzyk.png.back().shwlid.lidE.shwE;
2323  if (srVtx.fuzzyk.png.back().shwlid.shwE > highE) {
2324  highE = srVtx.fuzzyk.png.back().shwlid.shwE;
2325  // FillNueEnergyVars(slice, *slid,
2326  //rec.energy.nue);
2327  }
2328  }
2329  } // end if LID assn valid
2330  }
2331  } // end if showerlid is valid
2332 
2333  if (fmBPF.isValid()) {
2334 
2335  // BPF attempts (and rarely fails) to fit each prong under 3 different
2336  // particle assumptions: muon, pion, proton. So for each prong, there
2337  // will be at most 3 tracks. Each track _should_ have one and only one
2338  // FitSum, OverlapE, and TrkCntObj associated with it.
2339 
2340  // So... ...get all of these objects etc.
2341  std::vector<art::Ptr<rb::Track>> bpfTracks = fmBPF.at(pngId);
2342  art::FindManyP<rb::FitSum> fmBPFFitSums(
2343  bpfTracks, evt, fParams.BreakPointTrackLabel());
2344  art::FindManyP<rb::Energy> fmBPFOverlapE(
2345  bpfTracks, evt,fParams.OverlapEBPFLabel());
2346  art::FindManyP<cosrej::TrkCntObj> fmTrkCntBPF =
2347  FindManyPStrict<cosrej::TrkCntObj>(
2348  bpfTracks, evt, fParams.CosRejBPFLabel());
2349  art::FindManyP<ndreco::NDRecoBpfObj> fmTrkNDRecoBPF =
2350  FindManyPStrict<ndreco::NDRecoBpfObj>(
2351  bpfTracks, evt, fParams.NDRecoLabel());
2353  FindManyPStrict<trackinfo::TrackInfoObj>(
2354  bpfTracks, evt, fParams.TrackInfoBPF());
2355 
2356  for (size_t trackId = 0; trackId < bpfTracks.size();
2357  ++trackId) {
2358  if (fmBPFFitSums.isValid()) {
2359  std::vector<art::Ptr<rb::FitSum>> bpfFitSums =
2360  fmBPFFitSums.at(trackId);
2361  if (bpfFitSums.size() < 1) {
2362  std::cerr << "bpfFitSums.size() < 1, this should be 1, "
2363  "since each track has a fitsum!"
2364  << std::endl;
2365  continue;
2366  }
2367 
2368  //
2369  // Here we should fill the track and fitsum info
2370  //
2371 
2372  // For easy reference, pull out this specific prong.
2373  SRFuzzyKProng& srPng = srVtx.fuzzyk.png.back();
2374 
2375 
2376 
2377  // Sorry for the cut and paste, but it is only ~20 lines of actual code.
2378  // Fill the track if it is for the muon, then repeat for pion and proton.
2379 
2380  // If this is the muon, fill the muon track info...
2381  if(abs(bpfFitSums[0]->PDG()) == 13) {
2382 
2383  srPng.bpf.muon.IsValid = true;
2384 
2385  // Fill the Track Vars with overlapping E if that info is
2386  // present, otherwise use the default FillTrackVars
2387  // function.
2388  if (fmBPFOverlapE.isValid()) {
2389  std::vector<art::Ptr<rb::Energy>> overlapEs =
2390  fmBPFOverlapE.at(trackId);
2391  if (!overlapEs.empty()) {
2393  (*bpfTracks[trackId]), (*overlapEs[0]), srPng.bpf.muon,
2394  *slices, allTracks, sliceTracks, allTracksBirks,
2395  sliceTracksBirks, sliceId);
2396  } else {
2397  FillTrackVars((*bpfTracks[trackId]), srPng.bpf.muon,
2398  *slices, allTracks, sliceTracks,
2399  allTracksBirks, sliceTracksBirks,
2400  sliceId);
2401  }
2402  } else {
2403  FillTrackVars((*bpfTracks[trackId]), srPng.bpf.muon, *slices,
2404  allTracks, sliceTracks, allTracksBirks,
2405  sliceTracksBirks, sliceId);
2406  }
2407 
2408  // Fill the FitSum info
2409  FillTrackVarsBpfFitSum((*bpfFitSums[0]), srPng.bpf.muon);
2410  if (fmTrkCntBPF.isValid()) {
2411  // Get track containment info from TrkCntObj
2412  std::vector<art::Ptr<cosrej::TrkCntObj>> trkcnts =
2413  fmTrkCntBPF.at(trackId);
2414 
2415  if (!trkcnts.empty())
2416  FillTrackContainmentVars(*trkcnts[0], srPng.bpf.muon);
2417  }
2418 
2419  if (fmTrkNDRecoBPF.isValid()) {
2420  std::vector<art::Ptr<ndreco::NDRecoBpfObj>> ndrecos =
2421  fmTrkNDRecoBPF.at(trackId);
2422 
2423  if (!ndrecos.empty())
2424  FillNDRecoBPFVars(*ndrecos[0], srPng.bpf.muon);
2425  }
2426 
2427  // Fill the TrackInfo info (yes I said info twice...)
2428  if(fmTrkInfoBPF.isValid()){
2429  std::vector<art::Ptr<trackinfo::TrackInfoObj> > trkinfos = fmTrkInfoBPF.at(trackId);
2430  if(!trkinfos.empty())
2431  FillTrackInfoVars(*trkinfos[0], srPng.bpf.muon);
2432  }
2433 
2434  } // end if(PDG == 13)
2435 
2436  // If this is the pion, fill the pion track info...
2437  else if(abs(bpfFitSums[0]->PDG()) == 211) {
2438 
2439  srPng.bpf.pion.IsValid = true;
2440 
2441  // Fill the Track Vars with overlapping E if that info is
2442  // present, otherwise use the default FillTrackVars
2443  // function.
2444  if (fmBPFOverlapE.isValid()) {
2445  std::vector<art::Ptr<rb::Energy>> overlapEs =
2446  fmBPFOverlapE.at(trackId);
2447  if (!overlapEs.empty()) {
2449  (*bpfTracks[trackId]), (*overlapEs[0]), srPng.bpf.pion,
2450  *slices, allTracks, sliceTracks, allTracksBirks,
2451  sliceTracksBirks, sliceId);
2452  } else {
2453  FillTrackVars((*bpfTracks[trackId]), srPng.bpf.pion,
2454  *slices, allTracks, sliceTracks,
2455  allTracksBirks, sliceTracksBirks,
2456  sliceId);
2457  }
2458  } else {
2459  FillTrackVars((*bpfTracks[trackId]), srPng.bpf.pion, *slices,
2460  allTracks, sliceTracks, allTracksBirks,
2461  sliceTracksBirks, sliceId);
2462  }
2463 
2464  // Fill the FitSum info
2465  FillTrackVarsBpfFitSum((*bpfFitSums[0]), srPng.bpf.pion);
2466  if (fmTrkCntBPF.isValid()) {
2467  // Get track containment info from TrkCntObj
2468  std::vector<art::Ptr<cosrej::TrkCntObj>> trkcnts =
2469  fmTrkCntBPF.at(trackId);
2470 
2471  if (!trkcnts.empty())
2472  FillTrackContainmentVars(*trkcnts[0], srPng.bpf.pion);
2473  }
2474 
2475  if (fmTrkNDRecoBPF.isValid()) {
2476  std::vector<art::Ptr<ndreco::NDRecoBpfObj>> ndrecos =
2477  fmTrkNDRecoBPF.at(trackId);
2478 
2479  if (!ndrecos.empty())
2480  FillNDRecoBPFVars(*ndrecos[0], srPng.bpf.pion);
2481  }
2482 
2483  // Fill the TrackInfo info (yes I said info twice...)
2484  if(fmTrkInfoBPF.isValid()){
2485  std::vector<art::Ptr<trackinfo::TrackInfoObj> > trkinfos = fmTrkInfoBPF.at(trackId);
2486  if(!trkinfos.empty())
2487  FillTrackInfoVars(*trkinfos[0], srPng.bpf.pion);
2488  }
2489 
2490  } // end if(PDG == 211)
2491 
2492  // If this is the proton, fill the proton track info...
2493  else if(abs(bpfFitSums[0]->PDG()) == 2212) {
2494 
2495  srPng.bpf.proton.IsValid = true;
2496 
2497  // Fill the Track Vars with overlapping E if that info is
2498  // present, otherwise use the default FillTrackVars
2499  // function.
2500  if (fmBPFOverlapE.isValid()) {
2501  std::vector<art::Ptr<rb::Energy>> overlapEs =
2502  fmBPFOverlapE.at(trackId);
2503  if (!overlapEs.empty()) {
2505  (*bpfTracks[trackId]), (*overlapEs[0]), srPng.bpf.proton,
2506  *slices, allTracks, sliceTracks, allTracksBirks,
2507  sliceTracksBirks, sliceId);
2508  } else {
2509  FillTrackVars((*bpfTracks[trackId]), srPng.bpf.proton,
2510  *slices, allTracks, sliceTracks,
2511  allTracksBirks, sliceTracksBirks,
2512  sliceId);
2513  }
2514  } else {
2515  FillTrackVars((*bpfTracks[trackId]), srPng.bpf.proton, *slices,
2516  allTracks, sliceTracks, allTracksBirks,
2517  sliceTracksBirks, sliceId);
2518  }
2519 
2520  // Fill the FitSum info
2521  FillTrackVarsBpfFitSum((*bpfFitSums[0]), srPng.bpf.proton);
2522  if (fmTrkCntBPF.isValid()) {
2523  // Get track containment info from TrkCntObj
2524  std::vector<art::Ptr<cosrej::TrkCntObj>> trkcnts =
2525  fmTrkCntBPF.at(trackId);
2526 
2527  if (!trkcnts.empty())
2528  FillTrackContainmentVars(*trkcnts[0], srPng.bpf.proton);
2529  }
2530 
2531  if (fmTrkNDRecoBPF.isValid()) {
2532  std::vector<art::Ptr<ndreco::NDRecoBpfObj>> ndrecos =
2533  fmTrkNDRecoBPF.at(trackId);
2534 
2535  if (!ndrecos.empty())
2536  FillNDRecoBPFVars(*ndrecos[0], srPng.bpf.proton);
2537  }
2538 
2539  // Fill the TrackInfo info (yes I said info twice...)
2540  if(fmTrkInfoBPF.isValid()){
2541  std::vector<art::Ptr<trackinfo::TrackInfoObj> > trkinfos = fmTrkInfoBPF.at(trackId);
2542  if(!trkinfos.empty())
2543  FillTrackInfoVars(*trkinfos[0], srPng.bpf.proton);
2544  }
2545 
2546  } // end if(PDG == 2212)
2547 
2548  // otherwise, the PDG code was not one of the three track assumptions
2549  else {
2550  // throw a warning but continue...
2551  std::cout << "\n\nWARNING: BPF FitSum object contained non-standard PDG code "
2552  << bpfFitSums[0]->PDG()
2553  << ". Continuing but not filling info for this track.\n\n";
2554  } // end else
2555 
2556  } // fmBPFFitSums.isValid()
2557  else std::cerr << "fmBPFFitSums is not valid" << std::endl;
2558 
2559  } // loop trackId
2560 
2561  } // fmBPF.isValid()
2562 
2563  pngId++;
2564 
2565  } // end loop over 3D prongs
2566  } // there are non-zero 3D prongs
2567  rec.vtx.elastic.fuzzyk.nshwlid = nshwlids;
2568  rec.vtx.elastic.fuzzyk.fillSizes();
2569 
2570  // energy from prong orphaned hits
2571  art::FindManyP<rb::Prong> fmProngOrph(
2572  elastics, evt,fParams.FuzzyKOrphLabel());
2573  if (fmProngOrph.isValid()) {
2574  std::vector<art::Ptr<rb::Prong>> prongsOrph = fmProngOrph.at(vtxId);
2575  if (prongsOrph.size() > 0) {
2576  rec.vtx.elastic.fuzzyk.orphCalE =
2577  prongsOrph[0]->CalorimetricEnergy();
2578  }
2579  }
2580 
2581  std::sort(srVtx.fuzzyk.png.begin(), srVtx.fuzzyk.png.end(), sortProngCalE); // Use prong calE to sort
2582 
2583  float longest = 0;
2584  for (unsigned int pngIdx = 0; pngIdx < srVtx.fuzzyk.png.size();
2585  pngIdx++) {
2586  if (srVtx.fuzzyk.png[pngIdx].len > longest) {
2587  longest = srVtx.fuzzyk.png[pngIdx].len;
2588  srVtx.fuzzyk.longestidx = pngIdx; // note the longest prong index
2589  }
2590  }
2591 
2592  if (nshwlids >
2593  0) { // fill nue energy for primary shower if there were showers
2594  rec.energy.nue.lid.E =
2595  rec.vtx.elastic.fuzzyk.png[0].shwlid.lidE.E;
2596  rec.energy.nue.lid.depE =
2597  rec.vtx.elastic.fuzzyk.png[0].shwlid.lidE.depE;
2598  rec.energy.nue.lid.shwE =
2599  rec.vtx.elastic.fuzzyk.png[0].shwlid.lidE.shwE;
2600  rec.energy.nue.lid.hadE =
2601  rec.vtx.elastic.fuzzyk.png[0].shwlid.lidE.hadE;
2602  }
2603 
2604  } // end if 3D prongs are valid for this vertex
2605  } // end loop over vertices
2606 
2607  } // end ElasticArms
2608 
2609  FillNueEnergyVars(slice, rec.vtx, rec.energy.nue);
2610 
2611  // HoughVertex
2612 
2613  if (fmHoughVertex.isValid()) {
2614  std::vector<art::Ptr<rb::Vertex>> houghs = fmHoughVertex.at(sliceId);
2615 
2616  for (unsigned int vtxId = 0; vtxId < houghs.size(); ++vtxId) {
2617  rec.vtx.hough.push_back(SRHoughVertex());
2618  SRHoughVertex& srVtx = rec.vtx.hough.back();
2619 
2620  srVtx.time = houghs[vtxId]->GetT();
2621  srVtx.vtx = houghs[vtxId]->GetXYZ();
2622 
2623  /* if(fmFuzzyProng3D.isValid()){
2624  std::vector< art::Ptr<rb::Prong> > prongs = fmFuzzyProng3D.at(vtxId);
2625  // Sort prongs in descending length order
2626  std::sort(prongs.begin(), prongs.end(),sortRBProngLength);
2627  for(art::Ptr<rb::Prong> prong: prongs){
2628  if(prong->NXCell() == 0 || prong->NYCell() == 0){
2629  srVtx.fuzzyk.png2d.push_back(SRProng());
2630  FillProngVars(*prong, srVtx.fuzzyk.png2d.back(), slicess,
2631  sliceId);
2632  }// end 2d prong case
2633  else{
2634  srVtx.fuzzyk.png.push_back(SRProng());
2635  FillProngVars(*prong, srVtx.fuzzyk.png.back(), *slices, sliceId);
2636  }// end 3d prong case
2637  }// end loop over prongs
2638  }// end if prongs are valid for this vertex
2639  // Check if there is a valid instance label for 2D prongs
2640  if(fmFuzzyProng2D.isValid() && fFuzzyK2dAssnLabel != ""){
2641  std::vector< art::Ptr<rb::Prong> > prongs = fmFuzzyProng2D.at(vtxId);
2642  // Sort prongs in descending length order
2643  std::sort(prongs.begin(), prongs.end(),sortRBProngLength);
2644  for(art::Ptr<rb::Prong> prong: prongs){
2645  if(prong->NXCell() == 0 || prong->NYCell() == 0){
2646  srVtx.fuzzyk.png2d.push_back(SRProng());
2647  FillProngVars(*prong, srVtx.fuzzyk.png2d.back(), *slices,
2648  sliceId);
2649  }// end 2d prong case
2650  else{
2651  srVtx.fuzzyk.png.push_back(SRProng());
2652  FillProngVars(*prong, srVtx.fuzzyk.png.back(), *slices, sliceId);
2653  }// end 3d prong case
2654  }// end loop over prongs
2655  }// end if prongs are valid for this vertex
2656  */
2657  } // end for vtxId
2658  // Sort the vector
2659  // std::sort(rec.vtx.hough.begin(), rec.vtx.hough.end(),
2660  // sortHoughProng);
2661  }
2662 
2663  // Get all the VertexDT vertices associated with the slice
2664  if (fmVertexDT.isValid()) {
2665  std::vector<art::Ptr<rb::Vertex>> vtxs = fmVertexDT.at(sliceId);
2666  for (unsigned int vtxId = 0; vtxId < vtxs.size(); ++vtxId) {
2667  rec.vtx.vdt.push_back(SRVertexDT());
2668  SRVertexDT& srVtx = rec.vtx.vdt.back();
2669  srVtx.time = vtxs[vtxId]->GetT();
2670  srVtx.vtx = vtxs[vtxId]->GetXYZ();
2671  }
2672  }
2673  rec.vtx.fillSizes();
2674 
2675  // Get all MEFs aassociated with slice
2676  if (fmSlcME.isValid()) {
2677  std::vector<art::Ptr<me::SlcME>> michels = fmSlcME.at(sliceId);
2678  for (art::Ptr<me::SlcME> michel : michels) {
2679  AddSlcMEToVec(*michel, &rec.me.slc, *slices, allTracks, sliceTracks,
2680  allTracksBirks, sliceTracksBirks, sliceId);
2681  }
2682  }
2683  rec.me.fillSizes();
2684 
2685  // Get slice level PID from LIDBuilder
2686  if (foLID.isValid()) {
2687  cet::maybe_ref<art::Ptr<slid::EventLID> const> rlid(foLID.at(sliceId));
2688  art::Ptr<slid::EventLID> elid = rlid.ref();
2689  // If the pointer is valid
2690  if (elid) {
2691  rec.sel.lid.ann = elid->Value();
2692  rec.sel.lid.anne = elid->ValueWithE();
2693  rec.sel.lid.annepi0 = elid->ValueEPi0();
2694  rec.sel.lid.annecos = elid->ValueECos();
2695  } else {
2696  rec.sel.lid.setDefault();
2697  }
2698  }
2699 
2700  if (foSlcLID.isValid()) {
2702  foSlcLID.at(sliceId));
2703  art::Ptr<slid::SliceLID> slclid = rslclid.ref();
2704  if (slclid) {
2705  rec.sel.lid.rnn = slclid->Value();
2706  } else {
2707  rec.sel.lid.setDefault();
2708  }
2709  }
2710 
2711  /*
2712  if(fmShwLID.isValid()){
2713  std::vector<const slid::ShowerLID* > shwlids = fmShwLID.at(sliceId);
2714  // Sort here since the vector may not be
2715  sort(shwlids.begin(),shwlids.end(), slid::CompareByEnergy);
2716  art::FindOneP<rb::Shower> fos(shwlids, evt, fLIDLabel());
2717  for(unsigned int shwId = 0; shwId < shwlids.size(); ++shwId){
2718  cet::maybe_ref<art::Ptr<rb::Shower> const> rshw(fos.at(shwId));
2719  art::Ptr<rb::Shower> shw = rshw.ref();
2720  if (shw){
2721  rec.shw.shwlid.push_back(SRShower());
2722  FillShowerVars(*shw, rec.shw.shwlid.back(), *slices, sliceId);
2723  }
2724 
2725  rec.sel.elecid.shwlid.push_back(SRSLid());
2726  FillSlidVars(*shwlids[shwId], rec.sel.elecid.shwlid.back().shwlid);
2727 
2728  // Fill energy information for leading energy shower
2729  // which also serves as a slice pid
2730  if (shwId == 0){
2731  FillLIDEnergyVars(*shwlids[shwId], rec.energy.elecid);
2732  }
2733 
2734  // Fill energy vector with the other shower energies
2735  rec.energy.shwlid.push_back(SRSLidEnergy());
2736  FillLIDEnergyVars(*shwlids[shwId], rec.energy.shwlid.back());
2737 
2738  if (shw && shwId == 0){
2739  FillNueEnergyVars(slice, *shw, rec.energy.nue);
2740  }
2741  } // end for loop on showerlids
2742  } // end if showerlids are valid
2743 
2744  if(fmShwPID.isValid()){
2745  std::vector<const slid::ShowerPID* > shwpids = fmShwPID.at(sliceId);
2746  // Sort here since the vector may not be
2747  sort(shwpids.begin(),shwpids.end(), slid::CompareByEnergy);
2748  art::FindOneP<rb::Shower> fos(shwpids, evt, fSPIDLabel());
2749  for(unsigned int shwId = 0; shwId < shwpids.size(); ++shwId){
2750  rec.sel.elecid.shwpid.push_back(SRSPid());
2751  FillSpidVars(*shwpids[shwId], rec.sel.elecid.shwpid.back());
2752  } // end for loop on showerpids
2753  } // end if showerpids are valid
2754 
2755  } // end if event id is valid
2756  else{
2757  rec.sel.lid .setDefault();
2758  rec.energy.nue.lid .setDefault();
2759  }
2760 } // end if slice lids are valid
2761 else{
2762  rec.sel.lid .setDefault();
2763  rec.energy.lid .setDefault();
2764 }
2765 
2766 rec.shw.fillSizes();
2767 rec.energy.fillSizes();
2768 rec.sel.elecid.fillSizes();
2769  */
2771  bool isLemmed = GetAssociatedProduct(fmLem, sliceId, lem);
2772  std::vector<std::string> vec = {fParams.LEMNuePresel()};
2773 
2774  if (isLemmed && !(fParams.ApplyLEMNuePresel() &&
2775  rb::IsFiltered(evt, slices, sliceId, vec))) {
2776  FillLEMVars(lem, rec.sel.lem);
2777  } else
2778  rec.sel.lem.setDefault();
2779 
2780  rvp::RVP rvp;
2781  if (GetAssociatedProduct(fmRvp, sliceId, rvp))
2782  FillRVPVars(rvp, rec.sel.rvp);
2783  else
2784  rec.sel.rvp.setDefault();
2785 
2786  xnue::Xnue xnue;
2787  if (GetAssociatedProduct(fmXnue, sliceId, xnue)) {
2788  FillXnueVars(xnue, rec.sel.xnuepid);
2789  } else
2790  rec.sel.xnuepid.setDefault();
2791 
2792  numue::NumuE numuE;
2793  if (GetAssociatedProduct(fmNumuEnergy, sliceId, numuE)) {
2794  FillNumuEnergyVars(numuE, rec.energy.numu);
2795  } else {
2796  rec.energy.numu.setDefault();
2797  }
2798 
2799 
2800  cvn::RegNuResult evtregcvn;
2801  if (GetAssociatedProduct(fmEvtRegCVN, sliceId, evtregcvn)) {
2802  rec.energy.nue.regcvnEvtE = evtregcvn.Output();
2803  } else {
2804  rec.energy.nue.regcvnEvtE = -5.;
2805  }
2806 
2807  cvn::RegNuResult evtregcvn_opphorn;
2808  if (GetAssociatedProduct(fmEvtRegCVN_opphorn, sliceId, evtregcvn_opphorn)) {
2809  rec.energy.nue.regcvnEvtE_opphorn = evtregcvn.Output();
2810  } else {
2811  rec.energy.nue.regcvnEvtE_opphorn = -5.;
2812  }
2813 
2814 
2815 
2816  LSTME::LSTMEnergy numulstmE;
2817  if (GetAssociatedProduct(fmNumuLSTMEnergy, sliceId, numulstmE)) {
2818  FillNumuLSTMEnergyVars(numulstmE, rec.energy.numu.lstmmuon, rec.energy.numu.lstmnu);
2819  } else {
2820  rec.energy.numu.setLSTMDefault();
2821  }
2822 
2823  LSTME::LSTMEnergy numulstmE_opphorn;
2824  if (GetAssociatedProduct(fmNumuLSTMEnergy_opphorn, sliceId, numulstmE_opphorn)) {
2825  FillNumuLSTMEnergyVars(numulstmE_opphorn, rec.energy.numu.lstmmuon_opphorn, rec.energy.numu.lstmnu_opphorn);
2826  } else {
2827  rec.energy.numu.setLSTMDefault(true);
2828  }
2829 
2830  cvn::RegHadronResult hadregcvn;
2831  if (GetAssociatedProduct(fmHadRegCVN, sliceId, hadregcvn)) {
2832  rec.energy.numu.regcvnhadE = hadregcvn.Output();
2833  } else {
2834  rec.energy.numu.regcvnhadE = -5.;
2835  }
2836 
2837  cvn::RegHadronResult hadregcvn_opphorn;
2838  if (GetAssociatedProduct(fmHadRegCVN_opphorn, sliceId, hadregcvn_opphorn)) {
2839  rec.energy.numu.regcvnhadE_opphorn = hadregcvn.Output();
2840  } else {
2841  rec.energy.numu.regcvnhadE_opphorn = -5.;
2842  }
2843 
2844  SliceLID::Prediction sliceLID;
2845  if (GetAssociatedProduct(fmSliceLID, sliceId, sliceLID)) {
2846  FillSliceLID(sliceLID, rec.sel.slicelid);
2847  } else {
2848  rec.sel.slicelid.setDefault();
2849  }
2850 
2851  SliceLID::Prediction sliceLID_opphorn;
2852  if (GetAssociatedProduct(fmSliceLID_opphorn, sliceId, sliceLID_opphorn)) {
2853  FillSliceLID(sliceLID_opphorn, rec.sel.slicelid_opphorn);
2854  } else {
2855  rec.sel.slicelid_opphorn.setDefault();
2856  }
2857 
2859  if (GetAssociatedProduct(fmCosRej, sliceId, cosrej)) {
2860  FillCosRejVars(cosrej, rec.sel.cosrej, rec.sel.contain);
2861 
2862  float kDistAllTop = rec.sel.nuecosrej.distallpngtop;
2863  if (std::isnan(kDistAllTop)) kDistAllTop = -1000.0;
2864 
2865  float kDistAllBottom = rec.sel.nuecosrej.distallpngbottom;
2866  if (std::isnan(kDistAllBottom)) kDistAllBottom = -1000.0;
2867 
2868  float kDistAllEast = rec.sel.nuecosrej.distallpngeast;
2869  if (std::isnan(kDistAllEast)) kDistAllEast = -1000.0;
2870 
2871  float kDistAllWest = rec.sel.nuecosrej.distallpngwest;
2872  if (std::isnan(kDistAllWest)) kDistAllWest = -1000.0;
2873 
2874  float kDistAllBack = rec.sel.nuecosrej.distallpngback;
2875  if (std::isnan(kDistAllBack)) kDistAllBack = -1000.0;
2876 
2877  float kDistAllFront = rec.sel.nuecosrej.distallpngfront;
2878  if (std::isnan(kDistAllFront)) kDistAllFront = -1000.0;
2879 
2880  if (fDet ==
2881  kFARDET) { // adopted from CAFAna/Cuts/NumuCuts2017.h on 8.16.17
2882  std::pair<int, int> planes = calcFirstLastLivePlane(
2883  rec.slc.firstplane, std::bitset<14>(rec.hdr.dibmask));
2884  int planestofront = rec.slc.firstplane - planes.first;
2885  int planestoback = planes.second - rec.slc.lastplane;
2886  bool kcut1 =
2887  (rec.sel.contain.kalfwdcell > 6 && rec.sel.contain.kalbakcell > 6 &&
2888  rec.sel.contain.cosfwdcell > 0 && rec.sel.contain.cosbakcell > 7 &&
2889  planestofront > 1 && planestoback > 1);
2890  bool kcut2 = kDistAllTop > 60 && kDistAllBottom > 12 &&
2891  kDistAllEast > 16 && kDistAllWest > 12 &&
2892  kDistAllFront > 18 && kDistAllBack > 18;
2893  rec.sel.contain.numucontain = (kcut1 && kcut2);
2894 
2895  if (cosrej.CosFwdCell() > 0 && cosrej.CosBakCell() > 0 &&
2896  cosrej.KalFwdCell() > 10 && cosrej.KalBakCell() > 10 &&
2897  rec.slc.ncellsfromedge > 1 && rec.sel.contain.nplanestofront > 1 &&
2898  rec.sel.contain.nplanestoback > 1) {
2899  rec.sel.contain.numucontainSA = true;
2900  } else
2901  rec.sel.contain.numucontainSA = false;
2902  }
2903 
2904  else if (fDet == kNEARDET) { // adopted from CAFAna/Cuts/NumuCuts2017.h
2905  // on 8.16.17
2906  if (rec.vtx.elastic.IsValid < 1)
2907  rec.sel.contain.numucontain = false;
2908  else {
2909  for (unsigned int i = 0; i < rec.vtx.elastic.fuzzyk.nshwlid; ++i) {
2910  TVector3 start = rec.vtx.elastic.fuzzyk.png[i].shwlid.start;
2911  TVector3 stop = rec.vtx.elastic.fuzzyk.png[i].shwlid.stop;
2912  if (std::min(start.X(), stop.X()) < -180.0)
2913  rec.sel.contain.numucontain = false;
2914  if (std::max(start.X(), stop.X()) > 180.0)
2915  rec.sel.contain.numucontain = false;
2916  if (std::min(start.Y(), stop.Y()) < -180.0)
2917  rec.sel.contain.numucontain = false;
2918  if (std::max(start.Y(), stop.Y()) > 180.0)
2919  rec.sel.contain.numucontain = false;
2920  if (std::min(start.Z(), stop.Z()) < 20.0)
2921  rec.sel.contain.numucontain = false;
2922  if (std::max(start.Z(), stop.Z()) > 1525.0)
2923  rec.sel.contain.numucontain = false;
2924  }
2925  if (rec.trk.kalman.ntracks < 1) rec.sel.contain.numucontain = false;
2926  for (unsigned int i = 0; i < rec.trk.kalman.ntracks; ++i) {
2927  if (i == rec.trk.kalman.idxremid)
2928  continue;
2929  else if (rec.trk.kalman.tracks[i].start.Z() > 1275 ||
2930  rec.trk.kalman.tracks[i].stop.Z() > 1275)
2931  rec.sel.contain.numucontain = false;
2932  }
2933  }
2934 
2935  rec.sel.contain.numucontain =
2936  (rec.trk.kalman.ntracks > rec.trk.kalman.idxremid &&
2937  rec.slc.firstplane > 1 && rec.slc.lastplane < 212 &&
2938  rec.trk.kalman.tracks[0].start.Z() < 1100 &&
2939  (rec.trk.kalman.tracks[0].stop.Z() < 1275 ||
2940  rec.sel.contain.kalyposattrans < 55) &&
2941  rec.sel.contain.kalfwdcellnd > 5 &&
2942  rec.sel.contain.kalbakcellnd > 10);
2943 
2944  if (cosrej.KalFwdCellND() > 4 && cosrej.KalBakCellND() > 8 &&
2945  rec.slc.ncellsfromedge > 1 && rec.slc.firstplane > 1 &&
2946  rec.slc.lastplane < 212 &&
2947  rec.trk.kalman.tracks[0].start.Z() < 1150 &&
2948  (rec.trk.kalman.tracks[0].stop.Z() < 1275 ||
2949  rec.sel.contain.kalyposattrans < 55) &&
2950  (rec.energy.numu.ndhadcaltranE + rec.energy.numu.ndhadcalcatE) <
2951  0.03) {
2952  rec.sel.contain.numucontainSA = true;
2953  } else
2954  rec.sel.contain.numucontainSA = false;
2955  }
2956  }
2957 
2958  cosrej::NueCosRej nuecosrej;
2959  if (GetAssociatedProduct(fmNueCosRej, sliceId, nuecosrej)) {
2960  FillNueCosRejVars(nuecosrej, rec.sel.nuecosrej);
2961  /*
2962  // nue containment boolean - nue group please feel free to update.
2963  // If variables move from sandbox this needs to be changed too
2964  if(fDet == kFARDET){
2965  if(nuecosrej.ProngMaxX() < 750 && nuecosrej.ProngMinX() > -750 &&
2966  nuecosrej.ProngMaxY() < 750 && nuecosrej.ProngMinY() > -750 &&
2967  std::min(nuecosrej.StartDistToFront(),
2968  nuecosrej.StopDistToFront()) > 50 &&
2969  std::min(nuecosrej.StartDistToBack(),
2970  nuecosrej.StopDistToBack()) > 50){
2971  rec.sel.contain.nuecontain = true;
2972  }
2973  }
2974 
2975  if(fDet == kNEARDET){
2976  if(nuecosrej.ProngMaxX() < 185 && nuecosrej.ProngMinX() > -185 &&
2977  nuecosrej.ProngMaxY() < 185 && nuecosrej.ProngMinY() > -185 &&
2978  std::min(nuecosrej.StartDistToFront(),
2979  nuecosrej.StopDistToFront()) > 50 &&
2980  std::min(nuecosrej.StartDistToBack(),
2981  nuecosrej.StopDistToBack()) > 50){
2982  rec.sel.contain.nuecontain = true;
2983  }
2984  }
2985  */
2986  } else {
2987  rec.sel.nuecosrej.setDefault();
2988  }
2989 
2990  ncid::NCCosRej nccosrej; ///////////// NC Cosmic Rejection
2991  if (GetAssociatedProduct(fmNCCosRej, sliceId, nccosrej)) {
2992  FillNCCosRejVars(nccosrej, rec.sel.nccosrej);
2993  } else {
2994  rec.sel.nccosrej.setDefault();
2995  }
2996 
2997  ncpi0::NCPi0BkgRej ncpi0bkgrej;
2998  if (GetAssociatedProduct(fmNCPi0BkgRej, sliceId, ncpi0bkgrej)) {
2999  FillNCPi0BkgRejVars(ncpi0bkgrej, rec.sel.ncpi0bkgrej);
3000  } else {
3001  rec.sel.ncpi0bkgrej.setDefault();
3002  }
3003 
3004  presel::PreselObj nuepre;
3005  if (GetAssociatedProduct(fmNuePresel, sliceId, nuepre)) {
3006  FillNuePreselVars(nuepre, rec.sel.nuepre);
3007  } else {
3008  rec.sel.nuepre.setDefault();
3009  }
3010 
3011  presel::PreselObj rockpre;
3012  if (GetAssociatedProduct(fmRockPresel, sliceId, rockpre)) {
3013  FillRockPreselVars(rockpre, rec.sel.rockpre);
3014  } else {
3015  rec.sel.rockpre.setDefault();
3016  }
3017 
3018  presel::Veto veto, nueveto;
3019  if (GetAssociatedProduct(fmVeto, sliceId, veto) ||
3020  GetAssociatedProduct(fmNueVeto, sliceId, nueveto)) {
3021  FillVetoVars(veto, nueveto, rec.sel.veto);
3022  } else {
3023  rec.sel.veto.setDefault();
3024  }
3025 
3026  cvn::Result nuoneResult;
3027  if (GetAssociatedProduct(fmNuonE, sliceId, nuoneResult)) {
3028  FillNuonEResultVars(nuoneResult, rec.sel.nuone);
3029  } else {
3030  rec.sel.nuone.setDefault();
3031  }
3032 
3033  cvn::Result nuoneResult_opphorn;
3034  if (GetAssociatedProduct(fmNuonE_opphorn, sliceId, nuoneResult_opphorn)) {
3035  FillNuonEResultVars(nuoneResult_opphorn, rec.sel.nuone_opphorn);
3036  } else {
3037  rec.sel.nuone_opphorn.setDefault();
3038  }
3039 
3040  // --- Horrible hack to appease CAFAna.
3041  rec.sel.cvn.setDefault();
3042 
3043  // --- Loose Presel PtP Cut
3044  cvn::Result cvnResult_LoosePreselPtP;
3045  if (GetAssociatedProduct(fmCVN_LoosePreselPtP, sliceId, cvnResult_LoosePreselPtP)) {
3046  FillCVNResultVars(cvnResult_LoosePreselPtP, rec.sel.cvnloosepreselptp);
3047  } else {
3048  rec.sel.cvnloosepreselptp.setDefault();
3049  }
3050  cvn::Result cvnResult_LoosePreselPtP_opphorn;
3051  if (GetAssociatedProduct(fmCVN_LoosePreselPtP_opphorn, sliceId, cvnResult_LoosePreselPtP_opphorn)) {
3052  FillCVNResultVars(cvnResult_LoosePreselPtP_opphorn, rec.sel.cvnloosepreselptp_opphorn);
3053  } else {
3054  rec.sel.cvnloosepreselptp_opphorn.setDefault();
3055  }
3056  // --- Old Presel
3057  cvn::Result cvnResult_OldPresel;
3058  if (GetAssociatedProduct(fmCVN_OldPresel, sliceId, cvnResult_OldPresel)) {
3059  FillCVNResultVars(cvnResult_OldPresel, rec.sel.cvnoldpresel);
3060  } else {
3061  rec.sel.cvnoldpresel.setDefault();
3062  }
3063  cvn::Result cvnResult_OldPresel_opphorn;
3064  if (GetAssociatedProduct(fmCVN_OldPresel_opphorn, sliceId, cvnResult_OldPresel_opphorn)) {
3065  FillCVNResultVars(cvnResult_OldPresel_opphorn, rec.sel.cvnoldpresel_opphorn);
3066  } else {
3067  rec.sel.cvnoldpresel_opphorn.setDefault();
3068  }
3069  // --- No Cosmics used in training
3070  cvn::Result cvnResult_NoCosmics;
3071  if (GetAssociatedProduct(fmCVN_NoCosmics, sliceId, cvnResult_NoCosmics)) {
3072  FillCVNResultVars(cvnResult_NoCosmics, rec.sel.cvnnocosmics);
3073  } else {
3074  rec.sel.cvnnocosmics.setDefault();
3075  }
3076  cvn::Result cvnResult_NoCosmics_opphorn;
3077  if (GetAssociatedProduct(fmCVN_NoCosmics_opphorn, sliceId, cvnResult_NoCosmics_opphorn)) {
3078  FillCVNResultVars(cvnResult_NoCosmics_opphorn, rec.sel.cvnnocosmics_opphorn);
3079  } else {
3080  rec.sel.cvnnocosmics_opphorn.setDefault();
3081  }
3082  // --- Wrong Sign
3083  cvn::Result cvnResult_ws;
3084  if (GetAssociatedProduct(fmWS, sliceId, cvnResult_ws)) {
3085  rec.sel.wsid = cvnResult_ws.fOutput[0];
3086  }
3087 
3088  //cvn::Features cvnFeatures;
3089  //if (GetAssociatedProduct(fmCVNFeatures, sliceId, cvnFeatures) &&
3090  // fParams.FullTrainingInfo()) {
3091  // FillCVNFeaturesVars(cvnFeatures, rec.training.cvnfeatures,
3092  // fParams.NumCVNPrincipalComponents());
3093  //}
3094 
3095  if (fmCVNPixelMaps.isValid()){
3096  std::vector<art::Ptr<cvn::PixelMap>> cvnPixelMaps =
3097  fmCVNPixelMaps.at(sliceId);
3098  if (cvnPixelMaps.size())
3099  rec.training.cvnmapsfilled = true;
3100  else
3101  rec.training.cvnmapsfilled = false;
3102  if (fParams.FillPixelMaps()) {
3103  for (unsigned int ip = 0; ip < cvnPixelMaps.size(); ++ip) {
3104  rec.training.cvnmaps.push_back(SRPixelObjMap());
3105  SRPixelObjMap& pmap = rec.training.cvnmaps.back();
3106  FillCVNPixelMaps(*cvnPixelMaps[ip], pmap, fParams.UseGeVPixelMaps());
3107  }
3108  }
3109  } else {
3110  rec.training.cvnmapsfilled = false;
3111  }
3112 
3113  cvn::TrainingData cvnTrainingData;
3114  if (GetAssociatedProduct(fmCVNTrainingData, sliceId, cvnTrainingData) &&
3116  rec.training.trainingdata.push_back(SRTrainingData());
3117  SRTrainingData& srTrainingData = rec.training.trainingdata.back();
3118  FillCVNTrainingData(cvnTrainingData, srTrainingData);
3119  }
3120 
3121  // also store every hit in the slice. Doing this in a fixed size array the
3122  // size of the detector since the zeros compress better
3123  if (fParams.FillPixelMaps()) {
3124  rec.training.slicemaps.push_back(SRSliceMap());
3125  SRSliceMap& smap = rec.training.slicemaps.back();
3126  float corr = 15123.7; // Rough PE to GeV? Rev 17083
3127  float peCorrChunk;
3128  // This value is 2.0 / (PE/GeV)
3129  if (fParams.UseGeVPixelMaps())
3130  peCorrChunk = 1.322e-4 / 255.0;
3131  else
3132  peCorrChunk = 2.0 / 255.0;
3133  for (unsigned int ih = 0; ih < slice.NCell(); ih++) {
3134  rb::RecoHit rh = slice.RecoHit(ih);
3135  art::Ptr<rb::CellHit> ch = slice.Cell(ih);
3136  // convert planes to 0 to 447 in each view
3137  unsigned int pl = (ch->Plane() / 2);
3138  unsigned int chan = ch->Cell() + 384 * (pl + 448 * ch->View());
3139  // convert cell energy to an unsigned char. Scales same as in CVN
3140  float hitPE = ch->PE();
3141  float gev = (rh.IsCalibrated()) ? rh.GeV() : hitPE / corr;
3142  float pe = (rh.IsCalibrated()) ? rh.PECorr() : hitPE;
3143  float tv = (fParams.UseGeVPixelMaps()) ? gev : pe;
3144  tv /= 2100.0; // scale factor from CVN Pixel maps
3145  float truncateCorr = ceil(tv / peCorrChunk);
3146  unsigned char value;
3147  if (!fParams.UseEnergy()) {
3148  if (truncateCorr > 0)
3149  truncateCorr = 1;
3150  else
3151  truncateCorr = 0;
3152  }
3153  if (truncateCorr > 255)
3154  value = 255;
3155  else
3156  value = (unsigned char)truncateCorr;
3157  smap.slicemap[chan] = value;
3158  }
3159  }
3160 
3162  if (GetAssociatedProduct(fmMrccPar, sliceId, parent))
3163  FillMRCCParentInfo(parent, rec.parent.mrccpar, evt);
3164  else
3165  rec.parent.mrccpar.setDefault();
3166 
3167  //#######################################################
3168  // Fill truth information
3169  //#######################################################
3170 
3171  // Set mc branch values to default
3172  rec.mc.setDefault();
3173  if (bt->HaveTruthInfo()) {
3174  // Get all truth track hits associated with slice
3175 
3176  int currentNuIdx = 0;
3177  unsigned int nNuGlobal =
3178  0; // count of neutrino interactions in this slice
3179 
3180  for (unsigned int iEffPur = 0; iEffPur < sliceTruthTable[sliceId].size();
3181  ++iEffPur) {
3182  // continue if truth left no hits in slice
3183  if (!(sliceTruthTable[sliceId][iEffPur].nSliceHits > 0)) continue;
3184 
3185  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
3186  faveIdsEnergy[sliceId]) {
3187  cheat::NeutrinoEffPur* faveNu = &sliceTruthTable[sliceId][iEffPur];
3188 
3189  // If NOT Neutrinos
3190  if (!faveNu->neutrinoInt->NeutrinoSet()) {
3191  // Add faveNu to allcosmics
3192  if (fParams.FullTruthInfo())
3193  AddCosmicTruthToVec(faveNu, allTracks, sliceTracks,
3194  allTracksBirks, sliceTracksBirks,
3195  &rec.mc.allcosmics);
3196 
3197  // Add faveNu to cosmic
3198  AddCosmicTruthToVec(faveNu, allTracks, sliceTracks, allTracksBirks,
3199  sliceTracksBirks, &rec.mc.cosmic);
3200  }
3201  // Otherwise we are a neutrino
3202  else {
3203  // Add faveNu to allnus
3204  if (fParams.FullTruthInfo())
3205  AddMCTruthToVec(faveNu, allTracks, sliceTracks, allTracksBirks,
3206  sliceTracksBirks, &rec.mc.allnus, evt, spill);
3207 
3208  // Add faveNu to nu
3209  AddMCTruthToVec(faveNu, allTracks, sliceTracks, allTracksBirks,
3210  sliceTracksBirks, &rec.mc.nu, evt, spill);
3211  rec.mc.faveidxenergy = currentNuIdx;
3212 
3214  pv.push_back(faveNu->neutrinoInt);
3215 
3216  if (rec.mc.nu.size() > 0) {
3217  // if a neutrino is matched to the slice, fill the mc reweighting
3218  // object.
3219 
3221  pv, evt, fParams.ReweightLabel());
3223  if (GetAssociatedProduct(fmRW, 0, rw))
3224  rec.mc.nu[0].rwgt.genie = GenieReweightTable(rw);
3225 
3227  pv, evt, fParams.FluxReweightLabel());
3228  fxwgt::FluxWeights flxrw;
3229  if (GetAssociatedProduct(fmFLXRW, 0, flxrw))
3230  rec.mc.nu[0].rwgt.ppfx = FluxReweights(flxrw);
3231  else
3232  rec.mc.nu[0].rwgt.ppfx.setDefault();
3233 
3234 
3235 
3237  pv, evt, fParams.G4ReweightLabel());
3239  if (GetAssociatedProduct(fmG4RW, 0, g4w))
3240  rec.mc.nu[0].rwgt.geant4 = Geant4Reweights(g4w);
3241  else
3242  rec.mc.nu[0].rwgt.geant4.setDefault();
3243 
3244 
3245 
3246  }
3247 
3248 
3249  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
3250  faveIdsEff[sliceId])
3251  rec.mc.faveidxeff = currentNuIdx;
3252  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
3253  faveIdsPur[sliceId])
3254  rec.mc.faveidxpur = currentNuIdx;
3255  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
3256  faveIdsEffPur[sliceId])
3257  rec.mc.faveidxeffpur = currentNuIdx;
3258  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
3259  faveIdsEffThenPur[sliceId])
3260  rec.mc.faveidxeffthenpur = currentNuIdx;
3261 
3262  currentNuIdx++;
3263  }
3264  }
3265  }
3266 
3267  for (unsigned int iEffPur = 0; iEffPur < sliceTruthTable[sliceId].size();
3268  ++iEffPur) {
3269  if (sliceTruthTable[sliceId][iEffPur].neutrinoInt &&
3270  sliceTruthTable[sliceId][iEffPur].neutrinoInt->NeutrinoSet())
3271  nNuGlobal++;
3272  // continue if truth left no hits in slice
3273  if (!(sliceTruthTable[sliceId][iEffPur].nSliceHits > 0)) continue;
3274 
3275  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex !=
3276  faveIdsEnergy[sliceId]) {
3277  // Cosmics
3278  if (!sliceTruthTable[sliceId][iEffPur].neutrinoInt->NeutrinoSet()) {
3279  if (fParams.FullTruthInfo())
3280  AddCosmicTruthToVec(&sliceTruthTable[sliceId][iEffPur], allTracks,
3281  sliceTracks, allTracksBirks, sliceTracksBirks,
3282  &rec.mc.allcosmics);
3283  }
3284  // Neutrino
3285  else {
3286  if (fParams.FullTruthInfo())
3287  AddMCTruthToVec(&sliceTruthTable[sliceId][iEffPur], allTracks,
3288  sliceTracks, allTracksBirks, sliceTracksBirks,
3289  &rec.mc.allnus, evt, spill);
3290 
3291  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
3292  faveIdsEff[sliceId])
3293  rec.mc.faveidxeff = currentNuIdx;
3294  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
3295  faveIdsPur[sliceId])
3296  rec.mc.faveidxpur = currentNuIdx;
3297  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
3298  faveIdsEffPur[sliceId])
3299  rec.mc.faveidxeffpur = currentNuIdx;
3300  if ((int)sliceTruthTable[sliceId][iEffPur].truthColIndex ==
3301  faveIdsEffThenPur[sliceId])
3302  rec.mc.faveidxeffthenpur = currentNuIdx;
3303 
3304  currentNuIdx++;
3305  }
3306  }
3307  } // end loop to fill all* vectors
3308 
3309  // fill size fields for mc branch
3310  rec.mc.nnu = rec.mc.nu.size();
3311  rec.mc.ncosmic = rec.mc.cosmic.size();
3312  if (fParams.FullTruthInfo()) {
3313  rec.mc.nallnus = rec.mc.allnus.size();
3314  rec.mc.nallcosmics = rec.mc.allcosmics.size();
3315  }
3316  rec.mc.global.nnu = nNuGlobal;
3317  } // end MC block
3318 
3320 
3321  fRecTree->Fill();
3322  srcol->push_back(rec);
3323  util::CreateAssn(*this, evt, *srcol, art::Ptr<rb::Cluster>(slices, sliceId),
3324  *srAssn);
3325  } // end loop over slices
3326 
3327  *spillcol = spill;
3328 
3329  evt.put(std::move(srcol));
3330  evt.put(std::move(spillcol));
3331  if (fParams.FillNuTree()) evt.put(std::move(spilltruthcol));
3332  evt.put(std::move(srAssn));
3333 }
art::ServiceHandle< geo::LiveGeometry > livegeom
Near Detector underground.
Definition: SREnums.h:10
Atom< string > NumuEnergyLabel
Atom< string > DiscreteTrackLabel
T max(const caf::Proxy< T > &a, T b)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
Atom< string > MuonIDLabel
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
static MetadataManager & getInstance()
art::ServiceHandle< cheat::BackTracker > bt
Atom< InputTag > FuzzyK2DLabel
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
void FillCVNNeutronDaughterResultVars(const std::vector< art::Ptr< rb::PID > > &cvnneutrons, caf::SRCVNNeutronDaughterResult &cvnneutron)
Definition: FillPIDs.cxx:95
double HalfD() const
Definition: CellGeo.cxx:205
Atom< string > ClusterLabel
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
Far Detector at Ash River.
Definition: SREnums.h:11
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)
Atom< InputTag > OverlapEKalLabel
Atom< string > DiFShowerCVNLabelOldPresel
const Var kDistAllBottom([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngbottom)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngbottom);})
Distance of all showers in slice from the bottom edge of detector.
Definition: NueVars.h:33
Atom< string > NueSelLabel
art::ServiceHandle< geo::Geometry > geom
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
int NAnalysisChannels(int sr=0)
Definition: RunHistory.cxx:285
Neutrino energy output from RegCVN.
Definition: RegResult.h:29
Atom< string > RockPreselLabel
const Var kDistAllWest([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngwest)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngwest);})
Distance of all showers in slice from the west edge of detector.
Definition: NueVars.h:36
void FillVetoVars(const presel::Veto &veto, const presel::Veto &nueveto, caf::SRVeto &srveto)
Definition: FillPIDs.cxx:554
Atom< string > CosmicCVNLabel
Atom< bool > UseGeVPixelMaps
Atom< string > CVNParticleLabel
Atom< string > NueVetoLabel
X or Y views.
Definition: PlaneGeo.h:30
const Var kDistAllTop([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngtop)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngtop);})
Distance of all showers in slice from the top edge of detector.
Definition: NueVars.h:30
Atom< string > G4ReweightLabel
Atom< string > VetoLabel
unsigned short Plane() const
Definition: CellHit.h:39
Atom< string > KalmanTrackLabel
Atom< bool > FullTruthInfo
float ncVal
Cosmic CVN nc score for each time slice.
Atom< string > BreakPointTrackLabel
void FillMuonIDVars(const muonid::MuonID &muid, SRKalmanTrack &srTrk)
Definition: FillPIDs.cxx:48
void GetByLabelStrict(const art::Event &evt, const std::string &label, art::Handle< T > &handle) const
geo::View_t View() const
Definition: CellHit.h:41
Atom< string > DiFShowerLIDLabel
void FillTrackVars(const rb::Track &trk, SRTrack &srTrk, const std::vector< rb::Cluster > &sliceList, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, int sliceIdx)
Definition: FillReco.cxx:281
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const Var kDistAllBack([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngback)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngback);})
Distance of all showers in slice from the back edge of detector.
Definition: NueVars.h:42
bool DetFineTimingSetting() const
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
Atom< string > SlcLIDLabel
std::vector< NeutrinoWithIndex > allMCTruth() const
static bool compMuonID(const SRKalmanTrack &a, const SRKalmanTrack &b)
void FillTrackContainmentVars(const cosrej::TrkCntObj &trkcnt, SRTrack &srTrk)
Function to fill SRTrack containment information.
Definition: FillReco.cxx:555
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr std::uint32_t timeHigh() const
Definition: Timestamp.h:35
void CopyMuonIDVars(const SRKalmanTrack &srTrk, caf::SRMuonID &muid)
Definition: FillPIDs.cxx:723
Atom< bool > ApplyLEMNuePresel
float Output() const
Definition: RegResult.h:22
float cosmicVal
Cosmic CVN cosmic score for each time slice.
SRMCReweight rwgt
Definition: SRNeutrino.h:91
Atom< string > CVNLabelNoCosmicsOppHorn
void FillNDRecoPngVars(const ndreco::NDRecoPngObj &pionreco, SRFuzzyKProng &srPng)
Definition: FillReco.cxx:608
double HalfW() const
Definition: CellGeo.cxx:191
Truth info for all neutrinos in the spill.
A collection of associated CellHits.
Definition: Cluster.h:47
Atom< string > CVNTrainingDataLabel
Atom< string > NuMIBeamLabel
static bool sortProngCalE(const SRFuzzyKProng &a, const SRFuzzyKProng &b)
std::pair< int, int > calcFirstLastLivePlane(int plane, std::bitset< 14 > binary, caf::Det_t det=caf::kFARDET)
Atom< InputTag > FuzzyKWeight3DLabel
Atom< string > ReweightLabel
Atom< bool > IsSinglesOverlay
OStream cerr
Definition: OStream.cxx:7
Atom< bool > FillNuTree
Atom< bool > FillPixelMaps
std::map< std::string, std::string > & GetMetadata()
caf::StandardRecord * rec
Definition: tutCAFMacro.C:20
Atom< string > SliceLIDLabel
void abs(TH1 *hist)
Atom< string > CVNLabelLoosePreselPtPOppHorn
int KalBakCellND() const
Definition: CosRejObj.cxx:670
Atom< string > DiscreteTrack2dLabel
Atom< string > DiFShowerCVNLabelNoCosmics
Atom< string > SpillQualLabel
Atom< string > GeneratorLabel
TString ip
Definition: loadincs.C:5
double corr(TGraph *g, double thresh)
Atom< bool > UseEnergy
Atom< string > CosRejLabel
double DistToFront(TVector3 vertex)
const PlaneGeo * Plane(unsigned int i) const
double DistToBack(TVector3 vertex)
Atom< string > CVNLabelLoosePreselPtP
float nutauVal
Cosmic CVN nutau score for each time slice.
bool isRealData() const
Definition: Event.h:83
SRGeant4Weights Geant4Reweights(const g4rwgt::G4WeightTable &g4wgts)
Definition: FillTruth.cxx:880
Atom< string > MRParentSliceLabel
Defines an enumeration for prong classification.
void FillNDRecoTrkVars(const ndreco::NDRecoTrkObj &pionreco, SRKalmanTrack &srTrk)
Functions to fill SRKalmanTrack pion energy reconstruction information.
Definition: FillReco.cxx:576
Atom< string > CVNParticleOppHornLabel
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
int KalFwdCell() const
Definition: CosRejObj.cxx:652
SRFluxWeights FluxReweights(const fxwgt::FluxWeights &flxwgts)
Definition: FillTruth.cxx:870
Atom< bool > ApplyingFilter
std::vector< float > fOutput
Vector of outputs from neural net.
Definition: Result.h:30
Det_t fDet
Detector ID in caf namespace typedef.
void FillProngVars(const rb::Prong &prng, T &srPrng, const std::vector< rb::Cluster > &sliceList, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, int sliceIdx)
Definition: FillReco.cxx:179
void FillSlidVars(const slid::ShowerLID &slid, SRShowerLID &shwlid)
Definition: FillPIDs.cxx:202
Atom< InputTag > FuzzyKWeight2DLabel
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
void FillSliceLID(const SliceLID::Prediction &artSliceLID, caf::SRSliceLID &cafSliceLID)
Definition: FillPIDs.cxx:728
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
std::vector< int > SliceToOrderedNuIds(const std::vector< cheat::NeutrinoWithIndex > &nusWithIdx, const std::vector< std::vector< cheat::NeutrinoEffPur >> &slTruthTable, std::function< double(const cheat::NeutrinoEffPur &)> slMetric, std::function< double(const cheat::NeutrinoEffPur &)> nuMetric) const
Given a vector of indexed neutrino interaction and a vector of slice truth tables, returns a vector of neutrino interaction indices ordered by best match to the corresponding slice. Here, best match is determined according to the given cheat::NeutrinoEffPur functions for the slice and the nu.
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
Atom< string > SPIDLabel
TVector3 MaxXYZ() const
Definition: Cluster.cxx:492
Atom< string > RegCVNLabel
Atom< string > CVNPixelMapLabel
const Var kDistAllEast([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.distallpngeast)) return-1000.0f;return float(sr->sel.nuecosrej.distallpngeast);})
Distance of all showers in slice from the east edge of detector.
Definition: NueVars.h:39
Atom< string > SingleMixerLabel
float nueVal
Cosmic CVN nue score for each time slice.
void CopyRemidVars(const SRKalmanTrack &srTrk, caf::SRRemid &remid)
Definition: FillPIDs.cxx:712
PID
Definition: FillPIDs.h:14
void FillNCPi0BkgRejVars(const ncpi0::NCPi0BkgRej &ncpi0bkgrej, caf::SRNCPi0BkgRej &srncpi0bkgrej)
Definition: FillPIDs.cxx:530
bool isValid() const
Definition: Handle.h:189