CosmicRayTrackRecoveryAlgorithm.cc

Go to the documentation of this file.

#include "Pandora/AlgorithmHeaders.h"

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

#include "larpandoracontent/LArHelpers/LArGeometryHelper.h"
#include "larpandoracontent/LArHelpers/LArClusterHelper.h"
#include "larpandoracontent/LArHelpers/LArPointingClusterHelper.h"

using namespace pandora;

namespace lar_content
{

CosmicRayTrackRecoveryAlgorithm::CosmicRayTrackRecoveryAlgorithm() :
    m_clusterMinLength(10.f),
    m_clusterMinSpanZ(2.f),
    m_clusterMinOverlapX(6.f),
    m_clusterMaxDeltaX(3.f)
{
}

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

StatusCode CosmicRayTrackRecoveryAlgorithm::Run()
{
    // Get the available clusters for each view
    ClusterVector availableClustersU, availableClustersV, availableClustersW;
    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetAvailableClusters(m_inputClusterListNameU, availableClustersU));
    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetAvailableClusters(m_inputClusterListNameV, availableClustersV));
    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->GetAvailableClusters(m_inputClusterListNameW, availableClustersW));

    // Select clean clusters in each view
    ClusterVector cleanClustersU, cleanClustersV, cleanClustersW;
    this->SelectCleanClusters(availableClustersU, cleanClustersU);
    this->SelectCleanClusters(availableClustersV, cleanClustersV);
    this->SelectCleanClusters(availableClustersW, cleanClustersW);

    // Calculate sliding fit results for clean clusters
    TwoDSlidingFitResultMap slidingFitResultMap;
    this->BuildSlidingFitResultMap(cleanClustersU, slidingFitResultMap);
    this->BuildSlidingFitResultMap(cleanClustersV, slidingFitResultMap);
    this->BuildSlidingFitResultMap(cleanClustersW, slidingFitResultMap);

    // Match clusters between pairs of views (using start/end information)
    ClusterAssociationMap matchedClustersUV, matchedClustersVW, matchedClustersWU;
    this->MatchViews(cleanClustersU, cleanClustersV, slidingFitResultMap, matchedClustersUV);
    this->MatchViews(cleanClustersV, cleanClustersW, slidingFitResultMap, matchedClustersVW);
    this->MatchViews(cleanClustersW, cleanClustersU, slidingFitResultMap, matchedClustersWU);

    // Create candidate particles using one, two and three primary views
    ParticleList candidateParticles;

    this->MatchThreeViews(cleanClustersU, cleanClustersV, cleanClustersW,
        matchedClustersUV, matchedClustersVW, matchedClustersWU, candidateParticles);
    
    this->MatchTwoViews(cleanClustersU, cleanClustersV, cleanClustersW,
        matchedClustersUV, matchedClustersVW, matchedClustersWU, candidateParticles);
   
    this->MatchOneView(cleanClustersU, cleanClustersV, cleanClustersW,
        matchedClustersUV, matchedClustersVW, matchedClustersWU, candidateParticles);
 
    // Build particle flow objects from candidate particles
    this->BuildParticles(candidateParticles);

    return STATUS_CODE_SUCCESS;
}

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

StatusCode CosmicRayTrackRecoveryAlgorithm::GetAvailableClusters(const std::string &inputClusterListName, ClusterVector &clusterVector) const
{
    const ClusterList *pClusterList = NULL;
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*this,
        inputClusterListName, pClusterList))

    if (!pClusterList || pClusterList->empty())
    { 
        if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
            std::cout << "CosmicRayTrackRecoveryAlgorithm: unable to find cluster list " << inputClusterListName << std::endl;

        return STATUS_CODE_SUCCESS;
    }

    for (const Cluster *const pCluster : *pClusterList)
    {
        if (PandoraContentApi::IsAvailable(*this, pCluster))
            clusterVector.push_back(pCluster);
    }

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

    return STATUS_CODE_SUCCESS;
}

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

void CosmicRayTrackRecoveryAlgorithm::SelectCleanClusters(const ClusterVector &inputVector, ClusterVector &outputVector) const
{
    for (ClusterVector::const_iterator iter = inputVector.begin(), iterEnd = inputVector.end(); iter != iterEnd; ++iter)
    {
        const Cluster *const pCluster = *iter;

        // Remove clusters below a minimum length
        if (LArClusterHelper::GetLengthSquared(pCluster) < m_clusterMinLength * m_clusterMinLength)
            continue;

        // Remove clusters nearly parallel to Z or X
        CartesianVector minCoordinate(0.f, 0.f, 0.f);
        CartesianVector maxCoordinate(0.f, 0.f, 0.f);
        LArClusterHelper::GetClusterBoundingBox(pCluster, minCoordinate, maxCoordinate);

        const CartesianVector deltaCoordinate(maxCoordinate - minCoordinate);
        if (std::fabs(deltaCoordinate.GetZ()) < m_clusterMinSpanZ || std::fabs(deltaCoordinate.GetX()) < m_clusterMinOverlapX)
            continue;

        outputVector.push_back(pCluster);
    }
}

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

void CosmicRayTrackRecoveryAlgorithm::BuildSlidingFitResultMap(const ClusterVector &clusterVector, TwoDSlidingFitResultMap &slidingFitResultMap) const
{
    const unsigned int m_halfWindowLayers(25);
    const float slidingFitPitch(LArGeometryHelper::GetWireZPitch(this->GetPandora()));

    for (ClusterVector::const_iterator iter = clusterVector.begin(), iterEnd = clusterVector.end(); iter != iterEnd; ++iter)
    {
        if (slidingFitResultMap.end() == slidingFitResultMap.find(*iter))
        {
            try
            {
                const TwoDSlidingFitResult slidingFitResult(*iter, m_halfWindowLayers, slidingFitPitch);

                if (!slidingFitResultMap.insert(TwoDSlidingFitResultMap::value_type(*iter, slidingFitResult)).second)
                    throw StatusCodeException(STATUS_CODE_FAILURE);
            }
            catch (StatusCodeException &statusCodeException)
            {
                if (STATUS_CODE_FAILURE == statusCodeException.GetStatusCode())
                    throw statusCodeException;
            }
        }
    }
}

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

void CosmicRayTrackRecoveryAlgorithm::MatchViews(const ClusterVector &clusterVector1, const ClusterVector &clusterVector2, 
    const TwoDSlidingFitResultMap &slidingFitResultMap, ClusterAssociationMap &clusterAssociationMap) const
{
    for (ClusterVector::const_iterator iter1 = clusterVector1.begin(), iterEnd1 = clusterVector1.end(); iter1 != iterEnd1; ++iter1)
        this->MatchClusters(*iter1, clusterVector2, slidingFitResultMap, clusterAssociationMap);

    for (ClusterVector::const_iterator iter2 = clusterVector2.begin(), iterEnd2 = clusterVector2.end(); iter2 != iterEnd2; ++iter2)
        this->MatchClusters(*iter2, clusterVector1, slidingFitResultMap, clusterAssociationMap);
}

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

void CosmicRayTrackRecoveryAlgorithm::MatchClusters(const Cluster* const pSeedCluster, const ClusterVector &targetClusters, 
    const TwoDSlidingFitResultMap &slidingFitResultMap, ClusterAssociationMap &clusterAssociationMap) const
{
    // Match seed cluster to target clusters according to alignment in X position of start/end positions
    // Two possible matches: (a) one-to-one associations where both the track start and end positions match up
    //                       (b) one-to-two associations where the track is split into two clusters in one view
    // Require overlap in X (according to clusterMinOverlapX) and alignment in X (according to clusterMaxDeltaX)

    TwoDSlidingFitResultMap::const_iterator fsIter = slidingFitResultMap.find(pSeedCluster);

    if (slidingFitResultMap.end() == fsIter)
        throw StatusCodeException(STATUS_CODE_FAILURE);

    const TwoDSlidingFitResult &slidingFitResult1(fsIter->second);
    const CartesianVector &innerVertex1(slidingFitResult1.GetGlobalMinLayerPosition());
    const CartesianVector &outerVertex1(slidingFitResult1.GetGlobalMaxLayerPosition());
    const float xSpan1(std::fabs(outerVertex1.GetX() - innerVertex1.GetX()));

    const Cluster *pBestClusterInner(NULL);
    const Cluster *pBestClusterOuter(NULL);
    const Cluster *pBestCluster(NULL);

    float bestDisplacementInner(m_clusterMaxDeltaX);
    float bestDisplacementOuter(m_clusterMaxDeltaX);
    float bestDisplacement(2.f * m_clusterMaxDeltaX);

    for (ClusterVector::const_iterator tIter = targetClusters.begin(), tIterEnd = targetClusters.end(); tIter != tIterEnd; ++tIter)
    {
        const Cluster *const pTargetCluster = *tIter;

        if (LArClusterHelper::GetClusterHitType(pSeedCluster) == LArClusterHelper::GetClusterHitType(pTargetCluster))
            throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);

        TwoDSlidingFitResultMap::const_iterator ftIter = slidingFitResultMap.find(*tIter);

        if (slidingFitResultMap.end() == ftIter)
            throw StatusCodeException(STATUS_CODE_FAILURE);

        const TwoDSlidingFitResult &slidingFitResult2(ftIter->second);
        const CartesianVector &innerVertex2(slidingFitResult2.GetGlobalMinLayerPosition());
        const CartesianVector &outerVertex2(slidingFitResult2.GetGlobalMaxLayerPosition());
        const float xSpan2(std::fabs(outerVertex2.GetX() - innerVertex2.GetX()));

        if (xSpan2 > 1.5f * xSpan1)
            continue;

        const float xMin1(std::min(innerVertex1.GetX(), outerVertex1.GetX()));
        const float xMax1(std::max(innerVertex1.GetX(), outerVertex1.GetX()));
        const float xMin2(std::min(innerVertex2.GetX(), outerVertex2.GetX()));
        const float xMax2(std::max(innerVertex2.GetX(), outerVertex2.GetX()));
        const float xOverlap(std::min(xMax1,xMax2) - std::max(xMin1,xMin2));

        if (xOverlap < m_clusterMinOverlapX)
            continue;
            
        const float dxMin(std::fabs(xMin2 - xMin1));
        const float dxMax(std::fabs(xMax2 - xMax1));
           
        if (dxMin < bestDisplacementInner)
        {
            pBestClusterInner = pTargetCluster;
            bestDisplacementInner = dxMin;
        }

        if (dxMax < bestDisplacementOuter)
        {
            pBestClusterOuter = pTargetCluster;
            bestDisplacementOuter = dxMax;
        }

        if (dxMin + dxMax < bestDisplacement)
        {
            pBestCluster = pTargetCluster;
            bestDisplacement = dxMin + dxMax;
        }
    }

    if (pBestCluster)
    {
        ClusterList &seedList(clusterAssociationMap[pSeedCluster]);

        if (seedList.end() == std::find(seedList.begin(), seedList.end(), pBestCluster))
            seedList.push_back(pBestCluster);

        ClusterList &bestList(clusterAssociationMap[pBestCluster]);

        if (bestList.end() == std::find(bestList.begin(), bestList.end(), pSeedCluster))
            bestList.push_back(pSeedCluster);
    }
    else if (pBestClusterInner && pBestClusterOuter)
    {
        TwoDSlidingFitResultMap::const_iterator iterInner = slidingFitResultMap.find(pBestClusterInner);  
        TwoDSlidingFitResultMap::const_iterator iterOuter = slidingFitResultMap.find(pBestClusterOuter);

        if (slidingFitResultMap.end() == iterInner || slidingFitResultMap.end() == iterOuter)
            throw StatusCodeException(STATUS_CODE_FAILURE);

        const LArPointingCluster pointingClusterInner(iterInner->second);
        const LArPointingCluster pointingClusterOuter(iterOuter->second);

        LArPointingCluster::Vertex pointingVertexInner, pointingVertexOuter;

        try
        {
            LArPointingClusterHelper::GetClosestVertices(pointingClusterInner, pointingClusterOuter, pointingVertexInner, pointingVertexOuter);
        }
        catch (StatusCodeException &)
        {
            return;
        }

        const LArPointingCluster::Vertex pointingEndInner(pointingVertexInner.IsInnerVertex() ? pointingClusterInner.GetOuterVertex() : pointingClusterInner.GetInnerVertex());
        const LArPointingCluster::Vertex pointingEndOuter(pointingVertexOuter.IsInnerVertex() ? pointingClusterOuter.GetOuterVertex() : pointingClusterOuter.GetInnerVertex());

        const float rSpan((pointingEndInner.GetPosition() - pointingEndOuter.GetPosition()).GetMagnitude());

        if (LArPointingClusterHelper::IsEmission(pointingVertexInner.GetPosition(), pointingVertexOuter, -1.f, 0.75f * rSpan, 5.f, 10.f) &&
            LArPointingClusterHelper::IsEmission(pointingVertexOuter.GetPosition(), pointingVertexInner, -1.f, 0.75f * rSpan, 5.f, 10.f))
        {
            ClusterList &bestInnerList(clusterAssociationMap[pBestClusterInner]);

            if (bestInnerList.end() == std::find(bestInnerList.begin(), bestInnerList.end(), pSeedCluster))
                bestInnerList.push_back(pSeedCluster);

            ClusterList &bestOuterList(clusterAssociationMap[pBestClusterOuter]);

            if (bestOuterList.end() == std::find(bestOuterList.begin(), bestOuterList.end(), pSeedCluster))
                bestOuterList.push_back(pSeedCluster);
        }
    }
}

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

void CosmicRayTrackRecoveryAlgorithm::MatchThreeViews(const ClusterVector &clusterVectorU,  const ClusterVector &clusterVectorV,
    const ClusterVector &clusterVectorW, const ClusterAssociationMap &matchedClustersUV, const ClusterAssociationMap &matchedClustersVW, 
    const ClusterAssociationMap &matchedClustersWU, ParticleList &particleList) const
{  
    ClusterSet vetoList;
    this->BuildVetoList(particleList, vetoList);

    ParticleList newParticleList;

    const ClusterVector &clusterVector1(clusterVectorU);
    const ClusterVector &clusterVector2(clusterVectorV);
    const ClusterVector &clusterVector3(clusterVectorW);
 
    const ClusterAssociationMap &matchedClusters12(matchedClustersUV);
    const ClusterAssociationMap &matchedClusters23(matchedClustersVW);
    const ClusterAssociationMap &matchedClusters31(matchedClustersWU);

    for (ClusterVector::const_iterator iter1 = clusterVector1.begin(), iterEnd1 = clusterVector1.end(); iter1 != iterEnd1; ++iter1)
    {
        const Cluster *const pCluster1 = *iter1;

        if (vetoList.count(pCluster1))
            continue;

        const ClusterAssociationMap::const_iterator iter311 = matchedClusters31.find(pCluster1);
        const ClusterList matchedClusters31_pCluster1(iter311 != matchedClusters31.end() ? iter311->second : ClusterList());

        const ClusterAssociationMap::const_iterator iter121 = matchedClusters12.find(pCluster1);
        const ClusterList matchedClusters12_pCluster1(iter121 != matchedClusters12.end() ? iter121->second : ClusterList());

        for (ClusterVector::const_iterator iter2 = clusterVector2.begin(), iterEndV = clusterVector2.end(); iter2 != iterEndV; ++iter2)
        { 
            const Cluster *const pCluster2 = *iter2;

            if (vetoList.count(pCluster2))
                continue;

            const ClusterAssociationMap::const_iterator iter122 = matchedClusters12.find(pCluster2);
            const ClusterList matchedClusters12_pCluster2(iter122 != matchedClusters12.end() ? iter122->second : ClusterList());

            const ClusterAssociationMap::const_iterator iter232 = matchedClusters23.find(pCluster2);
            const ClusterList matchedClusters23_pCluster2(iter232 != matchedClusters23.end() ? iter232->second : ClusterList());

            for (ClusterVector::const_iterator iter3 = clusterVector3.begin(), iterEnd3 = clusterVector3.end(); iter3 != iterEnd3; ++iter3)
            { 
                const Cluster *const pCluster3 = *iter3;

                if (vetoList.count(pCluster3))
                    continue;

                const ClusterAssociationMap::const_iterator iter233 = matchedClusters23.find(pCluster3);
                const ClusterList matchedClusters23_pCluster3(iter233 != matchedClusters23.end() ? iter233->second : ClusterList());

                const ClusterAssociationMap::const_iterator iter313 = matchedClusters31.find(pCluster3);
                const ClusterList matchedClusters31_pCluster3(iter313 != matchedClusters31.end() ? iter313->second : ClusterList());

                const bool match12((matchedClusters12_pCluster1.size() +  matchedClusters12_pCluster2.size() > 0) &&
                    ((matchedClusters12_pCluster1.size() == 1 && std::find(matchedClusters12_pCluster1.begin(), matchedClusters12_pCluster1.end(), pCluster2) != matchedClusters12_pCluster1.end()) || 
                     (matchedClusters12_pCluster1.size() == 0)) &&
                    ((matchedClusters12_pCluster2.size() == 1 && std::find(matchedClusters12_pCluster2.begin(), matchedClusters12_pCluster2.end(), pCluster1) != matchedClusters12_pCluster2.end()) || 
                     (matchedClusters12_pCluster2.size() == 0)));

                const bool match23((matchedClusters23_pCluster2.size() +  matchedClusters23_pCluster3.size() > 0) &&
                    ((matchedClusters23_pCluster2.size() == 1 && std::find(matchedClusters23_pCluster2.begin(), matchedClusters23_pCluster2.end(), pCluster3) != matchedClusters23_pCluster2.end()) || 
                     (matchedClusters23_pCluster2.size() == 0)) &&
                    ((matchedClusters23_pCluster3.size() == 1 && std::find(matchedClusters23_pCluster3.begin(), matchedClusters23_pCluster3.end(), pCluster2) != matchedClusters23_pCluster3.end()) || 
                    (matchedClusters23_pCluster3.size() == 0)));

                const bool match31((matchedClusters31_pCluster3.size() +  matchedClusters31_pCluster1.size() > 0) &&
                    ((matchedClusters31_pCluster3.size() == 1 && std::find(matchedClusters31_pCluster3.begin(), matchedClusters31_pCluster3.end(), pCluster1) != matchedClusters31_pCluster3.end()) || 
                    (matchedClusters31_pCluster3.size() == 0)) &&
                    ((matchedClusters31_pCluster1.size() == 1 && std::find(matchedClusters31_pCluster1.begin(), matchedClusters31_pCluster1.end(), pCluster3) != matchedClusters31_pCluster1.end()) || 
                    (matchedClusters31_pCluster1.size() == 0)));

                if (match12 && match23 && match31)
                {
                    Particle newParticle;
                    newParticle.m_clusterList.push_back(pCluster1);
                    newParticle.m_clusterList.push_back(pCluster2);
                    newParticle.m_clusterList.push_back(pCluster3);
                    newParticleList.push_back(newParticle);
                }
            }
        }
    }

    return this->RemoveAmbiguities(newParticleList, particleList);
}

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

void CosmicRayTrackRecoveryAlgorithm::MatchTwoViews(const ClusterVector &clusterVectorU,  const ClusterVector &clusterVectorV,
    const ClusterVector &clusterVectorW, const ClusterAssociationMap &matchedClustersUV, const ClusterAssociationMap &matchedClustersVW, 
    const ClusterAssociationMap &matchedClustersWU, ParticleList &particleList) const
{
    ClusterSet vetoList;
    this->BuildVetoList(particleList, vetoList);

    ParticleList newParticleList;

    for (unsigned int iView = 0; iView < 3; ++iView)
    {
        const ClusterVector &clusterVector1((0 == iView) ? clusterVectorU : (1 == iView) ? clusterVectorV : clusterVectorW);
        const ClusterVector &clusterVector2((0 == iView) ? clusterVectorV : (1 == iView) ? clusterVectorW : clusterVectorU);    
        const ClusterVector &clusterVector3((0 == iView) ? clusterVectorW : (1 == iView) ? clusterVectorU : clusterVectorV);
 
        const ClusterAssociationMap &matchedClusters12(((0 == iView) ? matchedClustersUV : (1 == iView) ? matchedClustersVW : matchedClustersWU));
        const ClusterAssociationMap &matchedClusters23(((0 == iView) ? matchedClustersVW : (1 == iView) ? matchedClustersWU : matchedClustersUV));
        const ClusterAssociationMap &matchedClusters31(((0 == iView) ? matchedClustersWU : (1 == iView) ? matchedClustersUV : matchedClustersVW));

        for (ClusterVector::const_iterator iter1 = clusterVector1.begin(), iterEnd1 = clusterVector1.end(); iter1 != iterEnd1; ++iter1)
        {
            const Cluster *const pCluster1 = *iter1;

            if (vetoList.count(pCluster1))
                continue;

            const ClusterAssociationMap::const_iterator iter311 = matchedClusters31.find(pCluster1);
            const ClusterList matchedClusters31_pCluster1(iter311 != matchedClusters31.end() ? iter311->second : ClusterList());

            const ClusterAssociationMap::const_iterator iter121 = matchedClusters12.find(pCluster1);
            const ClusterList matchedClusters12_pCluster1(iter121 != matchedClusters12.end() ? iter121->second : ClusterList());

            for (ClusterVector::const_iterator iter2 = clusterVector2.begin(), iterEnd2 = clusterVector2.end(); iter2 != iterEnd2; ++iter2)
            { 
                const Cluster *const pCluster2 = *iter2;

                if (vetoList.count(pCluster2))
                    continue;

                const ClusterAssociationMap::const_iterator iter122 = matchedClusters12.find(pCluster2);
                const ClusterList matchedClusters12_pCluster2(iter122 != matchedClusters12.end() ? iter122->second : ClusterList());

                const ClusterAssociationMap::const_iterator iter232 = matchedClusters23.find(pCluster2);
                const ClusterList matchedClusters23_pCluster2(iter232 != matchedClusters23.end() ? iter232->second : ClusterList());

                const bool match12((matchedClusters12_pCluster1.size() == 1 && std::find(matchedClusters12_pCluster1.begin(), matchedClusters12_pCluster1.end(), pCluster2) != matchedClusters12_pCluster1.end()) &&
                    (matchedClusters12_pCluster2.size() == 1 && std::find(matchedClusters12_pCluster2.begin(), matchedClusters12_pCluster2.end(), pCluster1) != matchedClusters12_pCluster2.end()) &&
                    (matchedClusters23_pCluster2.size() == 0 && matchedClusters31_pCluster1.size() == 0));

                if (!match12)
                continue;

                Particle newParticle;
                newParticle.m_clusterList.push_back(pCluster1);
                newParticle.m_clusterList.push_back(pCluster2);
                    
                for (ClusterVector::const_iterator iter3 = clusterVector3.begin(), iterEnd3 = clusterVector3.end(); iter3 != iterEnd3; ++iter3)
                { 
                    const Cluster *const pCluster3 = *iter3;
  
                    if (vetoList.count(pCluster3))
                        continue;
 
                    const ClusterAssociationMap::const_iterator iter233 = matchedClusters23.find(pCluster3);
                    const ClusterList matchedClusters23_pCluster3(iter233 != matchedClusters23.end() ? iter233->second : ClusterList());

                    const ClusterAssociationMap::const_iterator iter313 = matchedClusters31.find(pCluster3);
                    const ClusterList matchedClusters31_pCluster3(iter313 != matchedClusters31.end() ? iter313->second : ClusterList());

                    const bool match3((matchedClusters31_pCluster3.size() +  matchedClusters23_pCluster3.size() > 0) &&
                        ((matchedClusters31_pCluster3.size() == 1 && std::find(matchedClusters31_pCluster3.begin(), matchedClusters31_pCluster3.end(), pCluster1) != matchedClusters31_pCluster3.end()) || 
                        (matchedClusters31_pCluster3.size() == 0)) &&
                        ((matchedClusters23_pCluster3.size() == 1 && std::find(matchedClusters23_pCluster3.begin(), matchedClusters23_pCluster3.end(), pCluster2) != matchedClusters23_pCluster3.end()) || 
                        (matchedClusters23_pCluster3.size() == 0)));

                    if (match3)
                        newParticle.m_clusterList.push_back(pCluster3);
                }

                newParticleList.push_back(newParticle);
            }
        }
    }   

    return this->RemoveAmbiguities(newParticleList, particleList);
}

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

void CosmicRayTrackRecoveryAlgorithm::MatchOneView(const ClusterVector &clusterVectorU,  const ClusterVector &clusterVectorV,
    const ClusterVector &clusterVectorW, const ClusterAssociationMap &matchedClustersUV, const ClusterAssociationMap &matchedClustersVW, 
    const ClusterAssociationMap &matchedClustersWU, ParticleList &particleList) const
{
    ClusterSet vetoList;
    this->BuildVetoList(particleList, vetoList);

    ParticleList newParticleList;

    for (unsigned int iView = 0; iView < 3; ++iView)
    {
        const ClusterVector &clusterVector1((0 == iView) ? clusterVectorU : (1 == iView) ? clusterVectorV : clusterVectorW);
        const ClusterVector &clusterVector2((0 == iView) ? clusterVectorV : (1 == iView) ? clusterVectorW : clusterVectorU);
        const ClusterVector &clusterVector3((0 == iView) ? clusterVectorW : (1 == iView) ? clusterVectorU : clusterVectorV);
 
        const ClusterAssociationMap &matchedClusters12(((0 == iView) ? matchedClustersUV : (1 == iView) ? matchedClustersVW : matchedClustersWU));
        const ClusterAssociationMap &matchedClusters23(((0 == iView) ? matchedClustersVW : (1 == iView) ? matchedClustersWU : matchedClustersUV));
        const ClusterAssociationMap &matchedClusters31(((0 == iView) ? matchedClustersWU : (1 == iView) ? matchedClustersUV : matchedClustersVW));

        for (ClusterVector::const_iterator iter1 = clusterVector1.begin(), iterEnd1 = clusterVector1.end(); iter1 != iterEnd1; ++iter1)
        {
            const Cluster *const pCluster1 = *iter1;

            if (vetoList.count(pCluster1))
                continue;

            const ClusterAssociationMap::const_iterator iter311 = matchedClusters31.find(pCluster1);
            const ClusterList matchedClusters31_pCluster1(iter311 != matchedClusters31.end() ? iter311->second : ClusterList());

            const ClusterAssociationMap::const_iterator iter121 = matchedClusters12.find(pCluster1);
            const ClusterList matchedClusters12_pCluster1(iter121 != matchedClusters12.end() ? iter121->second : ClusterList());

            if (matchedClusters12_pCluster1.size() + matchedClusters31_pCluster1.size() > 0)
                continue;

            Particle newParticle;
            newParticle.m_clusterList.push_back(pCluster1);

            for (ClusterVector::const_iterator iter2 = clusterVector2.begin(), iterEnd2 = clusterVector2.end(); iter2 != iterEnd2; ++iter2)
            {
                const Cluster *const pCluster2 = *iter2;

                if (vetoList.count(pCluster2))
                    continue;

                const ClusterAssociationMap::const_iterator iter122 = matchedClusters12.find(pCluster2);
                const ClusterList matchedClusters12_pCluster2(iter122 != matchedClusters12.end() ? iter122->second : ClusterList());

                const ClusterAssociationMap::const_iterator iter232 = matchedClusters23.find(pCluster2);
                const ClusterList matchedClusters23_pCluster2(iter232 != matchedClusters23.end() ? iter232->second : ClusterList());

                if (matchedClusters12_pCluster2.size() == 1 && std::find(matchedClusters12_pCluster2.begin(), matchedClusters12_pCluster2.end(), pCluster1) != matchedClusters12_pCluster2.end() &&
                    matchedClusters23_pCluster2.size() == 0)
                    newParticle.m_clusterList.push_back(pCluster2);
            }

            for (ClusterVector::const_iterator iter3 = clusterVector3.begin(), iterEnd3 = clusterVector3.end(); iter3 != iterEnd3; ++iter3)
            {
                const Cluster *const pCluster3 = *iter3;
 
                if (vetoList.count(pCluster3))
                    continue;

                const ClusterAssociationMap::const_iterator iter233 = matchedClusters23.find(pCluster3);
                const ClusterList matchedClusters23_pCluster3(iter233 != matchedClusters23.end() ? iter233->second : ClusterList());

                const ClusterAssociationMap::const_iterator iter313 = matchedClusters31.find(pCluster3);
                const ClusterList matchedClusters31_pCluster3(iter313 != matchedClusters31.end() ? iter313->second : ClusterList());

                if (matchedClusters31_pCluster3.size() == 1 && std::find(matchedClusters31_pCluster3.begin(), matchedClusters31_pCluster3.end(), pCluster1) != matchedClusters31_pCluster3.end() &&
                    matchedClusters23_pCluster3.size() == 0)
                    newParticle.m_clusterList.push_back(pCluster3);
            }

            if (newParticle.m_clusterList.size() > 1)
                newParticleList.push_back(newParticle);
        }
    }

    return this->RemoveAmbiguities(newParticleList, particleList);
}
 
//------------------------------------------------------------------------------------------------------------------------------------------

void CosmicRayTrackRecoveryAlgorithm::BuildVetoList(const ParticleList &particleList, ClusterSet &vetoList) const
{
    for (ParticleList::const_iterator pIter = particleList.begin(), pIterEnd = particleList.end(); pIter != pIterEnd; ++pIter)
    {
        const Particle &particle = *pIter;
        vetoList.insert(particle.m_clusterList.begin(), particle.m_clusterList.end());
    }
}

//------------------------------------------------------------------------------------------------------------------------------------------
    
void CosmicRayTrackRecoveryAlgorithm::RemoveAmbiguities(const ParticleList &inputParticleList, ParticleList &outputParticleList) const
{
    for (ParticleList::const_iterator pIter1 = inputParticleList.begin(), pIterEnd1 = inputParticleList.end(); pIter1 != pIterEnd1; ++pIter1)
    {
        const Particle &particle1 = *pIter1;
        const ClusterList &clusterList1 = particle1.m_clusterList;

        try
        {
            bool isUnique(true);

            for (ParticleList::const_iterator pIter2 = pIter1, pIterEnd2 = inputParticleList.end(); pIter2 != pIterEnd2; ++pIter2)
            {
                const Particle &particle2 = *pIter2;
                const ClusterList &clusterList2 = particle2.m_clusterList;

                ClusterSet duplicateSet;

                for (ClusterList::const_iterator cIter1 = clusterList1.begin(), cIterEnd1 = clusterList1.end(); cIter1 != cIterEnd1; ++cIter1)
                {
                    const Cluster *pCluster = *cIter1;

                    if (clusterList2.end() != std::find(clusterList2.begin(), clusterList2.end(), pCluster))
                        duplicateSet.insert(pCluster);
                }

                if (duplicateSet.size() > 0 && clusterList1.size() != clusterList2.size())
                {
                    isUnique = false;
                    break;
                }
            }

            if (!isUnique)
                throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);

            outputParticleList.push_back(particle1);
        }
        catch (StatusCodeException &)
        {
            std::cout << " Warning in CosmicRayTrackRecoveryAlgorithm: found duplicate particles in candidate list " << std::endl;
        }
    }
}

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

void CosmicRayTrackRecoveryAlgorithm::BuildParticles(const ParticleList &particleList)
{
    if (particleList.empty())
        return;

    const PfoList *pPfoList(nullptr); std::string pfoListName;
    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*this, pPfoList, pfoListName));

    for (const Particle &particle : particleList)
    {
        if (!PandoraContentApi::IsAvailable(*this, &particle.m_clusterList))
            continue;

        ClusterList clusterList;
        this->MergeClusters(particle.m_clusterList, clusterList);

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

        PandoraContentApi::ParticleFlowObject::Parameters pfoParameters;
        pfoParameters.m_particleId = MU_MINUS; // TRACK
        pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
        pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
        pfoParameters.m_energy = 0.f;
        pfoParameters.m_momentum = CartesianVector(0.f, 0.f, 0.f);
        pfoParameters.m_clusterList = clusterList;

        const ParticleFlowObject *pPfo(nullptr);
        PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*this, pfoParameters, pPfo));
    }

    if (!pPfoList->empty())
        PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*this, m_outputPfoListName));
}

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

void CosmicRayTrackRecoveryAlgorithm::MergeClusters(const ClusterList &inputClusterList, ClusterList &outputClusterList) const
{
    ClusterList clusterListU, clusterListV, clusterListW;
    LArClusterHelper::GetClustersByHitType(inputClusterList, TPC_VIEW_U, clusterListU);
    LArClusterHelper::GetClustersByHitType(inputClusterList, TPC_VIEW_V, clusterListV);
    LArClusterHelper::GetClustersByHitType(inputClusterList, TPC_VIEW_W, clusterListW);

    for (unsigned int iView = 0; iView < 3; ++iView)
    {
        const ClusterList clusterList((0 == iView) ? clusterListU : (1 == iView) ?  clusterListV : clusterListW);
        const std::string inputClusterListName((0 == iView) ? m_inputClusterListNameU : (1 == iView) ? m_inputClusterListNameV : m_inputClusterListNameW); 

        if (clusterList.empty())
            continue;

        const Cluster *const pSeedCluster(clusterList.front());

        if (!PandoraContentApi::IsAvailable(*this, pSeedCluster))
            throw StatusCodeException(STATUS_CODE_FAILURE);

        for (const Cluster *const pAssociatedCluster : clusterList)
        {
            if (pAssociatedCluster == pSeedCluster)
                continue;

            if (!PandoraContentApi::IsAvailable(*this, pAssociatedCluster))
                continue;

            PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::MergeAndDeleteClusters(*this, pSeedCluster, pAssociatedCluster,
                inputClusterListName, inputClusterListName));
        }

        outputClusterList.push_back(pSeedCluster);
    }
}

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

StatusCode CosmicRayTrackRecoveryAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
{
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
        "ClusterMinLength", m_clusterMinLength));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
        "ClusterMinSpanZ", m_clusterMinSpanZ));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
        "ClusterMinOverlapX", m_clusterMinOverlapX));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
        "ClusterMaxDeltaX", m_clusterMaxDeltaX));

    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));

    return STATUS_CODE_SUCCESS;
}

} // namespace lar_content