BranchGrowingAlgorithm.cc

Go to the documentation of this file.

#include "Pandora/AlgorithmHeaders.h"

#include "larpandoracontent/LArHelpers/LArClusterHelper.h"

#include "larpandoracontent/LArTrackShowerId/BranchGrowingAlgorithm.h"

using namespace pandora;

namespace lar_content
{

void BranchGrowingAlgorithm::FindAssociatedClusters(const Cluster *const pParticleSeed, ClusterVector &candidateClusters,
    ClusterUsageMap &forwardUsageMap, ClusterUsageMap &backwardUsageMap) const
{
    ClusterVector currentSeedAssociations, newSeedAssociations;
    currentSeedAssociations.push_back(pParticleSeed);

    unsigned int associationOrder(1);

    while (!currentSeedAssociations.empty())
    {
        for (ClusterVector::iterator iterI = candidateClusters.begin(), iterIEnd = candidateClusters.end(); iterI != iterIEnd; ++iterI)
        {
            const Cluster *const pCandidateCluster = *iterI;

            if (NULL == pCandidateCluster)
                continue;

            for (ClusterVector::iterator iterJ = currentSeedAssociations.begin(), iterJEnd = currentSeedAssociations.end(); iterJ != iterJEnd; ++iterJ)
            {
                const Cluster *const pAssociatedCluster = *iterJ;

                const AssociationType associationType(this->AreClustersAssociated(pAssociatedCluster, pCandidateCluster));

                if (NONE == associationType)
                    continue;

                // Check we store best association between this seed and candidate
                Association association(associationOrder, associationType);
                const Association &existingAssociation = forwardUsageMap[pParticleSeed][pCandidateCluster];

                if (association.GetType() > existingAssociation.GetType())
                {
                    // If not first association, check strength of previous association in chain
                    if (pParticleSeed != pAssociatedCluster)
                        association.SetType(std::min(association.GetType(), backwardUsageMap[pAssociatedCluster][pParticleSeed].GetType()));

                    forwardUsageMap[pParticleSeed][pCandidateCluster] = association;
                    backwardUsageMap[pCandidateCluster][pParticleSeed] = association;
                }

                newSeedAssociations.push_back(pCandidateCluster);
                *iterI = NULL;
            }
        }

        currentSeedAssociations = newSeedAssociations;
        newSeedAssociations.clear();
        ++associationOrder;
    }
}

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

void BranchGrowingAlgorithm::IdentifyClusterMerges(const ClusterVector &particleSeedVector, const ClusterUsageMap &backwardUsageMap,
    SeedAssociationList &seedAssociationList) const
{
    ClusterVector sortedCandidates;
    for (const auto &mapEntry : backwardUsageMap) sortedCandidates.push_back(mapEntry.first);
    std::sort(sortedCandidates.begin(), sortedCandidates.end(), LArClusterHelper::SortByNHits);

    for (const Cluster *const pCluster : sortedCandidates)
    {
        const ClusterAssociationMap &particleSeedUsageMap(backwardUsageMap.at(pCluster));

        if (particleSeedUsageMap.empty())
            throw StatusCodeException(STATUS_CODE_FAILURE);

        ClusterVector sortedSeeds;
        for (const auto &mapEntry : particleSeedUsageMap) sortedSeeds.push_back(mapEntry.first);
        std::sort(sortedSeeds.begin(), sortedSeeds.end(), LArClusterHelper::SortByNHits);

        const Cluster *pBestParticleSeed = NULL;
        AssociationType bestType(NONE);
        unsigned int bestOrder(std::numeric_limits<unsigned int>::max());

        for (const Cluster *const pParticleSeed : sortedSeeds)
        {
            const Association &association(particleSeedUsageMap.at(pParticleSeed));

            if ((association.GetType() > bestType) || ((association.GetType() == bestType) && (association.GetOrder() < bestOrder)))
            {
                // Break-out condition for single order associations
                if ((SINGLE_ORDER == association.GetType()) && (association.GetOrder() > 1))
                    continue;

                // Type is primary consideration; order breaks ties
                pBestParticleSeed = pParticleSeed;
                bestType = association.GetType();
                bestOrder = association.GetOrder();
            }
            else if ((association.GetType() == bestType) && (association.GetOrder() == bestOrder))
            {
                // Remove ambiguous cluster from algorithm
                pBestParticleSeed = NULL;
            }
        }

        if (NULL == pBestParticleSeed)
            continue;

        seedAssociationList[pBestParticleSeed].push_back(pCluster);
    }

    // Now deal with seeds that have no associations
    for (ClusterVector::const_iterator iter = particleSeedVector.begin(), iterEnd = particleSeedVector.end(); iter != iterEnd; ++iter)
    {
        const Cluster *const pParticleSeed = *iter;

        if (seedAssociationList.end() == seedAssociationList.find(pParticleSeed))
            seedAssociationList[pParticleSeed] = ClusterVector();
    }
}

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

StatusCode BranchGrowingAlgorithm::ReadSettings(const TiXmlHandle /*xmlHandle*/)
{
    return STATUS_CODE_SUCCESS;
}

} // namespace lar_content