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

Public Member Functions

 ClusterPathFinder (const fhicl::ParameterSet &)
 Constructor. More...
 
 ~ClusterPathFinder ()
 Destructor. More...
 
void configure (fhicl::ParameterSet const &pset) override
 
void initializeHistograms (art::TFileDirectory &) override
 Interface for initializing histograms if they are desired Note that the idea is to put hisgtograms in a subfolder. More...
 
void ModifyClusters (reco::ClusterParametersList &) const override
 Scan an input collection of clusters and modify those according to the specific implementing algorithm. More...
 
float getTimeToExecute () const override
 If monitoring, recover the time to execute a particular function. More...
 
virtual void configure (const fhicl::ParameterSet &)=0
 Interface for configuring the particular algorithm tool. More...
 

Private Types

using Point = std::tuple< float, float, const reco::ClusterHit3D * >
 
using PointList = std::list< Point >
 
using MinMaxPoints = std::pair< Point, Point >
 
using MinMaxPointPair = std::pair< MinMaxPoints, MinMaxPoints >
 

Private Member Functions

reco::ClusterParametersList::iterator breakIntoTinyBits (reco::ClusterParameters &cluster, reco::ClusterParametersList::iterator positionItr, reco::ClusterParametersList &outputClusterList, int level=0) const
 Use PCA to try to find path in cluster. More...
 
float closestApproach (const Eigen::Vector3f &, const Eigen::Vector3f &, const Eigen::Vector3f &, const Eigen::Vector3f &, Eigen::Vector3f &, Eigen::Vector3f &) const
 
void buildConvexHull (reco::ClusterParameters &clusterParameters, int level=0) const
 
void buildVoronoiDiagram (reco::ClusterParameters &clusterParameters) const
 

Private Attributes

bool m_enableMonitoring
 Data members to follow. More...
 
size_t m_minTinyClusterSize
 Minimum size for a "tiny" cluster. More...
 
float m_timeToProcess
 
std::unique_ptr< lar_cluster3d::IClusterAlgm_clusterAlg
 Algorithm to do 3D space point clustering. More...
 
PrincipalComponentsAlg m_pcaAlg
 

Detailed Description

Definition at line 36 of file ClusterPathFinder_tool.cc.

Member Typedef Documentation

Definition at line 95 of file ClusterPathFinder_tool.cc.

using lar_cluster3d::ClusterPathFinder::Point = std::tuple<float, float, const reco::ClusterHit3D*>
private

Definition at line 93 of file ClusterPathFinder_tool.cc.

Definition at line 94 of file ClusterPathFinder_tool.cc.

Constructor & Destructor Documentation

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

Constructor.

Parameters
pset

Definition at line 112 of file ClusterPathFinder_tool.cc.

References configure().

113  : m_pcaAlg(pset.get<fhicl::ParameterSet>("PrincipalComponentsAlg"))
114  {
115  this->configure(pset);
116  }
void configure(fhicl::ParameterSet const &pset) override
lar_cluster3d::ClusterPathFinder::~ClusterPathFinder ( )

Destructor.

Definition at line 120 of file ClusterPathFinder_tool.cc.

120 {}

Member Function Documentation

reco::ClusterParametersList::iterator lar_cluster3d::ClusterPathFinder::breakIntoTinyBits ( reco::ClusterParameters cluster,
reco::ClusterParametersList::iterator  positionItr,
reco::ClusterParametersList outputClusterList,
int  level = 0 
) const
private

Use PCA to try to find path in cluster.

Parameters
clusterParametersThe given cluster parameters object to try to split
clusterParametersListThe list of clusters

Definition at line 222 of file ClusterPathFinder_tool.cc.

References buildConvexHull(), reco::PrincipalComponents::flipAxis(), reco::ClusterParameters::getBestEdgeList(), reco::PrincipalComponents::getEigenVectors(), reco::ClusterParameters::getFullPCA(), reco::ClusterParameters::getHitPairListPtr(), reco::ClusterParameters::getSkeletonPCA(), reco::PrincipalComponents::getSvdOK(), art::detail::indent(), art::left(), m_minTinyClusterSize, m_pcaAlg, lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_3D(), lar_cluster3d::PrincipalComponentsAlg::PCAAnalysis_calc3DDocas(), art::right(), reco::ClusterParameters::UpdateParameters(), and reco::ClusterHit2D::USED.

Referenced by getTimeToExecute(), and ModifyClusters().

227  {
228  // This needs to be a recursive routine...
229  // Idea is to take the input cluster and order 3D hits by arclength along PCA primary axis
230  // If the cluster is above the minimum number of hits then we divide into two and call ourself
231  // with the two halves. This means we form a new cluster with hits and PCA and then call ourself
232  // If the cluster is below the minimum then we can't break any more, simply add this cluster to
233  // the new output list.
234 
235  // set an indention string
236  std::string pluses(level / 2, '+');
237  std::string indent(level / 2, ' ');
238 
239  indent += pluses;
240 
241  reco::ClusterParametersList::iterator inputPositionItr = positionItr;
242 
243  // To make best use of this we'll also want the PCA for this cluster... so...
244  // Recover the prime ingredients
245  reco::PrincipalComponents& fullPCA = clusterToBreak.getFullPCA();
246 
247  std::cout << indent << ">>> breakIntoTinyBits with "
248  << clusterToBreak.getHitPairListPtr().size() << " input hits " << std::endl;
249 
250  // It turns out that computing the convex hull surrounding the points in the 2D projection onto the
251  // plane of largest spread in the PCA is a good way to break up the cluster... and we do it here since
252  // we (currently) want this to be part of the standard output
253  buildConvexHull(clusterToBreak, level + 2);
254 
255  bool storeCurrentCluster(true);
256  int minimumClusterSize(m_minTinyClusterSize);
257 
258  // Create a rough cut intended to tell us when we've reached the land of diminishing returns
259  if (clusterToBreak.getBestEdgeList().size() > 6 &&
260  clusterToBreak.getHitPairListPtr().size() > size_t(2 * minimumClusterSize)) {
261  // Recover the list of 3D hits associated to this cluster
262  reco::HitPairListPtr& clusHitPairVector = clusterToBreak.getHitPairListPtr();
263 
264  // Calculate the doca to the PCA primary axis for each 3D hit
265  // Importantly, this also gives us the arclength along that axis to the hit
266  m_pcaAlg.PCAAnalysis_calc3DDocas(clusHitPairVector, fullPCA);
267 
268  // Sort the hits along the PCA
269  clusHitPairVector.sort([](const auto& left, const auto& right) {
270  return left->getArclenToPoca() < right->getArclenToPoca();
271  });
272 
273  // Now we use the convex hull vertex points to form split points for breaking up the incoming cluster
274  reco::EdgeList& bestEdgeList = clusterToBreak.getBestEdgeList();
275  std::vector<const reco::ClusterHit3D*> vertexHitVec;
276 
277  std::cout << indent << "+> Breaking cluster, convex hull has " << bestEdgeList.size()
278  << " edges to work with" << std::endl;
279 
280  for (const auto& edge : bestEdgeList) {
281  vertexHitVec.push_back(std::get<0>(edge));
282  vertexHitVec.push_back(std::get<1>(edge));
283  }
284 
285  // Sort this vector, we aren't worried about duplicates right now...
286  std::sort(vertexHitVec.begin(), vertexHitVec.end(), [](const auto& left, const auto& right) {
287  return left->getArclenToPoca() < right->getArclenToPoca();
288  });
289 
290  // Now we create a list of pairs of iterators to the start and end of each subcluster
291  using Hit3DItrPair =
292  std::pair<reco::HitPairListPtr::iterator, reco::HitPairListPtr::iterator>;
293  using VertexPairList = std::list<Hit3DItrPair>;
294 
295  VertexPairList vertexPairList;
296  reco::HitPairListPtr::iterator firstHitItr = clusHitPairVector.begin();
297 
298  for (const auto& hit3D : vertexHitVec) {
300  std::find(firstHitItr, clusHitPairVector.end(), hit3D);
301 
302  if (vertexItr == clusHitPairVector.end()) {
303  std::cout
304  << indent
305  << ">>>>>>>>>>>>>>>>> Hit not found in input list, cannot happen? <<<<<<<<<<<<<<<<<<<"
306  << std::endl;
307  break;
308  }
309 
310  std::cout << indent << "+> -- Distance from first to current vertex point: "
311  << std::distance(firstHitItr, vertexItr) << " first: " << *firstHitItr
312  << ", vertex: " << *vertexItr;
313 
314  // Require a minimum number of points...
315  if (std::distance(firstHitItr, vertexItr) > minimumClusterSize) {
316  vertexPairList.emplace_back(Hit3DItrPair(firstHitItr, vertexItr));
317  firstHitItr = vertexItr;
318 
319  std::cout << " ++ made pair ";
320  }
321 
322  std::cout << std::endl;
323  }
324 
325  // Not done if there is distance from first to end of list
326  if (std::distance(firstHitItr, clusHitPairVector.end()) > 0) {
327  std::cout << indent << "+> loop over vertices done, remant distance: "
328  << std::distance(firstHitItr, clusHitPairVector.end()) << std::endl;
329 
330  // In the event we don't have the minimum number of hits we simply extend the last pair
331  if (!vertexPairList.empty() &&
332  std::distance(firstHitItr, clusHitPairVector.end()) < minimumClusterSize)
333  vertexPairList.back().second = clusHitPairVector.end();
334  else
335  vertexPairList.emplace_back(Hit3DItrPair(firstHitItr, clusHitPairVector.end()));
336  }
337 
338  std::cout << indent << "+> ---> breaking cluster into " << vertexPairList.size()
339  << " subclusters" << std::endl;
340 
341  if (vertexPairList.size() > 1) {
342  storeCurrentCluster = false;
343 
344  // Ok, now loop through our pairs
345  for (auto& hit3DItrPair : vertexPairList) {
346  reco::ClusterParameters clusterParams;
347  reco::HitPairListPtr& hitPairListPtr = clusterParams.getHitPairListPtr();
348 
349  std::cout << indent << "+> -- building new cluster, size: "
350  << std::distance(hit3DItrPair.first, hit3DItrPair.second) << std::endl;
351 
352  // size the container...
353  hitPairListPtr.resize(std::distance(hit3DItrPair.first, hit3DItrPair.second));
354 
355  // and copy the hits into the container
356  std::copy(hit3DItrPair.first, hit3DItrPair.second, hitPairListPtr.begin());
357 
358  // First stage of feature extraction runs here
359  m_pcaAlg.PCAAnalysis_3D(hitPairListPtr, clusterParams.getFullPCA());
360 
361  // Recover the new fullPCA
362  reco::PrincipalComponents& newFullPCA = clusterParams.getFullPCA();
363 
364  // Must have a valid pca
365  if (newFullPCA.getSvdOK()) {
366  std::cout << indent << "+> -- >> cluster has a valid Full PCA" << std::endl;
367 
368  // Need to check if the PCA direction has been reversed
369  Eigen::Vector3f fullPrimaryVec(fullPCA.getEigenVectors().row(0));
370  Eigen::Vector3f newPrimaryVec(newFullPCA.getEigenVectors().row(0));
371 
372  // If the PCA's are opposite the flip the axes
373  if (fullPrimaryVec.dot(newPrimaryVec) < 0.) {
374  for (size_t vecIdx = 0; vecIdx < 3; vecIdx++)
375  newFullPCA.flipAxis(vecIdx);
376  }
377 
378  // Set the skeleton PCA to make sure it has some value
379  clusterParams.getSkeletonPCA() = clusterParams.getFullPCA();
380 
381  positionItr =
382  breakIntoTinyBits(clusterParams, positionItr, outputClusterList, level + 4);
383  }
384  }
385  }
386  }
387 
388  // First question, are we done breaking?
389  if (storeCurrentCluster) {
390  // I think this is where we fill out the rest of the parameters?
391  // Start by adding the 2D hits...
392  // See if we can avoid duplicates by temporarily transferring to a set
393  std::set<const reco::ClusterHit2D*> hitSet;
394 
395  // Loop through 3D hits to get a set of unique 2D hits
396  for (const auto& hit3D : clusterToBreak.getHitPairListPtr()) {
397  for (const auto& hit2D : hit3D->getHits()) {
398  if (hit2D) hitSet.insert(hit2D);
399  }
400  }
401 
402  // Now add these to the new cluster
403  for (const auto& hit2D : hitSet) {
404  hit2D->setStatusBit(reco::ClusterHit2D::USED);
405  clusterToBreak.UpdateParameters(hit2D);
406  }
407 
408  std::cout << indent << "*********>>> storing new subcluster of size "
409  << clusterToBreak.getHitPairListPtr().size() << std::endl;
410 
411  positionItr = outputClusterList.insert(positionItr, clusterToBreak);
412 
413  // The above points to the element, want the next element
414  positionItr++;
415  }
416  else if (inputPositionItr == positionItr)
417  std::cout << indent << "***** DID NOT STORE A CLUSTER *****" << std::endl;
418 
419  return positionItr;
420  }
intermediate_table::iterator iterator
void PCAAnalysis_calc3DDocas(const reco::HitPairListPtr &hitPairVector, const reco::PrincipalComponents &pca) const
bool getSvdOK() const
Definition: Cluster3D.h:240
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
void flipAxis(size_t axis)
Definition: Cluster3D.cxx:232
reco::PrincipalComponents & getSkeletonPCA()
Definition: Cluster3D.h:465
reco::HitPairListPtr & getHitPairListPtr()
Definition: Cluster3D.h:463
std::list< EdgeTuple > EdgeList
Definition: Cluster3D.h:337
void buildConvexHull(reco::ClusterParameters &clusterParameters, int level=0) const
reco::PrincipalComponents & getFullPCA()
Definition: Cluster3D.h:464
std::string indent(std::size_t const i)
size_t m_minTinyClusterSize
Minimum size for a "tiny" cluster.
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
reco::ClusterParametersList::iterator breakIntoTinyBits(reco::ClusterParameters &cluster, reco::ClusterParametersList::iterator positionItr, reco::ClusterParametersList &outputClusterList, int level=0) const
Use PCA to try to find path in cluster.
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:243
void lar_cluster3d::ClusterPathFinder::buildConvexHull ( reco::ClusterParameters clusterParameters,
int  level = 0 
) const
private

Definition at line 422 of file ClusterPathFinder_tool.cc.

References util::abs(), reco::PrincipalComponents::getAvePosition(), reco::ClusterParameters::getBestEdgeList(), lar_cluster3d::ConvexHull::getConvexHull(), lar_cluster3d::ConvexHull::getConvexHullArea(), reco::PrincipalComponents::getEigenVectors(), reco::ClusterParameters::getFullPCA(), reco::ClusterParameters::getHit3DToEdgeMap(), reco::ClusterParameters::getHitPairListPtr(), reco::ClusterHit3D::getPosition(), art::detail::indent(), art::left(), and art::right().

Referenced by breakIntoTinyBits().

424  {
425  // set an indention string
426  std::string minuses(level / 2, '-');
427  std::string indent(level / 2, ' ');
428 
429  indent += minuses;
430 
431  // The plan is to build the enclosing 2D polygon around the points in the PCA plane of most spread for this cluster
432  // To do so we need to start by building a list of 2D projections onto the plane of most spread...
433  reco::PrincipalComponents& pca = clusterParameters.getFullPCA();
434 
435  // Recover the parameters from the Principal Components Analysis that we need to project and accumulate
436  const Eigen::Vector3f& pcaCenter = pca.getAvePosition();
437 
438  //dcel2d::PointList pointList;
439  using Point = std::tuple<float, float, const reco::ClusterHit3D*>;
440  using PointList = std::list<Point>;
441  PointList pointList;
442 
443  // Loop through hits and do projection to plane
444  for (const auto& hit3D : clusterParameters.getHitPairListPtr()) {
445  Eigen::Vector3f pcaToHitVec(hit3D->getPosition()[0] - pcaCenter(0),
446  hit3D->getPosition()[1] - pcaCenter(1),
447  hit3D->getPosition()[2] - pcaCenter(2));
448  Eigen::Vector3f pcaToHit = pca.getEigenVectors() * pcaToHitVec;
449 
450  pointList.emplace_back(dcel2d::Point(pcaToHit(1), pcaToHit(2), hit3D));
451  }
452 
453  // Sort the point vec by increasing x, then increase y
454  pointList.sort([](const auto& left, const auto& right) {
455  return (std::abs(std::get<0>(left) - std::get<0>(right)) >
456  std::numeric_limits<float>::epsilon()) ?
457  std::get<0>(left) < std::get<0>(right) :
458  std::get<1>(left) < std::get<1>(right);
459  });
460 
461  // containers for finding the "best" hull...
462  std::vector<ConvexHull> convexHullVec;
463  std::vector<PointList> rejectedListVec;
464  bool increaseDepth(pointList.size() > 5);
465  float lastArea(std::numeric_limits<float>::max());
466 
467  while (increaseDepth) {
468  // Get another convexHull container
469  convexHullVec.push_back(ConvexHull(pointList));
470  rejectedListVec.push_back(PointList());
471 
472  const ConvexHull& convexHull = convexHullVec.back();
473  PointList& rejectedList = rejectedListVec.back();
474  const PointList& convexHullPoints = convexHull.getConvexHull();
475 
476  increaseDepth = false;
477 
478  if (convexHull.getConvexHullArea() > 0.) {
479  std::cout << indent << "-> built convex hull, 3D hits: " << pointList.size() << " with "
480  << convexHullPoints.size() << " vertices"
481  << ", area: " << convexHull.getConvexHullArea() << std::endl;
482  std::cout << indent << "-> -Points:";
483  for (const auto& point : convexHullPoints)
484  std::cout << " (" << std::get<0>(point) << "," << std::get<1>(point) << ")";
485  std::cout << std::endl;
486 
487  if (convexHullVec.size() < 2 || convexHull.getConvexHullArea() < 0.8 * lastArea) {
488  for (auto& point : convexHullPoints) {
489  pointList.remove(point);
490  rejectedList.emplace_back(point);
491  }
492  lastArea = convexHull.getConvexHullArea();
493  // increaseDepth = true;
494  }
495  }
496  }
497 
498  // do we have a valid convexHull?
499  while (!convexHullVec.empty() && convexHullVec.back().getConvexHullArea() < 0.5) {
500  convexHullVec.pop_back();
501  rejectedListVec.pop_back();
502  }
503 
504  // If we found the convex hull then build edges around the region
505  if (!convexHullVec.empty()) {
506  size_t nRejectedTotal(0);
507  reco::HitPairListPtr hitPairListPtr = clusterParameters.getHitPairListPtr();
508 
509  for (const auto& rejectedList : rejectedListVec) {
510  nRejectedTotal += rejectedList.size();
511 
512  for (const auto& rejectedPoint : rejectedList) {
513  std::cout << indent << "-> -- Point is "
514  << convexHullVec.back().findNearestDistance(rejectedPoint)
515  << " from nearest edge" << std::endl;
516 
517  if (convexHullVec.back().findNearestDistance(rejectedPoint) > 0.5)
518  hitPairListPtr.remove(std::get<2>(rejectedPoint));
519  }
520  }
521 
522  std::cout << indent << "-> Removed " << nRejectedTotal << " leaving " << pointList.size()
523  << "/" << hitPairListPtr.size() << " points" << std::endl;
524 
525  // Now add "edges" to the cluster to describe the convex hull (for the display)
526  reco::Hit3DToEdgeMap& edgeMap = clusterParameters.getHit3DToEdgeMap();
527  reco::EdgeList& bestEdgeList = clusterParameters.getBestEdgeList();
528 
529  Point lastPoint = convexHullVec.back().getConvexHull().front();
530 
531  for (auto& curPoint : convexHullVec.back().getConvexHull()) {
532  if (curPoint == lastPoint) continue;
533 
534  const reco::ClusterHit3D* lastPoint3D = std::get<2>(lastPoint);
535  const reco::ClusterHit3D* curPoint3D = std::get<2>(curPoint);
536 
537  float distBetweenPoints = (curPoint3D->getPosition()[0] - lastPoint3D->getPosition()[0]) *
538  (curPoint3D->getPosition()[0] - lastPoint3D->getPosition()[0]) +
539  (curPoint3D->getPosition()[1] - lastPoint3D->getPosition()[1]) *
540  (curPoint3D->getPosition()[1] - lastPoint3D->getPosition()[1]) +
541  (curPoint3D->getPosition()[2] - lastPoint3D->getPosition()[2]) *
542  (curPoint3D->getPosition()[2] - lastPoint3D->getPosition()[2]);
543 
544  distBetweenPoints = std::sqrt(distBetweenPoints);
545 
546  reco::EdgeTuple edge(lastPoint3D, curPoint3D, distBetweenPoints);
547 
548  edgeMap[lastPoint3D].push_back(edge);
549  edgeMap[curPoint3D].push_back(edge);
550  bestEdgeList.emplace_back(edge);
551 
552  lastPoint = curPoint;
553  }
554  }
555 
556  return;
557  }
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
const Eigen::Vector3f getPosition() const
Definition: Cluster3D.h:155
std::tuple< float, float, const reco::ClusterHit3D * > Point
constexpr auto abs(T v)
Returns the absolute value of the argument.
reco::EdgeList & getBestEdgeList()
Definition: Cluster3D.h:468
reco::Hit3DToEdgeMap & getHit3DToEdgeMap()
Definition: Cluster3D.h:466
reco::HitPairListPtr & getHitPairListPtr()
Definition: Cluster3D.h:463
std::list< EdgeTuple > EdgeList
Definition: Cluster3D.h:337
reco::PrincipalComponents & getFullPCA()
Definition: Cluster3D.h:464
std::tuple< const reco::ClusterHit3D *, const reco::ClusterHit3D *, double > EdgeTuple
Definition: Cluster3D.h:336
std::string indent(std::size_t const i)
const Eigen::Vector3f & getAvePosition() const
Definition: Cluster3D.h:244
std::list< const reco::ClusterHit3D * > HitPairListPtr
Definition: Cluster3D.h:326
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:42
std::unordered_map< const reco::ClusterHit3D *, reco::EdgeList > Hit3DToEdgeMap
Definition: Cluster3D.h:339
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:243
void lar_cluster3d::ClusterPathFinder::buildVoronoiDiagram ( reco::ClusterParameters clusterParameters) const
private

Definition at line 559 of file ClusterPathFinder_tool.cc.

References util::abs(), voronoi2d::VoronoiDiagram::buildVoronoiDiagram(), reco::PrincipalComponents::getAvePosition(), reco::ClusterParameters::getBestEdgeList(), reco::PrincipalComponents::getEigenVectors(), reco::ClusterParameters::getFaceList(), reco::ClusterParameters::getFullPCA(), reco::ClusterParameters::getHalfEdgeList(), reco::ClusterParameters::getHit3DToEdgeMap(), reco::ClusterParameters::getHitPairListPtr(), reco::ClusterHit3D::getPosition(), reco::ClusterParameters::getVertexList(), art::left(), and art::right().

Referenced by ModifyClusters().

560  {
561  // The plan is to build the enclosing 2D polygon around the points in the PCA plane of most spread for this cluster
562  // To do so we need to start by building a list of 2D projections onto the plane of most spread...
563  reco::PrincipalComponents& pca = clusterParameters.getFullPCA();
564 
565  // Recover the parameters from the Principal Components Analysis that we need to project and accumulate
566  const Eigen::Vector3f& pcaCenter = pca.getAvePosition();
567  dcel2d::PointList pointList;
568 
569  // Loop through hits and do projection to plane
570  for (const auto& hit3D : clusterParameters.getHitPairListPtr()) {
571  Eigen::Vector3f pcaToHitVec(hit3D->getPosition()[0] - pcaCenter(0),
572  hit3D->getPosition()[1] - pcaCenter(1),
573  hit3D->getPosition()[2] - pcaCenter(2));
574  Eigen::Vector3f pcaToHit = pca.getEigenVectors() * pcaToHitVec;
575 
576  pointList.emplace_back(dcel2d::Point(pcaToHit(1), pcaToHit(2), hit3D));
577  }
578 
579  // Sort the point vec by increasing x, then increase y
580  pointList.sort([](const auto& left, const auto& right) {
581  return (std::abs(std::get<0>(left) - std::get<0>(right)) >
582  std::numeric_limits<float>::epsilon()) ?
583  std::get<0>(left) < std::get<0>(right) :
584  std::get<1>(left) < std::get<1>(right);
585  });
586 
587  // Set up the voronoi diagram builder
588  voronoi2d::VoronoiDiagram voronoiDiagram(clusterParameters.getHalfEdgeList(),
589  clusterParameters.getVertexList(),
590  clusterParameters.getFaceList());
591 
592  // And make the diagram
593  voronoiDiagram.buildVoronoiDiagram(pointList);
594 
595  // Recover the voronoi diagram vertex list and the container to store them in
596  dcel2d::VertexList& vertexList = clusterParameters.getVertexList();
597 
598  // Now get the inverse of the rotation matrix so we can get the vertex positions,
599  // which lie in the plane of the two largest PCA axes, in the standard coordinate system
600  Eigen::Matrix3f rotationMatrixInv = pca.getEigenVectors().inverse();
601 
602  // Translate and fill
603  for (auto& vertex : vertexList) {
604  Eigen::Vector3f coords = rotationMatrixInv * vertex.getCoords();
605 
606  coords += pcaCenter;
607 
608  vertex.setCoords(coords);
609  }
610 
611  // Now do the Convex Hull
612  // Now add "edges" to the cluster to describe the convex hull (for the display)
613  reco::Hit3DToEdgeMap& edgeMap = clusterParameters.getHit3DToEdgeMap();
614  reco::EdgeList& bestEdgeList = clusterParameters.getBestEdgeList();
615 
616  // const dcel2d::PointList& edgePoints = voronoiDiagram.getConvexHull();
617  // PointList localList;
618  //
619  // for(const auto& edgePoint : edgePoints) localList.emplace_back(std::get<0>(edgePoint),std::get<1>(edgePoint),std::get<2>(edgePoint));
620  //
621  // // Sort the point vec by increasing x, then increase y
622  // localList.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);});
623  //
624  // // Why do I need to do this?
625  // ConvexHull convexHull(localList);
626  //
627  // Point lastPoint = convexHull.getConvexHull().front();
628  dcel2d::Point lastPoint = voronoiDiagram.getConvexHull().front();
629 
630  // std::cout << "@@@@>> Build convex hull, voronoi handed " << edgePoints.size() << " points, convexHull cut to " << convexHull.getConvexHull().size() << std::endl;
631 
632  // for(auto& curPoint : convexHull.getConvexHull())
633  for (auto& curPoint : voronoiDiagram.getConvexHull()) {
634  if (curPoint == lastPoint) continue;
635 
636  const reco::ClusterHit3D* lastPoint3D = std::get<2>(lastPoint);
637  const reco::ClusterHit3D* curPoint3D = std::get<2>(curPoint);
638 
639  float distBetweenPoints = (curPoint3D->getPosition()[0] - lastPoint3D->getPosition()[0]) *
640  (curPoint3D->getPosition()[0] - lastPoint3D->getPosition()[0]) +
641  (curPoint3D->getPosition()[1] - lastPoint3D->getPosition()[1]) *
642  (curPoint3D->getPosition()[1] - lastPoint3D->getPosition()[1]) +
643  (curPoint3D->getPosition()[2] - lastPoint3D->getPosition()[2]) *
644  (curPoint3D->getPosition()[2] - lastPoint3D->getPosition()[2]);
645 
646  distBetweenPoints = std::sqrt(distBetweenPoints);
647 
648  reco::EdgeTuple edge(lastPoint3D, curPoint3D, distBetweenPoints);
649 
650  edgeMap[lastPoint3D].push_back(edge);
651  edgeMap[curPoint3D].push_back(edge);
652  bestEdgeList.emplace_back(edge);
653 
654  lastPoint = curPoint;
655  }
656 
657  std::cout << "****> vertexList containted " << vertexList.size() << " vertices" << std::endl;
658 
659  return;
660  }
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:102
const Eigen::Vector3f getPosition() const
Definition: Cluster3D.h:155
constexpr auto abs(T v)
Returns the absolute value of the argument.
reco::EdgeList & getBestEdgeList()
Definition: Cluster3D.h:468
reco::Hit3DToEdgeMap & getHit3DToEdgeMap()
Definition: Cluster3D.h:466
reco::HitPairListPtr & getHitPairListPtr()
Definition: Cluster3D.h:463
std::list< EdgeTuple > EdgeList
Definition: Cluster3D.h:337
VoronoiDiagram class definiton.
Definition: Voronoi.h:30
reco::PrincipalComponents & getFullPCA()
Definition: Cluster3D.h:464
std::tuple< const reco::ClusterHit3D *, const reco::ClusterHit3D *, double > EdgeTuple
Definition: Cluster3D.h:336
void buildVoronoiDiagram(const dcel2d::PointList &)
Given an input set of 2D points construct a 2D voronoi diagram.
Definition: Voronoi.cxx:135
const Eigen::Vector3f & getAvePosition() const
Definition: Cluster3D.h:244
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
Definition: AssnsIter.h:94
dcel2d::FaceList & getFaceList()
Definition: Cluster3D.h:470
std::tuple< double, double, const reco::ClusterHit3D * > Point
Definitions used by the VoronoiDiagram algorithm.
Definition: DCEL.h:42
std::list< Point > PointList
Definition: DCEL.h:43
std::unordered_map< const reco::ClusterHit3D *, reco::EdgeList > Hit3DToEdgeMap
Definition: Cluster3D.h:339
dcel2d::HalfEdgeList & getHalfEdgeList()
Definition: Cluster3D.h:472
std::list< Vertex > VertexList
Definition: DCEL.h:169
const EigenVectors & getEigenVectors() const
Definition: Cluster3D.h:243
dcel2d::VertexList & getVertexList()
Definition: Cluster3D.h:471
vertex reconstruction
float lar_cluster3d::ClusterPathFinder::closestApproach ( const Eigen::Vector3f &  P0,
const Eigen::Vector3f &  u0,
const Eigen::Vector3f &  P1,
const Eigen::Vector3f &  u1,
Eigen::Vector3f &  poca0,
Eigen::Vector3f &  poca1 
) const
private

Definition at line 662 of file ClusterPathFinder_tool.cc.

References d, DEFINE_ART_CLASS_TOOL, den, and e.

Referenced by getTimeToExecute().

668  {
669  // Technique is to compute the arclength to each point of closest approach
670  Eigen::Vector3f w0 = P0 - P1;
671  float a(1.);
672  float b(u0.dot(u1));
673  float c(1.);
674  float d(u0.dot(w0));
675  float e(u1.dot(w0));
676  float den(a * c - b * b);
677 
678  float arcLen0 = (b * e - c * d) / den;
679  float arcLen1 = (a * e - b * d) / den;
680 
681  poca0 = P0 + arcLen0 * u0;
682  poca1 = P1 + arcLen1 * u1;
683 
684  return (poca0 - poca1).norm();
685  }
Float_t den
Definition: plot.C:35
Float_t d
Definition: plot.C:235
Float_t e
Definition: plot.C:35
virtual void lar_cluster3d::IClusterModAlg::configure ( const fhicl::ParameterSet )
pure virtualinherited

Interface for configuring the particular algorithm tool.

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

Definition at line 124 of file ClusterPathFinder_tool.cc.

References fhicl::ParameterSet::get(), m_clusterAlg, m_enableMonitoring, m_minTinyClusterSize, and m_timeToProcess.

Referenced by ClusterPathFinder().

125  {
126  m_enableMonitoring = pset.get<bool>("EnableMonitoring", true);
127  m_minTinyClusterSize = pset.get<size_t>("MinTinyClusterSize", 40);
128  m_clusterAlg =
129  art::make_tool<lar_cluster3d::IClusterAlg>(pset.get<fhicl::ParameterSet>("ClusterAlg"));
130 
131  m_timeToProcess = 0.;
132 
133  return;
134  }
std::unique_ptr< lar_cluster3d::IClusterAlg > m_clusterAlg
Algorithm to do 3D space point clustering.
size_t m_minTinyClusterSize
Minimum size for a "tiny" cluster.
bool m_enableMonitoring
Data members to follow.
float lar_cluster3d::ClusterPathFinder::getTimeToExecute ( ) const
inlineoverridevirtual

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

Implements lar_cluster3d::IClusterModAlg.

Definition at line 71 of file ClusterPathFinder_tool.cc.

References breakIntoTinyBits(), closestApproach(), and m_timeToProcess.

void lar_cluster3d::ClusterPathFinder::initializeHistograms ( art::TFileDirectory &  )
overridevirtual

Interface for initializing histograms if they are desired Note that the idea is to put hisgtograms in a subfolder.

Parameters
TFileDirectory- the folder to store the hists in

Implements lar_cluster3d::IClusterModAlg.

Definition at line 136 of file ClusterPathFinder_tool.cc.

137  {
138  return;
139  }
void lar_cluster3d::ClusterPathFinder::ModifyClusters ( reco::ClusterParametersList clusterParametersList) const
overridevirtual

Scan an input collection of clusters and modify those according to the specific implementing algorithm.

Parameters
clusterParametersListA list of cluster objects (parameters from associated hits)

Top level interface for algorithm to consider pairs of clusters from the input list and determine if they are consistent with each other and, therefore, should be merged. This is done by looking at the PCA for each cluster and looking at the projection of the primary axis along the vector connecting their centers.

Implements lar_cluster3d::IClusterModAlg.

Definition at line 141 of file ClusterPathFinder_tool.cc.

References breakIntoTinyBits(), buildVoronoiDiagram(), reco::ClusterParameters::daughterList(), reco::ClusterParameters::getHitPairListPtr(), m_clusterAlg, m_enableMonitoring, m_minTinyClusterSize, and m_timeToProcess.

142  {
150  // Initial clustering is done, now trim the list and get output parameters
151  cet::cpu_timer theClockBuildClusters;
152 
153  // Start clocks if requested
154  if (m_enableMonitoring) theClockBuildClusters.start();
155 
156  int countClusters(0);
157 
158  // This is the loop over candidate 3D clusters
159  // Note that it might be that the list of candidate clusters is modified by splitting
160  // So we use the following construct to make sure we get all of them
161  reco::ClusterParametersList::iterator clusterParametersListItr = clusterParametersList.begin();
162 
163  while (clusterParametersListItr != clusterParametersList.end()) {
164  // Dereference to get the cluster paramters
165  reco::ClusterParameters& clusterParameters = *clusterParametersListItr;
166 
167  std::cout << "**> Looking at Cluster " << countClusters++ << std::endl;
168 
169  // It turns out that computing the convex hull surrounding the points in the 2D projection onto the
170  // plane of largest spread in the PCA is a good way to break up the cluster... and we do it here since
171  // we (currently) want this to be part of the standard output
172  buildVoronoiDiagram(clusterParameters);
173 
174  // Make sure our cluster has enough hits...
175  if (clusterParameters.getHitPairListPtr().size() > m_minTinyClusterSize) {
176  // Get an interim cluster list
177  reco::ClusterParametersList reclusteredParameters;
178 
179  // Call the main workhorse algorithm for building the local version of candidate 3D clusters
180  m_clusterAlg->Cluster3DHits(clusterParameters.getHitPairListPtr(), reclusteredParameters);
181 
182  std::cout << ">>>>>>>>>>> Reclustered to " << reclusteredParameters.size()
183  << " Clusters <<<<<<<<<<<<<<<" << std::endl;
184 
185  // Only process non-empty results
186  if (!reclusteredParameters.empty()) {
187  // Loop over the reclustered set
188  for (auto& cluster : reclusteredParameters) {
189  std::cout << "****> Calling breakIntoTinyBits" << std::endl;
190 
191  // Break our cluster into smaller elements...
192  breakIntoTinyBits(cluster, cluster.daughterList().end(), cluster.daughterList(), 4);
193 
194  std::cout << "****> Broke Cluster with " << cluster.getHitPairListPtr().size()
195  << " into " << cluster.daughterList().size() << " sub clusters";
196  for (auto& clus : cluster.daughterList())
197  std::cout << ", " << clus.getHitPairListPtr().size();
198  std::cout << std::endl;
199 
200  // Add the daughters to the cluster
201  clusterParameters.daughterList().insert(clusterParameters.daughterList().end(),
202  cluster);
203  }
204  }
205  }
206 
207  // Go to next cluster parameters object
208  clusterParametersListItr++;
209  }
210 
211  if (m_enableMonitoring) {
212  theClockBuildClusters.stop();
213 
214  m_timeToProcess = theClockBuildClusters.accumulated_real_time();
215  }
216 
217  mf::LogDebug("Cluster3D") << ">>>>> Cluster Path finding done" << std::endl;
218 
219  return;
220  }
intermediate_table::iterator iterator
std::unique_ptr< lar_cluster3d::IClusterAlg > m_clusterAlg
Algorithm to do 3D space point clustering.
void buildVoronoiDiagram(reco::ClusterParameters &clusterParameters) const
reco::HitPairListPtr & getHitPairListPtr()
Definition: Cluster3D.h:463
Cluster finding and building.
ClusterParametersList & daughterList()
Definition: Cluster3D.h:438
size_t m_minTinyClusterSize
Minimum size for a "tiny" cluster.
reco::ClusterParametersList::iterator breakIntoTinyBits(reco::ClusterParameters &cluster, reco::ClusterParametersList::iterator positionItr, reco::ClusterParametersList &outputClusterList, int level=0) const
Use PCA to try to find path in cluster.
bool m_enableMonitoring
Data members to follow.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
std::list< ClusterParameters > ClusterParametersList
Definition: Cluster3D.h:393

Member Data Documentation

std::unique_ptr<lar_cluster3d::IClusterAlg> lar_cluster3d::ClusterPathFinder::m_clusterAlg
private

Algorithm to do 3D space point clustering.

Definition at line 108 of file ClusterPathFinder_tool.cc.

Referenced by configure(), and ModifyClusters().

bool lar_cluster3d::ClusterPathFinder::m_enableMonitoring
private

Data members to follow.

Definition at line 103 of file ClusterPathFinder_tool.cc.

Referenced by configure(), and ModifyClusters().

size_t lar_cluster3d::ClusterPathFinder::m_minTinyClusterSize
private

Minimum size for a "tiny" cluster.

Definition at line 104 of file ClusterPathFinder_tool.cc.

Referenced by breakIntoTinyBits(), configure(), and ModifyClusters().

PrincipalComponentsAlg lar_cluster3d::ClusterPathFinder::m_pcaAlg
private

Definition at line 109 of file ClusterPathFinder_tool.cc.

Referenced by breakIntoTinyBits().

float lar_cluster3d::ClusterPathFinder::m_timeToProcess
mutableprivate

Definition at line 105 of file ClusterPathFinder_tool.cc.

Referenced by configure(), getTimeToExecute(), and ModifyClusters().


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