10 #include "Pandora/AlgorithmHeaders.h"
22 TwoViewCosmicRayRemovalTool::TwoViewCosmicRayRemovalTool() :
24 m_slidingFitWindow(10000),
25 m_distanceToLine(0.5f),
26 m_minContaminationLength(3.f),
27 m_maxDistanceToHit(1.f),
28 m_minRemnantClusterSize(3),
29 m_maxDistanceToTrack(2.f)
39 if (PandoraContentApi::GetSettings(*m_pParentAlgorithm)->ShouldDisplayAlgorithmInfo())
40 std::cout <<
"----> Running Algorithm Tool: " << this->GetInstanceName() <<
", " << this->GetType() << std::endl;
42 bool changesMade(
false);
47 ClusterSet usedKeyClusters;
48 for (
const Cluster *
const pKeyCluster : sortedKeyClusters)
50 if (usedKeyClusters.count(pKeyCluster))
56 for (
const MatrixType::Element &
element : elementList)
57 usedKeyClusters.insert(
element.GetCluster1());
61 changesMade = (changesMade ? changesMade : changesMadeInIteration);
71 ClusterSet modifiedClusters, checkedClusters;
73 for (
const MatrixType::Element &
element : elementList)
79 if (checkedClusters.count(pDeltaRayCluster))
83 if ((modifiedClusters.count(
element.GetCluster1())) || (modifiedClusters.count(
element.GetCluster2())))
95 checkedClusters.insert(pDeltaRayCluster);
98 CaloHitList deltaRayHits;
101 if (deltaRayHits.empty())
105 CaloHitList deltaRayRemnantHits;
106 if (this->
GrowSeed(
element, hitType, deltaRayHits, deltaRayRemnantHits) != STATUS_CODE_SUCCESS)
109 if (deltaRayHits.size() == pDeltaRayCluster->GetNCaloHits())
112 modifiedClusters.insert(pDeltaRayCluster);
118 return !modifiedClusters.empty();
145 if (otherHitType == hitType)
153 float xMinDR(+std::numeric_limits<float>::max()), xMaxDR(-std::numeric_limits<float>::max());
154 pDeltaRayCluster->GetClusterSpanX(xMinDR, xMaxDR);
156 float xMinCR(-std::numeric_limits<float>::max()), xMaxCR(+std::numeric_limits<float>::max());
157 pMuonCluster->GetClusterSpanX(xMinCR, xMaxCR);
159 if ((xMinDR < xMinCR) || (xMaxDR > xMaxCR))
162 float zMinDR(+std::numeric_limits<float>::max()), zMaxDR(-std::numeric_limits<float>::max());
163 pDeltaRayCluster->GetClusterSpanZ(xMinDR, xMaxDR, zMinDR, zMaxDR);
165 float zMinCR(-std::numeric_limits<float>::max()), zMaxCR(+std::numeric_limits<float>::max());
166 pMuonCluster->GetClusterSpanZ(xMinCR, xMaxCR, zMinCR, zMaxCR);
168 if ((zMinDR < zMinCR) || (zMaxDR > zMaxCR))
187 CartesianVector muonDirection(0.f, 0.f, 0.f);
190 CartesianVector deltaRayVertex(0.f, 0.f, 0.f), muonVertex(0.f, 0.f, 0.f);
194 float furthestSeparation(0.f);
195 CartesianVector extendedPoint(0.f, 0.f, 0.f);
197 CaloHitList deltaRayHitList;
198 pDeltaRayCluster->GetOrderedCaloHitList().FillCaloHitList(deltaRayHitList);
200 for (
const CaloHit *
const pCaloHit : deltaRayHitList)
202 const CartesianVector &position(pCaloHit->GetPositionVector());
203 const float separation((position - muonVertex).GetMagnitude());
205 if (separation > furthestSeparation)
209 furthestSeparation = separation;
210 extendedPoint = position;
220 CaloHitList muonHitList;
221 pMuonCluster->GetOrderedCaloHitList().FillCaloHitList(muonHitList);
223 for (
const CaloHit *
const pCaloHit : muonHitList)
225 if (this->
IsInLineSegment(deltaRayVertex, extendedPoint, pCaloHit->GetPositionVector()))
235 const CartesianVector &hitPosition,
const CartesianVector &lineStart,
const CartesianVector &lineEnd,
const float distanceToLine)
const
237 CartesianVector lineDirection(lineStart - lineEnd);
238 lineDirection = lineDirection.GetUnitVector();
240 const float transverseDistanceFromLine(lineDirection.GetCrossProduct(hitPosition - lineStart).GetMagnitude());
242 if (transverseDistanceFromLine > distanceToLine)
251 const CartesianVector &lowerBoundary,
const CartesianVector &upperBoundary,
const CartesianVector &point)
const
253 const float segmentBoundaryGradient = (-1.f) * (upperBoundary.GetX() - lowerBoundary.GetX()) / (upperBoundary.GetZ() - lowerBoundary.GetZ());
254 const float xPointOnUpperLine((point.GetZ() - upperBoundary.GetZ()) / segmentBoundaryGradient + upperBoundary.GetX());
255 const float xPointOnLowerLine((point.GetZ() - lowerBoundary.GetZ()) / segmentBoundaryGradient + lowerBoundary.GetX());
257 if (std::fabs(xPointOnUpperLine - point.GetX()) < std::numeric_limits<float>::epsilon())
260 if (std::fabs(xPointOnLowerLine - point.GetX()) < std::numeric_limits<float>::epsilon())
263 if ((point.GetX() > xPointOnUpperLine) && (point.GetX() > xPointOnLowerLine))
266 if ((point.GetX() < xPointOnUpperLine) && (point.GetX() < xPointOnLowerLine))
276 const float chiSquared(element.GetOverlapResult().GetReducedChiSquared());
277 const unsigned int hitSum(element.GetCluster1()->GetNCaloHits() + element.GetCluster2()->GetNCaloHits());
279 for (
const MatrixType::Element &testElement : elementList)
284 if ((testElement.GetCluster1() == element.GetCluster1()) && (testElement.GetCluster2() == element.GetCluster2()))
287 const float testChiSquared(testElement.GetOverlapResult().GetReducedChiSquared());
288 const unsigned int testHitSum(testElement.GetCluster1()->GetNCaloHits() + testElement.GetCluster2()->GetNCaloHits());
290 if ((testHitSum < hitSum) || ((testHitSum == hitSum) && (testChiSquared > chiSquared)))
305 CartesianVector positionOnMuon(0.f, 0.f, 0.f), muonDirection(0.f, 0.f, 0.f);
310 CartesianPointVector deltaRayProjectedPositions;
314 CaloHitList deltaRayHitList;
317 CaloHitVector deltaRayHitVector(deltaRayHitList.begin(), deltaRayHitList.end());
320 for (
const CaloHit *
const pCaloHit : deltaRayHitVector)
322 const CartesianVector &position(pCaloHit->GetPositionVector());
324 for (
const CartesianVector &projectedPosition : deltaRayProjectedPositions)
326 const float distanceToProjectionSquared((position - projectedPosition).GetMagnitudeSquared());
334 const float distanceToMuonHitsSquared(muonDirection.GetCrossProduct(position - positionOnMuon).GetMagnitudeSquared());
336 if (distanceToMuonHitsSquared < m_maxDistanceToHit * m_maxDistanceToHit)
339 collectedHits.push_back(pCaloHit);
350 const MatrixType::Element &
element,
const HitType hitType, CaloHitList &deltaRayHits, CaloHitList &remnantHits)
const
355 return STATUS_CODE_NOT_FOUND;
358 CartesianVector positionOnMuon(0.f, 0.f, 0.f), muonDirection(0.f, 0.f, 0.f);
360 muonDirection) != STATUS_CODE_SUCCESS)
361 return STATUS_CODE_NOT_FOUND;
369 return STATUS_CODE_SUCCESS;
375 const Cluster *
const pDeltaRayCluster,
const float &minDistanceFromMuon,
const bool demandCloserToCollected,
376 const CaloHitList &protectedHits, CaloHitList &collectedHits)
const
378 CaloHitList deltaRayHitList;
379 pDeltaRayCluster->GetOrderedCaloHitList().FillCaloHitList(deltaRayHitList);
381 bool hitsAdded(
false);
387 for (
const CaloHit *
const pCaloHit : deltaRayHitList)
389 if (std::find(protectedHits.begin(), protectedHits.end(), pCaloHit) != protectedHits.end())
392 if (std::find(collectedHits.begin(), collectedHits.end(), pCaloHit) != collectedHits.end())
395 const float distanceToCollectedHits(collectedHits.empty()
396 ? std::numeric_limits<float>::max()
398 const float distanceToMuonHits(muonDirection.GetCrossProduct(pCaloHit->GetPositionVector() - positionOnMuon).GetMagnitude());
400 if ((distanceToMuonHits < minDistanceFromMuon) || (demandCloserToCollected && (distanceToCollectedHits > distanceToMuonHits)))
403 collectedHits.push_back(pCaloHit);
412 const MatrixType::Element &
element,
const HitType hitType, CaloHitList &collectedHits, CaloHitList &deltaRayRemnantHits)
const
422 std::string originalListName, fragmentListName;
423 ClusterList originalClusterList(1, pDeltaRayCluster);
425 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
427 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
428 PandoraContentApi::InitializeFragmentation(*
m_pParentAlgorithm, originalClusterList, originalListName, fragmentListName));
430 CaloHitList deltaRayHitList;
431 pDeltaRayCluster->GetOrderedCaloHitList().FillCaloHitList(deltaRayHitList);
433 const Cluster *pDeltaRay(
nullptr), *pDeltaRayRemnant(
nullptr);
435 for (
const CaloHit *
const pCaloHit : deltaRayHitList)
437 const bool isDeltaRay(std::find(collectedHits.begin(), collectedHits.end(), pCaloHit) != collectedHits.end());
438 const bool isDeltaRayRemnant(
439 !isDeltaRay && (std::find(deltaRayRemnantHits.begin(), deltaRayRemnantHits.end(), pCaloHit) != deltaRayRemnantHits.end()));
440 const Cluster *&pCluster(isDeltaRay ? pDeltaRay : isDeltaRayRemnant ? pDeltaRayRemnant : pMuonCluster);
445 parameters.m_caloHitList.push_back(pCaloHit);
446 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Cluster::Create(*
m_pParentAlgorithm, parameters, pCluster));
450 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToCluster(*
m_pParentAlgorithm, pCluster, pCaloHit));
454 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::EndFragmentation(*
m_pParentAlgorithm, fragmentListName, originalListName));
458 if (pDeltaRayRemnant)
459 this->
ReclusterRemnant(hitType, pMuonCluster, pDeltaRayRemnant, clusterVector, pfoVector);
461 clusterVector.push_back(pMuonCluster);
462 pfoVector.push_back(element.GetOverlapResult().GetCommonMuonPfoList().front());
463 clusterVector.push_back(pDeltaRay);
464 pfoVector.push_back(
nullptr);
472 const Cluster *
const pDeltaRayRemnant,
ClusterVector &clusterVector, PfoVector &pfoVector)
const
474 std::string caloHitListName(hitType == TPC_VIEW_U ?
"CaloHitListU" : hitType == TPC_VIEW_V ?
"CaloHitListV" :
"CaloHitListW");
476 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ReplaceCurrentList<CaloHit>(*
m_pParentAlgorithm, caloHitListName));
477 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
479 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Delete<Cluster>(*
m_pParentAlgorithm, pDeltaRayRemnant));
481 const ClusterList *pClusterList(
nullptr);
482 std::string newClusterListName;
483 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
486 const ClusterList remnantClusters(pClusterList->begin(), pClusterList->end());
488 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
490 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
493 for (
const Cluster *
const pRemnant : remnantClusters)
499 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::MergeAndDeleteClusters(*
m_pParentAlgorithm, pMuonCluster, pRemnant));
504 clusterVector.push_back(pRemnant);
505 pfoVector.push_back(
nullptr);
512 const MatrixType::Element &
element,
const HitType hitType, CartesianPointVector &projectedPositions)
const
514 HitTypeVector hitTypes({TPC_VIEW_U, TPC_VIEW_V, TPC_VIEW_W});
516 hitTypes.erase(std::find(hitTypes.begin(), hitTypes.end(), hitType));
521 if (!pCluster1 || !pCluster2)
522 return STATUS_CODE_NOT_FOUND;
531 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinSeparation",
m_minSeparation));
533 PANDORA_RETURN_RESULT_IF_AND_IF(
534 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"SlidingFitWindow",
m_slidingFitWindow));
536 PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"DistanceToLine",
m_distanceToLine));
538 PANDORA_RETURN_RESULT_IF_AND_IF(
539 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinContaminationLength",
m_minContaminationLength));
541 PANDORA_RETURN_RESULT_IF_AND_IF(
542 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MaxDistanceToHit",
m_maxDistanceToHit));
544 PANDORA_RETURN_RESULT_IF_AND_IF(
545 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinRemnantClusterSize",
m_minRemnantClusterSize));
547 PANDORA_RETURN_RESULT_IF_AND_IF(
548 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MaxDistanceToTrack",
m_maxDistanceToTrack));
550 return STATUS_CODE_SUCCESS;
void SplitDeltaRayCluster(const MatrixType::Element &element, const pandora::HitType hitType, pandora::CaloHitList &collectedHits, pandora::CaloHitList &deltaRayRemnantHits) const
Fragment the delta ray cluster, refining the matched delta ray cluster perhaps creating significant r...
void CollectHitsFromDeltaRay(const pandora::CartesianVector &positionOnMuon, const pandora::CartesianVector &muonDirection, const pandora::Cluster *const pDeltaRayCluster, const float &minDistanceFromMuon, const bool demandCloserToCollected, const pandora::CaloHitList &protectedHits, pandora::CaloHitList &collectedHits) const
Collect hits from the delta ray cluster to either keep (demandCloserToCollected == true) or separate ...
unsigned int m_slidingFitWindow
The sliding fit window used in cosmic ray parameterisations.
Header file for the pfo helper class.
pandora::StatusCode GetProjectedPositions(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2, pandora::CartesianPointVector &projectedPositions) const
Use two clusters from different views to calculate projected positions in the remaining third view...
bool RemoveCosmicRayHits(const MatrixType::ElementList &elementList) const
Remove hits from delta ray clusters that belong to the parent muon.
void CreateSeed(const MatrixType::Element &element, const pandora::HitType hitType, pandora::CaloHitList &collectedHits) const
Collect hits in the delta ray cluster that lie close to calculated projected positions forming a seed...
TwoViewDeltaRayMatchingAlgorithm class.
void UpdateForNewClusters(const pandora::ClusterVector &newClusterVector, const pandora::PfoVector &pfoVector)
Add a new cluster to algorithm ownership maps and, if it a delta ray cluster, to the underlying match...
const pandora::Cluster * GetCluster(const MatrixType::Element &element, const pandora::HitType hitType)
Get the address of the given hit type cluster.
std::vector< Element > ElementList
static void GetClosestPositions(const pandora::Cluster *const pCluster1, const pandora::Cluster *const pCluster2, pandora::CartesianVector &position1, pandora::CartesianVector &position2)
Get pair of closest positions for a pair of clusters.
void GetSortedKeyClusters(pandora::ClusterVector &sortedKeyClusters) const
Get a sorted vector of key clusters (view 1 clusters with current implementation) ...
unsigned int m_minRemnantClusterSize
The minimum hit number of a remnant cluster for it to be considered significant.
static float GetWireZPitch(const pandora::Pandora &pandora, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
pandora::StatusCode ProjectDeltaRayPositions(const MatrixType::Element &element, const pandora::HitType hitType, pandora::CartesianPointVector &projectedPositions) const
Use two views of a delta ray pfo to calculate projected positions in a given third view...
void UpdateUponDeletion(const pandora::Cluster *const pDeletedCluster)
Update to reflect cluster deletion.
HitTypeVector GetHitTypeVector()
Obtain the HitTypeVector of input views.
Header file for the geometry helper class.
bool PassElementChecks(const MatrixType::Element &element, const pandora::HitType hitType) const
Determine whether element satifies simple checks.
float m_distanceToLine
The maximum perpendicular distance of a position to a line for it to be considered close...
bool IsInLineSegment(const pandora::CartesianVector &lowerBoundary, const pandora::CartesianVector &upperBoundary, const pandora::CartesianVector &point) const
Whether the projection of a given position lies on a defined line.
static bool SortHitsByPosition(const pandora::CaloHit *const pLhs, const pandora::CaloHit *const pRhs)
Sort calo hits by their position (use Z, followed by X, followed by Y)
Header file for the cluster helper class.
pandora::StatusCode GrowSeed(const MatrixType::Element &element, const pandora::HitType hitType, pandora::CaloHitList &collectedHits, pandora::CaloHitList &deltaRayRemantHits) const
Examine remaining hits in the delta ray cluster adding them to the delta ray seed or parent muon if a...
pandora::StatusCode ParameteriseMuon(const pandora::ParticleFlowObject *const pParentMuon, const pandora::Cluster *const pDeltaRayCluster, pandora::CartesianVector &positionOnMuon, pandora::CartesianVector &muonDirection) const
Parameterise the projection of a cosmic ray track in order to avoid poor/sparse projections.
Header file for the lar two dimensional sliding fit result class.
void ReclusterRemnant(const pandora::HitType hitType, const pandora::Cluster *const pMuonCluster, const pandora::Cluster *const pDeltaRayRemnant, pandora::ClusterVector &clusterVector, pandora::PfoVector &pfoVector) const
Reculster the remnant cluster, merging insignificant created clusters into the parent muon (if proxim...
bool IsBestElement(const MatrixType::Element &element, const pandora::HitType hitType, const MatrixType::ElementList &elementList) const
Determine whether the input element is the best to use to modify the contaminated cluster (best is de...
float m_maxDistanceToTrack
The minimum distance of an insignificant remnant cluster from the cosmic ray track for it to be merge...
TwoViewDeltaRayMatchingAlgorithm * m_pParentAlgorithm
Address of the parent matching algorithm.
float m_minContaminationLength
The minimum projected length of a delta ray cluster onto the muon track for it to be considered conta...
bool Run(TwoViewDeltaRayMatchingAlgorithm *const pAlgorithm, MatrixType &overlapMatrix)
Run the algorithm tool.
bool IsContaminated(const MatrixType::Element &element, const pandora::HitType hitType) const
Determine whether the cluster under investigation has muon contamination.
void GetGlobalDirection(const float dTdL, pandora::CartesianVector &direction) const
Get global direction coordinates for given sliding linear fit gradient.
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
const LayerFitResultMap & GetLayerFitResultMap() const
Get the layer fit result map.
bool IsCloseToLine(const pandora::CartesianVector &hitPosition, const pandora::CartesianVector &lineStart, const pandora::CartesianVector &lineEnd, const float distanceToLine) const
Whether a given position is close to a defined line.
fhicl::Table< sbnd::crt::CRTDetSimParams > Parameters
bool IsMuonEndpoint(const MatrixType::Element &element, const pandora::HitType hitType) const
Determine whether the matched clusters suggest that the delta ray is at the endpoint of the cosmic ra...
required by fuzzyCluster table::sbnd_g4_services gaushitTruthMatch pandora
std::vector< pandora::HitType > HitTypeVector
float m_minSeparation
The minimum delta ray - parent muon cluster separation required to investigate delta ray cluster...
void GetConnectedElements(const pandora::Cluster *const pCluster, const bool ignoreUnavailable, ElementList &elementList) const
Get a list of elements connected to a specified cluster.
float m_maxDistanceToHit
The maximum distance of a hit from the cosmic ray track for it to be added into the CR...
const std::string & GetClusterListName(const pandora::HitType hitType) const
Get the cluster list name corresponding to a specified hit type.
std::vector< art::Ptr< recob::Cluster > > ClusterVector
pandora::StatusCode GetMuonCluster(const pandora::PfoList &commonMuonPfoList, const pandora::HitType hitType, const pandora::Cluster *&pMuonCluster) const
Return the cluster of the common cosmic ray pfo in a given view (function demands there to be only on...
const std::string & GetClusteringAlgName() const
Get the name of the clustering algorithm to be used to recluster created delta ray remnants...
TwoDSlidingFitResult class.
BEGIN_PROLOG could also be cout
static float GetClosestDistance(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2)
Get closest distance between clusters in a pair of cluster lists.