BackTracker.h

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file  BackTracker.h
/// \brief back track the reconstruction to the simulation
///
/// \version $Id: BackTracker.h,v 1.37 2012-12-06 15:03:39 lein Exp $
/// \author  brebel@fnal.gov
////////////////////////////////////////////////////////////////////////


#ifndef CHEAT_BACKTRACKER_H
#define CHEAT_BACKTRACKER_H

#include <map>
#include <unordered_map>
#include <vector>

#include "GeometryObjects/OfflineChan.h"
#include "GeometryObjects/PlaneGeo.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/WeightedHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "Simulation/ParticleNavigator.h"
#include "Simulation/Particle.h"
#include "Simulation/TrueEnergy.h"
#include "Simulation/FLSHit.h"
#include "Simulation/PhotonSignal.h"
#include "nusimdata/SimulationBase/MCTruth.h"

#include "art/Framework/Services/Registry/ActivityRegistry.h"
#include "art/Framework/Services/Registry/ServiceMacros.h"
#include "fhiclcpp/types/Atom.h"

///code to link reconstructed objects back to the MC truth information
namespace cheat{

  struct TrackIDE
  {
    int trackID;      ///< Geant4 supplied trackID
    float energyFrac; ///< fraction of hit energy from the particle with this trackID
    float energy;     ///< energy from the particle with this trackID
  };

  struct NeutrinoEffPur
  {
    art::Ptr<simb::MCTruth> neutrinoInt;   ///< Truth information about neutrino interaction
    double efficiency;                     ///< Efficiency (based on FLS energy) of neutrino interaction relative to slice
    double purity;                         ///< Purity (based on FLS energy) of neutrino interaction relative to slice
    unsigned int    truthColIndex;         ///< Position in column vector of MCTruths from the event
    double energySlice;                    ///< Sum of FLS hits from the neutrino contributing to hits included in the slice
    double energyTotal;                    ///< Sum of FLS hits from the neutrino contributing to any hit in the event
    int nSliceHits;                        ///< Number of hits from this neutrino in this slice
    int nTotalHits;                        ///< Total number of hits the neutrino left in the detector

  };

  struct NeutrinoWithIndex
  {
    art::Ptr<simb::MCTruth> neutrinoInt;   ///< Truth information about neutrino interaction
    unsigned int    truthColIndex;         ///< Position in column vector of MCTruths from the event
  };

  struct ParticleEffPur
  {
    int trackID;                           ///< Truth particle track ID
    double efficiency;                     ///< Efficiency of particle relative to track
    double purity;                         ///< Putity of particle relative to track
    int pdg;                               ///< Pdg of particle
  };

  /// Does \a a have less energy than \a b?
  bool CompareByEnergy(const TrackIDE& a, const TrackIDE& b);

  /// Function for NeutrinoEffPur's nu interaction to slice efficiency
  double EffMetric(const cheat::NeutrinoEffPur& ep);

  /// Function for NeutrinoEffPur's nu interaction to slice purity
  double PurMetric(const cheat::NeutrinoEffPur& ep);

  /// Function for NeutrinoEffPur's nu interaction to slice efficiency * purity
  double EffPurMetric(const cheat::NeutrinoEffPur& ep);

  /// Function for NeutrinoEffPur's nu interaction to slice energy
  double EnergyMetric(const cheat::NeutrinoEffPur& ep);

  /// Many functions come in four variants. The core implementation is in terms
  /// of vector<const T*>, which is the common basis all the others can be
  /// converted to. If you have an art::Handle<vector<T>> you can just
  /// dereference it to get this type. If you happen to have a vector<T> or
  /// vector<Ptr<T>> or PtrVector<T>, overloads are provided so that those are
  /// just as easy to work with.
  class BackTracker
  {
  public:
    struct Params
    {
      template<class T> using Atom = fhicl::Atom<T>;
      using Comment = fhicl::Comment;
      using Name = fhicl::Name;

      Atom<std::string> GeantModuleLabel{Name("GeantModuleLabel"), Comment("Module that produced the sim::Particle objects")};
      Atom<std::string> PhotonModuleLabel{Name("PhotonModuleLabel"), Comment("Module that produced the sim::PhotonSignal objects")};
      Atom<bool> DoHitsCheck{Name("DoHitsCheck"), Comment("Do hit by hit check in HitCollectionEfficiency?")};
      Atom<double> MinContribFrac{Name("MinContribFrac"), Comment("Minimum fractional contribution a particle must make to a CellHit's energy to be counted"), 0.1};
      Atom<int> MinPhysicsHits{Name("MinPhysicsHits"), Comment("Minimum number of hits in a cluster that are not noise for the function IsHitsNoise to return false"), 4};
      Atom<double> MinPhysicsFrac{Name("MinPhysicsFrac"), Comment("If cluster fails to have MinPhysicsHits, can still have function IsHitsNoise return false if fraction of physics hits is at least MinPhysicsFrac"), 0.5};
      Atom<bool> MRMode{Name("MRMode"), Comment("Are we running over a muon removed file?")};
      Atom<bool> MCOverlayMode{Name("MCOverlayMode"), Comment("If MC is overlaid on real data, still backtrack")};
      Atom<std::string> TrueEnergyModuleLabel{Name("TrueEnergyModuleLabel"), Comment("Module that produced the sim::TrueEnergy objects")};
    };

    // Allows 'nova --print-description' to work
    using Parameters = art::ServiceTable<Params>;

    BackTracker(const Parameters& params, art::ActivityRegistry& reg);
    ~BackTracker();

    /// This function rebuilds the various maps we need to answer backtracking queries.
    /// It is called automatically before each event is processed. For jobs involving
    /// Monte Carlo generation, this is too soon. So, ReadoutSim also calls this function
    /// when it knows all the relevant products have been produced. Requiring this kind
    /// of cooperation from other modules for a service to function properly is not
    /// recommended in general. In this case it seems like the least-bad option, but
    /// think very carefully before adding any other examples. "Do as I say, not as I do".
    void Rebuild(const art::Event& evt);

    /// \brief Is this a file with truth info in? (Is BackTracker going to be
    /// any use to you?)
    ///
    /// True if we succeeded in retrieving FLSHits, ie in general for MC
    /// files. Expected false in data, and mock data etc.
    bool HaveTruthInfo() const {return fHaveTruthInfo;}

    /// \brief Returns the PhotonSignals contributing the signal in the
    /// specified cell. WARNING: Use with extreme caution! Most users should use HitToPhotonSignal
    /// instead. Only use the cell version if you really want to know ALL the photon signals
    /// in the cell for the entire time window, NOT time-associated with the hit.
    ///
    /// Sorted from highest signal to lowest (which is probably not much use,
    /// since energy from one particle can be split over several signals).
    const std::vector<sim::PhotonSignal> CellToPhotonSignal(unsigned int const& plane,
                                                            unsigned int const& cell) const;

    /// \brief Returns the PhotonSignals contributing the signal in the
    /// specified hit.
    ///
    /// Sorted from highest signal to lowest (which is probably not much use,
    /// since energy from one particle can be split over several signals).
    const std::vector<sim::PhotonSignal> HitToPhotonSignal(const art::Ptr<rb::CellHit>& hit) const;

    /// \brief Returns the FLSHits contributing the signal in the specified
    /// cell. WARNING: Use with extreme caution! Most users should use HitToFLSHit
    /// instead. Only use the cell version if you really want to know ALL the FLSHits
    /// in the cell for the entire time window, NOT time-associated with the hit.
    ///
    /// The vector is sorted to return in order of highest contributing Geant4
    /// track to lowest.
    const std::vector<sim::FLSHit>    CellToFLSHit(unsigned int const& plane,
                                                   unsigned int const& cell) const;

    /// \brief All the FLSHits that contributed to this hit, sorted from most
    /// to least light
    std::vector<sim::FLSHit> HitToFLSHit(const rb::CellHit& hit) const;

    std::vector<sim::FLSHit> HitToFLSHit(const art::Ptr<rb::CellHit>& hit) const
    {
      return HitToFLSHit(*hit);
    }

    /// \brief All the FLSHits that were created by the track id trackID,
    /// sorted from most to least light
    std::vector<sim::FLSHit> ParticleToFLSHit(const int& trackID) const;
 
    /// \brief All the FLSHits that were created by the track id trackID,
    /// which also have a particular pdg. Most users should use the other
    /// version of this function. However, sometimes when our truth is 
    /// condensed, the same mother track ID is assigned to FLS hits made 
    /// by non-important daughters, who can have different PDG values. An
    /// example of this is a muon track, whose delta rays often have the track
    /// id of the mother muon. By using this version of the function, you can 
    /// ensure you are looking only at the FLS hits deposited by the muon.
    std::vector<sim::FLSHit> ParticleToFLSHit(const int trackID, const int pdg) const;

    /// \brief The total energy deposited in this cell by all FLSHits with this
    /// trackId
    const double EnergyFromTrackId(const art::Ptr<rb::CellHit>& hit, int trackId, bool useBirksE = false) const;

    /// \brief Returns the sim::Particle object contributing the most light to
    /// the specified rb::CellHit
    const sim::Particle*              HitToParticle(art::Ptr<rb::CellHit> const& hit,
                                                    bool quiet = false) const;

    /// \brief Returns the sim::Particle object contributing the most light to
    /// the specified rb::CellHit
    const sim::Particle*              HitToParticle(const rb::CellHit& hit,
                                                    bool quiet = false) const;

    /// Is this hit not associated with any particles?
    bool                              IsNoise(const art::Ptr<rb::CellHit>& hit) const;
    /// Is this hit not associated with any particles?
    bool                              IsNoise(const rb::CellHit& hit) const;

    /// Of this collection of hits, is it mostly noise or physics hits? True indicates noise.
    bool                              IsHitsNoise(std::vector< art::Ptr<rb::CellHit> > const& hits) const;

    /// \brief Returns a pointer to the sim::Particle object corresponding to
    /// the given TrackID
    const sim::Particle*              TrackIDToParticle(int const& id) const;
    const sim::Particle*              TrackIDToMotherParticle(int const& id) const;

    const art::Ptr<simb::MCTruth>&    TrackIDToMCTruth(int const& id) const;
    const art::Ptr<simb::MCTruth>&    ParticleToMCTruth(const sim::Particle* p) const;
    std::vector<const sim::Particle*> MCTruthToParticles(art::Ptr<simb::MCTruth> const& mct) const;

    /// \brief Returns the XYZ position of the energy deposition for a given
    /// hit.
    ///
    /// Note the return value is in cell coordinates, not detector
    TVector3 HitToXYZ(art::Ptr<rb::CellHit> const& hit, bool useBirksE = false) const;

    /// \brief Returns vector of sim::Particle objects contributing to the
    /// given collection of hits.
    ///
    /// Orders the returned vector with the particle contributing the most
    /// light in the collection first.
    std::vector<const sim::Particle*> HitsToParticle(const std::vector<const rb::CellHit*>& hits) const;
    std::vector<const sim::Particle*> HitsToParticle(const std::vector<rb::CellHit>& hits) const
    {
      return HitsToParticle(VecToVecRawPtr(hits));
    }
    std::vector<const sim::Particle*> HitsToParticle(const std::vector<art::Ptr<rb::CellHit>>& hits) const
    {
      return HitsToParticle(VecPtrToVecRawPtr(hits));
    }
    std::vector<const sim::Particle*> HitsToParticle(const art::PtrVector<rb::CellHit>& hits) const
    {
      return HitsToParticle(PtrVecToVecRawPtr(hits));
    }

    /// \brief Returns vector of TrackIDE structs contributing to the given
    /// collection of hits.
    ///
    /// Orders the returned vector with the trackID contributing the most
    /// energy in the collection first.
    std::vector<TrackIDE> HitsToTrackIDE(const std::vector<const rb::CellHit*>& hits, bool useBirksE = false) const;

    std::vector<TrackIDE> HitsToTrackIDE(const std::vector<rb::CellHit>& hits, bool useBirksE = false) const
    {
      return HitsToTrackIDE(VecToVecRawPtr(hits), useBirksE);
    }

    std::vector<TrackIDE> HitsToTrackIDE(const std::vector<art::Ptr<rb::CellHit>>& hits, bool useBirksE = false) const
    {
      return HitsToTrackIDE(VecPtrToVecRawPtr(hits), useBirksE);
    }

    std::vector<TrackIDE> HitsToTrackIDE(const art::PtrVector<rb::CellHit>& hits, bool useBirksE = false) const
    {
      return HitsToTrackIDE(PtrVecToVecRawPtr(hits), useBirksE);
    }

    /// Convenience function. \ref HitsToTrackIDE but for a single hit.
    std::vector<TrackIDE> HitToTrackIDE(const rb::CellHit& hit, bool useBirksE = false) const;

    std::vector<TrackIDE> HitToTrackIDE(const art::Ptr<rb::CellHit>& hit, bool useBirksE = false) const;


    /// \brief Returns the fraction of hits in a collection that come from the
    /// specified Geant4 track ids.
    ///
    double HitCollectionPurity(const std::set<int>&                trackIDs,
                               const std::vector<rb::WeightedHit>& whits,
                               std::map<int, double>*              purMap = 0,
                               std::map<int, int>*                 parents = 0,
                               bool                                energyPur = false) const;

    /// \brief Returns the fraction of hits in a collection that come from the
    /// specified Geant4 track ids.
    ///
    /// The optional \a purMap parameter returns the purities for each of the
    /// track IDs individually.
    ///
    /// The optional \a parents parameter maps particle IDs to the IDs they
    ///
    /// The optional \a energyPur parameter allows the purity to be calculated based on FLS energy instead of number of hits
    double HitCollectionPurity(const std::set<int>&               trackIDs,
                               const std::vector<const rb::CellHit*>& hits,
                               std::map<int, double>*             purMap = 0,
                               std::map<int, int>*                parents = 0,
                               bool                               energyPur = false) const;

    double HitCollectionPurity(const std::set<int>&            trackIDs,
                               const std::vector<rb::CellHit>& hits,
                               std::map<int, double>*          purMap = 0,
                               std::map<int, int>*             parents = 0,
                               bool                            energyPur = false) const
    {
      return HitCollectionPurity(trackIDs, VecToVecRawPtr(hits), purMap, parents, energyPur);
    }

    double HitCollectionPurity(const std::set<int>&               trackIDs,
                               const std::vector<art::Ptr<rb::CellHit>>& hits,
                               std::map<int, double>*             purMap = 0,
                               std::map<int, int>*                parents = 0,
                               bool                               energyPur = false) const
    {
      return HitCollectionPurity(trackIDs, VecPtrToVecRawPtr(hits), purMap, parents, energyPur);
    }

    double HitCollectionPurity(const std::set<int>&               trackIDs,
                               const art::PtrVector<rb::CellHit>& hits,
                               std::map<int, double>*             purMap = 0,
                               std::map<int, int>*                parents = 0,
                               bool                               energyPur = false) const
    {
      return HitCollectionPurity(trackIDs, PtrVecToVecRawPtr(hits), purMap, parents, energyPur);
    }

    /// \brief Returns the fraction of all energy in an event from a specific set
    /// of Geant4 track IDs that are represented in a collection of weighted hits
    ///
    /// \a desiredHits and 
    /// \a totalHits are also optional parameters that are the number of hits from the given set of 
    /// trackIds in the given slice and the total number of hits from these trackId's in the detector respectively
    double HitCollectionEfficiency(const std::set<int>&   trackIDs,
                                   const std::vector<rb::WeightedHit>& whits,
                                   const std::vector<const rb::CellHit*>& allhits,
                                   const geo::View_t&     view,
                                   std::map<int, double>* effMap = 0,
                                   bool                   energyEff = false,
                                   double*                desiredEnergy = 0,
                                   double*                totalEnergy = 0,
                                   int*                   desiredHits = 0,
                                   int*                   totalHits = 0) const;

    double HitCollectionEfficiency(const std::set<int>&   trackIDs,
                                   const std::vector<rb::WeightedHit>& whits,
                                   const std::vector<rb::CellHit>& allhits,
                                   const geo::View_t&     view,
                                   std::map<int, double>* effMap = 0,
                                   bool                   energyEff = false,
                                   double*                desiredEnergy = 0,
                                   double*                totalEnergy = 0,
                                   int*                   desiredHits = 0,
                                   int*                   totalHits = 0) const
    {
      return HitCollectionEfficiency(trackIDs, whits, VecToVecRawPtr(allhits), view, effMap, energyEff, desiredEnergy, totalEnergy, desiredHits, totalHits);
    }

    double HitCollectionEfficiency(const std::set<int>&   trackIDs,
                                   const std::vector<rb::WeightedHit>& whits,
                                   const std::vector<art::Ptr<rb::CellHit>>& allhits,
                                   const geo::View_t&     view,
                                   std::map<int, double>* effMap = 0,
                                   bool                   energyEff = false,
                                   double*                desiredEnergy = 0,
                                   double*                totalEnergy = 0,
                                   int*                   desiredHits = 0,
                                   int*                   totalHits = 0) const
    {
      return HitCollectionEfficiency(trackIDs, whits, VecPtrToVecRawPtr(allhits), view, effMap, energyEff, desiredEnergy, totalEnergy, desiredHits, totalHits);
    }

    double HitCollectionEfficiency(const std::set<int>&   trackIDs,
                                   const std::vector<rb::WeightedHit>& whits,
                                   const art::PtrVector<rb::CellHit>&  allhits,
                                   const geo::View_t&     view,
                                   std::map<int, double>* effMap = 0,
                                   bool                   energyEff = false,
                                   double*                desiredEnergy = 0,
                                   double*                totalEnergy = 0,
                                   int*                   desiredHits = 0,
                                   int*                   totalHits = 0) const
    {
      return HitCollectionEfficiency(trackIDs, whits, PtrVecToVecRawPtr(allhits), view, effMap, energyEff, desiredEnergy, totalEnergy, desiredHits, totalHits);
    }



    /// \brief Returns the fraction of all hits in an event from a specific set
    /// of Geant4 track IDs that are represented in a collection of hits
    ///
    /// The optional \a effMap parameter returns the efficiencies for each of
    /// the track IDs individually.
    ///
    /// The optional \a energyPur parameter allows the efficiency to be calculated based on FLS energy instead of number of hits
    /// \a desiredHits and \a totalHits are also optional parameters that are the number of hits from the given set of trackIds 
    /// in the given slice and the total number of hits from these trackId's in the detector respectively
    double HitCollectionEfficiency(const std::set<int>&   trackIDs,
                                   const std::vector<const rb::CellHit*>& hits,
                                   const std::vector<const rb::CellHit*>& allhits,
                                   const geo::View_t&     view,
                                   std::map<int, double>* effMap = 0,
                                   bool                   energyEff = false,
                                   double*                desiredEnergy = 0,
                                   double*                totalEnergy = 0,
                                   int*                   desiredHits = 0,
                                   int*                   totalHits = 0) const;

    double HitCollectionEfficiency(const std::set<int>&   trackIDs,
                                   const std::vector<rb::CellHit>& hits,
                                   const std::vector<rb::CellHit>& allhits,
                                   const geo::View_t&     view,
                                   std::map<int, double>* effMap = 0,
                                   bool                   energyEff = false,
                                   double*                desiredEnergy = 0,
                                   double*                totalEnergy = 0,
                                   int*                   desiredHits = 0,
                                   int*                   totalHits = 0) const
    {
      return HitCollectionEfficiency(trackIDs,
                                     VecToVecRawPtr(hits),
                                     VecToVecRawPtr(allhits),
                                     view,
                                     effMap, energyEff, desiredEnergy,
                                     totalEnergy, desiredHits, totalHits);
    }

    double HitCollectionEfficiency(const std::set<int>&   trackIDs,
                                   const std::vector<art::Ptr<rb::CellHit>>& hits,
                                   const std::vector<art::Ptr<rb::CellHit>>& allhits,
                                   const geo::View_t&     view,
                                   std::map<int, double>* effMap = 0,
                                   bool                   energyEff = false,
                                   double*                desiredEnergy = 0,
                                   double*                totalEnergy = 0,
                                   int*                   desiredHits = 0,
                                   int*                   totalHits = 0) const
    {
      return HitCollectionEfficiency(trackIDs,
                                     VecPtrToVecRawPtr(hits),
                                     VecPtrToVecRawPtr(allhits),
                                     view,
                                     effMap, energyEff, desiredEnergy,
                                     totalEnergy, desiredHits, totalHits);
    }

    double HitCollectionEfficiency(std::set<int>                      trackIDs,
                                   const art::PtrVector<rb::CellHit>& hits,
                                   const art::PtrVector<rb::CellHit>& allhits,
                                   const geo::View_t&                 view,
                                   std::map<int, double>*             effMap = 0,
                                   bool                               energyEff     = false,
                                   double*                            desiredEnergy = 0,
                                   double*                            totalEnergy   = 0,
                                   int*                               desiredHits   = 0,
                                   int*                               totalHits     = 0) const
    {
      return HitCollectionEfficiency(trackIDs,
                                     PtrVecToVecRawPtr(hits),
                                     PtrVecToVecRawPtr(allhits),
                                     view,
                                     effMap, energyEff, desiredEnergy,
                                     totalEnergy, desiredHits, totalHits);
    }

    /// \brief Tool for NumuEAna that allows one to see if primary muon (or any track ID you feed the function)
    /// contributes to enough hits in the slice to pass the hit cut. It is currently set at 4 hits in each view.
    bool PassMuonCuts(int trackID, art::PtrVector<rb::CellHit> const& hits) const;

    /// \brief Tool for function SliceToNeutrinoInteractions, tells how to sort by putting structures with highest efficiency first.
    static bool sort_eff(const NeutrinoEffPur& a, const NeutrinoEffPur& b);

    /// \brief Tool for function SliceToNeutrinoInteractions, tells how to sort by putting structures with highest purity first.
    static bool sort_pur(const NeutrinoEffPur& a, const NeutrinoEffPur& b);

    /// \brief Tool for function SliceToNeutrinoIndex, tells how to sort by putting structures with highest efficiency first.
    static bool sort_idxeff(const std::pair<int,NeutrinoEffPur>& a, const std::pair<int,NeutrinoEffPur>& b);

    /// \brief Tool for function SliceToNeutrinoIndex, tells how to sort by putting structures with highest purity first.
    static bool sort_idxpur(const std::pair<int,NeutrinoEffPur>& a, const std::pair<int,NeutrinoEffPur>& b);


    /// \brief Given a collection of hits (often a slice), returns vector of neutrino indices 
    /// corresponding to the vector from SliceToNeutrinoInteractions.
    std::vector<int>
    SliceToNeutrinoIndex(const std::vector<const rb::CellHit*>& sliceHits,
                         const std::vector<const rb::CellHit*>& allHits,
                         bool sortPur = false) const;

    /// \brief Given a collection of hits (often a slice), returns vector of structures of neutrino interactions,
    /// efficiency, and purity of that neutrino interaction relative to the slice. Efficiency is defined as FLS energy
    /// from neutrino interaction in slice / total FLS energy from neutrino interaction in event.
    /// This vector is sorted from the highest efficiency interaction to lowest.
    std::vector<NeutrinoEffPur>
    SliceToNeutrinoInteractions(const std::vector<const rb::CellHit*>& sliceHits,
                                const std::vector<const rb::CellHit*>& allHits,
                                bool sortPur = false) const;

    std::vector<NeutrinoEffPur>
    SliceToNeutrinoInteractions(const std::vector<rb::CellHit>& sliceHits,
                                const std::vector<rb::CellHit>& allHits,
                                bool sortPur = false) const
    {
      return SliceToNeutrinoInteractions(VecToVecRawPtr(sliceHits),
                                         VecToVecRawPtr(allHits),
                                         sortPur);
    }

    std::vector<NeutrinoEffPur>
    SliceToNeutrinoInteractions(const std::vector<art::Ptr<rb::CellHit>>& sliceHits,
                                const std::vector<art::Ptr<rb::CellHit>>& allHits,
                                bool sortPur = false) const
    {
      return SliceToNeutrinoInteractions(VecPtrToVecRawPtr(sliceHits),
                                         VecPtrToVecRawPtr(allHits),
                                         sortPur);
    }

    std::vector<NeutrinoEffPur>
    SliceToNeutrinoInteractions(const art::PtrVector<rb::CellHit>& sliceHits,
                                const art::PtrVector<rb::CellHit>& allHits,
                                bool sortPur = false) const
    {
      return SliceToNeutrinoInteractions(PtrVecToVecRawPtr(sliceHits),
                                         PtrVecToVecRawPtr(allHits),
                                         sortPur);
    }

    std::vector<NeutrinoEffPur>
    SliceToNeutrinoInteractions(const art::Ptr<rb::Cluster>& sliceCluster,
                                const std::vector<const rb::CellHit*>& allHits,
                                bool sortPur = false) const;

    std::vector<NeutrinoEffPur>
    SliceToNeutrinoInteractions(const art::Ptr<rb::Cluster>& sliceCluster,
                                const std::vector<rb::CellHit>& allHits,
                                bool sortPur = false) const
    {
      return SliceToNeutrinoInteractions(sliceCluster,
                                         VecToVecRawPtr(allHits),
                                         sortPur);
    }

    std::vector<NeutrinoEffPur>
    SliceToNeutrinoInteractions(const art::Ptr<rb::Cluster>& sliceCluster,
                                const std::vector<art::Ptr<rb::CellHit>>& allHits,
                                bool sortPur = false) const
    {
      return SliceToNeutrinoInteractions(sliceCluster,
                                         VecPtrToVecRawPtr(allHits),
                                         sortPur);
    }

    std::vector<NeutrinoEffPur>
    SliceToNeutrinoInteractions(const art::Ptr<rb::Cluster>& sliceCluster,
                                const art::PtrVector<rb::CellHit>& allHits,
                                bool sortPur = false) const
    {
      return SliceToNeutrinoInteractions(sliceCluster,
                                         PtrVecToVecRawPtr(allHits),
                                         sortPur);
    }


    ///\brief 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. 
    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;

    std::vector <int> SliceToOrderedNuIdsByEff(const std::vector<cheat::NeutrinoWithIndex>&           nusWithIdx,
                                               const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable) const;

    std::vector <int> SliceToOrderedNuIdsByPur(const std::vector<cheat::NeutrinoWithIndex>&           nusWithIdx,
                                               const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable) const;

    std::vector <int> SliceToOrderedNuIdsByEffPur(const std::vector<cheat::NeutrinoWithIndex>&           nusWithIdx,
                                                  const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable) const;

    std::vector <int> SliceToOrderedNuIdsByEnergy(const std::vector<cheat::NeutrinoWithIndex>&           nusWithIdx,
                                                  const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable) const;

    /// \brief Given a collection of hits (often a slice), returns vector of structures of MCTruth,
    /// efficiency, and purity of that neutrino interaction relative to the slice. Efficiency is defined as FLS energy
    /// from neutrino interaction in slice / total FLS energy from neutrino interaction in event.
    /// This vector is sorted from the highest efficiency interaction to lowest.  This function returns all MCTruth,
    /// including those without neutrino interactions, i.e. cosmics.
    std::vector<NeutrinoEffPur>
    SliceToMCTruth(const std::vector<const rb::CellHit*>& sliceHits,
                   const std::vector<const rb::CellHit*>& allHits,
                   bool sortPur = false) const;

    std::vector<NeutrinoEffPur>
    SliceToMCTruth(const std::vector<rb::CellHit>& sliceHits,
                   const std::vector<rb::CellHit>& allHits,
                   bool sortPur = false) const
    {
      return SliceToMCTruth(VecToVecRawPtr(sliceHits),
                            VecToVecRawPtr(allHits),
                            sortPur);
    }

    std::vector<NeutrinoEffPur>
    SliceToMCTruth(const std::vector<art::Ptr<rb::CellHit>>& sliceHits,
                   const std::vector<art::Ptr<rb::CellHit>>& allHits,
                   bool sortPur = false) const
    {
      return SliceToMCTruth(VecPtrToVecRawPtr(sliceHits),
                            VecPtrToVecRawPtr(allHits),
                            sortPur);
    }

    std::vector<NeutrinoEffPur>
    SliceToMCTruth(const art::PtrVector<rb::CellHit>& sliceHits,
                   const art::PtrVector<rb::CellHit>& allHits,
                   bool sortPur = false) const
    {
      return SliceToMCTruth(PtrVecToVecRawPtr(sliceHits),
                            PtrVecToVecRawPtr(allHits),
                            sortPur);
    }

    std::vector<NeutrinoEffPur>
    SliceToMCTruth(const art::Ptr<rb::Cluster>& sliceCluster,
                   const std::vector<const rb::CellHit*>& allHits,
                   bool sortPur = false) const;

    std::vector<NeutrinoEffPur>
    SliceToMCTruth(const art::Ptr<rb::Cluster>& sliceCluster,
                   const std::vector<rb::CellHit>& allHits,
                   bool sortPur = false) const
    {
      return SliceToMCTruth(sliceCluster, VecToVecRawPtr(allHits), sortPur);
    }

    std::vector<NeutrinoEffPur>
    SliceToMCTruth(const art::Ptr<rb::Cluster>& sliceCluster,
                   const std::vector<art::Ptr<rb::CellHit>>& allHits,
                   bool sortPur = false) const
    {
      return SliceToMCTruth(sliceCluster, VecPtrToVecRawPtr(allHits), sortPur);
    }

    std::vector<NeutrinoEffPur>
    SliceToMCTruth(const art::Ptr<rb::Cluster>& sliceCluster,
                   const art::PtrVector<rb::CellHit>& allHits,
                   bool sortPur = false) const
    {
      return SliceToMCTruth(sliceCluster, PtrVecToVecRawPtr(allHits), sortPur);
    }

    /// \brief Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures of MCTruth,
    /// efficiency, and purity of that MCTruth interaction relative to the slice. 
    /// The returned vector is the same size as the total number of slices, including the noise slice, which is treated the same
    /// as any other slice. It is sorted in the same order as the slice list.
    /// For each entry, there is a vector of NeutrinoEffPurs. This vector is the same size as the MCTruth list. It is sorted in
    /// the same order as the MCTruth list.
    /// Every slice/MCTruth pair shows up, even if the contribution is zero.
    /// Efficiency is defined as FLS energy from a MCTruth in a slice / total FLS energy from a MCTruth in event.
    /// This function is agnostic to MCTruth and treates them all the same, including those without neutrino interactions, i.e. cosmics.
    std::vector<std::vector<cheat::NeutrinoEffPur>>
    SlicesToMCTruthsTable(const std::vector<const rb::Cluster*>& sliceList) const;

    std::vector<std::vector<cheat::NeutrinoEffPur>>
    SlicesToMCTruthsTable(const std::vector<rb::Cluster>& sliceList) const
    {
      return SlicesToMCTruthsTable(VecToVecRawPtr(sliceList));
    }

    std::vector<std::vector<cheat::NeutrinoEffPur>>
    SlicesToMCTruthsTable(const std::vector<art::Ptr<rb::Cluster>>& sliceList) const
    {
      return SlicesToMCTruthsTable(VecPtrToVecRawPtr(sliceList));
    }

    std::vector<std::vector<cheat::NeutrinoEffPur>>
    SlicesToMCTruthsTable(const art::PtrVector<rb::Cluster>& sliceList) const
    {
      return SlicesToMCTruthsTable(PtrVecToVecRawPtr(sliceList));
    }

    /// \brief Given a collection of reconstructed tracks, this method finds the best match for each reco track to the
    /// input track id particles. 
    /// 
    /// For each reconstructed track, the true particle that best matches the track by
    /// purity from the input trkIDs is determined. Each input track ID is only allowed to match at most one input
    /// track and any track IDs specified in the input allowedDaughters are not counted against a tracks purity.
    /// The tracks are matched to particles in the order that they appear in the input vector. Each entry in the
    /// returned vector of ParticleEffPur corresponds to the track in the same position as the input track vector and
    /// holds the matched track id, efficiency, purity, and pdg of the best matched particle for the given track. If
    /// no track id is found to match the track, then the purity, efficiency, and pdg will be zero. You need to check
    /// that a match exists before trying to use the matched track id in other functions
    std::vector<ParticleEffPur>
    TracksToParticles(const std::vector<const rb::Track*>& tracks,
                      const std::set<int>& trkIDs,
                      const std::set<int>& allowedDaughters,
                      const std::vector<const rb::CellHit*>& allHits) const;

    std::vector<ParticleEffPur>
    TracksToParticles(const std::vector<rb::Track>& tracks,
                      const std::set<int>& trkIDs,
                      const std::set<int>& allowedDaughters,
                      const std::vector<rb::CellHit>& allHits) const
    {
      return TracksToParticles(VecToVecRawPtr(tracks),
                               trkIDs, allowedDaughters,
                               VecToVecRawPtr(allHits));
    }

    std::vector<ParticleEffPur>
    TracksToParticles(const std::vector<art::Ptr<rb::Track>>& tracks,
                      const std::set<int>& trkIDs,
                      const std::set<int>& allowedDaughters,
                      const std::vector<art::Ptr<rb::CellHit>>& allHits) const
    {
      return TracksToParticles(VecPtrToVecRawPtr(tracks),
                               trkIDs, allowedDaughters,
                               VecPtrToVecRawPtr(allHits));
    }

    std::vector<ParticleEffPur>
    TracksToParticles(const art::PtrVector<rb::Track>& tracks,
                      const std::set<int>& trkIDs,
                      const std::set<int>& allowedDaughters,
                      const art::PtrVector<rb::CellHit>& allHits) const
    {
      return TracksToParticles(PtrVecToVecRawPtr(tracks),
                               trkIDs, allowedDaughters,
                               PtrVecToVecRawPtr(allHits));
    }

    std::vector<ParticleEffPur>
    TracksToParticles(const std::vector<art::Ptr<rb::Track>>& tracks,
                      const std::set<int>& trkIDs,
                      const std::set<int>& allowedDaughters,
                      const art::PtrVector<rb::CellHit>& allHits) const
    {
      return TracksToParticles(VecPtrToVecRawPtr(tracks),
                               trkIDs, allowedDaughters,
                               PtrVecToVecRawPtr(allHits));
    }


    /// \brief Function primarily for CAFMaker - returns a vector of all MCTruth interactions along with their truth vector position
    std::vector <NeutrinoWithIndex> allNeutrinoInteractions() const;

    std::vector <NeutrinoWithIndex> allMCTruth() const;

    /// \brief Get a reference to the ParticleNavigator
    const sim::ParticleNavigator&    ParticleNavigator() const { return fParticleNav; }

    /// \brief Get all G4 track ids present in the event.
    ///
    /// These do not include ids from EM shower secondaries, tertiaries, etc as
    /// they are not saved as particles.
    const std::set<int>               GetTrackIDList() const { return fTrackIDs; }

    /// \brief:
    //  Return a vector of clusters where the hits are separated by MCTruth.  A hit that
    //  had contributions from more than one MCTruth is clustered with the MCTruth that
    //  contributed the most energy.  Noise hits are put into the last cluster.
    //
    //  NOTE: It is not always the case that all MCTruths contribute hits to the event record,
    //  so some of the clusters returned might be empty!
    std::vector<rb::Cluster> ClusterByTruth(std::vector< art::Ptr< rb::CellHit > > const &hits);

    /// \brief:
    //  Given a list of hits (usually all hits) this returns a cluster of only the PRIMARY hits
    //  associated with a specific MCTruth.
    //
    // NOTE: A hit will only be added to this cluster if its primary energy contribution came from this MCTruth.
    rb::Cluster MCTruthToCluster(art::Ptr<simb::MCTruth>                const& truth,
                                 std::vector< art::Ptr< rb::CellHit > > const& hits);

    void HitToParticlesMap(const std::vector<const rb::Cluster*>& sliceList,
                           const int& sliceIdx);

    void HitToParticlesMap(const std::vector<rb::Cluster>& sliceList,
                           const int& sliceIdx)
    {
      return HitToParticlesMap(VecToVecRawPtr(sliceList), sliceIdx);
    }

    void HitToParticlesMap(const std::vector<art::Ptr<rb::Cluster>>& sliceList,
                           const int& sliceIdx)
    {
      return HitToParticlesMap(VecPtrToVecRawPtr(sliceList), sliceIdx);
    }

    void HitToParticlesMap(const art::PtrVector<rb::Cluster>& sliceList,
                           const int& sliceIdx)
    {
      return HitToParticlesMap(PtrVecToVecRawPtr(sliceList), sliceIdx);
    }

    cheat::ParticleEffPur ClusterToParticle(const std::vector<const rb::Cluster*>& sliceList,
                                            const std::vector<const rb::CellHit*>& hits,
                                            const int& sliceIdx);

    cheat::ParticleEffPur ClusterToParticle(const std::vector<rb::Cluster>& sliceList,
                                            const std::vector<rb::CellHit>& hits,
                                            const int& sliceIdx)
    {
      return ClusterToParticle(VecToVecRawPtr(sliceList),
                               VecToVecRawPtr(hits),
                               sliceIdx);
    }

    cheat::ParticleEffPur ClusterToParticle(const std::vector<art::Ptr<rb::Cluster>>& sliceList,
                                            const std::vector<art::Ptr<rb::CellHit>>& hits,
                                            const int& sliceIdx)
    {
      return ClusterToParticle(VecPtrToVecRawPtr(sliceList),
                               VecPtrToVecRawPtr(hits),
                               sliceIdx);
    }

    cheat::ParticleEffPur ClusterToParticle(const art::PtrVector<rb::Cluster>& sliceList,
                                            const art::PtrVector< rb::CellHit >& hits,
                                            const int& sliceIdx)
    {
      return ClusterToParticle(PtrVecToVecRawPtr(sliceList),
                               PtrVecToVecRawPtr(hits),
                               sliceIdx);
    }

    /// \brief
    /// A function to calculate the amount of energy from a given MCParticle which is not contained in the detector.
    /// Will return the energy the particle had when it exited the detector. 
    /// Should a particle stop in the detector then the "Escaping energy" is 0.
    /// It has the option of returning the energy that a particle deposited energy in the detector.
    /// This is calclated from EnEnt - EnEsc, this has the deficiency that depositions from daughters may be doubly counted,
    ///   however this is what the function CalcDaughterEscapingEnergy is for.
    /// There is a flag "UseMCPart", this controls whether to use the saved MCParticles, or to use the more accurate GEANT4
    ///   information, which takes into account all daughter particles.
    ///   Currently it is defaulted to true, ie use the MCParticles, as the GEANT4 info is not usually in the event record yet.
    float CalcEscapingEnergy( const sim::Particle &Par, bool UseMCPart = true, std::string G4Label = "trueens" );
    float CalcEscapingEnergy( const sim::Particle &Par, float& EnDep, bool UseMCPart = true, std::string G4Label = "trueens" );

    /// \brief
    /// A function to be used in conjunction with CalcTotalEscapingEnergy to calculate the energy which the daughters of a 
    ///   a given sim::Particle deposited in the detector.
    float CalcDaughterEscapingEnergy( const sim::Particle &Par );

    /// \brief
    /// An extension of CalcEscapingEnergy which also takes into account the energy which the daughters of given sim::Particle's
    ///   deposit in the detector.
    /// To do this it calls CalcDaughterEscapingEnergy, which in turn calls CalcEscapingEnergy for all of the daughters.
    float CalcTotalEscapingEnergy( const sim::Particle &Par, bool UseMCPart = true, std::string G4Label = "trueens" );
    float CalcTotalEscapingEnergy( const sim::Particle &Par, float& ParEnDep, bool UseMCPart = true, std::string G4Label = "trueens" );

    /// \brief
    /// A function which returns the trajectory point at which the particle first enters and leave the detector volume.
    /// The trajectory point is returned so that the user can determine the energies etc which the particle had at this
    ///  point.
    ///  Should the particle have been generated in the detector then TrajEnt == 0 will be returned.
    ///  Should the particle stop in the detector then TrajEx == NumberTrajectoryPoints()-1 will be returned.
    ///  Should the particle never enter the detector then holder values TrajEnter == TrajEx == -5 will be returned.
    bool   InterceptsDetector( const sim::Particle &Par, int &TrajEnt, int &TrajEx );

  private:
    /// \brief Helper function for implementing overloads
    ///
    /// NB the pointers are only valid so long as the parent \a xs is still
    /// around
    template<class T> std::vector<const T*>
    PtrVecToVecRawPtr(const art::PtrVector<T>& xs) const
    {
      std::vector<const T*> ret;
      ret.reserve(xs.size());
      for(const art::Ptr<T>& x: xs) ret.push_back(x.get());
      return ret;
    }

    /// Helper function for implementing overloads
    ///
    /// NB the pointers are only valid so long as the parent \a xs is still
    /// around
    template<class T> std::vector<const T*>
    VecPtrToVecRawPtr(const std::vector<art::Ptr<T>>& xs) const
    {
      std::vector<const T*> ret;
      ret.reserve(xs.size());
      for(const art::Ptr<T>& x: xs) ret.push_back(x.get());
      return ret;
    }

    /// Helper function for implementing overloads
    ///
    /// NB the pointers are only valid so long as the parent \a xs is still
    /// around
    template<class T> std::vector<const T*>
    VecToVecRawPtr(const std::vector<T>& xs) const
    {
      std::vector<const T*> ret;
      ret.reserve(xs.size());
      for(const T& x: xs) ret.push_back(&x);
      return ret;
    }

    /// Helper for HitCollectionEfficiency and HitCollectionPurity
    class RawWeightedHit
    {
    public:
      RawWeightedHit(const rb::WeightedHit& h) : hit(h.hit.get()), weight(h.weight) {}
      RawWeightedHit(const rb::CellHit* h, double w) : hit(h), weight(w) {}
      const rb::CellHit* hit;
      double weight;
    };

    /// Internal implementation of public \ref HitCollectionPurity
    double HitCollectionPurity(const std::set<int>&               trackIDs,
                               const std::vector<RawWeightedHit>& whits,
                               std::map<int, double>*             purMap = 0,
                               std::map<int, int>*                parents = 0,
                               bool                               energyPur = false) const;

    /// Internal implementation of public \ref HitCollectionEfficiency
    double HitCollectionEfficiency(const std::set<int>&   trackIDs,
                                   const std::vector<RawWeightedHit>& whits,
                                   const std::vector<const rb::CellHit*>& allhits,
                                   const geo::View_t&     view,
                                   std::map<int, double>* effMap = 0,
                                   bool                   energyEff = false,
                                   double*                desiredEnergy = 0,
                                   double*                totalEnergy = 0,
                                   int*                   desiredHits = 0,
                                   int*                   totalHits = 0) const;

    /// \brief Internal helper, decides if a PhotonSignal and CellHit are close
    /// enough to match.
    bool CloseInTime(const rb::CellHit& chit, const sim::PhotonSignal& phot) const;

    /// \brief Internal helper, adds the contributions from each particle in
    /// this cell to the map
    void AccumulateHitContributions(rb::CellHit const& hit,
                                    std::map<int, int>         & idToNPhot) const;

    Params fParams;

    bool fHaveTruthInfo; ///< Set by \ref Rebuild

    novadaq::cnv::DetId                    fDetID;                  ///< detector id

    std::unordered_map< geo::OfflineChan, std::list<sim::FLSHit> > fCellToFLSHit;    ///< map of FLSHits in each cell in event
    std::unordered_map< int,              std::list<sim::FLSHit> > fTrackIDToFLSHit; ///< map of G4 track ids to FLSHits
    std::unordered_map< geo::OfflineChan, std::vector<sim::PhotonSignal> > fCellToPhotonSignal; ///< map of PhotonSignals in each cell in event
    sim::ParticleNavigator                 fParticleNav;            ///< Particle navigator to map track ID to Particle
    std::set<int>                          fTrackIDs;               ///< G4 track ids for all particles in the events
    std::vector< art::Ptr<simb::MCTruth> > fMCTruthList;            ///< all the MCTruths for this event
    std::map<int, int>                     fTrackIDToMCTruthIndex;  ///< map of track id values to MCTruthList entry
    std::vector< art::Ptr<sim::TrueEnergy> > fTrueEnergyList;       ///< all the TrueEnergy's for this event
    std::map<int, int>                       fTrackIDToTrueEIndex;  ///< map of track id values to TrueEnergy entry

    /// map of cellHit to trackIDEs that contributed energy to the hit. 
    /// cellhit information is encoded as TDC*1000,000+Plane*1000+Cell
    std::map<double, std::vector<cheat::TrackIDE> > fCellToIDs;
    std::map<int, float > fIdToEnergy;

    int      fSliceIDForHitParticleMap;
    int      fNumTrueEnergyWarnings;
  };
} // namespace

DECLARE_ART_SERVICE(cheat::BackTracker, LEGACY)
#endif // CHEAT_BACKTRACKER_H