Cluster3D.h

Go to the documentation of this file.

#ifndef RECO_CLUSTER3D_H
#define RECO_CLUSTER3D_H

#include <iosfwd>
#include <vector>
#include <list>
#include <set>
#include <map>
#include <unordered_map>
#include <memory>

#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "lardataobj/RecoBase/Hit.h"
#include "larreco/RecoAlg/Cluster3DAlgs/Voronoi/DCEL.h"

// Eigen
#include <Eigen/Dense>

namespace reco {
    
// First define a container object to augment the sparse 2D hit information
class ClusterHit2D
{
public:
    
    ClusterHit2D();   // Default constructor
    
private:
    
    mutable unsigned  m_statusBits;   
    mutable float     m_docaToAxis;   
    mutable float     m_arcLenToPoca; 
    float             m_xPosition;    
    float             m_timeTicks;    
    const recob::Hit* m_hit;          
    
public:
    
    enum StatusBits { SHAREDINPAIR    = 0x00080000,
                      SHAREDINTRIPLET = 0x00040000,
                      USEDINPAIR      = 0x00008000,
                      USEDINTRIPLET   = 0x00004000,
                      SHAREDINCLUSTER = 0x00000200,
                      USEDINCLUSTER   = 0x00000100,
                      USED            = 0x00000001
                    };
    
    ClusterHit2D(unsigned          statusBits,
                 float             doca,
                 float             poca,
                 float             xPosition,
                 float             timeTicks,
                 const recob::Hit& recobHit);
    
    unsigned          getStatusBits()   const {return m_statusBits;}
    float             getDocaToAxis()   const {return m_docaToAxis;}
    float             getArcLenToPoca() const {return m_arcLenToPoca;}
    float             getXPosition()    const {return m_xPosition;}
    float             getTimeTicks()    const {return m_timeTicks;}
    const recob::Hit& getHit()          const {return *m_hit;}
    
    void setStatusBit(unsigned bits)    const {m_statusBits   |= bits;}
    void clearStatusBits(unsigned bits) const {m_statusBits   &= ~bits;}
    void setDocaToAxis(float doca)      const {m_docaToAxis    = doca;}
    void setArcLenToPoca(float poca)    const {m_arcLenToPoca  = poca;}
    
    friend std::ostream& operator << (std::ostream& o, const ClusterHit2D& c);
    friend bool          operator <  (const ClusterHit2D & a, const ClusterHit2D & b);
    
};

using ClusterHit2DVec = std::vector<const reco::ClusterHit2D*>;

// Now define an object with the recob::Hit information that will comprise the 3D cluster
class ClusterHit3D
{
public:
    
    enum StatusBits { REJECTEDHIT     = 0x80000000,           
                      SKELETONHIT     = 0x10000000,           
                      EDGEHIT         = 0x20000000,           
                      SEEDHIT         = 0x40000000,           
                      MADESPACEPOINT  = 0x08000000,           
                      CONVEXHULLVTX   = 0x04000000,           
                      EXTREMEPOINT    = 0x02000000,           
                      SKELETONPOSAVE  = 0x00100000,           
                      CLUSTERVISITED  = 0x00008000,           
                      CLUSTERNOISE    = 0x00004000,           
                      CLUSTERATTACHED = 0x00002000,           
                      CLUSTERSHARED   = 0x00001000,           
                      PATHCHECKED     = 0x00000800,           
                      SELECTEDBYMST   = 0x00000100,           
                      PCAOUTLIER      = 0x00000080,           
                      HITINVIEW0      = 0x00000001,           
                      HITINVIEW1      = 0x00000002,           
                      HITINVIEW2      = 0x00000004            
                    };
        
    
    ClusterHit3D();   // Default constructor

    ClusterHit3D(size_t                          id,
                 unsigned int                    statusBits,
                 const float*                    position,
                 float                           totalCharge,
                 float                           avePeakTime,
                 float                           deltaPeakTime,
                 float                           sigmaPeakTime,
                 float                           hitChiSquare,
                 float                           docaToAxis,
                 float                           arclenToPoca,
                 const ClusterHit2DVec&          hitVec,
                 const std::vector<float>&       hitDelTSigVec,
                 const std::vector<geo::WireID>& wireIDVec);
    
    ClusterHit3D(const ClusterHit3D&);

    void initialize(size_t                          id,
                    unsigned int                    statusBits,
                    const float*                    position,
                    float                           totalCharge,
                    float                           avePeakTime,
                    float                           deltaPeakTime,
                    float                           sigmaPeakTime,
                    float                           hitChiSquare,
                    float                           docaToAxis,
                    float                           arclenToPoca,
                    const ClusterHit2DVec&          hitVec,
                    const std::vector<float>&       hitDelTSigVec,
                    const std::vector<geo::WireID>& wireIDVec);

    size_t                              getID()            const {return m_id;}
    unsigned int                        getStatusBits()    const {return m_statusBits;}
    const float*                        getPosition()      const {return m_position;}
    float                               getX()             const {return m_position[0];}
    float                               getY()             const {return m_position[1];}
    float                               getZ()             const {return m_position[2];}
    float                               getTotalCharge()   const {return m_totalCharge;}
    float                               getAvePeakTime()   const {return m_avePeakTime;}
    float                               getDeltaPeakTime() const {return m_deltaPeakTime;}
    float                               getSigmaPeakTime() const {return m_sigmaPeakTime;}
    float                               getHitChiSquare()  const {return m_hitChiSquare;}
    float                               getDocaToAxis()    const {return m_docaToAxis;}
    float                               getArclenToPoca()  const {return m_arclenToPoca;}
    const ClusterHit2DVec&              getHits()          const {return m_hitVector;}
    const std::vector<float>            getHitDelTSigVec() const {return m_hitDelTSigVec;}
    const std::vector<geo::WireID>&     getWireIDs()       const {return m_wireIDVector;}
    
    bool bitsAreSet(const unsigned int& bitsToCheck)       const {return m_statusBits & bitsToCheck;}

    void setID(const size_t& id)           const {m_id            = id;}
    void setStatusBit(unsigned bits)       const {m_statusBits   |= bits;}
    void clearStatusBits(unsigned bits)    const {m_statusBits   &= ~bits;}
    void setDocaToAxis(double doca)        const {m_docaToAxis    = doca;}
    void setArclenToPoca(double poca)      const {m_arclenToPoca  = poca;}
    void setWireID(const geo::WireID& wid) const;
    
    void setPosition(const float* pos) const {m_position[0] = pos[0]; m_position[1] = pos[1]; m_position[2] = pos[2];}

    const bool operator<(const reco::ClusterHit3D& other) const
    {
        if (m_position[2] != other.m_position[2]) return m_position[2] < other.m_position[2];
        else return m_position[0] < other.m_position[0];
    }

    const bool operator==(const reco::ClusterHit3D& other) const
    {
        return m_id == other.m_id;
    }
    
    friend std::ostream& operator << (std::ostream& o, const ClusterHit3D& c);
    //friend bool          operator <  (const ClusterHit3D & a, const ClusterHit3D & b);
    
private:
    
    mutable size_t                   m_id;                 
    mutable unsigned int             m_statusBits;         
    mutable float                    m_position[3];        
    float                            m_totalCharge;        
    float                            m_avePeakTime;        
    float                            m_deltaPeakTime;      
    float                            m_sigmaPeakTime;      
    float                            m_hitChiSquare;       
    mutable float                    m_docaToAxis;         
    mutable float                    m_arclenToPoca;       
    ClusterHit2DVec                  m_hitVector;          
    mutable std::vector<float>       m_hitDelTSigVec;      
    mutable std::vector<geo::WireID> m_wireIDVector;       
};
    
// We also need to define a container for the output of the PCA Analysis
class PrincipalComponents
{
public:

    typedef std::vector<std::vector<float>> EigenVectors;
    
    PrincipalComponents();
    
private:
    
    bool           m_svdOK;              
    int            m_numHitsUsed;        
    float          m_eigenValues[3];     
    EigenVectors   m_eigenVectors;       
    float          m_avePosition[3];     
    mutable double m_aveHitDoca;         
    
public:
    
    PrincipalComponents(bool ok, int nHits, const float* eigenValues, const EigenVectors& eigenVecs, const float* avePos, const float aveHitDoca = 9999.);
    
    bool                getSvdOK()            const {return m_svdOK;}
    int                 getNumHitsUsed()      const {return m_numHitsUsed;}
    const float*        getEigenValues()      const {return m_eigenValues;}
    const EigenVectors& getEigenVectors()     const {return m_eigenVectors;}
    const float*        getAvePosition()      const {return m_avePosition;}
    const float         getAveHitDoca()       const {return m_aveHitDoca;}
    
    void                flipAxis(size_t axis);
    void                setAveHitDoca(double doca) const {m_aveHitDoca = doca;}
    
    friend std::ostream&  operator << (std::ostream & o, const PrincipalComponents& a);
    friend bool operator < (const PrincipalComponents& a, const PrincipalComponents& b);
    
};

class Cluster3D
{
public:

    Cluster3D();  

private:

    mutable unsigned    m_statusBits;       
    PrincipalComponents m_pcaResults;       
    float               m_totalCharge;      
    float               m_startPosition[3]; 
    float               m_endPosition[3];   
    int                 m_clusterIdx;       
    
public:
    Cluster3D(unsigned                   statusBits,
              const PrincipalComponents& pcaResults,
              float                      totalCharge,
              const float*               startPosition,
              const float*               endPosition,
              int                        idx);
    
    unsigned                   getStatusBits()    const {return m_statusBits;}
    const PrincipalComponents& getPcaResults()    const {return m_pcaResults;}
    float                      getTotalCharge()   const {return m_totalCharge;}
    const float*               getStartPosition() const {return m_startPosition;}
    const float*               getEndPosition()   const {return m_endPosition;}
    int                        getClusterIdx()    const {return m_clusterIdx;}
    
    void setStatusBit(unsigned bits)    const {m_statusBits |= bits;}
    void clearStatusBits(unsigned bits) const {m_statusBits &= ~bits;}
    
    Cluster3D            operator +  (Cluster3D);
    friend std::ostream& operator << (std::ostream& o, const Cluster3D& c);
    friend bool          operator <  (const Cluster3D & a, const Cluster3D & b);
};

class RecobClusterParameters
{
public:
    RecobClusterParameters() : m_startTime(999999.),
    m_sigmaStartTime(1.),
    m_endTime(0.),
    m_sigmaEndTime(1.),
    m_totalCharge(0.),
    m_startWire(9999999),
    m_endWire(0),
    m_plane(100),
    m_view(geo::kUnknown)
    {
        m_hitVector.clear();
    }
    
    void UpdateParameters(const reco::ClusterHit2D* hit);
    
    float           m_startTime;
    float           m_sigmaStartTime;
    float           m_endTime;
    float           m_sigmaEndTime;
    float           m_totalCharge;
    unsigned int    m_startWire;
    unsigned int    m_endWire;
    unsigned int    m_plane;
    geo::View_t     m_view;
    ClusterHit2DVec m_hitVector;
};
    
    
using Hit2DListPtr             = std::list<const reco::ClusterHit2D*>;
using HitPairListPtr           = std::list<const reco::ClusterHit3D*>;
using HitPairSetPtr            = std::set<const reco::ClusterHit3D*>;
using HitPairListPtrList       = std::list<HitPairListPtr>;
using HitPairClusterMap        = std::map<int, HitPairListPtr>;
using HitPairList              = std::list<std::unique_ptr<reco::ClusterHit3D>>;
    
using PCAHitPairClusterMapPair = std::pair<reco::PrincipalComponents, reco::HitPairClusterMap::iterator>;
using PlaneToClusterParamsMap  = std::map<size_t, RecobClusterParameters>;
using EdgeTuple                = std::tuple<const reco::ClusterHit3D*,const reco::ClusterHit3D*,double>;
using EdgeList                 = std::list<EdgeTuple>;
using Hit3DToEdgePair          = std::pair<const reco::ClusterHit3D*, reco::EdgeList>;
using Hit3DToEdgeMap           = std::unordered_map<const reco::ClusterHit3D*, reco::EdgeList>;
using Hit2DToHit3DListMap      = std::unordered_map<const reco::ClusterHit2D*, reco::HitPairListPtr>;
//using VertexPoint              = Eigen::Vector3f;
//using VertexPointList          = std::list<Eigen::Vector3f>;
    
using ProjectedPoint           = std::tuple<float, float, const reco::ClusterHit3D*>;          
using ProjectedPointList       = std::list<ProjectedPoint>;
using ConvexHullKinkTuple      = std::tuple<ProjectedPoint, Eigen::Vector2f, Eigen::Vector2f>; 
using ConvexHullKinkTupleList  = std::list<ConvexHullKinkTuple>;
    
class ConvexHull
{
public:
    ConvexHull()
    {
        fProjectedPointList.clear(),
        fConvexHullEdgeMap.clear(),
        fConvexHullEdgeList.clear(),
        fConvexHullExtremePoints.clear(),
        fConvexHullKinkPoints.clear();
    }
    
    void clear()
    {
        fProjectedPointList.clear(),
        fConvexHullEdgeMap.clear(),
        fConvexHullEdgeList.clear(),
        fConvexHullExtremePoints.clear(),
        fConvexHullKinkPoints.clear();
    }
    
    reco::ProjectedPointList&      getProjectedPointList()      {return fProjectedPointList;}
    reco::Hit3DToEdgeMap&          getConvexHullEdgeMap()       {return fConvexHullEdgeMap;}
    reco::EdgeList&                getConvexHullEdgeList()      {return fConvexHullEdgeList;}
    reco::ProjectedPointList&      getConvexHullExtremePoints() {return fConvexHullExtremePoints;}
    reco::ConvexHullKinkTupleList& getConvexHullKinkPoints()    {return fConvexHullKinkPoints;}
    
    
private:
    reco::ProjectedPointList      fProjectedPointList;      
    reco::Hit3DToEdgeMap          fConvexHullEdgeMap;       
    reco::EdgeList                fConvexHullEdgeList;      
    reco::ProjectedPointList      fConvexHullExtremePoints; 
    reco::ConvexHullKinkTupleList fConvexHullKinkPoints;    
};

class ClusterParameters;
using ClusterParametersList = std::list<ClusterParameters>;
    
class ClusterParameters
{
public:
    ClusterParameters()
    {
        fClusterParams.clear();
        fHitPairListPtr.clear();
        fHit2DToHit3DListMap.clear();
        fHit3DToEdgeMap.clear();
        fBestHitPairListPtr.clear();
        fBestEdgeList.clear();
        fConvexHull.clear();
        fFaceList.clear();
        fVertexList.clear();
        fHalfEdgeList.clear();
        fClusterParameters.clear();
    }
    
    ClusterParameters(reco::HitPairClusterMap::iterator& mapItr) : fHitPairListPtr(mapItr->second)
    {
        fClusterParams.clear();
        fHit2DToHit3DListMap.clear();
        fHit3DToEdgeMap.clear();
        fBestHitPairListPtr.clear();
        fBestEdgeList.clear();
        fConvexHull.clear();
        fFaceList.clear();
        fVertexList.clear();
        fHalfEdgeList.clear();
    }
    
    ClusterParameters(reco::HitPairListPtr& hitList) : fHitPairListPtr(hitList)
    {
        fClusterParams.clear();
        fHit2DToHit3DListMap.clear();
        fHit3DToEdgeMap.clear();
        fBestHitPairListPtr.clear();
        fBestEdgeList.clear();
        fConvexHull.clear();
        fFaceList.clear();
        fVertexList.clear();
        fHalfEdgeList.clear();
    }
    
    ClusterParametersList& daughterList() {return fClusterParameters;}
    
    void UpdateParameters(const reco::ClusterHit2D* hit)
    {
        fClusterParams[hit->getHit().WireID().Plane].UpdateParameters(hit);
    }
    
    void addHit3D(const reco::ClusterHit3D* hit3D)
    {
        fHitPairListPtr.emplace_back(hit3D);
        
        for(const auto& hit2D : hit3D->getHits())
            if (hit2D) fHit2DToHit3DListMap[hit2D].emplace_back(hit3D);
    }
    
    void fillHit2DToHit3DListMap()
    {
        for(const auto& hit3D : fHitPairListPtr)
        {
            for(const auto& hit2D : hit3D->getHits())
                if (hit2D) fHit2DToHit3DListMap[hit2D].emplace_back(hit3D);
        }
    }
    
    reco::PlaneToClusterParamsMap& getClusterParams()       {return fClusterParams;}
    reco::Hit2DToHit3DListMap&     getHit2DToHit3DListMap() {return fHit2DToHit3DListMap;}
    reco::HitPairListPtr&          getHitPairListPtr()      {return fHitPairListPtr;}
    reco::PrincipalComponents&     getFullPCA()             {return fFullPCA;}
    reco::PrincipalComponents&     getSkeletonPCA()         {return fSkeletonPCA;}
    reco::Hit3DToEdgeMap&          getHit3DToEdgeMap()      {return fHit3DToEdgeMap;}
    reco::HitPairListPtr&          getBestHitPairListPtr()  {return fBestHitPairListPtr;}
    reco::EdgeList&                getBestEdgeList()        {return fBestEdgeList;}
    reco::ConvexHull&              getConvexHull()          {return fConvexHull;}
    dcel2d::FaceList&              getFaceList()            {return fFaceList;}
    dcel2d::VertexList&            getVertexList()          {return fVertexList;}
    dcel2d::HalfEdgeList&          getHalfEdgeList()        {return fHalfEdgeList;}

    friend bool operator < (const ClusterParameters &a, const ClusterParameters& b)
    {
        return a.fHitPairListPtr.size() > b.fHitPairListPtr.size();
    }

private:
    PlaneToClusterParamsMap       fClusterParams;
    reco::HitPairListPtr          fHitPairListPtr;      // This contains the list of 3D hits in the cluster
    reco::Hit2DToHit3DListMap     fHit2DToHit3DListMap; // Provides a mapping between 2D hits and 3D hits they make
    reco::PrincipalComponents     fFullPCA;             // PCA run over full set of 3D hits
    reco::PrincipalComponents     fSkeletonPCA;         // PCA run over just the "skeleton" 3D hits
    reco::Hit3DToEdgeMap          fHit3DToEdgeMap;
    reco::HitPairListPtr          fBestHitPairListPtr;
    reco::EdgeList                fBestEdgeList;        // This has become multiuse... really need to split it up
    reco::ConvexHull              fConvexHull;          // Convex hull object
    dcel2d::FaceList              fFaceList;            // Keeps track of "faces" from Voronoi Diagram
    dcel2d::VertexList            fVertexList;          // Keeps track of "vertices" from Voronoi Diagram
    dcel2d::HalfEdgeList          fHalfEdgeList;        // Keeps track of "halfedges" from Voronoi Diagram
    ClusterParametersList         fClusterParameters;   // For possible daughter clusters
};

using ClusterToHitPairSetPair = std::pair<reco::ClusterParameters*,HitPairSetPtr>;
using ClusterToHitPairSetMap  = std::unordered_map<reco::ClusterParameters*,HitPairSetPtr>;
using Hit2DToHit3DSetMap      = std::unordered_map<const reco::ClusterHit2D*,HitPairSetPtr>;
using Hit2DToClusterMap       = std::unordered_map<const reco::ClusterHit2D*,ClusterToHitPairSetMap>;
    
}

#endif //RECO_CLUSTER3D_H