MSTPathFinder_tool.cc
Go to the documentation of this file.
1 /**
2  * @file MSTPathFinder_tool.cc
3  *
4  * @brief art Tool for comparing clusters and merging those that are consistent
5  *
6  */
7 
8 // Framework Includes
11 #include "art_root_io/TFileService.h"
12 #include "cetlib/cpu_timer.h"
13 #include "cetlib/search_path.h"
14 
17 
18 // LArSoft includes
27 
28 // Eigen
29 #include <Eigen/Dense>
30 
31 // Root histograms
32 #include "TH1F.h"
33 #include "TH2F.h"
34 #include "TProfile.h"
35 
36 // std includes
37 #include <functional>
38 #include <iostream>
39 #include <memory>
40 #include <string>
41 #include <numeric> // std::accumulate
42 
43 //------------------------------------------------------------------------------------------------------------------------------------------
44 // implementation follows
45 
46 namespace lar_cluster3d {
47 
48  class MSTPathFinder : virtual public IClusterModAlg {
49  public:
50  /**
51  * @brief Constructor
52  *
53  * @param pset
54  */
55  explicit MSTPathFinder(const fhicl::ParameterSet&);
56 
57  /**
58  * @brief Destructor
59  */
61 
62  void configure(fhicl::ParameterSet const& pset) override;
63 
64  /**
65  * @brief Interface for initializing histograms if they are desired
66  * Note that the idea is to put hisgtograms in a subfolder
67  *
68  * @param TFileDirectory - the folder to store the hists in
69  */
70  void initializeHistograms(art::TFileDirectory&) override;
71 
72  /**
73  * @brief Scan an input collection of clusters and modify those according
74  * to the specific implementing algorithm
75  *
76  * @param clusterParametersList A list of cluster objects (parameters from associated hits)
77  */
78  void ModifyClusters(reco::ClusterParametersList&) const override;
79 
80  /**
81  * @brief If monitoring, recover the time to execute a particular function
82  */
83  float
84  getTimeToExecute() const override
85  {
86  return std::accumulate(fTimeVector.begin(), fTimeVector.end(), 0.);
87  }
88 
89  private:
90  /**
91  * @brief enumerate the possible values for time checking if monitoring timing
92  */
93  enum TimeValues {
96  RUNDBSCAN = 2,
100  };
101 
102  /**
103  * @brief Driver for Prim's algorithm
104  */
108 
109  /**
110  * @brief Prune the obvious ambiguous hits
111  */
113 
114  /**
115  * @brief Algorithm to find the best path through the given cluster
116  */
118 
119  /**
120  * @brief a depth first search to find longest branches
121  */
123  const reco::Hit3DToEdgeMap&,
124  float&) const;
125 
126  /**
127  * @brief Alternative version of FindBestPathInCluster utilizing an A* algorithm
128  */
130 
131  /**
132  * @brief Algorithm to find shortest path between two 3D hits
133  */
134  void AStar(const reco::ClusterHit3D*,
135  const reco::ClusterHit3D*,
136  float alpha,
138  reco::ClusterParameters&) const;
139 
140  using BestNodeTuple = std::tuple<const reco::ClusterHit3D*, float, float>;
141  using BestNodeMap = std::unordered_map<const reco::ClusterHit3D*, BestNodeTuple>;
142 
144  BestNodeMap&,
146  reco::EdgeList&) const;
147 
148  float DistanceBetweenNodes(const reco::ClusterHit3D*, const reco::ClusterHit3D*) const;
149 
150  /**
151  * @brief Find the lowest cost path between two nodes using MST edges
152  */
153  void LeastCostPath(const reco::EdgeTuple&,
154  const reco::ClusterHit3D*,
156  float&) const;
157 
158  /**
159  * @brief Add ability to build the convex hull (these needs to be split out! )
160  */
161  using MinMaxPoints = std::pair<reco::ProjectedPoint, reco::ProjectedPoint>;
162  using MinMaxPointPair = std::pair<MinMaxPoints, MinMaxPoints>;
163 
165 
167  const reco::ClusterHit3D*,
168  const reco::ClusterHit3D*) const;
169 
170  /**
171  * @brief Data members to follow
172  */
173  bool fEnableMonitoring; ///<
174  size_t fMinTinyClusterSize; ///< Minimum size for a "tiny" cluster
175  float fConvexHullKinkAngle; ///< Angle to declare a kink in convex hull calc
176  float fConvexHullMinSep; ///< Min hit separation to conisder in convex hull
177 
178  mutable std::vector<float> fTimeVector; ///<
179 
180  geo::Geometry const* fGeometry; //< pointer to the Geometry service
181 
182  PrincipalComponentsAlg fPCAAlg; // For running Principal Components Analysis
183  kdTree fkdTree; // For the kdTree
184 
185  std::unique_ptr<lar_cluster3d::IClusterParametersBuilder>
186  fClusterBuilder; ///< Common cluster builder tool
187  };
188 
190  : fPCAAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
191  , fkdTree(pset.get<fhicl::ParameterSet>("kdTree"))
192  {
193  this->configure(pset);
194  }
195 
196  //------------------------------------------------------------------------------------------------------------------------------------------
197 
199 
200  //------------------------------------------------------------------------------------------------------------------------------------------
201 
202  void
204  {
205  fEnableMonitoring = pset.get<bool>("EnableMonitoring", true);
206  fMinTinyClusterSize = pset.get<size_t>("MinTinyClusterSize", 40);
207  fConvexHullKinkAngle = pset.get<float>("ConvexHullKinkAgle", 0.95);
208  fConvexHullMinSep = pset.get<float>("ConvexHullMinSep", 0.65);
209 
211 
212  fGeometry = &*geometry;
213 
214  fTimeVector.resize(NUMTIMEVALUES, 0.);
215 
216  fClusterBuilder = art::make_tool<lar_cluster3d::IClusterParametersBuilder>(
217  pset.get<fhicl::ParameterSet>("ClusterParamsBuilder"));
218 
219  return;
220  }
221 
222  void
223  MSTPathFinder::initializeHistograms(art::TFileDirectory& histDir)
224  {
225  // It is assumed that the input TFileDirectory has been set up to group histograms into a common
226  // folder at the calling routine's level. Here we create one more level of indirection to keep
227  // histograms made by this tool separate.
228  // fFillHistograms = true;
229  //
230  // std::string dirName = "ConvexHullPath";
231  //
232  // art::TFileDirectory dir = histDir.mkdir(dirName.c_str());
233  //
234  // // Divide into two sets of hists... those for the overall cluster and
235  // // those for the subclusters
236  // fTopNum3DHits = dir.make<TH1F>("TopNum3DHits", "Number 3D Hits", 200, 0., 200.);
237 
238  return;
239  }
240 
241  void
243  {
244  /**
245  * @brief Top level interface tool for performing deghosting and primary path finding
246  * using a minimum spanning tree approach. This is a shell tool, it actually uses the
247  * Minimum Spanning Tree clusering tool...
248  */
249 
250  // Initial clustering is done, now trim the list and get output parameters
251  cet::cpu_timer theClockBuildClusters;
252 
253  // Start clocks if requested
254  if (fEnableMonitoring) theClockBuildClusters.start();
255 
256  // Ok, the idea here is to loop over the input clusters and the process one at a time and then use the MST algorithm
257  // to deghost and try to find the best path.
258  for (auto& clusterParams : clusterParametersList) {
259  // It turns out that computing the convex hull surrounding the points in the 2D projection onto the
260  // plane of largest spread in the PCA is a good way to break up the cluster... and we do it here since
261  // we (currently) want this to be part of the standard output
262  buildConvexHull(clusterParams, clusterParams.getHitPairListPtr());
263 
264  // Make sure our cluster has enough hits...
265  if (clusterParams.getHitPairListPtr().size() > fMinTinyClusterSize) {
266  // DBScan is driven of its "epsilon neighborhood". Computing adjacency within DBScan can be time
267  // consuming so the idea is the prebuild the adjaceny map and then run DBScan.
268  // The following call does this work
269  kdTree::KdTreeNodeList kdTreeNodeContainer;
270  kdTree::KdTreeNode topNode =
271  fkdTree.BuildKdTree(clusterParams.getHitPairListPtr(), kdTreeNodeContainer);
272 
274 
275  // We are making subclusters
276  reco::ClusterParametersList& daughterParametersList = clusterParams.daughterList();
277 
278  // Run DBScan to get candidate clusters
279  RunPrimsAlgorithm(clusterParams.getHitPairListPtr(), topNode, daughterParametersList);
280 
281  // Initial clustering is done, now trim the list and get output parameters
282  cet::cpu_timer theClockBuildClusters;
283 
284  // Start clocks if requested
285  if (fEnableMonitoring) theClockBuildClusters.start();
286 
287  fClusterBuilder->BuildClusterInfo(daughterParametersList);
288 
289  if (fEnableMonitoring) {
290  theClockBuildClusters.stop();
291 
292  fTimeVector[BUILDCLUSTERINFO] = theClockBuildClusters.accumulated_real_time();
293  }
294 
295  // Test run the path finding algorithm
296  for (auto& daughterParams : daughterParametersList)
297  FindBestPathInCluster(daughterParams, topNode);
298  }
299  }
300 
301  if (fEnableMonitoring) {
302  theClockBuildClusters.stop();
303 
304  fTimeVector[BUILDCLUSTERINFO] = theClockBuildClusters.accumulated_real_time();
305  }
306 
307  mf::LogDebug("MSTPathFinder") << ">>>>> Cluster Path finding done" << std::endl;
308 
309  return;
310  }
311 
312  //------------------------------------------------------------------------------------------------------------------------------------------
313  void
315  kdTree::KdTreeNode& topNode,
316  reco::ClusterParametersList& clusterParametersList) const
317  {
318  // If no hits then no work
319  if (hitPairList.empty()) return;
320 
321  // Now proceed with building the clusters
322  cet::cpu_timer theClockDBScan;
323 
324  // Start clocks if requested
325  if (fEnableMonitoring) theClockDBScan.start();
326 
327  // Initialization
328  size_t clusterIdx(0);
329 
330  // This will contain our list of edges
331  reco::EdgeList curEdgeList;
332 
333  // Get the first point
334  reco::HitPairListPtr::const_iterator freeHitItr = hitPairList.begin();
335  const reco::ClusterHit3D* lastAddedHit = *freeHitItr++;
336 
338 
339  // Make a cluster...
340  clusterParametersList.push_back(reco::ClusterParameters());
341 
342  // Get an iterator to the first cluster
343  reco::ClusterParameters* curCluster = &clusterParametersList.back();
344 
345  // We use pointers here because the objects they point to will change in the loop below
346  reco::Hit3DToEdgeMap* curEdgeMap = &curCluster->getHit3DToEdgeMap();
347  reco::HitPairListPtr* curClusterHitList = &curCluster->getHitPairListPtr();
348 
349  // Loop until all hits have been associated to a cluster
350  while (1) {
351  // and the 3D hit status bits
353 
354  // Purge the current list to get rid of edges which point to hits already in the cluster
355  for (reco::EdgeList::iterator curEdgeItr = curEdgeList.begin();
356  curEdgeItr != curEdgeList.end();) {
357  if (std::get<1>(*curEdgeItr)->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED)
358  curEdgeItr = curEdgeList.erase(curEdgeItr);
359  else
360  curEdgeItr++;
361  }
362 
363  // Add the lastUsedHit to the current cluster
364  curClusterHitList->push_back(lastAddedHit);
365 
366  // Set up to find the list of nearest neighbors to the last used hit...
367  kdTree::CandPairList CandPairList;
368  float bestDistance(1.5); //std::numeric_limits<float>::max());
369 
370  // And find them... result will be an unordered list of neigbors
371  fkdTree.FindNearestNeighbors(lastAddedHit, topNode, CandPairList, bestDistance);
372 
373  // Copy edges to the current list (but only for hits not already in a cluster)
374  // for(auto& pair : CandPairList)
375  // if (!(pair.second->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED)) curEdgeList.push_back(reco::EdgeTuple(lastAddedHit,pair.second,pair.first));
376  for (auto& pair : CandPairList) {
377  if (!(pair.second->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED)) {
378  double edgeWeight =
379  pair.first * lastAddedHit->getHitChiSquare() * pair.second->getHitChiSquare();
380 
381  curEdgeList.push_back(reco::EdgeTuple(lastAddedHit, pair.second, edgeWeight));
382  }
383  }
384 
385  // If the edge list is empty then we have a complete cluster
386  if (curEdgeList.empty()) {
387  std::cout << "-----------------------------------------------------------------------------"
388  "------------"
389  << std::endl;
390  std::cout << "**> Cluster idx: " << clusterIdx++ << " has " << curClusterHitList->size()
391  << " hits" << std::endl;
392 
393  // Look for the next "free" hit
394  freeHitItr = std::find_if(freeHitItr, hitPairList.end(), [](const auto& hit) {
395  return !(hit->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED);
396  });
397 
398  // If at end of input list we are done with all hits
399  if (freeHitItr == hitPairList.end()) break;
400 
401  std::cout << "##################################################################>"
402  "Processing another cluster"
403  << std::endl;
404 
405  // Otherwise, get a new cluster and set up
406  clusterParametersList.push_back(reco::ClusterParameters());
407 
408  curCluster = &clusterParametersList.back();
409 
410  curEdgeMap = &curCluster->getHit3DToEdgeMap();
411  curClusterHitList = &curCluster->getHitPairListPtr();
412  lastAddedHit = *freeHitItr++;
413  }
414  // Otherwise we are still processing the current cluster
415  else {
416  // Sort the list of edges by distance
417  curEdgeList.sort([](const auto& left, const auto& right) {
418  return std::get<2>(left) < std::get<2>(right);
419  });
420 
421  // Populate the map with the edges...
422  reco::EdgeTuple& curEdge = curEdgeList.front();
423 
424  (*curEdgeMap)[std::get<0>(curEdge)].push_back(curEdge);
425  (*curEdgeMap)[std::get<1>(curEdge)].push_back(
426  reco::EdgeTuple(std::get<1>(curEdge), std::get<0>(curEdge), std::get<2>(curEdge)));
427 
428  // Update the last hit to be added to the collection
429  lastAddedHit = std::get<1>(curEdge);
430  }
431  }
432 
433  if (fEnableMonitoring) {
434  theClockDBScan.stop();
435 
436  fTimeVector[RUNDBSCAN] = theClockDBScan.accumulated_real_time();
437  }
438 
439  return;
440  }
441 
442  void
444  {
445  reco::HitPairListPtr longestCluster;
446  float bestQuality(0.);
447  float aveNumEdges(0.);
448  size_t maxNumEdges(0);
449  size_t nIsolatedHits(0);
450 
451  // Now proceed with building the clusters
452  cet::cpu_timer theClockPathFinding;
453 
454  // Start clocks if requested
455  if (fEnableMonitoring) theClockPathFinding.start();
456 
457  reco::HitPairListPtr& hitPairList = curCluster.getHitPairListPtr();
458  reco::Hit3DToEdgeMap& curEdgeMap = curCluster.getHit3DToEdgeMap();
459  reco::EdgeList& bestEdgeList = curCluster.getBestEdgeList();
460 
461  // Do some spelunking...
462  for (const auto& hit : hitPairList) {
463  if (!curEdgeMap[hit].empty() && curEdgeMap[hit].size() == 1) {
464  float quality(0.);
465 
466  reco::HitPairListPtr tempList =
467  DepthFirstSearch(curEdgeMap[hit].front(), curEdgeMap, quality);
468 
469  tempList.push_front(std::get<0>(curEdgeMap[hit].front()));
470 
471  if (quality > bestQuality) {
472  longestCluster = tempList;
473  bestQuality = quality;
474  }
475 
476  nIsolatedHits++;
477  }
478 
479  aveNumEdges += float(curEdgeMap[hit].size());
480  maxNumEdges = std::max(maxNumEdges, curEdgeMap[hit].size());
481  }
482 
483  aveNumEdges /= float(hitPairList.size());
484  std::cout << "----> # isolated hits: " << nIsolatedHits
485  << ", longest branch: " << longestCluster.size()
486  << ", cluster size: " << hitPairList.size() << ", ave # edges: " << aveNumEdges
487  << ", max: " << maxNumEdges << std::endl;
488 
489  if (!longestCluster.empty()) {
490  hitPairList = longestCluster;
491  for (const auto& hit : hitPairList) {
492  for (const auto& edge : curEdgeMap[hit])
493  bestEdgeList.emplace_back(edge);
494  }
495 
496  std::cout << " ====> new cluster size: " << hitPairList.size() << std::endl;
497  }
498 
499  if (fEnableMonitoring) {
500  theClockPathFinding.stop();
501 
502  fTimeVector[PATHFINDING] += theClockPathFinding.accumulated_real_time();
503  }
504 
505  return;
506  }
507 
508  void
510  kdTree::KdTreeNode& topNode) const
511  {
512  // Set up for timing the function
513  cet::cpu_timer theClockPathFinding;
514 
515  // Start clocks if requested
516  if (fEnableMonitoring) theClockPathFinding.start();
517 
518  // Trial A* here
519  if (clusterParams.getHitPairListPtr().size() > 2) {
520  // Do a quick PCA to determine our parameter "alpha"
521  fPCAAlg.PCAAnalysis_3D(clusterParams.getHitPairListPtr(), clusterParams.getFullPCA());
522 
523  // Recover the new fullPCA
524  reco::PrincipalComponents& pca = clusterParams.getFullPCA();
525 
526  // The chances of a failure are remote, still we should check
527  if (pca.getSvdOK()) {
528  // Get references to what we need....
529  reco::HitPairListPtr& curCluster = clusterParams.getHitPairListPtr();
530  reco::Hit3DToEdgeMap& curEdgeMap = clusterParams.getHit3DToEdgeMap();
531  float pcaLen = 3.0 * sqrt(pca.getEigenValues()[2]);
532  float pcaWidth = 3.0 * sqrt(pca.getEigenValues()[1]);
533  float pcaHeight = 3.0 * sqrt(pca.getEigenValues()[0]);
534  const Eigen::Vector3f& pcaCenter = pca.getAvePosition();
535  float alpha = std::min(float(1.), std::max(float(0.001), pcaWidth / pcaLen));
536 
537  // Create a temporary container for the isolated points
538  reco::ProjectedPointList isolatedPointList;
539 
540  // Go through and find the isolated points, for those get the projection to the plane of maximum spread
541  for (const auto& hit3D : curCluster) {
542  // the definition of an isolated hit is that it only has one associated edge
543  if (!curEdgeMap[hit3D].empty() && curEdgeMap[hit3D].size() == 1) {
544  Eigen::Vector3f pcaToHitVec(hit3D->getPosition()[0] - pcaCenter(0),
545  hit3D->getPosition()[1] - pcaCenter(1),
546  hit3D->getPosition()[2] - pcaCenter(2));
547  Eigen::Vector3f pcaToHit = pca.getEigenVectors() * pcaToHitVec;
548 
549  // This sets x,y where x is the longer spread, y the shorter
550  isolatedPointList.emplace_back(pcaToHit(2), pcaToHit(1), hit3D);
551  }
552  }
553 
554  std::cout << "************* Finding best path with A* in cluster *****************"
555  << std::endl;
556  std::cout << "**> There are " << curCluster.size() << " hits, " << isolatedPointList.size()
557  << " isolated hits, the alpha parameter is " << alpha << std::endl;
558  std::cout << "**> PCA len: " << pcaLen << ", wid: " << pcaWidth << ", height: " << pcaHeight
559  << ", ratio: " << pcaHeight / pcaWidth << std::endl;
560 
561  // If no isolated points then nothing to do...
562  if (isolatedPointList.size() > 1) {
563  // Sort the point vec by increasing x, if same then by increasing y.
564  isolatedPointList.sort([](const auto& left, const auto& right) {
565  return (std::abs(std::get<0>(left) - std::get<0>(right)) >
566  std::numeric_limits<float>::epsilon()) ?
567  std::get<0>(left) < std::get<0>(right) :
568  std::get<1>(left) < std::get<1>(right);
569  });
570 
571  // Ok, get the two most distance points...
572  const reco::ClusterHit3D* startHit = std::get<2>(isolatedPointList.front());
573  const reco::ClusterHit3D* stopHit = std::get<2>(isolatedPointList.back());
574 
575  std::cout << "**> Sorted " << isolatedPointList.size()
576  << " hits, longest distance: " << DistanceBetweenNodes(startHit, stopHit)
577  << std::endl;
578 
579  // Call the AStar function to try to find the best path...
580  // AStar(startHit,stopHit,alpha,topNode,clusterParams);
581 
582  float cost(std::numeric_limits<float>::max());
583 
584  LeastCostPath(curEdgeMap[startHit].front(), stopHit, clusterParams, cost);
585 
586  clusterParams.getBestHitPairListPtr().push_front(startHit);
587 
588  std::cout << "**> Best path has " << clusterParams.getBestHitPairListPtr().size()
589  << " hits, " << clusterParams.getBestEdgeList().size() << " edges" << std::endl;
590  }
591 
592  // Recalculate the PCA based on the hits comprisig the path
593  fPCAAlg.PCAAnalysis_3D(clusterParams.getBestHitPairListPtr(), pca);
594 
595  // And now compute the convex hull
596  buildConvexHull(clusterParams, clusterParams.getBestHitPairListPtr());
597  }
598  else {
599  std::cout << "++++++>>> PCA failure! # hits: " << clusterParams.getHitPairListPtr().size()
600  << std::endl;
601  }
602  }
603 
604  if (fEnableMonitoring) {
605  theClockPathFinding.stop();
606 
607  fTimeVector[PATHFINDING] += theClockPathFinding.accumulated_real_time();
608  }
609 
610  return;
611  }
612 
613  void
615  const reco::ClusterHit3D* goalNode,
616  float alpha,
617  kdTree::KdTreeNode& topNode,
618  reco::ClusterParameters& clusterParams) const
619  {
620  // Recover the list of hits and edges
621  reco::HitPairListPtr& pathNodeList = clusterParams.getBestHitPairListPtr();
622  reco::EdgeList& bestEdgeList = clusterParams.getBestEdgeList();
623  reco::Hit3DToEdgeMap& curEdgeMap = clusterParams.getHit3DToEdgeMap();
624 
625  // Find the shortest path from start to goal using an A* algorithm
626  // Keep track of the nodes which have already been evaluated
627  reco::HitPairListPtr closedList;
628 
629  // Keep track of the nodes that have been "discovered" but yet to be evaluated
630  reco::HitPairListPtr openList = {startNode};
631 
632  // Create a map which, for each node, will tell us the node it can be most efficiencly reached from.
633  BestNodeMap bestNodeMap;
634 
635  bestNodeMap[startNode] =
636  BestNodeTuple(startNode, 0., DistanceBetweenNodes(startNode, goalNode));
637 
638  alpha = 1.; //std::max(0.5,alpha);
639 
640  while (!openList.empty()) {
641  // The list is not empty so by def we will return something
642  reco::HitPairListPtr::iterator currentNodeItr = openList.begin();
643 
644  // If the list contains more than one element then we need to find the one with the smallest total estimated cost to the end
645  if (openList.size() > 1)
646  currentNodeItr = std::min_element(
647  openList.begin(), openList.end(), [bestNodeMap](const auto& next, const auto& best) {
648  return std::get<2>(bestNodeMap.at(next)) < std::get<2>(bestNodeMap.at(best));
649  });
650 
651  // Easier to deal directly with the pointer to the node
652  const reco::ClusterHit3D* currentNode = *currentNodeItr;
653 
654  // Check to see if we have reached the goal and need to evaluate the path
655  if (currentNode == goalNode) {
656  // The path reconstruction will
657  ReconstructBestPath(goalNode, bestNodeMap, pathNodeList, bestEdgeList);
658 
659  break;
660  }
661 
662  // Otherwise need to keep evaluating
663  else {
664  openList.erase(currentNodeItr);
666 
667  // Get tuple values for the current node
668  const BestNodeTuple& currentNodeTuple = bestNodeMap.at(currentNode);
669  float currentNodeScore = std::get<1>(currentNodeTuple);
670 
671  // Recover the edges associated to the current point
672  const reco::EdgeList& curEdgeList = curEdgeMap[currentNode];
673 
674  for (const auto& curEdge : curEdgeList) {
675  const reco::ClusterHit3D* candHit3D = std::get<1>(curEdge);
676 
677  if (candHit3D->getStatusBits() & reco::ClusterHit3D::PATHCHECKED) continue;
678 
679  float tentative_gScore = currentNodeScore + std::get<2>(curEdge);
680 
681  // Have we seen the candidate node before?
682  BestNodeMap::iterator candNodeItr = bestNodeMap.find(candHit3D);
683 
684  if (candNodeItr == bestNodeMap.end()) { openList.push_back(candHit3D); }
685  else if (tentative_gScore > std::get<1>(candNodeItr->second))
686  continue;
687 
688  // Make a guess at score to get to target...
689  float guessToTarget = DistanceBetweenNodes(candHit3D, goalNode) / 0.3;
690 
691  bestNodeMap[candHit3D] =
692  BestNodeTuple(currentNode, tentative_gScore, tentative_gScore + guessToTarget);
693  }
694  }
695  }
696 
697  return;
698  }
699 
700  void
702  BestNodeMap& bestNodeMap,
703  reco::HitPairListPtr& pathNodeList,
704  reco::EdgeList& bestEdgeList) const
705  {
706  while (std::get<0>(bestNodeMap.at(goalNode)) != goalNode) {
707  const reco::ClusterHit3D* nextNode = std::get<0>(bestNodeMap[goalNode]);
708  reco::EdgeTuple bestEdge =
709  reco::EdgeTuple(goalNode, nextNode, DistanceBetweenNodes(goalNode, nextNode));
710 
711  pathNodeList.push_front(goalNode);
712  bestEdgeList.push_front(bestEdge);
713 
714  goalNode = nextNode;
715  }
716 
717  pathNodeList.push_front(goalNode);
718 
719  return;
720  }
721 
722  void
724  const reco::ClusterHit3D* goalNode,
725  reco::ClusterParameters& clusterParams,
726  float& showMeTheMoney) const
727  {
728  // Recover the mapping between hits and edges
729  reco::Hit3DToEdgeMap& curEdgeMap = clusterParams.getHit3DToEdgeMap();
730 
731  reco::Hit3DToEdgeMap::const_iterator edgeListItr = curEdgeMap.find(std::get<1>(curEdge));
732 
733  showMeTheMoney = std::numeric_limits<float>::max();
734 
735  if (edgeListItr != curEdgeMap.end() && !edgeListItr->second.empty()) {
736  reco::HitPairListPtr& bestNodeList = clusterParams.getBestHitPairListPtr();
737  reco::EdgeList& bestEdgeList = clusterParams.getBestEdgeList();
738 
739  for (const auto& edge : edgeListItr->second) {
740  // skip the self reference
741  if (std::get<1>(edge) == std::get<0>(curEdge)) continue;
742 
743  // Have we found the droid we are looking for?
744  if (std::get<1>(edge) == goalNode) {
745  bestNodeList.push_back(goalNode);
746  bestEdgeList.push_back(edge);
747  showMeTheMoney = std::get<2>(edge);
748  break;
749  }
750 
751  // Keep searching, it is out there somewhere...
752  float currentCost(0.);
753 
754  LeastCostPath(edge, goalNode, clusterParams, currentCost);
755 
756  if (currentCost < std::numeric_limits<float>::max()) {
757  showMeTheMoney = std::get<2>(edge) + currentCost;
758  break;
759  }
760  }
761  }
762 
763  if (showMeTheMoney < std::numeric_limits<float>::max()) {
764  clusterParams.getBestHitPairListPtr().push_front(std::get<1>(curEdge));
765  clusterParams.getBestEdgeList().push_front(curEdge);
766  }
767 
768  return;
769  }
770 
771  float
773  const reco::ClusterHit3D* node2) const
774  {
775  const Eigen::Vector3f& node1Pos = node1->getPosition();
776  const Eigen::Vector3f& node2Pos = node2->getPosition();
777  float deltaNode[] = {
778  node1Pos[0] - node2Pos[0], node1Pos[1] - node2Pos[1], node1Pos[2] - node2Pos[2]};
779 
780  // Standard euclidean distance
781  return std::sqrt(deltaNode[0] * deltaNode[0] + deltaNode[1] * deltaNode[1] +
782  deltaNode[2] * deltaNode[2]);
783 
784  // Manhatten distance
785  //return std::fabs(deltaNode[0]) + std::fabs(deltaNode[1]) + std::fabs(deltaNode[2]);
786  /*
787  // Chebyshev distance
788  // We look for maximum distance by wires...
789 
790  // Now go through the hits and compare view by view to look for delta wire and tigher constraint on delta t
791  int wireDeltas[] = {0,0,0};
792 
793  for(size_t idx = 0; idx < 3; idx++)
794  wireDeltas[idx] = std::abs(int(node1->getWireIDs()[idx].Wire) - int(node2->getWireIDs()[idx].Wire));
795 
796  // put wire deltas in order...
797  std::sort(wireDeltas, wireDeltas + 3);
798 
799  return std::sqrt(deltaNode[0]*deltaNode[0] + 0.09 * float(wireDeltas[2]*wireDeltas[2]));
800  */
801  }
802 
805  const reco::Hit3DToEdgeMap& hitToEdgeMap,
806  float& bestTreeQuality) const
807  {
808  reco::HitPairListPtr hitPairListPtr;
809  float bestQuality(0.);
810  float curEdgeWeight = std::max(0.3, std::get<2>(curEdge));
811  float curEdgeProj(1. / curEdgeWeight);
812 
813  reco::Hit3DToEdgeMap::const_iterator edgeListItr = hitToEdgeMap.find(std::get<1>(curEdge));
814 
815  if (edgeListItr != hitToEdgeMap.end()) {
816  // The input edge weight has quality factors applied, recalculate just the position difference
817  const Eigen::Vector3f& firstHitPos = std::get<0>(curEdge)->getPosition();
818  const Eigen::Vector3f& secondHitPos = std::get<1>(curEdge)->getPosition();
819  float curEdgeVec[] = {secondHitPos[0] - firstHitPos[0],
820  secondHitPos[1] - firstHitPos[1],
821  secondHitPos[2] - firstHitPos[2]};
822  float curEdgeMag = std::sqrt(curEdgeVec[0] * curEdgeVec[0] + curEdgeVec[1] * curEdgeVec[1] +
823  curEdgeVec[2] * curEdgeVec[2]);
824 
825  curEdgeMag = std::max(float(0.1), curEdgeMag);
826 
827  for (const auto& edge : edgeListItr->second) {
828  // skip the self reference
829  if (std::get<1>(edge) == std::get<0>(curEdge)) continue;
830 
831  float quality(0.);
832 
833  reco::HitPairListPtr tempList = DepthFirstSearch(edge, hitToEdgeMap, quality);
834 
835  if (quality > bestQuality) {
836  hitPairListPtr = tempList;
837  bestQuality = quality;
838  curEdgeProj = 1. / curEdgeMag;
839  }
840  }
841  }
842 
843  hitPairListPtr.push_front(std::get<1>(curEdge));
844 
845  bestTreeQuality += bestQuality + curEdgeProj;
846 
847  return hitPairListPtr;
848  }
849 
850  void
852  reco::Hit2DToClusterMap& hit2DToClusterMap) const
853  {
854 
855  // Recover the HitPairListPtr from the input clusterParams (which will be the
856  // only thing that has been provided)
857  reco::HitPairListPtr& hitPairVector = clusterParams.getHitPairListPtr();
858 
859  size_t nStartedWith(hitPairVector.size());
860  size_t nRejectedHits(0);
861 
862  reco::HitPairListPtr goodHits;
863 
864  // Loop through the hits and try to week out the clearly ambiguous ones
865  for (const auto& hit3D : hitPairVector) {
866  // Loop to try to remove ambiguous hits
867  size_t n2DHitsIn3DHit(0);
868  size_t nThisClusterOnly(0);
869  size_t nOtherCluster(0);
870 
871  // reco::ClusterParameters* otherCluster;
872  const std::set<const reco::ClusterHit3D*>* otherClusterHits = 0;
873 
874  for (const auto& hit2D : hit3D->getHits()) {
875  if (!hit2D) continue;
876 
877  n2DHitsIn3DHit++;
878 
879  if (hit2DToClusterMap[hit2D].size() < 2)
880  nThisClusterOnly = hit2DToClusterMap[hit2D][&clusterParams].size();
881  else {
882  for (const auto& clusterHitMap : hit2DToClusterMap[hit2D]) {
883  if (clusterHitMap.first == &clusterParams) continue;
884 
885  if (clusterHitMap.second.size() > nOtherCluster) {
886  nOtherCluster = clusterHitMap.second.size();
887  // otherCluster = clusterHitMap.first;
888  otherClusterHits = &clusterHitMap.second;
889  }
890  }
891  }
892  }
893 
894  if (n2DHitsIn3DHit < 3 && nThisClusterOnly > 1 && nOtherCluster > 0) {
895  bool skip3DHit(false);
896 
897  for (const auto& otherHit3D : *otherClusterHits) {
898  size_t nOther2DHits(0);
899 
900  for (const auto& otherHit2D : otherHit3D->getHits()) {
901  if (!otherHit2D) continue;
902 
903  nOther2DHits++;
904  }
905 
906  if (nOther2DHits > 2) {
907  skip3DHit = true;
908  nRejectedHits++;
909  break;
910  }
911  }
912 
913  if (skip3DHit) continue;
914  }
915 
916  goodHits.emplace_back(hit3D);
917  }
918 
919  std::cout << "###>> Input " << nStartedWith << " hits, rejected: " << nRejectedHits
920  << std::endl;
921 
922  hitPairVector.resize(goodHits.size());
923  std::copy(goodHits.begin(), goodHits.end(), hitPairVector.begin());
924 
925  return;
926  }
927 
928  void
930  reco::HitPairListPtr& hitPairListPtr,
931  int level) const
932  {
933  // set an indention string
934  std::string minuses(level / 2, '-');
935  std::string indent(level / 2, ' ');
936 
937  indent += minuses;
938 
939  // The plan is to build the enclosing 2D polygon around the points in the PCA plane of most spread for this cluster
940  // To do so we need to start by building a list of 2D projections onto the plane of most spread...
941  reco::PrincipalComponents& pca = clusterParameters.getFullPCA();
942 
943  // Recover the parameters from the Principal Components Analysis that we need to project and accumulate
944  const Eigen::Vector3f& pcaCenter = pca.getAvePosition();
945  reco::ConvexHull& convexHull = clusterParameters.getConvexHull();
946  reco::ProjectedPointList& pointList = convexHull.getProjectedPointList();
947 
948  // Loop through hits and do projection to plane
949  for (const auto& hit3D : hitPairListPtr) {
950  Eigen::Vector3f pcaToHitVec(hit3D->getPosition()[0] - pcaCenter(0),
951  hit3D->getPosition()[1] - pcaCenter(1),
952  hit3D->getPosition()[2] - pcaCenter(2));
953  Eigen::Vector3f pcaToHit = pca.getEigenVectors() * pcaToHitVec;
954 
955  pointList.emplace_back(pcaToHit(1), pcaToHit(2), hit3D);
956  }
957 
958  // Sort the point vec by increasing x, then increase y
959  pointList.sort([](const auto& left, const auto& right) {
960  return (std::abs(std::get<0>(left) - std::get<0>(right)) >
961  std::numeric_limits<float>::epsilon()) ?
962  std::get<0>(left) < std::get<0>(right) :
963  std::get<1>(left) < std::get<1>(right);
964  });
965 
966  // containers for finding the "best" hull...
967  std::vector<ConvexHull> convexHullVec;
968  std::vector<reco::ProjectedPointList> rejectedListVec;
969  bool increaseDepth(pointList.size() > 3);
970  float lastArea(std::numeric_limits<float>::max());
971 
972  while (increaseDepth) {
973  // Get another convexHull container
974  convexHullVec.push_back(ConvexHull(pointList, fConvexHullKinkAngle, fConvexHullMinSep));
975  rejectedListVec.push_back(reco::ProjectedPointList());
976 
977  const ConvexHull& convexHull = convexHullVec.back();
978  reco::ProjectedPointList& rejectedList = rejectedListVec.back();
979  const reco::ProjectedPointList& convexHullPoints = convexHull.getConvexHull();
980 
981  increaseDepth = false;
982 
983  if (convexHull.getConvexHullArea() > 0.) {
984  if (convexHullVec.size() < 2 || convexHull.getConvexHullArea() < 0.8 * lastArea) {
985  for (auto& point : convexHullPoints) {
986  pointList.remove(point);
987  rejectedList.emplace_back(point);
988  }
989  lastArea = convexHull.getConvexHullArea();
990  // increaseDepth = true;
991  }
992  }
993  }
994 
995  // do we have a valid convexHull?
996  while (!convexHullVec.empty() && convexHullVec.back().getConvexHullArea() < 0.5) {
997  convexHullVec.pop_back();
998  rejectedListVec.pop_back();
999  }
1000 
1001  // If we found the convex hull then build edges around the region
1002  if (!convexHullVec.empty()) {
1003  size_t nRejectedTotal(0);
1004  reco::HitPairListPtr locHitPairListPtr = hitPairListPtr;
1005 
1006  for (const auto& rejectedList : rejectedListVec) {
1007  nRejectedTotal += rejectedList.size();
1008 
1009  for (const auto& rejectedPoint : rejectedList) {
1010  if (convexHullVec.back().findNearestDistance(rejectedPoint) > 0.5)
1011  locHitPairListPtr.remove(std::get<2>(rejectedPoint));
1012  }
1013  }
1014 
1015  // Now add "edges" to the cluster to describe the convex hull (for the display)
1016  reco::ProjectedPointList& convexHullPointList = convexHull.getConvexHullPointList();
1017  reco::Hit3DToEdgeMap& edgeMap = convexHull.getConvexHullEdgeMap();
1018  reco::EdgeList& edgeList = convexHull.getConvexHullEdgeList();
1019 
1020  reco::ProjectedPoint lastPoint = convexHullVec.back().getConvexHull().front();
1021 
1022  for (auto& curPoint : convexHullVec.back().getConvexHull()) {
1023  if (curPoint == lastPoint) continue;
1024 
1025  const reco::ClusterHit3D* lastPoint3D = std::get<2>(lastPoint);
1026  const reco::ClusterHit3D* curPoint3D = std::get<2>(curPoint);
1027 
1028  float distBetweenPoints = (curPoint3D->getPosition()[0] - lastPoint3D->getPosition()[0]) *
1029  (curPoint3D->getPosition()[0] - lastPoint3D->getPosition()[0]) +
1030  (curPoint3D->getPosition()[1] - lastPoint3D->getPosition()[1]) *
1031  (curPoint3D->getPosition()[1] - lastPoint3D->getPosition()[1]) +
1032  (curPoint3D->getPosition()[2] - lastPoint3D->getPosition()[2]) *
1033  (curPoint3D->getPosition()[2] - lastPoint3D->getPosition()[2]);
1034 
1035  distBetweenPoints = std::sqrt(distBetweenPoints);
1036 
1037  reco::EdgeTuple edge(lastPoint3D, curPoint3D, distBetweenPoints);
1038 
1039  convexHullPointList.push_back(curPoint);
1040  edgeMap[lastPoint3D].push_back(edge);
1041  edgeMap[curPoint3D].push_back(edge);
1042  edgeList.emplace_back(edge);
1043 
1044  lastPoint = curPoint;
1045  }
1046 
1047  // Store the "extreme" points
1048  const ConvexHull::PointList& extremePoints = convexHullVec.back().getExtremePoints();
1049  reco::ProjectedPointList& extremePointList = convexHull.getConvexHullExtremePoints();
1050 
1051  for (const auto& point : extremePoints)
1052  extremePointList.push_back(point);
1053 
1054  // Store the "kink" points
1055  const reco::ConvexHullKinkTupleList& kinkPoints = convexHullVec.back().getKinkPoints();
1056  reco::ConvexHullKinkTupleList& kinkPointList = convexHull.getConvexHullKinkPoints();
1057 
1058  for (const auto& kink : kinkPoints)
1059  kinkPointList.push_back(kink);
1060  }
1061 
1062  return;
1063  }
1064 
1065  float
1067  const reco::ClusterHit3D* first3D,
1068  const reco::ClusterHit3D* last3D) const
1069  {
1070  float largestDistance(0.);
1071 
1072  // Note that edges are vectors and that the convex hull edge list will be ordered
1073  // The idea is that the maximum distance from a given edge is to the edge just before the edge that "turns back" towards the current edge
1074  // meaning that the dot product of the two edges becomes negative.
1075  reco::EdgeList::const_iterator firstEdgeItr = convexHull.begin();
1076 
1077  while (firstEdgeItr != convexHull.end()) {
1078  reco::EdgeList::const_iterator nextEdgeItr = firstEdgeItr;
1079 
1080  // Eigen::Vector2f firstEdgeVec(std::get<3>(*firstEdgeItr),std::get<);
1081  // Eigen::Vector2f lastPrimaryVec(lastPCA.getEigenVectors()[0][0],lastPCA.getEigenVectors()[0][1],lastPCA.getEigenVectors()[0][2]);
1082  // float cosToLast = newPrimaryVec.dot(lastPrimaryVec);
1083 
1084  while (++nextEdgeItr != convexHull.end()) {}
1085  }
1086 
1087  return largestDistance;
1088  }
1089 
1091 } // namespace lar_cluster3d
reco::HitPairListPtr & getBestHitPairListPtr()
Definition: Cluster3D.h:480
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:244
#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:30
void FindBestPathInCluster(reco::ClusterParameters &) const
Algorithm to find the best path through the given cluster.
std::string string
Definition: nybbler.cc:12
const Eigen::Vector3f getPosition() const
Definition: Cluster3D.h:158
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:353
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:34
std::list< KdTreeNode > KdTreeNodeList
Definition: kdTree.h:67
Implements a kdTree for use in clustering.
size_t FindNearestNeighbors(const reco::ClusterHit3D *, const KdTreeNode &, CandPairList &, float &) const
Definition: kdTree.cxx:170
reco::EdgeList & getBestEdgeList()
Definition: Cluster3D.h:481
reco::Hit3DToEdgeMap & getHit3DToEdgeMap()
Definition: Cluster3D.h:479
intermediate_table::const_iterator const_iterator
reco::HitPairListPtr & getHitPairListPtr()
Definition: Cluster3D.h:476
float fConvexHullMinSep
Min hit separation to conisder in convex hull.
MSTPathFinder(const fhicl::ParameterSet &)
Constructor.
unsigned int getStatusBits() const
Definition: Cluster3D.h:157
art framework interface to geometry description
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:92
Implements a ConvexHull for use in clustering.
T abs(T value)
std::list< EdgeTuple > EdgeList
Definition: Cluster3D.h:345
define a kd tree node
Definition: kdTree.h:118
fInnerVessel push_back(Point(-578.400000, 0.000000, 0.000000))
const EigenValues & getEigenValues() const
Definition: Cluster3D.h:246
reco::PrincipalComponents & getFullPCA()
Definition: Cluster3D.h:477
reco::Hit3DToEdgeMap & getConvexHullEdgeMap()
Definition: Cluster3D.h:385
void LeastCostPath(const reco::EdgeTuple &, const reco::ClusterHit3D *, reco::ClusterParameters &, float &) const
Find the lowest cost path between two nodes using MST edges.
reco::ProjectedPointList & getConvexHullExtremePoints()
Definition: Cluster3D.h:387
reco::ConvexHullKinkTupleList & getConvexHullKinkPoints()
Definition: Cluster3D.h:388
Define a container for working with the convex hull.
Definition: Cluster3D.h:360
T get(std::string const &key) const
Definition: ParameterSet.h:271
std::unique_ptr< lar_cluster3d::IClusterParametersBuilder > fClusterBuilder
Common cluster builder tool.
const Eigen::Vector3f & getAvePosition() const
Definition: Cluster3D.h:248
double alpha
Definition: doAna.cpp:15
Path checking algorithm has seen this hit.
Definition: Cluster3D.h:111
void ReconstructBestPath(const reco::ClusterHit3D *, BestNodeMap &, reco::HitPairListPtr &, reco::EdgeList &) const
std::list< ConvexHullKinkTuple > ConvexHullKinkTupleList
Definition: Cluster3D.h:355
std::tuple< const reco::ClusterHit3D *, const reco::ClusterHit3D *, double > EdgeTuple
Definition: Cluster3D.h:344
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:335
float getConvexHullArea() const
recover the area of the convex hull
Definition: ConvexHull.h:78
static int max(int a, int b)
The geometry of one entire detector, as served by art.
Definition: Geometry.h:196
const PointList & getConvexHull() const
recover the list of convex hull vertices
Definition: ConvexHull.h:58
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:79
ConvexHull class definiton.
Definition: ConvexHull.h:27
Encapsulate the geometry of a wire.
geo::Geometry const * fGeometry
float findConvexHullEndPoints(const reco::EdgeList &, const reco::ClusterHit3D *, const reco::ClusterHit3D *) const
reco::ConvexHull & getConvexHull()
Definition: Cluster3D.h:482
Declaration of signal hit object.
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
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
void AStar(const reco::ClusterHit3D *, const reco::ClusterHit3D *, float alpha, kdTree::KdTreeNode &, reco::ClusterParameters &) const
Algorithm to find shortest path between two 3D hits.
T copy(T const &v)
std::pair< reco::ProjectedPoint, reco::ProjectedPoint > MinMaxPoints
Add ability to build the convex hull (these needs to be split out! )
std::unordered_map< const reco::ClusterHit2D *, ClusterToHitPairSetMap > Hit2DToClusterMap
Definition: Cluster3D.h:511
float getHitChiSquare() const
Definition: Cluster3D.h:166
std::unordered_map< const reco::ClusterHit3D *, reco::EdgeList > Hit3DToEdgeMap
Definition: Cluster3D.h:347
double accumulated_real_time() const
Definition: cpu_timer.cc:66
void configure(fhicl::ParameterSet const &pset) override
float getTimeToExecute() const
Definition: kdTree.h:95
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
reco::ProjectedPointList & getConvexHullPointList()
Definition: Cluster3D.h:384
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:383
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:112
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:386
void start()
Definition: cpu_timer.cc:83
std::tuple< float, float, const reco::ClusterHit3D * > ProjectedPoint
Projected coordinates and pointer to hit.
Definition: Cluster3D.h:352
std::list< ClusterParameters > ClusterParametersList
Definition: Cluster3D.h:404
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:97
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:247
attached to a cluster
Definition: Cluster3D.h:109
QTextStream & endl(QTextStream &s)
void setStatusBit(unsigned bits) const
Definition: Cluster3D.h:180