LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
lar_cluster3d::MSTPathFinder Class Referenceabstract
Inheritance diagram for lar_cluster3d::MSTPathFinder:
lar_cluster3d::IClusterModAlg

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...
 
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 *, 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
 

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::IClusterParametersBuilderfClusterBuilder
 Common cluster builder tool. More...
 

Detailed Description

Definition at line 50 of file MSTPathFinder_tool.cc.

Member Typedef Documentation

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

Definition at line 140 of file MSTPathFinder_tool.cc.

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

Definition at line 139 of file MSTPathFinder_tool.cc.

Definition at line 161 of file MSTPathFinder_tool.cc.

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

Definition at line 160 of file MSTPathFinder_tool.cc.

Member Enumeration Documentation

enumerate the possible values for time checking if monitoring timing

Enumerator
BUILDTHREEDHITS 
BUILDHITTOHITMAP 
RUNDBSCAN 
BUILDCLUSTERINFO 
PATHFINDING 
NUMTIMEVALUES 

Definition at line 94 of file MSTPathFinder_tool.cc.

Constructor & Destructor Documentation

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

Constructor.

Parameters
pset

Definition at line 184 of file MSTPathFinder_tool.cc.

References configure().

185  : fPCAAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
186  , fkdTree(pset.get<fhicl::ParameterSet>("kdTree"))
187  {
188  this->configure(pset);
189  }
PrincipalComponentsAlg fPCAAlg
void configure(fhicl::ParameterSet const &pset) override
lar_cluster3d::MSTPathFinder::~MSTPathFinder ( )

Destructor.

Definition at line 193 of file MSTPathFinder_tool.cc.

193 {}

Member Function Documentation

void lar_cluster3d::MSTPathFinder::AStar ( const reco::ClusterHit3D startNode,
const reco::ClusterHit3D goalNode,
reco::ClusterParameters clusterParams 
) const
private

Algorithm to find shortest path between two 3D hits.

Definition at line 602 of file MSTPathFinder_tool.cc.

References DistanceBetweenNodes(), reco::ClusterParameters::getBestEdgeList(), reco::ClusterParameters::getBestHitPairListPtr(), reco::ClusterParameters::getHit3DToEdgeMap(), reco::ClusterHit3D::getStatusBits(), reco::ClusterHit3D::PATHCHECKED, ReconstructBestPath(), and reco::ClusterHit3D::setStatusBit().

605  {
606  // Recover the list of hits and edges
607  reco::HitPairListPtr& pathNodeList = clusterParams.getBestHitPairListPtr();
608  reco::EdgeList& bestEdgeList = clusterParams.getBestEdgeList();
609  reco::Hit3DToEdgeMap& curEdgeMap = clusterParams.getHit3DToEdgeMap();
610 
611  // Find the shortest path from start to goal using an A* algorithm
612  // Keep track of the nodes which have already been evaluated
613  reco::HitPairListPtr closedList;
614 
615  // Keep track of the nodes that have been "discovered" but yet to be evaluated
616  reco::HitPairListPtr openList = {startNode};
617 
618  // Create a map which, for each node, will tell us the node it can be most efficiencly reached from.
619  BestNodeMap bestNodeMap;
620 
621  bestNodeMap[startNode] =
622  BestNodeTuple(startNode, 0., DistanceBetweenNodes(startNode, goalNode));
623 
624  while (!openList.empty()) {
625  // The list is not empty so by def we will return something
626  reco::HitPairListPtr::iterator currentNodeItr = openList.begin();
627 
628  // If the list contains more than one element then we need to find the one with the smallest total estimated cost to the end
629  if (openList.size() > 1)
630  currentNodeItr = std::min_element(
631  openList.begin(), openList.end(), [bestNodeMap](const auto& next, const auto& best) {
632  return std::get<2>(bestNodeMap.at(next)) < std::get<2>(bestNodeMap.at(best));
633  });
634 
635  // Easier to deal directly with the pointer to the node
636  const reco::ClusterHit3D* currentNode = *currentNodeItr;
637 
638  // Check to see if we have reached the goal and need to evaluate the path
639  if (currentNode == goalNode) {
640  // The path reconstruction will
641  ReconstructBestPath(goalNode, bestNodeMap, pathNodeList, bestEdgeList);
642 
643  break;
644  }
645 
646  // Otherwise need to keep evaluating
647  else {
648  openList.erase(currentNodeItr);
650 
651  // Get tuple values for the current node
652  const BestNodeTuple& currentNodeTuple = bestNodeMap.at(currentNode);
653  float currentNodeScore = std::get<1>(currentNodeTuple);
654 
655  // Recover the edges associated to the current point
656  const reco::EdgeList& curEdgeList = curEdgeMap[currentNode];
657 
658  for (const auto& curEdge : curEdgeList) {
659  const reco::ClusterHit3D* candHit3D = std::get<1>(curEdge);
660 
661  if (candHit3D->getStatusBits() & reco::ClusterHit3D::PATHCHECKED) continue;
662 
663  float tentative_gScore = currentNodeScore + std::get<2>(curEdge);
664 
665  // Have we seen the candidate node before?
666  BestNodeMap::iterator candNodeItr = bestNodeMap.find(candHit3D);
667 
668  if (candNodeItr == bestNodeMap.end()) { openList.push_back(candHit3D); }
669  else if (tentative_gScore > std::get<1>(candNodeItr->second))
670  continue;
671 
672  // Make a guess at score to get to target...
673  float guessToTarget = DistanceBetweenNodes(candHit3D, goalNode) / 0.3;
674 
675  bestNodeMap[candHit3D] =
676  BestNodeTuple(currentNode, tentative_gScore, tentative_gScore + guessToTarget);
677  }
678  }
679  }
680 
681  return;
682  }
reco::HitPairListPtr & getBestHitPairListPtr()
Definition: Cluster3D.h:467
intermediate_table::iterator iterator
std::unordered_map< const reco::ClusterHit3D *, BestNodeTuple > BestNodeMap
std::tuple< const reco::ClusterHit3D *, float, float > BestNodeTuple
reco::EdgeList & getBestEdgeList()
Definition: Cluster3D.h:468
reco::Hit3DToEdgeMap & getHit3DToEdgeMap()
Definition: Cluster3D.h:466
unsigned int getStatusBits() const
Definition: Cluster3D.h:154
std::list< EdgeTuple > EdgeList
Definition: Cluster3D.h:337
Path checking algorithm has seen this hit.
Definition: Cluster3D.h:108
void ReconstructBestPath(const reco::ClusterHit3D *, BestNodeMap &, reco::HitPairListPtr &, reco::EdgeList &) const
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
float DistanceBetweenNodes(const reco::ClusterHit3D *, const reco::ClusterHit3D *) const
std::unordered_map< const reco::ClusterHit3D *, reco::EdgeList > Hit3DToEdgeMap
Definition: Cluster3D.h:339
void setStatusBit(unsigned bits) const
Definition: Cluster3D.h:177
void lar_cluster3d::MSTPathFinder::buildConvexHull ( reco::ClusterParameters clusterParameters,
reco::HitPairListPtr hitPairListPtr,
int  level = 0 
) const
private

Definition at line 907 of file MSTPathFinder_tool.cc.

References util::abs(), DEFINE_ART_CLASS_TOOL, fConvexHullKinkAngle, fConvexHullMinSep, reco::PrincipalComponents::getAvePosition(), lar_cluster3d::ConvexHull::getConvexHull(), reco::ClusterParameters::getConvexHull(), lar_cluster3d::ConvexHull::getConvexHullArea(), reco::ConvexHull::getConvexHullEdgeList(), reco::ConvexHull::getConvexHullEdgeMap(), reco::ConvexHull::getConvexHullExtremePoints(), reco::ConvexHull::getConvexHullKinkPoints(), reco::ConvexHull::getConvexHullPointList(), reco::PrincipalComponents::getEigenVectors(), reco::ClusterParameters::getFullPCA(), reco::ClusterHit3D::getPosition(), reco::ConvexHull::getProjectedPointList(), art::detail::indent(), art::left(), and art::right().

Referenced by FindBestPathInCluster(), and ModifyClusters().

910  {
911  // set an indention string
912  std::string minuses(level / 2, '-');
913  std::string indent(level / 2, ' ');
914 
915  indent += minuses;
916 
917  // The plan is to build the enclosing 2D polygon around the points in the PCA plane of most spread for this cluster
918  // To do so we need to start by building a list of 2D projections onto the plane of most spread...
919  reco::PrincipalComponents& pca = clusterParameters.getFullPCA();
920 
921  // Recover the parameters from the Principal Components Analysis that we need to project and accumulate
922  const Eigen::Vector3f& pcaCenter = pca.getAvePosition();
923  reco::ConvexHull& convexHull = clusterParameters.getConvexHull();
924  reco::ProjectedPointList& pointList = convexHull.getProjectedPointList();
925 
926  // Loop through hits and do projection to plane
927  for (const auto& hit3D : hitPairListPtr) {
928  Eigen::Vector3f pcaToHitVec(hit3D->getPosition()[0] - pcaCenter(0),
929  hit3D->getPosition()[1] - pcaCenter(1),
930  hit3D->getPosition()[2] - pcaCenter(2));
931  Eigen::Vector3f pcaToHit = pca.getEigenVectors() * pcaToHitVec;
932 
933  pointList.emplace_back(pcaToHit(1), pcaToHit(2), hit3D);
934  }
935 
936  // Sort the point vec by increasing x, then increase y
937  pointList.sort([](const auto& left, const auto& right) {
938  return (std::abs(std::get<0>(left) - std::get<0>(right)) >
939  std::numeric_limits<float>::epsilon()) ?
940  std::get<0>(left) < std::get<0>(right) :
941  std::get<1>(left) < std::get<1>(right);
942  });
943 
944  // containers for finding the "best" hull...
945  std::vector<ConvexHull> convexHullVec;
946  std::vector<reco::ProjectedPointList> rejectedListVec;
947  bool increaseDepth(pointList.size() > 3);
948  float lastArea(std::numeric_limits<float>::max());
949 
950  while (increaseDepth) {
951  // Get another convexHull container
952  convexHullVec.push_back(ConvexHull(pointList, fConvexHullKinkAngle, fConvexHullMinSep));
953  rejectedListVec.push_back(reco::ProjectedPointList());
954 
955  const ConvexHull& convexHull = convexHullVec.back();
956  reco::ProjectedPointList& rejectedList = rejectedListVec.back();
957  const reco::ProjectedPointList& convexHullPoints = convexHull.getConvexHull();
958 
959  increaseDepth = false;
960 
961  if (convexHull.getConvexHullArea() > 0.) {
962  if (convexHullVec.size() < 2 || convexHull.getConvexHullArea() < 0.8 * lastArea) {
963  for (auto& point : convexHullPoints) {
964  pointList.remove(point);
965  rejectedList.emplace_back(point);
966  }
967  lastArea = convexHull.getConvexHullArea();
968  // increaseDepth = true;
969  }
970  }
971  }
972 
973  // do we have a valid convexHull?
974  while (!convexHullVec.empty() && convexHullVec.back().getConvexHullArea() < 0.5) {
975  convexHullVec.pop_back();
976  rejectedListVec.pop_back();
977  }
978 
979  // If we found the convex hull then build edges around the region
980  if (!convexHullVec.empty()) {
981  reco::HitPairListPtr locHitPairListPtr = hitPairListPtr;
982 
983  for (const auto& rejectedList : rejectedListVec) {
984  for (const auto& rejectedPoint : rejectedList) {
985  if (convexHullVec.back().findNearestDistance(rejectedPoint) > 0.5)
986  locHitPairListPtr.remove(std::get<2>(rejectedPoint));
987  }
988  }
989 
990  // Now add "edges" to the cluster to describe the convex hull (for the display)
991  reco::ProjectedPointList& convexHullPointList = convexHull.getConvexHullPointList();
992  reco::Hit3DToEdgeMap& edgeMap = convexHull.getConvexHullEdgeMap();
993  reco::EdgeList& edgeList = convexHull.getConvexHullEdgeList();
994 
995  reco::ProjectedPoint lastPoint = convexHullVec.back().getConvexHull().front();
996 
997  for (auto& curPoint : convexHullVec.back().getConvexHull()) {
998  if (curPoint == lastPoint) continue;
999 
1000  const reco::ClusterHit3D* lastPoint3D = std::get<2>(lastPoint);
1001  const reco::ClusterHit3D* curPoint3D = std::get<2>(curPoint);
1002 
1003  float distBetweenPoints = (curPoint3D->getPosition()[0] - lastPoint3D->getPosition()[0]) *
1004  (curPoint3D->getPosition()[0] - lastPoint3D->getPosition()[0]) +
1005  (curPoint3D->getPosition()[1] - lastPoint3D->getPosition()[1]) *
1006  (curPoint3D->getPosition()[1] - lastPoint3D->getPosition()[1]) +
1007  (curPoint3D->getPosition()[2] - lastPoint3D->getPosition()[2]) *
1008  (curPoint3D->getPosition()[2] - lastPoint3D->getPosition()[2]);
1009 
1010  distBetweenPoints = std::sqrt(distBetweenPoints);
1011 
1012  reco::EdgeTuple edge(lastPoint3D, curPoint3D, distBetweenPoints);
1013 
1014  convexHullPointList.push_back(curPoint);
1015  edgeMap[lastPoint3D].push_back(edge);
1016  edgeMap[curPoint3D].push_back(edge);
1017  edgeList.emplace_back(edge);
1018 
1019  lastPoint = curPoint;
1020  }
1021 
1022  // Store the "extreme" points
1023  const ConvexHull::PointList& extremePoints = convexHullVec.back().getExtremePoints();
1024  reco::ProjectedPointList& extremePointList = convexHull.getConvexHullExtremePoints();
1025 
1026  for (const auto& point : extremePoints)
1027  extremePointList.push_back(point);
1028 
1029  // Store the "kink" points
1030  const reco::ConvexHullKinkTupleList& kinkPoints = convexHullVec.back().getKinkPoints();
1031  reco::ConvexHullKinkTupleList& kinkPointList = convexHull.getConvexHullKinkPoints();
1032 
1033  for (const auto& kink : kinkPoints)
1034  kinkPointList.push_back(kink);
1035  }
1036 
1037  return;
1038  }
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
std::tuple< float, float, const reco::ClusterHit3D * > ProjectedPoint
Projected coordinates and pointer to hit.
Definition: Cluster3D.h:345
const Eigen::Vector3f getPosition() const
Definition: Cluster3D.h:155
std::list< ProjectedPoint > ProjectedPointList
Definition: Cluster3D.h:346
std::list< Point > PointList
The list of the projected points.
Definition: ConvexHull.h:33
constexpr auto abs(T v)
Returns the absolute value of the argument.
float fConvexHullMinSep
Min hit separation to conisder in convex hull.
std::list< EdgeTuple > EdgeList
Definition: Cluster3D.h:337
reco::PrincipalComponents & getFullPCA()
Definition: Cluster3D.h:464
std::tuple< const reco::ClusterHit3D *, const reco::ClusterHit3D *, double > EdgeTuple
Definition: Cluster3D.h:336
Define a container for working with the convex hull.
Definition: Cluster3D.h:354
std::string indent(std::size_t const i)
const Eigen::Vector3f & getAvePosition() const
Definition: Cluster3D.h:244
std::list< ConvexHullKinkTuple > ConvexHullKinkTupleList
Definition: Cluster3D.h:349
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
float fConvexHullKinkAngle
Angle to declare a kink in convex hull calc.
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
reco::ConvexHull & getConvexHull()
Definition: Cluster3D.h:469
std::unordered_map< const reco::ClusterHit3D *, reco::EdgeList > Hit3DToEdgeMap
Definition: Cluster3D.h:339
reco::ProjectedPointList & getProjectedPointList()
Definition: Cluster3D.h:370
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:243
virtual void lar_cluster3d::IClusterModAlg::configure ( const fhicl::ParameterSet )
pure virtualinherited

Interface for configuring the particular algorithm tool.

Parameters
ParameterSetThe input set of parameters for configuration
void lar_cluster3d::MSTPathFinder::configure ( fhicl::ParameterSet const &  pset)
override

Definition at line 197 of file MSTPathFinder_tool.cc.

References fClusterBuilder, fConvexHullKinkAngle, fConvexHullMinSep, fEnableMonitoring, fGeometry, fMinTinyClusterSize, fTimeVector, fhicl::ParameterSet::get(), and NUMTIMEVALUES.

Referenced by MSTPathFinder().

198  {
199  fEnableMonitoring = pset.get<bool>("EnableMonitoring", true);
200  fMinTinyClusterSize = pset.get<size_t>("MinTinyClusterSize", 40);
201  fConvexHullKinkAngle = pset.get<float>("ConvexHullKinkAgle", 0.95);
202  fConvexHullMinSep = pset.get<float>("ConvexHullMinSep", 0.65);
203 
205 
206  fGeometry = &*geometry;
207 
208  fTimeVector.resize(NUMTIMEVALUES, 0.);
209 
210  fClusterBuilder = art::make_tool<lar_cluster3d::IClusterParametersBuilder>(
211  pset.get<fhicl::ParameterSet>("ClusterParamsBuilder"));
212 
213  return;
214  }
size_t fMinTinyClusterSize
Minimum size for a "tiny" cluster.
float fConvexHullMinSep
Min hit separation to conisder in convex hull.
std::unique_ptr< lar_cluster3d::IClusterParametersBuilder > fClusterBuilder
Common cluster builder tool.
float fConvexHullKinkAngle
Angle to declare a kink in convex hull calc.
geo::Geometry const * fGeometry
bool fEnableMonitoring
Data members to follow.
std::vector< float > fTimeVector
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 784 of file MSTPathFinder_tool.cc.

Referenced by FindBestPathInCluster().

787  {
788  reco::HitPairListPtr hitPairListPtr;
789  float bestQuality(0.);
790  float curEdgeWeight = std::max(0.3, std::get<2>(curEdge));
791  float curEdgeProj(1. / curEdgeWeight);
792 
793  reco::Hit3DToEdgeMap::const_iterator edgeListItr = hitToEdgeMap.find(std::get<1>(curEdge));
794 
795  if (edgeListItr != hitToEdgeMap.end()) {
796  // The input edge weight has quality factors applied, recalculate just the position difference
797  const Eigen::Vector3f& firstHitPos = std::get<0>(curEdge)->getPosition();
798  const Eigen::Vector3f& secondHitPos = std::get<1>(curEdge)->getPosition();
799  float curEdgeVec[] = {secondHitPos[0] - firstHitPos[0],
800  secondHitPos[1] - firstHitPos[1],
801  secondHitPos[2] - firstHitPos[2]};
802  float curEdgeMag = std::sqrt(curEdgeVec[0] * curEdgeVec[0] + curEdgeVec[1] * curEdgeVec[1] +
803  curEdgeVec[2] * curEdgeVec[2]);
804 
805  curEdgeMag = std::max(float(0.1), curEdgeMag);
806 
807  for (const auto& edge : edgeListItr->second) {
808  // skip the self reference
809  if (std::get<1>(edge) == std::get<0>(curEdge)) continue;
810 
811  float quality(0.);
812 
813  reco::HitPairListPtr tempList = DepthFirstSearch(edge, hitToEdgeMap, quality);
814 
815  if (quality > bestQuality) {
816  hitPairListPtr = tempList;
817  bestQuality = quality;
818  curEdgeProj = 1. / curEdgeMag;
819  }
820  }
821  }
822 
823  hitPairListPtr.push_front(std::get<1>(curEdge));
824 
825  bestTreeQuality += bestQuality + curEdgeProj;
826 
827  return hitPairListPtr;
828  }
reco::HitPairListPtr DepthFirstSearch(const reco::EdgeTuple &, const reco::Hit3DToEdgeMap &, float &) const
a depth first search to find longest branches
intermediate_table::const_iterator const_iterator
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
float lar_cluster3d::MSTPathFinder::DistanceBetweenNodes ( const reco::ClusterHit3D node1,
const reco::ClusterHit3D node2 
) const
private

Definition at line 753 of file MSTPathFinder_tool.cc.

References reco::ClusterHit3D::getPosition().

Referenced by AStar(), FindBestPathInCluster(), and ReconstructBestPath().

755  {
756  const Eigen::Vector3f& node1Pos = node1->getPosition();
757  const Eigen::Vector3f& node2Pos = node2->getPosition();
758  float deltaNode[] = {
759  node1Pos[0] - node2Pos[0], node1Pos[1] - node2Pos[1], node1Pos[2] - node2Pos[2]};
760 
761  // Standard euclidean distance
762  return std::sqrt(deltaNode[0] * deltaNode[0] + deltaNode[1] * deltaNode[1] +
763  deltaNode[2] * deltaNode[2]);
764 
765  // Manhatten distance
766  //return std::fabs(deltaNode[0]) + std::fabs(deltaNode[1]) + std::fabs(deltaNode[2]);
767  /*
768  // Chebyshev distance
769  // We look for maximum distance by wires...
770 
771  // Now go through the hits and compare view by view to look for delta wire and tigher constraint on delta t
772  int wireDeltas[] = {0,0,0};
773 
774  for(size_t idx = 0; idx < 3; idx++)
775  wireDeltas[idx] = std::abs(int(node1->getWireIDs()[idx].Wire) - int(node2->getWireIDs()[idx].Wire));
776 
777  // put wire deltas in order...
778  std::sort(wireDeltas, wireDeltas + 3);
779 
780  return std::sqrt(deltaNode[0]*deltaNode[0] + 0.09 * float(wireDeltas[2]*wireDeltas[2]));
781  */
782  }
const Eigen::Vector3f getPosition() const
Definition: Cluster3D.h:155
void lar_cluster3d::MSTPathFinder::FindBestPathInCluster ( reco::ClusterParameters curCluster) const
private

Algorithm to find the best path through the given cluster.

Definition at line 433 of file MSTPathFinder_tool.cc.

References DepthFirstSearch(), util::empty(), fEnableMonitoring, fTimeVector, reco::ClusterParameters::getBestEdgeList(), reco::ClusterParameters::getHit3DToEdgeMap(), reco::ClusterParameters::getHitPairListPtr(), PATHFINDING, and util::size().

Referenced by ModifyClusters().

434  {
435  reco::HitPairListPtr longestCluster;
436  float bestQuality(0.);
437  float aveNumEdges(0.);
438  size_t maxNumEdges(0);
439  size_t nIsolatedHits(0);
440 
441  // Now proceed with building the clusters
442  cet::cpu_timer theClockPathFinding;
443 
444  // Start clocks if requested
445  if (fEnableMonitoring) theClockPathFinding.start();
446 
447  reco::HitPairListPtr& hitPairList = curCluster.getHitPairListPtr();
448  reco::Hit3DToEdgeMap& curEdgeMap = curCluster.getHit3DToEdgeMap();
449  reco::EdgeList& bestEdgeList = curCluster.getBestEdgeList();
450 
451  // Do some spelunking...
452  for (const auto& hit : hitPairList) {
453  if (!curEdgeMap[hit].empty() && curEdgeMap[hit].size() == 1) {
454  float quality(0.);
455 
456  reco::HitPairListPtr tempList =
457  DepthFirstSearch(curEdgeMap[hit].front(), curEdgeMap, quality);
458 
459  tempList.push_front(std::get<0>(curEdgeMap[hit].front()));
460 
461  if (quality > bestQuality) {
462  longestCluster = tempList;
463  bestQuality = quality;
464  }
465 
466  nIsolatedHits++;
467  }
468 
469  aveNumEdges += float(curEdgeMap[hit].size());
470  maxNumEdges = std::max(maxNumEdges, curEdgeMap[hit].size());
471  }
472 
473  aveNumEdges /= float(hitPairList.size());
474  std::cout << "----> # isolated hits: " << nIsolatedHits
475  << ", longest branch: " << longestCluster.size()
476  << ", cluster size: " << hitPairList.size() << ", ave # edges: " << aveNumEdges
477  << ", max: " << maxNumEdges << std::endl;
478 
479  if (!longestCluster.empty()) {
480  hitPairList = longestCluster;
481  for (const auto& hit : hitPairList) {
482  for (const auto& edge : curEdgeMap[hit])
483  bestEdgeList.emplace_back(edge);
484  }
485 
486  std::cout << " ====> new cluster size: " << hitPairList.size() << std::endl;
487  }
488 
489  if (fEnableMonitoring) {
490  theClockPathFinding.stop();
491 
492  fTimeVector[PATHFINDING] += theClockPathFinding.accumulated_real_time();
493  }
494 
495  return;
496  }
reco::HitPairListPtr DepthFirstSearch(const reco::EdgeTuple &, const reco::Hit3DToEdgeMap &, float &) const
a depth first search to find longest branches
reco::EdgeList & getBestEdgeList()
Definition: Cluster3D.h:468
reco::Hit3DToEdgeMap & getHit3DToEdgeMap()
Definition: Cluster3D.h:466
reco::HitPairListPtr & getHitPairListPtr()
Definition: Cluster3D.h:463
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
std::list< EdgeTuple > EdgeList
Definition: Cluster3D.h:337
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
Detector simulation of raw signals on wires.
bool fEnableMonitoring
Data members to follow.
std::vector< float > fTimeVector
std::unordered_map< const reco::ClusterHit3D *, reco::EdgeList > Hit3DToEdgeMap
Definition: Cluster3D.h:339
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
void lar_cluster3d::MSTPathFinder::FindBestPathInCluster ( reco::ClusterParameters clusterParams,
kdTree::KdTreeNode  
) const
private

Alternative version of FindBestPathInCluster utilizing an A* algorithm.

Definition at line 498 of file MSTPathFinder_tool.cc.

References util::abs(), buildConvexHull(), DistanceBetweenNodes(), util::empty(), fEnableMonitoring, fPCAAlg, fTimeVector, reco::PrincipalComponents::getAvePosition(), reco::ClusterParameters::getBestEdgeList(), reco::ClusterParameters::getBestHitPairListPtr(), reco::PrincipalComponents::getEigenValues(), reco::PrincipalComponents::getEigenVectors(), reco::ClusterParameters::getFullPCA(), reco::ClusterParameters::getHit3DToEdgeMap(), reco::ClusterParameters::getHitPairListPtr(), reco::PrincipalComponents::getSvdOK(), LeastCostPath(), art::left(), PATHFINDING, lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_3D(), art::right(), and util::size().

500  {
501  // Set up for timing the function
502  cet::cpu_timer theClockPathFinding;
503 
504  // Start clocks if requested
505  if (fEnableMonitoring) theClockPathFinding.start();
506 
507  // Trial A* here
508  if (clusterParams.getHitPairListPtr().size() > 2) {
509  // Do a quick PCA to determine our parameter "alpha"
510  fPCAAlg.PCAAnalysis_3D(clusterParams.getHitPairListPtr(), clusterParams.getFullPCA());
511 
512  // Recover the new fullPCA
513  reco::PrincipalComponents& pca = clusterParams.getFullPCA();
514 
515  // The chances of a failure are remote, still we should check
516  if (pca.getSvdOK()) {
517  // Get references to what we need....
518  reco::HitPairListPtr& curCluster = clusterParams.getHitPairListPtr();
519  reco::Hit3DToEdgeMap& curEdgeMap = clusterParams.getHit3DToEdgeMap();
520  float pcaLen = 3.0 * sqrt(pca.getEigenValues()[2]);
521  float pcaWidth = 3.0 * sqrt(pca.getEigenValues()[1]);
522  float pcaHeight = 3.0 * sqrt(pca.getEigenValues()[0]);
523  const Eigen::Vector3f& pcaCenter = pca.getAvePosition();
524  float alpha = std::min(float(1.), std::max(float(0.001), pcaWidth / pcaLen));
525 
526  // Create a temporary container for the isolated points
527  reco::ProjectedPointList isolatedPointList;
528 
529  // Go through and find the isolated points, for those get the projection to the plane of maximum spread
530  for (const auto& hit3D : curCluster) {
531  // the definition of an isolated hit is that it only has one associated edge
532  if (!curEdgeMap[hit3D].empty() && curEdgeMap[hit3D].size() == 1) {
533  Eigen::Vector3f pcaToHitVec(hit3D->getPosition()[0] - pcaCenter(0),
534  hit3D->getPosition()[1] - pcaCenter(1),
535  hit3D->getPosition()[2] - pcaCenter(2));
536  Eigen::Vector3f pcaToHit = pca.getEigenVectors() * pcaToHitVec;
537 
538  // This sets x,y where x is the longer spread, y the shorter
539  isolatedPointList.emplace_back(pcaToHit(2), pcaToHit(1), hit3D);
540  }
541  }
542 
543  std::cout << "************* Finding best path with A* in cluster *****************"
544  << std::endl;
545  std::cout << "**> There are " << curCluster.size() << " hits, " << isolatedPointList.size()
546  << " isolated hits, the alpha parameter is " << alpha << std::endl;
547  std::cout << "**> PCA len: " << pcaLen << ", wid: " << pcaWidth << ", height: " << pcaHeight
548  << ", ratio: " << pcaHeight / pcaWidth << std::endl;
549 
550  // If no isolated points then nothing to do...
551  if (isolatedPointList.size() > 1) {
552  // Sort the point vec by increasing x, if same then by increasing y.
553  isolatedPointList.sort([](const auto& left, const auto& right) {
554  return (std::abs(std::get<0>(left) - std::get<0>(right)) >
555  std::numeric_limits<float>::epsilon()) ?
556  std::get<0>(left) < std::get<0>(right) :
557  std::get<1>(left) < std::get<1>(right);
558  });
559 
560  // Ok, get the two most distance points...
561  const reco::ClusterHit3D* startHit = std::get<2>(isolatedPointList.front());
562  const reco::ClusterHit3D* stopHit = std::get<2>(isolatedPointList.back());
563 
564  std::cout << "**> Sorted " << isolatedPointList.size()
565  << " hits, longest distance: " << DistanceBetweenNodes(startHit, stopHit)
566  << std::endl;
567 
568  // Call the AStar function to try to find the best path...
569  // AStar(startHit,stopHit,clusterParams);
570 
571  float cost(std::numeric_limits<float>::max());
572 
573  LeastCostPath(curEdgeMap[startHit].front(), stopHit, clusterParams, cost);
574 
575  clusterParams.getBestHitPairListPtr().push_front(startHit);
576 
577  std::cout << "**> Best path has " << clusterParams.getBestHitPairListPtr().size()
578  << " hits, " << clusterParams.getBestEdgeList().size() << " edges" << std::endl;
579  }
580 
581  // Recalculate the PCA based on the hits comprisig the path
582  fPCAAlg.PCAAnalysis_3D(clusterParams.getBestHitPairListPtr(), pca);
583 
584  // And now compute the convex hull
585  buildConvexHull(clusterParams, clusterParams.getBestHitPairListPtr());
586  }
587  else {
588  std::cout << "++++++>>> PCA failure! # hits: " << clusterParams.getHitPairListPtr().size()
589  << std::endl;
590  }
591  }
592 
593  if (fEnableMonitoring) {
594  theClockPathFinding.stop();
595 
596  fTimeVector[PATHFINDING] += theClockPathFinding.accumulated_real_time();
597  }
598 
599  return;
600  }
reco::HitPairListPtr & getBestHitPairListPtr()
Definition: Cluster3D.h:467
bool getSvdOK() const
Definition: Cluster3D.h:240
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
std::list< ProjectedPoint > ProjectedPointList
Definition: Cluster3D.h:346
constexpr auto abs(T v)
Returns the absolute value of the argument.
reco::EdgeList & getBestEdgeList()
Definition: Cluster3D.h:468
reco::Hit3DToEdgeMap & getHit3DToEdgeMap()
Definition: Cluster3D.h:466
reco::HitPairListPtr & getHitPairListPtr()
Definition: Cluster3D.h:463
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
const EigenValues & getEigenValues() const
Definition: Cluster3D.h:242
reco::PrincipalComponents & getFullPCA()
Definition: Cluster3D.h:464
void LeastCostPath(const reco::EdgeTuple &, const reco::ClusterHit3D *, reco::ClusterParameters &, float &) const
Find the lowest cost path between two nodes using MST edges.
const Eigen::Vector3f & getAvePosition() const
Definition: Cluster3D.h:244
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
PrincipalComponentsAlg fPCAAlg
float DistanceBetweenNodes(const reco::ClusterHit3D *, const reco::ClusterHit3D *) const
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
bool fEnableMonitoring
Data members to follow.
std::vector< float > fTimeVector
std::unordered_map< const reco::ClusterHit3D *, reco::EdgeList > Hit3DToEdgeMap
Definition: Cluster3D.h:339
void buildConvexHull(reco::ClusterParameters &, reco::HitPairListPtr &, int level=0) const
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:243
float lar_cluster3d::MSTPathFinder::getTimeToExecute ( ) const
inlineoverridevirtual

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

Implements lar_cluster3d::IClusterModAlg.

Definition at line 85 of file MSTPathFinder_tool.cc.

References fTimeVector.

86  {
87  return std::accumulate(fTimeVector.begin(), fTimeVector.end(), 0.);
88  }
std::vector< float > fTimeVector
void lar_cluster3d::MSTPathFinder::initializeHistograms ( art::TFileDirectory &  )
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 216 of file MSTPathFinder_tool.cc.

217  {
218  // It is assumed that the input TFileDirectory has been set up to group histograms into a common
219  // folder at the calling routine's level. Here we create one more level of indirection to keep
220  // histograms made by this tool separate.
221  // fFillHistograms = true;
222  //
223  // std::string dirName = "ConvexHullPath";
224  //
225  // art::TFileDirectory dir = histDir.mkdir(dirName.c_str());
226  //
227  // // Divide into two sets of hists... those for the overall cluster and
228  // // those for the subclusters
229  // fTopNum3DHits = dir.make<TH1F>("TopNum3DHits", "Number 3D Hits", 200, 0., 200.);
230 
231  return;
232  }
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 705 of file MSTPathFinder_tool.cc.

References reco::ClusterParameters::getBestEdgeList(), reco::ClusterParameters::getBestHitPairListPtr(), and reco::ClusterParameters::getHit3DToEdgeMap().

Referenced by FindBestPathInCluster().

709  {
710  // Recover the mapping between hits and edges
711  reco::Hit3DToEdgeMap& curEdgeMap = clusterParams.getHit3DToEdgeMap();
712 
713  reco::Hit3DToEdgeMap::const_iterator edgeListItr = curEdgeMap.find(std::get<1>(curEdge));
714 
715  showMeTheMoney = std::numeric_limits<float>::max();
716 
717  if (edgeListItr != curEdgeMap.end() && !edgeListItr->second.empty()) {
718  reco::HitPairListPtr& bestNodeList = clusterParams.getBestHitPairListPtr();
719  reco::EdgeList& bestEdgeList = clusterParams.getBestEdgeList();
720 
721  for (const auto& edge : edgeListItr->second) {
722  // skip the self reference
723  if (std::get<1>(edge) == std::get<0>(curEdge)) continue;
724 
725  // Have we found the droid we are looking for?
726  if (std::get<1>(edge) == goalNode) {
727  bestNodeList.push_back(goalNode);
728  bestEdgeList.push_back(edge);
729  showMeTheMoney = std::get<2>(edge);
730  break;
731  }
732 
733  // Keep searching, it is out there somewhere...
734  float currentCost(0.);
735 
736  LeastCostPath(edge, goalNode, clusterParams, currentCost);
737 
738  if (currentCost < std::numeric_limits<float>::max()) {
739  showMeTheMoney = std::get<2>(edge) + currentCost;
740  break;
741  }
742  }
743  }
744 
745  if (showMeTheMoney < std::numeric_limits<float>::max()) {
746  clusterParams.getBestHitPairListPtr().push_front(std::get<1>(curEdge));
747  clusterParams.getBestEdgeList().push_front(curEdge);
748  }
749 
750  return;
751  }
reco::HitPairListPtr & getBestHitPairListPtr()
Definition: Cluster3D.h:467
reco::EdgeList & getBestEdgeList()
Definition: Cluster3D.h:468
reco::Hit3DToEdgeMap & getHit3DToEdgeMap()
Definition: Cluster3D.h:466
intermediate_table::const_iterator const_iterator
std::list< EdgeTuple > EdgeList
Definition: Cluster3D.h:337
void LeastCostPath(const reco::EdgeTuple &, const reco::ClusterHit3D *, reco::ClusterParameters &, float &) const
Find the lowest cost path between two nodes using MST edges.
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
std::unordered_map< const reco::ClusterHit3D *, reco::EdgeList > Hit3DToEdgeMap
Definition: Cluster3D.h:339
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
clusterParametersListA 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 234 of file MSTPathFinder_tool.cc.

References BUILDCLUSTERINFO, buildConvexHull(), BUILDHITTOHITMAP, lar_cluster3d::kdTree::BuildKdTree(), fClusterBuilder, fEnableMonitoring, FindBestPathInCluster(), fkdTree, fMinTinyClusterSize, fTimeVector, lar_cluster3d::kdTree::getTimeToExecute(), and RunPrimsAlgorithm().

235  {
242  // Initial clustering is done, now trim the list and get output parameters
243  cet::cpu_timer theClockBuildClusters;
244 
245  // Start clocks if requested
246  if (fEnableMonitoring) theClockBuildClusters.start();
247 
248  // Ok, the idea here is to loop over the input clusters and the process one at a time and then use the MST algorithm
249  // to deghost and try to find the best path.
250  for (auto& clusterParams : clusterParametersList) {
251  // It turns out that computing the convex hull surrounding the points in the 2D projection onto the
252  // plane of largest spread in the PCA is a good way to break up the cluster... and we do it here since
253  // we (currently) want this to be part of the standard output
254  buildConvexHull(clusterParams, clusterParams.getHitPairListPtr());
255 
256  // Make sure our cluster has enough hits...
257  if (clusterParams.getHitPairListPtr().size() > fMinTinyClusterSize) {
258  // DBScan is driven of its "epsilon neighborhood". Computing adjacency within DBScan can be time
259  // consuming so the idea is the prebuild the adjaceny map and then run DBScan.
260  // The following call does this work
261  kdTree::KdTreeNodeList kdTreeNodeContainer;
262  kdTree::KdTreeNode topNode =
263  fkdTree.BuildKdTree(clusterParams.getHitPairListPtr(), kdTreeNodeContainer);
264 
266 
267  // We are making subclusters
268  reco::ClusterParametersList& daughterParametersList = clusterParams.daughterList();
269 
270  // Run DBScan to get candidate clusters
271  RunPrimsAlgorithm(clusterParams.getHitPairListPtr(), topNode, daughterParametersList);
272 
273  // Initial clustering is done, now trim the list and get output parameters
274  cet::cpu_timer theClockBuildClusters;
275 
276  // Start clocks if requested
277  if (fEnableMonitoring) theClockBuildClusters.start();
278 
279  fClusterBuilder->BuildClusterInfo(daughterParametersList);
280 
281  if (fEnableMonitoring) {
282  theClockBuildClusters.stop();
283 
284  fTimeVector[BUILDCLUSTERINFO] = theClockBuildClusters.accumulated_real_time();
285  }
286 
287  // Test run the path finding algorithm
288  for (auto& daughterParams : daughterParametersList)
289  FindBestPathInCluster(daughterParams, topNode);
290  }
291  }
292 
293  if (fEnableMonitoring) {
294  theClockBuildClusters.stop();
295 
296  fTimeVector[BUILDCLUSTERINFO] = theClockBuildClusters.accumulated_real_time();
297  }
298 
299  mf::LogDebug("MSTPathFinder") << ">>>>> Cluster Path finding done" << std::endl;
300 
301  return;
302  }
void RunPrimsAlgorithm(const reco::HitPairListPtr &, kdTree::KdTreeNode &, reco::ClusterParametersList &) const
Driver for Prim&#39;s algorithm.
size_t fMinTinyClusterSize
Minimum size for a "tiny" cluster.
void FindBestPathInCluster(reco::ClusterParameters &) const
Algorithm to find the best path through the given cluster.
std::list< KdTreeNode > KdTreeNodeList
Definition: kdTree.h:69
std::unique_ptr< lar_cluster3d::IClusterParametersBuilder > fClusterBuilder
Common cluster builder tool.
bool fEnableMonitoring
Data members to follow.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
std::vector< float > fTimeVector
float getTimeToExecute() const
Definition: kdTree.h:104
void buildConvexHull(reco::ClusterParameters &, reco::HitPairListPtr &, int level=0) const
KdTreeNode & BuildKdTree(Hit3DVec::iterator, Hit3DVec::iterator, KdTreeNodeList &, int depth=0) const
Given an input set of ClusterHit3D objects, build a kd tree structure.
Definition: kdTree.cxx:109
std::list< ClusterParameters > ClusterParametersList
Definition: Cluster3D.h:393
void lar_cluster3d::MSTPathFinder::PruneAmbiguousHits ( reco::ClusterParameters clusterParams,
reco::Hit2DToClusterMap hit2DToClusterMap 
) const
private

Prune the obvious ambiguous hits.

Definition at line 830 of file MSTPathFinder_tool.cc.

References reco::ClusterParameters::getHitPairListPtr(), and util::size().

832  {
833 
834  // Recover the HitPairListPtr from the input clusterParams (which will be the
835  // only thing that has been provided)
836  reco::HitPairListPtr& hitPairVector = clusterParams.getHitPairListPtr();
837 
838  size_t nStartedWith(hitPairVector.size());
839  size_t nRejectedHits(0);
840 
841  reco::HitPairListPtr goodHits;
842 
843  // Loop through the hits and try to week out the clearly ambiguous ones
844  for (const auto& hit3D : hitPairVector) {
845  // Loop to try to remove ambiguous hits
846  size_t n2DHitsIn3DHit(0);
847  size_t nThisClusterOnly(0);
848  size_t nOtherCluster(0);
849 
850  // reco::ClusterParameters* otherCluster;
851  const std::set<const reco::ClusterHit3D*>* otherClusterHits = 0;
852 
853  for (const auto& hit2D : hit3D->getHits()) {
854  if (!hit2D) continue;
855 
856  n2DHitsIn3DHit++;
857 
858  if (hit2DToClusterMap[hit2D].size() < 2)
859  nThisClusterOnly = hit2DToClusterMap[hit2D][&clusterParams].size();
860  else {
861  for (const auto& clusterHitMap : hit2DToClusterMap[hit2D]) {
862  if (clusterHitMap.first == &clusterParams) continue;
863 
864  if (clusterHitMap.second.size() > nOtherCluster) {
865  nOtherCluster = clusterHitMap.second.size();
866  // otherCluster = clusterHitMap.first;
867  otherClusterHits = &clusterHitMap.second;
868  }
869  }
870  }
871  }
872 
873  if (n2DHitsIn3DHit < 3 && nThisClusterOnly > 1 && nOtherCluster > 0) {
874  bool skip3DHit(false);
875 
876  for (const auto& otherHit3D : *otherClusterHits) {
877  size_t nOther2DHits(0);
878 
879  for (const auto& otherHit2D : otherHit3D->getHits()) {
880  if (!otherHit2D) continue;
881 
882  nOther2DHits++;
883  }
884 
885  if (nOther2DHits > 2) {
886  skip3DHit = true;
887  nRejectedHits++;
888  break;
889  }
890  }
891 
892  if (skip3DHit) continue;
893  }
894 
895  goodHits.emplace_back(hit3D);
896  }
897 
898  std::cout << "###>> Input " << nStartedWith << " hits, rejected: " << nRejectedHits
899  << std::endl;
900 
901  hitPairVector.resize(goodHits.size());
902  std::copy(goodHits.begin(), goodHits.end(), hitPairVector.begin());
903 
904  return;
905  }
reco::HitPairListPtr & getHitPairListPtr()
Definition: Cluster3D.h:463
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
void lar_cluster3d::MSTPathFinder::ReconstructBestPath ( const reco::ClusterHit3D goalNode,
BestNodeMap bestNodeMap,
reco::HitPairListPtr pathNodeList,
reco::EdgeList bestEdgeList 
) const
private

Definition at line 684 of file MSTPathFinder_tool.cc.

References DistanceBetweenNodes().

Referenced by AStar().

688  {
689  while (std::get<0>(bestNodeMap.at(goalNode)) != goalNode) {
690  const reco::ClusterHit3D* nextNode = std::get<0>(bestNodeMap[goalNode]);
691  reco::EdgeTuple bestEdge =
692  reco::EdgeTuple(goalNode, nextNode, DistanceBetweenNodes(goalNode, nextNode));
693 
694  pathNodeList.push_front(goalNode);
695  bestEdgeList.push_front(bestEdge);
696 
697  goalNode = nextNode;
698  }
699 
700  pathNodeList.push_front(goalNode);
701 
702  return;
703  }
std::tuple< const reco::ClusterHit3D *, const reco::ClusterHit3D *, double > EdgeTuple
Definition: Cluster3D.h:336
float DistanceBetweenNodes(const reco::ClusterHit3D *, const reco::ClusterHit3D *) const
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 305 of file MSTPathFinder_tool.cc.

References reco::ClusterHit3D::CLUSTERATTACHED, fEnableMonitoring, lar_cluster3d::kdTree::FindNearestNeighbors(), fkdTree, fTimeVector, reco::ClusterParameters::getHit3DToEdgeMap(), reco::ClusterHit3D::getHitChiSquare(), reco::ClusterParameters::getHitPairListPtr(), art::left(), art::right(), RUNDBSCAN, and reco::ClusterHit3D::setStatusBit().

Referenced by ModifyClusters().

308  {
309  // If no hits then no work
310  if (hitPairList.empty()) return;
311 
312  // Now proceed with building the clusters
313  cet::cpu_timer theClockDBScan;
314 
315  // Start clocks if requested
316  if (fEnableMonitoring) theClockDBScan.start();
317 
318  // Initialization
319  size_t clusterIdx(0);
320 
321  // This will contain our list of edges
322  reco::EdgeList curEdgeList;
323 
324  // Get the first point
325  reco::HitPairListPtr::const_iterator freeHitItr = hitPairList.begin();
326  const reco::ClusterHit3D* lastAddedHit = *freeHitItr++;
327 
329 
330  // Make a cluster...
331  clusterParametersList.push_back(reco::ClusterParameters());
332 
333  // Get an iterator to the first cluster
334  reco::ClusterParameters* curCluster = &clusterParametersList.back();
335 
336  // We use pointers here because the objects they point to will change in the loop below
337  reco::Hit3DToEdgeMap* curEdgeMap = &curCluster->getHit3DToEdgeMap();
338  reco::HitPairListPtr* curClusterHitList = &curCluster->getHitPairListPtr();
339 
340  // Loop until all hits have been associated to a cluster
341  while (1) {
342  // and the 3D hit status bits
344 
345  // Purge the current list to get rid of edges which point to hits already in the cluster
346  for (reco::EdgeList::iterator curEdgeItr = curEdgeList.begin();
347  curEdgeItr != curEdgeList.end();) {
348  if (std::get<1>(*curEdgeItr)->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED)
349  curEdgeItr = curEdgeList.erase(curEdgeItr);
350  else
351  curEdgeItr++;
352  }
353 
354  // Add the lastUsedHit to the current cluster
355  curClusterHitList->push_back(lastAddedHit);
356 
357  // Set up to find the list of nearest neighbors to the last used hit...
358  kdTree::CandPairList CandPairList;
359  float bestDistance(1.5); //std::numeric_limits<float>::max());
360 
361  // And find them... result will be an unordered list of neigbors
362  fkdTree.FindNearestNeighbors(lastAddedHit, topNode, CandPairList, bestDistance);
363 
364  // Copy edges to the current list (but only for hits not already in a cluster)
365  // for(auto& pair : CandPairList)
366  // if (!(pair.second->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED)) curEdgeList.push_back(reco::EdgeTuple(lastAddedHit,pair.second,pair.first));
367  for (auto& pair : CandPairList) {
368  if (!(pair.second->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED)) {
369  double edgeWeight =
370  pair.first * lastAddedHit->getHitChiSquare() * pair.second->getHitChiSquare();
371 
372  curEdgeList.push_back(reco::EdgeTuple(lastAddedHit, pair.second, edgeWeight));
373  }
374  }
375 
376  // If the edge list is empty then we have a complete cluster
377  if (curEdgeList.empty()) {
378  std::cout << "-----------------------------------------------------------------------------"
379  "------------"
380  << std::endl;
381  std::cout << "**> Cluster idx: " << clusterIdx++ << " has " << curClusterHitList->size()
382  << " hits" << std::endl;
383 
384  // Look for the next "free" hit
385  freeHitItr = std::find_if(freeHitItr, hitPairList.end(), [](const auto& hit) {
386  return !(hit->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED);
387  });
388 
389  // If at end of input list we are done with all hits
390  if (freeHitItr == hitPairList.end()) break;
391 
392  std::cout << "##################################################################>"
393  "Processing another cluster"
394  << std::endl;
395 
396  // Otherwise, get a new cluster and set up
397  clusterParametersList.push_back(reco::ClusterParameters());
398 
399  curCluster = &clusterParametersList.back();
400 
401  curEdgeMap = &curCluster->getHit3DToEdgeMap();
402  curClusterHitList = &curCluster->getHitPairListPtr();
403  lastAddedHit = *freeHitItr++;
404  }
405  // Otherwise we are still processing the current cluster
406  else {
407  // Sort the list of edges by distance
408  curEdgeList.sort([](const auto& left, const auto& right) {
409  return std::get<2>(left) < std::get<2>(right);
410  });
411 
412  // Populate the map with the edges...
413  reco::EdgeTuple& curEdge = curEdgeList.front();
414 
415  (*curEdgeMap)[std::get<0>(curEdge)].push_back(curEdge);
416  (*curEdgeMap)[std::get<1>(curEdge)].push_back(
417  reco::EdgeTuple(std::get<1>(curEdge), std::get<0>(curEdge), std::get<2>(curEdge)));
418 
419  // Update the last hit to be added to the collection
420  lastAddedHit = std::get<1>(curEdge);
421  }
422  }
423 
424  if (fEnableMonitoring) {
425  theClockDBScan.stop();
426 
427  fTimeVector[RUNDBSCAN] = theClockDBScan.accumulated_real_time();
428  }
429 
430  return;
431  }
intermediate_table::iterator iterator
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
size_t FindNearestNeighbors(const reco::ClusterHit3D *, const KdTreeNode &, CandPairList &, float &) const
Definition: kdTree.cxx:183
reco::Hit3DToEdgeMap & getHit3DToEdgeMap()
Definition: Cluster3D.h:466
intermediate_table::const_iterator const_iterator
reco::HitPairListPtr & getHitPairListPtr()
Definition: Cluster3D.h:463
std::list< EdgeTuple > EdgeList
Definition: Cluster3D.h:337
std::tuple< const reco::ClusterHit3D *, const reco::ClusterHit3D *, double > EdgeTuple
Definition: Cluster3D.h:336
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
Detector simulation of raw signals on wires.
std::list< CandPair > CandPairList
Definition: kdTree.h:84
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
bool fEnableMonitoring
Data members to follow.
std::vector< float > fTimeVector
float getHitChiSquare() const
Definition: Cluster3D.h:163
std::unordered_map< const reco::ClusterHit3D *, reco::EdgeList > Hit3DToEdgeMap
Definition: Cluster3D.h:339
attached to a cluster
Definition: Cluster3D.h:106
void setStatusBit(unsigned bits) const
Definition: Cluster3D.h:177

Member Data Documentation

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

Common cluster builder tool.

Definition at line 181 of file MSTPathFinder_tool.cc.

Referenced by configure(), and ModifyClusters().

float lar_cluster3d::MSTPathFinder::fConvexHullKinkAngle
private

Angle to declare a kink in convex hull calc.

Definition at line 170 of file MSTPathFinder_tool.cc.

Referenced by buildConvexHull(), and configure().

float lar_cluster3d::MSTPathFinder::fConvexHullMinSep
private

Min hit separation to conisder in convex hull.

Definition at line 171 of file MSTPathFinder_tool.cc.

Referenced by buildConvexHull(), and configure().

bool lar_cluster3d::MSTPathFinder::fEnableMonitoring
private

Data members to follow.

Definition at line 168 of file MSTPathFinder_tool.cc.

Referenced by configure(), FindBestPathInCluster(), ModifyClusters(), and RunPrimsAlgorithm().

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

Definition at line 175 of file MSTPathFinder_tool.cc.

Referenced by configure().

kdTree lar_cluster3d::MSTPathFinder::fkdTree
private

Definition at line 178 of file MSTPathFinder_tool.cc.

Referenced by ModifyClusters(), and RunPrimsAlgorithm().

size_t lar_cluster3d::MSTPathFinder::fMinTinyClusterSize
private

Minimum size for a "tiny" cluster.

Definition at line 169 of file MSTPathFinder_tool.cc.

Referenced by configure(), and ModifyClusters().

PrincipalComponentsAlg lar_cluster3d::MSTPathFinder::fPCAAlg
private

Definition at line 177 of file MSTPathFinder_tool.cc.

Referenced by FindBestPathInCluster().

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

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