LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
MSTPathFinder_tool.cc
Go to the documentation of this file.
1 
8 // Framework Includes
11 #include "art_root_io/TFileService.h"
12 #include "cetlib/cpu_timer.h"
13 #include "cetlib/search_path.h"
14 #include "fhiclcpp/ParameterSet.h"
16 
19 
20 // LArSoft includes
29 
30 // Eigen
31 #include <Eigen/Dense>
32 
33 // Root histograms
34 #include "TH1F.h"
35 #include "TH2F.h"
36 #include "TProfile.h"
37 
38 // std includes
39 #include <functional>
40 #include <iostream>
41 #include <memory>
42 #include <numeric> // std::accumulate
43 #include <string>
44 
45 //------------------------------------------------------------------------------------------------------------------------------------------
46 // implementation follows
47 
48 namespace lar_cluster3d {
49 
50  class MSTPathFinder : virtual public IClusterModAlg {
51  public:
57  explicit MSTPathFinder(const fhicl::ParameterSet&);
58 
63 
64  void configure(fhicl::ParameterSet const& pset) override;
65 
72  void initializeHistograms(art::TFileDirectory&) override;
73 
80  void ModifyClusters(reco::ClusterParametersList&) const override;
81 
85  float getTimeToExecute() const override
86  {
87  return std::accumulate(fTimeVector.begin(), fTimeVector.end(), 0.);
88  }
89 
90  private:
94  enum TimeValues {
97  RUNDBSCAN = 2,
101  };
102 
109 
114 
119 
124  const reco::Hit3DToEdgeMap&,
125  float&) const;
126 
131 
135  void AStar(const reco::ClusterHit3D*,
136  const reco::ClusterHit3D*,
137  reco::ClusterParameters&) const;
138 
139  using BestNodeTuple = std::tuple<const reco::ClusterHit3D*, float, float>;
140  using BestNodeMap = std::unordered_map<const reco::ClusterHit3D*, BestNodeTuple>;
141 
143  BestNodeMap&,
145  reco::EdgeList&) const;
146 
147  float DistanceBetweenNodes(const reco::ClusterHit3D*, const reco::ClusterHit3D*) const;
148 
152  void LeastCostPath(const reco::EdgeTuple&,
153  const reco::ClusterHit3D*,
155  float&) const;
156 
160  using MinMaxPoints = std::pair<reco::ProjectedPoint, reco::ProjectedPoint>;
161  using MinMaxPointPair = std::pair<MinMaxPoints, MinMaxPoints>;
162 
163  void buildConvexHull(reco::ClusterParameters&, reco::HitPairListPtr&, int level = 0) const;
164 
172 
173  mutable std::vector<float> fTimeVector;
174 
175  geo::Geometry const* fGeometry; //< pointer to the Geometry service
176 
177  PrincipalComponentsAlg fPCAAlg; // For running Principal Components Analysis
178  kdTree fkdTree; // For the kdTree
179 
180  std::unique_ptr<lar_cluster3d::IClusterParametersBuilder>
182  };
183 
185  : fPCAAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
186  , fkdTree(pset.get<fhicl::ParameterSet>("kdTree"))
187  {
188  this->configure(pset);
189  }
190 
191  //------------------------------------------------------------------------------------------------------------------------------------------
192 
194 
195  //------------------------------------------------------------------------------------------------------------------------------------------
196 
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  }
215 
216  void MSTPathFinder::initializeHistograms(art::TFileDirectory& /* histDir */)
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  }
233 
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  }
303 
304  //------------------------------------------------------------------------------------------------------------------------------------------
306  kdTree::KdTreeNode& topNode,
307  reco::ClusterParametersList& clusterParametersList) const
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  }
432 
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  }
497 
499  kdTree::KdTreeNode& /* topNode */) const
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  }
601 
603  const reco::ClusterHit3D* goalNode,
604  reco::ClusterParameters& clusterParams) const
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  }
683 
685  BestNodeMap& bestNodeMap,
686  reco::HitPairListPtr& pathNodeList,
687  reco::EdgeList& bestEdgeList) const
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  }
704 
706  const reco::ClusterHit3D* goalNode,
707  reco::ClusterParameters& clusterParams,
708  float& showMeTheMoney) const
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  }
752 
754  const reco::ClusterHit3D* node2) const
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  }
783 
785  const reco::Hit3DToEdgeMap& hitToEdgeMap,
786  float& bestTreeQuality) const
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  }
829 
831  reco::Hit2DToClusterMap& hit2DToClusterMap) const
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  }
906 
908  reco::HitPairListPtr& hitPairListPtr,
909  int level) const
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  }
1039 
1041 } // namespace lar_cluster3d
reco::HitPairListPtr & getBestHitPairListPtr()
Definition: Cluster3D.h:467
intermediate_table::iterator iterator
void RunPrimsAlgorithm(const reco::HitPairListPtr &, kdTree::KdTreeNode &, reco::ClusterParametersList &) const
Driver for Prim&#39;s algorithm.
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
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
std::unordered_map< const reco::ClusterHit3D *, BestNodeTuple > BestNodeMap
std::tuple< const reco::ClusterHit3D *, float, float > BestNodeTuple
size_t fMinTinyClusterSize
Minimum size for a "tiny" cluster.
kdTree class definiton
Definition: kdTree.h:31
void FindBestPathInCluster(reco::ClusterParameters &) const
Algorithm to find the best path through the given cluster.
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
void initializeHistograms(art::TFileDirectory &) override
Interface for initializing histograms if they are desired Note that the idea is to put hisgtograms in...
std::list< ProjectedPoint > ProjectedPointList
Definition: Cluster3D.h:346
Declaration of signal hit object.
reco::HitPairListPtr DepthFirstSearch(const reco::EdgeTuple &, const reco::Hit3DToEdgeMap &, float &) const
a depth first search to find longest branches
std::pair< MinMaxPoints, MinMaxPoints > MinMaxPointPair
std::list< Point > PointList
The list of the projected points.
Definition: ConvexHull.h:33
std::list< KdTreeNode > KdTreeNodeList
Definition: kdTree.h:69
constexpr auto abs(T v)
Returns the absolute value of the argument.
Implements a kdTree for use in clustering.
size_t FindNearestNeighbors(const reco::ClusterHit3D *, const KdTreeNode &, CandPairList &, float &) const
Definition: kdTree.cxx:183
reco::EdgeList & getBestEdgeList()
Definition: Cluster3D.h:468
reco::Hit3DToEdgeMap & getHit3DToEdgeMap()
Definition: Cluster3D.h:466
intermediate_table::const_iterator const_iterator
reco::HitPairListPtr & getHitPairListPtr()
Definition: Cluster3D.h:463
float fConvexHullMinSep
Min hit separation to conisder in convex hull.
MSTPathFinder(const fhicl::ParameterSet &)
Constructor.
unsigned int getStatusBits() const
Definition: Cluster3D.h:154
TimeValues
enumerate the possible values for time checking if monitoring timing
IClusterModAlg interface class definiton.
void PruneAmbiguousHits(reco::ClusterParameters &, reco::Hit2DToClusterMap &) const
Prune the obvious ambiguous hits.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
Implements a ConvexHull for use in clustering.
std::list< EdgeTuple > EdgeList
Definition: Cluster3D.h:337
define a kd tree node
Definition: kdTree.h:127
const EigenValues & getEigenValues() const
Definition: Cluster3D.h:242
reco::PrincipalComponents & getFullPCA()
Definition: Cluster3D.h:464
reco::Hit3DToEdgeMap & getConvexHullEdgeMap()
Definition: Cluster3D.h:372
void LeastCostPath(const reco::EdgeTuple &, const reco::ClusterHit3D *, reco::ClusterParameters &, float &) const
Find the lowest cost path between two nodes using MST edges.
std::tuple< const reco::ClusterHit3D *, const reco::ClusterHit3D *, double > EdgeTuple
Definition: Cluster3D.h:336
reco::ProjectedPointList & getConvexHullExtremePoints()
Definition: Cluster3D.h:374
parameter set interface
std::unordered_map< const reco::ClusterHit2D *, ClusterToHitPairSetMap > Hit2DToClusterMap
Definition: Cluster3D.h:499
reco::ConvexHullKinkTupleList & getConvexHullKinkPoints()
Definition: Cluster3D.h:375
Define a container for working with the convex hull.
Definition: Cluster3D.h:354
T get(std::string const &key) const
Definition: ParameterSet.h:314
std::unique_ptr< lar_cluster3d::IClusterParametersBuilder > fClusterBuilder
Common cluster builder tool.
std::string indent(std::size_t const i)
const Eigen::Vector3f & getAvePosition() const
Definition: Cluster3D.h:244
Path checking algorithm has seen this hit.
Definition: Cluster3D.h:108
void ReconstructBestPath(const reco::ClusterHit3D *, BestNodeMap &, reco::HitPairListPtr &, reco::EdgeList &) const
std::list< ConvexHullKinkTuple > ConvexHullKinkTupleList
Definition: Cluster3D.h:349
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
float getConvexHullArea() const
recover the area of the convex hull
Definition: ConvexHull.h:77
The geometry of one entire detector, as served by art.
Definition: Geometry.h:181
const PointList & getConvexHull() const
recover the list of convex hull vertices
Definition: ConvexHull.h:57
PrincipalComponentsAlg fPCAAlg
float DistanceBetweenNodes(const reco::ClusterHit3D *, const reco::ClusterHit3D *) const
float fConvexHullKinkAngle
Angle to declare a kink in convex hull calc.
Detector simulation of raw signals on wires.
std::list< CandPair > CandPairList
Definition: kdTree.h:84
ConvexHull class definiton.
Definition: ConvexHull.h:26
This header file defines the interface to a principal components analysis designed to be used within ...
Encapsulate the geometry of a wire.
geo::Geometry const * fGeometry
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
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
Utility object to perform functions of association.
bool fEnableMonitoring
Data members to follow.
Encapsulate the construction of a single detector plane.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
This provides an art tool interface definition for 3D Cluster algorithms.
std::vector< float > fTimeVector
std::pair< reco::ProjectedPoint, reco::ProjectedPoint > MinMaxPoints
Add ability to build the convex hull (these needs to be split out! )
float getHitChiSquare() const
Definition: Cluster3D.h:163
std::unordered_map< const reco::ClusterHit3D *, reco::EdgeList > Hit3DToEdgeMap
Definition: Cluster3D.h:339
void configure(fhicl::ParameterSet const &pset) override
float getTimeToExecute() const
Definition: kdTree.h:104
reco::ProjectedPointList & getConvexHullPointList()
Definition: Cluster3D.h:371
void buildConvexHull(reco::ClusterParameters &, reco::HitPairListPtr &, int level=0) const
float getTimeToExecute() const override
If monitoring, recover the time to execute a particular function.
reco::ProjectedPointList & getProjectedPointList()
Definition: Cluster3D.h:370
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
void ModifyClusters(reco::ClusterParametersList &) const override
Scan an input collection of clusters and modify those according to the specific implementing algorith...
reco::EdgeList & getConvexHullEdgeList()
Definition: Cluster3D.h:373
art framework interface to geometry description
std::list< ClusterParameters > ClusterParametersList
Definition: Cluster3D.h:393
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:243
attached to a cluster
Definition: Cluster3D.h:106
void AStar(const reco::ClusterHit3D *, const reco::ClusterHit3D *, reco::ClusterParameters &) const
Algorithm to find shortest path between two 3D hits.
void setStatusBit(unsigned bits) const
Definition: Cluster3D.h:177