ClusterCrawlerAlg.h

Go to the documentation of this file.

// ClusterCrawlerAlg.h
//
// ClusterCrawlerAlg class
//
// Bruce Baller
//
#ifndef CLUSTERCRAWLERALG_H
#define CLUSTERCRAWLERALG_H


// C/C++ standard libraries
#include <array>
#include <map>
#include <vector>
#include <memory> // std::move()
#include <utility> // std::pair<>

// framework libraries
#include "fhiclcpp/ParameterSet.h" 
#include "art/Framework/Services/Registry/ServiceHandle.h" 

// LArSoft libraries
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "larcore/Geometry/Geometry.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardata/DetectorInfoServices/LArPropertiesService.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "larreco/RecoAlg/CCHitFinderAlg.h"
#include "larreco/RecoAlg/LinFitAlg.h"


namespace cluster {
  
  class ClusterCrawlerAlg {
    public:
    
    // some functions to handle the CTP_t type
    typedef unsigned int CTP_t;
    
    static constexpr unsigned int CTPpad = 1000; // alignment for CTP sub-items
    static CTP_t EncodeCTP
      (unsigned int cryo, unsigned int tpc, unsigned int plane)
      { return cryo * CTPpad*CTPpad + tpc * CTPpad + plane; }
    static CTP_t EncodeCTP(const geo::PlaneID& planeID)
      { return EncodeCTP(planeID.Cryostat, planeID.TPC, planeID.Plane); }
    static geo::PlaneID DecodeCTP(CTP_t CTP)
      { return { CTP / (CTPpad*CTPpad), CTP / CTPpad % CTPpad, CTP % CTPpad }; }
    
    
    struct ClusterStore {
      short ID;         // Cluster ID. ID < 0 = abandoned cluster
      short ProcCode;   // Processor code for debugging
      short StopCode;   // code for the reason for stopping cluster tracking
      CTP_t CTP;        // Cryostat/TPC/Plane code
      float BeginSlp;   // beginning slope (= DS end = high wire number)
      float BeginSlpErr; // error
      float BeginAng;
      unsigned int BeginWir;   // begin wire
      float BeginTim;   // begin time
      float BeginChg;   // beginning average charge
                        float BeginChgNear; // charge near the cluster Begin
                        short BeginVtx;         // ID of the Begin vertex
      float EndSlp;     // end slope (= US end = low  wire number)
      float EndAng;
      float EndSlpErr;  // error
      unsigned int EndWir;     // end wire
      float EndTim;     // ending time
      float EndChg;     // ending average charge
                        float EndChgNear; // charge near the cluster End
      short EndVtx;     // ID of the end vertex
      std::vector<unsigned int> tclhits; // hits on the cluster
    }; // ClusterStore

    struct VtxStore {
      float Wire;
                        float WireErr;
      float Time;
                        float TimeErr;
      unsigned short NClusters;  // = 0 for abandoned vertices
                        float ChiDOF;
      short Topo;                       // 1 = US-US, 2 = US-DS, 3 = DS-US, 4 = DS-DS, 5 = Star,
                                                                                                // 6 = hammer, 7 = vtx3clustermatch, 8 = vtx3clustersplit, 9 = FindTrajVertices
      CTP_t CTP;
      bool Fixed;                 // Vertex position fixed (should not be re-fit)
    };
    
    struct Vtx3Store {
      std::array<short, 3> Ptr2D; // pointers to 2D vertices in each plane
      float X;                    // x position
      float XErr;                 // x position error
      float Y;                    // y position
      float YErr;                 // y position error
      float Z;                    // z position
      float ZErr;                 // z position error
      short Wire;                 // wire number for an incomplete 3D vertex
      unsigned short CStat;
      unsigned short TPC;
      unsigned short ProcCode;
    };
    
    ClusterCrawlerAlg(fhicl::ParameterSet const& pset);

    virtual void reconfigure(fhicl::ParameterSet const& pset);
    void RunCrawler(std::vector<recob::Hit> const& srchits);

    
    std::vector<short> const& GetinClus() const {return inClus; }
    
    std::vector<recob::Hit>&& YieldHits() { return std::move(fHits); }
    
    std::vector<ClusterStore> const& GetClusters() const { return tcl; }
    
    std::vector<VtxStore> const& GetEndPoints() const { return vtx; }
    
    std::vector<Vtx3Store> const& GetVertices() const { return vtx3; }
    
    
    void ClearResults();
    
        
    static bool SortByMultiplet(recob::Hit const& a, recob::Hit const& b);
    

    private:
    
    unsigned short fNumPass;                 
    std::vector<unsigned short> fMaxHitsFit; 
    std::vector<unsigned short> fMinHits;    
    std::vector<unsigned short> fNHitsAve;   
    std::vector<float> fChiCut;     
    std::vector<float> fKinkChiRat; 
    std::vector<float> fKinkAngCut; 
    std::vector<float> fChgCut;     
    std::vector<unsigned short> fMaxWirSkip; 
    std::vector<unsigned short> fMinWirAfterSkip; 
    std::vector<bool> fDoMerge;     
    std::vector<float> fTimeDelta;  
    std::vector<float> fMergeChgCut;  
    std::vector<bool> fFindVertices;    
    std::vector<bool> fLACrawl;    
                bool fFindHammerClusters;                                        
  //  bool fFindVLAClusters;                                     ///< look for Very Large Angle clusters
    bool fRefineVertexClusters;
                
    std::vector<float> fMinAmp;                                                                 

    float fKillGarbageClusters;
    bool fChkClusterDS;
    bool fVtxClusterSplit;
    bool fFindStarVertices;
  //  bool fFindTrajVertices;

    // global cuts and parameters 
    float fHitErrFac;   
    float fHitMinAmp;   
    float fClProjErrFac;   
    float fMinHitFrac;
    float fLAClusAngleCut;  
                unsigned short fLAClusMaxHitsFit; 
    float fLAClusSlopeCut;
    bool fMergeAllHits;
    float fHitMergeChiCut; 
    float fMergeOverlapAngCut;   
    unsigned short fAllowNoHitWire;
                float fVertex2DCut;     
    float fVertex2DWireErrCut;
    float fVertex3DCut;   

    int fDebugPlane;
    int fDebugWire;  
    int fDebugHit;   
    
    // Wires that have been determined by some filter (e.g. NoiseFilter) to be good
    std::vector<geo::WireID> fFilteredWires;
 
    // these variables define the cluster used during crawling
    float clpar[3];     
    float clparerr[2];  
    float clChisq;     
    float fAveChg;  
  //  float fChgRMS;  ///< average charge RMS at leading edge of cluster
    float fChgSlp;  
    float fAveHitWidth; 
    
    bool prt;
    bool vtxprt;
    unsigned short NClusters;
    
    art::ServiceHandle<geo::Geometry> geom;
    
    std::vector<recob::Hit> fHits; 
    std::vector<short> inClus;    
    std::vector<bool> mergeAvailable; 
    std::vector< ClusterStore > tcl; 
    std::vector< VtxStore > vtx; 
    std::vector< Vtx3Store > vtx3; 
    
    trkf::LinFitAlg fLinFitAlg;

    float clBeginSlp;  
    float clBeginAng;
    float clBeginSlpErr; 
    unsigned int clBeginWir;  
    float clBeginTim;  
    float clBeginChg;  
                float clBeginChgNear; 
    float clEndSlp;    
    float clEndAng;
    float clEndSlpErr; 
    unsigned int clEndWir;    
    float clEndTim;    
    float clEndChg;    
                float clEndChgNear;  
    short clStopCode;     
    short clProcCode;     
    CTP_t clCTP;        
                bool clLA;                                      
    
    unsigned int fFirstWire;    
    unsigned int fFirstHit;     
    unsigned int fLastWire;      
    unsigned int cstat;         // the current cryostat
    unsigned int tpc;         // the current TPC
    unsigned int plane;         // the current plane
    unsigned int fNumWires;   // number of wires in the current plane
    unsigned int fMaxTime;    // number of time samples in the current plane
    float fScaleF;     
    
    unsigned short pass;

    // vector of pairs of first (.first) and last+1 (.second) hit on each wire
    // in the range fFirstWire to fLastWire. A value of -2 indicates that there
    // are no hits on the wire. A value of -1 indicates that the wire is dead
    std::vector< std::pair<int, int> > WireHitRange;
    
    std::vector<unsigned int> fcl2hits;  
    std::vector<float> chifits;   
    std::vector<short> hitNear;   
    
    std::vector<float> chgNear; 
                float fChgNearWindow;           
                float fChgNearCut;                              

    std::string fhitsModuleLabel;
    // ******** crawling routines *****************

    // Loops over wires looking for seed clusters
    void ClusterLoop();
    // Returns true if the hits on a cluster have a consistent width
    bool ClusterHitsOK(short nHitChk);
    // Finds a hit on wire kwire, adds it to the cluster and re-fits it
    void AddHit(unsigned int kwire, bool& HitOK, bool& SigOK);
    // Finds a hit on wire kwire, adds it to a LargeAngle cluster and re-fits it
    void AddLAHit(unsigned int kwire, bool& ChkCharge, bool& HitOK, bool& SigOK);
    // Fits the cluster hits in fcl2hits to a straight line
    void FitCluster();
    // Fits the charge of the cluster hits in fcl2hits
    void FitClusterChg();
     // Fits the middle of a temporary cluster it1 using hits iht to iht + nhit
    void FitClusterMid(unsigned short it1, unsigned int iht, short nhit);
    void FitClusterMid(std::vector<unsigned int>& hitVec, unsigned int iht, short nhit);
    // Crawls along a trail of hits UpStream
    void CrawlUS();
    // Crawls along a trail of hits UpStream - Large Angle version
    void LACrawlUS();
    
    void KillGarbageClusters();

    // ************** cluster merging routines *******************

    // Compares two cluster combinations to see if they should be merged
    void ChkMerge();
    // Checks merge for cluster cl2 within the bounds of cl1
    void ChkMerge12(unsigned short it1, unsigned short it2, bool& didit);
    // Merges clusters cl1 and cl2
    void DoMerge(unsigned short it1, unsigned short it2, short ProcCode);
      
    // ************** hit merging routines *******************
    // check the number of nearby hits to the ones added to the current cluster.
    // If there are too many, merge the hits and re-fit
    void ChkClusterNearbyHits(bool prt);
    // merge the hits in a multiplet into one hit
    void MergeHits(const unsigned int theHit, bool& didMerge);
    void FixMultipletLocalIndices
      (size_t begin, size_t end, short int multiplicity = -1);

    // ************** cluster finish routines *******************

    // Check the fraction of wires with hits
    void CheckClusterHitFrac(bool prt);

    // Try to extend clusters downstream
    void ChkClusterDS();

    // Try to merge overlapping clusters
    void MergeOverlap();
    
    void MakeClusterObsolete(unsigned short icl);
    void RestoreObsoleteCluster(unsigned short icl);
    
    void RemoveObsoleteHits();

    // ************** 2D vertex routines *******************

    // Find 2D vertices
    void FindVertices();
    // Find 2D star topology vertices
    void FindStarVertices();
    // check a vertex (vw, fvt) made with clusters it1, and it2 against the
    // vector of existing clusters
    void ChkVertex(float fvw, float fvt, unsigned short it1, unsigned short it2, short topo);
    // try to attach a cluster to an existing vertex
    void ClusterVertex(unsigned short it2);
    // try to attach a cluster to a specified vertex
    void VertexCluster(unsigned short ivx);
    // Refine cluster ends near vertices
    void RefineVertexClusters(unsigned short ivx);
    // Split clusters that cross a vertex
    bool VtxClusterSplit();
    // returns true if a vertex is encountered while crawling
    bool CrawlVtxChk(unsigned int kwire);
    // returns true if this cluster is between a vertex and another
    // cluster that is associated with the vertex
    bool CrawlVtxChk2();
    // use a vertex constraint to start a cluster
    void VtxConstraint(unsigned int iwire, unsigned int ihit, unsigned int jwire, unsigned int& useHit, bool& doConstrain);
    // fit the vertex position
    void FitVtx(unsigned short iv);
    // weight and fit all vertices
    void FitAllVtx(CTP_t inCTP);

    // ************** 3D vertex routines *******************

    // match vertices between planes
    void VtxMatch(geo::TPCID const& tpcid);
    // Match clusters to endpoints using 3D vertex information
    void Vtx3ClusterMatch(geo::TPCID const& tpcid);
    // split clusters using 3D vertex information
    void Vtx3ClusterSplit(geo::TPCID const& tpcid);
                // look for a long cluster that stops at a short cluster in two views
                void FindHammerClusters();

    // ************** utility routines *******************

    // inits everything
    void CrawlInit();
    // inits the cluster stuff
    void ClusterInit();
    // fills the wirehitrange vector for the supplied Cryostat/TPC/Plane code
    void GetHitRange(CTP_t CTP);
    // Stores cluster information in a temporary vector
    bool TmpStore();
    // Gets a temp cluster and puts it into the working cluster variables
    void TmpGet(unsigned short it1);
    // Does just what it says
    void CalculateAveHitWidth();
    // Shortens the fcl2hits, chifits, etc vectors by the specified amount
    void FclTrimUS(unsigned short nTrim);
    // Splits a cluster into two clusters at position pos. Associates the
    // new clusters with a vertex
    bool SplitCluster(unsigned short icl, unsigned short pos, unsigned short ivx);
    // Counts the number of dead wires in the range spanned by fcl2hits
    unsigned int DeadWireCount();
    unsigned int DeadWireCount(unsigned int inWire1, unsigned int inWire2);
    // return true if the pre-merged it1 and it2 clusters will meet the quality requirement
    bool ChkMergedClusterHitFrac(unsigned short it1, unsigned short it2);
    // Prints cluster information to the screen
    void PrintClusters();
    void PrintVertices();
    // check for a signal on all wires between two points
    bool ChkSignal(unsigned int iht, unsigned int jht);
    bool ChkSignal(unsigned int wire1, float time1, unsigned int wire2, float time2);
    // returns an angle-dependent scale factor for weighting fits, etc
    float AngleFactor(float slope);
    // calculate the kink angle between hits 0-2 and 3 - 5 on the leading edge of
    // the cluster under construction
    float EndKinkAngle();
    bool CheckHitDuplicates(std::string location, std::string marker = "") const;
    // Find the distance of closest approach between the end of a cluster and a (wire,tick) position
    float DoCA(short icl, unsigned short end, float vwire, float vtick);
    // Find the Chisq/DOF between the end of a cluster and a (wire,tick) position
    float ClusterVertexChi(short icl, unsigned short end, unsigned short ivx);
    // Find the Chisq/DOF between a point and a vertex
    float PointVertexChi(float wire, float tick, unsigned short ivx);
    std::string PrintHit(unsigned int iht);
    
    std::pair<size_t, size_t> FindHitMultiplet(size_t iHit) const;

    void CheckHitClusterAssociations();
    
    static bool areInSameMultiplet(recob::Hit const& first_hit, recob::Hit const& second_hit);

    
  }; // class ClusterCrawlerAlg


} // namespace cluster

#endif // ifndef CLUSTERCRAWLERALG_H