Recluster_module.cc

Go to the documentation of this file.

#include <fstream>
#include <cstdlib>
#include <vector>
#include <string>
#include "TVector3.h"

// ART includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Utilities/InputTag.h"

// NOvA includes
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Shower.h"
#include "RecoBase/FilterList.h"
#include "Utilities/AssociationUtil.h"
#include "ShowerLID/DeconvolveAlg.h"
#include "RecoBase/Vertex.h"

namespace slid{

  class Recluster : public art::EDProducer
  {

  public:

    Recluster( fhicl::ParameterSet const& pset );
    ~Recluster();

    void produce(art::Event& evt);
    void reconfigure(const fhicl::ParameterSet& pset);
    void beginRun(art::Run& run);

  private:

    std::string fSliceLabel;         ///< Label of slicer module
    bool        fProngInput;         ///< Lable of prong producing module
    std::string fProngLabel;         ///< True = Prong input, False = Track input 
    std::string fVertexLabel;        ///< Lable of vertex producing module
    std::string fProng3DInstLabel;     ///< Label of instance label of 3d prongs
    std::string fProng2DInstLabel;     ///< Label of instance label of 2d prongs
    std::string fTrack3DLabel;     ///< Label of 3d tracks
    std::string fTrack2DLabel;     ///< Label of 2d tracks
    std::vector<std::string> fFilterLabels; ///< Labels of filter lists to obey

    bool        fRespectFilter;      ///< Skip slices that failed filter?
    bool        fRecluster;          ///< Should shower reclustering be performed?
    bool        fDeconvolve;         ///< Should cell-energy deconvolution be performed?
    bool        fConsider2D;         ///< Should 2-D prongs be considered for reclustering/deconvolution?
    bool        fApplyPEThreshold;   ///< If this parameter is true, cellhits with PE below fPEThreshold will not be reclustered
    float       fMaxPlane;           ///< The extent of the reclustering cylinder, beyond prong stop point
    float       fMIPPlane;           ///< The extent of the MIP cylinder from the prong start point
    float       fMaxRadius;          ///< The radius of the reclustering cylinder, in cell widths
    float       fMIPRadius;          ///< The radius of the MIP cylinder, in cell widths

    float       fPEThreshold;        ///< If a cell PE is less than this value, it will not be added to the shower.


    DeconvolveAlg*       fDeconvolveAlg;     ///< An object of DeconvolveAlg
    fhicl::ParameterSet  fDeconvolveAlgPSet; ///< ParameterSet of DecolvolveAlg
    art::ServiceHandle<geo::Geometry> fGeom; ///< A private handle to geometry service

  };// end Recluster class definition


  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////


  Recluster::Recluster(fhicl::ParameterSet const& pset)
    : fDeconvolveAlg(0)
    , fDeconvolveAlgPSet(pset.get< fhicl::ParameterSet >("DeconvolveAlgPSet") )
  {
    this->reconfigure(pset);

    // Produces a vector of reclustered showers and their
    // association with slices
    produces< std::vector <rb::Shower> >();
    produces< art::Assns  <rb::Shower, rb::Cluster> >();
    produces< art::Assns  <rb::Shower, rb::Prong> >();

  }// End of constructor


  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  Recluster::~Recluster()
  {
    if( fDeconvolveAlg )
      delete fDeconvolveAlg;
  }

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  void Recluster::reconfigure( const fhicl::ParameterSet& pset )
  {
    fSliceLabel       = pset.get< std::string >("SliceLabel");
    fProngInput       = pset.get< bool >("ProngInput");
    fProngLabel       = pset.get< std::string >("ProngLabel");
    fVertexLabel      = pset.get< std::string >("VertexLabel");
    fProng3DInstLabel   = pset.get< std::string >("Prong3DInstLabel");
    fProng2DInstLabel   = pset.get< std::string >("Prong2DInstLabel");
    fTrack3DLabel   = pset.get< std::string >("Track3DLabel");
    fTrack2DLabel   = pset.get< std::string >("Track2DLabel");


    fMaxPlane         = pset.get< float >("MaxPlane");
    fMIPPlane         = pset.get< float >("MIPPlane");
    fMaxRadius        = pset.get< float >("MaxRadius");
    fMIPRadius        = pset.get< float >("MIPRadius");
    fPEThreshold      = pset.get< float >("PEThreshold");

    fRespectFilter    = pset.get< bool  >("RespectFilter");
    fRecluster        = pset.get< bool  >("Recluster");
    fDeconvolve       = pset.get< bool  >("Deconvolve");
    fConsider2D       = pset.get< bool  >("Consider2D");
    fApplyPEThreshold = pset.get< bool  >("ApplyPEThreshold");
    fFilterLabels     = pset.get< std::vector<std::string> >("FilterLabels");
  }

  /*********************************************************************/
   //////////////////////////////////////////////////////////////////////

  void Recluster::beginRun(art::Run& run)
  {
    fDeconvolveAlg = new DeconvolveAlg(fDeconvolveAlgPSet);
  }

  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  void Recluster::produce( art::Event& evt)
  {

    std::unique_ptr< std::vector< rb::Shower > >
      reclustCol( new std::vector< rb::Shower > );
    std::unique_ptr< art::Assns< rb::Shower, rb::Cluster > >
      sliceAssassin( new art::Assns< rb::Shower, rb::Cluster > );
    // TODO: This is strictly temporary. Once validation
    // is over, we won't need prong-shower association any
    // longer and shouldn't be produced.
    std::unique_ptr< art::Assns< rb::Shower, rb::Prong > >
      prongAssassin( new art::Assns< rb::Shower, rb::Prong > );

    // grab the slices in this event
    art::Handle< std::vector <rb::Cluster> > slices;
    evt.getByLabel( fSliceLabel, slices );

    std::vector< rb::Shower > showerVec;
    // grab the prong-slice associations in this event
    art::FindManyP<rb::Prong> 
      fmprong3D(slices, evt, art::InputTag(fProngLabel, fProng3DInstLabel));
    art::FindManyP<rb::Prong> 
      fmprong2D(slices, evt, art::InputTag(fProngLabel, fProng2DInstLabel));

    art::FindManyP<rb::Track>
      fmtrack3D(slices, evt, fTrack3DLabel);
    art::FindManyP<rb::Track>
      fmtrack2D(slices, evt, fTrack2DLabel);



    // grab the vertex-slice associations in this event
    art::FindManyP<rb::Vertex> fmvertex(slices, evt, fVertexLabel);

    // If there are no vertices or no prongs in this event
    // move on.

    if( ((fProngInput==true&&fmprong3D.isValid())||(fProngInput==false&&fmtrack3D.isValid())) || fmvertex.isValid() ) {

      int nSli = slices->size();
      for(int iSli = 0; iSli < nSli; ++iSli){
        //if there was nue preselection, skip slice if it failed
        if(rb::IsFiltered(evt,slices,iSli, fFilterLabels) 
           && fRespectFilter )
          continue;

        art::Ptr<rb::Cluster> thisSlice(slices, iSli);

        // if this slice is noise skip it
        if( thisSlice->IsNoise())
          continue;

//      std::vector<art::Ptr<rb::Prong>  >  prongs = fmprong3D.at(iSli);

        std::vector<art::Ptr<rb::Track>  >  tracks3D; 
        tracks3D.clear();
        if(fProngInput==false) tracks3D = fmtrack3D.at(iSli);


        std::vector<art::Ptr<rb::Prong>  >  prongs;
        prongs.clear();
        if(fProngInput==true){//prong input
          std::vector<art::Ptr<rb::Prong>  >  prongs3D = fmprong3D.at(iSli);
          for (unsigned int ip=0; ip<prongs3D.size(); ++ip){
              prongs.push_back(prongs3D[ip]);
          }
          //if instance labels are in use, combine 2d and 3d
          if ((fProng3DInstLabel != "")&&(fProng2DInstLabel != "")&&fmprong2D.isValid()){
            std::vector<art::Ptr<rb::Prong> >  prongs2D  = fmprong2D.at(iSli);
            for (unsigned int ip=0; ip<prongs2D.size(); ++ip){ 
              prongs.push_back(prongs2D[ip]);
            }
          }
        }

        if(fProngInput==false){// track input, convert tracks into prongs
          std::vector<art::Ptr<rb::Track>  >  tracks3D = fmtrack3D.at(iSli);
          for (unsigned int ip=0; ip<tracks3D.size(); ++ip){
              prongs.push_back(tracks3D[ip]);
          }
          if(fmtrack2D.isValid()){ 
            std::vector<art::Ptr<rb::Track>  >  tracks2D = fmtrack2D.at(iSli);
            for (unsigned int ip=0; ip<tracks2D.size(); ++ip){
                prongs.push_back(tracks2D[ip]);
            }
          }
        }

        int nPro = prongs.size();

        std::vector< rb::Prong > newProngs;
        std::vector< int >  assProng;

        std::vector<art::Ptr<rb::Vertex>  >  vertices   = fmvertex.at(iSli);
        if( vertices.empty() )
          continue;
        // Right now we only use fuzzyk for vertices which makes
        // one vertex per slice. If there comes along another
        // vertex finder that makes more than one per slice,
        // the hope is that the 0th one will be the primary vertex.
        art::Ptr<rb::Vertex> vtx = vertices[0];

        // Make a list of hits in slice that were not in
        // any prong. These will be subject to a PE cut
        // since hits not on prongs are likely to be noise

        rb::Cluster orphans = *(thisSlice);
        if( fApplyPEThreshold ){
          for( int iPro = 0; iPro < nPro; ++iPro){
            rb::Prong prong = *(prongs[iPro]);
            if( prong.Is2D() && !fConsider2D)
              continue;
            int nc = prong.NCell();
            for( int ic = 0; ic < nc; ++ic)
              orphans.RemoveHit( prong.Cell(ic));
          }

          for(int io = 0; io < (int)orphans.NCell(); ++io)
            if( orphans.Cell(io)->PE() > fPEThreshold ){
              orphans.RemoveHit( io );
              io = io -1;
            }
        }
        // We should now have a list of below threshold orphan hits

        for( int iPro = 0; iPro < nPro; ++iPro){

          rb::Prong thisProng  = *(prongs[iPro]);
          if( thisProng.Is2D() && !fConsider2D )
            continue;

          if( fRecluster ){

            const rb::Cluster showerClust = (rb::Cluster)thisProng;

            // Get the cluster that is the slice minus this prong
            rb::Cluster lesserSlice = (showerClust.NCell() < thisSlice->NCell()) ? thisSlice->Exclude( &showerClust ) : rb::Cluster();
            if( fApplyPEThreshold ){
              int norphans = orphans.NCell();
              for(int io = 0; io < norphans; ++io)
                lesserSlice.RemoveHit( orphans.Cell(io) );
            }

            // Get all the hits in this cluster and sort by plane-
            art::PtrVector<rb::CellHit> sliceHits = lesserSlice.AllCells();

            rb::SortByPlaneAndCell( sliceHits );

            // Reclustering is done to grab electron hits that may
            // have been missed by the prong making algorithm.
            // To recluster, a series of two collinear cylinders is
            // defined. The first cylinder captures the MIP behav-
            // -iour of the electron; the second cylinder captures the
            // EM shower. Their radii and lengths are set in .fcl file


            TVector3 dir   = thisProng.Dir();
            TVector3 start = thisProng.Start();
            TVector3 stop  = start + thisProng.TotalLength()*dir;
            int stopPlane, startPlane, vtxPlane;

            if( (start - vtx->GetXYZ()).Mag() > (stop - vtx->GetXYZ()).Mag() ){
              // vertex is closer to the stop point. Reverse
              // the shower direction.
              TVector3 temp = start;
              start = stop;
              stop  = temp;
              thisProng.SetDir( -1*dir);
              thisProng.SetStart( start );
            }

            startPlane = thisProng.MinPlane();
            stopPlane  = thisProng.MaxPlane();


            // There is a better way of determining shower start and stop planes
            // than just using the cluster min and max planes.
            // And that is to use geo::CellID as follows. However, the following
            // is being kept commented, to agree with RecoJMShower for now. (2014/07/11)

            // const geo::CellUniqueId cida = fGeom->CellId( start.X(), start.Y(),
            //                                          start.Z(), 1, 1, 0, 1 );
            // if( cida )
            //   startPlane = fGeom->getPlaneID( cida );
            // else
            //   startPlane = thisProng.MinPlane();

            // const geo::CellUniqueId cidb = fGeom->CellId( stop.X(), stop.Y(),
            //                                          stop.Z(), 1, 1, 0, 1 );
            // if( cidb )
            //   stopPlane = fGeom->getPlaneID( cidb );
            // else
            //   stopPlane = thisProng.MaxPlane();

            const geo::CellUniqueId cidv = fGeom->CellId( vtx->GetX(), 
                                                          vtx->GetY(),
                                                          vtx->GetZ(), 
                                                          1, 1, 0, 1 );

            if( cidv )
              vtxPlane = fGeom->getPlaneID( cidv );
            else
              vtxPlane = thisProng.MinPlane();

            if( abs(vtxPlane - startPlane) > abs(vtxPlane-stopPlane))
              std::swap(startPlane, stopPlane);


            float cellW = 2*fGeom->Plane(1)->Cell(1)->HalfW();

            for( auto hit : sliceHits ){

              int thisPlane = hit->Plane();
              int diffFromStart = thisPlane - startPlane;

              if( startPlane > stopPlane ){
                diffFromStart = -diffFromStart;
              }

              if( thisPlane < std::min(startPlane, stopPlane) ||
                  thisPlane > std::max(startPlane, stopPlane) )
                continue;

              float distToCore = fDeconvolveAlg->DistanceToCore( hit, 
                                                                 thisProng );

              if( distToCore == 9999 )
                continue;

              if( (diffFromStart <  fMIPPlane  &&
                   distToCore <= fMIPRadius*cellW) ||
                  (diffFromStart >= fMIPPlane && 
                   distToCore <= fMaxRadius*cellW))
                thisProng.Add( hit );
              
            }// end loop over relevant cell hits
          }// end if recluster
          newProngs.push_back( thisProng );
          assProng.push_back( iPro );
        }// end loop over prongs

        std::vector< float > totpe;
        for( int iPro = 0; iPro < nPro; ++iPro){
          totpe.push_back( prongs[iPro]->TotalPE() );
        }
        if( fDeconvolve )
          fDeconvolveAlg->GetWeights( newProngs, totpe );

        int nNewShower = newProngs.size();

        for( int i = 0; i < nNewShower; ++i){
          rb::Shower newShower(newProngs[i]);
          newShower.SetID(i + (100*iSli) );
          newShower.SetTotalLength( newProngs[i].TotalLength() );
          reclustCol->push_back( newShower );
          util::CreateAssn( *this, evt, *reclustCol, 
                            thisSlice ,*sliceAssassin);
          util::CreateAssn( *this, evt, *reclustCol, 
                            prongs[ assProng[i] ] ,*prongAssassin);
        }

      }// end loop over slices

    }// end if slice-vertex or slice-prong associations are invalid


    evt.put( std::move(reclustCol));
    evt.put( std::move(sliceAssassin));
    evt.put( std::move(prongAssassin));
  }// End of producer


  /**********************************************************************/
  ////////////////////////////////////////////////////////////////////////

  DEFINE_ART_MODULE(Recluster)

}// end of namespace slid