9 #include "Pandora/AlgorithmHeaders.h" 22 CosmicRaySplittingAlgorithm::CosmicRaySplittingAlgorithm() :
23 m_clusterMinLength(10.
f),
24 m_halfWindowLayers(30),
26 m_maxCosSplittingAngle(0.9925
f),
27 m_minCosMergingAngle(0.94
f),
28 m_maxTransverseDisplacement(5.
f),
29 m_maxLongitudinalDisplacement(30.
f),
30 m_maxLongitudinalDisplacementSquared(m_maxLongitudinalDisplacement * m_maxLongitudinalDisplacement)
38 const ClusterList *pClusterList = NULL;
39 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*
this, pClusterList));
50 ClusterSet splitClusters;
54 if (splitClusters.count(*bIter) > 0)
59 if (slidingFitResultMap.end() == bFitIter)
65 CartesianVector splitPosition(0.
f, 0.
f, 0.
f);
66 CartesianVector splitDirection1(0.
f, 0.
f, 0.
f);
67 CartesianVector splitDirection2(0.
f, 0.
f, 0.
f);
69 if (STATUS_CODE_SUCCESS != this->
FindBestSplitPosition(branchSlidingFitResult, splitPosition, splitDirection1, splitDirection2))
81 if (splitClusters.count(*rIter) > 0)
86 if (slidingFitResultMap.end() == rFitIter)
94 float lengthSquared1(std::numeric_limits<float>::max());
95 float lengthSquared2(std::numeric_limits<float>::max());
97 if (STATUS_CODE_SUCCESS !=
98 this->
ConfirmSplitPosition(branchSlidingFitResult, replacementSlidingFitResult, splitPosition, splitDirection1,
99 splitDirection2, lengthSquared1, lengthSquared2))
102 if (lengthSquared1 < bestLengthSquared1)
104 bestLengthSquared1 = lengthSquared1;
105 bestReplacementIter1 = rFitIter;
108 if (lengthSquared2 < bestLengthSquared2)
110 bestLengthSquared2 = lengthSquared2;
111 bestReplacementIter2 = rFitIter;
115 const Cluster *pReplacementCluster1 = NULL;
116 const Cluster *pReplacementCluster2 = NULL;
118 if (slidingFitResultMap.end() != bestReplacementIter1)
119 pReplacementCluster1 = bestReplacementIter1->first;
121 if (slidingFitResultMap.end() != bestReplacementIter2)
122 pReplacementCluster2 = bestReplacementIter2->first;
124 if (NULL == pReplacementCluster1 && NULL == pReplacementCluster2)
127 const Cluster *
const pBranchCluster = branchSlidingFitResult.
GetCluster();
130 if (pReplacementCluster1 && pReplacementCluster2)
132 if (!(this->
IdentifyCrossedTracks(pBranchCluster, pReplacementCluster1, pReplacementCluster2, splitPosition)))
135 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=,
136 this->
PerformDoubleSplit(pBranchCluster, pReplacementCluster1, pReplacementCluster2, splitPosition, splitDirection1, splitDirection2));
139 else if (pReplacementCluster1)
141 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=,
142 this->
PerformSingleSplit(pBranchCluster, pReplacementCluster1, splitPosition, splitDirection1, splitDirection2));
144 else if (pReplacementCluster2)
146 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=,
147 this->
PerformSingleSplit(pBranchCluster, pReplacementCluster2, splitPosition, splitDirection2, splitDirection1));
151 if (pReplacementCluster1)
152 splitClusters.insert(pReplacementCluster1);
154 if (pReplacementCluster2)
155 splitClusters.insert(pReplacementCluster2);
157 splitClusters.insert(pBranchCluster);
160 return STATUS_CODE_SUCCESS;
169 const Cluster *
const pCluster = *iter;
174 clusterVector.push_back(pCluster);
186 if (slidingFitResultMap.end() == slidingFitResultMap.find(*iter))
193 if (!slidingFitResultMap.insert(TwoDSlidingFitResultMap::value_type(*iter, slidingFitResult)).second)
194 throw StatusCodeException(STATUS_CODE_FAILURE);
196 catch (StatusCodeException &statusCodeException)
198 if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
199 throw statusCodeException;
208 CartesianVector &splitPosition, CartesianVector &splitDirection1, CartesianVector &splitDirection2)
const 212 bool foundSplit(
false);
218 float minL(0.
f), maxL(0.
f), minT(0.
f), maxT(0.
f);
222 const unsigned int nSamplingPoints =
static_cast<unsigned int>((maxL - minL) /
m_samplingPitch);
224 for (
unsigned int n = 0;
n < nSamplingPoints; ++
n)
226 const float alpha((0.5
f + static_cast<float>(
n)) / static_cast<float>(nSamplingPoints));
227 const float rL(minL + (maxL - minL) * alpha);
229 CartesianVector centralPosition(0.
f, 0.
f, 0.
f), forwardDirection(0.
f, 0.
f, 0.
f), backwardDirection(0.
f, 0.
f, 0.
f);
231 if ((STATUS_CODE_SUCCESS != branchSlidingFitResult.
GetGlobalFitPosition(rL, centralPosition)) ||
232 (STATUS_CODE_SUCCESS != branchSlidingFitResult.
GetGlobalFitDirection(rL + halfWindowLength, forwardDirection)) ||
233 (STATUS_CODE_SUCCESS != branchSlidingFitResult.
GetGlobalFitDirection(rL - halfWindowLength, backwardDirection)))
238 const float cosTheta(forwardDirection.GetDotProduct(backwardDirection));
240 if (cosTheta < splitCosTheta)
242 CartesianVector forwardPosition(0.
f, 0.
f, 0.
f);
243 CartesianVector backwardPosition(0.
f, 0.
f, 0.
f);
245 if ((STATUS_CODE_SUCCESS != branchSlidingFitResult.
GetGlobalFitPosition(rL + halfWindowLength, forwardPosition)) ||
246 (STATUS_CODE_SUCCESS != branchSlidingFitResult.
GetGlobalFitPosition(rL - halfWindowLength, backwardPosition)))
251 CartesianVector forwardVectorOutwards(forwardPosition - centralPosition);
252 CartesianVector backwardVectorOutwards(backwardPosition - centralPosition);
254 splitPosition = centralPosition;
255 splitDirection1 = (forwardDirection.GetDotProduct(forwardVectorOutwards) > 0.f) ? forwardDirection : forwardDirection * -1.
f;
256 splitDirection2 = (backwardDirection.GetDotProduct(backwardVectorOutwards) > 0.f) ? backwardDirection : backwardDirection * -1.
f;
257 splitCosTheta = cosTheta;
263 return STATUS_CODE_NOT_FOUND;
265 return STATUS_CODE_SUCCESS;
271 const TwoDSlidingFitResult &replacementSlidingFitResult,
const CartesianVector &splitPosition,
const CartesianVector &splitDirection1,
272 const CartesianVector &splitDirection2,
float &bestLengthSquared1,
float &bestLengthSquared2)
const 275 bestLengthSquared1 = std::numeric_limits<float>::max();
276 bestLengthSquared2 = std::numeric_limits<float>::max();
278 bool foundMatch(
false);
287 const CartesianVector branchDirection1(splitDirection1 * -1.
f);
288 const CartesianVector branchDirection2(splitDirection2 * -1.
f);
290 const float cosSplittingAngle(-splitDirection1.GetDotProduct(splitDirection2));
291 const float branchLengthSquared((branchVertex2 - branchVertex1).GetMagnitudeSquared());
292 const float replacementLengthSquared((maxPosition - minPosition).GetMagnitudeSquared());
295 for (
unsigned int iFwd = 0; iFwd < 2; ++iFwd)
297 const CartesianVector pAxis((0 == iFwd) ? (maxPosition - minPosition) : (minPosition - maxPosition));
298 const CartesianVector vtxPosition((0 == iFwd) ? minPosition : maxPosition);
299 const CartesianVector endPosition((0 == iFwd) ? maxPosition : minPosition);
300 const CartesianVector vtxDirection((0 == iFwd) ? (pAxis.GetDotProduct(minDirection) > 0 ? minDirection : minDirection * -1.f)
301 : (pAxis.GetDotProduct(maxDirection) > 0 ? maxDirection : maxDirection * -1.f));
307 const float lengthSquared((vtxPosition - splitPosition).GetMagnitudeSquared());
308 const float lengthSquared1((vtxPosition - branchVertex1).GetMagnitudeSquared());
309 const float lengthSquared2((vtxPosition - branchVertex2).GetMagnitudeSquared());
311 if (vtxDisplacement > endDisplacement)
314 if (std::min(lengthSquared, std::min(lengthSquared1, lengthSquared2)) > std::min(branchLengthSquared, replacementLengthSquared))
318 float impactL(0.
f), impactT(0.
f);
322 impactT > std::max(1.5
f, 0.577
f * (1.
f + impactL)))
326 if (vtxDirection.GetDotProduct(branchDirection1) > std::max(
m_minCosMergingAngle, cosSplittingAngle))
328 if (lengthSquared < bestLengthSquared1)
330 bestLengthSquared1 = lengthSquared;
336 if (vtxDirection.GetDotProduct(branchDirection2) > std::max(
m_minCosMergingAngle, cosSplittingAngle))
338 if (lengthSquared < bestLengthSquared2)
340 bestLengthSquared2 = lengthSquared;
347 return STATUS_CODE_NOT_FOUND;
349 return STATUS_CODE_SUCCESS;
355 const CartesianVector &splitPosition,
const CartesianVector &forwardDirection,
const CartesianVector &backwardDirection)
const 357 CaloHitList caloHitListToMove;
358 this->
GetCaloHitListToMove(pBranchCluster, pReplacementCluster, splitPosition, forwardDirection, backwardDirection, caloHitListToMove);
360 CaloHitList caloHitListToKeep;
363 if (caloHitListToKeep.empty())
364 return PandoraContentApi::MergeAndDeleteClusters(*
this, pReplacementCluster, pBranchCluster);
366 return this->
SplitCluster(pBranchCluster, pReplacementCluster, caloHitListToMove);
372 const Cluster *
const pReplacementCluster2,
const CartesianVector &splitPosition,
const CartesianVector &splitDirection1,
373 const CartesianVector &splitDirection2)
const 375 CaloHitList caloHitListToMove1, caloHitListToMove2;
376 this->
GetCaloHitListsToMove(pBranchCluster, splitPosition, splitDirection1, splitDirection2, caloHitListToMove1, caloHitListToMove2);
378 CaloHitList caloHitListToKeep1;
381 if (caloHitListToKeep1.empty())
382 return STATUS_CODE_FAILURE;
384 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->
SplitCluster(pBranchCluster, pReplacementCluster1, caloHitListToMove1));
386 CaloHitList caloHitListToKeep2;
389 if (caloHitListToKeep2.empty())
390 return PandoraContentApi::MergeAndDeleteClusters(*
this, pReplacementCluster2, pBranchCluster);
392 return this->
SplitCluster(pBranchCluster, pReplacementCluster2, caloHitListToMove2);
398 const CartesianVector &splitPosition,
const CartesianVector &forwardDirection,
const CartesianVector &backwardDirection,
399 CaloHitList &caloHitListToMove)
const 403 const CartesianVector vtxDirection((forwardPosition - vtxPosition).GetUnitVector());
405 CaloHitList caloHitListToSort, caloHitListToCheck;
406 pBranchCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitListToSort);
410 const CaloHit *
const pCaloHit = *iter;
412 if (forwardDirection.GetDotProduct(pCaloHit->GetPositionVector() - forwardPosition) > 0.
f)
414 caloHitListToMove.push_back(pCaloHit);
416 else if (forwardDirection.GetDotProduct(pCaloHit->GetPositionVector() - vtxPosition) > -1.25
f)
418 caloHitListToCheck.push_back(pCaloHit);
422 float closestLengthSquared(std::numeric_limits<float>::max());
426 const CaloHit *
const pCaloHit = *iter;
428 if (vtxDirection.GetCrossProduct(pCaloHit->GetPositionVector() - forwardPosition).GetMagnitude() >
429 backwardDirection.GetCrossProduct(pCaloHit->GetPositionVector() - forwardPosition).GetMagnitude())
432 const float lengthSquared((pCaloHit->GetPositionVector() - vtxPosition).GetMagnitudeSquared());
434 if (lengthSquared < closestLengthSquared)
435 closestLengthSquared = lengthSquared;
440 const CaloHit *
const pCaloHit = *iter;
442 if ((pCaloHit->GetPositionVector() - vtxPosition).GetMagnitudeSquared() >= closestLengthSquared)
443 caloHitListToMove.push_back(pCaloHit);
450 const CartesianVector &splitDirection1,
const CartesianVector &splitDirection2, CaloHitList &caloHitListToMove1, CaloHitList &caloHitListToMove2)
const 452 CaloHitList caloHitListToSort;
453 pBranchCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitListToSort);
457 const CaloHit *
const pCaloHit = *iter;
459 const float displacement1(splitDirection1.GetDotProduct(pCaloHit->GetPositionVector() - splitPosition));
460 const float displacement2(splitDirection2.GetDotProduct(pCaloHit->GetPositionVector() - splitPosition));
462 if (displacement1 > displacement2)
464 caloHitListToMove1.push_back(pCaloHit);
468 caloHitListToMove2.push_back(pCaloHit);
475 const Cluster *
const pReplacementCluster2,
const pandora::CartesianVector &splitPosition)
const 477 CartesianVector branchVertex1(0.
f, 0.
f, 0.
f), branchVertex2(0.
f, 0.
f, 0.
f);
483 if ((replacementVertex2 - replacementVertex1).GetMagnitudeSquared() > (branchVertex2 - branchVertex1).GetMagnitudeSquared())
492 const Cluster *
const pBranchCluster,
const CaloHitList &caloHitListToMove, CaloHitList &caloHitListToKeep)
const 494 if (caloHitListToMove.empty())
495 return STATUS_CODE_FAILURE;
497 CaloHitList caloHitList;
498 pBranchCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitList);
502 const CaloHit *
const pCaloHit = *iter;
504 if (caloHitListToMove.end() == std::find(caloHitListToMove.begin(), caloHitListToMove.end(), pCaloHit))
505 caloHitListToKeep.push_back(pCaloHit);
508 return STATUS_CODE_SUCCESS;
514 const Cluster *
const pBranchCluster,
const Cluster *
const pReplacementCluster,
const CaloHitList &caloHitListToMove)
const 516 if (caloHitListToMove.empty())
517 return STATUS_CODE_FAILURE;
521 const CaloHit *
const pCaloHit = *iter;
522 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::RemoveFromCluster(*
this, pBranchCluster, pCaloHit));
523 PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToCluster(*
this, pReplacementCluster, pCaloHit));
526 return STATUS_CODE_SUCCESS;
533 PANDORA_RETURN_RESULT_IF_AND_IF(
534 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"ClusterMinLength",
m_clusterMinLength));
536 PANDORA_RETURN_RESULT_IF_AND_IF(
537 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"SlidingFitHalfWindow",
m_halfWindowLayers));
539 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"SamplingPitch",
m_samplingPitch));
542 return STATUS_CODE_INVALID_PARAMETER;
544 PANDORA_RETURN_RESULT_IF_AND_IF(
545 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MaxCosSplittingAngle",
m_maxCosSplittingAngle));
547 PANDORA_RETURN_RESULT_IF_AND_IF(
548 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinCosMergingAngle",
m_minCosMergingAngle));
550 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
553 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
557 return STATUS_CODE_SUCCESS;
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.
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...
void GetCaloHitListsToMove(const pandora::Cluster *const pBranchCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection1, const pandora::CartesianVector &splitDirection2, pandora::CaloHitList &caloHitList1, pandora::CaloHitList &caloHitList2) const
Get lists of calo hits to move in order to split a branch cluster into two segments for case of two r...
bool IdentifyCrossedTracks(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster1, const pandora::Cluster *const pReplacementCluster2, const pandora::CartesianVector &splitPosition) const
Identify crossed tracks formed from the branch cluster and its replacement cluster.
Header file for the cosmic ray splitting algorithm class.
float m_maxLongitudinalDisplacement
maximum longitudinal displacement of associated clusters
void GetCaloHitListToMove(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &forwardDirection, const pandora::CartesianVector &backwardDirection, pandora::CaloHitList &caloHitList) const
Get list of calo hits to move in order to split a branch cluster into two segments for case of one re...
static void GetImpactParameters(const LArPointingCluster::Vertex &pointingVertex, const LArPointingCluster::Vertex &targetVertex, float &longitudinal, float &transverse)
Calculate impact parameters between a pair of pointing vertices.
float GetLayerFitHalfWindowLength() const
Get the layer fit half window length.
static pandora::HitType GetClusterHitType(const pandora::Cluster *const pCluster)
Get the hit type associated with a two dimensional cluster.
pandora::StatusCode Run()
pandora::CartesianVector GetGlobalMinLayerDirection() const
Get global direction corresponding to the fit result in minimum fit layer.
float m_clusterMinLength
minimum length of clusters for this algorithm
pandora::StatusCode SplitCluster(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster, const pandora::CaloHitList &caloHitListToMove) const
Split the branch cluster and add hits to the replacement cluster.
Header file for the geometry helper class.
pandora::StatusCode PerformSingleSplit(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &forwardDirection, const pandora::CartesianVector &backwardDirection) const
Split a branch cluster for case of one replacement cluster.
float m_minCosMergingAngle
largest relative angle allowed when merging clusters
Header file for the cluster helper class.
pandora::StatusCode GetCaloHitListToKeep(const pandora::Cluster *const pBranchCluster, const pandora::CaloHitList &caloHitListToMove, pandora::CaloHitList &caloHitListToKeep) const
Split the branch cluster and add hits to the replacement cluster.
pandora::StatusCode GetGlobalFitPosition(const float rL, pandora::CartesianVector &position) const
Get global fit position for a given longitudinal coordinate.
float m_maxCosSplittingAngle
smallest scatter angle allowed when splitting cluster
float m_maxLongitudinalDisplacementSquared
longitudinal displacement squared
pandora::StatusCode PerformDoubleSplit(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster1, const pandora::Cluster *const pReplacementCluster2, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection1, const pandora::CartesianVector &splitDirection2) const
Split a branch cluster for case of two replacement clusters.
pandora::CartesianVector GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.
static float GetWirePitch(const pandora::Pandora &pandora, const pandora::HitType view, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
std::unordered_map< const pandora::Cluster *, TwoDSlidingFitResult > TwoDSlidingFitResultMap
static void GetExtremalCoordinates(const pandora::ClusterList &clusterList, pandora::CartesianVector &innerCoordinate, pandora::CartesianVector &outerCoordinate)
Get positions of the two most distant calo hits in a list of cluster (ordered by Z) ...
void BuildSlidingFitResultMap(const pandora::ClusterVector &clusterVector, TwoDSlidingFitResultMap &slidingFitResultMap) const
Build the map of sliding fit results.
float m_maxTransverseDisplacement
maximum transverse displacement of associated clusters
const pandora::Cluster * GetCluster() const
Get the address of the cluster, if originally provided.
pandora::StatusCode GetGlobalFitDirection(const float rL, pandora::CartesianVector &direction) const
Get global fit direction for a given longitudinal coordinate.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
pandora::StatusCode ConfirmSplitPosition(const TwoDSlidingFitResult &branchSlidingFitResult, const TwoDSlidingFitResult &replacementSlidingFitResult, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection1, const pandora::CartesianVector &splitDirection2, float &lengthSquared1, float &lengthSquared2) const
Find a second replacement cluster that aligns with the scatter of the first branch cluster...
pandora::CartesianVector GetGlobalMaxLayerDirection() const
Get global direction corresponding to the fit result in maximum fit layer.
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
pandora::StatusCode FindBestSplitPosition(const TwoDSlidingFitResult &slidingFitResult, pandora::CartesianVector &splitPosition, pandora::CartesianVector &splitDirection1, pandora::CartesianVector &splitDirection2) const
Find the position of greatest scatter along a sliding linear fit.
unsigned int m_halfWindowLayers
number of layers to use for half-window of sliding fit
static pandora::CartesianVector GetClosestPosition(const pandora::CartesianVector &position, const pandora::ClusterList &clusterList)
Get closest position in a list of clusters to a specified input position vector.
void GetLocalPosition(const pandora::CartesianVector &position, float &rL, float &rT) const
Get local sliding fit coordinates for a given global position.
pandora::CartesianVector GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.
TwoDSlidingFitResult class.
float m_samplingPitch
sampling pitch for walking along sliding linear fit
static float GetClosestDistance(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2)
Get closest distance between clusters in a pair of cluster lists.