TransverseExtensionAlgorithm.cc

Go to the documentation of this file.

 
#include "Pandora/AlgorithmHeaders.h"
 
#include "larpandoracontent/LArHelpers/LArClusterHelper.h"
#include "larpandoracontent/LArHelpers/LArPointingClusterHelper.h"
 
#include "larpandoracontent/LArTwoDReco/LArClusterAssociation/TransverseExtensionAlgorithm.h"
 
using namespace pandora;
 
namespace lar_content
{
 
TransverseExtensionAlgorithm::TransverseExtensionAlgorithm() :
    m_minClusterLength(5.f),
    m_maxLongitudinalDisplacement(10.f),
    m_maxTransverseDisplacement(1.f)
{
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void TransverseExtensionAlgorithm::GetListOfCleanClusters(const ClusterList *const pClusterList, ClusterVector &clusterVector) const
{
    for (ClusterList::const_iterator iter = pClusterList->begin(), iterEnd = pClusterList->end(); iter != iterEnd; ++iter)
        clusterVector.push_back(*iter);
 
    std::sort(clusterVector.begin(), clusterVector.end(), LArClusterHelper::SortByNHits);
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void TransverseExtensionAlgorithm::FillClusterAssociationMatrix(const ClusterVector &clusterVector, ClusterAssociationMatrix &clusterAssociationMatrix) const
{
    // Convert long clusters into pointing clusters
    LArPointingClusterList pointingClusterList;
 
    for (ClusterVector::const_iterator iter = clusterVector.begin(), iterEnd = clusterVector.end(); iter != iterEnd; ++iter)
    {
        const Cluster *const pCluster(*iter);
 
        if (LArClusterHelper::GetLengthSquared(pCluster) < m_minClusterLength * m_minClusterLength)
            continue;
 
        try
        {
            pointingClusterList.push_back(LArPointingCluster(*iter));
        }
        catch (StatusCodeException &)
        {
        }
    }
 
    // Form associations between clusters
    for (LArPointingClusterList::const_iterator iter1 = pointingClusterList.begin(), iterEnd1 = pointingClusterList.end(); iter1 != iterEnd1; ++iter1)
    {
        const LArPointingCluster &parentCluster = *iter1;
 
        for (ClusterVector::const_iterator iter2 = clusterVector.begin(), iterEnd2 = clusterVector.end(); iter2 != iterEnd2; ++iter2)
        {
            const Cluster *const pDaughterCluster(*iter2);
 
            if (parentCluster.GetCluster() == pDaughterCluster)
                continue;
 
            this->FillClusterAssociationMatrix(parentCluster, pDaughterCluster, clusterAssociationMatrix);
        }
    }
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void TransverseExtensionAlgorithm::FillClusterAssociationMatrix(const LArPointingCluster &pointingCluster,
    const Cluster *const pDaughterCluster, ClusterAssociationMatrix &clusterAssociationMatrix) const
{
    const Cluster *const pParentCluster(pointingCluster.GetCluster());
 
    if (pParentCluster == pDaughterCluster)
        return;
 
    for (unsigned int useInner = 0; useInner < 2; ++useInner)
    {
        const LArPointingCluster::Vertex &pointingVertex(useInner == 1 ? pointingCluster.GetInnerVertex() : pointingCluster.GetOuterVertex());
        const ClusterAssociation::VertexType vertexType(useInner == 1 ? ClusterAssociation::INNER : ClusterAssociation::OUTER);
 
        if (pointingVertex.GetRms() > 0.5f)
            continue;
 
        const float projectedDisplacement(LArClusterHelper::GetClosestDistance(pointingVertex.GetPosition(), pDaughterCluster));
 
        if (projectedDisplacement > m_maxLongitudinalDisplacement)
            continue;
 
        ClusterAssociation::AssociationType associationType(ClusterAssociation::WEAK);
        float figureOfMerit(projectedDisplacement);
 
        CartesianVector firstCoordinate(0.f, 0.f, 0.f);
        CartesianVector secondCoordinate(0.f, 0.f, 0.f);
        LArClusterHelper::GetExtremalCoordinates(pDaughterCluster, firstCoordinate, secondCoordinate);
 
        float firstL(0.f), firstT(0.f), secondT(0.f), secondL(0.f);
        LArPointingClusterHelper::GetImpactParameters(pointingVertex, firstCoordinate, firstL, firstT);
        LArPointingClusterHelper::GetImpactParameters(pointingVertex, secondCoordinate, secondL, secondT);
 
        const float innerL(firstL < secondL ? firstL : secondL);
        const float innerT(firstL < secondL ? firstT : secondT);
        const float outerL(firstL > secondL ? firstL : secondL);
        const float outerT(firstL > secondL ? firstT : secondT);
 
        if (innerL > 0.f && innerL < 2.5f && outerL < m_maxLongitudinalDisplacement && innerT < m_maxTransverseDisplacement &&
            outerT < 1.5f * m_maxTransverseDisplacement)
        {
            associationType = ClusterAssociation::STRONG;
            figureOfMerit = outerL;
        }
 
        (void)clusterAssociationMatrix[pParentCluster].insert(
            ClusterAssociationMap::value_type(pDaughterCluster, ClusterAssociation(vertexType, vertexType, associationType, figureOfMerit)));
    }
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void TransverseExtensionAlgorithm::FillClusterMergeMap(const ClusterAssociationMatrix &parentToDaughterMatrix, ClusterMergeMap &clusterMergeMap) const
{
    ClusterAssociationMatrix daughterToParentMatrix;
 
    // Loop over parent clusters and select nearby daughter clusters that are closer than another parent cluster
    ClusterVector sortedParentClusters;
    for (const auto &mapEntry : parentToDaughterMatrix)
        sortedParentClusters.push_back(mapEntry.first);
    std::sort(sortedParentClusters.begin(), sortedParentClusters.end(), LArClusterHelper::SortByNHits);
 
    for (const Cluster *const pParentCluster : sortedParentClusters)
    {
        const ClusterAssociationMap &daughterToAssociationMap(parentToDaughterMatrix.at(pParentCluster));
 
        float maxDisplacementInner(std::numeric_limits<float>::max());
        float maxDisplacementOuter(std::numeric_limits<float>::max());
 
        // Find the nearest parent cluster
        ClusterVector sortedLocalDaughterClusters;
        for (const auto &mapEntry : daughterToAssociationMap)
            sortedLocalDaughterClusters.push_back(mapEntry.first);
        std::sort(sortedLocalDaughterClusters.begin(), sortedLocalDaughterClusters.end(), LArClusterHelper::SortByNHits);
 
        for (const Cluster *const pDaughterCluster : sortedLocalDaughterClusters)
        {
            const ClusterAssociation &clusterAssociation(daughterToAssociationMap.at(pDaughterCluster));
 
            if (clusterAssociation.GetAssociation() == ClusterAssociation::WEAK)
            {
                if (clusterAssociation.GetParent() == ClusterAssociation::INNER && clusterAssociation.GetFigureOfMerit() < maxDisplacementInner)
                    maxDisplacementInner = clusterAssociation.GetFigureOfMerit();
 
                if (clusterAssociation.GetParent() == ClusterAssociation::OUTER && clusterAssociation.GetFigureOfMerit() < maxDisplacementOuter)
                    maxDisplacementOuter = clusterAssociation.GetFigureOfMerit();
            }
        }
 
        // Select daughter clusters that are closer than the nearest parent cluster
        for (const Cluster *const pDaughterCluster : sortedLocalDaughterClusters)
        {
            const ClusterAssociation &clusterAssociation(daughterToAssociationMap.at(pDaughterCluster));
 
            if (clusterAssociation.GetAssociation() == ClusterAssociation::STRONG)
            {
                if (clusterAssociation.GetParent() == ClusterAssociation::INNER && clusterAssociation.GetFigureOfMerit() < maxDisplacementInner)
                    (void)daughterToParentMatrix[pDaughterCluster].insert(ClusterAssociationMap::value_type(pParentCluster, clusterAssociation));
 
                if (clusterAssociation.GetParent() == ClusterAssociation::OUTER && clusterAssociation.GetFigureOfMerit() < maxDisplacementOuter)
                    (void)daughterToParentMatrix[pDaughterCluster].insert(ClusterAssociationMap::value_type(pParentCluster, clusterAssociation));
            }
        }
    }
 
    // Loop over daughter clusters and select the nearest parent clusters
    ClusterVector sortedDaughterClusters;
    for (const auto &mapEntry : daughterToParentMatrix)
        sortedDaughterClusters.push_back(mapEntry.first);
    std::sort(sortedDaughterClusters.begin(), sortedDaughterClusters.end(), LArClusterHelper::SortByNHits);
 
    // Loop over parent clusters and select nearby daughter clusters that are closer than another parent cluster
    for (const Cluster *const pDaughterCluster : sortedDaughterClusters)
    {
        const ClusterAssociationMap &parentToAssociationMap(daughterToParentMatrix.at(pDaughterCluster));
 
        const Cluster *pParentCluster(nullptr);
        float minDisplacement(std::numeric_limits<float>::max());
 
        ClusterVector sortedLocalParentClusters;
        for (const auto &mapEntry : parentToAssociationMap)
            sortedLocalParentClusters.push_back(mapEntry.first);
        std::sort(sortedLocalParentClusters.begin(), sortedLocalParentClusters.end(), LArClusterHelper::SortByNHits);
 
        for (const Cluster *const pCandidateParentCluster : sortedLocalParentClusters)
        {
            const ClusterAssociation &clusterAssociation(parentToAssociationMap.at(pCandidateParentCluster));
 
            if (clusterAssociation.GetFigureOfMerit() < minDisplacement)
            {
                minDisplacement = clusterAssociation.GetFigureOfMerit();
                pParentCluster = pCandidateParentCluster;
            }
        }
 
        if (pParentCluster)
        {
            ClusterList &parentList(clusterMergeMap[pParentCluster]);
 
            if (parentList.end() == std::find(parentList.begin(), parentList.end(), pDaughterCluster))
                parentList.push_back(pDaughterCluster);
 
            ClusterList &daughterList(clusterMergeMap[pDaughterCluster]);
 
            if (daughterList.end() == std::find(daughterList.begin(), daughterList.end(), pParentCluster))
                daughterList.push_back(pParentCluster);
        }
    }
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode TransverseExtensionAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
{
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinClusterLength", m_minClusterLength));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
        XmlHelper::ReadValue(xmlHandle, "MaxLongitudinalDisplacement", m_maxLongitudinalDisplacement));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
        XmlHelper::ReadValue(xmlHandle, "MaxTransverseDisplacement", m_maxTransverseDisplacement));
 
    return ClusterExtensionAlgorithm::ReadSettings(xmlHandle);
}
 
} // namespace lar_content