8 #include "Pandora/AlgorithmHeaders.h" 21 HitWidthClusterMergingAlgorithm::HitWidthClusterMergingAlgorithm() :
22 m_maxConstituentHitWidth(0.5
f),
23 m_hitWidthScalingFactor(1.
f),
24 m_fittingWeight(20.
f),
25 m_minClusterWeight(0.5
f),
26 m_maxXMergeDistance(5.
f),
27 m_maxZMergeDistance(2.
f),
28 m_minMergeCosOpeningAngle(0.97
f),
29 m_minDirectionDeviationCosAngle(0.9
f),
30 m_minClusterSparseness(0.3
f)
42 for (
const Cluster *
const pCluster : *pClusterList)
48 const unsigned int numberOfProposedConstituentHits(
51 if (numberOfProposedConstituentHits == 0)
55 const float clusterSparseness(1.
f - (static_cast<float>(pCluster->GetNCaloHits()) / static_cast<float>(numberOfProposedConstituentHits)));
63 clusterVector.push_back(pCluster);
74 for (
ClusterVector::const_iterator iterCurrentCluster = clusterVector.begin(); iterCurrentCluster != clusterVector.end(); ++iterCurrentCluster)
76 const Cluster *
const pCurrentCluster = *iterCurrentCluster;
82 if (iterCurrentCluster == iterTestCluster)
85 const Cluster *
const pTestCluster = *iterTestCluster;
92 clusterAssociationMap[pCurrentCluster].m_forwardAssociations.insert(pTestCluster);
93 clusterAssociationMap[pTestCluster].m_backwardAssociations.insert(pCurrentCluster);
111 float currentMaxX(currentHigherXExtrema.GetX()), testMaxX(testHigherXExtrema.GetX());
115 if (std::fabs(testMaxX - currentMaxX) > std::numeric_limits<float>::epsilon())
116 return (testMaxX > currentMaxX);
120 if (std::fabs(testMaxX - currentMaxX) > std::numeric_limits<float>::epsilon())
121 return (testMaxX < currentMaxX);
144 CartesianVector testMergePoint(0.
f, 0.
f, 0.
f);
161 CartesianVector currentClusterDirection(0.
f, 0.
f, 0.
f), testClusterDirection(0.
f, 0.
f, 0.
f);
169 catch (
const StatusCodeException &)
180 newConstituentHitVector.insert(newConstituentHitVector.end(), testFitParameters.
GetConstituentHitVector().begin(),
183 const CartesianVector midpoint((currentFitParameters.
GetHigherXExtrema() + testMergePoint) * 0.5);
184 CartesianVector newClusterDirection(0.
f, 0.
f, 0.
f);
190 catch (
const StatusCodeException &)
209 float minDistanceSquared(std::numeric_limits<float>::max());
212 const CartesianVector &hitPosition(constituentHit.GetPositionVector());
213 const float separationDistanceSquared(hitPosition.GetDistanceSquared(position));
215 if (separationDistanceSquared < minDistanceSquared)
217 minDistanceSquared = separationDistanceSquared;
218 closestPoint = hitPosition;
226 CartesianVector &direction,
const CartesianVector &fitReferencePoint,
const float fittingWeight)
const 228 float weightSum(0.
f), weightedLSum(0.
f), weightedTSum(0.
f);
229 bool isLConstant(
true), isTConstant(
true);
236 CartesianVector axisDirection(0.
f, 0.
f, 0.
f), orthoDirection(0.
f, 0.
f, 0.
f);
237 this->
GetFittingAxes(constituentHitSubsetVector, axisDirection, orthoDirection);
240 float firstHitL(0.
f), firstHitT(0.
f);
241 this->
GetFittingCoordinates(axisDirection, constituentHitSubsetVector.begin()->GetPositionVector(), firstHitL, firstHitT);
245 const float hitWeight(constituentHit.GetHitWidth());
246 float rL(0.
f), rT(0.
f);
249 if (std::fabs(firstHitL - rL) > std::numeric_limits<float>::epsilon())
252 if (std::fabs(firstHitT - rT) > std::numeric_limits<float>::epsilon())
255 weightedLSum += rL * hitWeight;
256 weightedTSum += rT * hitWeight;
257 weightSum += hitWeight;
260 if (weightSum < std::numeric_limits<float>::epsilon())
262 std::cout <<
"HitWidthClusterMergingAlgorithm::GetWeightedGradient - hit weight in fit is negative or equivalent to zero" << std::endl;
263 throw StatusCodeException(STATUS_CODE_NOT_ALLOWED);
270 direction = (orthoDirection.GetX() < 0.f) ? orthoDirection * (-1.
f) : orthoDirection;
278 direction = (axisDirection.GetX() < 0.f) ? axisDirection * (-1.
f) : axisDirection;
282 const float weightedLMean(weightedLSum / weightSum), weightedTMean(weightedTSum / weightSum);
283 float numerator(0.
f), denominator(0.
f);
287 const float hitWeight(constituentHit.GetHitWidth());
288 float rL(0.
f), rT(0.
f);
291 numerator += hitWeight * (rL - weightedLMean) * (rT - weightedTMean);
292 denominator += hitWeight * pow(rL - weightedLMean, 2);
295 const float gradient(numerator / denominator);
301 if (direction.GetX() < 0.f)
313 std::sort(sortedConstituentHitVector.begin(), sortedConstituentHitVector.end(),
316 float weightCount(0.
f);
319 constituentHitSubsetVector.push_back(constituentHit);
321 weightCount += constituentHit.GetHitWidth();
323 if (weightCount > fittingWeight)
331 CartesianVector &axisDirection, CartesianVector &orthoDirection)
const 335 if (constituentHitSubsetPositionVector.size() < 2)
336 throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);
338 CartesianVector centroid(0.
f, 0.
f, 0.
f);
343 axisDirection = eigenVecs.at(0);
346 if (axisDirection.GetZ() < 0.f)
347 axisDirection *= -1.
f;
350 const CartesianVector yAxis(0.
f, 1.
f, 0.
f);
351 orthoDirection = yAxis.GetCrossProduct(axisDirection).GetUnitVector();
357 const CartesianVector &axisDirection,
const CartesianVector &constituentHitPosition,
float &rL,
float &rT)
const 360 const CartesianVector xAxis(1.
f, 0.
f, 0.
f);
361 const float openingAngle(axisDirection.GetOpeningAngle(xAxis)), c(std::cos(openingAngle)), s(std::sin(openingAngle));
363 rL = (c * constituentHitPosition.GetX()) + (s * constituentHitPosition.GetZ());
364 rT = (c * constituentHitPosition.GetZ()) - (s * constituentHitPosition.GetX());
371 const CartesianVector xAxis(1.
f, 0.
f, 0.
f);
372 const float openingAngle(axisDirection.GetOpeningAngle(xAxis)), c(std::cos(openingAngle)), s(std::sin(openingAngle));
373 const float deltaL(1.
f), deltaT(gradient);
375 const float x = (c * deltaL) - (s * deltaT);
376 const float z = (c * deltaT) + (s * deltaL);
378 globalDirection.SetValues(x, 0.
f, z);
379 globalDirection = globalDirection.GetUnitVector();
389 for (
const Cluster *
const pCluster : clusterVector)
393 if (primaryMapIter == clusterAssociationMap.end())
398 primaryMapIter->second.m_forwardAssociations.begin(), primaryMapIter->second.m_forwardAssociations.end());
402 for (
const Cluster *
const pConsideredCluster : primaryForwardAssociations)
404 for (
const Cluster *
const pPrimaryCluster : primaryForwardAssociations)
406 if (pConsideredCluster == pPrimaryCluster)
412 if (secondaryMapIter == clusterAssociationMap.end())
415 const ClusterSet &secondaryForwardAssociations(secondaryMapIter->second.m_forwardAssociations);
417 if (secondaryForwardAssociations.find(pConsideredCluster) != secondaryForwardAssociations.end())
419 ClusterSet &tempPrimaryForwardAssociations(tempMap.find(pCluster)->second.m_forwardAssociations);
423 if (forwardAssociationToRemove == tempPrimaryForwardAssociations.end())
426 ClusterSet &tempPrimaryBackwardAssociations(tempMap.find(pConsideredCluster)->second.m_backwardAssociations);
430 if (backwardAssociationToRemove == tempPrimaryBackwardAssociations.end())
433 tempPrimaryForwardAssociations.erase(forwardAssociationToRemove);
434 tempPrimaryBackwardAssociations.erase(backwardAssociationToRemove);
440 clusterAssociationMap = tempMap;
447 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"FittingWeight",
m_fittingWeight));
449 PANDORA_RETURN_RESULT_IF_AND_IF(
450 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinClusterWeight",
m_minClusterWeight));
452 PANDORA_RETURN_RESULT_IF_AND_IF(
453 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MaxXMergeDistance",
m_maxXMergeDistance));
455 PANDORA_RETURN_RESULT_IF_AND_IF(
456 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MaxZMergeDistance",
m_maxZMergeDistance));
458 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
461 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
464 PANDORA_RETURN_RESULT_IF_AND_IF(
465 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MaxConstituentHitWidth",
m_maxConstituentHitWidth));
467 PANDORA_RETURN_RESULT_IF_AND_IF(
468 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"HitWidthScalingFactor",
m_hitWidthScalingFactor));
470 PANDORA_RETURN_RESULT_IF_AND_IF(
471 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinClusterSparseness",
m_minClusterSparseness));
float m_minClusterWeight
The threshold hit weight of the original, unscaled cluster to be considered in the merging process...
static bool SortByNHits(const pandora::Cluster *const pLhs, const pandora::Cluster *const pRhs)
Sort clusters by number of hits, then layer span, then inner layer, then position, then pulse-height.
const pandora::CartesianVector & GetLowerXExtrema() const
Returns the lower x extremal point of the constituent hits.
static float GetOriginalTotalClusterWeight(const pandora::Cluster *const pCluster)
Sum the widths of the original, unscaled hits contained within a cluster.
static const ClusterParameters & GetClusterParameters(const pandora::Cluster *const pCluster, const ClusterToParametersMap &clusterToParametersMap)
Return the cluster parameters of a given cluster, exception thrown if not found in map [cluster -> cl...
bool IsExtremalCluster(const bool isForward, const pandora::Cluster *const pCurrentCluster, const pandora::Cluster *const pTestCluster) const
Determine which of two clusters is extremal.
pandora::CartesianVector EigenValues
float m_maxXMergeDistance
The maximum x distance between merging points of associated clusters, units cm.
Header file for the hit width cluster merging algorithm class.
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
void GetFittingAxes(const LArHitWidthHelper::ConstituentHitVector &constituentHitSubsetVector, pandora::CartesianVector &axisDirection, pandora::CartesianVector &orthoDirection) const
Obtain the axes of the fitting frame.
std::unordered_map< const pandora::Cluster *, ClusterAssociation > ClusterAssociationMap
SortByDistanceToPoint class.
static ConstituentHitVector GetConstituentHits(const pandora::Cluster *const pCluster, const float maxConstituentHitWidth, const float hitWidthScalingFactor, const bool isUniform)
Break up the cluster hits into constituent hits.
Header file for the principal curve analysis helper class.
void PopulateClusterAssociationMap(const pandora::ClusterVector &clusterVector, ClusterAssociationMap &clusterAssociationMap) const
Populate the cluster association map.
void GetGlobalDirection(const pandora::CartesianVector &axisDirection, const float gradient, pandora::CartesianVector &globalDirection) const
Translate a gradient in the fitting coordinate frame to a direction vector in the detector frame...
float m_minDirectionDeviationCosAngle
The minimum cosine opening angle of the direction of and associated cluster before and after merge...
float m_minClusterSparseness
The threshold sparseness of a cluster to be considered in the merging process.
SortByHigherExtrema class.
void FindClosestPointToPosition(const pandora::CartesianVector &position, const LArHitWidthHelper::ConstituentHitVector &constituentHitVector, pandora::CartesianVector &closestPoint) const
Determine the position of the constituent hit that lies closest to a specified position.
float m_minMergeCosOpeningAngle
The minimum cosine opening angle of the directions of associated clusters.
Header file for the geometry helper class.
Header file for the cluster helper class.
const pandora::CartesianVector & GetHigherXExtrema() const
Returns the higher x extremal point of the constituent hits.
std::vector< ConstituentHit > ConstituentHitVector
bool AreClustersAssociated(const LArHitWidthHelper::ClusterParameters ¤tClusterParameters, const LArHitWidthHelper::ClusterParameters &testClusterParameters) const
Determine whether two clusters are associated.
static unsigned int GetNProposedConstituentHits(const pandora::Cluster *const pCluster, const float maxConstituentHitWidth, const float hitWidthScalingFactor)
Return the number of constituent hits that a given cluster would be broken into.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
static pandora::CartesianPointVector GetConstituentHitPositionVector(const ConstituentHitVector &constituentHitVector)
Obtain a vector of the contituent hit central positions.
void RemoveShortcutAssociations(const pandora::ClusterVector &clusterVector, ClusterAssociationMap &clusterAssociationMap) const
Remove 'shortcut' associations from the cluster association map.
float m_maxConstituentHitWidth
The maximum hit width of a constituent hit of broken up hit, units cm.
static pandora::CartesianVector GetExtremalCoordinatesHigherX(const ConstituentHitVector &constituentHitVector)
Return the higher x extremal point of the constituent hits.
static void RunPca(const T &t, pandora::CartesianVector ¢roid, EigenValues &outputEigenValues, EigenVectors &outputEigenVectors)
Run principal component analysis using input calo hits (TPC_VIEW_U,V,W or TPC_3D; all treated as 3D p...
std::vector< pandora::CartesianVector > EigenVectors
float m_hitWidthScalingFactor
The scaling factor of the hit widths.
void GetFittingCoordinates(const pandora::CartesianVector &axisDirection, const pandora::CartesianVector &constituentHitPosition, float &rL, float &rT) const
Translate from (x, y, z) coordinates to (rL, rT) coordinates.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
void GetClusterDirection(const LArHitWidthHelper::ConstituentHitVector &constituentHitVector, pandora::CartesianVector &direction, const pandora::CartesianVector &fitReferencePoint, const float fittingWeight) const
Determine the cluster direction at a reference point by performing a weighted least squared fit to th...
float m_fittingWeight
The maximum hit weight considered in the least squared fit.
const ConstituentHitVector & GetConstituentHitVector() const
Returns the vector of constituent hits.
LArHitWidthHelper::ClusterToParametersMap m_clusterToParametersMap
The map [cluster -> cluster parameters].
float m_maxZMergeDistance
The maximum z distance between merging points of associated clusters, units cm.
void GetConstituentHitSubsetVector(const LArHitWidthHelper::ConstituentHitVector &constituentHitVector, const pandora::CartesianVector &fitReferencePoint, const float fittingWeight, LArHitWidthHelper::ConstituentHitVector &constituentHitSubsetVector) const
Obtain a vector of the minimum number of hits closest to a reference point that exceed a given weight...
void GetListOfCleanClusters(const pandora::ClusterList *const pClusterList, pandora::ClusterVector &clusterVector) const
Populate cluster vector with subset of cluster list, containing clusters judged to be clean...