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