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 <utility> // std::pair<>
#include <vector>

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

// LArSoft libraries
#include "larcore/Geometry/Geometry.h"
#include "larcore/Geometry/WireReadout.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "lardataobj/RecoBase/Hit.h"
#include "larreco/RecoAlg/LinFitAlg.h"
namespace detinfo {
  class DetectorClocksData;
  class DetectorPropertiesData;
}

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;
    };

    explicit ClusterCrawlerAlg(fhicl::ParameterSet const& pset);

    void RunCrawler(detinfo::DetectorClocksData const& clock_data,
                    detinfo::DetectorPropertiesData const& det_prop,
                    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<recob::Hit> GetHits() { return 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 const> geom;
    geo::WireReadoutGeom const* wireReadoutGeom{
      &art::ServiceHandle<geo::WireReadout const>()->Get()};

    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(detinfo::DetectorPropertiesData const& det_prop, geo::TPCID const& tpcid);
    // Match clusters to endpoints using 3D vertex information
    void Vtx3ClusterMatch(detinfo::DetectorPropertiesData const& det_prop, geo::TPCID const& tpcid);
    // split clusters using 3D vertex information
    void Vtx3ClusterSplit(detinfo::DetectorPropertiesData const& det_prop, geo::TPCID const& tpcid);
    // look for a long cluster that stops at a short cluster in two views
    void FindHammerClusters(detinfo::DetectorPropertiesData const& det_prop);

    // ************** 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