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.