11 #include "art_root_io/TFileService.h" 12 #include "cetlib/cpu_timer.h" 13 #include "cetlib/search_path.h" 31 #include <Eigen/Dense> 140 using BestNodeMap = std::unordered_map<const reco::ClusterHit3D*, BestNodeTuple>;
160 using MinMaxPoints = std::pair<reco::ProjectedPoint, reco::ProjectedPoint>;
180 std::unique_ptr<lar_cluster3d::IClusterParametersBuilder>
210 fClusterBuilder = art::make_tool<lar_cluster3d::IClusterParametersBuilder>(
243 cet::cpu_timer theClockBuildClusters;
250 for (
auto& clusterParams : clusterParametersList) {
271 RunPrimsAlgorithm(clusterParams.getHitPairListPtr(), topNode, daughterParametersList);
274 cet::cpu_timer theClockBuildClusters;
282 theClockBuildClusters.stop();
288 for (
auto& daughterParams : daughterParametersList)
294 theClockBuildClusters.stop();
299 mf::LogDebug(
"MSTPathFinder") <<
">>>>> Cluster Path finding done" << std::endl;
310 if (hitPairList.empty())
return;
313 cet::cpu_timer theClockDBScan;
319 size_t clusterIdx(0);
347 curEdgeItr != curEdgeList.end();) {
349 curEdgeItr = curEdgeList.erase(curEdgeItr);
355 curClusterHitList->push_back(lastAddedHit);
359 float bestDistance(1.5);
367 for (
auto& pair : CandPairList) {
370 pair.first * lastAddedHit->
getHitChiSquare() * pair.second->getHitChiSquare();
372 curEdgeList.push_back(
reco::EdgeTuple(lastAddedHit, pair.second, edgeWeight));
377 if (curEdgeList.empty()) {
378 std::cout <<
"-----------------------------------------------------------------------------" 381 std::cout <<
"**> Cluster idx: " << clusterIdx++ <<
" has " << curClusterHitList->size()
382 <<
" hits" << std::endl;
385 freeHitItr = std::find_if(freeHitItr, hitPairList.end(), [](
const auto&
hit) {
390 if (freeHitItr == hitPairList.end())
break;
392 std::cout <<
"##################################################################>" 393 "Processing another cluster" 399 curCluster = &clusterParametersList.back();
403 lastAddedHit = *freeHitItr++;
408 curEdgeList.sort([](
const auto&
left,
const auto&
right) {
409 return std::get<2>(
left) < std::get<2>(
right);
415 (*curEdgeMap)[std::get<0>(curEdge)].push_back(curEdge);
416 (*curEdgeMap)[std::get<1>(curEdge)].push_back(
417 reco::EdgeTuple(std::get<1>(curEdge), std::get<0>(curEdge), std::get<2>(curEdge)));
420 lastAddedHit = std::get<1>(curEdge);
425 theClockDBScan.stop();
436 float bestQuality(0.);
437 float aveNumEdges(0.);
438 size_t maxNumEdges(0);
439 size_t nIsolatedHits(0);
442 cet::cpu_timer theClockPathFinding;
452 for (
const auto&
hit : hitPairList) {
459 tempList.push_front(std::get<0>(curEdgeMap[
hit].front()));
461 if (quality > bestQuality) {
462 longestCluster = tempList;
463 bestQuality = quality;
469 aveNumEdges += float(curEdgeMap[
hit].
size());
470 maxNumEdges = std::max(maxNumEdges, curEdgeMap[
hit].
size());
473 aveNumEdges /= float(hitPairList.size());
474 std::cout <<
"----> # isolated hits: " << nIsolatedHits
475 <<
", longest branch: " << longestCluster.size()
476 <<
", cluster size: " << hitPairList.size() <<
", ave # edges: " << aveNumEdges
477 <<
", max: " << maxNumEdges << std::endl;
479 if (!longestCluster.empty()) {
480 hitPairList = longestCluster;
481 for (
const auto&
hit : hitPairList) {
482 for (
const auto& edge : curEdgeMap[
hit])
483 bestEdgeList.emplace_back(edge);
486 std::cout <<
" ====> new cluster size: " << hitPairList.size() << std::endl;
490 theClockPathFinding.stop();
502 cet::cpu_timer theClockPathFinding;
524 float alpha = std::min(
float(1.), std::max(
float(0.001), pcaWidth / pcaLen));
530 for (
const auto& hit3D : curCluster) {
532 if (!curEdgeMap[hit3D].
empty() && curEdgeMap[hit3D].
size() == 1) {
533 Eigen::Vector3f pcaToHitVec(hit3D->getPosition()[0] - pcaCenter(0),
534 hit3D->getPosition()[1] - pcaCenter(1),
535 hit3D->getPosition()[2] - pcaCenter(2));
539 isolatedPointList.emplace_back(pcaToHit(2), pcaToHit(1), hit3D);
543 std::cout <<
"************* Finding best path with A* in cluster *****************" 545 std::cout <<
"**> There are " << curCluster.size() <<
" hits, " << isolatedPointList.size()
546 <<
" isolated hits, the alpha parameter is " << alpha << std::endl;
547 std::cout <<
"**> PCA len: " << pcaLen <<
", wid: " << pcaWidth <<
", height: " << pcaHeight
548 <<
", ratio: " << pcaHeight / pcaWidth << std::endl;
551 if (isolatedPointList.size() > 1) {
553 isolatedPointList.sort([](
const auto&
left,
const auto&
right) {
555 std::numeric_limits<float>::epsilon()) ?
564 std::cout <<
"**> Sorted " << isolatedPointList.size()
571 float cost(std::numeric_limits<float>::max());
573 LeastCostPath(curEdgeMap[startHit].front(), stopHit, clusterParams, cost);
578 <<
" hits, " << clusterParams.
getBestEdgeList().size() <<
" edges" << std::endl;
588 std::cout <<
"++++++>>> PCA failure! # hits: " << clusterParams.
getHitPairListPtr().size()
594 theClockPathFinding.stop();
621 bestNodeMap[startNode] =
624 while (!openList.empty()) {
629 if (openList.size() > 1)
630 currentNodeItr = std::min_element(
631 openList.begin(), openList.end(), [bestNodeMap](
const auto& next,
const auto& best) {
632 return std::get<2>(bestNodeMap.at(next)) < std::get<2>(bestNodeMap.at(best));
639 if (currentNode == goalNode) {
648 openList.erase(currentNodeItr);
652 const BestNodeTuple& currentNodeTuple = bestNodeMap.at(currentNode);
653 float currentNodeScore = std::get<1>(currentNodeTuple);
658 for (
const auto& curEdge : curEdgeList) {
663 float tentative_gScore = currentNodeScore + std::get<2>(curEdge);
668 if (candNodeItr == bestNodeMap.end()) { openList.push_back(candHit3D); }
669 else if (tentative_gScore > std::get<1>(candNodeItr->second))
675 bestNodeMap[candHit3D] =
676 BestNodeTuple(currentNode, tentative_gScore, tentative_gScore + guessToTarget);
689 while (std::get<0>(bestNodeMap.at(goalNode)) != goalNode) {
694 pathNodeList.push_front(goalNode);
695 bestEdgeList.push_front(bestEdge);
700 pathNodeList.push_front(goalNode);
708 float& showMeTheMoney)
const 715 showMeTheMoney = std::numeric_limits<float>::max();
717 if (edgeListItr != curEdgeMap.end() && !edgeListItr->second.empty()) {
721 for (
const auto& edge : edgeListItr->second) {
723 if (std::get<1>(edge) == std::get<0>(curEdge))
continue;
726 if (std::get<1>(edge) == goalNode) {
727 bestNodeList.push_back(goalNode);
728 bestEdgeList.push_back(edge);
729 showMeTheMoney = std::get<2>(edge);
734 float currentCost(0.);
738 if (currentCost < std::numeric_limits<float>::max()) {
739 showMeTheMoney = std::get<2>(edge) + currentCost;
745 if (showMeTheMoney < std::numeric_limits<float>::max()) {
756 const Eigen::Vector3f& node1Pos = node1->
getPosition();
757 const Eigen::Vector3f& node2Pos = node2->
getPosition();
758 float deltaNode[] = {
759 node1Pos[0] - node2Pos[0], node1Pos[1] - node2Pos[1], node1Pos[2] - node2Pos[2]};
762 return std::sqrt(deltaNode[0] * deltaNode[0] + deltaNode[1] * deltaNode[1] +
763 deltaNode[2] * deltaNode[2]);
786 float& bestTreeQuality)
const 789 float bestQuality(0.);
790 float curEdgeWeight = std::max(0.3, std::get<2>(curEdge));
791 float curEdgeProj(1. / curEdgeWeight);
795 if (edgeListItr != hitToEdgeMap.end()) {
797 const Eigen::Vector3f& firstHitPos = std::get<0>(curEdge)->getPosition();
798 const Eigen::Vector3f& secondHitPos = std::get<1>(curEdge)->getPosition();
799 float curEdgeVec[] = {secondHitPos[0] - firstHitPos[0],
800 secondHitPos[1] - firstHitPos[1],
801 secondHitPos[2] - firstHitPos[2]};
802 float curEdgeMag = std::sqrt(curEdgeVec[0] * curEdgeVec[0] + curEdgeVec[1] * curEdgeVec[1] +
803 curEdgeVec[2] * curEdgeVec[2]);
805 curEdgeMag = std::max(
float(0.1), curEdgeMag);
807 for (
const auto& edge : edgeListItr->second) {
809 if (std::get<1>(edge) == std::get<0>(curEdge))
continue;
815 if (quality > bestQuality) {
816 hitPairListPtr = tempList;
817 bestQuality = quality;
818 curEdgeProj = 1. / curEdgeMag;
823 hitPairListPtr.push_front(std::get<1>(curEdge));
825 bestTreeQuality += bestQuality + curEdgeProj;
827 return hitPairListPtr;
838 size_t nStartedWith(hitPairVector.size());
839 size_t nRejectedHits(0);
844 for (
const auto& hit3D : hitPairVector) {
846 size_t n2DHitsIn3DHit(0);
847 size_t nThisClusterOnly(0);
848 size_t nOtherCluster(0);
851 const std::set<const reco::ClusterHit3D*>* otherClusterHits = 0;
853 for (
const auto& hit2D : hit3D->getHits()) {
854 if (!hit2D)
continue;
858 if (hit2DToClusterMap[hit2D].
size() < 2)
859 nThisClusterOnly = hit2DToClusterMap[hit2D][&clusterParams].
size();
861 for (
const auto& clusterHitMap : hit2DToClusterMap[hit2D]) {
862 if (clusterHitMap.first == &clusterParams)
continue;
864 if (clusterHitMap.second.size() > nOtherCluster) {
865 nOtherCluster = clusterHitMap.second.size();
867 otherClusterHits = &clusterHitMap.second;
873 if (n2DHitsIn3DHit < 3 && nThisClusterOnly > 1 && nOtherCluster > 0) {
874 bool skip3DHit(
false);
876 for (
const auto& otherHit3D : *otherClusterHits) {
877 size_t nOther2DHits(0);
879 for (
const auto& otherHit2D : otherHit3D->getHits()) {
880 if (!otherHit2D)
continue;
885 if (nOther2DHits > 2) {
892 if (skip3DHit)
continue;
895 goodHits.emplace_back(hit3D);
898 std::cout <<
"###>> Input " << nStartedWith <<
" hits, rejected: " << nRejectedHits
901 hitPairVector.resize(goodHits.size());
902 std::copy(goodHits.begin(), goodHits.end(), hitPairVector.begin());
912 std::string minuses(level / 2,
'-');
913 std::string
indent(level / 2,
' ');
927 for (
const auto& hit3D : hitPairListPtr) {
928 Eigen::Vector3f pcaToHitVec(hit3D->getPosition()[0] - pcaCenter(0),
929 hit3D->getPosition()[1] - pcaCenter(1),
930 hit3D->getPosition()[2] - pcaCenter(2));
933 pointList.emplace_back(pcaToHit(1), pcaToHit(2), hit3D);
937 pointList.sort([](
const auto&
left,
const auto&
right) {
939 std::numeric_limits<float>::epsilon()) ?
945 std::vector<ConvexHull> convexHullVec;
946 std::vector<reco::ProjectedPointList> rejectedListVec;
947 bool increaseDepth(pointList.size() > 3);
948 float lastArea(std::numeric_limits<float>::max());
950 while (increaseDepth) {
955 const ConvexHull& convexHull = convexHullVec.back();
959 increaseDepth =
false;
962 if (convexHullVec.size() < 2 || convexHull.
getConvexHullArea() < 0.8 * lastArea) {
963 for (
auto& point : convexHullPoints) {
964 pointList.remove(point);
965 rejectedList.emplace_back(point);
974 while (!convexHullVec.empty() && convexHullVec.back().getConvexHullArea() < 0.5) {
975 convexHullVec.pop_back();
976 rejectedListVec.pop_back();
980 if (!convexHullVec.empty()) {
983 for (
const auto& rejectedList : rejectedListVec) {
984 for (
const auto& rejectedPoint : rejectedList) {
985 if (convexHullVec.back().findNearestDistance(rejectedPoint) > 0.5)
986 locHitPairListPtr.remove(std::get<2>(rejectedPoint));
997 for (
auto& curPoint : convexHullVec.back().getConvexHull()) {
998 if (curPoint == lastPoint)
continue;
1003 float distBetweenPoints = (curPoint3D->getPosition()[0] - lastPoint3D->
getPosition()[0]) *
1004 (curPoint3D->getPosition()[0] - lastPoint3D->
getPosition()[0]) +
1005 (curPoint3D->getPosition()[1] - lastPoint3D->
getPosition()[1]) *
1006 (curPoint3D->getPosition()[1] - lastPoint3D->
getPosition()[1]) +
1007 (curPoint3D->getPosition()[2] - lastPoint3D->
getPosition()[2]) *
1008 (curPoint3D->getPosition()[2] - lastPoint3D->
getPosition()[2]);
1010 distBetweenPoints = std::sqrt(distBetweenPoints);
1014 convexHullPointList.push_back(curPoint);
1015 edgeMap[lastPoint3D].push_back(edge);
1016 edgeMap[curPoint3D].push_back(edge);
1017 edgeList.emplace_back(edge);
1019 lastPoint = curPoint;
1026 for (
const auto& point : extremePoints)
1027 extremePointList.push_back(point);
1033 for (
const auto& kink : kinkPoints)
1034 kinkPointList.push_back(kink);
reco::HitPairListPtr & getBestHitPairListPtr()
void RunPrimsAlgorithm(const reco::HitPairListPtr &, kdTree::KdTreeNode &, reco::ClusterParametersList &) const
Driver for Prim's algorithm.
constexpr auto const & right(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
void PCAAnalysis_3D(const reco::HitPairListPtr &hitPairList, reco::PrincipalComponents &pca, bool skeletonOnly=false) const
std::unordered_map< const reco::ClusterHit3D *, BestNodeTuple > BestNodeMap
std::tuple< const reco::ClusterHit3D *, float, float > BestNodeTuple
size_t fMinTinyClusterSize
Minimum size for a "tiny" cluster.
void FindBestPathInCluster(reco::ClusterParameters &) const
Algorithm to find the best path through the given cluster.
std::tuple< float, float, const reco::ClusterHit3D * > ProjectedPoint
Projected coordinates and pointer to hit.
const Eigen::Vector3f getPosition() const
void initializeHistograms(art::TFileDirectory &) override
Interface for initializing histograms if they are desired Note that the idea is to put hisgtograms in...
std::list< ProjectedPoint > ProjectedPointList
Declaration of signal hit object.
reco::HitPairListPtr DepthFirstSearch(const reco::EdgeTuple &, const reco::Hit3DToEdgeMap &, float &) const
a depth first search to find longest branches
std::pair< MinMaxPoints, MinMaxPoints > MinMaxPointPair
~MSTPathFinder()
Destructor.
std::list< Point > PointList
The list of the projected points.
std::list< KdTreeNode > KdTreeNodeList
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
reco::EdgeList & getBestEdgeList()
reco::Hit3DToEdgeMap & getHit3DToEdgeMap()
reco::HitPairListPtr & getHitPairListPtr()
float fConvexHullMinSep
Min hit separation to conisder in convex hull.
MSTPathFinder(const fhicl::ParameterSet &)
Constructor.
unsigned int getStatusBits() const
TimeValues
enumerate the possible values for time checking if monitoring timing
IClusterModAlg interface class definiton.
void PruneAmbiguousHits(reco::ClusterParameters &, reco::Hit2DToClusterMap &) const
Prune the obvious ambiguous hits.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Implements a ConvexHull for use in clustering.
std::list< EdgeTuple > EdgeList
const EigenValues & getEigenValues() const
reco::PrincipalComponents & getFullPCA()
reco::Hit3DToEdgeMap & getConvexHullEdgeMap()
void LeastCostPath(const reco::EdgeTuple &, const reco::ClusterHit3D *, reco::ClusterParameters &, float &) const
Find the lowest cost path between two nodes using MST edges.
std::tuple< const reco::ClusterHit3D *, const reco::ClusterHit3D *, double > EdgeTuple
reco::ProjectedPointList & getConvexHullExtremePoints()
std::unordered_map< const reco::ClusterHit2D *, ClusterToHitPairSetMap > Hit2DToClusterMap
reco::ConvexHullKinkTupleList & getConvexHullKinkPoints()
Define a container for working with the convex hull.
T get(std::string const &key) const
std::unique_ptr< lar_cluster3d::IClusterParametersBuilder > fClusterBuilder
Common cluster builder tool.
std::string indent(std::size_t const i)
const Eigen::Vector3f & getAvePosition() const
Path checking algorithm has seen this hit.
void ReconstructBestPath(const reco::ClusterHit3D *, BestNodeMap &, reco::HitPairListPtr &, reco::EdgeList &) const
std::list< ConvexHullKinkTuple > ConvexHullKinkTupleList
std::list< const reco::ClusterHit3D * > HitPairListPtr
float getConvexHullArea() const
recover the area of the convex hull
The geometry of one entire detector, as served by art.
const PointList & getConvexHull() const
recover the list of convex hull vertices
PrincipalComponentsAlg fPCAAlg
float DistanceBetweenNodes(const reco::ClusterHit3D *, const reco::ClusterHit3D *) const
float fConvexHullKinkAngle
Angle to declare a kink in convex hull calc.
Detector simulation of raw signals on wires.
std::list< CandPair > CandPairList
ConvexHull class definiton.
This header file defines the interface to a principal components analysis designed to be used within ...
Encapsulate the geometry of a wire.
geo::Geometry const * fGeometry
constexpr auto const & left(const_AssnsIter< L, R, D, Dir > const &a, const_AssnsIter< L, R, D, Dir > const &b)
reco::ConvexHull & getConvexHull()
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Utility object to perform functions of association.
bool fEnableMonitoring
Data members to follow.
Encapsulate the construction of a single detector plane.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
This provides an art tool interface definition for 3D Cluster algorithms.
std::vector< float > fTimeVector
std::pair< reco::ProjectedPoint, reco::ProjectedPoint > MinMaxPoints
Add ability to build the convex hull (these needs to be split out! )
float getHitChiSquare() const
std::unordered_map< const reco::ClusterHit3D *, reco::EdgeList > Hit3DToEdgeMap
void configure(fhicl::ParameterSet const &pset) override
float getTimeToExecute() const
reco::ProjectedPointList & getConvexHullPointList()
void buildConvexHull(reco::ClusterParameters &, reco::HitPairListPtr &, int level=0) const
float getTimeToExecute() const override
If monitoring, recover the time to execute a particular function.
reco::ProjectedPointList & getProjectedPointList()
KdTreeNode & BuildKdTree(Hit3DVec::iterator, Hit3DVec::iterator, KdTreeNodeList &, int depth=0) const
Given an input set of ClusterHit3D objects, build a kd tree structure.
void ModifyClusters(reco::ClusterParametersList &) const override
Scan an input collection of clusters and modify those according to the specific implementing algorith...
reco::EdgeList & getConvexHullEdgeList()
art framework interface to geometry description
std::list< ClusterParameters > ClusterParametersList
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
const EigenVectors & getEigenVectors() const
void AStar(const reco::ClusterHit3D *, const reco::ClusterHit3D *, reco::ClusterParameters &) const
Algorithm to find shortest path between two 3D hits.
void setStatusBit(unsigned bits) const