9 #include "Pandora/AlgorithmHeaders.h" 23 CosmicRayRemovalTool::CosmicRayRemovalTool() :
24 m_slidingFitWindow(10000),
25 m_minContaminationLength(3.
f),
26 m_maxDistanceToHit(1.
f),
27 m_minRemnantClusterSize(3),
28 m_maxDistanceToTrack(2.
f)
38 if (PandoraContentApi::GetSettings(*m_pParentAlgorithm)->ShouldDisplayAlgorithmInfo())
39 std::cout <<
"----> Running Algorithm Tool: " << this->GetInstanceName() <<
", " << this->GetType() << std::endl;
41 bool changesMade(
false);
46 ClusterSet usedKeyClusters;
47 for (
const Cluster *
const pKeyCluster : sortedKeyClusters)
49 if (usedKeyClusters.count(pKeyCluster))
55 for (
const TensorType::Element &element : elementList)
56 usedKeyClusters.insert(element.GetClusterU());
60 changesMade = (changesMade ? changesMade : changesMadeInIteration);
70 ClusterSet modifiedClusters, checkedClusters;
72 for (
const TensorType::Element &element : elementList)
74 for (
const HitType hitType : {TPC_VIEW_U, TPC_VIEW_V, TPC_VIEW_W})
76 const Cluster *
const pDeltaRayCluster(element.GetCluster(hitType));
78 if (checkedClusters.count(pDeltaRayCluster))
82 if ((modifiedClusters.count(element.GetClusterU())) || (modifiedClusters.count(element.GetClusterV())) ||
83 (modifiedClusters.count(element.GetClusterW())))
92 if (!this->
IsBestElement(element, hitType, elementList, modifiedClusters))
95 checkedClusters.insert(pDeltaRayCluster);
97 CaloHitList deltaRayHits;
98 this->
CreateSeed(element, hitType, deltaRayHits);
100 if (deltaRayHits.empty())
104 CaloHitList remnantHits;
105 if (this->
GrowSeed(element, hitType, deltaRayHits, remnantHits) != STATUS_CODE_SUCCESS)
109 if (deltaRayHits.size() == pDeltaRayCluster->GetNCaloHits())
112 modifiedClusters.insert(pDeltaRayCluster);
118 return !modifiedClusters.empty();
136 const Cluster *pMuonCluster(
nullptr), *
const pDeltaRayCluster(element.GetCluster(hitType));
144 CartesianVector muonDirection(0.
f, 0.
f, 0.
f);
147 CartesianVector deltaRayVertex(0.
f, 0.
f, 0.
f), muonVertex(0.
f, 0.
f, 0.
f);
151 float furthestSeparation(0.
f);
152 CartesianVector extendedPoint(0.
f, 0.
f, 0.
f);
154 CaloHitList deltaRayHitList;
155 pDeltaRayCluster->GetOrderedCaloHitList().FillCaloHitList(deltaRayHitList);
157 for (
const CaloHit *
const pCaloHit : deltaRayHitList)
159 const CartesianVector &position(pCaloHit->GetPositionVector());
160 const float separation((position - muonVertex).GetMagnitude());
162 if (separation > furthestSeparation)
166 furthestSeparation = separation;
167 extendedPoint = position;
177 CaloHitList muonHitList;
178 pMuonCluster->GetOrderedCaloHitList().FillCaloHitList(muonHitList);
180 for (
const CaloHit *
const pCaloHit : muonHitList)
182 if (this->
IsInLineSegment(deltaRayVertex, extendedPoint, pCaloHit->GetPositionVector()))
194 CartesianVector positionOnMuon(0.
f, 0.
f, 0.
f), muonDirection(0.
f, 0.
f, 0.
f);
196 positionOnMuon, muonDirection) != STATUS_CODE_SUCCESS)
199 CartesianPointVector deltaRayProjectedPositions;
203 CaloHitList deltaRayHitList;
204 element.GetCluster(hitType)->GetOrderedCaloHitList().FillCaloHitList(deltaRayHitList);
206 CaloHitVector deltaRayHitVector(deltaRayHitList.begin(), deltaRayHitList.end());
209 for (
const CaloHit *
const pCaloHit : deltaRayHitVector)
211 const CartesianVector &position(pCaloHit->GetPositionVector());
213 for (
const CartesianVector &projectedPosition : deltaRayProjectedPositions)
215 const float distanceToProjectionSquared((position - projectedPosition).GetMagnitudeSquared());
223 const float distanceToMuonHitsSquared(muonDirection.GetCrossProduct(position - positionOnMuon).GetMagnitudeSquared());
225 if (distanceToMuonHitsSquared < m_maxDistanceToHit * m_maxDistanceToHit)
228 collectedHits.push_back(pCaloHit);
240 const Cluster *
const pDeltaRayCluster(element.GetCluster(hitType)), *pMuonCluster(
nullptr);
243 return STATUS_CODE_NOT_FOUND;
246 CartesianVector positionOnMuon(0.
f, 0.
f, 0.
f), muonDirection(0.
f, 0.
f, 0.
f);
248 muonDirection) != STATUS_CODE_SUCCESS)
249 return STATUS_CODE_NOT_FOUND;
257 return STATUS_CODE_SUCCESS;
263 const Cluster *
const pDeltaRayCluster,
const float &minDistanceFromMuon,
const bool demandCloserToCollected,
264 const CaloHitList &protectedHits, CaloHitList &collectedHits)
const 266 CaloHitList deltaRayHitList;
267 pDeltaRayCluster->GetOrderedCaloHitList().FillCaloHitList(deltaRayHitList);
269 bool hitsAdded(
false);
275 for (
const CaloHit *
const pCaloHit : deltaRayHitList)
277 if (std::find(protectedHits.begin(), protectedHits.end(), pCaloHit) != protectedHits.end())
280 if (std::find(collectedHits.begin(), collectedHits.end(), pCaloHit) != collectedHits.end())
283 const float distanceToCollectedHits(collectedHits.empty()
284 ? std::numeric_limits<float>::max()
286 const float distanceToMuonHits(muonDirection.GetCrossProduct(pCaloHit->GetPositionVector() - positionOnMuon).GetMagnitude());
288 if ((distanceToMuonHits < minDistanceFromMuon) || (demandCloserToCollected && (distanceToCollectedHits > distanceToMuonHits)))
291 collectedHits.push_back(pCaloHit);
300 const TensorType::Element &element,
const HitType hitType, CaloHitList &collectedHits, CaloHitList &deltaRayRemnantHits)
const 302 const Cluster *
const pDeltaRayCluster(element.GetCluster(hitType)), *pMuonCluster(
nullptr);
310 std::string originalListName, fragmentListName;
311 ClusterList originalClusterList(1, pDeltaRayCluster);
313 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
315 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
316 PandoraContentApi::InitializeFragmentation(*
m_pParentAlgorithm, originalClusterList, originalListName, fragmentListName));
318 CaloHitList deltaRayHitList;
319 pDeltaRayCluster->GetOrderedCaloHitList().FillCaloHitList(deltaRayHitList);
321 const Cluster *pDeltaRay(
nullptr), *pDeltaRayRemnant(
nullptr);
323 for (
const CaloHit *
const pCaloHit : deltaRayHitList)
325 const bool isDeltaRay(std::find(collectedHits.begin(), collectedHits.end(), pCaloHit) != collectedHits.end());
326 const bool isDeltaRayRemnant(
327 !isDeltaRay && (std::find(deltaRayRemnantHits.begin(), deltaRayRemnantHits.end(), pCaloHit) != deltaRayRemnantHits.end()));
328 const Cluster *&pCluster(isDeltaRay ? pDeltaRay : isDeltaRayRemnant ? pDeltaRayRemnant : pMuonCluster);
332 PandoraContentApi::Cluster::Parameters parameters;
333 parameters.m_caloHitList.push_back(pCaloHit);
334 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Cluster::Create(*
m_pParentAlgorithm, parameters, pCluster));
338 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToCluster(*
m_pParentAlgorithm, pCluster, pCaloHit));
342 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::EndFragmentation(*
m_pParentAlgorithm, fragmentListName, originalListName));
346 if (pDeltaRayRemnant)
347 this->
ReclusterRemnant(hitType, pMuonCluster, pDeltaRayRemnant, clusterVector, pfoVector);
349 clusterVector.push_back(pMuonCluster);
350 pfoVector.push_back(element.GetOverlapResult().GetCommonMuonPfoList().front());
351 clusterVector.push_back(pDeltaRay);
352 pfoVector.push_back(
nullptr);
362 std::string caloHitListName(hitType == TPC_VIEW_U ?
"CaloHitListU" : hitType == TPC_VIEW_V ?
"CaloHitListV" :
"CaloHitListW");
364 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ReplaceCurrentList<CaloHit>(*
m_pParentAlgorithm, caloHitListName));
365 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
367 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Delete<Cluster>(*
m_pParentAlgorithm, pDeltaRayRemnant));
369 const ClusterList *pClusterList(
nullptr);
370 std::string newClusterListName;
371 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
374 const ClusterList remnantClusters(pClusterList->begin(), pClusterList->end());
376 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
378 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=,
381 for (
const Cluster *
const pRemnant : remnantClusters)
387 PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::MergeAndDeleteClusters(*
m_pParentAlgorithm, pMuonCluster, pRemnant));
392 clusterVector.push_back(pRemnant);
393 pfoVector.push_back(
nullptr);
401 PANDORA_RETURN_RESULT_IF_AND_IF(
402 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"SlidingFitWindow",
m_slidingFitWindow));
404 PANDORA_RETURN_RESULT_IF_AND_IF(
405 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinContaminationLength",
m_minContaminationLength));
407 PANDORA_RETURN_RESULT_IF_AND_IF(
408 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MaxDistanceToHit",
m_maxDistanceToHit));
410 PANDORA_RETURN_RESULT_IF_AND_IF(
411 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MinRemnantClusterSize",
m_minRemnantClusterSize));
413 PANDORA_RETURN_RESULT_IF_AND_IF(
414 STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
"MaxDistanceToTrack",
m_maxDistanceToTrack));
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.
void GetConnectedElements(const pandora::Cluster *const pCluster, const bool ignoreUnavailable, ElementList &elementList) const
Get a list of elements connected to a specified cluster.
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...
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
float m_distanceToLine
The maximum perpendicular distance of a position to a line for it to be considered close...
bool RemoveCosmicRayHits(const TensorType::ElementList &elementList) const
Remove hits from delta ray clusters that belong to the parent muon.
ThreeViewDeltaRayMatchingAlgorithm * m_pParentAlgorithm
Address of the parent matching algorithm.
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.
float m_maxDistanceToHit
The maximum distance of a hit from the cosmic ray track for it to be added into the CR...
void UpdateUponDeletion(const pandora::Cluster *const pDeletedCluster)
Update to reflect cluster deletion.
Header file for the geometry helper class.
float m_maxDistanceToTrack
The minimum distance of an insignificant remnant cluster from the cosmic ray track for it to be merge...
bool IsContaminated(const TensorType::Element &element, const pandora::HitType hitType) const
Determine whether the cluster under investigation has muon contamination.
std::string GetClusteringAlgName() const
Get the name of the clustering algorithm to be used to recluster created delta ray remnants...
std::vector< Element > ElementList
bool IsBestElement(const TensorType::Element &element, const pandora::HitType hitType, const TensorType::ElementList &elementList, const pandora::ClusterSet &modifiedClusters) const
Determine whether the input element is the best to use to modify the contaminated cluster (best is de...
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.
virtual bool PassElementChecks(const TensorType::Element &element, const pandora::HitType hitType) const
Determine whether element satifies simple checks.
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.
virtual bool PassElementChecks(const TensorType::Element &element, const pandora::HitType hitType) const =0
Determine whether element satifies simple checks.
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.
ThreeViewDeltaRayMatchingAlgorithm class.
static float GetWirePitch(const pandora::Pandora &pandora, const pandora::HitType view, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
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.
void GetSortedKeyClusters(pandora::ClusterVector &sortedKeyClusters) const
Get a sorted vector of key clusters (U clusters with current implementation)
pandora::StatusCode ProjectDeltaRayPositions(const TensorType::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 GetGlobalDirection(const float dTdL, pandora::CartesianVector &direction) const
Get global direction coordinates for given sliding linear fit gradient.
const LayerFitResultMap & GetLayerFitResultMap() const
Get the layer fit result map.
void SplitDeltaRayCluster(const TensorType::Element &element, const pandora::HitType hitType, pandora::CaloHitList &collectedHits, pandora::CaloHitList &deltaRayRemnantHits) const
Fragment the delta ray cluster adding hits to the muon track and perhaps creating significant remnant...
pandora::StatusCode GrowSeed(const TensorType::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...
unsigned int m_minRemnantClusterSize
The minimum hit number of a remnant cluster for it to be considered significant.
bool Run(ThreeViewDeltaRayMatchingAlgorithm *const pAlgorithm, TensorType &overlapTensor)
Run the algorithm tool.
void ReclusterRemnant(const pandora::HitType hitType, const pandora::Cluster *const pMuonCluster, const pandora::Cluster *const pDeltaRayRemnant, pandora::ClusterVector &clusterVector, pandora::PfoVector &pfoVector) const
Reculster a given remnant cluster, merging insignificant created clusters into the parent muon (if pr...
const std::string & GetClusterListName(const pandora::HitType hitType) const
Get the cluster list name corresponding to a specified hit type.
float m_minContaminationLength
The minimum projected length of a delta ray cluster onto the muon track for it to be considered conta...
std::vector< art::Ptr< recob::Cluster > > ClusterVector
bool IsMuonEndpoint(const TensorType::Element &element, const bool ignoreHitType, const pandora::HitType hitTypeToIgnore=pandora::TPC_VIEW_U) const
Determine whether the matched clusters suggest that the delta ray is at the endpoint of the cosmic ra...
void CreateSeed(const TensorType::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...
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...
TwoDSlidingFitResult class.
static float GetClosestDistance(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2)
Get closest distance between clusters in a pair of cluster lists.
virtual pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)=0