All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
MinSpanTreeAlg_tool.cc
Go to the documentation of this file.
1 /**
2  * @file MinSpanTreeAlg.cxx
3  *
4  * @brief Producer module to create 3D clusters from input hits
5  *
6  */
7 
8 // Framework Includes
12 #include "cetlib/cpu_timer.h"
13 #include "fhiclcpp/ParameterSet.h"
15 
16 // LArSoft includes
27 
28 // std includes
29 #include <iostream>
30 #include <memory>
31 #include <unordered_map>
32 
33 // Eigen includes
34 #include <Eigen/Core>
35 
36 #include "TVector3.h"
37 
38 //------------------------------------------------------------------------------------------------------------------------------------------
39 // implementation follows
40 
41 namespace lar_cluster3d {
42 
43 class MinSpanTreeAlg : virtual public IClusterAlg
44 {
45 public:
46  /**
47  * @brief Constructor
48  *
49  * @param pset
50  */
51  explicit MinSpanTreeAlg(const fhicl::ParameterSet&);
52 
53  /**
54  * @brief Destructor
55  */
57 
58  void configure(fhicl::ParameterSet const &pset) override;
59 
60  /**
61  * @brief Given a set of recob hits, run DBscan to form 3D clusters
62  *
63  * @param hitPairList The input list of 3D hits to run clustering on
64  * @param hitPairClusterMap A map of hits that have been clustered
65  * @param clusterParametersList A list of cluster objects (parameters from associated hits)
66  */
67  void Cluster3DHits(reco::HitPairList& hitPairList,
68  reco::ClusterParametersList& clusterParametersList) const override;
69 
71  reco::ClusterParametersList& clusterParametersList) const override {return;}
72 
73  /**
74  * @brief If monitoring, recover the time to execute a particular function
75  */
76  float getTimeToExecute(TimeValues index) const override {return m_timeVector.at(index);}
77 
78 private:
79 
80  /**
81  * @brief Driver for Prim's algorithm
82  */
84 
85  /**
86  * @brief Prune the obvious ambiguous hits
87  */
89 
90  /**
91  * @brief Algorithm to find the best path through the given cluster
92  */
94 
95  /**
96  * @brief a depth first search to find longest branches
97  */
99 
100  /**
101  * @brief Alternative version of FindBestPathInCluster utilizing an A* algorithm
102  */
104 
105  /**
106  * @brief Algorithm to find shortest path between two 3D hits
107  */
109 
110  using BestNodeTuple = std::tuple<const reco::ClusterHit3D*,float,float>;
111  using BestNodeMap = std::unordered_map<const reco::ClusterHit3D*,BestNodeTuple>;
112 
114 
116 
117  /**
118  * @brief Find the lowest cost path between two nodes using MST edges
119  */
120  void LeastCostPath(const reco::EdgeTuple&,
121  const reco::ClusterHit3D*,
123  float&) const;
124 
125  void CheckHitSorting(reco::ClusterParameters& clusterParams) const;
126 
127  /**
128  * @brief define data structure for keeping track of channel status
129  */
130  using ChannelStatusVec = std::vector<size_t>;
131  using ChannelStatusByPlaneVec = std::vector<ChannelStatusVec>;
132 
133  /**
134  * @brief Data members to follow
135  */
137  mutable std::vector<float> m_timeVector; ///<
138  std::vector<std::vector<float>> m_wireDir; ///<
139 
140  geo::Geometry const* m_geometry; //< pointer to the Geometry service
141 
142  PrincipalComponentsAlg m_pcaAlg; // For running Principal Components Analysis
143  kdTree m_kdTree; // For the kdTree
144 
145  std::unique_ptr<lar_cluster3d::IClusterParametersBuilder> m_clusterBuilder; ///< Common cluster builder tool
146 };
147 
149  m_pcaAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg")),
150  m_kdTree(pset.get<fhicl::ParameterSet>("kdTree"))
151 {
152  this->configure(pset);
153 }
154 
155 //------------------------------------------------------------------------------------------------------------------------------------------
156 
158 {
159 }
160 
161 //------------------------------------------------------------------------------------------------------------------------------------------
162 
164 {
165  m_enableMonitoring = pset.get<bool> ("EnableMonitoring", true );
166 
168 
169  m_geometry = &*geometry;
170 
171  m_timeVector.resize(NUMTIMEVALUES, 0.);
172 
173  // Determine the unit directon and normal vectors to the wires
174  m_wireDir.resize(3);
175 
176  raw::ChannelID_t uChannel(0);
177  std::vector<geo::WireID> uWireID = m_geometry->ChannelToWire(uChannel);
178  const geo::WireGeo* uWireGeo = m_geometry->WirePtr(uWireID[0]);
179 
180  TVector3 uWireDir = uWireGeo->Direction();
181 
182  m_wireDir[0].resize(3);
183  m_wireDir[0][0] = uWireDir[0];
184  m_wireDir[0][1] = -uWireDir[2];
185  m_wireDir[0][2] = uWireDir[1];
186 
187  raw::ChannelID_t vChannel(2400);
188  std::vector<geo::WireID> vWireID = m_geometry->ChannelToWire(vChannel);
189  const geo::WireGeo* vWireGeo = m_geometry->WirePtr(vWireID[0]);
190 
191  TVector3 vWireDir = vWireGeo->Direction();
192 
193  m_wireDir[1].resize(3);
194  m_wireDir[1][0] = vWireDir[0];
195  m_wireDir[1][1] = -vWireDir[2];
196  m_wireDir[1][2] = vWireDir[1];
197 
198  m_wireDir[2].resize(3);
199  m_wireDir[2][0] = 0.;
200  m_wireDir[2][1] = 0.;
201  m_wireDir[2][2] = 1.;
202 
203  m_clusterBuilder = art::make_tool<lar_cluster3d::IClusterParametersBuilder>(pset.get<fhicl::ParameterSet>("ClusterParamsBuilder"));
204 
205  return;
206 }
207 
209  reco::ClusterParametersList& clusterParametersList) const
210 {
211  /**
212  * @brief Driver for processing input 2D hits, transforming to 3D hits and building lists
213  * of associated 3D hits (candidate 3D clusters)
214  */
215 
216  // Zero the time vector
217  if (m_enableMonitoring) std::fill(m_timeVector.begin(),m_timeVector.end(),0.);
218 
219  // DBScan is driven of its "epsilon neighborhood". Computing adjacency within DBScan can be time
220  // consuming so the idea is the prebuild the adjaceny map and then run DBScan.
221  // The following call does this work
222  kdTree::KdTreeNodeList kdTreeNodeContainer;
223  kdTree::KdTreeNode topNode = m_kdTree.BuildKdTree(hitPairList, kdTreeNodeContainer);
224 
226 
227  // Run DBScan to get candidate clusters
228  RunPrimsAlgorithm(hitPairList, topNode, clusterParametersList);
229 
230  // Initial clustering is done, now trim the list and get output parameters
231  cet::cpu_timer theClockBuildClusters;
232 
233  // Start clocks if requested
234  if (m_enableMonitoring) theClockBuildClusters.start();
235 
236  m_clusterBuilder->BuildClusterInfo(clusterParametersList);
237 
238  if (m_enableMonitoring)
239  {
240  theClockBuildClusters.stop();
241 
242  m_timeVector[BUILDCLUSTERINFO] = theClockBuildClusters.accumulated_real_time();
243  }
244 
245  // Test run the path finding algorithm
246  for (auto& clusterParams : clusterParametersList) FindBestPathInCluster(clusterParams, topNode);
247 
248  mf::LogDebug("MinSpanTreeAlg") << ">>>>> Cluster3DHits done, found " << clusterParametersList.size() << " clusters" << std::endl;
249 
250  return;
251 }
252 
253 //------------------------------------------------------------------------------------------------------------------------------------------
255  kdTree::KdTreeNode& topNode,
256  reco::ClusterParametersList& clusterParametersList) const
257 {
258  // If no hits then no work
259  if (hitPairList.empty()) return;
260 
261  // Now proceed with building the clusters
262  cet::cpu_timer theClockDBScan;
263 
264  // Start clocks if requested
265  if (m_enableMonitoring) theClockDBScan.start();
266 
267  // Initialization
268  size_t clusterIdx(0);
269 
270  // This will contain our list of edges
271  reco::EdgeList curEdgeList;
272 
273  // Get the first point
274  reco::HitPairList::iterator freeHitItr = hitPairList.begin();
275  const reco::ClusterHit3D* lastAddedHit = &(*freeHitItr++);
276 
278 
279  // Make a cluster...
280  clusterParametersList.push_back(reco::ClusterParameters());
281 
282  // Get an iterator to the first cluster
283  reco::ClusterParametersList::iterator curClusterItr = --clusterParametersList.end();
284 
285  // We use pointers here because the objects they point to will change in the loop below
286  reco::Hit3DToEdgeMap* curEdgeMap = &(*curClusterItr).getHit3DToEdgeMap();
287  reco::HitPairListPtr* curCluster = &(*curClusterItr).getHitPairListPtr();
288 
289  // Loop until all hits have been associated to a cluster
290  while(1)
291  {
292  // and the 3D hit status bits
294 
295  // Purge the current list to get rid of edges which point to hits already in the cluster
296  for(reco::EdgeList::iterator curEdgeItr = curEdgeList.begin(); curEdgeItr != curEdgeList.end();)
297  {
298  if (std::get<1>(*curEdgeItr)->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED)
299  curEdgeItr = curEdgeList.erase(curEdgeItr);
300  else curEdgeItr++;
301  }
302 
303  // Add the lastUsedHit to the current cluster
304  curCluster->push_back(lastAddedHit);
305 
306  // Set up to find the list of nearest neighbors to the last used hit...
307  kdTree::CandPairList CandPairList;
308  float bestDistance(1.5); //std::numeric_limits<float>::max());
309 
310  // And find them... result will be an unordered list of neigbors
311  m_kdTree.FindNearestNeighbors(lastAddedHit, topNode, CandPairList, bestDistance);
312 
313  // Copy edges to the current list (but only for hits not already in a cluster)
314 // for(auto& pair : CandPairList)
315 // if (!(pair.second->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED)) curEdgeList.push_back(reco::EdgeTuple(lastAddedHit,pair.second,pair.first));
316  for(auto& pair : CandPairList)
317  {
318  if (!(pair.second->getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED))
319  {
320  double edgeWeight = lastAddedHit->getHitChiSquare() * pair.second->getHitChiSquare();
321 
322  curEdgeList.push_back(reco::EdgeTuple(lastAddedHit,pair.second,edgeWeight));
323  }
324  }
325 
326  // If the edge list is empty then we have a complete cluster
327  if (curEdgeList.empty())
328  {
329  std::cout << "-----------------------------------------------------------------------------------------" << std::endl;
330  std::cout << "**> Cluster idx: " << clusterIdx++ << " has " << curCluster->size() << " hits" << std::endl;
331 
332  // Look for the next "free" hit
333  freeHitItr = std::find_if(freeHitItr,hitPairList.end(),[](const auto& hit){return !(hit.getStatusBits() & reco::ClusterHit3D::CLUSTERATTACHED);});
334 
335  // If at end of input list we are done with all hits
336  if (freeHitItr == hitPairList.end()) break;
337 
338  std::cout << "##################################################################>Processing another cluster" << std::endl;
339 
340  // Otherwise, get a new cluster and set up
341  clusterParametersList.push_back(reco::ClusterParameters());
342 
343  curClusterItr = --clusterParametersList.end();
344 
345  curEdgeMap = &(*curClusterItr).getHit3DToEdgeMap();
346  curCluster = &(*curClusterItr).getHitPairListPtr();
347  lastAddedHit = &(*freeHitItr++);
348  }
349  // Otherwise we are still processing the current cluster
350  else
351  {
352  // Sort the list of edges by distance
353  curEdgeList.sort([](const auto& left,const auto& right){return std::get<2>(left) < std::get<2>(right);});
354 
355  // Populate the map with the edges...
356  reco::EdgeTuple& curEdge = curEdgeList.front();
357 
358  (*curEdgeMap)[std::get<0>(curEdge)].push_back(curEdge);
359  (*curEdgeMap)[std::get<1>(curEdge)].push_back(reco::EdgeTuple(std::get<1>(curEdge),std::get<0>(curEdge),std::get<2>(curEdge)));
360 
361  // Update the last hit to be added to the collection
362  lastAddedHit = std::get<1>(curEdge);
363  }
364  }
365 
366  if (m_enableMonitoring)
367  {
368  theClockDBScan.stop();
369 
370  m_timeVector[RUNDBSCAN] = theClockDBScan.accumulated_real_time();
371  }
372 
373  return;
374 }
375 
377 {
378  reco::HitPairListPtr longestCluster;
379  float bestQuality(0.);
380  float aveNumEdges(0.);
381  size_t maxNumEdges(0);
382  size_t nIsolatedHits(0);
383 
384  // Now proceed with building the clusters
385  cet::cpu_timer theClockPathFinding;
386 
387  // Start clocks if requested
388  if (m_enableMonitoring) theClockPathFinding.start();
389 
390  reco::HitPairListPtr& hitPairList = curCluster.getHitPairListPtr();
391  reco::Hit3DToEdgeMap& curEdgeMap = curCluster.getHit3DToEdgeMap();
392  reco::EdgeList& bestEdgeList = curCluster.getBestEdgeList();
393 
394  // Do some spelunking...
395  for(const auto& hit : hitPairList)
396  {
397  if (!curEdgeMap[hit].empty() && curEdgeMap[hit].size() == 1)
398  {
399  float quality(0.);
400 
401  reco::HitPairListPtr tempList = DepthFirstSearch(curEdgeMap[hit].front(), curEdgeMap, quality);
402 
403  tempList.push_front(std::get<0>(curEdgeMap[hit].front()));
404 
405  if (quality > bestQuality)
406  {
407  longestCluster = tempList;
408  bestQuality = quality;
409  }
410 
411  nIsolatedHits++;
412  }
413 
414  aveNumEdges += float(curEdgeMap[hit].size());
415  maxNumEdges = std::max(maxNumEdges,curEdgeMap[hit].size());
416  }
417 
418  aveNumEdges /= float(hitPairList.size());
419  std::cout << "----> # isolated hits: " << nIsolatedHits << ", longest branch: " << longestCluster.size() << ", cluster size: " << hitPairList.size() << ", ave # edges: " << aveNumEdges << ", max: " << maxNumEdges << std::endl;
420 
421  if (!longestCluster.empty())
422  {
423  hitPairList = longestCluster;
424  for(const auto& hit : hitPairList)
425  {
426  for(const auto& edge : curEdgeMap[hit]) bestEdgeList.emplace_back(edge);
427  }
428 
429  std::cout << " ====> new cluster size: " << hitPairList.size() << std::endl;
430  }
431 
432  if (m_enableMonitoring)
433  {
434  theClockPathFinding.stop();
435 
436  m_timeVector[PATHFINDING] += theClockPathFinding.accumulated_real_time();
437  }
438 
439  return;
440 }
441 
443 {
444  // Set up for timing the function
445  cet::cpu_timer theClockPathFinding;
446 
447  // Start clocks if requested
448  if (m_enableMonitoring) theClockPathFinding.start();
449 
450  // Trial A* here
451  if (clusterParams.getHitPairListPtr().size() > 2)
452  {
453  // Get references to what we need....
454  reco::HitPairListPtr& curCluster = clusterParams.getHitPairListPtr();
455  reco::Hit3DToEdgeMap& curEdgeMap = clusterParams.getHit3DToEdgeMap();
456 
457  // Do a quick PCA to determine our parameter "alpha"
459  m_pcaAlg.PCAAnalysis_3D(curCluster, pca);
460 
461  // The chances of a failure are remote, still we should check
462  if (pca.getSvdOK())
463  {
464  float pcaLen = 3.0*sqrt(pca.getEigenValues()[2]);
465  float pcaWidth = 3.0*sqrt(pca.getEigenValues()[1]);
466  float pcaHeight = 3.0*sqrt(pca.getEigenValues()[0]);
467  const Eigen::Vector3f& pcaCenter = pca.getAvePosition();
468  float alpha = std::min(float(1.),std::max(float(0.001),pcaWidth/pcaLen));
469 
470  // Create a temporary container for the isolated points
471  reco::ProjectedPointList isolatedPointList;
472 
473  // Go through and find the isolated points, for those get the projection to the plane of maximum spread
474  for(const auto& hit3D : curCluster)
475  {
476  // the definition of an isolated hit is that it only has one associated edge
477  if (!curEdgeMap[hit3D].empty() && curEdgeMap[hit3D].size() == 1)
478  {
479  Eigen::Vector3f pcaToHitVec(hit3D->getPosition()[0] - pcaCenter(0),
480  hit3D->getPosition()[1] - pcaCenter(1),
481  hit3D->getPosition()[2] - pcaCenter(2));
482  Eigen::Vector3f pcaToHit = pca.getEigenVectors() * pcaToHitVec;
483 
484  // This sets x,y where x is the longer spread, y the shorter
485  isolatedPointList.emplace_back(pcaToHit(2),pcaToHit(1),hit3D);
486  }
487  }
488 
489  std::cout << "************* Finding best path with A* in cluster *****************" << std::endl;
490  std::cout << "**> There are " << curCluster.size() << " hits, " << isolatedPointList.size() << " isolated hits, the alpha parameter is " << alpha << std::endl;
491  std::cout << "**> PCA len: " << pcaLen << ", wid: " << pcaWidth << ", height: " << pcaHeight << ", ratio: " << pcaHeight/pcaWidth << std::endl;
492 
493  // If no isolated points then nothing to do...
494  if (isolatedPointList.size() > 1)
495  {
496  // Sort the point vec by increasing x, if same then by increasing y.
497  isolatedPointList.sort([](const auto& left, const auto& right){return (std::abs(std::get<0>(left) - std::get<0>(right)) > std::numeric_limits<float>::epsilon()) ? std::get<0>(left) < std::get<0>(right) : std::get<1>(left) < std::get<1>(right);});
498 
499  // Ok, get the two most distance points...
500  const reco::ClusterHit3D* startHit = std::get<2>(isolatedPointList.front());
501  const reco::ClusterHit3D* stopHit = std::get<2>(isolatedPointList.back());
502 
503  std::cout << "**> Sorted " << isolatedPointList.size() << " hits, longest distance: " << DistanceBetweenNodes(startHit,stopHit) << std::endl;
504 
505  // Call the AStar function to try to find the best path...
506 // AStar(startHit,stopHit,alpha,topNode,clusterParams);
507 
508  float cost(std::numeric_limits<float>::max());
509 
510  LeastCostPath(curEdgeMap[startHit].front(),stopHit,clusterParams,cost);
511 
512  clusterParams.getBestHitPairListPtr().push_front(startHit);
513 
514  std::cout << "**> Best path has " << clusterParams.getBestHitPairListPtr().size() << " hits, " << clusterParams.getBestEdgeList().size() << " edges" << std::endl;
515  }
516  }
517  else
518  {
519  std::cout << "++++++>>> PCA failure! # hits: " << curCluster.size() << std::endl;
520  }
521  }
522 
523  if (m_enableMonitoring)
524  {
525  theClockPathFinding.stop();
526 
527  m_timeVector[PATHFINDING] += theClockPathFinding.accumulated_real_time();
528  }
529 
530  return;
531 }
532 
534  const reco::ClusterHit3D* goalNode,
535  float alpha,
536  kdTree::KdTreeNode& topNode,
537  reco::ClusterParameters& clusterParams) const
538 {
539  // Recover the list of hits and edges
540  reco::HitPairListPtr& pathNodeList = clusterParams.getBestHitPairListPtr();
541  reco::EdgeList& bestEdgeList = clusterParams.getBestEdgeList();
542  reco::Hit3DToEdgeMap& curEdgeMap = clusterParams.getHit3DToEdgeMap();
543 
544  // Find the shortest path from start to goal using an A* algorithm
545  // Keep track of the nodes which have already been evaluated
546  reco::HitPairListPtr closedList;
547 
548  // Keep track of the nodes that have been "discovered" but yet to be evaluated
549  reco::HitPairListPtr openList = {startNode};
550 
551  // Create a map which, for each node, will tell us the node it can be most efficiencly reached from.
552  BestNodeMap bestNodeMap;
553 
554  bestNodeMap[startNode] = BestNodeTuple(startNode,0.,DistanceBetweenNodes(startNode,goalNode));
555 
556  alpha = 1.; //std::max(0.5,alpha);
557 
558  while(!openList.empty())
559  {
560  // The list is not empty so by def we will return something
561  reco::HitPairListPtr::iterator currentNodeItr = openList.begin();
562 
563  // If the list contains more than one element then we need to find the one with the smallest total estimated cost to the end
564  if (openList.size() > 1)
565  currentNodeItr = std::min_element(openList.begin(),openList.end(),[bestNodeMap](const auto& next, const auto& best){return std::get<2>(bestNodeMap.at(next)) < std::get<2>(bestNodeMap.at(best));});
566 
567  // Easier to deal directly with the pointer to the node
568  const reco::ClusterHit3D* currentNode = *currentNodeItr;
569 
570  // Check to see if we have reached the goal and need to evaluate the path
571  if (currentNode == goalNode)
572  {
573  // The path reconstruction will
574  ReconstructBestPath(goalNode, bestNodeMap, pathNodeList, bestEdgeList);
575 
576  break;
577  }
578 
579  // Otherwise need to keep evaluating
580  else
581  {
582  openList.erase(currentNodeItr);
584 
585  // Get tuple values for the current node
586  const BestNodeTuple& currentNodeTuple = bestNodeMap.at(currentNode);
587  float currentNodeScore = std::get<1>(currentNodeTuple);
588 
589  // Recover the edges associated to the current point
590  const reco::EdgeList& curEdgeList = curEdgeMap[currentNode];
591 
592  for(const auto& curEdge : curEdgeList)
593  {
594  const reco::ClusterHit3D* candHit3D = std::get<1>(curEdge);
595 
596  if (candHit3D->getStatusBits() & reco::ClusterHit3D::PATHCHECKED) continue;
597 
598  float tentative_gScore = currentNodeScore + std::get<2>(curEdge);
599 
600  // Have we seen the candidate node before?
601  BestNodeMap::iterator candNodeItr = bestNodeMap.find(candHit3D);
602 
603  if (candNodeItr == bestNodeMap.end())
604  {
605  openList.push_back(candHit3D);
606  }
607  else if (tentative_gScore > std::get<1>(candNodeItr->second)) continue;
608 
609  // Make a guess at score to get to target...
610  float guessToTarget = DistanceBetweenNodes(candHit3D,goalNode) / 0.3;
611 
612  bestNodeMap[candHit3D] = BestNodeTuple(currentNode,tentative_gScore, tentative_gScore + guessToTarget);
613  }
614  }
615  }
616 
617  return;
618 }
619 
621  BestNodeMap& bestNodeMap,
622  reco::HitPairListPtr& pathNodeList,
623  reco::EdgeList& bestEdgeList) const
624 {
625  while(std::get<0>(bestNodeMap.at(goalNode)) != goalNode)
626  {
627  const reco::ClusterHit3D* nextNode = std::get<0>(bestNodeMap[goalNode]);
628  reco::EdgeTuple bestEdge = reco::EdgeTuple(goalNode,nextNode,DistanceBetweenNodes(goalNode,nextNode));
629 
630  pathNodeList.push_front(goalNode);
631  bestEdgeList.push_front(bestEdge);
632 
633  goalNode = nextNode;
634  }
635 
636  pathNodeList.push_front(goalNode);
637 
638  return;
639 }
640 
642  const reco::ClusterHit3D* goalNode,
643  reco::ClusterParameters& clusterParams,
644  float& showMeTheMoney) const
645 {
646  // Recover the mapping between hits and edges
647  reco::Hit3DToEdgeMap& curEdgeMap = clusterParams.getHit3DToEdgeMap();
648 
649  reco::Hit3DToEdgeMap::const_iterator edgeListItr = curEdgeMap.find(std::get<1>(curEdge));
650 
651  showMeTheMoney = std::numeric_limits<float>::max();
652 
653  if (edgeListItr != curEdgeMap.end() && !edgeListItr->second.empty())
654  {
655  reco::HitPairListPtr& bestNodeList = clusterParams.getBestHitPairListPtr();
656  reco::EdgeList& bestEdgeList = clusterParams.getBestEdgeList();
657 
658  for(const auto& edge : edgeListItr->second)
659  {
660  // skip the self reference
661  if (std::get<1>(edge) == std::get<0>(curEdge)) continue;
662 
663  // Have we found the droid we are looking for?
664  if (std::get<1>(edge) == goalNode)
665  {
666  bestNodeList.push_back(goalNode);
667  bestEdgeList.push_back(edge);
668  showMeTheMoney = std::get<2>(edge);
669  break;
670  }
671 
672  // Keep searching, it is out there somewhere...
673  float currentCost(0.);
674 
675  LeastCostPath(edge,goalNode,clusterParams,currentCost);
676 
677  if (currentCost < std::numeric_limits<float>::max())
678  {
679  showMeTheMoney = std::get<2>(edge) + currentCost;
680  break;
681  }
682  }
683  }
684 
685  if (showMeTheMoney < std::numeric_limits<float>::max())
686  {
687  clusterParams.getBestHitPairListPtr().push_front(std::get<1>(curEdge));
688  clusterParams.getBestEdgeList().push_front(curEdge);
689  }
690 
691  return;
692 }
693 
695 {
696  const Eigen::Vector3f& node1Pos = node1->getPosition();
697  const Eigen::Vector3f& node2Pos = node2->getPosition();
698  float deltaNode[] = {node1Pos[0]-node2Pos[0], node1Pos[1]-node2Pos[1], node1Pos[2]-node2Pos[2]};
699 
700  // Standard euclidean distance
701  return std::sqrt(deltaNode[0]*deltaNode[0]+deltaNode[1]*deltaNode[1]+deltaNode[2]*deltaNode[2]);
702 
703  // Manhatten distance
704  //return std::fabs(deltaNode[0]) + std::fabs(deltaNode[1]) + std::fabs(deltaNode[2]);
705 /*
706  // Chebyshev distance
707  // We look for maximum distance by wires...
708 
709  // Now go through the hits and compare view by view to look for delta wire and tigher constraint on delta t
710  int wireDeltas[] = {0,0,0};
711 
712  for(size_t idx = 0; idx < 3; idx++)
713  wireDeltas[idx] = std::abs(int(node1->getWireIDs()[idx].Wire) - int(node2->getWireIDs()[idx].Wire));
714 
715  // put wire deltas in order...
716  std::sort(wireDeltas, wireDeltas + 3);
717 
718  return std::sqrt(deltaNode[0]*deltaNode[0] + 0.09 * float(wireDeltas[2]*wireDeltas[2]));
719  */
720 }
721 
723  const reco::Hit3DToEdgeMap& hitToEdgeMap,
724  float& bestTreeQuality) const
725 {
726  reco::HitPairListPtr hitPairListPtr;
727  float bestQuality(0.);
728  float curEdgeWeight = std::max(0.3,std::get<2>(curEdge));
729  float curEdgeProj(1./curEdgeWeight);
730 
731  reco::Hit3DToEdgeMap::const_iterator edgeListItr = hitToEdgeMap.find(std::get<1>(curEdge));
732 
733  if (edgeListItr != hitToEdgeMap.end())
734  {
735  // The input edge weight has quality factors applied, recalculate just the position difference
736  const Eigen::Vector3f& firstHitPos = std::get<0>(curEdge)->getPosition();
737  const Eigen::Vector3f& secondHitPos = std::get<1>(curEdge)->getPosition();
738  float curEdgeVec[] = {secondHitPos[0]-firstHitPos[0],secondHitPos[1]-firstHitPos[1],secondHitPos[2]-firstHitPos[2]};
739  float curEdgeMag = std::sqrt(curEdgeVec[0]*curEdgeVec[0]+curEdgeVec[1]*curEdgeVec[1]+curEdgeVec[2]*curEdgeVec[2]);
740 
741  curEdgeMag = std::max(float(0.1),curEdgeMag);
742 
743  for(const auto& edge : edgeListItr->second)
744  {
745  // skip the self reference
746  if (std::get<1>(edge) == std::get<0>(curEdge)) continue;
747 
748  float quality(0.);
749 
750  reco::HitPairListPtr tempList = DepthFirstSearch(edge,hitToEdgeMap,quality);
751 
752  if (quality > bestQuality)
753  {
754  hitPairListPtr = tempList;
755  bestQuality = quality;
756  curEdgeProj = 1./curEdgeMag;
757  }
758  }
759  }
760 
761  hitPairListPtr.push_front(std::get<1>(curEdge));
762 
763  bestTreeQuality += bestQuality + curEdgeProj;
764 
765  return hitPairListPtr;
766 }
767 
769 {
770 
771  // Recover the HitPairListPtr from the input clusterParams (which will be the
772  // only thing that has been provided)
773  reco::HitPairListPtr& hitPairVector = clusterParams.getHitPairListPtr();
774 
775  size_t nStartedWith(hitPairVector.size());
776  size_t nRejectedHits(0);
777 
778  reco::HitPairListPtr goodHits;
779 
780  // Loop through the hits and try to week out the clearly ambiguous ones
781  for(const auto& hit3D : hitPairVector)
782  {
783  // Loop to try to remove ambiguous hits
784  size_t n2DHitsIn3DHit(0);
785  size_t nThisClusterOnly(0);
786  size_t nOtherCluster(0);
787 
788  // reco::ClusterParameters* otherCluster;
789  const std::set<const reco::ClusterHit3D*>* otherClusterHits = 0;
790 
791  for(const auto& hit2D : hit3D->getHits())
792  {
793  if (!hit2D) continue;
794 
795  n2DHitsIn3DHit++;
796 
797  if (hit2DToClusterMap[hit2D].size() < 2) nThisClusterOnly = hit2DToClusterMap[hit2D][&clusterParams].size();
798  else
799  {
800  for(const auto& clusterHitMap : hit2DToClusterMap[hit2D])
801  {
802  if (clusterHitMap.first == &clusterParams) continue;
803 
804  if (clusterHitMap.second.size() > nOtherCluster)
805  {
806  nOtherCluster = clusterHitMap.second.size();
807  // otherCluster = clusterHitMap.first;
808  otherClusterHits = &clusterHitMap.second;
809  }
810  }
811  }
812  }
813 
814  if (n2DHitsIn3DHit < 3 && nThisClusterOnly > 1 && nOtherCluster > 0)
815  {
816  bool skip3DHit(false);
817 
818  for(const auto& otherHit3D : *otherClusterHits)
819  {
820  size_t nOther2DHits(0);
821 
822  for(const auto& otherHit2D : otherHit3D->getHits())
823  {
824  if (!otherHit2D) continue;
825 
826  nOther2DHits++;
827  }
828 
829  if (nOther2DHits > 2)
830  {
831  skip3DHit = true;
832  nRejectedHits++;
833  break;
834  }
835  }
836 
837  if (skip3DHit) continue;
838 
839  }
840 
841  goodHits.emplace_back(hit3D);
842  }
843 
844  std::cout << "###>> Input " << nStartedWith << " hits, rejected: " << nRejectedHits << std::endl;
845 
846  hitPairVector.resize(goodHits.size());
847  std::copy(goodHits.begin(),goodHits.end(),hitPairVector.begin());
848 
849  return;
850 }
851 
852 struct HitPairClusterOrder
853 {
855  {
856  // Watch out for the case where two clusters can have the same number of hits!
857  return (*left).getHitPairListPtr().size() > (*right).getHitPairListPtr().size();
858  }
859 };
860 
862 {
863 public:
864  SetCheckHitOrder(const std::vector<size_t>& plane) : m_plane(plane) {}
865 
867  {
868  // Check if primary view's hit is on the same wire
869  if (left->getWireIDs()[m_plane[0]] == right->getWireIDs()[m_plane[0]])
870  {
871  // Same wire but not same hit, order by primary hit time
872  if (left->getHits()[m_plane[0]] && right->getHits()[m_plane[0]] && left->getHits()[m_plane[0]] != right->getHits()[m_plane[0]])
873  {
874  return left->getHits()[m_plane[0]]->getHit()->PeakTime() < right->getHits()[m_plane[0]]->getHit()->PeakTime();
875  }
876 
877  // Primary view is same hit, look at next view's wire
878  if (left->getWireIDs()[m_plane[1]] == right->getWireIDs()[m_plane[1]])
879  {
880  // Same wire but not same hit, order by secondary hit time
881  if (left->getHits()[m_plane[1]] && right->getHits()[m_plane[1]] && left->getHits()[m_plane[1]] != right->getHits()[m_plane[1]])
882  {
883  return left->getHits()[m_plane[1]]->getHit()->PeakTime() < right->getHits()[m_plane[1]]->getHit()->PeakTime();
884  }
885 
886  // All that is left is the final view... and this can't be the same hit... (else it is the same 3D hit)
887  return left->getWireIDs()[m_plane[2]] < right->getWireIDs()[m_plane[2]];
888  }
889 
890  return left->getWireIDs()[m_plane[1]] < right->getWireIDs()[m_plane[1]];
891  }
892 
893  // Order by primary view's wire number
894  return left->getWireIDs()[m_plane[0]] < right->getWireIDs()[m_plane[0]];
895  }
896 
897 private:
898  const std::vector<size_t>& m_plane;
899 };
900 
902 {
903  reco::HitPairListPtr& curCluster = clusterParams.getHitPairListPtr();
904 
905  // Trial A* here
906  if (curCluster.size() > 2)
907  {
908  // Do a quick PCA to determine our parameter "alpha"
910  m_pcaAlg.PCAAnalysis_3D(curCluster, pca);
911 
912  if (pca.getSvdOK())
913  {
914  const Eigen::Vector3f& pcaAxis = pca.getEigenVectors().row(2);
915 
916  std::vector<size_t> closestPlane = {0, 0, 0 };
917  std::vector<float> bestAngle = {0.,0.,0.};
918 
919  for(size_t plane = 0; plane < 3; plane++)
920  {
921  const std::vector<float>& wireDir = m_wireDir[plane];
922 
923  float dotProd = std::fabs(pcaAxis[0]*wireDir[0] + pcaAxis[1]*wireDir[1] + pcaAxis[2]*wireDir[2]);
924 
925  if (dotProd > bestAngle[0])
926  {
927  bestAngle[2] = bestAngle[1];
928  closestPlane[2] = closestPlane[1];
929  bestAngle[1] = bestAngle[0];
930  closestPlane[1] = closestPlane[0];
931  closestPlane[0] = plane;
932  bestAngle[0] = dotProd;
933  }
934  else if (dotProd > bestAngle[1])
935  {
936  bestAngle[2] = bestAngle[1];
937  closestPlane[2] = closestPlane[1];
938  closestPlane[1] = plane;
939  bestAngle[1] = dotProd;
940  }
941  else
942  {
943  closestPlane[2] = plane;
944  bestAngle[2] = dotProd;
945  }
946  }
947 
948  // Get a copy of our 3D hits
949  reco::HitPairListPtr localHitList = curCluster;
950 
951  // Sort the hits
952  localHitList.sort(SetCheckHitOrder(closestPlane));
953 
954  // Ok, let's print it all and take a look
955  std::cout << "********************************************************************************************" << std::endl;
956  std::cout << "**>>>>> longest axis: " << closestPlane[0] << ", best angle: " << bestAngle[0] << std::endl;
957  std::cout << "**>>>>> second axis: " << closestPlane[1] << ", best angle: " << bestAngle[1] << std::endl;
958  std::cout << " " << std::endl;
959 
960  reco::HitPairListPtr::iterator firstHitItr = localHitList.begin();
961  reco::HitPairListPtr::iterator lastHitItr = localHitList.begin();
962 
963  size_t bestPlane = closestPlane[0];
964 
965  reco::HitPairListPtr testList;
966 
967  while(firstHitItr != localHitList.end())
968  {
969  const reco::ClusterHit3D* currentHit = *firstHitItr;
970 
971  // Search for the last matching best view hit
972  while(lastHitItr != localHitList.end())
973  {
974  // If a different wire on the best view then we're certainly done
975  if (currentHit->getWireIDs()[bestPlane] != (*lastHitItr)->getWireIDs()[bestPlane]) break;
976 
977  // More subtle test to see if same wire but different hit (being careful of case of no hit)
978  if (currentHit->getHits()[bestPlane] && (*lastHitItr)->getHits()[bestPlane] && currentHit->getHits()[bestPlane] != (*lastHitItr)->getHits()[bestPlane]) break;
979 
980  // Yet event more subtle test...
981  if ((!(currentHit->getHits()[bestPlane]) && (*lastHitItr)->getHits()[bestPlane]) || (currentHit->getHits()[bestPlane] && !((*lastHitItr)->getHits()[bestPlane]))) break;
982 
983  // Not there yet...
984  lastHitItr++;
985  }
986 
987  // How many hits in this chain?
988 // size_t numHits(std::distance(firstHitItr,lastHitItr));
989 // float minOverlapFraction(0.);
990 
991 // if (numHits > 1)
992 // {
993 // reco::HitPairListPtr::iterator bestMinOverlapItr = std::max_element(firstHitItr,lastHitItr,[](const auto& left, const auto& right){return left->getMinOverlapFraction() < right->getMinOverlapFraction();});
994 //
995 // minOverlapFraction = std::min(0.999*(*bestMinOverlapItr)->getMinOverlapFraction(),0.90);
996 // }
997 
998  while(firstHitItr != lastHitItr)
999  {
1000 // if (currentHit->getMinOverlapFraction() > minOverlapFraction) testList.push_back(currentHit); //currentHit->setStatusBit(reco::ClusterHit3D::SKELETONHIT);
1001 
1002  currentHit = *++firstHitItr;
1003  }
1004 
1005  firstHitItr = lastHitItr;
1006  }
1007 /*
1008  for(const auto& hit : localHitList)
1009  {
1010  std::cout << "- wires: ";
1011 
1012  for(size_t idx = 0; idx < 3; idx++)
1013  {
1014  float viewTime = -1.;
1015 
1016  if (hit->getHits()[closestView[idx]]) viewTime = hit->getHits()[closestView[idx]]->getTimeTicks();
1017 
1018  std::cout << closestView[idx] << ":" << hit->getWireIDs()[closestView[idx]].Wire << " - " << viewTime << ", ";
1019  }
1020 
1021  bool isSkeleton = hit->getStatusBits() & reco::ClusterHit3D::SKELETONHIT;
1022 
1023  std::cout << "ave time: " << hit->getAvePeakTime() << ", min/max overlap: " << hit->getMinOverlapFraction() << "/" << hit->getMaxOverlapFraction() << ", tagged: " << isSkeleton << std::endl;
1024 
1025  if (isSkeleton) testList.push_back(hit);
1026  }
1027 */
1028  curCluster = testList;
1029  }
1030  }
1031 
1032  return;
1033 }
1034 
1036 } // namespace lar_cluster3d
bool operator()(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right) const
reco::HitPairListPtr & getBestHitPairListPtr()
Definition: Cluster3D.h:480
intermediate_table::iterator iterator
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
void configure(fhicl::ParameterSet const &pset) override
std::list< reco::ClusterHit3D > HitPairList
Definition: Cluster3D.h:339
bool getSvdOK() const
Definition: Cluster3D.h:244
std::vector< std::vector< float > > m_wireDir
#define DEFINE_ART_CLASS_TOOL(tool)
Definition: ToolMacros.h:42
void LeastCostPath(const reco::EdgeTuple &, const reco::ClusterHit3D *, reco::ClusterParameters &, float &) const
Find the lowest cost path between two nodes using MST edges.
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
float getTimeToExecute(TimeValues index) const override
If monitoring, recover the time to execute a particular function.
std::tuple< const reco::ClusterHit3D *, float, float > BestNodeTuple
kdTree class definiton
Definition: kdTree.h:30
bool m_enableMonitoring
Data members to follow.
const Eigen::Vector3f getPosition() const
Definition: Cluster3D.h:158
std::list< ProjectedPoint > ProjectedPointList
Definition: Cluster3D.h:353
void Cluster3DHits(reco::HitPairListPtr &hitPairList, reco::ClusterParametersList &clusterParametersList) const override
Given a set of recob hits, run DBscan to form 3D clusters.
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
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
const std::vector< size_t > & m_plane
reco::EdgeList & getBestEdgeList()
Definition: Cluster3D.h:481
MinSpanTreeAlg(const fhicl::ParameterSet &)
Constructor.
reco::Hit3DToEdgeMap & getHit3DToEdgeMap()
Definition: Cluster3D.h:479
intermediate_table::const_iterator const_iterator
reco::HitPairListPtr & getHitPairListPtr()
Definition: Cluster3D.h:476
IClusterAlg interface class definiton.
Definition: IClusterAlg.h:25
std::vector< size_t > ChannelStatusVec
define data structure for keeping track of channel status
float DistanceBetweenNodes(const reco::ClusterHit3D *, const reco::ClusterHit3D *) const
unsigned int getStatusBits() const
Definition: Cluster3D.h:157
void AStar(const reco::ClusterHit3D *, const reco::ClusterHit3D *, float alpha, kdTree::KdTreeNode &, reco::ClusterParameters &) const
Algorithm to find shortest path between two 3D hits.
art framework interface to geometry description
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
T abs(T value)
std::list< EdgeTuple > EdgeList
Definition: Cluster3D.h:345
define a kd tree node
Definition: kdTree.h:118
void PruneAmbiguousHits(reco::ClusterParameters &, reco::Hit2DToClusterMap &) const
Prune the obvious ambiguous hits.
fInnerVessel push_back(Point(-578.400000, 0.000000, 0.000000))
const EigenValues & getEigenValues() const
Definition: Cluster3D.h:246
std::vector< ChannelStatusVec > ChannelStatusByPlaneVec
std::unordered_map< const reco::ClusterHit3D *, BestNodeTuple > BestNodeMap
T get(std::string const &key) const
Definition: ParameterSet.h:271
const Eigen::Vector3f & getAvePosition() const
Definition: Cluster3D.h:248
void RunPrimsAlgorithm(reco::HitPairList &, kdTree::KdTreeNode &, reco::ClusterParametersList &) const
Driver for Prim&#39;s algorithm.
double alpha
Definition: doAna.cpp:15
Path checking algorithm has seen this hit.
Definition: Cluster3D.h:111
TimeValues
enumerate the possible values for time checking if monitoring timing
Definition: IClusterAlg.h:61
void CheckHitSorting(reco::ClusterParameters &clusterParams) const
void Cluster3DHits(reco::HitPairList &hitPairList, reco::ClusterParametersList &clusterParametersList) const override
Given a set of recob hits, run DBscan to form 3D clusters.
std::tuple< const reco::ClusterHit3D *, const reco::ClusterHit3D *, double > EdgeTuple
Definition: Cluster3D.h:344
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:335
This provides an art tool interface definition for 3D Cluster algorithms.
static int max(int a, int b)
The geometry of one entire detector, as served by art.
Definition: Geometry.h:196
Definition of data types for geometry description.
std::unique_ptr< lar_cluster3d::IClusterParametersBuilder > m_clusterBuilder
Common cluster builder tool.
Vector Direction() const
Definition: WireGeo.h:587
Detector simulation of raw signals on wires.
std::list< CandPair > CandPairList
Definition: kdTree.h:79
Encapsulate the geometry of a wire.
Declaration of signal hit object.
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
def fill(s)
Definition: translator.py:93
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
const std::vector< geo::WireID > & getWireIDs() const
Definition: Cluster3D.h:173
T copy(T const &v)
std::unordered_map< const reco::ClusterHit2D *, ClusterToHitPairSetMap > Hit2DToClusterMap
Definition: Cluster3D.h:511
bool operator()(const reco::ClusterParametersList::iterator &left, const reco::ClusterParametersList::iterator &right)
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
unsigned int ChannelID_t
Type representing the ID of a readout channel.
Definition: RawTypes.h:28
float getTimeToExecute() const
Definition: kdTree.h:95
auto const & get(AssnsNode< L, R, D > const &r)
Definition: AssnsNode.h:115
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
const ClusterHit2DVec & getHits() const
Definition: Cluster3D.h:171
void start()
Definition: cpu_timer.cc:83
void ReconstructBestPath(const reco::ClusterHit3D *, BestNodeMap &, reco::HitPairListPtr &, reco::EdgeList &) const
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)
SetCheckHitOrder(const std::vector< size_t > &plane)
void FindBestPathInCluster(reco::ClusterParameters &) const
Algorithm to find the best path through the given cluster.
void setStatusBit(unsigned bits) const
Definition: Cluster3D.h:180
WireGeo const * WirePtr(geo::WireID const &wireid) const
Returns the specified wire.
reco::HitPairListPtr DepthFirstSearch(const reco::EdgeTuple &, const reco::Hit3DToEdgeMap &, float &) const
a depth first search to find longest branches