DeltaRayParentAlgorithm.cc

Go to the documentation of this file.

 
#include "Pandora/AlgorithmHeaders.h"
 
#include "larpandoracontent/LArHelpers/LArClusterHelper.h"
#include "larpandoracontent/LArHelpers/LArPfoHelper.h"
 
#include "larpandoracontent/LArThreeDReco/LArCosmicRay/DeltaRayParentAlgorithm.h"
 
using namespace pandora;
 
namespace lar_content
{
 
DeltaRayParentAlgorithm::DeltaRayParentAlgorithm() : m_distanceForMatching(5.f)
{
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode DeltaRayParentAlgorithm::Run()
{
    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 || pMuonPfoList->empty())
        return STATUS_CODE_SUCCESS;
 
    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 || pDeltaRayPfoList->empty())
        return STATUS_CODE_SUCCESS;
 
    PfoLengthMap pfoLengthMap;
    this->InitialisePfoLengthMap(pMuonPfoList, pDeltaRayPfoList, pfoLengthMap);
 
    PfoVector deltaRayPfoVector(pDeltaRayPfoList->begin(), pDeltaRayPfoList->end());
    std::sort(deltaRayPfoVector.begin(), deltaRayPfoVector.end(), LArPfoHelper::SortByNHits);
 
    for (const ParticleFlowObject *const pPfo : deltaRayPfoVector)
    {
        if (!pPfo->GetParentPfoList().empty())
            continue;
 
        const ParticleFlowObject *pParentPfo(nullptr);
        this->FindParentPfo(pfoLengthMap, pPfo, pParentPfo);
 
        if (!pParentPfo)
            continue;
 
        this->AssignToParentPfo(pMuonPfoList, pDeltaRayPfoList, pPfo, pParentPfo, pfoLengthMap);
    }
 
    return STATUS_CODE_SUCCESS;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void DeltaRayParentAlgorithm::InitialisePfoLengthMap(const PfoList *const muonPfoList, const PfoList *const deltaRayPfoList, PfoLengthMap &pfoLengthMap) const
{
    for (const ParticleFlowObject *const pPfo : *muonPfoList)
        pfoLengthMap[pPfo] = LArPfoHelper::GetTwoDLengthSquared(pPfo);
 
    for (const ParticleFlowObject *const pPfo : *deltaRayPfoList)
        pfoLengthMap[pPfo] = LArPfoHelper::GetTwoDLengthSquared(pPfo);
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void DeltaRayParentAlgorithm::FindParentPfo(const PfoLengthMap &pfoLengthMap, const ParticleFlowObject *const pPfo, const ParticleFlowObject *&pParentPfo) const
{
    const PfoLengthMap::const_iterator currentIter(pfoLengthMap.find(pPfo));
 
    if (currentIter == pfoLengthMap.end())
        throw StatusCodeException(STATUS_CODE_FAILURE);
 
    PfoVector allPfoVector;
 
    for (auto &entry : pfoLengthMap)
        allPfoVector.push_back(entry.first);
 
    const float lengthSquared(currentIter->second);
    float bestDistance(m_distanceForMatching);
 
    for (const ParticleFlowObject *const pTestParent : allPfoVector)
    {
        if (pTestParent == pPfo)
            continue;
 
        const PfoLengthMap::const_iterator testIter(pfoLengthMap.find(pTestParent));
 
        if (testIter == pfoLengthMap.end())
            throw StatusCodeException(STATUS_CODE_FAILURE);
 
        if (testIter->second < lengthSquared)
            continue;
 
        float distance(std::numeric_limits<float>::max());
 
        if (this->GetTwoDSeparation(pPfo, pTestParent, distance) == STATUS_CODE_NOT_FOUND)
            continue;
 
        if (distance < bestDistance)
        {
            pParentPfo = pTestParent;
            bestDistance = distance;
        }
    }
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode DeltaRayParentAlgorithm::GetTwoDSeparation(const ParticleFlowObject *const pPfo1, const ParticleFlowObject *const pPfo2, float &separation) const
{
    ClusterList clusterListU1, clusterListV1, clusterListW1;
    ClusterList clusterListU2, clusterListV2, clusterListW2;
 
    LArPfoHelper::GetClusters(pPfo1, TPC_VIEW_U, clusterListU1);
    LArPfoHelper::GetClusters(pPfo1, TPC_VIEW_V, clusterListV1);
    LArPfoHelper::GetClusters(pPfo1, TPC_VIEW_W, clusterListW1);
 
    LArPfoHelper::GetClusters(pPfo2, TPC_VIEW_U, clusterListU2);
    LArPfoHelper::GetClusters(pPfo2, TPC_VIEW_V, clusterListV2);
    LArPfoHelper::GetClusters(pPfo2, TPC_VIEW_W, clusterListW2);
 
    float numViews(0.f), distance(0.f);
 
    if (!clusterListU1.empty() && !clusterListU2.empty())
    {
        distance += LArClusterHelper::GetClosestDistance(clusterListU1, clusterListU2);
        numViews += 1.f;
    }
 
    if (!clusterListV1.empty() && !clusterListV2.empty())
    {
        distance += LArClusterHelper::GetClosestDistance(clusterListV1, clusterListV2);
        numViews += 1.f;
    }
 
    if (!clusterListW1.empty() && !clusterListW2.empty())
    {
        distance += LArClusterHelper::GetClosestDistance(clusterListW1, clusterListW2);
        numViews += 1.f;
    }
 
    if (numViews < std::numeric_limits<float>::epsilon())
        return STATUS_CODE_NOT_FOUND;
 
    separation = distance / numViews;
 
    return STATUS_CODE_SUCCESS;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void DeltaRayParentAlgorithm::AssignToParentPfo(const PfoList *const pMuonPfoList, const PfoList *const pDeltaRayPfoList,
    const ParticleFlowObject *const pPfo, const ParticleFlowObject *const pParentPfo, PfoLengthMap &pfoLengthMap) const
{
    if (std::find(pMuonPfoList->begin(), pMuonPfoList->end(), pParentPfo) != pMuonPfoList->end())
        PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SetPfoParentDaughterRelationship(*this, pParentPfo, pPfo));
 
    if (std::find(pDeltaRayPfoList->begin(), pDeltaRayPfoList->end(), pParentPfo) != pDeltaRayPfoList->end())
    {
        ClusterList pfoClusters;
        LArPfoHelper::GetTwoDClusterList(pPfo, pfoClusters);
 
        PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Delete(*this, pPfo, m_deltaRayPfoListName));
 
        for (const Cluster *const pCluster : pfoClusters)
            PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToPfo(*this, pParentPfo, pCluster));
 
        this->UpdatePfoLengthMap({pPfo, pParentPfo}, pParentPfo, pfoLengthMap);
    }
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void DeltaRayParentAlgorithm::UpdatePfoLengthMap(const PfoList &pfosToRemove, const ParticleFlowObject *const pPfoToAdd, PfoLengthMap &pfoLengthMap) const
{
    for (const ParticleFlowObject *const pPfoToRemove : pfosToRemove)
    {
        const PfoLengthMap::const_iterator iter(pfoLengthMap.find(pPfoToRemove));
 
        if (iter == pfoLengthMap.end())
            throw StatusCodeException(STATUS_CODE_FAILURE);
 
        pfoLengthMap.erase(iter);
    }
 
    pfoLengthMap[pPfoToAdd] = LArPfoHelper::GetTwoDLengthSquared(pPfoToAdd);
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode DeltaRayParentAlgorithm::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_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "DistanceForMatching", m_distanceForMatching));
 
    return STATUS_CODE_SUCCESS;
}
 
} // namespace lar_content