SliceNavigator_service.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file   SliceNavigator_service.cc
// \brief  Look up products dependent on the current slice
// \author Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

#include "EventDisplayServices/SliceNavigator.h"

#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Run.h"

#include "EventDisplayServices/GeometryDrawingOptions.h"
#include "nuevdb/EventDisplayBase/NavState.h"
#include "nuevdb/EventDisplayBase/ServiceTable.h"
#include "nuevdb/EventDisplayBase/Printable.h"
#include "nuevdb/EventDisplayBase/EventHolder.h"

#include "Utilities/AssociationUtil.h"

#include "MEFinder/MEClusters.h"
#include "GeometryObjects/OfflineChan.h"
#include "GeometryObjects/CellGeo.h"
#include "Geometry/Geometry.h"
#include "MCCheater/BackTracker.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/HoughResult.h"
#include "RecoBase/Prong.h"
#include "RecoBase/Track.h"
#include "RecoBase/Vertex.h"
#include "RecoBase/FilterList.h"
#include "nusimdata/SimulationBase/MCTruth.h"

#include "TTimer.h"

namespace evd
{
  //---------------------------------------------------------------
  SliceNavigator::SliceNavigator(const fhicl::ParameterSet& pset,
                                 art::ActivityRegistry& reg)
    : evdb::Reconfigurable{pset},
      fAutoPrintCount(0),
      fRun(0),
      fFullEvent(true),
      fCurrentSlice(0),
      fMaxSliceId(0),
      fInitialSlice(-1),
      fTimer(0)
  {
    reconfigure(pset);
    
    evdb::ServiceTable::OverrideCategory("SliceNavigator",
                                         evdb::kDRAWING_SERVICE);

    reg.sPostProcessEvent.watch(this, &SliceNavigator::postProcessEvent);
    reg.sPreBeginRun.watch(this, &SliceNavigator::preBeginRun);
  }

  //---------------------------------------------------------------
  SliceNavigator::~SliceNavigator()
  {
    fSlicesToNavigate.clear();
    delete fTimer;
  }

  //---------------------------------------------------------------
  void SliceNavigator::reconfigure(fhicl::ParameterSet const& pset)
  {
    fSliceLabel          = pset.get<std::string>("SliceLabel");
    fCellHitLabel        = pset.get<std::string>("CellHitLabel");
    fAutoAdvanceInterval = pset.get<int>("AutoAdvanceInterval");
    fAutoPrintPattern    = pset.get<std::string>("AutoPrintPattern");
    fAutoPrintMax        = pset.get<int>("AutoPrintMax");
    fFilterResults       = (pset.get<int>("SliceOpt.val") & 1);
    fDimHits             = (pset.get<int>("SliceOpt.val") & 2);
    fDimOthers           = (pset.get<int>("SliceOpt.val") & 4);

    fFullEvent = !pset.get<bool>("InitiallyEnabled");
    fInitialSlice = pset.get<int>("InitialSlice");

    this->ParseModuleInstanceString(fSliceLabel, fSliceInstance);
    this->ParseModuleInstanceString(fCellHitLabel, fCellHitInstance);
    
  }
  
  //---------------------------------------------------------------
  void SliceNavigator::ParseModuleInstanceString(std::string & mod,
                                                 std::string & inst)
  {
    const size_t splitPos = mod.find_first_of('_');
    if(splitPos == std::string::npos){
      inst = "";
      return;
    }
    inst = mod.substr(splitPos+1);
    mod.resize(splitPos);
  }

  //---------------------------------------------------------------
  void SliceNavigator::preBeginRun(art::Run const & r)
  {
    fRun = r.run();
  }

  //---------------------------------------------------------------
  void SliceNavigator::postProcessEvent(const art::Event& evt)
  {
    // Clear the slices to navigate vector
    fSlicesToNavigate.clear();

    art::Handle<std::vector<rb::Cluster>> slices;
    
    art::InputTag sliceTag(fSliceLabel, fSliceInstance);
    
    evt.getByLabel( sliceTag, slices);

    if(slices.failedToGet()){
      std::cerr << "Couldn't find any slices under label "
                << fSliceLabel << std::endl;
      fMaxSliceId = 0;
      return;
    }

    int nSli = slices->size();
    for( int iSli = 0; iSli < nSli; ++iSli){
      // if we only want to look at slices that
      // survived filtering, throw the rest away
      if(!fFullEvent && fFilterResults && rb::IsFiltered(evt,slices,iSli) )
      continue;
      if(!fFullEvent && (*slices)[iSli].IsNoise())
      continue;

      fSlicesToNavigate.push_back( iSli );
      if(iSli == fInitialSlice){
        fCurrentSlice = fSlicesToNavigate.size()-1;
        fInitialSlice = -1; // Don't do this again
      }
    }
    if (fSlicesToNavigate.size() > 0) fMaxSliceId = fSlicesToNavigate.back();

    // Protecting against segfaults:  If we've advanced to the next event
    // with fCurrentSlice set to the max number of slices, and the next event
    // has fewer slices than the previous, we need to make sure we don't 
    // suffer from an index out of range problem.
    if( (fCurrentSlice + 1 > int(fSlicesToNavigate.size())) && !fSlicesToNavigate.empty() ) 
      fCurrentSlice = fSlicesToNavigate.size() - 1;

    if(fCurrentSlice == -1){
      if(!fSlicesToNavigate.empty())
        // We pressed the prev slice button
        fCurrentSlice = fSlicesToNavigate.size()-1;
      else
        evdb::NavState::Set(evdb::kPREV_EVENT);
    }

    if(fCurrentSlice == -2 ){
      if(!fSlicesToNavigate.empty())
        // We pressed the next button, but there
        // no more slices left in the previous event
        fCurrentSlice = 0;
      else
        evdb::NavState::Set(evdb::kNEXT_EVENT);
    }

    art::Handle<std::vector<rb::Cluster>> slicecol;
    evt.getByLabel(sliceTag, slicecol);
    fSlices.clear();
    for(unsigned int i = 0; i < slicecol->size(); ++i)
      fSlices.push_back(art::Ptr<rb::Cluster>(slicecol, i));

    art::Handle<std::vector<rb::CellHit>> hitcol;
    art::InputTag cellTag(fCellHitLabel, fCellHitInstance);
    evt.getByLabel(cellTag, hitcol);
    fAllHits.clear();
    for(unsigned int i = 0; i < hitcol->size(); ++i)
      fAllHits.push_back(art::Ptr<rb::CellHit>(hitcol, i));
      
    // Auto-printing
    if (fAutoPrintMax > 0 || fAutoPrintMax == -1){
      ++fAutoPrintCount;
      std::map<std::string, evdb::Printable*>& ps = evdb::Printable::GetPrintables();
      for(std::map<std::string, evdb::Printable*>::iterator it = ps.begin(); it != ps.end(); ++it){
        evdb::Printable* p = it->second;
        
        std::string filename = fAutoPrintPattern;
        // std::string filename = "evd-slice_r%r_e%e_s%s.gif";
        
        // Insert instances of run variable into filename
        if (filename.find("%r") != std::string::npos)
          stringReplace(filename, "%r", std::to_string(fRun));
        
        // Insert instances of event variable into filename
        if (filename.find("%e") != std::string::npos)
          stringReplace(filename, "%e", std::to_string(evt.event()));
        
        // Insert instances of slice variable into filename
        if (filename.find("%s") != std::string::npos)
          stringReplace(filename, "%s", std::to_string(this->CurrentSlice()));
        
        // Print file
        p->Print(TString(filename));
      }
      if (fAutoPrintCount >= (fAutoPrintMax - 1) && fAutoPrintMax != -1)
      {
        std::cout << "Hit autoprint max and exiting." << std::endl;
        exit(0);
      }
    }

    if(fTimer){
      // Auto-advance is active. Fire off one more shot
      fTimer->Start(fAutoAdvanceInterval, true);
    }
  }

  //---------------------------------------------------------------
  template<class T> size_t SliceNavigator::
  GetProducts(const art::Event&      evt,
              const std::string&     which,
              const std::string&     instance,
              std::vector<const T*>& prods)
  {
    return GetProductsImpl(evt, which, instance, prods);
  }

  //---------------------------------------------------------------
  template<class T> size_t SliceNavigator::
  GetProductsImpl(const art::Event&      evt,
                  const std::string&     which,
                  const std::string&     instance,
                  std::vector<const T*>& prods)
  {
    if(fFullEvent)
      return GetProductsImplFull(evt, which, instance, prods);
    else
      return GetProductsImplPart(evt, which, instance, prods);
  }


  //---------------------------------------------------------------
  // Need a special function for clusts again
  // because cluster-cluster assns are not allowed.
  // See generic GetProductsImplFull for the story
  template<> size_t SliceNavigator::
  GetProductsImplFull(const art::Event&      evt,
                      const std::string&     which,
                      const std::string&     instance,
                      std::vector<const rb::Cluster*>& clusts)
  {
    std::vector<const rb::Cluster*> temp(clusts);
    try{
      evt.getView(which, instance, temp);
      temp.swap(clusts);
    }
    catch(cet::exception& e){
      std::cerr << "Failed to get products of type " << typeid(rb::Cluster).name()
                << " under label \"" << which << "\": "
                << e << std::endl;
    }
    return clusts.size();
  }

  //---------------------------------------------------------------
  template<class T> size_t SliceNavigator::
  GetProductsImplFull(const art::Event&      evt,
                      const std::string&     which,
                      const std::string&     instance,
                      std::vector<const T*>& prods)
  {

    art::Handle<std::vector<T> > tempHandle;

    std::vector<art::Ptr<T> > temp;

    try{
      evt.getByLabel(art::InputTag(which, instance), tempHandle);
      art::fill_ptr_vector( temp,tempHandle);
    }
    catch(cet::exception& e){
      std::cerr << "Failed to get products of type " << typeid(T).name()
                << " under label \"" << which << "\": "
                << e << std::endl;
    }

    std::vector<art::Ptr<T> > tempcopy;
    // If fFilterResults is turned on, we don't want to
    // draw anything from a slice that was filtered out
    // So each product has to be mapped to its associated
    // slice cluster and check if the slice was filtered out.
    if( fFilterResults &&
        typeid(T) != typeid(rb::CellHit) &&
        typeid(T) != typeid(simb::MCTruth)){

      art::Handle<std::vector<rb::Cluster> > slicesHandle;
      evt.getByLabel( art::InputTag(fSliceLabel, fSliceInstance), slicesHandle);
      std::vector<art::Ptr<rb::Cluster> >    slices;
      art::fill_ptr_vector( slices, slicesHandle);
      art::FindOneP<rb::Cluster> fot(temp, evt, which);

      if(!fot.isValid()){
        std::cerr << "Nothing of type " << typeid(T).name()
        << " with label \"" << which << "\" and instance \"" << instance
        << "\" associated with slices" << std::endl;
        return 0;
      }

      int nT = temp.size();
      for( int i = 0; i < nT; i++){

        art::Ptr<rb::Cluster> slice = fot.at( i );
        std::vector< art::Ptr< rb::Cluster> >::iterator slicePos =
        std::find(slices.begin(), slices.end(), slice) ;
        int iSlice = std::distance( slices.begin(), slicePos );
        
        if( fFilterResults && rb::IsFiltered(evt,slicesHandle, iSlice ) )
        continue;
        tempcopy.push_back( temp[i] );
        
      }// end loop over temp
    }
    else
      tempcopy.swap( temp );

    for( art::Ptr<T> p: tempcopy)
      prods.push_back( p.get());

    return prods.size();
  }


  //---------------------------------------------------------------
  template<class T> size_t SliceNavigator::
  GetProductsImplPart(const art::Event&      evt,
                      const std::string&     which,
                      const std::string&     instance,
                      std::vector<const T*>& prods)
  {
    art::Handle<std::vector<rb::Cluster>> slices;
    evt.getByLabel(art::InputTag(fSliceLabel, fSliceInstance), slices);

    if(slices.failedToGet()){
      std::cerr << "Couldn't find any slices under label "
                << fSliceLabel << std::endl;
      return 0;
    }

    art::FindManyP<T> fmt(slices, evt, art::InputTag(which, instance));
    if(!fmt.isValid()){
      std::cerr << "Nothing of type " << typeid(T).name()
                << " with label \"" << which << "\" and instance \"" << instance
                << "\" associated with slices" << std::endl;
      return 0;
    }

    std::vector<art::Ptr<T>> prodPtrs = fmt.at( fSlicesToNavigate[fCurrentSlice] );
    prods.clear();
    for(art::Ptr<T> p: prodPtrs) prods.push_back(p.get());
    return prods.size();
  }

  //---------------------------------------------------------------
  // Special case for Tracks because there's this MichelE trick
  template<> size_t SliceNavigator::
  GetProducts(const art::Event&              evt,
              const std::string&             whichTrack,
              const std::string&             whichMichel,
              const std::string&             instance,
              std::vector<const rb::Track*>& tracks)
  {
    std::vector<const me::TrkME*> clusts;
    GetProductsImpl(evt, whichTrack, instance, tracks);
    GetProductsImpl(evt, whichMichel, instance, clusts);

    for(const me::TrkME* c: clusts)
      tracks.push_back(c->ParentTrk().get());

    return tracks.size();
  }

  //---------------------------------------------------------------
  // Special case clusters because they *are* the slice, not just associated
  template<> size_t SliceNavigator::
  GetProducts(const art::Event&                evt,
              const std::string&               which,
              const std::string&               instance,
              std::vector<const rb::Cluster*>& clusts)
  {
    // If it's the actual slices, look them up by index.
    // Otherwise we'll get nothing
    // (not possible to associate Cluster with Cluster)
    if(which == fSliceLabel && instance == ""){
      art::Handle<std::vector<rb::Cluster>> slices;
      evt.getByLabel(fSliceLabel, slices);
      if(slices.failedToGet()){
        std::cerr << "Couldn't find any slices under label "
                  << fSliceLabel << std::endl;
        return 0;
      }
      if( fFullEvent ){
        int nSli = slices->size();
        for( int iSli=0; iSli < nSli; ++iSli){
          if ( fFilterResults && rb::IsFiltered(evt,slices, iSli ) )
          continue;
          clusts.push_back((art::Ptr<rb::Cluster>(slices, iSli).get())) ;
        }
      }// if full event
      else{
        int idx = fSlicesToNavigate[ fCurrentSlice];
        if( fCurrentSlice >=0 )
        clusts.push_back( (art::Ptr<rb::Cluster>(slices, idx)).get() );
      }
    }
    else
      GetProductsImplFull(evt, which, instance, clusts);

    return clusts.size();
  }

  //---------------------------------------------------------------
  // Special case for MCTruth because we have to go via backtracker
  template<> size_t SliceNavigator::
  GetProducts(const art::Event&              evt,
              const std::string&             which,
              const std::string&             instance,
              std::vector<const simb::MCTruth*>& truths)
  {
    truths.clear();

    if(fFullEvent){
      GetProductsImpl(evt, which, instance, truths);
    }
    else{
      art::ServiceHandle<cheat::BackTracker> bt;

      std::vector<cheat::NeutrinoEffPur> eps =
        bt->SliceToMCTruth(fSlices[fSlicesToNavigate[fCurrentSlice]],
                           fAllHits);

      if(!eps.empty()) truths = {eps.front().neutrinoInt.get()};
    }

    return truths.size();
  }

  //---------------------------------------------------------------
  bool SliceNavigator::FullEventMode() const
  {
    return fFullEvent;
  }

  //---------------------------------------------------------------
  void SliceNavigator::SetFullEventMode(bool mode)
  {
    if(mode == fFullEvent) return;

    fFullEvent = mode;
    evdb::NavState::Set(evdb::kRELOAD_EVENT);
  }

  //---------------------------------------------------------------
  void SliceNavigator::Previous()
  {
    --fCurrentSlice;

    if(fCurrentSlice < 0){
      fCurrentSlice = -1;
      evdb::NavState::Set(evdb::kPREV_EVENT);
      fSlicesToNavigate.clear();
    }
    else{
      // Go to the slice requested
      evdb::NavState::Set(evdb::kRELOAD_EVENT);
    }
  }

  //---------------------------------------------------------------
  void SliceNavigator::Next()
  {
    ++fCurrentSlice;

    if(fCurrentSlice >= int(fSlicesToNavigate.size())){
      fCurrentSlice = -2;
      fSlicesToNavigate.clear();
      evdb::NavState::Set(evdb::kNEXT_EVENT);
    }
    else{
      // Go to the slice requested
      evdb::NavState::Set(evdb::kRELOAD_EVENT);
    }
  }

  //---------------------------------------------------------------
  void SliceNavigator::AutoAdvance()
  {
    delete fTimer;
    fTimer = new TTimer;
    fTimer->SetObject(this);
    // One-shot timer, will be re-primed by postProcessEvent()
    fTimer->Start(fAutoAdvanceInterval, true);
  }

  //---------------------------------------------------------------
  void SliceNavigator::StopAutoAdvance()
  {
    if(!fTimer) return;
    fTimer->Stop();
    delete fTimer;
    fTimer = 0;
  }

  //---------------------------------------------------------------
  bool SliceNavigator::AutoAdvancing()
  {
    return fTimer != 0;
  }

  //---------------------------------------------------------------
  int SliceNavigator::CurrentSlice() const
  {
    if(fSlicesToNavigate.empty() || fCurrentSlice < 0) return 0;

    return fSlicesToNavigate [fCurrentSlice ];
  }

  //---------------------------------------------------------------
  int SliceNavigator::CurrentSliceIndex() const
  {
    return fCurrentSlice;
  }

  //---------------------------------------------------------------
  void SliceNavigator::SetCurrentSlice(int slice)
  {
    auto it = std::find(fSlicesToNavigate.begin(),
                        fSlicesToNavigate.end(),
                        slice);
    if(it != fSlicesToNavigate.end()){
      fCurrentSlice = it-fSlicesToNavigate.begin();
      evdb::NavState::Set(evdb::kRELOAD_EVENT);
    }
    else{
      // What else can we do?
      fCurrentSlice = 0;
      evdb::NavState::Set(evdb::kRELOAD_EVENT);
    }
  }

  //---------------------------------------------------------------
  int SliceNavigator::MaxSlice() const
  {
    return fMaxSliceId;
  }

  //---------------------------------------------------------------
  void SliceNavigator::GetBounds(const art::Event& evt,
                                 double& xmin, double& xmax,
                                 double& ymin, double& ymax,
                                 double& zmin, double& zmax,
                                 double trimfrac)
  {
    // Fallback to the whole detector in case we don't get values below
    art::ServiceHandle<geo::Geometry> geom;
    xmin = -geom->DetHalfWidth();
    xmax = +geom->DetHalfWidth();
    ymin = -geom->DetHalfHeight();
    ymax = +geom->DetHalfHeight();
    zmin = 0;
    zmax = geom->DetLength();

    art::Handle<std::vector<rb::Cluster>> slices;
    evt.getByLabel( art::InputTag(fSliceLabel, fSliceInstance), slices);

    if(slices.failedToGet()){
      std::cerr << "Couldn't find any slices under label "
                << fSliceLabel << std::endl;
      return;
    }

    if(fCurrentSlice >= 0 && fCurrentSlice < int(slices->size())){
      const rb::Cluster& slice = (*slices)[fSlicesToNavigate[fCurrentSlice]];

      if(trimfrac == 0){
        const TVector3 r0 = slice.MinXYZ();
        const TVector3 r1 = slice.MaxXYZ();
        xmin = r0.X(); xmax = r1.X();
        ymin = r0.Y(); ymax = r1.Y();
        zmin = r0.Z(); zmax = r1.Z();
        return;
      }

      std::vector<double> xs, ys, zs;

      for(unsigned int i = 0; i < slice.NCell(); ++i){
        const art::Ptr<rb::CellHit>& chit = slice.Cell(i);
        double xyz[3];
        geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
        zs.push_back(xyz[2]);
        if(chit->View() == geo::kX)
          xs.push_back(xyz[0]);
        else
          ys.push_back(xyz[1]);
      }

      std::sort(xs.begin(), xs.end());
      std::sort(ys.begin(), ys.end());
      std::sort(zs.begin(), zs.end());

      xmin = xs[xs.size()*trimfrac];
      ymin = ys[ys.size()*trimfrac];
      zmin = zs[zs.size()*trimfrac];

      // For small events we want to pick up the last hit, not the N-1th hit
      xmax = xs[(xs.size()-1)*(1-trimfrac)+1];
      ymax = ys[(ys.size()-1)*(1-trimfrac)+1];
      zmax = zs[(zs.size()-1)*(1-trimfrac)+1];
    }
  }

  //---------------------------------------------------------------
  void SliceNavigator::GetBox(const art::Event& evt,
                              double& xmin, double& xmax,
                              double& ymin, double& ymax,
                              double& zmin, double& zmax,
                              double trimfrac)
  {
    // Fallback to the whole detector in case we don't get values below
    art::ServiceHandle<geo::Geometry> geom;
    xmin = -geom->DetHalfWidth();
    xmax = +geom->DetHalfWidth();
    ymin = -geom->DetHalfHeight();
    ymax = +geom->DetHalfHeight();
    zmin = 0;
    zmax = geom->DetLength();

    // Use Fiducial dimensions if smaller than slice
    art::ServiceHandle<evd::GeometryDrawingOptions> drawopt;

    float xbox=0;
    float ybox=0;
    float zbox=0;
    if ( (drawopt->fSliceBoxSize).size() == 3 ) {
      std::vector< float > sizes = drawopt->fSliceBoxSize;
      xbox = sizes[0];
      ybox = sizes[1];
      zbox = sizes[2];
    }
    art::Handle<std::vector<rb::Cluster>> slices;
    evt.getByLabel( art::InputTag(fSliceLabel, fSliceInstance), slices);

    if(slices.failedToGet()){
      std::cerr << "Couldn't find any slices under label "
                << fSliceLabel << std::endl;
      return;
    }

    if(fCurrentSlice >= 0 && fCurrentSlice < int(slices->size())){
      const rb::Cluster& slice = (*slices)[fSlicesToNavigate[fCurrentSlice]];

      if(trimfrac == 0){
        const TVector3 r0 = slice.MinXYZ();
        const TVector3 r1 = slice.MaxXYZ();
        xmin = r0.X(); xmax = r1.X();
        ymin = r0.Y(); ymax = r1.Y();
        zmin = r0.Z(); zmax = r1.Z();
        return;
      }

      std::vector<double> xs, ys, zs;

      for(unsigned int i = 0; i < slice.NCell(); ++i){
        const art::Ptr<rb::CellHit>& chit = slice.Cell(i);
        double xyz[3];
        geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
        zs.push_back(xyz[2]);
        if(chit->View() == geo::kX)
          xs.push_back(xyz[0]);
        else
          ys.push_back(xyz[1]);
      }
      std::sort(xs.begin(), xs.end());
      std::sort(ys.begin(), ys.end());
      std::sort(zs.begin(), zs.end());

      xmin = xs[xs.size()*trimfrac];
      ymin = ys[ys.size()*trimfrac];
      zmin = zs[zs.size()*trimfrac];

      // For small events we want to pick up the last hit, not the N-1th hit
      xmax = xs[(xs.size()-1)*(1-trimfrac)+1];
      ymax = ys[(ys.size()-1)*(1-trimfrac)+1];
      zmax = zs[(zs.size()-1)*(1-trimfrac)+1];

      double xmid=(xmax+xmin)/2;
      double ymid=(ymax+ymin)/2;
      double zmid=(zmax+zmin)/2;
      xmax= ( xbox==0 ) ? xmax : xmid+xbox/2;
      ymax= ( ybox==0 ) ? ymax : ymid+ybox/2;
      zmax= ( zbox==0 ) ? zmax : zmid+zbox/2;
      xmin= ( xbox==0 ) ? xmin : xmid-xbox/2;
      ymin= ( ybox==0 ) ? ymin : ymid-ybox/2;
      zmin= ( zbox==0 ) ? zmin : zmid-zbox/2;
    }
  }

  //---------------------------------------------------------------
  std::vector<const rb::Cluster*> SliceNavigator::SlicesToNavigate(const art::Event& evt)
  {
    std::vector<const rb::Cluster*> clusts;
    art::Handle<std::vector<rb::Cluster>> slices;
    evt.getByLabel(art::InputTag(fSliceLabel, fSliceInstance), slices);
    if(slices.failedToGet()) return {};

    for(unsigned int i = 0; i < fSlicesToNavigate.size(); ++i)
      clusts.push_back(&(*slices)[fSlicesToNavigate[i]]);

    return clusts;
  }// End of SlicesToNavigate

  //---------------------------------------------------------------
  std::string SliceNavigator::stringReplace(std::string &s,
                               const std::string &toReplace,
                               const std::string &replaceWith)
  {
      return(s.replace(s.find(toReplace), toReplace.length(), replaceWith));
  }

  //---------------------------------------------------------------
  Bool_t SliceNavigator::HandleTimer(TTimer* t)
  {
    // This is the auto-advance timer
    Next();
    return true;
  }

  // Need to explicitly instantiate every product that can be got. It's a lot
  // of typing, so use a macro
#define INSTANTIATE(T) template size_t SliceNavigator::GetProducts<T>(const art::Event&, const std::string&, const std::string&, std::vector<const T*>&)

  INSTANTIATE(geo::OfflineChan);
  INSTANTIATE(rb::CellHit);
  INSTANTIATE(rb::HoughResult);
  INSTANTIATE(rb::Prong);
  INSTANTIATE(rb::Vertex);
  INSTANTIATE(rb::Track);
  INSTANTIATE(simb::MCTruth);

  DEFINE_ART_SERVICE(SliceNavigator)

} // end namespace evd
////////////////////////////////////////////////////////////////////////