OneViewDeltaRayMatchingAlgorithm.cc

Go to the documentation of this file.

#include "Pandora/AlgorithmHeaders.h"

#include "larpandoracontent/LArHelpers/LArClusterHelper.h"
#include "larpandoracontent/LArHelpers/LArPfoHelper.h"

#include "larpandoracontent/LArUtility/KDTreeLinkerAlgoT.h"

#include "larpandoracontent/LArThreeDReco/LArCosmicRay/OneViewDeltaRayMatchingAlgorithm.h"

using namespace pandora;

namespace lar_content
{

OneViewDeltaRayMatchingAlgorithm::OneViewDeltaRayMatchingAlgorithm() :
    m_overlapExtension(1.f),
    m_minClusterHits(3)
{
}

//------------------------------------------------------------------------------------------------------------------------------------------

StatusCode OneViewDeltaRayMatchingAlgorithm::Run()
{
    const PfoList muonPfoList(this->GetMuonPfoList());
    PfoList allPfoList(this->GetDeltaRayPfoList());

    allPfoList.insert(allPfoList.end(), muonPfoList.begin(), muonPfoList.end());

    m_deltaRayMatchingContainers.FillContainers(
        allPfoList, this->GetInputClusterList(TPC_VIEW_U), this->GetInputClusterList(TPC_VIEW_V), this->GetInputClusterList(TPC_VIEW_W));

    for (const HitType hitType : {TPC_VIEW_U, TPC_VIEW_V, TPC_VIEW_W})
        this->PerformOneViewMatching(hitType);

    for (const HitType hitType : {TPC_VIEW_U, TPC_VIEW_V, TPC_VIEW_W})
        this->PerformRecovery(hitType);

    m_deltaRayMatchingContainers.ClearContainers();

    return STATUS_CODE_SUCCESS;
}

//------------------------------------------------------------------------------------------------------------------------------------------

const ClusterList OneViewDeltaRayMatchingAlgorithm::GetInputClusterList(const HitType hitType)
{
    const std::string inputClusterListName((hitType == TPC_VIEW_U) ? m_inputClusterListNameU
            : (hitType == TPC_VIEW_V)                              ? m_inputClusterListNameV
                                                                   : m_inputClusterListNameW);

    const ClusterList *pInputClusterList(nullptr);

    PANDORA_THROW_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*this, inputClusterListName, pInputClusterList));

    if (!pInputClusterList)
        return ClusterList();

    return *pInputClusterList;
}

//------------------------------------------------------------------------------------------------------------------------------------------

const PfoList OneViewDeltaRayMatchingAlgorithm::GetMuonPfoList()
{
    const PfoList *pMuonPfoList(nullptr);

    PANDORA_THROW_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*this, m_muonPfoListName, pMuonPfoList));

    if (!pMuonPfoList)
        return PfoList();

    return *pMuonPfoList;
}

//------------------------------------------------------------------------------------------------------------------------------------------

const PfoList OneViewDeltaRayMatchingAlgorithm::GetDeltaRayPfoList()
{
    const PfoList *pDeltaRayPfoList(nullptr);

    PANDORA_THROW_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*this, m_deltaRayPfoListName, pDeltaRayPfoList));

    if (!pDeltaRayPfoList)
        return PfoList();

    return *pDeltaRayPfoList;
}

//------------------------------------------------------------------------------------------------------------------------------------------

void OneViewDeltaRayMatchingAlgorithm::PerformOneViewMatching(const HitType hitType)
{
    const DeltaRayMatchingContainers::ClusterProximityMap &clusterProximityMap(m_deltaRayMatchingContainers.GetClusterProximityMap(hitType));
    const DeltaRayMatchingContainers::ClusterToPfoMap &clusterToPfoMap(m_deltaRayMatchingContainers.GetClusterToPfoMap(hitType));
    ClusterVector availableClusterList;

    for (auto &entry : clusterProximityMap)
    {
        if (entry.first->IsAvailable())
            availableClusterList.push_back(entry.first);
    }

    std::sort(availableClusterList.begin(), availableClusterList.end(), LArClusterHelper::SortByNHits);

    ClusterSet modifiedClusters;

    for (const Cluster *const pAvailableCluster : availableClusterList)
    {
        if (modifiedClusters.count(pAvailableCluster))
            continue;

        const DeltaRayMatchingContainers::ClusterProximityMap::const_iterator iter(clusterProximityMap.find(pAvailableCluster));

        // ATTN: Map will update during loop
        if (iter == clusterProximityMap.end())
            continue;

        bool found(false);
        const ClusterList &nearbyClusters(clusterProximityMap.at(pAvailableCluster));
        PfoVector nearbyMuonPfoVector;

        do
        {
            found = false;

            const ParticleFlowObject *pClosestMuonPfo(nullptr);
            float closestDistance(std::numeric_limits<float>::max());

            for (const Cluster *const pNearbyCluster : nearbyClusters)
            {
                if (!this->IsMuonPfo(pNearbyCluster))
                    continue;

                if (std::find(nearbyMuonPfoVector.begin(), nearbyMuonPfoVector.end(), clusterToPfoMap.at(pNearbyCluster)) !=
                    nearbyMuonPfoVector.end())
                    continue;

                found = true;

                const float separation(LArClusterHelper::GetClosestDistance(pNearbyCluster, pAvailableCluster));

                if (separation < closestDistance)
                {
                    closestDistance = separation;
                    pClosestMuonPfo = clusterToPfoMap.at(pNearbyCluster);
                }
            }

            if (pClosestMuonPfo)
                nearbyMuonPfoVector.push_back(pClosestMuonPfo);
        } while (found);

        if (nearbyMuonPfoVector.empty())
            continue;

        if (this->AddIntoExistingDeltaRay(pAvailableCluster, nearbyMuonPfoVector))
            continue;

        this->CreateDeltaRay(pAvailableCluster, nearbyMuonPfoVector, modifiedClusters);
    }
}

//------------------------------------------------------------------------------------------------------------------------------------------

bool OneViewDeltaRayMatchingAlgorithm::IsMuonPfo(const Cluster *const pCluster)
{
    const HitType hitType(LArClusterHelper::GetClusterHitType(pCluster));
    const DeltaRayMatchingContainers::ClusterToPfoMap &clusterToPfoMap(m_deltaRayMatchingContainers.GetClusterToPfoMap(hitType));
    const DeltaRayMatchingContainers::ClusterToPfoMap::const_iterator iter(clusterToPfoMap.find(pCluster));

    if (iter == clusterToPfoMap.end())
        return false;

    const ParticleFlowObject *const pPfo(iter->second);
    const PfoList &muonPfoList(this->GetMuonPfoList());

    return (std::find(muonPfoList.begin(), muonPfoList.end(), pPfo) != muonPfoList.end());
}

//------------------------------------------------------------------------------------------------------------------------------------------

bool OneViewDeltaRayMatchingAlgorithm::AddIntoExistingDeltaRay(const Cluster *const pAvailableCluster, const PfoVector &nearbyMuonPfoVector)
{
    const HitType hitType(LArClusterHelper::GetClusterHitType(pAvailableCluster));
    const HitType projectedHitType1((hitType == TPC_VIEW_U) ? TPC_VIEW_V : (hitType == TPC_VIEW_V) ? TPC_VIEW_W : TPC_VIEW_U);
    const HitType projectedHitType2((projectedHitType1 == TPC_VIEW_U) ? TPC_VIEW_V
            : (projectedHitType1 == TPC_VIEW_V)                       ? TPC_VIEW_W
                                                                      : TPC_VIEW_U);
    const DeltaRayMatchingContainers::ClusterToPfoMap &clusterToPfoMap1(m_deltaRayMatchingContainers.GetClusterToPfoMap(projectedHitType1));
    const DeltaRayMatchingContainers::ClusterToPfoMap &clusterToPfoMap2(m_deltaRayMatchingContainers.GetClusterToPfoMap(projectedHitType2));

    for (const ParticleFlowObject *const pNearbyMuonPfo : nearbyMuonPfoVector)
    {
        const Cluster *const pProjectedCluster1(this->GetBestProjectedCluster({pAvailableCluster}, pNearbyMuonPfo, projectedHitType1, false));
        const Cluster *const pProjectedCluster2(this->GetBestProjectedCluster({pAvailableCluster}, pNearbyMuonPfo, projectedHitType2, false));

        if ((!pProjectedCluster1) || (!pProjectedCluster2))
            continue;

        const ParticleFlowObject *const pPfo1(clusterToPfoMap1.at(pProjectedCluster1));
        const ParticleFlowObject *const pPfo2(clusterToPfoMap2.at(pProjectedCluster2));

        if (pPfo1 == pPfo2)
        {
            PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToPfo(*this, pPfo1, pAvailableCluster));

            m_deltaRayMatchingContainers.AddClustersToPfoMaps(pPfo1);

            return true;
        }
    }

    return false;
}

//------------------------------------------------------------------------------------------------------------------------------------------

const Cluster *OneViewDeltaRayMatchingAlgorithm::GetBestProjectedCluster(
    const ClusterList &deltaRayClusterGroup, const ParticleFlowObject *const pNearbyMuonPfo, const HitType hitType, const bool findAvailable)
{
    ClusterList muonClusterList;
    LArPfoHelper::GetClusters(pNearbyMuonPfo, hitType, muonClusterList);

    if (muonClusterList.size() != 1)
        return nullptr;

    const DeltaRayMatchingContainers::ClusterProximityMap &clusterProximityMap(m_deltaRayMatchingContainers.GetClusterProximityMap(hitType));
    auto muonProximityIter(clusterProximityMap.find(muonClusterList.front()));

    if (muonProximityIter == clusterProximityMap.end())
        return nullptr;

    float spanMinX(0.f), spanMaxX(0.f);
    this->GetClusterSpanX(deltaRayClusterGroup, spanMinX, spanMaxX);

    unsigned int highestHit(0);
    const Cluster *pProjectedCluster(nullptr);

    for (const Cluster *const pNearbyCluster : muonProximityIter->second)
    {
        if (findAvailable && !pNearbyCluster->IsAvailable())
            continue;

        if (!findAvailable && !this->IsDeltaRayPfo(pNearbyCluster))
            continue;

        float minX(0.f), maxX(0.f);
        pNearbyCluster->GetClusterSpanX(minX, maxX);

        if ((maxX < (spanMinX - m_overlapExtension)) || (minX > (spanMaxX + m_overlapExtension)))
            continue;

        if (pNearbyCluster->GetNCaloHits() > highestHit)
        {
            highestHit = pNearbyCluster->GetNCaloHits();
            pProjectedCluster = pNearbyCluster;
        }
    }

    return pProjectedCluster;
}

//------------------------------------------------------------------------------------------------------------------------------------------

bool OneViewDeltaRayMatchingAlgorithm::IsDeltaRayPfo(const Cluster *const pCluster)
{
    const HitType hitType(LArClusterHelper::GetClusterHitType(pCluster));
    const DeltaRayMatchingContainers::ClusterToPfoMap &clusterToPfoMap(m_deltaRayMatchingContainers.GetClusterToPfoMap(hitType));

    const DeltaRayMatchingContainers::ClusterToPfoMap::const_iterator iter(clusterToPfoMap.find(pCluster));

    if (iter == clusterToPfoMap.end())
        return false;

    const ParticleFlowObject *const pDeltaRayPfo(iter->second);
    const PfoList &deltaRayPfoList(this->GetDeltaRayPfoList());

    return (std::find(deltaRayPfoList.begin(), deltaRayPfoList.end(), pDeltaRayPfo) != deltaRayPfoList.end());
}

//------------------------------------------------------------------------------------------------------------------------------------------

void OneViewDeltaRayMatchingAlgorithm::GetClusterSpanX(const ClusterList &clusterList, float &spanMinX, float &spanMaxX)
{
    spanMinX = std::numeric_limits<float>::max();
    spanMaxX = -std::numeric_limits<float>::max();

    for (const Cluster *const pCluster : clusterList)
    {
        float minX(0.f), maxX(0.f);
        pCluster->GetClusterSpanX(minX, maxX);

        if (minX < spanMinX)
            spanMinX = minX;

        if (maxX > spanMaxX)
            spanMaxX = maxX;
    }
}

//------------------------------------------------------------------------------------------------------------------------------------------

void OneViewDeltaRayMatchingAlgorithm::CreateDeltaRay(const Cluster *const pAvailableCluster, const PfoVector &nearbyMuonPfoVector, ClusterSet &modifiedClusters)
{
    ClusterList clusterGroup, consideredClusters;
    this->GetNearbyAvailableClusters(pAvailableCluster, consideredClusters, clusterGroup);

    for (const Cluster *const pModifiedCluster : clusterGroup)
        modifiedClusters.insert(pModifiedCluster);

    const HitType hitType(LArClusterHelper::GetClusterHitType(pAvailableCluster));
    const HitType projectedHitType1((hitType == TPC_VIEW_U) ? TPC_VIEW_V : (hitType == TPC_VIEW_V) ? TPC_VIEW_W : TPC_VIEW_U);
    const HitType projectedHitType2((projectedHitType1 == TPC_VIEW_U) ? TPC_VIEW_V
            : (projectedHitType1 == TPC_VIEW_V)                       ? TPC_VIEW_W
                                                                      : TPC_VIEW_U);
    ClusterList projectedClusters1, projectedClusters2;

    for (const ParticleFlowObject *const pNearbyMuonPfo : nearbyMuonPfoVector)
    {
        const Cluster *const pProjectedCluster1(this->GetBestProjectedCluster(clusterGroup, pNearbyMuonPfo, projectedHitType1, true));
        const Cluster *const pProjectedCluster2(this->GetBestProjectedCluster(clusterGroup, pNearbyMuonPfo, projectedHitType2, true));

        if ((!pProjectedCluster1) && (!pProjectedCluster2))
            continue;

        ClusterList consideredClusters1;
        if (pProjectedCluster1)
            this->GetNearbyAvailableClusters(pProjectedCluster1, consideredClusters1, projectedClusters1);

        ClusterList consideredClusters2;
        if (pProjectedCluster2)
            this->GetNearbyAvailableClusters(pProjectedCluster2, consideredClusters2, projectedClusters2);
    }

    const Cluster *const pCluster1(this->MergeClusterGroup(clusterGroup));
    const Cluster *const pCluster2(this->MergeClusterGroup(projectedClusters1));
    const Cluster *const pCluster3(this->MergeClusterGroup(projectedClusters2));

    this->CreatePfo(pCluster1, pCluster2, pCluster3);
}

//------------------------------------------------------------------------------------------------------------------------------------------

void OneViewDeltaRayMatchingAlgorithm::GetNearbyAvailableClusters(const Cluster *const pCluster, ClusterList &consideredClusters, ClusterList &foundClusters)
{
    const HitType hitType(LArClusterHelper::GetClusterHitType(pCluster));
    const DeltaRayMatchingContainers::ClusterProximityMap &clusterProximityMap(m_deltaRayMatchingContainers.GetClusterProximityMap(hitType));

    consideredClusters.push_back(pCluster);

    if (!pCluster->IsAvailable())
        return;

    if (std::find(foundClusters.begin(), foundClusters.end(), pCluster) == foundClusters.end())
        foundClusters.push_back(pCluster);

    const DeltaRayMatchingContainers::ClusterProximityMap::const_iterator proximityIter(clusterProximityMap.find(pCluster));

    if (proximityIter == clusterProximityMap.end())
        return;

    for (const Cluster *const pNearbyCluster : proximityIter->second)
    {
        if (!pNearbyCluster->IsAvailable())
            continue;

        if (std::find(consideredClusters.begin(), consideredClusters.end(), pNearbyCluster) != consideredClusters.end())
            continue;

        if (std::find(foundClusters.begin(), foundClusters.end(), pNearbyCluster) != foundClusters.end())
            continue;

        this->GetNearbyAvailableClusters(pNearbyCluster, consideredClusters, foundClusters);
    }
}

//------------------------------------------------------------------------------------------------------------------------------------------

const Cluster *OneViewDeltaRayMatchingAlgorithm::MergeClusterGroup(const ClusterList &clusterGroup)
{
    if (clusterGroup.empty())
        return nullptr;

    const Cluster *const pClusterToEnlarge(clusterGroup.front());

    if (clusterGroup.size() == 1)
        return pClusterToEnlarge;

    const HitType hitType(LArClusterHelper::GetClusterHitType(pClusterToEnlarge));
    const std::string inputClusterListName((hitType == TPC_VIEW_U) ? m_inputClusterListNameU
            : (hitType == TPC_VIEW_V)                              ? m_inputClusterListNameV
                                                                   : m_inputClusterListNameW);

    for (const Cluster *const pClusterToDelete : clusterGroup)
    {
        m_deltaRayMatchingContainers.RemoveClusterFromContainers(pClusterToDelete);

        if (pClusterToDelete != pClusterToEnlarge)
        {
            PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ReplaceCurrentList<Cluster>(*this, inputClusterListName));
            PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::MergeAndDeleteClusters(*this, pClusterToEnlarge, pClusterToDelete));
        }
    }

    m_deltaRayMatchingContainers.AddClustersToContainers({pClusterToEnlarge}, {nullptr});

    return pClusterToEnlarge;
}

//------------------------------------------------------------------------------------------------------------------------------------------

void OneViewDeltaRayMatchingAlgorithm::CreatePfo(const Cluster *const pCluster1, const Cluster *const pCluster2, const Cluster *const pCluster3)
{
    const PfoList *pPfoList(nullptr);
    std::string pfoListName;
    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*this, pPfoList, pfoListName));

    PandoraContentApi::ParticleFlowObject::Parameters pfoParameters;
    pfoParameters.m_particleId = E_MINUS;
    pfoParameters.m_charge = PdgTable::GetParticleCharge(E_MINUS);
    pfoParameters.m_mass = PdgTable::GetParticleMass(E_MINUS);
    pfoParameters.m_energy = 0.f;
    pfoParameters.m_momentum = CartesianVector(0.f, 0.f, 0.f);

    if (pCluster1)
        pfoParameters.m_clusterList.push_back(pCluster1);

    if (pCluster2)
        pfoParameters.m_clusterList.push_back(pCluster2);

    if (pCluster3)
        pfoParameters.m_clusterList.push_back(pCluster3);

    if (pfoParameters.m_clusterList.empty())
        throw StatusCodeException(STATUS_CODE_FAILURE);

    const ParticleFlowObject *pPfo(nullptr);

    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*this, pfoParameters, pPfo));
    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*this, m_outputPfoListName));
    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ReplaceCurrentList<Pfo>(*this, m_outputPfoListName));

    m_deltaRayMatchingContainers.AddClustersToPfoMaps(pPfo);
}

//------------------------------------------------------------------------------------------------------------------------------------------

void OneViewDeltaRayMatchingAlgorithm::PerformRecovery(const HitType hitType)
{
    const ClusterList &inputClusterList(this->GetInputClusterList(hitType));

    for (const Cluster *const pCluster : inputClusterList)
    {
        if (!pCluster->IsAvailable())
            continue;

        if (pCluster->GetNCaloHits() < m_minClusterHits)
            continue;

        this->CreatePfo(pCluster, nullptr, nullptr);
    }
}

//------------------------------------------------------------------------------------------------------------------------------------------

StatusCode OneViewDeltaRayMatchingAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
{
    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "MuonPfoListName", m_muonPfoListName));

    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "DeltaRayPfoListName", m_deltaRayPfoListName));

    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "InputClusterListNameU", m_inputClusterListNameU));

    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "InputClusterListNameV", m_inputClusterListNameV));

    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "InputClusterListNameW", m_inputClusterListNameW));

    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "OutputPfoListName", m_outputPfoListName));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
        XmlHelper::ReadValue(xmlHandle, "SearchRegion1D", m_deltaRayMatchingContainers.m_searchRegion1D));

    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "OverlapExtension", m_overlapExtension));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinClusterHits", m_minClusterHits));

    return STATUS_CODE_SUCCESS;
}

} // namespace lar_content