Inheritance diagram for lar_cluster3d::MSTPathFinder:

Public Member Functions
	MSTPathFinder (const fhicl::ParameterSet &)
	Constructor. More...

	~MSTPathFinder ()
	Destructor. More...

void	configure (fhicl::ParameterSet const &pset) override

void	initializeHistograms (art::TFileDirectory &) override
	Interface for initializing histograms if they are desired Note that the idea is to put hisgtograms in a subfolder. More...

void	ModifyClusters (reco::ClusterParametersList &) const override
	Scan an input collection of clusters and modify those according to the specific implementing algorithm. More...

float	getTimeToExecute () const override
	If monitoring, recover the time to execute a particular function. More...

Public Member Functions inherited from lar_cluster3d::IClusterModAlg
virtual	~IClusterModAlg () noexcept=default
	Virtual Destructor. More...

virtual void	configure (const fhicl::ParameterSet &)=0
	Interface for configuring the particular algorithm tool. More...

Private Types
enum	TimeValues { BUILDTHREEDHITS = 0, BUILDHITTOHITMAP = 1, RUNDBSCAN = 2, BUILDCLUSTERINFO = 3, PATHFINDING = 4, NUMTIMEVALUES }
	enumerate the possible values for time checking if monitoring timing More...

using	BestNodeTuple = std::tuple< const reco::ClusterHit3D *, float, float >

using	BestNodeMap = std::unordered_map< const reco::ClusterHit3D *, BestNodeTuple >

using	MinMaxPoints = std::pair< reco::ProjectedPoint, reco::ProjectedPoint >
	Add ability to build the convex hull (these needs to be split out! ) More...

using	MinMaxPointPair = std::pair< MinMaxPoints, MinMaxPoints >

Private Member Functions
void	RunPrimsAlgorithm (const reco::HitPairListPtr &, kdTree::KdTreeNode &, reco::ClusterParametersList &) const
	Driver for Prim's algorithm. More...

void	PruneAmbiguousHits (reco::ClusterParameters &, reco::Hit2DToClusterMap &) const
	Prune the obvious ambiguous hits. More...

void	FindBestPathInCluster (reco::ClusterParameters &) const
	Algorithm to find the best path through the given cluster. More...

reco::HitPairListPtr	DepthFirstSearch (const reco::EdgeTuple &, const reco::Hit3DToEdgeMap &, float &) const
	a depth first search to find longest branches More...

void	FindBestPathInCluster (reco::ClusterParameters &, kdTree::KdTreeNode &) const
	Alternative version of FindBestPathInCluster utilizing an A* algorithm. More...

void	AStar (const reco::ClusterHit3D , const reco::ClusterHit3D , float alpha, kdTree::KdTreeNode &, reco::ClusterParameters &) const
	Algorithm to find shortest path between two 3D hits. More...

void	ReconstructBestPath (const reco::ClusterHit3D *, BestNodeMap &, reco::HitPairListPtr &, reco::EdgeList &) const

float	DistanceBetweenNodes (const reco::ClusterHit3D , const reco::ClusterHit3D ) const

void	LeastCostPath (const reco::EdgeTuple &, const reco::ClusterHit3D *, reco::ClusterParameters &, float &) const
	Find the lowest cost path between two nodes using MST edges. More...

void	buildConvexHull (reco::ClusterParameters &, reco::HitPairListPtr &, int level=0) const

float	findConvexHullEndPoints (const reco::EdgeList &, const reco::ClusterHit3D , const reco::ClusterHit3D ) const

Private Attributes
bool	fEnableMonitoring
	Data members to follow. More...

size_t	fMinTinyClusterSize
	Minimum size for a "tiny" cluster. More...

float	fConvexHullKinkAngle
	Angle to declare a kink in convex hull calc. More...

float	fConvexHullMinSep
	Min hit separation to conisder in convex hull. More...

std::vector< float >	fTimeVector

geo::Geometry const *	fGeometry

PrincipalComponentsAlg	fPCAAlg

kdTree	fkdTree

std::unique_ptr< lar_cluster3d::IClusterParametersBuilder >	fClusterBuilder
	Common cluster builder tool. More...

Detailed Description

Definition at line 48 of file MSTPathFinder_tool.cc.

Member Typedef Documentation

using lar_cluster3d::MSTPathFinder::BestNodeMap = std::unordered_map<const reco::ClusterHit3D*, BestNodeTuple>

private

Definition at line 141 of file MSTPathFinder_tool.cc.

using lar_cluster3d::MSTPathFinder::BestNodeTuple = std::tuple<const reco::ClusterHit3D*, float, float>

private

Definition at line 140 of file MSTPathFinder_tool.cc.

using lar_cluster3d::MSTPathFinder::MinMaxPointPair = std::pair<MinMaxPoints, MinMaxPoints>

private

Definition at line 162 of file MSTPathFinder_tool.cc.

using lar_cluster3d::MSTPathFinder::MinMaxPoints = std::pair<reco::ProjectedPoint, reco::ProjectedPoint>

private

Add ability to build the convex hull (these needs to be split out! )

Definition at line 161 of file MSTPathFinder_tool.cc.

Member Enumeration Documentation

enum lar_cluster3d::MSTPathFinder::TimeValues

private

enumerate the possible values for time checking if monitoring timing

Enumerator
BUILDTHREEDHITS
BUILDHITTOHITMAP
RUNDBSCAN
BUILDCLUSTERINFO
PATHFINDING
NUMTIMEVALUES

Definition at line 93 of file MSTPathFinder_tool.cc.

                     {
       BUILDTHREEDHITS = 0,
       BUILDHITTOHITMAP = 1,
       RUNDBSCAN = 2,
       BUILDCLUSTERINFO = 3,
       PATHFINDING = 4,
       NUMTIMEVALUES
     };

Constructor & Destructor Documentation

lar_cluster3d::MSTPathFinder::MSTPathFinder ( const fhicl::ParameterSet & )

explicit

Constructor.

Parameters

pset

Definition at line 189 of file MSTPathFinder_tool.cc.

     : fPCAAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
     , fkdTree(pset.get<fhicl::ParameterSet>("kdTree"))
   {
     this->configure(pset);
   }

lar_cluster3d::MSTPathFinder::~MSTPathFinder ( )

Destructor.

Definition at line 198 of file MSTPathFinder_tool.cc.

198 {}

Member Function Documentation

void lar_cluster3d::MSTPathFinder::AStar	(	const reco::ClusterHit3D *	startNode,
		const reco::ClusterHit3D *	goalNode,
		float	alpha,
		kdTree::KdTreeNode &	topNode,
		reco::ClusterParameters &	clusterParams
	)		const

private

Algorithm to find shortest path between two 3D hits.

Definition at line 614 of file MSTPathFinder_tool.cc.

   {
     // Recover the list of hits and edges
     reco::HitPairListPtr& pathNodeList = clusterParams.getBestHitPairListPtr();
     reco::EdgeList& bestEdgeList = clusterParams.getBestEdgeList();
     reco::Hit3DToEdgeMap& curEdgeMap = clusterParams.getHit3DToEdgeMap();
 
     // Find the shortest path from start to goal using an A* algorithm
     // Keep track of the nodes which have already been evaluated
     reco::HitPairListPtr closedList;
 
     // Keep track of the nodes that have been "discovered" but yet to be evaluated
     reco::HitPairListPtr openList = {startNode};
 
     // Create a map which, for each node, will tell us the node it can be most efficiencly reached from.
     BestNodeMap bestNodeMap;
 
     bestNodeMap[startNode] =
       BestNodeTuple(startNode, 0., DistanceBetweenNodes(startNode, goalNode));
 
     alpha = 1.; //std::max(0.5,alpha);
 
     while (!openList.empty()) {
       // The list is not empty so by def we will return something
       reco::HitPairListPtr::iterator currentNodeItr = openList.begin();
 
       // If the list contains more than one element then we need to find the one with the smallest total estimated cost to the end
       if (openList.size() > 1)
         currentNodeItr = std::min_element(
           openList.begin(), openList.end(), [bestNodeMap](const auto& next, const auto& best) {
             return std::get<2>(bestNodeMap.at(next)) < std::get<2>(bestNodeMap.at(best));
           });
 
       // Easier to deal directly with the pointer to the node
       const reco::ClusterHit3D* currentNode = *currentNodeItr;
 
       // Check to see if we have reached the goal and need to evaluate the path
       if (currentNode == goalNode) {
         // The path reconstruction will
         ReconstructBestPath(goalNode, bestNodeMap, pathNodeList, bestEdgeList);
 
         break;
       }
 
       // Otherwise need to keep evaluating
       else {
         openList.erase(currentNodeItr);
         currentNode->setStatusBit(reco::ClusterHit3D::PATHCHECKED);
 
         // Get tuple values for the current node
         const BestNodeTuple& currentNodeTuple = bestNodeMap.at(currentNode);
         float currentNodeScore = std::get<1>(currentNodeTuple);
 
         // Recover the edges associated to the current point
         const reco::EdgeList& curEdgeList = curEdgeMap[currentNode];
 
         for (const auto& curEdge : curEdgeList) {
           const reco::ClusterHit3D* candHit3D = std::get<1>(curEdge);
 
           if (candHit3D->getStatusBits() & reco::ClusterHit3D::PATHCHECKED) continue;
 
           float tentative_gScore = currentNodeScore + std::get<2>(curEdge);
 
           // Have we seen the candidate node before?
           BestNodeMap::iterator candNodeItr = bestNodeMap.find(candHit3D);
 
           if (candNodeItr == bestNodeMap.end()) { openList.push_back(candHit3D); }
           else if (tentative_gScore > std::get<1>(candNodeItr->second))
             continue;
 
           // Make a guess at score to get to target...
           float guessToTarget = DistanceBetweenNodes(candHit3D, goalNode) / 0.3;
 
           bestNodeMap[candHit3D] =
             BestNodeTuple(currentNode, tentative_gScore, tentative_gScore + guessToTarget);
         }
       }
     }
 
     return;
   }

void lar_cluster3d::MSTPathFinder::buildConvexHull	(	reco::ClusterParameters &	clusterParameters,
		reco::HitPairListPtr &	hitPairListPtr,
		int	level = `0`
	)		const

private

Definition at line 929 of file MSTPathFinder_tool.cc.

   {
     // set an indention string
     std::string minuses(level / 2, '-');
     std::string indent(level / 2, ' ');
 
     indent += minuses;
 
     // The plan is to build the enclosing 2D polygon around the points in the PCA plane of most spread for this cluster
     // To do so we need to start by building a list of 2D projections onto the plane of most spread...
     reco::PrincipalComponents& pca = clusterParameters.getFullPCA();
 
     // Recover the parameters from the Principal Components Analysis that we need to project and accumulate
     const Eigen::Vector3f& pcaCenter = pca.getAvePosition();
     reco::ConvexHull& convexHull = clusterParameters.getConvexHull();
     reco::ProjectedPointList& pointList = convexHull.getProjectedPointList();
 
     // Loop through hits and do projection to plane
     for (const auto& hit3D : hitPairListPtr) {
       Eigen::Vector3f pcaToHitVec(hit3D->getPosition()[0] - pcaCenter(0),
                                   hit3D->getPosition()[1] - pcaCenter(1),
                                   hit3D->getPosition()[2] - pcaCenter(2));
       Eigen::Vector3f pcaToHit = pca.getEigenVectors() * pcaToHitVec;
 
       pointList.emplace_back(pcaToHit(1), pcaToHit(2), hit3D);
     }
 
     // Sort the point vec by increasing x, then increase y
     pointList.sort([](const auto& left, const auto& right) {
       return (std::abs(std::get<0>(left) - std::get<0>(right)) >
               std::numeric_limits<float>::epsilon()) ?
                std::get<0>(left) < std::get<0>(right) :
                std::get<1>(left) < std::get<1>(right);
     });
 
     // containers for finding the "best" hull...
     std::vector<ConvexHull> convexHullVec;
     std::vector<reco::ProjectedPointList> rejectedListVec;
     bool increaseDepth(pointList.size() > 3);
     float lastArea(std::numeric_limits<float>::max());
 
     while (increaseDepth) {
       // Get another convexHull container
       convexHullVec.push_back(ConvexHull(pointList, fConvexHullKinkAngle, fConvexHullMinSep));
       rejectedListVec.push_back(reco::ProjectedPointList());
 
       const ConvexHull& convexHull = convexHullVec.back();
       reco::ProjectedPointList& rejectedList = rejectedListVec.back();
       const reco::ProjectedPointList& convexHullPoints = convexHull.getConvexHull();
 
       increaseDepth = false;
 
       if (convexHull.getConvexHullArea() > 0.) {
         if (convexHullVec.size() < 2 || convexHull.getConvexHullArea() < 0.8 * lastArea) {
           for (auto& point : convexHullPoints) {
             pointList.remove(point);
             rejectedList.emplace_back(point);
           }
           lastArea = convexHull.getConvexHullArea();
           //                increaseDepth = true;
         }
       }
     }
 
     // do we have a valid convexHull?
     while (!convexHullVec.empty() && convexHullVec.back().getConvexHullArea() < 0.5) {
       convexHullVec.pop_back();
       rejectedListVec.pop_back();
     }
 
     // If we found the convex hull then build edges around the region
     if (!convexHullVec.empty()) {
       size_t nRejectedTotal(0);
       reco::HitPairListPtr locHitPairListPtr = hitPairListPtr;
 
       for (const auto& rejectedList : rejectedListVec) {
         nRejectedTotal += rejectedList.size();
 
         for (const auto& rejectedPoint : rejectedList) {
           if (convexHullVec.back().findNearestDistance(rejectedPoint) > 0.5)
             locHitPairListPtr.remove(std::get<2>(rejectedPoint));
         }
       }
 
       // Now add "edges" to the cluster to describe the convex hull (for the display)
       reco::ProjectedPointList& convexHullPointList = convexHull.getConvexHullPointList();
       reco::Hit3DToEdgeMap& edgeMap = convexHull.getConvexHullEdgeMap();
       reco::EdgeList& edgeList = convexHull.getConvexHullEdgeList();
 
       reco::ProjectedPoint lastPoint = convexHullVec.back().getConvexHull().front();
 
       for (auto& curPoint : convexHullVec.back().getConvexHull()) {
         if (curPoint == lastPoint) continue;
 
         const reco::ClusterHit3D* lastPoint3D = std::get<2>(lastPoint);
         const reco::ClusterHit3D* curPoint3D = std::get<2>(curPoint);
 
         float distBetweenPoints = (curPoint3D->getPosition()[0] - lastPoint3D->getPosition()[0]) *
                                     (curPoint3D->getPosition()[0] - lastPoint3D->getPosition()[0]) +
                                   (curPoint3D->getPosition()[1] - lastPoint3D->getPosition()[1]) *
                                     (curPoint3D->getPosition()[1] - lastPoint3D->getPosition()[1]) +
                                   (curPoint3D->getPosition()[2] - lastPoint3D->getPosition()[2]) *
                                     (curPoint3D->getPosition()[2] - lastPoint3D->getPosition()[2]);
 
         distBetweenPoints = std::sqrt(distBetweenPoints);
 
         reco::EdgeTuple edge(lastPoint3D, curPoint3D, distBetweenPoints);
 
         convexHullPointList.push_back(curPoint);
         edgeMap[lastPoint3D].push_back(edge);
         edgeMap[curPoint3D].push_back(edge);
         edgeList.emplace_back(edge);
 
         lastPoint = curPoint;
       }
 
       // Store the "extreme" points
       const ConvexHull::PointList& extremePoints = convexHullVec.back().getExtremePoints();
       reco::ProjectedPointList& extremePointList = convexHull.getConvexHullExtremePoints();
 
       for (const auto& point : extremePoints)
         extremePointList.push_back(point);
 
       // Store the "kink" points
       const reco::ConvexHullKinkTupleList& kinkPoints = convexHullVec.back().getKinkPoints();
       reco::ConvexHullKinkTupleList& kinkPointList = convexHull.getConvexHullKinkPoints();
 
       for (const auto& kink : kinkPoints)
         kinkPointList.push_back(kink);
     }
 
     return;
   }

void lar_cluster3d::MSTPathFinder::configure ( fhicl::ParameterSet const & pset )

override

Definition at line 203 of file MSTPathFinder_tool.cc.

   {
     fEnableMonitoring = pset.get<bool>("EnableMonitoring", true);
     fMinTinyClusterSize = pset.get<size_t>("MinTinyClusterSize", 40);
     fConvexHullKinkAngle = pset.get<float>("ConvexHullKinkAgle", 0.95);
     fConvexHullMinSep = pset.get<float>("ConvexHullMinSep", 0.65);
 
     art::ServiceHandle<geo::Geometry const> geometry;
 
     fGeometry = &*geometry;
 
     fTimeVector.resize(NUMTIMEVALUES, 0.);
 
     fClusterBuilder = art::make_tool<lar_cluster3d::IClusterParametersBuilder>(
       pset.get<fhicl::ParameterSet>("ClusterParamsBuilder"));
 
     return;
   }

reco::HitPairListPtr lar_cluster3d::MSTPathFinder::DepthFirstSearch	(	const reco::EdgeTuple &	curEdge,
		const reco::Hit3DToEdgeMap &	hitToEdgeMap,
		float &	bestTreeQuality
	)		const

private

a depth first search to find longest branches

Definition at line 804 of file MSTPathFinder_tool.cc.

   {
     reco::HitPairListPtr hitPairListPtr;
     float bestQuality(0.);
     float curEdgeWeight = std::max(0.3, std::get<2>(curEdge));
     float curEdgeProj(1. / curEdgeWeight);
 
     reco::Hit3DToEdgeMap::const_iterator edgeListItr = hitToEdgeMap.find(std::get<1>(curEdge));
 
     if (edgeListItr != hitToEdgeMap.end()) {
       // The input edge weight has quality factors applied, recalculate just the position difference
       const Eigen::Vector3f& firstHitPos = std::get<0>(curEdge)->getPosition();
       const Eigen::Vector3f& secondHitPos = std::get<1>(curEdge)->getPosition();
       float curEdgeVec[] = {secondHitPos[0] - firstHitPos[0],
                             secondHitPos[1] - firstHitPos[1],
                             secondHitPos[2] - firstHitPos[2]};
       float curEdgeMag = std::sqrt(curEdgeVec[0] * curEdgeVec[0] + curEdgeVec[1] * curEdgeVec[1] +
                                    curEdgeVec[2] * curEdgeVec[2]);
 
       curEdgeMag = std::max(float(0.1), curEdgeMag);
 
       for (const auto& edge : edgeListItr->second) {
         // skip the self reference
         if (std::get<1>(edge) == std::get<0>(curEdge)) continue;
 
         float quality(0.);
 
         reco::HitPairListPtr tempList = DepthFirstSearch(edge, hitToEdgeMap, quality);
 
         if (quality > bestQuality) {
           hitPairListPtr = tempList;
           bestQuality = quality;
           curEdgeProj = 1. / curEdgeMag;
         }
       }
     }
 
     hitPairListPtr.push_front(std::get<1>(curEdge));
 
     bestTreeQuality += bestQuality + curEdgeProj;
 
     return hitPairListPtr;
   }

float lar_cluster3d::MSTPathFinder::DistanceBetweenNodes	(	const reco::ClusterHit3D *	node1,
		const reco::ClusterHit3D *	node2
	)		const

private

Definition at line 772 of file MSTPathFinder_tool.cc.

   {
     const Eigen::Vector3f& node1Pos = node1->getPosition();
     const Eigen::Vector3f& node2Pos = node2->getPosition();
     float deltaNode[] = {
       node1Pos[0] - node2Pos[0], node1Pos[1] - node2Pos[1], node1Pos[2] - node2Pos[2]};
 
     // Standard euclidean distance
     return std::sqrt(deltaNode[0] * deltaNode[0] + deltaNode[1] * deltaNode[1] +
                      deltaNode[2] * deltaNode[2]);
 
     // Manhatten distance
     //return std::fabs(deltaNode[0]) + std::fabs(deltaNode[1]) + std::fabs(deltaNode[2]);
     /*
      // Chebyshev distance
      // We look for maximum distance by wires...
 
      // Now go through the hits and compare view by view to look for delta wire and tigher constraint on delta t
      int wireDeltas[] = {0,0,0};
 
      for(size_t idx = 0; idx < 3; idx++)
      wireDeltas[idx] = std::abs(int(node1->getWireIDs()[idx].Wire) - int(node2->getWireIDs()[idx].Wire));
 
      // put wire deltas in order...
      std::sort(wireDeltas, wireDeltas + 3);
 
      return std::sqrt(deltaNode[0]*deltaNode[0] + 0.09 * float(wireDeltas[2]*wireDeltas[2]));
      */
   }

void lar_cluster3d::MSTPathFinder::FindBestPathInCluster ( reco::ClusterParameters & curCluster ) const

private

Algorithm to find the best path through the given cluster.

Definition at line 443 of file MSTPathFinder_tool.cc.

   {
     reco::HitPairListPtr longestCluster;
     float bestQuality(0.);
     float aveNumEdges(0.);
     size_t maxNumEdges(0);
     size_t nIsolatedHits(0);
 
     // Now proceed with building the clusters
     cet::cpu_timer theClockPathFinding;
 
     // Start clocks if requested
     if (fEnableMonitoring) theClockPathFinding.start();
 
     reco::HitPairListPtr& hitPairList = curCluster.getHitPairListPtr();
     reco::Hit3DToEdgeMap& curEdgeMap = curCluster.getHit3DToEdgeMap();
     reco::EdgeList& bestEdgeList = curCluster.getBestEdgeList();
 
     // Do some spelunking...
     for (const auto& hit : hitPairList) {
       if (!curEdgeMap[hit].empty() && curEdgeMap[hit].size() == 1) {
         float quality(0.);
 
         reco::HitPairListPtr tempList =
           DepthFirstSearch(curEdgeMap[hit].front(), curEdgeMap, quality);
 
         tempList.push_front(std::get<0>(curEdgeMap[hit].front()));
 
         if (quality > bestQuality) {
           longestCluster = tempList;
           bestQuality = quality;
         }
 
         nIsolatedHits++;
       }
 
       aveNumEdges += float(curEdgeMap[hit].size());
       maxNumEdges = std::max(maxNumEdges, curEdgeMap[hit].size());
     }
 
     aveNumEdges /= float(hitPairList.size());
     std::cout << "----> # isolated hits: " << nIsolatedHits
               << ", longest branch: " << longestCluster.size()
               << ", cluster size: " << hitPairList.size() << ", ave # edges: " << aveNumEdges
               << ", max: " << maxNumEdges << std::endl;
 
     if (!longestCluster.empty()) {
       hitPairList = longestCluster;
       for (const auto& hit : hitPairList) {
         for (const auto& edge : curEdgeMap[hit])
           bestEdgeList.emplace_back(edge);
       }
 
       std::cout << "        ====> new cluster size: " << hitPairList.size() << std::endl;
     }
 
     if (fEnableMonitoring) {
       theClockPathFinding.stop();
 
       fTimeVector[PATHFINDING] += theClockPathFinding.accumulated_real_time();
     }
 
     return;
   }

void lar_cluster3d::MSTPathFinder::FindBestPathInCluster	(	reco::ClusterParameters &	clusterParams,
		kdTree::KdTreeNode &	topNode
	)		const

private

Alternative version of FindBestPathInCluster utilizing an A* algorithm.

Definition at line 509 of file MSTPathFinder_tool.cc.

   {
     // Set up for timing the function
     cet::cpu_timer theClockPathFinding;
 
     // Start clocks if requested
     if (fEnableMonitoring) theClockPathFinding.start();
 
     // Trial A* here
     if (clusterParams.getHitPairListPtr().size() > 2) {
       // Do a quick PCA to determine our parameter "alpha"
       fPCAAlg.PCAAnalysis_3D(clusterParams.getHitPairListPtr(), clusterParams.getFullPCA());
 
       // Recover the new fullPCA
       reco::PrincipalComponents& pca = clusterParams.getFullPCA();
 
       // The chances of a failure are remote, still we should check
       if (pca.getSvdOK()) {
         // Get references to what we need....
         reco::HitPairListPtr& curCluster = clusterParams.getHitPairListPtr();
         reco::Hit3DToEdgeMap& curEdgeMap = clusterParams.getHit3DToEdgeMap();
         float pcaLen = 3.0 * sqrt(pca.getEigenValues()[2]);
         float pcaWidth = 3.0 * sqrt(pca.getEigenValues()[1]);
         float pcaHeight = 3.0 * sqrt(pca.getEigenValues()[0]);
         const Eigen::Vector3f& pcaCenter = pca.getAvePosition();
         float alpha = std::min(float(1.), std::max(float(0.001), pcaWidth / pcaLen));
 
         // Create a temporary container for the isolated points
         reco::ProjectedPointList isolatedPointList;
 
         // Go through and find the isolated points, for those get the projection to the plane of maximum spread
         for (const auto& hit3D : curCluster) {
           // the definition of an isolated hit is that it only has one associated edge
           if (!curEdgeMap[hit3D].empty() && curEdgeMap[hit3D].size() == 1) {
             Eigen::Vector3f pcaToHitVec(hit3D->getPosition()[0] - pcaCenter(0),
                                         hit3D->getPosition()[1] - pcaCenter(1),
                                         hit3D->getPosition()[2] - pcaCenter(2));
             Eigen::Vector3f pcaToHit = pca.getEigenVectors() * pcaToHitVec;
 
             // This sets x,y where x is the longer spread, y the shorter
             isolatedPointList.emplace_back(pcaToHit(2), pcaToHit(1), hit3D);
           }
         }
 
         std::cout << "************* Finding best path with A* in cluster *****************"
                   << std::endl;
         std::cout << "**> There are " << curCluster.size() << " hits, " << isolatedPointList.size()
                   << " isolated hits, the alpha parameter is " << alpha << std::endl;
         std::cout << "**> PCA len: " << pcaLen << ", wid: " << pcaWidth << ", height: " << pcaHeight
                   << ", ratio: " << pcaHeight / pcaWidth << std::endl;
 
         // If no isolated points then nothing to do...
         if (isolatedPointList.size() > 1) {
           // Sort the point vec by increasing x, if same then by increasing y.
           isolatedPointList.sort([](const auto& left, const auto& right) {
             return (std::abs(std::get<0>(left) - std::get<0>(right)) >
                     std::numeric_limits<float>::epsilon()) ?
                      std::get<0>(left) < std::get<0>(right) :
                      std::get<1>(left) < std::get<1>(right);
           });
 
           // Ok, get the two most distance points...
           const reco::ClusterHit3D* startHit = std::get<2>(isolatedPointList.front());
           const reco::ClusterHit3D* stopHit = std::get<2>(isolatedPointList.back());
 
           std::cout << "**> Sorted " << isolatedPointList.size()
                     << " hits, longest distance: " << DistanceBetweenNodes(startHit, stopHit)
                     << std::endl;
 
           // Call the AStar function to try to find the best path...
           //                AStar(startHit,stopHit,alpha,topNode,clusterParams);
 
           float cost(std::numeric_limits<float>::max());
 
           LeastCostPath(curEdgeMap[startHit].front(), stopHit, clusterParams, cost);
 
           clusterParams.getBestHitPairListPtr().push_front(startHit);
 
           std::cout << "**> Best path has " << clusterParams.getBestHitPairListPtr().size()
                     << " hits, " << clusterParams.getBestEdgeList().size() << " edges" << std::endl;
         }
 
         // Recalculate the PCA based on the hits comprisig the path
         fPCAAlg.PCAAnalysis_3D(clusterParams.getBestHitPairListPtr(), pca);
 
         // And now compute the convex hull
         buildConvexHull(clusterParams, clusterParams.getBestHitPairListPtr());
       }
       else {
         std::cout << "++++++>>> PCA failure! # hits: " << clusterParams.getHitPairListPtr().size()
                   << std::endl;
       }
     }
 
     if (fEnableMonitoring) {
       theClockPathFinding.stop();
 
       fTimeVector[PATHFINDING] += theClockPathFinding.accumulated_real_time();
     }
 
     return;
   }

float lar_cluster3d::MSTPathFinder::findConvexHullEndPoints	(	const reco::EdgeList &	convexHull,
		const reco::ClusterHit3D *	first3D,
		const reco::ClusterHit3D *	last3D
	)		const

private

Definition at line 1066 of file MSTPathFinder_tool.cc.

   {
     float largestDistance(0.);
 
     // Note that edges are vectors and that the convex hull edge list will be ordered
     // The idea is that the maximum distance from a given edge is to the edge just before the edge that "turns back" towards the current edge
     // meaning that the dot product of the two edges becomes negative.
     reco::EdgeList::const_iterator firstEdgeItr = convexHull.begin();
 
     while (firstEdgeItr != convexHull.end()) {
       reco::EdgeList::const_iterator nextEdgeItr = firstEdgeItr;
 
       //        Eigen::Vector2f firstEdgeVec(std::get<3>(*firstEdgeItr),std::get<);
       //        Eigen::Vector2f lastPrimaryVec(lastPCA.getEigenVectors()[0][0],lastPCA.getEigenVectors()[0][1],lastPCA.getEigenVectors()[0][2]);
       //        float           cosToLast = newPrimaryVec.dot(lastPrimaryVec);
 
       while (++nextEdgeItr != convexHull.end()) {}
     }
 
     return largestDistance;
   }

float lar_cluster3d::MSTPathFinder::getTimeToExecute ( ) const

inlineoverridevirtual

If monitoring, recover the time to execute a particular function.

Implements lar_cluster3d::IClusterModAlg.

Definition at line 84 of file MSTPathFinder_tool.cc.

     {
       return std::accumulate(fTimeVector.begin(), fTimeVector.end(), 0.);
     }

void lar_cluster3d::MSTPathFinder::initializeHistograms ( art::TFileDirectory & histDir )

overridevirtual

Interface for initializing histograms if they are desired Note that the idea is to put hisgtograms in a subfolder.

Parameters

TFileDirectory - the folder to store the hists in

Implements lar_cluster3d::IClusterModAlg.

Definition at line 223 of file MSTPathFinder_tool.cc.

   {
     // It is assumed that the input TFileDirectory has been set up to group histograms into a common
     // folder at the calling routine's level. Here we create one more level of indirection to keep
     // histograms made by this tool separate.
     //    fFillHistograms = true;
     //
     //    std::string dirName = "ConvexHullPath";
     //
     //    art::TFileDirectory dir = histDir.mkdir(dirName.c_str());
     //
     //    // Divide into two sets of hists... those for the overall cluster and
     //    // those for the subclusters
     //    fTopNum3DHits = dir.make<TH1F>("TopNum3DHits",  "Number 3D Hits",  200,    0.,   200.);
 
     return;
   }

void lar_cluster3d::MSTPathFinder::LeastCostPath	(	const reco::EdgeTuple &	curEdge,
		const reco::ClusterHit3D *	goalNode,
		reco::ClusterParameters &	clusterParams,
		float &	showMeTheMoney
	)		const

private

Find the lowest cost path between two nodes using MST edges.

Definition at line 723 of file MSTPathFinder_tool.cc.

   {
     // Recover the mapping between hits and edges
     reco::Hit3DToEdgeMap& curEdgeMap = clusterParams.getHit3DToEdgeMap();
 
     reco::Hit3DToEdgeMap::const_iterator edgeListItr = curEdgeMap.find(std::get<1>(curEdge));
 
     showMeTheMoney = std::numeric_limits<float>::max();
 
     if (edgeListItr != curEdgeMap.end() && !edgeListItr->second.empty()) {
       reco::HitPairListPtr& bestNodeList = clusterParams.getBestHitPairListPtr();
       reco::EdgeList& bestEdgeList = clusterParams.getBestEdgeList();
 
       for (const auto& edge : edgeListItr->second) {
         // skip the self reference
         if (std::get<1>(edge) == std::get<0>(curEdge)) continue;
 
         // Have we found the droid we are looking for?
         if (std::get<1>(edge) == goalNode) {
           bestNodeList.push_back(goalNode);
           bestEdgeList.push_back(edge);
           showMeTheMoney = std::get<2>(edge);
           break;
         }
 
         // Keep searching, it is out there somewhere...
         float currentCost(0.);
 
         LeastCostPath(edge, goalNode, clusterParams, currentCost);
 
         if (currentCost < std::numeric_limits<float>::max()) {
           showMeTheMoney = std::get<2>(edge) + currentCost;
           break;
         }
       }
     }
 
     if (showMeTheMoney < std::numeric_limits<float>::max()) {
       clusterParams.getBestHitPairListPtr().push_front(std::get<1>(curEdge));
       clusterParams.getBestEdgeList().push_front(curEdge);
     }
 
     return;
   }

void lar_cluster3d::MSTPathFinder::ModifyClusters ( reco::ClusterParametersList & clusterParametersList ) const

overridevirtual

Scan an input collection of clusters and modify those according to the specific implementing algorithm.

Parameters

clusterParametersList A list of cluster objects (parameters from associated hits)

Top level interface tool for performing deghosting and primary path finding using a minimum spanning tree approach. This is a shell tool, it actually uses the Minimum Spanning Tree clusering tool...

Implements lar_cluster3d::IClusterModAlg.

Definition at line 242 of file MSTPathFinder_tool.cc.

   {
     /**
      *  @brief Top level interface tool for performing deghosting and primary path finding
      *  using a minimum spanning tree approach. This is a shell tool, it actually uses the
      *  Minimum Spanning Tree clusering tool...
      */
 
     // Initial clustering is done, now trim the list and get output parameters
     cet::cpu_timer theClockBuildClusters;
 
     // Start clocks if requested
     if (fEnableMonitoring) theClockBuildClusters.start();
 
     // Ok, the idea here is to loop over the input clusters and the process one at a time and then use the MST algorithm
     // to deghost and try to find the best path.
     for (auto& clusterParams : clusterParametersList) {
       // It turns out that computing the convex hull surrounding the points in the 2D projection onto the
       // plane of largest spread in the PCA is a good way to break up the cluster... and we do it here since
       // we (currently) want this to be part of the standard output
       buildConvexHull(clusterParams, clusterParams.getHitPairListPtr());
 
       // Make sure our cluster has enough hits...
       if (clusterParams.getHitPairListPtr().size() > fMinTinyClusterSize) {
         // DBScan is driven of its "epsilon neighborhood". Computing adjacency within DBScan can be time
         // consuming so the idea is the prebuild the adjaceny map and then run DBScan.
         // The following call does this work
         kdTree::KdTreeNodeList kdTreeNodeContainer;
         kdTree::KdTreeNode topNode =
           fkdTree.BuildKdTree(clusterParams.getHitPairListPtr(), kdTreeNodeContainer);
 
         if (fEnableMonitoring) fTimeVector.at(BUILDHITTOHITMAP) = fkdTree.getTimeToExecute();
 
         // We are making subclusters
         reco::ClusterParametersList& daughterParametersList = clusterParams.daughterList();
 
         // Run DBScan to get candidate clusters
         RunPrimsAlgorithm(clusterParams.getHitPairListPtr(), topNode, daughterParametersList);
 
         // Initial clustering is done, now trim the list and get output parameters
         cet::cpu_timer theClockBuildClusters;
 
         // Start clocks if requested
         if (fEnableMonitoring) theClockBuildClusters.start();
 
         fClusterBuilder->BuildClusterInfo(daughterParametersList);
 
         if (fEnableMonitoring) {
           theClockBuildClusters.stop();
 
           fTimeVector[BUILDCLUSTERINFO] = theClockBuildClusters.accumulated_real_time();
         }
 
         // Test run the path finding algorithm
         for (auto& daughterParams : daughterParametersList)
           FindBestPathInCluster(daughterParams, topNode);
       }
     }
 
     if (fEnableMonitoring) {
       theClockBuildClusters.stop();
 
       fTimeVector[BUILDCLUSTERINFO] = theClockBuildClusters.accumulated_real_time();
     }
 
     mf::LogDebug("MSTPathFinder") << ">>>>> Cluster Path finding done" << std::endl;
 
     return;
   }

void lar_cluster3d::MSTPathFinder::PruneAmbiguousHits	(	reco::ClusterParameters &	clusterParams,
		reco::Hit2DToClusterMap &	hit2DToClusterMap
	)		const

private

Prune the obvious ambiguous hits.

Definition at line 851 of file MSTPathFinder_tool.cc.

   {
 
     // Recover the HitPairListPtr from the input clusterParams (which will be the
     // only thing that has been provided)
     reco::HitPairListPtr& hitPairVector = clusterParams.getHitPairListPtr();
 
     size_t nStartedWith(hitPairVector.size());
     size_t nRejectedHits(0);
 
     reco::HitPairListPtr goodHits;
 
     // Loop through the hits and try to week out the clearly ambiguous ones
     for (const auto& hit3D : hitPairVector) {
       // Loop to try to remove ambiguous hits
       size_t n2DHitsIn3DHit(0);
       size_t nThisClusterOnly(0);
       size_t nOtherCluster(0);
 
       //        reco::ClusterParameters* otherCluster;
       const std::set<const reco::ClusterHit3D*>* otherClusterHits = 0;
 
       for (const auto& hit2D : hit3D->getHits()) {
         if (!hit2D) continue;
 
         n2DHitsIn3DHit++;
 
         if (hit2DToClusterMap[hit2D].size() < 2)
           nThisClusterOnly = hit2DToClusterMap[hit2D][&clusterParams].size();
         else {
           for (const auto& clusterHitMap : hit2DToClusterMap[hit2D]) {
             if (clusterHitMap.first == &clusterParams) continue;
 
             if (clusterHitMap.second.size() > nOtherCluster) {
               nOtherCluster = clusterHitMap.second.size();
               //                        otherCluster     = clusterHitMap.first;
               otherClusterHits = &clusterHitMap.second;
             }
           }
         }
       }
 
       if (n2DHitsIn3DHit < 3 && nThisClusterOnly > 1 && nOtherCluster > 0) {
         bool skip3DHit(false);
 
         for (const auto& otherHit3D : *otherClusterHits) {
           size_t nOther2DHits(0);
 
           for (const auto& otherHit2D : otherHit3D->getHits()) {
             if (!otherHit2D) continue;
 
             nOther2DHits++;
           }
 
           if (nOther2DHits > 2) {
             skip3DHit = true;
             nRejectedHits++;
             break;
           }
         }
 
         if (skip3DHit) continue;
       }
 
       goodHits.emplace_back(hit3D);
     }
 
     std::cout << "###>> Input " << nStartedWith << " hits, rejected: " << nRejectedHits
               << std::endl;
 
     hitPairVector.resize(goodHits.size());
     std::copy(goodHits.begin(), goodHits.end(), hitPairVector.begin());
 
     return;
   }

void lar_cluster3d::MSTPathFinder::ReconstructBestPath	(	const reco::ClusterHit3D *	goalNode,
		BestNodeMap &	bestNodeMap,
		reco::HitPairListPtr &	pathNodeList,
		reco::EdgeList &	bestEdgeList
	)		const

private

Definition at line 701 of file MSTPathFinder_tool.cc.

   {
     while (std::get<0>(bestNodeMap.at(goalNode)) != goalNode) {
       const reco::ClusterHit3D* nextNode = std::get<0>(bestNodeMap[goalNode]);
       reco::EdgeTuple bestEdge =
         reco::EdgeTuple(goalNode, nextNode, DistanceBetweenNodes(goalNode, nextNode));
 
       pathNodeList.push_front(goalNode);
       bestEdgeList.push_front(bestEdge);
 
       goalNode = nextNode;
     }
 
     pathNodeList.push_front(goalNode);
 
     return;
   }

void lar_cluster3d::MSTPathFinder::RunPrimsAlgorithm	(	const reco::HitPairListPtr &	hitPairList,
		kdTree::KdTreeNode &	topNode,
		reco::ClusterParametersList &	clusterParametersList
	)		const

private

Driver for Prim's algorithm.

Definition at line 314 of file MSTPathFinder_tool.cc.

   {
     // If no hits then no work
     if (hitPairList.empty()) return;
 
     // Now proceed with building the clusters
     cet::cpu_timer theClockDBScan;
 
     // Start clocks if requested
     if (fEnableMonitoring) theClockDBScan.start();
 
     // Initialization
     size_t clusterIdx(0);
 
     // This will contain our list of edges
     reco::EdgeList curEdgeList;
 
     // Get the first point
     reco::HitPairListPtr::const_iterator freeHitItr = hitPairList.begin();
     const reco::ClusterHit3D* lastAddedHit = *freeHitItr++;
 
     lastAddedHit->setStatusBit(reco::ClusterHit3D::CLUSTERATTACHED);
 
     // Make a cluster...
     clusterParametersList.push_back(reco::ClusterParameters());
 
     // Get an iterator to the first cluster
     reco::ClusterParameters* curCluster = &clusterParametersList.back();
 
     // We use pointers here because the objects they point to will change in the loop below
     reco::Hit3DToEdgeMap* curEdgeMap = &curCluster->getHit3DToEdgeMap();
     reco::HitPairListPtr* curClusterHitList = &curCluster->getHitPairListPtr();
 
     // Loop until all hits have been associated to a cluster
     while (1) {
       // and the 3D hit status bits
       lastAddedHit->setStatusBit(reco::ClusterHit3D::CLUSTERATTACHED);
 
       // Purge the current list to get rid of edges which point to hits already in the cluster
       for (reco::EdgeList::iterator curEdgeItr = curEdgeList.begin();
            curEdgeItr != curEdgeList.end();) {
         if (std::get<1>(*curEdgeItr)->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED)
           curEdgeItr = curEdgeList.erase(curEdgeItr);
         else
           curEdgeItr++;
       }
 
       // Add the lastUsedHit to the current cluster
       curClusterHitList->push_back(lastAddedHit);
 
       // Set up to find the list of nearest neighbors to the last used hit...
       kdTree::CandPairList CandPairList;
       float bestDistance(1.5); //std::numeric_limits<float>::max());
 
       // And find them... result will be an unordered list of neigbors
       fkdTree.FindNearestNeighbors(lastAddedHit, topNode, CandPairList, bestDistance);
 
       // Copy edges to the current list (but only for hits not already in a cluster)
       //        for(auto& pair : CandPairList)
       //            if (!(pair.second->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED)) curEdgeList.push_back(reco::EdgeTuple(lastAddedHit,pair.second,pair.first));
       for (auto& pair : CandPairList) {
         if (!(pair.second->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED)) {
           double edgeWeight =
             pair.first * lastAddedHit->getHitChiSquare() * pair.second->getHitChiSquare();
 
           curEdgeList.push_back(reco::EdgeTuple(lastAddedHit, pair.second, edgeWeight));
         }
       }
 
       // If the edge list is empty then we have a complete cluster
       if (curEdgeList.empty()) {
         std::cout << "-----------------------------------------------------------------------------"
                      "------------"
                   << std::endl;
         std::cout << "**> Cluster idx: " << clusterIdx++ << " has " << curClusterHitList->size()
                   << " hits" << std::endl;
 
         // Look for the next "free" hit
         freeHitItr = std::find_if(freeHitItr, hitPairList.end(), [](const auto& hit) {
           return !(hit->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED);
         });
 
         // If at end of input list we are done with all hits
         if (freeHitItr == hitPairList.end()) break;
 
         std::cout << "##################################################################>"
                      "Processing another cluster"
                   << std::endl;
 
         // Otherwise, get a new cluster and set up
         clusterParametersList.push_back(reco::ClusterParameters());
 
         curCluster = &clusterParametersList.back();
 
         curEdgeMap = &curCluster->getHit3DToEdgeMap();
         curClusterHitList = &curCluster->getHitPairListPtr();
         lastAddedHit = *freeHitItr++;
       }
       // Otherwise we are still processing the current cluster
       else {
         // Sort the list of edges by distance
         curEdgeList.sort([](const auto& left, const auto& right) {
           return std::get<2>(left) < std::get<2>(right);
         });
 
         // Populate the map with the edges...
         reco::EdgeTuple& curEdge = curEdgeList.front();
 
         (*curEdgeMap)[std::get<0>(curEdge)].push_back(curEdge);
         (*curEdgeMap)[std::get<1>(curEdge)].push_back(
           reco::EdgeTuple(std::get<1>(curEdge), std::get<0>(curEdge), std::get<2>(curEdge)));
 
         // Update the last hit to be added to the collection
         lastAddedHit = std::get<1>(curEdge);
       }
     }
 
     if (fEnableMonitoring) {
       theClockDBScan.stop();
 
       fTimeVector[RUNDBSCAN] = theClockDBScan.accumulated_real_time();
     }
 
     return;
   }

Member Data Documentation

std::unique_ptr<lar_cluster3d::IClusterParametersBuilder> lar_cluster3d::MSTPathFinder::fClusterBuilder

private

Common cluster builder tool.

Definition at line 186 of file MSTPathFinder_tool.cc.

float lar_cluster3d::MSTPathFinder::fConvexHullKinkAngle

private

Angle to declare a kink in convex hull calc.

Definition at line 175 of file MSTPathFinder_tool.cc.

float lar_cluster3d::MSTPathFinder::fConvexHullMinSep

private

Min hit separation to conisder in convex hull.

Definition at line 176 of file MSTPathFinder_tool.cc.

bool lar_cluster3d::MSTPathFinder::fEnableMonitoring

private

Data members to follow.

Definition at line 173 of file MSTPathFinder_tool.cc.

geo::Geometry const* lar_cluster3d::MSTPathFinder::fGeometry

private

Definition at line 180 of file MSTPathFinder_tool.cc.

kdTree lar_cluster3d::MSTPathFinder::fkdTree

private

Definition at line 183 of file MSTPathFinder_tool.cc.

size_t lar_cluster3d::MSTPathFinder::fMinTinyClusterSize

private

Minimum size for a "tiny" cluster.

Definition at line 174 of file MSTPathFinder_tool.cc.

PrincipalComponentsAlg lar_cluster3d::MSTPathFinder::fPCAAlg

private

Definition at line 182 of file MSTPathFinder_tool.cc.

std::vector<float> lar_cluster3d::MSTPathFinder::fTimeVector

mutableprivate

Definition at line 178 of file MSTPathFinder_tool.cc.

The documentation for this class was generated from the following file:

larreco/larreco/RecoAlg/Cluster3DAlgs/PathFinding/MSTPathFinder_tool.cc

Public Member Functions

Private Types

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation