12 #include "cetlib/cpu_timer.h" 32 #include <unordered_map> 110 using BestNodeMap = std::unordered_map<const reco::ClusterHit3D*, BestNodeTuple>;
148 std::unique_ptr<lar_cluster3d::IClusterParametersBuilder>
168 auto const uWireDir = uWireGeo->
Direction();
169 m_wireDir[0] = {(float)uWireDir.X(), (float)-uWireDir.Z(), (float)uWireDir.Y()};
175 auto const vWireDir = vWireGeo->
Direction();
176 m_wireDir[1] = {(float)vWireDir.X(), (float)-vWireDir.Z(), (float)vWireDir.Y()};
179 m_clusterBuilder = art::make_tool<lar_cluster3d::IClusterParametersBuilder>(
206 cet::cpu_timer theClockBuildClusters;
214 theClockBuildClusters.stop();
220 for (
auto& clusterParams : clusterParametersList)
223 mf::LogDebug(
"MinSpanTreeAlg") <<
">>>>> Cluster3DHits done, found " 224 << clusterParametersList.size() <<
" clusters" << std::endl;
233 if (hitPairList.empty())
return;
236 cet::cpu_timer theClockDBScan;
242 size_t clusterIdx(0);
270 curEdgeItr != curEdgeList.end();) {
272 curEdgeItr = curEdgeList.erase(curEdgeItr);
278 curCluster->push_back(lastAddedHit);
282 float bestDistance(1.5);
288 for (
auto& pair : CandPairList) {
290 double edgeWeight = lastAddedHit->
getHitChiSquare() * pair.second->getHitChiSquare();
292 curEdgeList.push_back(
reco::EdgeTuple(lastAddedHit, pair.second, edgeWeight));
297 if (curEdgeList.empty()) {
298 std::cout <<
"-----------------------------------------------------------------------------" 301 std::cout <<
"**> Cluster idx: " << clusterIdx++ <<
" has " << curCluster->size() <<
" hits" 305 freeHitItr = std::find_if(freeHitItr, hitPairList.end(), [](
const auto&
hit) {
310 if (freeHitItr == hitPairList.end())
break;
312 std::cout <<
"##################################################################>" 313 "Processing another cluster" 319 curClusterItr = --clusterParametersList.end();
321 curEdgeMap = &(*curClusterItr).getHit3DToEdgeMap();
322 curCluster = &(*curClusterItr).getHitPairListPtr();
323 lastAddedHit = &(*freeHitItr++);
328 curEdgeList.sort([](
const auto&
left,
const auto&
right) {
329 return std::get<2>(
left) < std::get<2>(
right);
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)));
340 lastAddedHit = std::get<1>(curEdge);
345 theClockDBScan.stop();
354 float bestQuality(0.);
355 float aveNumEdges(0.);
356 size_t maxNumEdges(0);
357 size_t nIsolatedHits(0);
360 cet::cpu_timer theClockPathFinding;
370 for (
const auto&
hit : hitPairList) {
377 tempList.push_front(std::get<0>(curEdgeMap[
hit].front()));
379 if (quality > bestQuality) {
380 longestCluster = tempList;
381 bestQuality = quality;
387 aveNumEdges += float(curEdgeMap[
hit].
size());
388 maxNumEdges = std::max(maxNumEdges, curEdgeMap[
hit].
size());
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;
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);
404 std::cout <<
" ====> new cluster size: " << hitPairList.size() << std::endl;
408 theClockPathFinding.stop();
418 cet::cpu_timer theClockPathFinding;
439 float alpha = std::min(
float(1.), std::max(
float(0.001), pcaWidth / pcaLen));
445 for (
const auto& hit3D : curCluster) {
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));
454 isolatedPointList.emplace_back(pcaToHit(2), pcaToHit(1), hit3D);
458 std::cout <<
"************* Finding best path with A* in cluster *****************" 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;
466 if (isolatedPointList.size() > 1) {
468 isolatedPointList.sort([](
const auto&
left,
const auto&
right) {
470 std::numeric_limits<float>::epsilon()) ?
479 std::cout <<
"**> Sorted " << isolatedPointList.size()
483 float cost(std::numeric_limits<float>::max());
485 LeastCostPath(curEdgeMap[startHit].front(), stopHit, clusterParams, cost);
490 <<
" hits, " << clusterParams.
getBestEdgeList().size() <<
" edges" << std::endl;
494 std::cout <<
"++++++>>> PCA failure! # hits: " << curCluster.size() << std::endl;
499 theClockPathFinding.stop();
524 bestNodeMap[startNode] =
527 while (!openList.empty()) {
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));
542 if (currentNode == goalNode) {
551 openList.erase(currentNodeItr);
555 const BestNodeTuple& currentNodeTuple = bestNodeMap.at(currentNode);
556 float currentNodeScore = std::get<1>(currentNodeTuple);
561 for (
const auto& curEdge : curEdgeList) {
566 float tentative_gScore = currentNodeScore + std::get<2>(curEdge);
571 if (candNodeItr == bestNodeMap.end()) { openList.push_back(candHit3D); }
572 else if (tentative_gScore > std::get<1>(candNodeItr->second))
578 bestNodeMap[candHit3D] =
579 BestNodeTuple(currentNode, tentative_gScore, tentative_gScore + guessToTarget);
590 while (std::get<0>(bestNodeMap.at(goalNode)) != goalNode) {
595 pathNodeList.push_front(goalNode);
596 bestEdgeList.push_front(bestEdge);
601 pathNodeList.push_front(goalNode);
607 float& showMeTheMoney)
const 614 showMeTheMoney = std::numeric_limits<float>::max();
616 if (edgeListItr != curEdgeMap.end() && !edgeListItr->second.empty()) {
620 for (
const auto& edge : edgeListItr->second) {
622 if (std::get<1>(edge) == std::get<0>(curEdge))
continue;
625 if (std::get<1>(edge) == goalNode) {
626 bestNodeList.push_back(goalNode);
627 bestEdgeList.push_back(edge);
628 showMeTheMoney = std::get<2>(edge);
633 float currentCost(0.);
637 if (currentCost < std::numeric_limits<float>::max()) {
638 showMeTheMoney = std::get<2>(edge) + currentCost;
644 if (showMeTheMoney < std::numeric_limits<float>::max()) {
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]};
659 return std::sqrt(deltaNode[0] * deltaNode[0] + deltaNode[1] * deltaNode[1] +
660 deltaNode[2] * deltaNode[2]);
665 float& bestTreeQuality)
const 668 float bestQuality(0.);
669 float curEdgeWeight = std::max(0.3, std::get<2>(curEdge));
670 float curEdgeProj(1. / curEdgeWeight);
674 if (edgeListItr != hitToEdgeMap.end()) {
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]);
684 curEdgeMag = std::max(
float(0.1), curEdgeMag);
686 for (
const auto& edge : edgeListItr->second) {
688 if (std::get<1>(edge) == std::get<0>(curEdge))
continue;
694 if (quality > bestQuality) {
695 hitPairListPtr = tempList;
696 bestQuality = quality;
697 curEdgeProj = 1. / curEdgeMag;
702 hitPairListPtr.push_front(std::get<1>(curEdge));
704 bestTreeQuality += bestQuality + curEdgeProj;
706 return hitPairListPtr;
717 size_t nStartedWith(hitPairVector.size());
718 size_t nRejectedHits(0);
723 for (
const auto& hit3D : hitPairVector) {
725 size_t n2DHitsIn3DHit(0);
726 size_t nThisClusterOnly(0);
727 size_t nOtherCluster(0);
729 const std::set<const reco::ClusterHit3D*>* otherClusterHits = 0;
731 for (
const auto& hit2D : hit3D->getHits()) {
732 if (!hit2D)
continue;
736 if (hit2DToClusterMap[hit2D].
size() < 2)
737 nThisClusterOnly = hit2DToClusterMap[hit2D][&clusterParams].
size();
739 for (
const auto& clusterHitMap : hit2DToClusterMap[hit2D]) {
740 if (clusterHitMap.first == &clusterParams)
continue;
742 if (clusterHitMap.second.size() > nOtherCluster) {
743 nOtherCluster = clusterHitMap.second.size();
744 otherClusterHits = &clusterHitMap.second;
750 if (n2DHitsIn3DHit < 3 && nThisClusterOnly > 1 && nOtherCluster > 0) {
751 bool skip3DHit(
false);
753 for (
const auto& otherHit3D : *otherClusterHits) {
754 size_t nOther2DHits(0);
756 for (
const auto& otherHit2D : otherHit3D->getHits()) {
757 if (!otherHit2D)
continue;
762 if (nOther2DHits > 2) {
769 if (skip3DHit)
continue;
772 goodHits.emplace_back(hit3D);
775 std::cout <<
"###>> Input " << nStartedWith <<
" hits, rejected: " << nRejectedHits
778 hitPairVector.resize(goodHits.size());
779 std::copy(goodHits.begin(), goodHits.end(), hitPairVector.begin());
787 return (*left).getHitPairListPtr().size() > (*right).getHitPairListPtr().size();
800 if (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();
809 if (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();
835 if (curCluster.size() > 2) {
843 std::vector<size_t> closestPlane = {0, 0, 0};
844 std::vector<float> bestAngle = {0., 0., 0.};
846 for (
size_t plane = 0; plane < 3; plane++) {
847 const std::vector<float>& wireDir =
m_wireDir[plane];
850 std::fabs(pcaAxis[0] * wireDir[0] + pcaAxis[1] * wireDir[1] + pcaAxis[2] * wireDir[2]);
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;
860 else if (dotProd > bestAngle[1]) {
861 bestAngle[2] = bestAngle[1];
862 closestPlane[2] = closestPlane[1];
863 closestPlane[1] = plane;
864 bestAngle[1] = dotProd;
867 closestPlane[2] = plane;
868 bestAngle[2] = dotProd;
879 std::cout <<
"*****************************************************************************" 882 std::cout <<
"**>>>>> longest axis: " << closestPlane[0] <<
", best angle: " << bestAngle[0]
884 std::cout <<
"**>>>>> second axis: " << closestPlane[1] <<
", best angle: " << bestAngle[1]
886 std::cout <<
" " << std::endl;
891 size_t bestPlane = closestPlane[0];
895 while (firstHitItr != localHitList.end()) {
899 while (lastHitItr != localHitList.end()) {
901 if (currentHit->
getWireIDs()[bestPlane] != (*lastHitItr)->getWireIDs()[bestPlane])
905 if (currentHit->
getHits()[bestPlane] && (*lastHitItr)->getHits()[bestPlane] &&
906 currentHit->
getHits()[bestPlane] != (*lastHitItr)->getHits()[bestPlane])
910 if ((!(currentHit->
getHits()[bestPlane]) && (*lastHitItr)->getHits()[bestPlane]) ||
911 (currentHit->
getHits()[bestPlane] && !((*lastHitItr)->getHits()[bestPlane])))
918 while (firstHitItr != lastHitItr) {
919 currentHit = *++firstHitItr;
922 firstHitItr = lastHitItr;
924 curCluster = testList;
bool operator()(const reco::ClusterHit3D *left, const reco::ClusterHit3D *right) const
reco::HitPairListPtr & getBestHitPairListPtr()
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
std::tuple< const reco::ClusterHit3D *, float, float > BestNodeTuple
std::list< reco::ClusterHit3D > HitPairList
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)
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.
bool m_enableMonitoring
Data members to follow.
const Eigen::Vector3f getPosition() const
std::list< ProjectedPoint > ProjectedPointList
Declaration of signal hit object.
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
const std::vector< size_t > & m_plane
WireGeo const * WirePtr(WireID const &wireid) const
Returns the specified wire.
reco::EdgeList & getBestEdgeList()
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
MinSpanTreeAlg(const fhicl::ParameterSet &)
Constructor.
reco::Hit3DToEdgeMap & getHit3DToEdgeMap()
reco::HitPairListPtr & getHitPairListPtr()
IClusterAlg interface class definiton.
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
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.
std::list< EdgeTuple > EdgeList
void PruneAmbiguousHits(reco::ClusterParameters &, reco::Hit2DToClusterMap &) const
Prune the obvious ambiguous hits.
const EigenValues & getEigenValues() const
std::vector< ChannelStatusVec > ChannelStatusByPlaneVec
std::unordered_map< const reco::ClusterHit3D *, BestNodeTuple > BestNodeMap
Vector_t Direction() const
std::tuple< const reco::ClusterHit3D *, const reco::ClusterHit3D *, double > EdgeTuple
Interface for a class providing readout channel mapping to geometry.
std::unordered_map< const reco::ClusterHit2D *, ClusterToHitPairSetMap > Hit2DToClusterMap
T get(std::string const &key) const
const Eigen::Vector3f & getAvePosition() const
void RunPrimsAlgorithm(reco::HitPairList &, kdTree::KdTreeNode &, reco::ClusterParametersList &) const
Driver for Prim's algorithm.
Path checking algorithm has seen this hit.
TimeValues
enumerate the possible values for time checking if monitoring timing
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
This provides an art tool interface definition for 3D Cluster algorithms.
Description of the physical geometry of one entire detector.
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.
std::vector< float > m_timeVector
Detector simulation of raw signals on wires.
std::list< CandPair > CandPairList
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)
virtual std::vector< WireID > ChannelToWire(raw::ChannelID_t channel) const =0
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
const std::vector< geo::WireID > & getWireIDs() const
geo::GeometryCore const * m_geometry
bool operator()(const reco::ClusterParametersList::iterator &left, const reco::ClusterParametersList::iterator &right)
float getHitChiSquare() const
std::unordered_map< const reco::ClusterHit3D *, reco::EdgeList > Hit3DToEdgeMap
unsigned int ChannelID_t
Type representing the ID of a readout channel.
float getTimeToExecute() const
KdTreeNode & BuildKdTree(Hit3DVec::iterator, Hit3DVec::iterator, KdTreeNodeList &, int depth=0) const
Given an input set of ClusterHit3D objects, build a kd tree structure.
geo::WireReadoutGeom const * m_wireReadoutGeom
const ClusterHit2DVec & getHits() const
art framework interface to geometry description
void ReconstructBestPath(const reco::ClusterHit3D *, BestNodeMap &, reco::HitPairListPtr &, reco::EdgeList &) const
std::list< ClusterParameters > ClusterParametersList
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
const EigenVectors & getEigenVectors() const
SetCheckHitOrder(const std::vector< size_t > &plane)
void FindBestPathInCluster(reco::ClusterParameters &) const
Algorithm to find the best path through the given cluster.
PrincipalComponentsAlg m_pcaAlg
void setStatusBit(unsigned bits) const
reco::HitPairListPtr DepthFirstSearch(const reco::EdgeTuple &, const reco::Hit3DToEdgeMap &, float &) const
a depth first search to find longest branches