Slicer_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////
// \file    Slicer.cxx
// \brief   TODO
// \version $Id: Slicer.cxx,v 1.44 2012-11-09 17:23:01 jpaley Exp $
// \authors rbpatter@caltech.edu, bckhouse@caltech.edu
///////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"

#include "CMap/service/CMapService.h"
#include "Geometry/Geometry.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"

#include <cmath>
#include <string>



// Slicer header
namespace slicer {
  
  class Slicer : public art::EDProducer {
  public:
    explicit Slicer(fhicl::ParameterSet const &pset);
    virtual ~Slicer();
    
    void produce(art::Event& evt);
    
  protected:
    // internal methods
    
    void NeighborSlice(std::vector<art::Ptr<rb::CellHit> > hitlist,
                       std::unique_ptr< std::vector<rb::Cluster> >& clustercol);
    
    void WindowSlice(std::vector<art::Ptr<rb::CellHit> > hitlist,
                     std::unique_ptr< std::vector<rb::Cluster> >& clustercol);
    
    // I might want to use RecoHit::T() or CellHit::TNS().  This method
    // will facilitate that choice.  For now, just CellHits...
    float HitTime(art::Ptr<rb::CellHit> hit) { return hit->TNS(); }
    
    // Helper for NeighborSlice
    bool AreNeighbors(art::Ptr<rb::CellHit> hit_i, art::Ptr<rb::CellHit> hit_j, bool weak = false);
    
    // Helper for WindowSlice
    std::vector<std::vector<art::Ptr<rb::CellHit> > >
    WindowSliceZ(std::vector<art::Ptr<rb::CellHit> >& hits, bool returnAllSlices = false) const;
    
    // Helper for WindowSlice, to make BlockSlices
    std::vector<std::vector<art::Ptr<rb::CellHit> > >
    BlockSlice(std::vector<art::Ptr<rb::CellHit> >& hits) const;
    
    // configuration settings
    std::string  fAlgorithm;          ///< Neighbor/Window
    bool         fSliceByBlock;       ///< If true, slices per block are written
    std::string  fCellHitInput;       ///< Read CellHits from Cal/[fCellHitFolderIn]
    bool         fDoNotSliceAnything; ///< If true, one big slice is written out
    unsigned int fMinHits;            ///< The minimum number of hits needed to create a slice (this is in any view)
    
    // Configuration specific to the Window slicer
    double fTimeWindow;  // ns
    bool   fSeparateZGaps; ///< If true, the neighbor algorithm will separate zGaps
    double fZGapAllowed; // planes
  };
  
}





namespace slicer
{
  //----------------------------------------------------------------------
  Slicer::Slicer(fhicl::ParameterSet const &pset) :
  fAlgorithm         (pset.get< std::string >("Algorithm")         ),
  fSliceByBlock      (pset.get< bool        >("SliceByBlock")      ),
  fCellHitInput      (pset.get< std::string >("CellHitInput")      ),
  fDoNotSliceAnything(pset.get< bool        >("DoNotSliceAnything")),
  fMinHits           (pset.get< int         >("MinHits")           ),
  fTimeWindow        (pset.get< double      >("TimeWindow")        ),
  fSeparateZGaps     (pset.get< bool        >("SeparateZGaps")      ),
  fZGapAllowed       (pset.get< int         >("ZGapAllowed")       )
  {
    assert(fAlgorithm == "Neighbor" || fAlgorithm == "Window");
    
    produces< std::vector<rb::Cluster> >();
  }
  
  //----------------------------------------------------------------------
  Slicer::~Slicer()
  {
  }
  
  //----------------------------------------------------------------------
  void Slicer::produce(art::Event& evt)
  {
    // Load hit list (CellHits)
    art::Handle< std::vector<rb::CellHit> > hitcol;
    evt.getByLabel(fCellHitInput, hitcol);
    
    std::vector<art::Ptr<rb::CellHit > > hitlist;
    for(unsigned int i = 0; i < hitcol->size(); ++i){
      art::Ptr<rb::CellHit> hit(hitcol, i);
      hitlist.push_back(hit);
    }
    
    std::unique_ptr< std::vector<rb::Cluster> > clustercol(new std::vector<rb::Cluster>);
    
    if(fAlgorithm == "Neighbor") NeighborSlice(hitlist, clustercol);
    if(fAlgorithm == "Window")   WindowSlice(hitlist, clustercol);
    
    evt.put(std::move(clustercol));
  }
  
  //----------------------------------------------------------------------
  void Slicer::
  NeighborSlice(std::vector<art::Ptr<rb::CellHit> > hitlist,
                std::unique_ptr< std::vector<rb::Cluster> >& clustercol)
  {
    // Slicer settings (some detector-specific).  Hard-coding these
    // until I have a better feel for the use-cases and the required
    // user interface(s).
    
    int maxNeighbors = 1;
    double bigGap    = 0.; // ns
    double lonerGap  = 0.; // ns
    // So that hits with the same exact time make it into a slice, independent
    // of what order they happen to be sorted in
    double graceGap  = 1; // ns
    int minNumWeakNeighbors = 0; // Don't require anything
    
    //get the geometry
    art::ServiceHandle<geo::Geometry> geo;
    
    // figure out which detector we are using
    int det = geo->DetId();
    
    switch (det) {
      case novadaq::cnv::kNEARDET:
        bigGap = 40;
        lonerGap = 100;
        //bigGap = 530; // With current DCS timing quality, we can only trust gaps longer
        //lonerGap = 530; // than 500 ns (so, using something between 500 and 500+62.5)
        minNumWeakNeighbors = 0;
        break;
      case novadaq::cnv::kNDOS:
        // bigGap = 80;
        // lonerGap = 250;
        bigGap = 530; // With current DCS timing quality, we can only trust gaps longer
        lonerGap = 530; // than 500 ns (so, using something between 500 and 500+62.5)
        // Empirically, have seen hits within one tick like this get missed. More
        // likely in NDOS due to the large number of bad channels. And not too
        // much danger of dragging the whole detector in, because it's smaller.
        graceGap = 62.5+1;
        minNumWeakNeighbors = 0;
        break;
      case novadaq::cnv::kFARDET:
      default: // use kFar setting as default
        bigGap = 100;
        lonerGap = 350;
        //bigGap = 530; // With current DCS timing quality, we can only trust gaps longer
        //lonerGap = 530; // than 500 ns (so, using something between 500 and 500+62.5)
        minNumWeakNeighbors = 1; // Require at least one neighbor
        break;
    }
    
    if (fDoNotSliceAnything) {
      // effectively turn off slicing
      bigGap = 9E10;
    }
    
    
    // Time sort the CellHits in the vector
    rb::SortByTime(hitlist);
    const int nHits = hitlist.size();
    
    // -----------------------------------------------------------------------
    // These lists are in step with the hit list, holding the chunk number
    // each hit is associated with and the number of "neighbors" the hit has
    std::vector<int> hit_chunkID;
    std::vector<int> hit_numNeighbors;
    
    
    // -----------------------------------------------------------------------
    // Count neighbors for each hit.  Hits with zero or one neighbor can't
    // bridge the rough-cut time slices
    
    for (int i=0;i<nHits;i++) {
      hit_numNeighbors.push_back(0);
    }
    for (int i=0;i<nHits;i++) {
      art::Ptr<rb::CellHit> hit_i = hitlist[i];
      for (int j=i+1;j<nHits;j++) {
        art::Ptr<rb::CellHit> hit_j = hitlist[j];
        if (AreNeighbors(hit_i, hit_j)) {
          hit_numNeighbors[i]++;
          hit_numNeighbors[j]++;
        }
      }
    }
    
    
    // -----------------------------------------------------------------------
    // Look for significant gaps in time, skipping any loner hits.  Assign
    // a number (chunkID) to each independent set of hits.
    
    int currChunkID = -1; // will increment
    double prevT = -9e9;
    
    if (fDoNotSliceAnything) currChunkID = 0;
    
    for (int i=0;i<nHits;i++) {
      art::Ptr<rb::CellHit> hit_i = hitlist[i];
      
      if (prevT+bigGap<HitTime(hit_i)) {
        // -- big gap ==> new chunk
        currChunkID++;
        hit_chunkID.push_back(currChunkID);
      }
      else {
        // -- little gap ==> gotta see if it's a loner hit (or loner hits)...
        // skip any loner hits that are present
        int j=i;
        while (j<nHits&&hit_numNeighbors[j]<=maxNeighbors) j++;
        
        // if not a loner, automatic bridging:
        bool bridge_them = (j==i)||fDoNotSliceAnything;
        if (!bridge_them) {
          // does the next non-loner hit bridge the gap anyway?
          if (j<nHits) {
            art::Ptr<rb::CellHit> hit_j = hitlist[j];
            if (!(prevT+bigGap<HitTime(hit_j))) {
              bridge_them = true;
            }
          }
        }
        for (;i<=j&&i<nHits;i++) {
          if (!bridge_them) {
            // new chunk for each loner
            currChunkID++;
          }
          hit_chunkID.push_back(currChunkID);
        }
        hit_i = hitlist[--i];
      }
      
      prevT = HitTime(hit_i);
    }
    const int nChunks = currChunkID+1;
    
    
    // -----------------------------------------------------------------------
    // Now look for chunks with too few hits in them.  Figure out if these
    // hits can be associated with a preceeding (preferred) or a following
    // chunk by seeing if they are neighbors with any of the hits in that
    // chunk AND if they are close enough in time.
    
    // In general, the definition of "neighbor" could be different here.
    // We'll have to tune this up as we learn more about our events.
    
    // Repeat this until nothing new gets saved.  (As hits get saved,
    // other hits may be saved in a second (third?) round.  Shouldn't
    // take more than one or two passes, and speed shouldn't be an issue
    // here, so I'm coding this up as transparently as I can.
    
    bool savedAnyHits = true;
    
    // for counting how many hits each chunk has
    unsigned int numHitsInChunk[nChunks];
    
    while (savedAnyHits) {
      
      savedAnyHits = false;
      
      // count how many hits each chunk has
      for (int iC=0;iC<nChunks;iC++) {
        numHitsInChunk[iC] = 0;
      }
      for (int i=0;i<nHits;i++) {
        numHitsInChunk[hit_chunkID[i]]++;
      }
      
      // go through trying to save the hits in the tiny chunks
      for (int i=0;i<nHits;i++) {
        
        art::Ptr<rb::CellHit> hit_i = hitlist[i];
        
        if (numHitsInChunk[hit_chunkID[i]]<fMinHits) {
          // too few hits in this chunk, so let's try to save the hit(s).
          
          // Have we done something with this hit yet?
          bool saved = false;
          
          // first, try to attach this hit to a preceeding (good) chunk
          // if that doesn't work, try to attach to a follow-on (good) chunk
          for(int step = -1; step <= +1; step += 2){
            // Either from this hit back to the start, or forwards to the end
            for (int j=i+step;j>=0 && j<nHits;j+=step) {
              // j hit isn't in a good chunk, no good to us
              if (numHitsInChunk[hit_chunkID[j]]<fMinHits) continue;
              
              // ok, this j-th hit is in a good chunk
              art::Ptr<rb::CellHit> hit_j = hitlist[j];
              
              const double dt = fabs(HitTime(hit_i)-HitTime(hit_j));
              
              // too far in time already; stop looking in this direction
              if (dt>lonerGap) break;
              
              if (dt<graceGap || AreNeighbors(hit_i,hit_j)) {
                // it's very close in time to the i-th hit or closeish in
                // time and close in space. Save it!
                hit_chunkID[i] = hit_chunkID[j];
                saved = true;
                savedAnyHits = true;
                break;
              }
            } // j loop
            if(saved) break;
          } // step loop
        } // if small chunk
      } // i loop
    } // while (savedAnyHits)
    
    // Note: as long as the loop above terminates when (!savedAnyHits),
    // then numHitsInChunk[] should be up to date (for use below).
    
    //  // RBP - diagnostics:
    //  for (int i=0;i<nHits;i++) {
    //    const rb::CellHit *hit_i = hitlist[i];
    //    cout << i << " " << hit_chunkID[i] << " "
    //   << hit_numNeighbors[i] << " " << hit_i->View() << " "
    //   << hit_i->Plane() << " " << hit_i->Cell() << " " << HitTime(hit_i) << endl;
    //  }
    //  cout << "----------------" << endl;
    
    // Count how many neighbours each hit has from its own chunk, but with
    // weak conditions.
    std::vector<int> numWeakNeighbors;
    numWeakNeighbors.resize(nHits);
    
    for (int i=0;i<nHits-1;i++) {
      currChunkID = hit_chunkID[i];
      for (int j=i+1;j<nHits;j++) {
        if(i == j) continue; // A hit doesn't count as its own neighbour
        if(hit_chunkID[j] != currChunkID) continue; // Must be in same chunk
        
        if(AreNeighbors(hitlist[i], hitlist[j], true)){
          numWeakNeighbors[i]++;
          numWeakNeighbors[j]++;
        }
      }
    }

    // Now we remove hits from chunks if they are really far away from all other hits
    // i.e. no weak neighbors 
    // These hits are relegated to the noise slice
    for (int i=0;i<nHits;i++) 
    {
      if(numWeakNeighbors[i] < minNumWeakNeighbors)
      {
        currChunkID = hit_chunkID[i];
        // Decrement the number of hits in that chunk because we are removing it 
        numHitsInChunk[currChunkID] -= 1;
        // Move that hit to the noise slice
        hit_chunkID[i] = nChunks;
      }
      
    }
    
    // We're going to create vectors of cell hits for our time slices 
    // and try to divide ones with large z-gaps, if desired 
    std::vector<std::vector<art::Ptr<rb::CellHit> > >* tSlices 
    = new std::vector<std::vector<art::Ptr<rb::CellHit> > >(nChunks+1);
    std::vector<std::vector<art::Ptr<rb::CellHit> > >* slices;


    for (int i=0;i<nHits;i++) {
      currChunkID = hit_chunkID[i];
      art::Ptr<rb::CellHit> hit = hitlist[i];
      if (currChunkID>=0&&numHitsInChunk[currChunkID]>=fMinHits) {
        // push to a good time slice
        (*tSlices)[currChunkID].push_back(hit);
      }
      else {
        // stray hit; push to noise cluster
        (*tSlices)[nChunks].push_back(hit);
      }
    }
    

    // Split z-gaps if desired 
    if(fSeparateZGaps){
      slices = new std::vector<std::vector<art::Ptr<rb::CellHit> > >;
      // Loop over time slices to separate z-gaps, excluding the noise slice 
      for(unsigned int ts = 0; ts < tSlices->size() - 1; ++ts)
      {
        int nHitsTime = (*tSlices)[ts].size();
        int nHitsZ = 0;
        // call the z-slicing function with returnAllSlices set to true
        std::vector<std::vector<art::Ptr<rb::CellHit> > > zSlices =
        WindowSliceZ((*tSlices)[ts], true);
        
        for(unsigned int zs = 0; zs < zSlices.size(); ++zs)
        {
          // Order hits by time
          rb::SortByTime(zSlices[zs]);
          // Add to slices vector
          slices->push_back(zSlices[zs]);
          nHitsZ += zSlices[zs].size();
          
        }
        assert(nHitsTime == nHitsZ);
      }
       
      slices->push_back((*tSlices)[nChunks]);
      delete tSlices;
    }else{
      // Slices are just the time slices 
      slices = tSlices;
    
    }
    
    // -----------------------------------------------------------------------
    // For now, call it quits here.  Any chunks with fewer than fMinHits hits
    // will be banished to a "noise" cluster, to be written out last.
    // Similarly for hits really far from all others in their chunk.
    

    // cluster storage:
    std::vector<rb::Cluster> clusters(slices->size());
    
    


    // index of the noise slice
    size_t ns = slices->size() - 1; 
    // Add all of the physics slices, the last in the vector is the noise slice
    for(size_t s = 0; s < ns; ++s)
    {

      if ((*slices)[s].size() >= fMinHits) 
      {
        // add hits from this slice to the cluster
        for(size_t h = 0; h < (*slices)[s].size(); ++h)
        {
          clusters[s].Add((*slices)[s][h]);
        }
      }
      else 
      {
        // slice is too small; move hits to noise slice
        for(size_t h = 0; h < (*slices)[s].size(); ++h)
        {
          
          (*slices)[ns].push_back((*slices)[s][h]);
        }
      }
    }

    // order the noise slice by time
    rb::SortByTime((*slices)[ns]);
    // Add the hits in the noise slice to its cluster 
    for(size_t h = 0; h < (*slices)[ns].size(); ++h)
    {      
      clusters[ns].Add((*slices)[ns][h]);
    }

    // tag the noise cluster as such:
    clusters[ns].SetNoise(true);
    
    if(slices)  delete slices;
    // We're going to do a sanity check to make sure slices contain all the hits
    int nHitsCalHit = hitlist.size();
    int nHitsSlicer = 0;

    
    // Write the clusters.
    // Always write the noise cluster (iC==ns), even if empty.
    //iterate over the vector and remove any that fail to have more than
    //the minimum hits
    for (size_t iC=0;iC < clusters.size(); ++iC) {
      if (iC==ns || clusters[iC].NCell()>=fMinHits) {
        nHitsSlicer += clusters[iC].NCell();
        clustercol->push_back(clusters[iC]);
      }
      else{
        assert(clusters[iC].NCell() == 0 && " skipping a cluster with a non-zero size" );
      }
    }
    assert(nHitsCalHit == nHitsSlicer);
    
    
  }
  
  //----------------------------------------------------------------------
  bool Slicer::AreNeighbors(art::Ptr<rb::CellHit> hit_i, art::Ptr<rb::CellHit> hit_j, bool weak) {
    
    // same view & small plane gap & small cell gap ==> neighbors
    
    const int allowedPlaneGap = weak ? 20 : 5;
    const int allowedCellGap = weak ? 24 : 6;
    
    return
    ( hit_i->View()==hit_j->View() &&
     TMath::Abs(hit_i->Plane()-hit_j->Plane()) <= allowedPlaneGap &&
     TMath::Abs(hit_i->Cell() -hit_j->Cell() ) <= allowedCellGap );
    
  }
  
  //----------------------------------------------------------------------
  void Slicer::
  WindowSlice(std::vector<art::Ptr<rb::CellHit> > hits,
              std::unique_ptr<std::vector<rb::Cluster> >& slices)
  {
    // A timeslice is at most fTimeWindow long, and must contain at least
    // fMinHits hits. Any fTimeWindow-sized subwindow drawn within it
    // will also contain at least fMinHits hits.
    
    slices->clear();
    if(hits.size() == 0) return;
    
    rb::Cluster noiseClust;
    noiseClust.SetNoise(true);
    
    if(hits.empty()){
      slices->push_back(noiseClust);
      return;
    }
    
    std::vector<std::vector<art::Ptr<rb::CellHit> > > tSlices;
    
    rb::SortByTime(hits);
    
    typedef std::vector<art::Ptr<rb::CellHit> >::iterator it_t;
    // The inclusive->exclusive range we're considering for the current slice
    it_t itStart = hits.begin();
    it_t itEnd   = hits.begin();
    
    // We use this point to test if we can expand the slice
    it_t itTest  = hits.begin();
    
    while(true){
      // Grow the slice by moving the end until it's fTimeWindow in size
      while(HitTime(*itEnd) - HitTime(*itStart) < fTimeWindow){
        ++itEnd;
        // We ran out of hits altogether
        if(itEnd == hits.end()) break;
      }
      // Break all the way out of the loop
      if(itEnd == hits.end()) break;
      
      // The window doesn't have enough charge, drop first hit and try again
      if(itEnd - itStart < fMinHits){
        noiseClust.Add(*itStart);
        ++itStart;
        assert(itStart <= itEnd);
        continue;
      }
      
      // The window is long enough and has enough charge
      // Catch itTest up until it defines the start of a fTimeWindow-sized
      // window ending on itEnd
      while(HitTime(*itEnd) - HitTime(*itTest) > fTimeWindow) ++itTest;
      
      // If growing the window more would cause the itTest-itEnd window
      // not to have enough hits, then write out what we have.
      if(itEnd+1 - itTest < fMinHits){
        std::vector<art::Ptr<rb::CellHit> > slice(itStart, itEnd);
        tSlices.push_back(slice);
        itStart = itEnd;
        continue;
      }
      
      // We can still legally grow the window. Do that and try again
      ++itEnd;
      if(itEnd == hits.end()) break;
    }
    
    // If there are enough hits in the final slice, make it one too
    if(itEnd - itStart > fMinHits){
      std::vector<art::Ptr<rb::CellHit> > slice(itStart, itEnd);
      tSlices.push_back(slice);
    }
    else{
      for(it_t it = itStart; it != itEnd; ++it) noiseClust.Add(*it);
    }
    
    // Split every time slice up based on z if it has gaps in
    for(unsigned int s = 0; s < tSlices.size(); ++s){
      std::vector<std::vector<art::Ptr<rb::CellHit> > > zSlices =
      WindowSliceZ(tSlices[s]);
      for(unsigned int n = 0; n < zSlices.size(); ++n){
        rb::Cluster clust;
        for(unsigned int h = 0; h < zSlices[n].size(); ++h){
          clust.Add(zSlices[n][h]);
        }
        slices->push_back(clust);
      }
    }
    
    // Noise cluster goes at the end
    slices->push_back(noiseClust);
  }
  
  ////////////////////////////////////////////////////////////////////////
  std::vector<std::vector<art::Ptr<rb::CellHit> > > Slicer::
  WindowSliceZ(std::vector<art::Ptr<rb::CellHit> >& hits, bool returnAllSlices) const
  {
    std::vector<std::vector<art::Ptr<rb::CellHit> > > zSlices;
    
    // Shouldn't happen, but still...
    if(hits.size() == 0) return zSlices;
    
    rb::SortByPlane(hits);
    
    // The extremes, inclusive and exclusive of our candidate slice
    typedef std::vector<art::Ptr<rb::CellHit> >::iterator it_t;
    it_t itBegin = hits.begin();
    it_t itEnd   = hits.begin();
    
    while(true){
      const int endPlane = (*itEnd)->Plane();
      ++itEnd;
      
      if(itEnd == hits.end()) break;
      
      const int nextPlane = (*itEnd)->Plane();
      
      if(nextPlane - endPlane > fZGapAllowed){
        std::vector<art::Ptr<rb::CellHit> > slice(itBegin, itEnd);
        if(slice.size() > fMinHits || returnAllSlices) zSlices.push_back(slice);
        // Start a new slice
        itBegin = itEnd;
      }
    } // end while
    
    // Whatever hits are left, give them their own slice
    if(itBegin != hits.end()){
      std::vector<art::Ptr<rb::CellHit> > slice(itBegin, itEnd);
      if(slice.size() > fMinHits|| returnAllSlices) zSlices.push_back(slice);
    }
    
    std::vector<std::vector<art::Ptr<rb::CellHit> > > bSlices;
    if(fSliceByBlock) {
      // Split every z slice up based on block position
      for(unsigned int s = 0; s < zSlices.size(); ++s){
        bSlices = BlockSlice(zSlices[s]);
      }
    }
    if(fSliceByBlock) return bSlices;
    else return zSlices;
  }
  
  ////////////////////////////////////////////////////////////////////////
  std::vector<std::vector<art::Ptr<rb::CellHit> > > Slicer::
  BlockSlice(std::vector<art::Ptr<rb::CellHit> >& hits) const
  {
    // Get the geometry
    art::ServiceHandle<geo::Geometry> geo;
    art::ServiceHandle<cmap::CMap>    cmap;
    
    // Figure out which detector we are using
    // We need to know in order to grab the correct block id
    const int det = geo->DetId();
    
    std::vector<std::vector<art::Ptr<rb::CellHit> > > bSlices;
    
    // Shouldn't happen, but still...
    if(hits.size() == 0) return bSlices;
    
    // The extremes, inclusive and exclusive of our candidate slice
    typedef std::vector<art::Ptr<rb::CellHit> >::iterator it_t;
    it_t itBegin = hits.begin();
    it_t itEnd   = hits.begin();
    
    while(true){
      daqchannelmap::lchan lchan = cmap->Map()->encodeLChan(det, (*itEnd)->Plane(),1);
      const int block = cmap->Map()->computeBlock(lchan);
      ++itEnd;
      
      // Check that we have reached end of hit list.
      // If we have, break out of loop
      if(itEnd == hits.end()) break;
      
      daqchannelmap::lchan lchan2 = cmap->Map()->encodeLChan(det, (*itEnd)->Plane(),1);
      const int nextBlock = cmap->Map()->computeBlock(lchan2);
      
      if(block != nextBlock){
        std::vector<art::Ptr<rb::CellHit> > slice(itBegin, itEnd);
        bSlices.push_back(slice);
        // Start a new slice
        itBegin = itEnd;
      }
    } // end while
    
    // Save all of the hits that remain into a block slice
    if(itBegin != hits.end()){
      std::vector<art::Ptr<rb::CellHit> > slice(itBegin,itEnd);
      bSlices.push_back(slice);
    }
    return bSlices;
  }
  
} // end namespace slicer
/////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/ModuleMacros.h"
namespace slicer
{
  DEFINE_ART_MODULE(Slicer)
}