BackTracker_service.cc

Go to the documentation of this file.

//////////////////////////////////////////////////////////////////////
// $Id: BackTracker_service.cc,v 1.7 2012-12-06 15:03:39 lein Exp $
//
//
// \file: BackTracker_service.cc
//
// brebel@fnal.gov
//
////////////////////////////////////////////////////////////////////////

#include "MCCheater/BackTracker.h"

#include <map>
#include <algorithm>

// Framework includes
#include "canvas/Persistency/Common/FindOneP.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "Simulation/FLSHitList.h"
#include "Geometry/Geometry.h"

// need a map of CellUniqueId to vector of FLSHits to make this work
// really should just have a CellHit/RawDigit to provide the Cell and Plane
// and then look through the event to gett the proper FLSHits
// have to know where the FLSHits are stored in the event if want to pass art::Event around

namespace cheat{

  //----------------------------------------------------------------------
  static bool sortFLSHit(sim::FLSHit const& a, sim::FLSHit const& b)
  {
    return a.GetEdep() > b.GetEdep();
  }

  //----------------------------------------------------------------------
  bool CompareByEnergy(const TrackIDE& a, const TrackIDE& b)
  {
    return a.energy < b.energy;
  }

  //----------------------------------------------------------------------
  double EffMetric(const cheat::NeutrinoEffPur& ep)
  {
    return ep.efficiency;
  }

  //----------------------------------------------------------------------        
  double PurMetric(const cheat::NeutrinoEffPur& ep)
  {
    return ep.purity;
  }
  
  //----------------------------------------------------------------------
  double EffPurMetric(const cheat::NeutrinoEffPur& ep)
  {
    return ep.efficiency*ep.purity;
  }
  
  //----------------------------------------------------------------------
  double EnergyMetric(const cheat::NeutrinoEffPur& ep)
  {
    return ep.energySlice;
  }

  //----------------------------------------------------------------------
  BackTracker::BackTracker(const Parameters& params,
                           art::ActivityRegistry& reg)
    : fParams(params()), fHaveTruthInfo(false)
  {
    reg.sPreProcessEvent.watch(this, &BackTracker::Rebuild);
    fNumTrueEnergyWarnings = 0;
  }

  //----------------------------------------------------------------------
  BackTracker::~BackTracker()
  {
  }

  //----------------------------------------------------------------------
  void BackTracker::Rebuild(const art::Event& evt)
  {
    // Clear out anything remaining from previous calls to Rebuild
    fCellToFLSHit.clear();
    fTrackIDToFLSHit.clear();
    fCellToPhotonSignal.clear();
    fTrackIDs.clear();
    fParticleNav.clear();
    fMCTruthList.clear();
    fSliceIDForHitParticleMap = -1;
    fTrueEnergyList.clear();
    fTrackIDToTrueEIndex.clear();
    fHaveTruthInfo = false;               // Haven't got anything yet
    fDetID = novadaq::cnv::kUNKNOWN_DET;  // Haven't found detector type yet

    // do nothing if this is data
    if(evt.isRealData() && !fParams.MRMode() 
                        && !fParams.MCOverlayMode()) return;

    // get detector for CloseInTime function
    art::ServiceHandle<geo::Geometry>  geom;
    fDetID = geom->DetId();
    
    // fill the map of OfflineChan to std::vector<FLSHit>
    // first get the FLSHitList from the event
    art::Handle< std::vector<sim::FLSHitList> > flslistHandle;
    evt.getByLabel(fParams.GeantModuleLabel(), flslistHandle);
    
    // Probably this is a MC generation job and this is the preEvent callback
    // which is too early. That's OK, we'll get called at the right time by
    // ReadoutSim in a moment. Or we're mock data or something.
    if(flslistHandle.failedToGet()) return;

    fHaveTruthInfo = true; // We got FLSHits, the rest should go well too

    // loop over the FLSHits and fill the std::vector for each 
    // channel
    for( size_t l = 0; l < flslistHandle->size(); ++l){
      const std::vector<sim::FLSHit> &hits = flslistHandle->at(l).fHits;
      for( size_t f = 0; f < hits.size(); ++f){
        const sim::FLSHit& hit = hits[f];
        const geo::OfflineChan chan(hit.GetPlaneID(), hit.GetCellID());
        fCellToFLSHit[chan].push_back(hit);
        fTrackIDToFLSHit[hit.GetTrackID()].push_back(hit);
      }// end loop over flshits
    }// end loop over flshit lists
    
    // loop over the map and sort the vectors of flshits
    for(auto itr = fCellToFLSHit.begin(); itr != fCellToFLSHit.end(); itr++){
      itr->second.sort(sortFLSHit);
    }

    // loop over the map and sort the vectors of flshits
    for(auto pfitr = fTrackIDToFLSHit.begin(); pfitr != fTrackIDToFLSHit.end(); pfitr++){
      pfitr->second.sort(sortFLSHit);
    }

    // Similarly, PhotonSignals
    art::Handle<std::vector<sim::PhotonSignal> > phots;
    evt.getByLabel(fParams.PhotonModuleLabel(), phots);

    if(!phots.failedToGet()){
      for(size_t i = 0; i < phots->size(); ++i){
        const sim::PhotonSignal& phot = (*phots)[i];
        const geo::OfflineChan chan(phot.Plane(), phot.Cell());
        fCellToPhotonSignal[chan].push_back(phot);
      } // end for i

      for(auto it = fCellToPhotonSignal.begin(); it != fCellToPhotonSignal.end(); ++it){
        // Sort in reverse order, ie from most photons to least
        std::sort(it->second.rbegin(), it->second.rend(), sim::CompareByNPhoton);
      }
    }
    else{
      mf::LogWarning("BackTracker") << "unable to find PhotonTransport "
                                    << "created in " << fParams.PhotonModuleLabel() << " "
                                    << "any attempt to get the PhotonTransport objects from "
                                    << "the backtracker will fail";
    }

    art::Handle< std::vector<sim::Particle> > pHandle;
    evt.getByLabel(fParams.GeantModuleLabel(), pHandle);
    if(pHandle.failedToGet()){
      mf::LogWarning("BackTracker") << "unable to find sim::Particles "
                                    << "created in " << fParams.GeantModuleLabel() << " "
                                    << "any attempt to get the Particles from "
                                    << "the backtracker will fail";
    }
    else{
      // Keep track of all the MCTruths we discover associated to the particles
      std::set<art::Ptr<simb::MCTruth>> truths;
      // And which truth each particle points at
      std::map<int, art::Ptr<simb::MCTruth>> trackIDToMCTruthPtr;

      art::FindOneP<simb::MCTruth> fo(pHandle, evt, fParams.GeantModuleLabel());
      for(size_t p = 0; p < pHandle->size(); ++p){

        sim::Particle *part = new sim::Particle(pHandle->at(p));
        fParticleNav.Add(part);

        // get the simb::MCTruth associated to this vector<Particle>
        art::Ptr<simb::MCTruth> mct = fo.at(p);
        truths.insert(mct);
        try{
          trackIDToMCTruthPtr.emplace(pHandle->at(p).TrackId(), mct);
        }
        catch(cet::exception &ex){
          mf::LogWarning("BackTracker") << "unable to find MCTruth from Particle list "
                                        << "created in " << fParams.GeantModuleLabel() << " "
                                        << "any attempt to get the MCTruth objects from "
                                        << "the backtracker will fail\n"
                                        << "message from caught exception:\n" << ex;
        }
      } // end loop over particles to get MCTruthList

      // Flatten the MCTruth set into a vector
      fMCTruthList.insert(fMCTruthList.end(), truths.begin(), truths.end());

      // Translate all the truth pointers into indices
      std::map<art::Ptr<simb::MCTruth>, int> truthIdxs;
      for(unsigned int i = 0; i < fMCTruthList.size(); ++i){
        truthIdxs.emplace(fMCTruthList[i], i);
      }

      // Convert trackIDToMCTruthPtr from Ptrs to indices
      for(auto it: trackIDToMCTruthPtr){
        fTrackIDToMCTruthIndex[it.first] = truthIdxs[it.second];
      }
    }


    // fill the set of track ids from the ParticleNavigator
    fTrackIDs.clear();
    for(auto& it: fParticleNav) fTrackIDs.insert(it.first);

    LOG_DEBUG("BackTracker") << fParticleNav;

    // Get the sim::TrueEnergy data products out of the event. I don't think that I'll be able to utilise the Assn
    //  between sim::Particle and sim::TrueEnergy here, as I would have to somehow have to pass the Assn around.
    art::Handle< std::vector<sim::TrueEnergy> > teHandle;
    evt.getByLabel(fParams.TrueEnergyModuleLabel(), teHandle);
    if(!teHandle.failedToGet()) {
      // Get the TrueEnergy Assn, and loop through the particles. 
      art::FindManyP<sim::TrueEnergy> teAssn(pHandle, evt, fParams.TrueEnergyModuleLabel());
      for(size_t p = 0; p < pHandle->size(); ++p){ 
        int parTrId = pHandle->at(p).TrackId();
        // There should just be a 1 to 1 correspondance, however for some reason
        //  the first few particles may have multiple assns, and the last few none.
        bool match_truep = false;
        for(size_t te = 0; te < teAssn.at(p).size(); ++te){
          if (parTrId == teAssn.at(p)[te]->TrackId()) {
            match_truep = true;
            fTrueEnergyList.push_back   ( teAssn.at(p)[te] );
            fTrackIDToTrueEIndex.emplace( parTrId, p );     
          }
        }
        if (!match_truep) {
          // If didn't match, then it's one of the last particles, and so has an assn with a particle at the. 
          //  start of the particle loop. Lets restart the loop over the particles, pulling those assns...
          for(size_t par2 = 0; par2 < pHandle->size(); ++par2){
            for(size_t te = 0; te < teAssn.at(par2).size(); ++te){
              if (parTrId == teAssn.at(par2)[te]->TrackId()) {
                fTrueEnergyList.push_back   ( teAssn.at(par2)[te] );
                fTrackIDToTrueEIndex.emplace( parTrId, p );
                match_truep = true;
              }
              if (match_truep) break;
            }
            if (match_truep) break;
          }
        }
      }

      /*
        // I think that if using the true energy module you may have to use this technique.
        // Totally my fault, I made the assns in different way for some reason...
        art::FindOneP<sim::TrueEnergy> teAssn(pHandle, evt, fParams.TrueEnergyModuleLabel());
        for(size_t p = 0; p < pHandle->size(); ++p){ 
        art::Ptr<sim::TrueEnergy > MyEns = teAssn.at( p );
        fTrueEnergyList.push_back( MyEns );
        fTrackIDToTrueEIndex.emplace( pHandle->at(p).TrackId(), p );
        }
      */
    } else if (fNumTrueEnergyWarnings<10) {
      mf::LogWarning("BackTracker") << "unable to find sim::TrueEnergy "
                                    << "created in " << fParams.TrueEnergyModuleLabel() << " "
                                    << "any attempt to get the sim::TrueEnergy's from "
                                    << "the backtracker will fail";
      ++fNumTrueEnergyWarnings;
    }
  }

  //----------------------------------------------------------------------
  const std::vector<sim::PhotonSignal> BackTracker::CellToPhotonSignal(unsigned int const& plane,
                                                                       unsigned int const& cell) const
  {
    auto it = fCellToPhotonSignal.find(geo::OfflineChan(plane, cell));

    if(it == fCellToPhotonSignal.end()){
      LOG_DEBUG("BackTracker") << "no collection of sim::PhotonSignal for "
                               << " plane = " << plane << " cell = " << cell
                               << " probably this was a noise hit, returning an"
                               << " empty vector";
      const std::vector<sim::PhotonSignal> empty;
      return empty;
    }
    return it->second;
  }

  //----------------------------------------------------------------------
  const std::vector<sim::PhotonSignal> BackTracker::HitToPhotonSignal(const art::Ptr<rb::CellHit>& hit) const
  {
    auto it = fCellToPhotonSignal.find(geo::OfflineChan(hit->Plane(), hit->Cell()));
    
    // Nothing to accumulate from this cell
    if(it == fCellToPhotonSignal.end()){
      LOG_DEBUG("BackTracker") << "no collection of sim::PhotonSignal for "
                               << " plane = " << hit->Plane() << " cell = " << hit->Cell()
                               << " probably this was a noise hit, returning an"
                               << " empty vector";
      const std::vector<sim::PhotonSignal> empty;
      return empty;
    }
    
    // Candidate photons, in that they're in the right cell.
    const std::vector<sim::PhotonSignal>& phots = it->second;

    std::vector<sim::PhotonSignal> hitPhots;
    
    // For every photon that's sufficiently in time with the hit,
    // accumulate its energy into the relevant trackID.
    for(unsigned int i = 0; i < phots.size(); ++i){
      const sim::PhotonSignal& phot = phots[i];
      if(CloseInTime(*hit.get(), phot)){
        hitPhots.push_back(phot);
      }
    }
    return hitPhots;
  }
  
  //----------------------------------------------------------------------
  const std::vector<sim::FLSHit> BackTracker::CellToFLSHit(unsigned int const& plane,
                                                           unsigned int const& cell) const
  {
    auto it = fCellToFLSHit.find(geo::OfflineChan(plane, cell));

    if(it == fCellToFLSHit.end()){
      LOG_DEBUG("BackTracker") << "no collection of sim::FLSHit for "
                               << " plane = " << plane << " cell = " << cell
                               << " probably this was a noise hit, returning an"
                               << " empty vector";
      const std::vector<sim::FLSHit> empty;
      return empty;
    }
    return std::vector<sim::FLSHit>(it->second.begin(), it->second.end());
  }

  //----------------------------------------------------------------------
  std::vector<sim::FLSHit> BackTracker::HitToFLSHit(const rb::CellHit& hit) const
  {
    std::map<int, int> idToNPhot;
    AccumulateHitContributions(hit, idToNPhot);

    // Next, sort by number of photons. This turns out to be a bit fiddly.

    std::vector<std::pair<int, int> > nphotToId;
    nphotToId.reserve(idToNPhot.size());
    for(std::map<int, int>::iterator it = idToNPhot.begin(); it != idToNPhot.end(); ++it){
      nphotToId.push_back(std::make_pair(it->second, it->first));
    }

    if(nphotToId.empty()){
      LOG_DEBUG("BackTracker") << "No matching photons, this is probably a noise hit. Returning empty vector";
      return std::vector<sim::FLSHit>();
    }

    // Sort in reverse order (ie decreasing number of photons).
    // Tuples compare by first element first, so this works
    std::sort(nphotToId.rbegin(), nphotToId.rend());

    auto it = fCellToFLSHit.find(geo::OfflineChan(hit.Plane(), hit.Cell()));
    if ( it == fCellToFLSHit.end() ) {
      throw cet::exception("BackTracker")
        << "PhotonSignal with no matching FLSHit. That shouldn't be possible\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }

    std::vector<sim::FLSHit> ret;
    const std::list<sim::FLSHit>& flshits = it->second;

    ret.reserve(nphotToId.size()*flshits.size()); // Overestimate

    for(unsigned int idIdx = 0; idIdx < nphotToId.size(); ++idIdx){
      const int trackid = nphotToId[idIdx].second;
      for(const sim::FLSHit& fls: flshits){
        if(fls.GetTrackID() == trackid){
          ret.push_back(fls);
        }
      } // end for flsIdx
    } // end for idIdx

    return ret;
  }

  //----------------------------------------------------------------------
  const double BackTracker::EnergyFromTrackId(const art::Ptr<rb::CellHit>& hit, int trackId, bool useBirksE) const
  {
    auto it = fCellToFLSHit.find(geo::OfflineChan(hit->Plane(), hit->Cell()));
    if(it == fCellToFLSHit.end()){
      mf::LogWarning("BackTracker") << "No FLSHits in this cell";
      return 0;
    }

    const std::list<sim::FLSHit>& flshits = it->second;

    double energy = 0;

    for(const sim::FLSHit& fls: flshits){
      if(fls.GetTrackID() == trackId){
        energy += ((useBirksE) ? fls.GetEdepBirks() : fls.GetEdep());
      }
    } // end for flsIdx

    return energy;
  }


  //----------------------------------------------------------------------
  const sim::Particle* BackTracker::HitToParticle(const art::Ptr<rb::CellHit>& hit,
                                                  bool quiet) const
  {
    return HitToParticle( *hit.get(), quiet);
  }

  //----------------------------------------------------------------------
  const sim::Particle* BackTracker::HitToParticle(const rb::CellHit& hit,
                                                  bool quiet) const
  {
    std::map<int, int> idToNPhot;
    AccumulateHitContributions(hit, idToNPhot);

    // Go through the map and find the particle ID that contributed the most light
    int bestId = -1;
    int mostPhot = -1;
    for(std::map<int, int>::iterator it = idToNPhot.begin(); it != idToNPhot.end(); ++it){
      if(it->second > mostPhot){
        mostPhot = it->second;
        bestId = it->first;
      }
    }

    if(bestId == -1){
      if(!quiet){
        mf::LogWarning("BackTracker") << "No photons in this cell are in time with the hit. "
                                      << "Returning null pointer";
      }
      return 0;
    }

    return TrackIDToParticle(bestId);
  }

  //----------------------------------------------------------------------
  bool BackTracker::IsNoise(const art::Ptr<rb::CellHit>& hit) const
  {
    return (HitToParticle(hit, true) == 0);
  }

  //----------------------------------------------------------------------
  bool BackTracker::IsNoise(const rb::CellHit& hit) const
  {
    return (HitToParticle(hit, true) == 0);
  }

  //----------------------------------------------------------------------
  bool BackTracker::IsHitsNoise(std::vector< art::Ptr<rb::CellHit> > const& hits) const
  {
    int totalNoise = 0;
    int total = hits.size();

    for(size_t h = 0; h < hits.size(); ++h){
      art::Ptr<rb::CellHit> hit = hits[h];
      if (IsNoise(hit)) totalNoise += 1;
    }

    int totalPhysics = total - totalNoise;
    double physicsFrac = (double)totalPhysics / (double)total;

    // If at least MinPhysicsHits physics hits in the collection, not a noise group. 
    // If at least MinPhysicsFrac of the hits are physics, not a noise group. 
    // Otherwise, call it a noise group.

    return((totalPhysics < fParams.MinPhysicsHits()) && (physicsFrac < fParams.MinPhysicsFrac()));
  }


  //----------------------------------------------------------------------
  const sim::Particle* BackTracker::TrackIDToParticle(int const& id) const
  {
    sim::ParticleNavigator::const_iterator part_it = fParticleNav.find(abs(id));

    if(part_it == fParticleNav.end()){
      mf::LogWarning("BackTracker") << "can't find particle with track id "
                                    << id << " in sim::ParticleNavigator"
                                    << " returning null pointer";
      return 0;
    }

    return part_it->second;
  }

  //----------------------------------------------------------------------
  const sim::Particle* BackTracker::TrackIDToMotherParticle(int const& id) const
  {
    // get the mother id from the particle navigator
    // the EveId was adopted in the Rebuild method
 
    return this->TrackIDToParticle(fParticleNav.EveId(abs(id)));
  }

  //----------------------------------------------------------------------
  const art::Ptr<simb::MCTruth>& BackTracker::TrackIDToMCTruth(int const& id) const
  {
    // find the entry in the MCTruth collection for this track id
    size_t mct = fTrackIDToMCTruthIndex.find(abs(id))->second;
    
    if (mct > fMCTruthList.size() ) {
      throw cet::exception("BackTracker")
        << "attempting to find MCTruth index for "
        << "out of range value: " << mct
        << "/" << fMCTruthList.size() << "\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }
    return fMCTruthList[mct];
  }

  //----------------------------------------------------------------------
  const art::Ptr<simb::MCTruth>& BackTracker::ParticleToMCTruth(const sim::Particle* p) const
  {
    return this->TrackIDToMCTruth(p->TrackId());
  }

  //----------------------------------------------------------------------
  std::vector<const sim::Particle*> BackTracker::MCTruthToParticles(art::Ptr<simb::MCTruth> const& mct) const
  {
    std::vector<const sim::Particle*> ret;

    for(std::set<int>::const_iterator itr = fTrackIDs.begin(); itr != fTrackIDs.end(); itr++){
      if( this->TrackIDToMCTruth(*itr) == mct ) ret.push_back(fParticleNav.find(*itr)->second);
    }

    return ret;
  }

  //----------------------------------------------------------------------
  // this method orders the vector of particles from contributing most
  // to the collection to least
  std::vector<const sim::Particle*> BackTracker::
  HitsToParticle(const std::vector<const rb::CellHit*>& hits) const
  {
    std::map<int, int> idToNPhot;
    for(size_t h = 0; h < hits.size(); ++h){
      AccumulateHitContributions(*hits[h], idToNPhot);
    }

    // Now "all" we need to do is sort by number of photons, and discard any IDs
    // that we can't find particles for. This turns out to be a bit fiddly.

    std::vector<std::pair<int, int> > nphotToId;
    for(std::map<int, int>::iterator it = idToNPhot.begin(); it != idToNPhot.end(); ++it){
      nphotToId.push_back(std::make_pair(it->second, it->first));
    }

    // Sort in reverse order (ie increasing number of photons).
    // Tuples compare by first element first, so this works
    std::sort(nphotToId.rbegin(), nphotToId.rend());

    std::vector<const sim::Particle*> ret;
    for(unsigned int n = 0; n < nphotToId.size(); ++n){
      const sim::Particle* part = TrackIDToParticle(nphotToId[n].second);
      if(part) ret.push_back(part);
    }

    return ret;
  }

  //----------------------------------------------------------------------
  // this method orders the vector of particles from contributing most
  // to the collection to least
  std::vector<TrackIDE> BackTracker::
  HitsToTrackIDE(const std::vector<const rb::CellHit*>& hits, bool useBirksE) const
  {
    
    std::map<int, double> sigtracks;

    double totalE = 0;
    for(size_t h = 0; h < hits.size(); ++h){
      const std::vector<sim::FLSHit> fls = HitToFLSHit(*hits[h]);
      for(size_t f = 0; f < fls.size(); ++f){
        if( useBirksE){
          sigtracks[fls[f].GetTrackID()] += fls[f].GetEdepBirks();
          totalE += fls[f].GetEdepBirks();
        }
        else{
          sigtracks[fls[f].GetTrackID()] += fls[f].GetEdep();
          totalE += fls[f].GetEdep();
        }
      }
    } // end for h

    // Make trackIDE objects from our collected energies
    std::vector<TrackIDE> trackIDE;
    for(std::map<int, double>::iterator it = sigtracks.begin(); it != sigtracks.end(); ++it){
      TrackIDE ide;
      ide.trackID = it->first;
      ide.energy = it->second;
      ide.energyFrac = it->second/totalE;
      trackIDE.push_back(ide);
    }

    // Sort TrackIDE objects from most energy contributed to least
    std::sort(trackIDE.rbegin(), trackIDE.rend(), CompareByEnergy);

    return trackIDE;
  }

  //----------------------------------------------------------------------
  std::vector<TrackIDE> BackTracker::
  HitToTrackIDE(const art::Ptr<rb::CellHit>& hit, bool useBirksE) const
  {
    return HitsToTrackIDE(std::vector<const rb::CellHit*>(1, hit.get()), useBirksE);
  }

  //----------------------------------------------------------------------
  std::vector<TrackIDE> BackTracker::
  HitToTrackIDE(const rb::CellHit& hit, bool useBirksE) const
  {
    return HitsToTrackIDE(std::vector<const rb::CellHit*>(1, &hit), useBirksE);
  }

  //----------------------------------------------------------------------
  TVector3 BackTracker::HitToXYZ(art::Ptr<rb::CellHit> const& hit, bool useBirksE) const
  {
    double x = 0.;
    double y = 0.;
    double z = 0.;
    double w = 0.;

    const std::vector<sim::FLSHit> fhits = this->HitToFLSHit(hit);
    
    // loop over the FLSHits for this cell and get the weighted central 
    // position of each track id in it.
    for(const sim::FLSHit& fhit: fhits){
      double weight = ((useBirksE) ? fhit.GetEdepBirks() : fhit.GetEdep());
      
      w += weight;
      x += weight * 0.5*(fhit.GetExitX() + fhit.GetEntryX());
      y += weight * 0.5*(fhit.GetExitY() + fhit.GetEntryY());
      z += weight * 0.5*(fhit.GetExitZ() + fhit.GetEntryZ());

    }// end loop over sim::FLSHits
        
    // if the sum of the weights is still 0, then return
    // the obviously stupid default values
    if(w < 1.e-5) return TVector3(-999, -999, -999);

    return TVector3(x/w, y/w, z/w);
  }

  //----------------------------------------------------------------------
  double BackTracker::
  HitCollectionPurity(const std::set<int>&                trackIDs,
                      const std::vector<rb::WeightedHit>& whits,
                      std::map<int, double>*              purMap,
                      std::map<int, int>*                 parents,
                      bool                                energyPur) const
  {
    const std::vector<RawWeightedHit> rwh(whits.begin(), whits.end());
    return HitCollectionPurity(trackIDs, rwh, purMap, parents, energyPur);
  }

  //----------------------------------------------------------------------
  double BackTracker::
  HitCollectionPurity(const std::set<int>&               trackIDs,
                      const std::vector<RawWeightedHit>& whits,
                      std::map<int, double>*             purMap,
                      std::map<int, int>*                parents,
                      bool                               energyPur) const

  {
    // get the list of EveIDs that correspond to the hits in this collection
    // if the EveID shows up in the input list of trackIDs, then it counts
    double total   = whits.size();
    double desired = 0;

    double desiredE = 0;
    double totalE   = 0;

    if(purMap) purMap->clear();

    for(size_t h = 0; h < whits.size(); ++h){

      double totalEInCell   = 0;
      double desiredEInCell = 0;

      const rb::CellHit* hit = whits[h].hit;
      double weight = whits[h].weight;

      assert( weight <= 1);

      if(IsNoise(*hit)) continue;

      std::vector<TrackIDE> tracks = HitToTrackIDE(*hit);

      // The set of trackIDs (or parents) that contributed to this hit
      std::set<int> tally;

      for(size_t t = 0; t < tracks.size(); ++t){
        const int trackID = tracks[t].trackID;
        if(trackIDs.find(trackID) != trackIDs.end()){
          tally.insert(trackID);
        }
        if(parents){
          std::map<int, int>::iterator parentsIt = parents->find(trackID);
          if(parentsIt != parents->end()){
            tally.insert(parentsIt->second);
          }
        }
      } // end for t


      if (energyPur){
        for(size_t t = 0; t < tracks.size(); ++t){
          totalEInCell += tracks[t].energy * weight;
          if(trackIDs.find(tracks[t].trackID) != trackIDs.end()){
            desiredE += tracks[t].energy;
            desiredEInCell += tracks[t].energy;
          }
        }
      }//End of loop for counting energy when doing purity based on FLS hits

      // Purity in a cell can sometimes be greater than one if the weight of
      // a hit is less than the fraction of energy from desired track IDs
      // in that cell.
      if( desiredEInCell >= totalEInCell )
        totalE += desiredEInCell;
      else
        totalE += totalEInCell;

      if(purMap){
        for(std::set<int>::iterator it = tally.begin(); it != tally.end(); ++it){
          ++(*purMap)[*it];
        }
      }

      if(!tally.empty()) ++desired;


    }// end loop over whits

    double purity = 0.;

    if((total  > 0) && !energyPur) purity = desired/total;
    if((totalE > 0) &&  energyPur) purity = desiredE/totalE;
    
    assert(purity >= 0);
    assert(purity <= 1);

    if(total > 0 && purMap){
      for(std::map<int, double>::iterator it = purMap->begin(); it != purMap->end(); ++it){
        it->second /= total;
        // This is a purity, it should be between 0-100%                                                               
        assert(it->second >= 0);
        assert(it->second <= 1);
      }
    }

    return purity;

  }
  

  //----------------------------------------------------------------------
  

  double BackTracker::
  HitCollectionPurity(const std::set<int>&               trackIDs,
                      const std::vector<const rb::CellHit*>& hits,
                      std::map<int, double>*             purMap,
                      std::map<int, int>*                parents,
                      bool                               energyPur) const
  {
    // convert all cell hits to weighted hits with weight 1.
    std::vector<RawWeightedHit> wHits;
    for(const rb::CellHit* h: hits) wHits.emplace_back(h, 1);

    return HitCollectionPurity(trackIDs, wHits, purMap, parents, energyPur);

  }

  //----------------------------------------------------------------------
  double BackTracker::
  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,
                          bool                   energyEff,
                          double*                desiredEnergy,
                          double*                totalEnergy,
                          int*                   desiredHits,
                          int*                   totHits) const
  {

    // convert cell hits to weighted hits with weight 1
    std::vector<RawWeightedHit> wHits;
    for(const rb::CellHit* h: hits) wHits.emplace_back(h, 1);

    return HitCollectionEfficiency(trackIDs, wHits,  allhits, view,
                                   effMap, energyEff, desiredEnergy,
                                   totalEnergy, desiredHits, totHits);
  }
  //----------------------------------------------------------------------
  double BackTracker::
  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,
                          bool                                energyEff,
                          double*                             desiredEnergy,
                          double*                             totalEnergy,
                          int*                                desiredHits,
                          int*                                totHits) const
  {
    const std::vector<RawWeightedHit> rwh(whits.begin(), whits.end());
    return HitCollectionEfficiency(trackIDs, rwh, allhits, view,
                                   effMap, energyEff, desiredEnergy, totalEnergy, desiredHits, totHits);
  }

  //----------------------------------------------------------------------
  double BackTracker::
  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,
                          bool                               energyEff,
                          double*                            desiredEnergy,
                          double*                            totalEnergy,
                          int*                               desiredHits,
                          int*                               totHits) const

  {

    // allhits should be a superset of hits
    assert(allhits.size() >= whits.size());

    // to check it's a superset more thoroughly, optional
    if(fParams.DoHitsCheck()){
      for(unsigned int n = 0; n < whits.size(); ++n){
        // TODO: should be checking actual equality not the pointers
        // assert(std::find(allhits.begin(), allhits.end(), whits[n].hit) != allhits.end());
      }
    }
    std::vector<const rb::CellHit*> hits_copy;
    for(int i = 0; i <  (int)whits.size(); i++)
      hits_copy.push_back(whits[i].hit);

    // Check for duplicates
    std::sort(hits_copy.begin(), hits_copy.end());
    for(unsigned int n = 0; n+1 < hits_copy.size(); ++n)
      assert(hits_copy[n] < hits_copy[n+1]);
    
    // get the list of EveIDs that correspond to the hits in this collection
    // and the energy associated with the desired trackID
    
    double desired = 0.;
    double total   = 0.;
    
    double desiredE = 0.;
    double totalE   = 0.;

    if(effMap) effMap->clear();
    std::map<int, double> totMap;

    for(size_t h = 0; h < whits.size(); ++h){

      double totalEInCell   = 0.;
      double desiredEInCell = 0.;

      const rb::CellHit* hit = whits[h].hit;
      double weight = whits[h].weight;
      assert( weight <= 1);
      
      if(IsNoise(*hit)) continue;
      
      bool hitCounted = false;
      
      std::vector<TrackIDE> tracks = HitToTrackIDE(*hit);
      
      assert(hit->View() == view || view == geo::kXorY);
      
      // don't double count if this hit has more than one of the
      // desired track IDs associated with it
      for(size_t t = 0; t < tracks.size(); ++t){
        totalEInCell += tracks[t].energy;
        if(trackIDs.find(tracks[t].trackID) != trackIDs.end() ){
          desiredEInCell += tracks[t].energy;
          if(tracks[t].energyFrac >= fParams.MinContribFrac()){
            if(!hitCounted){
              desired += 1.;   
              hitCounted = true;
            }
            if(effMap) ++(*effMap)[tracks[t].trackID];
          }
        }
      }// end loop over track IDEs

      if( (weight*totalEInCell) >= desiredEInCell)
        desiredE += desiredEInCell;
      else
        desiredE += (weight*totalEInCell);

    }// end loop over hits

    // now figure out how many hits in the whole collection are associated with this id
    for(size_t h = 0; h < allhits.size(); ++h){
      const rb::CellHit* hit = allhits[h];
      if(IsNoise(*hit)) continue;

      std::vector<TrackIDE> tracks = HitToTrackIDE(*hit);

      bool hitCounted = false;

      // check that we are looking at the appropriate view here
      // in the case of 3D objects we take all hits
      if(hit->View() != view && view != geo::kXorY) continue;
      
      for(size_t t = 0; t < tracks.size(); ++t){
        // don't worry about hits where the energy fraction for the chosen     
        // trackID is < MinContribFrac
        // also don't double count if this hit has more than one of the   
        // desired track IDs associated with it                           
        if(trackIDs.find(tracks[t].trackID) != trackIDs.end()){
          totalE += tracks[t].energy;
          if(tracks[t].energyFrac >= fParams.MinContribFrac()){
            if(!hitCounted){
              total += 1.;
              hitCounted = true;
            }
            if(effMap) ++totMap[tracks[t].trackID];
          }// end if energy fraction is above threshold
        }// end if the current track ID is in the list of ones to check
      }// end loop over tracks
    }// end loop over all hits

    double efficiency = 0.;
    if((total  > 0.) && !energyEff)     efficiency = desired/total;
    if((totalE > 0.) &&  energyEff)     efficiency = desiredE/totalE;

    if((total  > 0.) && desiredHits  )  *desiredHits   = desired;
    if((total  > 0.) && totHits      )  *totHits       = total;
    if((totalE > 0.) && desiredEnergy)  *desiredEnergy = desiredE;
    if((totalE > 0.) && totalEnergy  )  *totalEnergy   = totalE;

    if(effMap){
      for(std::map<int, double>::iterator it = effMap->begin(); it != effMap->end(); ++it){
        if(totMap[it->first] > 0) 
          it->second /= totMap[it->first];
        // This is an efficiency. It should be between 0-100%                                                           
        assert(it->second >= 0);
        assert(it->second <= 1);
      }
    }
    
    return efficiency;
  }


  //----------------------------------------------------------------------
  bool BackTracker::PassMuonCuts(int                                trackID,
                                 art::PtrVector<rb::CellHit> const& hits) const
  {
    int nxcells = 0;
    int nycells = 0;

    for(size_t h = 0; h < hits.size(); ++h){

      art::Ptr<rb::CellHit> hit = hits[h];
      std::vector<TrackIDE> tracks = HitToTrackIDE(hit);
      for(size_t t = 0; t < tracks.size(); ++t){
        if(tracks[t].trackID == trackID){

          if (hit->View() == geo::kX){
            nxcells++;
          }// end of loop over x hits

          if (hit->View() == geo::kY){
            nycells++;
          }// end of loop over y hits
          
        }// end of loop over hits from the muon
      }// end of loop over tracks
    }// end loop over hits


    if ((nxcells >= 3) && (nycells >= 3)) return 1;
    else return 0;
  }
    

  //----------------------------------------------------------------------
  bool BackTracker::sort_eff(const NeutrinoEffPur& a, const NeutrinoEffPur& b){
    return a.efficiency>b.efficiency;
  }

  //----------------------------------------------------------------------
  bool BackTracker::sort_pur(const NeutrinoEffPur& a, const NeutrinoEffPur& b){
    return a.purity>b.purity;
  }
  //----------------------------------------------------------------------
  bool BackTracker::sort_idxeff(const std::pair<int,NeutrinoEffPur>& a, const std::pair<int,NeutrinoEffPur>& b){
    return a.second.efficiency>b.second.efficiency;
  }
  //----------------------------------------------------------------------
  bool BackTracker::sort_idxpur(const std::pair<int,NeutrinoEffPur>& a, const std::pair<int,NeutrinoEffPur>& b){
    return a.second.purity>b.second.purity;
  }

  //----------------------------------------------------------------------
  std::vector<int> BackTracker::
  SliceToNeutrinoIndex(const std::vector<const rb::CellHit*>& sliceHits,
                       const std::vector<const rb::CellHit*>& allHits,
                       bool                               sortPur) const
   {
     std::vector < std::pair<unsigned int,NeutrinoEffPur> > sliceNeutrinoInteractions;

     for (unsigned int s = 0; s < fMCTruthList.size(); s++){
       const art::Ptr<simb::MCTruth> neutrinoInteraction = fMCTruthList[s];
       if(!neutrinoInteraction->NeutrinoSet())
         continue;
       
       std::vector<const sim::Particle*> particleList = MCTruthToParticles(neutrinoInteraction);
       std::set<int> neutrinoTrackIDs;
       for (unsigned int i = 0; i < particleList.size(); ++i){
         neutrinoTrackIDs.insert(particleList[i]->TrackId());
       }
       
       double energySlice = -1.0;
       double energyTotal = -1.0;
       int    nSliceHits  = 0;
       int    nTotalHits  = 0;
       
       double sliceEff = HitCollectionEfficiency(neutrinoTrackIDs, sliceHits, allHits, geo::kXorY, 0, true, &energySlice, &energyTotal, &nSliceHits, &nTotalHits);
       double slicePur = HitCollectionPurity(neutrinoTrackIDs, sliceHits, 0, 0, true);
       
       if (sliceEff > 0){
         NeutrinoEffPur sliceNeutrinoInteraction = {neutrinoInteraction, sliceEff, slicePur, s, energySlice, energyTotal, nSliceHits, nTotalHits};
         sliceNeutrinoInteractions.push_back( std::make_pair(s,sliceNeutrinoInteraction) );
       }
     }//End of loop over neutrino interactions 
     
     if (!sortPur) std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),sort_idxeff);
     else          std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),sort_idxpur);
     
     std::vector<int> nuIndex;

     for (unsigned int i = 0; i<sliceNeutrinoInteractions.size(); ++i ){
       nuIndex.push_back( int(sliceNeutrinoInteractions[i].first) );
     }
     
    return nuIndex;
  }
  

  //----------------------------------------------------------------------
  std::vector<NeutrinoEffPur> BackTracker::
  SliceToNeutrinoInteractions(const std::vector<const rb::CellHit*>& sliceHits,
                              const std::vector<const rb::CellHit*>& allHits,
                              bool                               sortPur) const
   {
    std::vector <NeutrinoEffPur> sliceNeutrinoInteractions;
    
    for (unsigned int s = 0; s < fMCTruthList.size(); s++){
      const art::Ptr<simb::MCTruth> neutrinoInteraction = fMCTruthList[s];
      if(!neutrinoInteraction->NeutrinoSet())
        continue;
      
      std::vector<const sim::Particle*> particleList = MCTruthToParticles(neutrinoInteraction);
      std::set<int> neutrinoTrackIDs;
      for (unsigned int i = 0; i < particleList.size(); ++i){
        neutrinoTrackIDs.insert(particleList[i]->TrackId());
      }

      double energySlice = -1.0;
      double energyTotal = -1.0;
      int    nSliceHits  = 0;
      int    nTotalHits  = 0;

      double sliceEff = HitCollectionEfficiency(neutrinoTrackIDs, sliceHits, allHits, geo::kXorY, 0, true, &energySlice, &energyTotal, &nSliceHits, &nTotalHits);
      double slicePur = HitCollectionPurity(neutrinoTrackIDs, sliceHits, 0, 0, true);

      if (sliceEff > 0){
        NeutrinoEffPur sliceNeutrinoInteraction = {neutrinoInteraction, sliceEff, slicePur, s, energySlice, energyTotal, nSliceHits, nTotalHits};
        sliceNeutrinoInteractions.push_back(sliceNeutrinoInteraction);
      }
    }//End of loop over neutrino interactions 

    if (!sortPur) std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),sort_eff);
    else          std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),sort_pur);

    return sliceNeutrinoInteractions;
  }
  
  //----------------------------------------------------------------------
  std::vector<NeutrinoEffPur> BackTracker::
  SliceToNeutrinoInteractions(const art::Ptr<rb::Cluster>& sliceCluster,
                              const std::vector<const rb::CellHit*>& allHits,
                              bool sortPur) const
  {
    return SliceToNeutrinoInteractions(PtrVecToVecRawPtr(sliceCluster->AllCells()),
                                       allHits, sortPur);
  }
  
  //......................................................................
  std::vector<int> BackTracker::
  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
  {

    for( auto slc: slTruthTable )
      assert( slc.size() == nusWithIdx.size() && "Badly shaped Slice/truth table!" );
    
    // Fill with one Nu Id per slice
    std::vector<int> nuIds( slTruthTable.size(), -1 );
    if( nusWithIdx.empty() ) return nuIds;

    // Find the best slice for each nu
    std::vector<int> bestSlices( nusWithIdx.size() );
    for( size_t nuIdx = 0; nuIdx < nusWithIdx.size(); ++nuIdx ){
      float maxValue = 0;
      int  bestSlice = -1;

      for( size_t sliceIdx = 0; sliceIdx < slTruthTable.size(); ++sliceIdx ){
        const cheat::NeutrinoEffPur& effpur = slTruthTable[sliceIdx][nuIdx];

        if( !effpur.neutrinoInt.isNull() ){
          if( slMetric(effpur) > maxValue ){
            maxValue  = slMetric(effpur);
            bestSlice = sliceIdx;
          }
        }

      } //sliceIdx

      bestSlices[nuIdx] = bestSlice;

    } // nuIdx

    // Find the best nu for each slice.
    for( size_t sliceIdx = 0; sliceIdx < slTruthTable.size(); ++sliceIdx ){
      float maxValue = 0;
      int   bestNu   = -1;

      for( size_t nuIdx = 0; nuIdx < slTruthTable[sliceIdx].size(); ++nuIdx ){
        const cheat::NeutrinoEffPur& effpur = slTruthTable[sliceIdx][nuIdx];

        if( !effpur.neutrinoInt.isNull() ){
          if( nuMetric(effpur) > maxValue ){
            maxValue = nuMetric(effpur);
            bestNu   = nuIdx;
          }
        }
      } //nuIdx

      if( bestNu >= 0 && bestSlices[bestNu] == (int)sliceIdx )
        nuIds[sliceIdx] = bestNu;

    }//sliceIdx

    return nuIds;
  }

  //----------------------------------------------------------------------
  std::vector<int> BackTracker::
  SliceToOrderedNuIdsByEff(const std::vector<cheat::NeutrinoWithIndex>& nusWithIdx,
                           const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable) const
  {

    return SliceToOrderedNuIds(nusWithIdx, slTruthTable, EffMetric, EffMetric);
 
  }

  //----------------------------------------------------------------------
  std::vector<int> BackTracker::
  SliceToOrderedNuIdsByPur(const std::vector<cheat::NeutrinoWithIndex>& nusWithIdx,
                           const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable) const
  {

    return SliceToOrderedNuIds(nusWithIdx, slTruthTable, PurMetric, PurMetric);
 
  }

  //----------------------------------------------------------------------
  std::vector<int> BackTracker::
  SliceToOrderedNuIdsByEffPur(const std::vector<cheat::NeutrinoWithIndex>& nusWithIdx,
                           const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable) const
  {

    return SliceToOrderedNuIds(nusWithIdx, slTruthTable, EffPurMetric, EffPurMetric);
 
  }

  //----------------------------------------------------------------------
  std::vector<int> BackTracker::
  SliceToOrderedNuIdsByEnergy(const std::vector<cheat::NeutrinoWithIndex>& nusWithIdx,
                              const std::vector<std::vector<cheat::NeutrinoEffPur>>& slTruthTable) const
  {

    return SliceToOrderedNuIds(nusWithIdx, slTruthTable, EnergyMetric, EnergyMetric);
 
  }

  //----------------------------------------------------------------------
  std::vector<NeutrinoEffPur> BackTracker::
  SliceToMCTruth(const std::vector<const rb::CellHit*>& sliceHits,
                 const std::vector<const rb::CellHit*>& allHits,
                 bool sortPur) const
  {
    std::vector <NeutrinoEffPur> sliceNeutrinoInteractions;
    
    for (unsigned int s = 0; s < fMCTruthList.size(); s++){
      const art::Ptr<simb::MCTruth> neutrinoInteraction = fMCTruthList[s];
      
      std::vector<const sim::Particle*> particleList = MCTruthToParticles(neutrinoInteraction);
      std::set<int> neutrinoTrackIDs;
      for (unsigned int i = 0; i < particleList.size(); ++i){
        neutrinoTrackIDs.insert(particleList[i]->TrackId());
      }
      
      double energySlice = -1.0;
      double energyTotal = -1.0;
      int    nSliceHits  = 0;
      int    nTotalHits  = 0;

      double sliceEff = HitCollectionEfficiency(neutrinoTrackIDs, sliceHits, allHits, geo::kXorY, 0, true, &energySlice, &energyTotal, &nSliceHits, &nTotalHits);
      double slicePur = HitCollectionPurity(neutrinoTrackIDs, sliceHits, 0, 0, true);
      
      if (sliceEff > 0){
        NeutrinoEffPur sliceNeutrinoInteraction = {neutrinoInteraction, sliceEff, slicePur, s, energySlice, energyTotal, nSliceHits, nTotalHits};
        sliceNeutrinoInteractions.push_back(sliceNeutrinoInteraction);
      }
    }//End of loop over neutrino interactions 
    
    if (!sortPur) std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),sort_eff);
    else          std::sort(sliceNeutrinoInteractions.begin(),sliceNeutrinoInteractions.end(),sort_pur);
    
    return sliceNeutrinoInteractions;
  }
  
  //----------------------------------------------------------------------
  std::vector<NeutrinoEffPur> BackTracker::
  SliceToMCTruth(const art::Ptr<rb::Cluster>& sliceCluster,
                 const std::vector<const rb::CellHit*>& allHits,
                 bool sortPur) const
  {
    return SliceToMCTruth(PtrVecToVecRawPtr(sliceCluster->AllCells()),
                          allHits, sortPur);
  }

  //----------------------------------------------------------------------
  std::vector< std::vector<cheat::NeutrinoEffPur> > BackTracker::
  SlicesToMCTruthsTable(const std::vector<const rb::Cluster*>& sliceList) const
  {

    unsigned int numSlices = sliceList.size();
    std::vector<NeutrinoWithIndex> truthList = allMCTruth();
    unsigned int numTruths = truthList.size();

    //This is what will eventually be returned
    std::vector< std::vector<cheat::NeutrinoEffPur> > table(numSlices,std::vector<cheat::NeutrinoEffPur>(numTruths));
    
    //Loop to fill NeutrinoEffPur's with proper truth information
    for (unsigned int s = 0; s < numSlices; s++){
      for (unsigned int t = 0; t < numTruths; t++){
        cheat::NeutrinoEffPur truthFill = {truthList[t].neutrinoInt,0.0,0.0,truthList[t].truthColIndex,0.0,0.0,0,0};
        table[s][t] = truthFill; 
      }
    }//End of defining basic container with truths filled in

    for (unsigned int s = 0; s < numSlices; s++){
      //Vector to hold the total energy contributions to slice from each MC Truth
      std::vector<double> sliceEnergies(numTruths,0.0);
      //Vector to hold the total number of hits in the slice from each MC Truth
      std::vector<int> sliceHits(numTruths,0);

      for (unsigned int h = 0; h < sliceList[s]->NCell(); h++){
        //Don't worry about noise hits - they can't count towards the energy an interaction deposited
        if (IsNoise(sliceList[s]->Cell(h))) continue;
        //Vector to hold the total energy contributions to hit from each MC Truth
        std::vector<double> hitEnergies(numTruths,0.0);

        const std::vector<TrackIDE> trackIDs = HitToTrackIDE(sliceList[s]->Cell(h));
        for (unsigned int id = 0; id < trackIDs.size(); id++){
          int trackID = trackIDs[id].trackID;
          //Gets the MCTruth index for the particle trackID
          int truthID = fTrackIDToMCTruthIndex.find(abs(trackID))->second;
          hitEnergies[truthID] += trackIDs[id].energy;
        }//End of loop over track ids for the hit   
        
        //Add hit values into the slice vectors
        for (unsigned int q = 0; q < numTruths; q++){  
          if (hitEnergies[q] > 0.0){
            sliceEnergies[q] += hitEnergies[q];
            sliceHits[q] += 1;
          }//Only do for MCTruths with nonzero contributions to hit
        }//End of loop adding hit values to slice vectors  
   
      }//End of loop over hits in slice

      double totalSliceE = std::accumulate(sliceEnergies.begin(),sliceEnergies.end(),0.0);
      for (unsigned int t = 0; t < numTruths; t++){

        double purity = 0.0;
        if (totalSliceE > 0){
          purity = (sliceEnergies[t])/totalSliceE;
        }

        if ((purity < 0) || (purity > 1)) {
        throw cet::exception("BackTrackerIllegalPurity")
          << "The purity value "<<purity<<" is not allowed. \n";
      }

        table[s][t].purity = purity;
        table[s][t].energySlice = sliceEnergies[t];
        table[s][t].nSliceHits = sliceHits[t];

        }//End of loop over MCTruths for the slice
    }//End of loop over slices

    for (unsigned int t = 0; t < numTruths; t++){

      double totalMCTruthE = 0.0;
      int totalMCTruthHits = 0;
      
      for (unsigned int s = 0; s < numSlices; s++){
        totalMCTruthE += table[s][t].energySlice;
        totalMCTruthHits += table[s][t].nSliceHits;
      }//Loop over slices to find summed values

      for (unsigned int s = 0; s < numSlices; s++){

        double efficiency = 0.0;
        if (totalMCTruthE > 0){
          efficiency = (table[s][t].energySlice)/totalMCTruthE;
        }

        if ((efficiency < 0) || (efficiency > 1)) {
        throw cet::exception("BackTrackerIllegalEfficiency")
          << "The efficiency value "<<efficiency<<" is not allowed. \n";
      }

        table[s][t].efficiency = efficiency;
        table[s][t].energyTotal = totalMCTruthE;
        table[s][t].nTotalHits = totalMCTruthHits;
      }//Loop over slices to insert final values
      
    }//End of final loop over MCTruths for summed values
    
    return table;
  }

  //----------------------------------------------------------------------
  std::vector<ParticleEffPur>  BackTracker::
  TracksToParticles(const std::vector<const rb::Track*>& tracks,
                    const std::set<int>& trackIDs,
                    const std::set<int>& allowedDaughterIDs,
                    const std::vector<const rb::CellHit*>& allHits) const
  {
    std::vector<ParticleEffPur> matches;

    std::set<int> trkIDs = trackIDs; 
    std::set<int> allowedDaughters = allowedDaughterIDs;

    std::map<int, int> parents;
    for(std::set<int>::iterator it = trkIDs.begin(); it != trkIDs.end(); ++it){
      const sim::Particle* mother = TrackIDToParticle(*it);
      if(!mother) continue;
      for(int i = 0; i < mother->NumberDaughters(); ++i){
        const int id = mother->Daughter(i);
        assert(TrackIDToParticle(id));
        if(allowedDaughters.count(TrackIDToParticle(id)->PdgCode())){
          parents[id] = *it;
        }
      } // end for i
    } // end for it

      // loop over all tracks
    for(unsigned int itrack = 0; itrack < tracks.size(); ++itrack){
      // get the purity of the track for all track IDs
      std::map<int, double> purs;
      HitCollectionPurity(trkIDs,tracks[itrack]->AllCells(),&purs,&parents);
      // find the most pure trackID for this track
      int bestid = -1;
      double bestpur = 0.0;
      for(std::set<int>::iterator it = trkIDs.begin(); it != trkIDs.end(); ++it){
        const double p = purs[*it];
        if(p > bestpur){
          bestpur = p;
          bestid = *it;
        }
      } // end of loop over track ids

      // update this tracks info for the matched true particle      
      double recoPur = bestpur;
      int recotrkid = bestid;
      int recopdg = 0;
      double recoEff = 0.0;
      if(bestpur > 0.0){
        // found a real particle that matched the track, update pdg and efficiency to match this
        recopdg = TrackIDToParticle(recotrkid)->PdgCode();
        std::set<int> ideff;
        ideff.insert(recotrkid);
        recoEff = HitCollectionEfficiency(ideff,
                                          PtrVecToVecRawPtr(tracks[itrack]->AllCells()),
                                          allHits, geo::kXorY);
        // remove this track id from the set of track ids so no other tracks can try to match to it
        trkIDs.erase(trkIDs.find(recotrkid));
        // remove any instance where this track id is the parent particle in the parents map
        std::map<int, int>::iterator parentsIt = parents.begin();
        while(parentsIt != parents.end()){
          if(parentsIt->second == recotrkid){ parents.erase(parentsIt++); }
          else{ ++parentsIt; }
        }
      }
      // record this tracks match
      ParticleEffPur match = {recotrkid, recoEff, recoPur, recopdg};
      matches.push_back(match);
    } // end of loop over tracks for matching to truth
    
    return matches;

  }

  //----------------------------------------------------------------------
  std::vector<NeutrinoWithIndex> BackTracker::allNeutrinoInteractions() const
   {
    std::vector <NeutrinoWithIndex> allNeutrinoInteractions;
    
    for (unsigned int s = 0; s < fMCTruthList.size(); s++){
      const art::Ptr<simb::MCTruth> neutrinoInteraction = fMCTruthList[s];
      if(!neutrinoInteraction->NeutrinoSet()) 
        continue;
      
      NeutrinoWithIndex neutrino = {neutrinoInteraction, s};
      allNeutrinoInteractions.push_back(neutrino);
    }

    return allNeutrinoInteractions;
  }

  //----------------------------------------------------------------------
  std::vector<NeutrinoWithIndex> BackTracker::allMCTruth() const
  {
    std::vector <NeutrinoWithIndex> allTruth;
    
    for (unsigned int s = 0; s < fMCTruthList.size(); s++){
      const art::Ptr<simb::MCTruth> neutrinoInteraction = fMCTruthList[s];
      NeutrinoWithIndex neutrino = {neutrinoInteraction, s};
      allTruth.push_back(neutrino);
    }
    
    return allTruth;
  }
  
  
  //----------------------------------------------------------------------
  bool BackTracker::CloseInTime(const rb::CellHit& chit, const sim::PhotonSignal& phot) const
  {
    //The constants used to determine this value change depending on detector
    //and hit timing reconstruction. See docdb 9092 to see how the current 
    //constants were determined (May 2013). To redo the study and see if new
    //constants are required, use MCCheckOut/CloseInTime_module.cc.

    // Empirically determined as an over-generous width of the timing reconstruction.
    double kWindowNs = 1000;
    // Empirically determined as the approximate mean bias of the timing reconstruction.
    // This is the amount the CellHit is expected to be earlier than the photons.
    double kOffsetNs = 350;

    if ((fDetID == novadaq::cnv::kFARDET)&&(chit.GoodTiming())){
      kWindowNs = 500;
      kOffsetNs = 10;
    }

    if ((fDetID == novadaq::cnv::kNEARDET)&&(chit.GoodTiming())){
      kWindowNs = 200;
      kOffsetNs = 10;
    }

    const double dt = chit.TNS()-phot.TimeMean();

    return fabs(dt+kOffsetNs) <= kWindowNs;
  }

  //----------------------------------------------------------------------
  void BackTracker::AccumulateHitContributions(const rb::CellHit& hit,
                                               std::map<int, int> & idToNPhot) const
  {
    auto it = fCellToPhotonSignal.find(geo::OfflineChan(hit.Plane(), hit.Cell()));
    // Nothing to accumulate from this cell
    if(it == fCellToPhotonSignal.end()) return;

    // Candidate photons, in that they're in the right cell.
    const std::vector<sim::PhotonSignal>& phots = it->second;

    // For every photon that's sufficiently in time with the hit,
    // accumulate its energy into the relevant trackID.
    for(unsigned int i = 0; i < phots.size(); ++i){
      const sim::PhotonSignal& phot = phots[i];
      if(CloseInTime(hit, phot)){
        idToNPhot[phot.TrackId()] += phot.NPhoton();
      }
    }
  }
  //----------------------------------------------------------------------

  std::vector<sim::FLSHit> BackTracker::ParticleToFLSHit(const int& trackID) const{
    std::unordered_map< int, std::list<sim::FLSHit> >::const_iterator it = fTrackIDToFLSHit.find(trackID);
    if( it == fTrackIDToFLSHit.end() ){
      std::vector<sim::FLSHit> empty;
      LOG_DEBUG("BackTracker") << "no collection of sim::FLSHit for "
                               <<"Track ID "<<trackID
                               << ". Returning empty vector";

      return empty;
    }
    else
      return std::vector<sim::FLSHit>(it->second.begin(), it->second.end());
  }


  //----------------------------------------------------------------------
  std::vector<sim::FLSHit> BackTracker::ParticleToFLSHit(const int trackID, const int pdg) const{

      const std::vector<sim::FLSHit>& flshits_with_same_track_id = ParticleToFLSHit(trackID);

      std::vector<sim::FLSHit> out_flshits;

      const int nflshits = flshits_with_same_track_id.size();

      /// Fill only FLSHits with a given PDG
      for(int i=0; i<nflshits; ++i){
          const sim::FLSHit& flshit= flshits_with_same_track_id[i];
          if(flshit.GetPDG() == pdg){
            out_flshits.push_back(flshit);
          }
      }

      return out_flshits;
  }

  //----------------------------------------------------------------------
  std::vector<rb::Cluster> BackTracker::ClusterByTruth(std::vector< art::Ptr< rb::CellHit > > const& hits)
  {
    std::vector<rb::Cluster> clusters(fMCTruthList.size()+1);

    for(unsigned int hiti = 0; hiti < hits.size(); ++hiti) {
      art::Ptr<rb::CellHit> hit = hits[hiti];

      // If the hit is noise, put it into the last cluster.
      if(this->IsNoise(hit)) {
        clusters[fMCTruthList.size()].Add(hit);
        continue;
      }

      const sim::Particle *part = this->HitToParticle(hit);
      unsigned int idx = fTrackIDToMCTruthIndex[part->TrackId()];
      clusters[idx].Add(hit);

    } // end loop over hits

    // Label the noise slice.
    clusters[fMCTruthList.size()].SetNoise(true);

    return clusters;
  }

  //----------------------------------------------------------------------
  rb::Cluster BackTracker::MCTruthToCluster(art::Ptr<simb::MCTruth>                const& mct,
                                            std::vector< art::Ptr< rb::CellHit > > const& hits)
  {
    rb::Cluster cluster;

    // figure out where in the MCTruth list mct belongs
    unsigned int MCTidx;
    for(MCTidx = 0; MCTidx < fMCTruthList.size(); ++MCTidx) {
      if(fMCTruthList[MCTidx] == mct) break;
    }

    // loop over hits and add to the cluster
    for(unsigned int hiti = 0; hiti < hits.size(); ++hiti) {
      art::Ptr<rb::CellHit> hit = hits[hiti];

      // If the hit is noise, skip it.
      if(this->IsNoise(hit)) continue;

      const sim::Particle *part = this->HitToParticle(hit);
      unsigned int idx = fTrackIDToMCTruthIndex[part->TrackId()];
      if(idx == MCTidx) cluster.Add(hit);

    } // end loop over hits
    
    return cluster;
  }                                                    

  //----------------------------------------------------------------------

  void BackTracker::HitToParticlesMap(const std::vector<const rb::Cluster*>& sliceList,
                                     const int& sliceIdx)
  {
    fCellToIDs.clear();
    fIdToEnergy.clear();
    
    unsigned int nsli = sliceList.size();

    for (unsigned int isli = 0; isli < nsli; isli++){
      unsigned int nh = sliceList[isli]->NCell();
      for (unsigned int ih = 0; ih < nh; ih++){

        art::Ptr<rb::CellHit> hit = sliceList[isli]->Cell(ih);
        
        if (IsNoise(hit)) continue;
        double chan = (hit->TDC()*1000000) + (hit->Plane()*1000) + hit->Cell();
        
        std::vector<cheat::TrackIDE> ides = HitToTrackIDE(hit);
        fCellToIDs[chan] = ides;

        for(auto ide : ides)
          fIdToEnergy[ide.trackID] += ide.energy;
      }// end of loop over slice hits
    }// end of loop over slices
    fSliceIDForHitParticleMap = sliceIdx;

  }// end of hits to particles table
  
  //----------------------------------------------------------------------

  cheat::ParticleEffPur BackTracker::ClusterToParticle(const std::vector<const rb::Cluster*>& sliceList,
                                                       const std::vector<const rb::CellHit*>& hits,
                                                       const int& sliceIdx)
  {
    
    if(sliceIdx != fSliceIDForHitParticleMap){
      HitToParticlesMap(sliceList, sliceIdx);
    }
    
    double totalE = 0.0;
    unsigned int nh = hits.size();

    // make a map of how much energy each track ID contributes to 
    // this cluster. Also keep track of how much total true energy 
    // there is in this cluster.

    std::map< int, double> idToEnergyInProng;
    for(unsigned int ih = 0; ih < nh; ih++){
      const rb::CellHit* hit = hits[ih];
      double chan = (hit->TDC()*1000000) + (hit->Plane()*1000) + hit->Cell();
      std::vector<cheat::TrackIDE> ides = fCellToIDs[chan];
      for(auto ide: ides){
        idToEnergyInProng[ide.trackID] += ide.energy;
        totalE += ide.energy;
      }// end loop over ides   
    }// end loop over cluster hits

    // now find which trackID contributed the most energy to 
    // the cluster
    int maxID = -1; 
    double maxE = -1.0;

    cheat::ParticleEffPur puff = {-1,-1,-1,-1};
    if(idToEnergyInProng.size() == 0)
      return puff;

    for(auto idtoe: idToEnergyInProng){
      if(idtoe.second > maxE ){
        maxE = idtoe.second;
        maxID = idtoe.first;
      }
    }// end loop over trackID to energy map

    puff.trackID    = maxID;
    puff.purity     = maxE/totalE; 
    puff.efficiency = maxE/fIdToEnergy[maxID];
    puff.pdg        = TrackIDToParticle(maxID)->PdgCode();

    return puff;
  }// end of ClusterToParticle
  
  //----------------------------------------------------------------------

  float BackTracker::CalcEscapingEnergy( const sim::Particle &Par, bool UseMCPart, std::string G4Label ) {
    float TempVar = -5;
    return CalcEscapingEnergy( Par, TempVar, UseMCPart, G4Label );
  }// end of CalcEscapingEnergy
  
  //----------------------------------------------------------------------

  float BackTracker::CalcEscapingEnergy( const sim::Particle &Par, float& EnDep, bool UseMCPart, std::string G4Label ) {
    /// 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.

    // Set to some default values.
    float EnEsc, EnEnt;
    EnEsc = EnEnt = EnDep = -5;

    // --- Do I want to use the information stored in MCParticle?
    if (UseMCPart) {    
      // Make variables to call InterceptsDetector.
      int    TrajEnter, TrajExit;
      TrajEnter = TrajExit = -1;
      // Call InterceptsDetector
      bool IntDet = InterceptsDetector( Par, TrajEnter, TrajExit );
      
      // If IntDet == false then the particle never entered the detector, so return 0's.
      if (!IntDet) {
        return EnEsc;
      } else {
        // Now work out the energy which the particle deposited...
        //  Set EnEnt
        EnEnt = (float)Par.E( TrajEnter ) - Par.Mass();
        //  If TrajExit isn't NTraj-1, the particle left the volume then some energy was lost!
        //  If TrajExit == NTraj-1 then we can just leave EnEsc as 0.
        if ( TrajExit != (int)Par.NumberTrajectoryPoints()-1 ) { 
          EnEsc = (float)Par.E( TrajExit ) - Par.Mass();
        }
        // The energy deposited is just the entering energy - the exiting energy.
        EnDep = EnEnt - EnEsc;
      }
      // now return
      return EnEsc;
    } else {
      std::map<int, int>::iterator TrEIt = fTrackIDToTrueEIndex.find( Par.TrackId() );
      if ( TrEIt != fTrackIDToTrueEIndex.end() ) {
        art::Ptr< sim::TrueEnergy > ThisEn = fTrueEnergyList[ TrEIt->second ];
        /*
          std::cout << "Looking at the sim::TrueEnergy from the data product it is, TrackID " << ThisEn->TrackID() << ", EntEn " << ThisEn->EnteringEnergy()
          << ", EscEn " << ThisEn->EscapingEnergy() << ", TotEDep " << ThisEn->TotalDepEnergy() << ", TotEscEn " << ThisEn->TotalEscEnergy()
          << std::endl;
        */
        EnDep = ThisEn->TotalDepEnergy();
        return ThisEn->TotalEscEnergy();
      } else {
        mf::LogWarning("BackTracker") << "Didn't find TrackId " << Par.TrackId() << " in my map of sim::TrueEnergy's, returning -5.";
        return -5;
      }
    }
  }// end of CalcEscapingEnergy

  //----------------------------------------------------------------------

  float BackTracker::CalcDaughterEscapingEnergy( const sim::Particle &Par ) {
    float DauEsc = 0.;
    int NumDaught = Par.NumberDaughters();
    for (int dau=0; dau<NumDaught; ++dau) {
      const sim::Particle* ThisDau = TrackIDToParticle( Par.Daughter( dau ) );
      float ThisEsc = CalcEscapingEnergy( (*ThisDau), true );
      DauEsc += std::max((float)0., ThisEsc);
      /*
      std::cout << "\tDaughter " << dau << " of " << NumDaught << " was a " << ThisDau->PdgCode() << ", TrackId " << ThisDau->TrackId()
                << ", it had " << ThisDau->NumberDaughters() << " daughters itself. "
                << " EscEn for this particle = " << ThisEsc << " ==> DauEsc = " << DauEsc
                << std::endl;
      */
      if ( ThisDau->NumberDaughters() ) {
        //std::cout << "\t\t****** This particle has daughters, so recall the function... ******" << std::endl;
        float DaughtEsc = CalcDaughterEscapingEnergy( (*ThisDau) );
        DauEsc += DaughtEsc;
        //std::cout << "\t\t****** The sum of all EEscs was " << DaughtEsc << " ==> DauEsc is " << DauEsc << std::endl;
      } // If this particle has daughters itself
    } // Loop through the daughters of my particle.
    //std::cout << "\n\tAt the end of all of that, the sum of all escaping daughters was " << DauEsc << std::endl;
    return DauEsc;
  } // end of CalcDaughterEscapingEnergy

  //----------------------------------------------------------------------

  float BackTracker::CalcTotalEscapingEnergy( const sim::Particle &Par, bool UseMCPart, std::string G4Label ) {
    float TempVar = 0.;
    return CalcTotalEscapingEnergy( Par, TempVar, UseMCPart, G4Label );
  }

  //----------------------------------------------------------------------

  float BackTracker::CalcTotalEscapingEnergy( const sim::Particle &Par, float& TotEnDep, bool UseMCPart, std::string G4Label ) {
    if ( UseMCPart ) {
      float ParEDep = 0;
      float ParEEsc = CalcEscapingEnergy( Par, ParEDep, UseMCPart, G4Label );
      float DauEEsc = CalcDaughterEscapingEnergy( Par );

      float TotEnEsc = std::max((float)0., ParEEsc) + std::max((float)0., DauEEsc);
      TotEnDep = ParEDep - DauEEsc;
      // --- If negative, then had particle had very little energy --> may need to consider rest mass energy.
      if ( TotEnDep < 0 ) {
        //std::cout << "\t\t ** TotEnDep is negative " << TotEnDep << ", initially set to 0, but will check Part type. " << std::endl;
        // --- Default TotEnDep to 0.
        TotEnDep = 0;
        // --- If this particle decayed in the detector (ParEdep != 0) then it's likely that the particles it produced left the detector,
        //     and the sum of their energies will be more than the initial KE of the particle.
        //     For now I will modify this for some choice particles, not sure how best to do this...
        if ( ParEDep != 0 &&
             ( abs(Par.PdgCode()) == 11  || // Electron
               abs(Par.PdgCode()) == 13  || // Muon
               abs(Par.PdgCode()) == 111 || // Pi0
               abs(Par.PdgCode()) == 211 || // Pion
               abs(Par.PdgCode()) == 321 )  // Kaon
             ) {
          /*
          std::cout << "\t\t ** This particle likely decayed, so adding RestMass of " << Par.Mass()
                    << " --- TotEnDep = " << ParEDep << " - " << DauEEsc << " + " << Par.Mass() << " = " << ParEDep - DauEEsc + Par.Mass()
                    << std::endl;
          */
          TotEnDep = std::max( 0., ParEDep - DauEEsc + Par.Mass() );
        }
      }
      /*
      std::cout << "\t***The end of CalcTotalEscapingEnergy -- ParEEsc " << ParEEsc << ", ParEDep " << ParEDep << ", DauEEsc " << DauEEsc << ".\n"
                << "\t TotEnEsc = " << ParEEsc << " + " << DauEEsc << " = " << TotEnEsc << ".\n"
                << "\t TotEnDep = " << ParEDep << " - " << DauEEsc << " = " << TotEnDep << "."
                << std::endl;
      */
      return TotEnEsc;
    } else {
      std::map<int, int>::iterator TrEIt = fTrackIDToTrueEIndex.find( Par.TrackId() );
      if ( TrEIt != fTrackIDToTrueEIndex.end() ) {
        art::Ptr< sim::TrueEnergy > ThisEn = fTrueEnergyList[ TrEIt->second ];
        /*
          std::cout << "Looking at the sim::TrueEnergy from the data product it is, TrackID " << ThisEn->TrackID() << ", EntEn " << ThisEn->EnteringEnergy()
          << ", EscEn " << ThisEn->EscapingEnergy() << ", TotEDep " << ThisEn->TotalDepEnergy() << ", TotEscEn " << ThisEn->TotalEscEnergy()
          << std::endl;
        */
        TotEnDep = ThisEn->TotalDepEnergy();
        return ThisEn->TotalEscEnergy();
      } else {
        mf::LogWarning("BackTracker") << "Didn't find TrackId " << Par.TrackId() << " in my map of sim::TrueEnergy's, returning -5.";
        return -5;
      }
    }
  }
  
  //----------------------------------------------------------------------

  bool BackTracker::InterceptsDetector( const sim::Particle &Par, int &TrajEnt, int &TrajEx ) {
    /// 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.
    
    // Initially set some holder values.
    TrajEnt = TrajEx = -5;
    // Make a handle to geometry service
    art::ServiceHandle<geo::Geometry> geom;
    // Get the Number of trajectory points.
    unsigned int NumTraj = Par.NumberTrajectoryPoints();
    bool EnteredVol = false;
    bool PrevHit    = false;
    for (unsigned int trj=0; trj<NumTraj; ++trj) {
      // Is the hit in the detector?
      bool InDet = geom->isInsideDetectorBigBox( Par.Vx(trj), Par.Vy(trj), Par.Vz(trj) );
      if ( InDet ) {
        // If this is the first hit in the volume.
        if ( !EnteredVol ) {
          TrajEnt    = trj;
          EnteredVol = true;
        }
        // If this is the last trajectory point.
        if ( EnteredVol && InDet && (trj==NumTraj-1) ) {
          TrajEx   = trj;
        }
      }
      // If last hit was in volume, and this hit isn't.
      if ( PrevHit && !InDet) {
        TrajEx = trj-1;
      }
      // Set the PrevHit variable.
      PrevHit = InDet;
    } // Loop over TrajPoints
    // Finally, return
    return EnteredVol;
  } //end InterceptsDetector

  //----------------------------------------------------------------------

} // namespace


namespace cheat{
  DEFINE_ART_SERVICE(BackTracker)
}