11 #include "art_root_io/TFileService.h" 29 #include <Eigen/Dense> 141 using BestNodeMap = std::unordered_map<const reco::ClusterHit3D*, BestNodeTuple>;
161 using MinMaxPoints = std::pair<reco::ProjectedPoint, reco::ProjectedPoint>;
185 std::unique_ptr<lar_cluster3d::IClusterParametersBuilder>
216 fClusterBuilder = art::make_tool<lar_cluster3d::IClusterParametersBuilder>(
258 for (
auto& clusterParams : clusterParametersList) {
279 RunPrimsAlgorithm(clusterParams.getHitPairListPtr(), topNode, daughterParametersList);
290 theClockBuildClusters.stop();
296 for (
auto& daughterParams : daughterParametersList)
302 theClockBuildClusters.
stop();
319 if (hitPairList.empty())
return;
328 size_t clusterIdx(0);
356 curEdgeItr != curEdgeList.end();) {
358 curEdgeItr = curEdgeList.erase(curEdgeItr);
364 curClusterHitList->push_back(lastAddedHit);
368 float bestDistance(1.5);
376 for (
auto& pair : CandPairList) {
379 pair.first * lastAddedHit->
getHitChiSquare() * pair.second->getHitChiSquare();
381 curEdgeList.push_back(
reco::EdgeTuple(lastAddedHit, pair.second, edgeWeight));
386 if (curEdgeList.empty()) {
387 std::cout <<
"-----------------------------------------------------------------------------" 390 std::cout <<
"**> Cluster idx: " << clusterIdx++ <<
" has " << curClusterHitList->size()
394 freeHitItr = std::find_if(freeHitItr, hitPairList.end(), [](
const auto&
hit) {
399 if (freeHitItr == hitPairList.end())
break;
401 std::cout <<
"##################################################################>" 402 "Processing another cluster" 408 curCluster = &clusterParametersList.back();
412 lastAddedHit = *freeHitItr++;
417 curEdgeList.sort([](
const auto&
left,
const auto&
right) {
418 return std::get<2>(
left) < std::get<2>(
right);
424 (*curEdgeMap)[std::get<0>(curEdge)].
push_back(curEdge);
425 (*curEdgeMap)[std::get<1>(curEdge)].
push_back(
426 reco::EdgeTuple(std::get<1>(curEdge), std::get<0>(curEdge), std::get<2>(curEdge)));
429 lastAddedHit = std::get<1>(curEdge);
434 theClockDBScan.
stop();
446 float bestQuality(0.);
447 float aveNumEdges(0.);
448 size_t maxNumEdges(0);
449 size_t nIsolatedHits(0);
462 for (
const auto&
hit : hitPairList) {
469 tempList.push_front(std::get<0>(curEdgeMap[
hit].front()));
471 if (quality > bestQuality) {
472 longestCluster = tempList;
473 bestQuality = quality;
483 aveNumEdges /=
float(hitPairList.size());
484 std::cout <<
"----> # isolated hits: " << nIsolatedHits
485 <<
", longest branch: " << longestCluster.size()
486 <<
", cluster size: " << hitPairList.size() <<
", ave # edges: " << aveNumEdges
487 <<
", max: " << maxNumEdges <<
std::endl;
489 if (!longestCluster.empty()) {
490 hitPairList = longestCluster;
491 for (
const auto&
hit : hitPairList) {
492 for (
const auto& edge : curEdgeMap[
hit])
493 bestEdgeList.emplace_back(edge);
496 std::cout <<
" ====> new cluster size: " << hitPairList.size() <<
std::endl;
500 theClockPathFinding.
stop();
541 for (
const auto& hit3D : curCluster) {
543 if (!curEdgeMap[hit3D].
empty() && curEdgeMap[hit3D].
size() == 1) {
544 Eigen::Vector3f pcaToHitVec(hit3D->getPosition()[0] - pcaCenter(0),
545 hit3D->getPosition()[1] - pcaCenter(1),
546 hit3D->getPosition()[2] - pcaCenter(2));
550 isolatedPointList.emplace_back(pcaToHit(2), pcaToHit(1), hit3D);
554 std::cout <<
"************* Finding best path with A* in cluster *****************" 556 std::cout <<
"**> There are " << curCluster.size() <<
" hits, " << isolatedPointList.size()
557 <<
" isolated hits, the alpha parameter is " << alpha <<
std::endl;
558 std::cout <<
"**> PCA len: " << pcaLen <<
", wid: " << pcaWidth <<
", height: " << pcaHeight
559 <<
", ratio: " << pcaHeight / pcaWidth <<
std::endl;
562 if (isolatedPointList.size() > 1) {
564 isolatedPointList.sort([](
const auto&
left,
const auto&
right) {
566 std::numeric_limits<float>::epsilon()) ?
575 std::cout <<
"**> Sorted " << isolatedPointList.size()
584 LeastCostPath(curEdgeMap[startHit].front(), stopHit, clusterParams, cost);
599 std::cout <<
"++++++>>> PCA failure! # hits: " << clusterParams.
getHitPairListPtr().size()
605 theClockPathFinding.
stop();
635 bestNodeMap[startNode] =
640 while (!openList.empty()) {
645 if (openList.size() > 1)
646 currentNodeItr = std::min_element(
647 openList.begin(), openList.end(), [bestNodeMap](
const auto& next,
const auto& best) {
648 return std::get<2>(bestNodeMap.at(next)) < std::get<2>(bestNodeMap.at(best));
655 if (currentNode == goalNode) {
664 openList.erase(currentNodeItr);
668 const BestNodeTuple& currentNodeTuple = bestNodeMap.at(currentNode);
669 float currentNodeScore = std::get<1>(currentNodeTuple);
674 for (
const auto& curEdge : curEdgeList) {
679 float tentative_gScore = currentNodeScore + std::get<2>(curEdge);
684 if (candNodeItr == bestNodeMap.end()) { openList.push_back(candHit3D); }
685 else if (tentative_gScore > std::get<1>(candNodeItr->second))
691 bestNodeMap[candHit3D] =
692 BestNodeTuple(currentNode, tentative_gScore, tentative_gScore + guessToTarget);
706 while (std::get<0>(bestNodeMap.at(goalNode)) != goalNode) {
711 pathNodeList.push_front(goalNode);
712 bestEdgeList.push_front(bestEdge);
717 pathNodeList.push_front(goalNode);
726 float& showMeTheMoney)
const 735 if (edgeListItr != curEdgeMap.end() && !edgeListItr->second.empty()) {
739 for (
const auto& edge : edgeListItr->second) {
741 if (std::get<1>(edge) == std::get<0>(curEdge))
continue;
744 if (std::get<1>(edge) == goalNode) {
745 bestNodeList.push_back(goalNode);
746 bestEdgeList.push_back(edge);
747 showMeTheMoney = std::get<2>(edge);
752 float currentCost(0.);
757 showMeTheMoney = std::get<2>(edge) + currentCost;
775 const Eigen::Vector3f& node1Pos = node1->
getPosition();
776 const Eigen::Vector3f& node2Pos = node2->
getPosition();
777 float deltaNode[] = {
778 node1Pos[0] - node2Pos[0], node1Pos[1] - node2Pos[1], node1Pos[2] - node2Pos[2]};
781 return std::sqrt(deltaNode[0] * deltaNode[0] + deltaNode[1] * deltaNode[1] +
782 deltaNode[2] * deltaNode[2]);
806 float& bestTreeQuality)
const 809 float bestQuality(0.);
810 float curEdgeWeight =
std::max(0.3, std::get<2>(curEdge));
811 float curEdgeProj(1. / curEdgeWeight);
815 if (edgeListItr != hitToEdgeMap.end()) {
817 const Eigen::Vector3f& firstHitPos = std::get<0>(curEdge)->getPosition();
818 const Eigen::Vector3f& secondHitPos = std::get<1>(curEdge)->getPosition();
819 float curEdgeVec[] = {secondHitPos[0] - firstHitPos[0],
820 secondHitPos[1] - firstHitPos[1],
821 secondHitPos[2] - firstHitPos[2]};
822 float curEdgeMag = std::sqrt(curEdgeVec[0] * curEdgeVec[0] + curEdgeVec[1] * curEdgeVec[1] +
823 curEdgeVec[2] * curEdgeVec[2]);
825 curEdgeMag =
std::max(
float(0.1), curEdgeMag);
827 for (
const auto& edge : edgeListItr->second) {
829 if (std::get<1>(edge) == std::get<0>(curEdge))
continue;
835 if (quality > bestQuality) {
836 hitPairListPtr = tempList;
837 bestQuality = quality;
838 curEdgeProj = 1. / curEdgeMag;
843 hitPairListPtr.push_front(std::get<1>(curEdge));
845 bestTreeQuality += bestQuality + curEdgeProj;
847 return hitPairListPtr;
859 size_t nStartedWith(hitPairVector.size());
860 size_t nRejectedHits(0);
865 for (
const auto& hit3D : hitPairVector) {
867 size_t n2DHitsIn3DHit(0);
868 size_t nThisClusterOnly(0);
869 size_t nOtherCluster(0);
872 const std::set<const reco::ClusterHit3D*>* otherClusterHits = 0;
874 for (
const auto& hit2D : hit3D->getHits()) {
875 if (!hit2D)
continue;
879 if (hit2DToClusterMap[hit2D].
size() < 2)
880 nThisClusterOnly = hit2DToClusterMap[hit2D][&clusterParams].
size();
882 for (
const auto& clusterHitMap : hit2DToClusterMap[hit2D]) {
883 if (clusterHitMap.first == &clusterParams)
continue;
885 if (clusterHitMap.second.size() > nOtherCluster) {
886 nOtherCluster = clusterHitMap.second.size();
888 otherClusterHits = &clusterHitMap.second;
894 if (n2DHitsIn3DHit < 3 && nThisClusterOnly > 1 && nOtherCluster > 0) {
895 bool skip3DHit(
false);
897 for (
const auto& otherHit3D : *otherClusterHits) {
898 size_t nOther2DHits(0);
900 for (
const auto& otherHit2D : otherHit3D->getHits()) {
901 if (!otherHit2D)
continue;
906 if (nOther2DHits > 2) {
913 if (skip3DHit)
continue;
916 goodHits.emplace_back(hit3D);
919 std::cout <<
"###>> Input " << nStartedWith <<
" hits, rejected: " << nRejectedHits
922 hitPairVector.resize(goodHits.size());
923 std::copy(goodHits.begin(), goodHits.end(), hitPairVector.begin());
949 for (
const auto& hit3D : hitPairListPtr) {
950 Eigen::Vector3f pcaToHitVec(hit3D->getPosition()[0] - pcaCenter(0),
951 hit3D->getPosition()[1] - pcaCenter(1),
952 hit3D->getPosition()[2] - pcaCenter(2));
955 pointList.emplace_back(pcaToHit(1), pcaToHit(2), hit3D);
959 pointList.sort([](
const auto&
left,
const auto&
right) {
961 std::numeric_limits<float>::epsilon()) ?
967 std::vector<ConvexHull> convexHullVec;
968 std::vector<reco::ProjectedPointList> rejectedListVec;
969 bool increaseDepth(pointList.size() > 3);
972 while (increaseDepth) {
977 const ConvexHull& convexHull = convexHullVec.back();
981 increaseDepth =
false;
984 if (convexHullVec.size() < 2 || convexHull.
getConvexHullArea() < 0.8 * lastArea) {
985 for (
auto& point : convexHullPoints) {
986 pointList.remove(point);
987 rejectedList.emplace_back(point);
996 while (!convexHullVec.empty() && convexHullVec.back().getConvexHullArea() < 0.5) {
997 convexHullVec.pop_back();
998 rejectedListVec.pop_back();
1002 if (!convexHullVec.empty()) {
1003 size_t nRejectedTotal(0);
1006 for (
const auto& rejectedList : rejectedListVec) {
1007 nRejectedTotal += rejectedList.size();
1009 for (
const auto& rejectedPoint : rejectedList) {
1010 if (convexHullVec.back().findNearestDistance(rejectedPoint) > 0.5)
1011 locHitPairListPtr.remove(std::get<2>(rejectedPoint));
1022 for (
auto& curPoint : convexHullVec.back().getConvexHull()) {
1023 if (curPoint == lastPoint)
continue;
1028 float distBetweenPoints = (curPoint3D->getPosition()[0] - lastPoint3D->
getPosition()[0]) *
1029 (curPoint3D->getPosition()[0] - lastPoint3D->
getPosition()[0]) +
1030 (curPoint3D->getPosition()[1] - lastPoint3D->
getPosition()[1]) *
1031 (curPoint3D->getPosition()[1] - lastPoint3D->
getPosition()[1]) +
1032 (curPoint3D->getPosition()[2] - lastPoint3D->
getPosition()[2]) *
1033 (curPoint3D->getPosition()[2] - lastPoint3D->
getPosition()[2]);
1035 distBetweenPoints = std::sqrt(distBetweenPoints);
1039 convexHullPointList.push_back(curPoint);
1040 edgeMap[lastPoint3D].push_back(edge);
1041 edgeMap[curPoint3D].push_back(edge);
1042 edgeList.emplace_back(edge);
1044 lastPoint = curPoint;
1051 for (
const auto& point : extremePoints)
1052 extremePointList.push_back(point);
1058 for (
const auto& kink : kinkPoints)
1059 kinkPointList.push_back(kink);
1070 float largestDistance(0.);
1077 while (firstEdgeItr != convexHull.end()) {
1084 while (++nextEdgeItr != convexHull.end()) {}
1087 return largestDistance;
reco::HitPairListPtr & getBestHitPairListPtr()
void RunPrimsAlgorithm(const reco::HitPairListPtr &, kdTree::KdTreeNode &, reco::ClusterParametersList &) const
Driver for Prim's algorithm.
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.
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
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
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
art framework interface to geometry description
TimeValues
enumerate the possible values for time checking if monitoring timing
IClusterModAlg interface class definiton.
void PruneAmbiguousHits(reco::ClusterParameters &, reco::Hit2DToClusterMap &) const
Prune the obvious ambiguous hits.
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Implements a ConvexHull for use in clustering.
std::list< EdgeTuple > EdgeList
fInnerVessel push_back(Point(-578.400000, 0.000000, 0.000000))
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.
reco::ProjectedPointList & getConvexHullExtremePoints()
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.
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::tuple< const reco::ClusterHit3D *, const reco::ClusterHit3D *, double > EdgeTuple
std::list< const reco::ClusterHit3D * > HitPairListPtr
float getConvexHullArea() const
recover the area of the convex hull
static int max(int a, int b)
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.
Encapsulate the geometry of a wire.
geo::Geometry const * fGeometry
float findConvexHullEndPoints(const reco::EdgeList &, const reco::ClusterHit3D *, const reco::ClusterHit3D *) const
reco::ConvexHull & getConvexHull()
Declaration of signal hit object.
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
bool fEnableMonitoring
Data members to follow.
Encapsulate the construction of a single detector plane.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
This provides an art tool interface definition for 3D Cluster algorithms.
std::vector< float > fTimeVector
void AStar(const reco::ClusterHit3D *, const reco::ClusterHit3D *, float alpha, kdTree::KdTreeNode &, reco::ClusterParameters &) const
Algorithm to find shortest path between two 3D hits.
std::pair< reco::ProjectedPoint, reco::ProjectedPoint > MinMaxPoints
Add ability to build the convex hull (these needs to be split out! )
std::unordered_map< const reco::ClusterHit2D *, ClusterToHitPairSetMap > Hit2DToClusterMap
float getHitChiSquare() const
std::unordered_map< const reco::ClusterHit3D *, reco::EdgeList > Hit3DToEdgeMap
double accumulated_real_time() const
void configure(fhicl::ParameterSet const &pset) override
float getTimeToExecute() const
auto const & get(AssnsNode< L, R, D > const &r)
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()
std::tuple< float, float, const reco::ClusterHit3D * > ProjectedPoint
Projected coordinates and pointer to hit.
std::list< ClusterParameters > ClusterParametersList
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
const EigenVectors & getEigenVectors() const
QTextStream & endl(QTextStream &s)
void setStatusBit(unsigned bits) const