MMSlicer_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////
// \file    MMSlicer.cxx
// \brief   Based on slicer algorithm, but with a ADC cut
// Reproducing the trigger process
// \authors ricken
///////////////////////////////////////////////////////////////////////////
#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 "art/Framework/Core/ModuleMacros.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 MMSlicer : public art::EDProducer {
    public:
      explicit MMSlicer(fhicl::ParameterSet const &pset);
      virtual ~MMSlicer();

      void produce(art::Event& evt);

    protected:
      // internal methods
      void WindowSlice(std::vector<art::Ptr<rb::CellHit> > hitlist,
          std::unique_ptr< std::vector<rb::Cluster> >& clustercol);
      float HitTime(art::Ptr<rb::CellHit> hit) { return hit->TNS(); }

      // Helper for WindowSlice
      std::vector<std::vector<art::Ptr<rb::CellHit> > >
        WindowSliceZ(std::vector<art::Ptr<rb::CellHit> >& hits) const;


      // configuration settings
      std::string  fCellHitInput;       ///< Read CellHits from Cal/[fCellHitFolderIn]
      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
      double fZGapAllowed; // planes
      int    fADCmin,fADCmax;

  };

}

namespace slicer
{
  //----------------------------------------------------------------------
  MMSlicer::MMSlicer(fhicl::ParameterSet const &pset) :
    fCellHitInput      (pset.get< std::string >("CellHitInput")      ),
    fMinHits           (pset.get< int         >("MinHits")           ),
    fTimeWindow        (pset.get< double      >("TimeWindow")        ),
    fZGapAllowed       (pset.get< int         >("ZGapAllowed")       ),
    fADCmin            (pset.get< int         >("ADCmin") ),
    fADCmax            (pset.get< int         >("ADCmax") )
  {
    produces< std::vector<rb::Cluster> >();
  }

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

  //----------------------------------------------------------------------
  void MMSlicer::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);
      if (hit->ADC() > fADCmin && hit->ADC() < fADCmax)
        hitlist.push_back(hit);
    }

    std::unique_ptr< std::vector<rb::Cluster> > clustercol(new std::vector<rb::Cluster>);
    WindowSlice(hitlist, clustercol);
    evt.put(std::move(clustercol));
  }


  //----------------------------------------------------------------------
  void MMSlicer::
    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() < fMinHits ) 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;
          if(itEnd == hits.end()) break;    //out of hits
        }
        if(itEnd == hits.end()) break;      // Break all the way out of the loop

        // The window doesn't have enough hits, drop first hit and try again
        if(itEnd - itStart < fMinHits){      
          ++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);
      }

      //    std::cout<<"the number of time slices: "<<tSlices.size()<<std::endl;

      // 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;
          //Start adding hits:
          for(unsigned int h = 0; h != zSlices[n].size(); ++h){
            clust.Add(zSlices[n][h]);
          }
          if (clust.NXCell()>=3 && clust.NYCell()>=3)
            slices->push_back(clust);
        }
      }
    }

  ////////////////////////////////////////////////////////////////////////
  std::vector<std::vector<art::Ptr<rb::CellHit> > > MMSlicer::
    WindowSliceZ(std::vector<art::Ptr<rb::CellHit> >& hits) const
    {
      std::vector<std::vector<art::Ptr<rb::CellHit> > > zSlices;

      rb::SortByPlane(hits);

      //  std::cout<<"total number hits in the time slice: "<<hits.size()<<std::endl;

      if (hits.size()<fMinHits+1) return zSlices;

      // 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(itEnd != hits.end()-1){
        const int endPlane = (*itEnd)->Plane();
        ++itEnd;
        const int nextPlane = (*itEnd)->Plane();


        if(nextPlane - endPlane > fZGapAllowed){
          std::vector<art::Ptr<rb::CellHit> > slice(itBegin, itEnd);
          if (itEnd-itBegin > fMinHits) zSlices.push_back(slice);

          //        if(slice.size() > fMinHits) zSlices.push_back(slice);
          // Start a new slice
          itBegin = itEnd;
        }
      } // end while


      // Whatever hits are left, give them their own slice
      //    std::cout<<"remaing hits number: "<<itEnd-itBegin<<std::endl;
      std::vector<art::Ptr<rb::CellHit> > slice(itBegin, itEnd);

      if (slice.size() > fMinHits) zSlices.push_back(slice);

      return zSlices;
    }

  ////////////////////////////////////////////////////////////////////////

} // end namespace slicer
/////////////////////////////////////////////////////////////////////////



DEFINE_ART_MODULE(slicer::MMSlicer)