Cluster3D.h

Go to the documentation of this file.

#ifndef RECO_CLUSTER3D_H
#define RECO_CLUSTER3D_H

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

#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "larreco/RecoAlg/Cluster3DAlgs/Voronoi/DCEL.h"
namespace recob {
  class Hit;
}

// Eigen
#include <Eigen/Core>

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;             
    geo::WireID m_wireID;          
    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 geo::WireID& wireID,
                 const recob::Hit* recobHit);

    ClusterHit2D(const ClusterHit2D&);
    ClusterHit2D& operator=(ClusterHit2D const&);

    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 geo::WireID& WireID() const { return m_wireID; }
    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; }

    void setHit(const recob::Hit* hit) { m_hit = hit; }

    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 Eigen::Vector3f& position,
                 float totalCharge,
                 float avePeakTime,
                 float deltaPeakTime,
                 float sigmaPeakTime,
                 float hitChiSquare,
                 float overlapFraction,
                 float chargeAsymmetry,
                 float docaToAxis,
                 float arclenToPoca,
                 const ClusterHit2DVec& hitVec,
                 const std::vector<float>& hitDelTSigVec,
                 const std::vector<geo::WireID>& wireIDVec);

    ClusterHit3D(const ClusterHit3D&);
    ClusterHit3D& operator=(ClusterHit3D const&);

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

    size_t getID() const { return fID; }
    unsigned int getStatusBits() const { return fStatusBits; }
    const Eigen::Vector3f getPosition() const { return fPosition; }
    float getX() const { return fPosition[0]; }
    float getY() const { return fPosition[1]; }
    float getZ() const { return fPosition[2]; }
    float getTotalCharge() const { return fTotalCharge; }
    float getAvePeakTime() const { return fAvePeakTime; }
    float getDeltaPeakTime() const { return fDeltaPeakTime; }
    float getSigmaPeakTime() const { return fSigmaPeakTime; }
    float getHitChiSquare() const { return fHitChiSquare; }
    float getOverlapFraction() const { return fOverlapFraction; }
    float getChargeAsymmetry() const { return fChargeAsymmetry; }
    float getDocaToAxis() const { return fDocaToAxis; }
    float getArclenToPoca() const { return fArclenToPoca; }
    const ClusterHit2DVec& getHits() const { return fHitVector; }
    const std::vector<float> getHitDelTSigVec() const { return fHitDelTSigVec; }
    const std::vector<geo::WireID>& getWireIDs() const { return fWireIDVector; }

    ClusterHit2DVec& getHits() { return fHitVector; }

    bool bitsAreSet(const unsigned int& bitsToCheck) const { return fStatusBits & bitsToCheck; }

    void setID(const size_t& id) const { fID = id; }
    void setStatusBit(unsigned bits) const { fStatusBits |= bits; }
    void clearStatusBits(unsigned bits) const { fStatusBits &= ~bits; }
    void setDocaToAxis(double doca) const { fDocaToAxis = doca; }
    void setArclenToPoca(double poca) const { fArclenToPoca = poca; }
    void setWireID(const geo::WireID& wid) const;

    void setPosition(const Eigen::Vector3f& pos) const { fPosition = pos; }

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

    bool operator==(const reco::ClusterHit3D& other) const { return fID == other.fID; }

    friend std::ostream& operator<<(std::ostream& o, const ClusterHit3D& c);
    //friend bool          operator <  (const ClusterHit3D & a, const ClusterHit3D & b);

  private:
    mutable size_t fID;                
    mutable unsigned int fStatusBits;  
    mutable Eigen::Vector3f fPosition; 
    float fTotalCharge;                
    float fAvePeakTime;                
    float fDeltaPeakTime;              
    float fSigmaPeakTime;              
    float fHitChiSquare;               
    float fOverlapFraction;            
    float fChargeAsymmetry;            
    mutable float fDocaToAxis;         
    mutable float fArclenToPoca;       
    ClusterHit2DVec fHitVector;        
    mutable std::vector<float> fHitDelTSigVec;      
    mutable std::vector<geo::WireID> fWireIDVector; 
  };

  // We also need to define a container for the output of the PCA Analysis
  class PrincipalComponents {
  public:
    using EigenValues = Eigen::Vector3f;
    using EigenVectors = Eigen::Matrix3f;

    PrincipalComponents();

  private:
    bool m_svdOK;                  
    int m_numHitsUsed;             
    EigenValues m_eigenValues;     
    EigenVectors m_eigenVectors;   
    Eigen::Vector3f m_avePosition; 
    mutable double m_aveHitDoca;   

  public:
    PrincipalComponents(bool ok,
                        int nHits,
                        const EigenValues& eigenValues,
                        const EigenVectors& eigenVecs,
                        const Eigen::Vector3f& avePos,
                        const float aveHitDoca = 9999.);

    bool getSvdOK() const { return m_svdOK; }
    int getNumHitsUsed() const { return m_numHitsUsed; }
    const EigenValues& getEigenValues() const { return m_eigenValues; }
    const EigenVectors& getEigenVectors() const { return m_eigenVectors; }
    const Eigen::Vector3f& getAvePosition() const { return m_avePosition; }
    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(geo::PlaneID())
      , 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;
    geo::PlaneID 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<reco::ClusterHit3D>;
  //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(), fConvexHullPointList.clear(), fConvexHullEdgeMap.clear(),
        fConvexHullEdgeList.clear(), fConvexHullExtremePoints.clear(),
        fConvexHullKinkPoints.clear();
    }

    void clear()
    {
      fProjectedPointList.clear(), fConvexHullPointList.clear(), fConvexHullEdgeMap.clear(),
        fConvexHullEdgeList.clear(), fConvexHullExtremePoints.clear(),
        fConvexHullKinkPoints.clear();
    }

    reco::ProjectedPointList& getProjectedPointList() { return fProjectedPointList; }
    reco::ProjectedPointList& getConvexHullPointList() { return fConvexHullPointList; }
    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::ProjectedPointList fConvexHullPointList; 
    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->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