DeltaRayRemovalTool.cc

Go to the documentation of this file.

 
#include "Pandora/AlgorithmHeaders.h"
 
#include "larpandoracontent/LArHelpers/LArClusterHelper.h"
#include "larpandoracontent/LArHelpers/LArGeometryHelper.h"
 
#include "larpandoracontent/LArObjects/LArTwoDSlidingFitResult.h"
 
#include "larpandoracontent/LArThreeDReco/LArCosmicRay/DeltaRayRemovalTool.h"
 
using namespace pandora;
 
namespace lar_content
{
 
DeltaRayRemovalTool::DeltaRayRemovalTool() :
    m_slidingFitWindow(10000),
    m_minDeviationFromTransverse(0.35f),
    m_contaminationWindow(5.f),
    m_significantHitThreshold(3),
    m_minDistanceFromMuon(1.f),
    m_maxDistanceToCollected(1.f)
{
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool DeltaRayRemovalTool::Run(ThreeViewDeltaRayMatchingAlgorithm *const pAlgorithm, TensorType &overlapTensor)
{
    m_pParentAlgorithm = pAlgorithm;
 
    if (PandoraContentApi::GetSettings(*m_pParentAlgorithm)->ShouldDisplayAlgorithmInfo())
        std::cout << "----> Running Algorithm Tool: " << this->GetInstanceName() << ", " << this->GetType() << std::endl;
 
    bool changesMade(false);
 
    ClusterVector sortedKeyClusters;
    overlapTensor.GetSortedKeyClusters(sortedKeyClusters);
 
    ClusterSet usedKeyClusters;
    for (const Cluster *const pKeyCluster : sortedKeyClusters)
    {
        if (usedKeyClusters.count(pKeyCluster))
            continue;
 
        TensorType::ElementList elementList;
        overlapTensor.GetConnectedElements(pKeyCluster, true, elementList);
 
        for (const TensorType::Element &element : elementList)
            usedKeyClusters.insert(element.GetClusterU());
 
        const bool changesMadeInIteration = this->RemoveDeltaRayHits(elementList);
 
        changesMade = (changesMade ? changesMade : changesMadeInIteration);
    }
 
    return changesMade;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool DeltaRayRemovalTool::RemoveDeltaRayHits(const TensorType::ElementList &elementList) const
{
    ClusterSet modifiedClusters, checkedClusters;
 
    for (const TensorType::Element &element : elementList)
    {
        for (const HitType hitType : {TPC_VIEW_U, TPC_VIEW_V, TPC_VIEW_W})
        {
            const Cluster *pDeltaRayCluster(element.GetCluster(hitType));
            const ParticleFlowObject *const pMuonPfo(element.GetOverlapResult().GetCommonMuonPfoList().front());
 
            if (checkedClusters.count(pDeltaRayCluster))
                continue;
 
            // ATTN: The underlying tensor will update during this loop, do not proceed if element has been modified
            if ((modifiedClusters.count(element.GetClusterU())) || (modifiedClusters.count(element.GetClusterV())) ||
                (modifiedClusters.count(element.GetClusterW())))
                continue;
 
            if (!this->PassElementChecks(element, hitType))
                continue;
 
            if (!this->IsContaminated(element, hitType))
                continue;
 
            if (!this->IsBestElement(element, hitType, elementList, modifiedClusters))
                continue;
 
            checkedClusters.insert(pDeltaRayCluster);
 
            CaloHitList deltaRayHits;
            if (m_pParentAlgorithm->CollectHitsFromMuon(nullptr, nullptr, pDeltaRayCluster, pMuonPfo, m_minDistanceFromMuon,
                    m_maxDistanceToCollected, deltaRayHits) != STATUS_CODE_SUCCESS)
                continue;
 
            modifiedClusters.insert(pDeltaRayCluster);
 
            this->SplitMuonCluster(element, hitType, deltaRayHits);
        }
    }
 
    return !modifiedClusters.empty();
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool DeltaRayRemovalTool::PassElementChecks(const TensorType::Element &element, const HitType hitType) const
{
    // ATTN: Avoid endpoints, topological michel reconstruction is very ambiguous
    if (this->IsMuonEndpoint(element, false))
        return false;
 
    return RemovalBaseTool::PassElementChecks(element, hitType);
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool DeltaRayRemovalTool::IsContaminated(const TensorType::Element &element, const HitType hitType) const
{
    const Cluster *pMuonCluster(nullptr), *const pDeltaRayCluster(element.GetCluster(hitType));
 
    if (m_pParentAlgorithm->GetMuonCluster(element.GetOverlapResult().GetCommonMuonPfoList(), hitType, pMuonCluster) != STATUS_CODE_SUCCESS)
        return false;
 
    const float slidingFitPitch(LArGeometryHelper::GetWireZPitch(this->GetPandora()));
    const TwoDSlidingFitResult slidingFitResult(pMuonCluster, m_slidingFitWindow, slidingFitPitch);
 
    CartesianVector muonDirection(0.f, 0.f, 0.f);
    slidingFitResult.GetGlobalDirection(slidingFitResult.GetLayerFitResultMap().begin()->second.GetGradient(), muonDirection);
 
    const CartesianVector xAxis(1.f, 0.f, 0.f);
    const bool isTransverse((muonDirection.GetOpeningAngle(xAxis) < m_minDeviationFromTransverse) ||
                            ((M_PI - muonDirection.GetOpeningAngle(xAxis)) < m_minDeviationFromTransverse));
 
    if (isTransverse)
        return false;
 
    CartesianVector deltaRayVertex(0.f, 0.f, 0.f), muonVertex(0.f, 0.f, 0.f);
    LArClusterHelper::GetClosestPositions(pDeltaRayCluster, pMuonCluster, deltaRayVertex, muonVertex);
 
    CaloHitList minusMuonHits, minusDeltaRayHits, plusMuonHits, plusDeltaRayHits;
    const CartesianVector minusPosition(muonVertex - (muonDirection * m_contaminationWindow));
    const CartesianVector plusPosition(muonVertex + (muonDirection * m_contaminationWindow));
 
    this->FindExtrapolatedHits(pMuonCluster, muonVertex, minusPosition, minusMuonHits);
    this->FindExtrapolatedHits(pMuonCluster, muonVertex, plusPosition, plusMuonHits);
    this->FindExtrapolatedHits(pDeltaRayCluster, muonVertex, minusPosition, minusDeltaRayHits);
    this->FindExtrapolatedHits(pDeltaRayCluster, muonVertex, plusPosition, plusDeltaRayHits);
 
    if ((minusMuonHits.size() < m_significantHitThreshold) && (plusMuonHits.size() < m_significantHitThreshold))
        return true;
 
    if ((minusMuonHits.size() < m_significantHitThreshold) && (minusDeltaRayHits.size() >= m_significantHitThreshold) &&
        (plusDeltaRayHits.size() < m_significantHitThreshold) && (plusMuonHits.size() >= m_significantHitThreshold))
    {
        return true;
    }
 
    if ((plusMuonHits.size() < m_significantHitThreshold) && (plusDeltaRayHits.size() >= m_significantHitThreshold) &&
        (minusDeltaRayHits.size() < m_significantHitThreshold) && (minusMuonHits.size() >= m_significantHitThreshold))
    {
        return true;
    }
 
    return false;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void DeltaRayRemovalTool::SplitMuonCluster(const TensorType::Element &element, const HitType hitType, const CaloHitList &deltaRayHits) const
{
    const Cluster *pDeltaRayCluster(element.GetCluster(hitType)), *pMuonCluster(nullptr);
 
    if (m_pParentAlgorithm->GetMuonCluster(element.GetOverlapResult().GetCommonMuonPfoList(), hitType, pMuonCluster) != STATUS_CODE_SUCCESS)
        throw StatusCodeException(STATUS_CODE_FAILURE);
 
    m_pParentAlgorithm->UpdateUponDeletion(pMuonCluster);
    m_pParentAlgorithm->UpdateUponDeletion(pDeltaRayCluster);
 
    m_pParentAlgorithm->SplitMuonCluster(m_pParentAlgorithm->GetClusterListName(hitType), pMuonCluster, deltaRayHits, pDeltaRayCluster);
 
    ClusterVector clusterVector;
    PfoVector pfoVector;
    clusterVector.push_back(pMuonCluster);
    pfoVector.push_back(element.GetOverlapResult().GetCommonMuonPfoList().front());
    clusterVector.push_back(pDeltaRayCluster);
    pfoVector.push_back(nullptr);
    m_pParentAlgorithm->UpdateForNewClusters(clusterVector, pfoVector);
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode DeltaRayRemovalTool::ReadSettings(const TiXmlHandle xmlHandle)
{
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SlidingFitWindow", m_slidingFitWindow));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
        XmlHelper::ReadValue(xmlHandle, "MinDeviationFromTransverse", m_minDeviationFromTransverse));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "ContaminationWindow", m_contaminationWindow));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
        XmlHelper::ReadValue(xmlHandle, "SignificantHitThreshold", m_significantHitThreshold));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinDistanceFromMuon", m_minDistanceFromMuon));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxDistanceToCollected", m_maxDistanceToCollected));
 
    return STATUS_CODE_SUCCESS;
}
 
} // namespace lar_content