pandora_srcnav/_cosmic_ray_splitting_algorithm_8cc_source.html

#include "Pandora/AlgorithmHeaders.h"


#include "larpandoracontent/LArTwoDReco/LArCosmicRay/CosmicRaySplittingAlgorithm.h"


#include "larpandoracontent/LArHelpers/LArClusterHelper.h"

#include "larpandoracontent/LArHelpers/LArGeometryHelper.h"

#include "larpandoracontent/LArHelpers/LArPointingClusterHelper.h"


using namespace pandora;


namespace lar_content

{


CosmicRaySplittingAlgorithm::CosmicRaySplittingAlgorithm() :

    m_clusterMinLength(10.f),

    m_halfWindowLayers(30),

    m_samplingPitch(1.f),

    m_maxCosSplittingAngle(0.9925f),

    m_minCosMergingAngle(0.94f),

    m_maxTransverseDisplacement(5.f),

    m_maxLongitudinalDisplacement(30.f),

    m_maxLongitudinalDisplacementSquared(m_maxLongitudinalDisplacement * m_maxLongitudinalDisplacement)

{

}


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


StatusCode CosmicRaySplittingAlgorithm::Run()

{

    const ClusterList *pClusterList = NULL;

    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pClusterList));


    // Get ordered list of clean clusters

    ClusterVector clusterVector;

    this->GetListOfCleanClusters(pClusterList, clusterVector);


    // Calculate sliding fit results for clean clusters

    TwoDSlidingFitResultMap slidingFitResultMap;

    this->BuildSlidingFitResultMap(clusterVector, slidingFitResultMap);


    // Loop over clusters, identify and perform splits

    ClusterSet splitClusters;


    for (ClusterVector::const_iterator bIter = clusterVector.begin(), bIterEnd1 = clusterVector.end(); bIter != bIterEnd1; ++bIter)

    {

        if (splitClusters.count(*bIter) > 0)

            continue;


        TwoDSlidingFitResultMap::const_iterator bFitIter = slidingFitResultMap.find(*bIter);


        if (slidingFitResultMap.end() == bFitIter)

            continue;


        const TwoDSlidingFitResult &branchSlidingFitResult(bFitIter->second);


        // Find best split position for candidate branch cluster

        CartesianVector splitPosition(0.f, 0.f, 0.f);

        CartesianVector splitDirection1(0.f, 0.f, 0.f);

        CartesianVector splitDirection2(0.f, 0.f, 0.f);


        if (STATUS_CODE_SUCCESS != this->FindBestSplitPosition(branchSlidingFitResult, splitPosition, splitDirection1, splitDirection2))

            continue;


        // Find candidate replacement clusters to merge into branch cluster at the split position

        TwoDSlidingFitResultMap::const_iterator bestReplacementIter1(slidingFitResultMap.end());

        TwoDSlidingFitResultMap::const_iterator bestReplacementIter2(slidingFitResultMap.end());


        float bestLengthSquared1(m_maxLongitudinalDisplacementSquared);

        float bestLengthSquared2(m_maxLongitudinalDisplacementSquared);


        for (ClusterVector::const_iterator rIter = clusterVector.begin(), rIterEnd = clusterVector.end(); rIter != rIterEnd; ++rIter)

        {

            if (splitClusters.count(*rIter) > 0)

                continue;


            TwoDSlidingFitResultMap::const_iterator rFitIter = slidingFitResultMap.find(*rIter);


            if (slidingFitResultMap.end() == rFitIter)

                continue;


            const TwoDSlidingFitResult &replacementSlidingFitResult(rFitIter->second);


            if (branchSlidingFitResult.GetCluster() == replacementSlidingFitResult.GetCluster())

                continue;


            float lengthSquared1(std::numeric_limits<float>::max());

            float lengthSquared2(std::numeric_limits<float>::max());


            if (STATUS_CODE_SUCCESS != this->ConfirmSplitPosition(branchSlidingFitResult, replacementSlidingFitResult, splitPosition,

                                           splitDirection1, splitDirection2, lengthSquared1, lengthSquared2))

                continue;


            if (lengthSquared1 < bestLengthSquared1)

            {

                bestLengthSquared1 = lengthSquared1;

                bestReplacementIter1 = rFitIter;

            }


            if (lengthSquared2 < bestLengthSquared2)

            {

                bestLengthSquared2 = lengthSquared2;

                bestReplacementIter2 = rFitIter;

            }

        }


        const Cluster *pReplacementCluster1 = NULL;

        const Cluster *pReplacementCluster2 = NULL;


        if (slidingFitResultMap.end() != bestReplacementIter1)

            pReplacementCluster1 = bestReplacementIter1->first;


        if (slidingFitResultMap.end() != bestReplacementIter2)

            pReplacementCluster2 = bestReplacementIter2->first;


        if (NULL == pReplacementCluster1 && NULL == pReplacementCluster2)

            continue;


        const Cluster *const pBranchCluster = branchSlidingFitResult.GetCluster();


        // Crossed tracks

        if (pReplacementCluster1 && pReplacementCluster2)

        {

            if (!(this->IdentifyCrossedTracks(pBranchCluster, pReplacementCluster1, pReplacementCluster2, splitPosition)))

                continue;


            PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=,

                this->PerformDoubleSplit(pBranchCluster, pReplacementCluster1, pReplacementCluster2, splitPosition, splitDirection1, splitDirection2));

        }

        // Single scatter

        else if (pReplacementCluster1)

        {

            PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=,

                this->PerformSingleSplit(pBranchCluster, pReplacementCluster1, splitPosition, splitDirection1, splitDirection2));

        }

        else if (pReplacementCluster2)

        {

            PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=,

                this->PerformSingleSplit(pBranchCluster, pReplacementCluster2, splitPosition, splitDirection2, splitDirection1));

        }


        // Choose not to re-use clusters (for now)

        if (pReplacementCluster1)

            splitClusters.insert(pReplacementCluster1);


        if (pReplacementCluster2)

            splitClusters.insert(pReplacementCluster2);


        splitClusters.insert(pBranchCluster);

    }


    return STATUS_CODE_SUCCESS;

}


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


void CosmicRaySplittingAlgorithm::GetListOfCleanClusters(const ClusterList *const pClusterList, ClusterVector &clusterVector) const

{

    for (ClusterList::const_iterator iter = pClusterList->begin(), iterEnd = pClusterList->end(); iter != iterEnd; ++iter)

    {

        const Cluster *const pCluster = *iter;


        if (LArClusterHelper::GetLengthSquared(pCluster) < m_clusterMinLength * m_clusterMinLength)

            continue;


        clusterVector.push_back(pCluster);

    }


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

}


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


void CosmicRaySplittingAlgorithm::BuildSlidingFitResultMap(const ClusterVector &clusterVector, TwoDSlidingFitResultMap &slidingFitResultMap) const

{

    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;

            }

        }

    }

}


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


StatusCode CosmicRaySplittingAlgorithm::FindBestSplitPosition(const TwoDSlidingFitResult &branchSlidingFitResult,

    CartesianVector &splitPosition, CartesianVector &splitDirection1, CartesianVector &splitDirection2) const

{

    // Find position of greatest scatter for this cluster

    float splitCosTheta(m_maxCosSplittingAngle);

    bool foundSplit(false);


    const CartesianVector &minPosition(branchSlidingFitResult.GetGlobalMinLayerPosition());

    const CartesianVector &maxPosition(branchSlidingFitResult.GetGlobalMaxLayerPosition());

    const float halfWindowLength(branchSlidingFitResult.GetLayerFitHalfWindowLength());


    float minL(0.f), maxL(0.f), minT(0.f), maxT(0.f);

    branchSlidingFitResult.GetLocalPosition(minPosition, minL, minT);

    branchSlidingFitResult.GetLocalPosition(maxPosition, maxL, maxT);


    const unsigned int nSamplingPoints = static_cast<unsigned int>((maxL - minL) / m_samplingPitch);


    for (unsigned int n = 0; n < nSamplingPoints; ++n)

    {

        const float alpha((0.5f + static_cast<float>(n)) / static_cast<float>(nSamplingPoints));

        const float rL(minL + (maxL - minL) * alpha);


        CartesianVector centralPosition(0.f, 0.f, 0.f), forwardDirection(0.f, 0.f, 0.f), backwardDirection(0.f, 0.f, 0.f);


        if ((STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitPosition(rL, centralPosition)) ||

            (STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitDirection(rL + halfWindowLength, forwardDirection)) ||

            (STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitDirection(rL - halfWindowLength, backwardDirection)))

        {

            continue;

        }


        const float cosTheta(forwardDirection.GetDotProduct(backwardDirection));


        if (cosTheta < splitCosTheta)

        {

            CartesianVector forwardPosition(0.f, 0.f, 0.f);

            CartesianVector backwardPosition(0.f, 0.f, 0.f);


            if ((STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitPosition(rL + halfWindowLength, forwardPosition)) ||

                (STATUS_CODE_SUCCESS != branchSlidingFitResult.GetGlobalFitPosition(rL - halfWindowLength, backwardPosition)))

            {

                continue;

            }


            CartesianVector forwardVectorOutwards(forwardPosition - centralPosition);

            CartesianVector backwardVectorOutwards(backwardPosition - centralPosition);


            splitPosition = centralPosition;

            splitDirection1 = (forwardDirection.GetDotProduct(forwardVectorOutwards) > 0.f) ? forwardDirection : forwardDirection * -1.f;

            splitDirection2 = (backwardDirection.GetDotProduct(backwardVectorOutwards) > 0.f) ? backwardDirection : backwardDirection * -1.f;

            splitCosTheta = cosTheta;

            foundSplit = true;

        }

    }


    if (!foundSplit)

        return STATUS_CODE_NOT_FOUND;


    return STATUS_CODE_SUCCESS;

}


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


StatusCode CosmicRaySplittingAlgorithm::ConfirmSplitPosition(const TwoDSlidingFitResult &branchSlidingFitResult,

    const TwoDSlidingFitResult &replacementSlidingFitResult, const CartesianVector &splitPosition, const CartesianVector &splitDirection1,

    const CartesianVector &splitDirection2, float &bestLengthSquared1, float &bestLengthSquared2) const

{

    // Check if the replacement cluster points to the split position on the branch cluster

    bestLengthSquared1 = std::numeric_limits<float>::max();

    bestLengthSquared2 = std::numeric_limits<float>::max();


    bool foundMatch(false);


    const CartesianVector minPosition(replacementSlidingFitResult.GetGlobalMinLayerPosition());

    const CartesianVector maxPosition(replacementSlidingFitResult.GetGlobalMaxLayerPosition());

    const CartesianVector minDirection(replacementSlidingFitResult.GetGlobalMinLayerDirection());

    const CartesianVector maxDirection(replacementSlidingFitResult.GetGlobalMaxLayerDirection());


    const CartesianVector branchVertex1(branchSlidingFitResult.GetGlobalMinLayerPosition());

    const CartesianVector branchVertex2(branchSlidingFitResult.GetGlobalMaxLayerPosition());

    const CartesianVector branchDirection1(splitDirection1 * -1.f);

    const CartesianVector branchDirection2(splitDirection2 * -1.f);


    const float cosSplittingAngle(-splitDirection1.GetDotProduct(splitDirection2));

    const float branchLengthSquared((branchVertex2 - branchVertex1).GetMagnitudeSquared());

    const float replacementLengthSquared((maxPosition - minPosition).GetMagnitudeSquared());


    // Loop over each end of the replacement cluster

    for (unsigned int iFwd = 0; iFwd < 2; ++iFwd)

    {

        const CartesianVector pAxis((0 == iFwd) ? (maxPosition - minPosition) : (minPosition - maxPosition));

        const CartesianVector vtxPosition((0 == iFwd) ? minPosition : maxPosition);

        const CartesianVector endPosition((0 == iFwd) ? maxPosition : minPosition);

        const CartesianVector vtxDirection((0 == iFwd) ? (pAxis.GetDotProduct(minDirection) > 0 ? minDirection : minDirection * -1.f)

                                                       : (pAxis.GetDotProduct(maxDirection) > 0 ? maxDirection : maxDirection * -1.f));


        // Choose the correct end of the replacement cluster and require proximity to the branch cluster

        const float vtxDisplacement(LArClusterHelper::GetClosestDistance(vtxPosition, branchSlidingFitResult.GetCluster()));

        const float endDisplacement(LArClusterHelper::GetClosestDistance(endPosition, branchSlidingFitResult.GetCluster()));


        const float lengthSquared((vtxPosition - splitPosition).GetMagnitudeSquared());

        const float lengthSquared1((vtxPosition - branchVertex1).GetMagnitudeSquared());

        const float lengthSquared2((vtxPosition - branchVertex2).GetMagnitudeSquared());


        if (vtxDisplacement > endDisplacement)

            continue;


        if (std::min(lengthSquared, std::min(lengthSquared1, lengthSquared2)) > std::min(branchLengthSquared, replacementLengthSquared))

            continue;


        // Require good pointing information between replacement cluster and candidate split position

        float impactL(0.f), impactT(0.f);

        LArPointingClusterHelper::GetImpactParameters(vtxPosition, vtxDirection, splitPosition, impactL, impactT);


        if (impactT > m_maxTransverseDisplacement || impactL > m_maxLongitudinalDisplacement || impactL < -1.f ||

            impactT > std::max(1.5f, 0.577f * (1.f + impactL)))

            continue;


        // Check the segment of the branch cluster above the split position

        if (vtxDirection.GetDotProduct(branchDirection1) > std::max(m_minCosMergingAngle, cosSplittingAngle))

        {

            if (lengthSquared < bestLengthSquared1)

            {

                bestLengthSquared1 = lengthSquared;

                foundMatch = true;

            }

        }


        // Check the segment of the branch cluster below the split position

        if (vtxDirection.GetDotProduct(branchDirection2) > std::max(m_minCosMergingAngle, cosSplittingAngle))

        {

            if (lengthSquared < bestLengthSquared2)

            {

                bestLengthSquared2 = lengthSquared;

                foundMatch = true;

            }

        }

    }


    if (!foundMatch)

        return STATUS_CODE_NOT_FOUND;


    return STATUS_CODE_SUCCESS;

}


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


StatusCode CosmicRaySplittingAlgorithm::PerformSingleSplit(const Cluster *const pBranchCluster, const Cluster *const pReplacementCluster,

    const CartesianVector &splitPosition, const CartesianVector &forwardDirection, const CartesianVector &backwardDirection) const

{

    CaloHitList caloHitListToMove;

    this->GetCaloHitListToMove(pBranchCluster, pReplacementCluster, splitPosition, forwardDirection, backwardDirection, caloHitListToMove);


    CaloHitList caloHitListToKeep;

    this->GetCaloHitListToKeep(pBranchCluster, caloHitListToMove, caloHitListToKeep);


    if (caloHitListToKeep.empty())

        return PandoraContentApi::MergeAndDeleteClusters(*this, pReplacementCluster, pBranchCluster);


    return this->SplitCluster(pBranchCluster, pReplacementCluster, caloHitListToMove);

}


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


StatusCode CosmicRaySplittingAlgorithm::PerformDoubleSplit(const Cluster *const pBranchCluster, const Cluster *const pReplacementCluster1,

    const Cluster *const pReplacementCluster2, const CartesianVector &splitPosition, const CartesianVector &splitDirection1,

    const CartesianVector &splitDirection2) const

{

    CaloHitList caloHitListToMove1, caloHitListToMove2;

    this->GetCaloHitListsToMove(pBranchCluster, splitPosition, splitDirection1, splitDirection2, caloHitListToMove1, caloHitListToMove2);


    CaloHitList caloHitListToKeep1;

    this->GetCaloHitListToKeep(pBranchCluster, caloHitListToMove1, caloHitListToKeep1);


    if (caloHitListToKeep1.empty())

        return STATUS_CODE_FAILURE;


    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, this->SplitCluster(pBranchCluster, pReplacementCluster1, caloHitListToMove1));


    CaloHitList caloHitListToKeep2;

    this->GetCaloHitListToKeep(pBranchCluster, caloHitListToMove2, caloHitListToKeep2);


    if (caloHitListToKeep2.empty())

        return PandoraContentApi::MergeAndDeleteClusters(*this, pReplacementCluster2, pBranchCluster);


    return this->SplitCluster(pBranchCluster, pReplacementCluster2, caloHitListToMove2);

}


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


void CosmicRaySplittingAlgorithm::GetCaloHitListToMove(const Cluster *const pBranchCluster, const Cluster *const pReplacementCluster,

    const CartesianVector &splitPosition, const CartesianVector &forwardDirection, const CartesianVector &backwardDirection,

    CaloHitList &caloHitListToMove) const

{

    const CartesianVector forwardPosition(LArClusterHelper::GetClosestPosition(splitPosition, pBranchCluster));

    const CartesianVector vtxPosition(LArClusterHelper::GetClosestPosition(splitPosition, pReplacementCluster));

    const CartesianVector vtxDirection((forwardPosition - vtxPosition).GetUnitVector());


    CaloHitList caloHitListToSort, caloHitListToCheck;

    pBranchCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitListToSort);


    for (CaloHitList::const_iterator iter = caloHitListToSort.begin(), iterEnd = caloHitListToSort.end(); iter != iterEnd; ++iter)

    {

        const CaloHit *const pCaloHit = *iter;


        if (forwardDirection.GetDotProduct(pCaloHit->GetPositionVector() - forwardPosition) > 0.f)

        {

            caloHitListToMove.push_back(pCaloHit);

        }

        else if (forwardDirection.GetDotProduct(pCaloHit->GetPositionVector() - vtxPosition) > -1.25f)

        {

            caloHitListToCheck.push_back(pCaloHit);

        }

    }


    float closestLengthSquared(std::numeric_limits<float>::max());


    for (CaloHitList::const_iterator iter = caloHitListToCheck.begin(), iterEnd = caloHitListToCheck.end(); iter != iterEnd; ++iter)

    {

        const CaloHit *const pCaloHit = *iter;


        if (vtxDirection.GetCrossProduct(pCaloHit->GetPositionVector() - forwardPosition).GetMagnitude() >

            backwardDirection.GetCrossProduct(pCaloHit->GetPositionVector() - forwardPosition).GetMagnitude())

            continue;


        const float lengthSquared((pCaloHit->GetPositionVector() - vtxPosition).GetMagnitudeSquared());


        if (lengthSquared < closestLengthSquared)

            closestLengthSquared = lengthSquared;

    }


    for (CaloHitList::const_iterator iter = caloHitListToCheck.begin(), iterEnd = caloHitListToCheck.end(); iter != iterEnd; ++iter)

    {

        const CaloHit *const pCaloHit = *iter;


        if ((pCaloHit->GetPositionVector() - vtxPosition).GetMagnitudeSquared() >= closestLengthSquared)

            caloHitListToMove.push_back(pCaloHit);

    }

}


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


void CosmicRaySplittingAlgorithm::GetCaloHitListsToMove(const Cluster *const pBranchCluster, const CartesianVector &splitPosition,

    const CartesianVector &splitDirection1, const CartesianVector &splitDirection2, CaloHitList &caloHitListToMove1, CaloHitList &caloHitListToMove2) const

{

    CaloHitList caloHitListToSort;

    pBranchCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitListToSort);


    for (CaloHitList::const_iterator iter = caloHitListToSort.begin(), iterEnd = caloHitListToSort.end(); iter != iterEnd; ++iter)

    {

        const CaloHit *const pCaloHit = *iter;


        const float displacement1(splitDirection1.GetDotProduct(pCaloHit->GetPositionVector() - splitPosition));

        const float displacement2(splitDirection2.GetDotProduct(pCaloHit->GetPositionVector() - splitPosition));


        if (displacement1 > displacement2)

        {

            caloHitListToMove1.push_back(pCaloHit);

        }

        else

        {

            caloHitListToMove2.push_back(pCaloHit);

        }

    }

}


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


bool CosmicRaySplittingAlgorithm::IdentifyCrossedTracks(const Cluster *const pBranchCluster, const Cluster *const pReplacementCluster1,

    const Cluster *const pReplacementCluster2, const pandora::CartesianVector &splitPosition) const

{

    CartesianVector branchVertex1(0.f, 0.f, 0.f), branchVertex2(0.f, 0.f, 0.f);

    LArClusterHelper::GetExtremalCoordinates(pBranchCluster, branchVertex1, branchVertex2);


    const CartesianVector replacementVertex1(LArClusterHelper::GetClosestPosition(splitPosition, pReplacementCluster1));

    const CartesianVector replacementVertex2(LArClusterHelper::GetClosestPosition(splitPosition, pReplacementCluster2));


    if ((replacementVertex2 - replacementVertex1).GetMagnitudeSquared() > (branchVertex2 - branchVertex1).GetMagnitudeSquared())

        return false;


    return true;

}


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


StatusCode CosmicRaySplittingAlgorithm::GetCaloHitListToKeep(

    const Cluster *const pBranchCluster, const CaloHitList &caloHitListToMove, CaloHitList &caloHitListToKeep) const

{

    if (caloHitListToMove.empty())

        return STATUS_CODE_FAILURE;


    CaloHitList caloHitList;

    pBranchCluster->GetOrderedCaloHitList().FillCaloHitList(caloHitList);


    for (CaloHitList::const_iterator iter = caloHitList.begin(), iterEnd = caloHitList.end(); iter != iterEnd; ++iter)

    {

        const CaloHit *const pCaloHit = *iter;


        if (caloHitListToMove.end() == std::find(caloHitListToMove.begin(), caloHitListToMove.end(), pCaloHit))

            caloHitListToKeep.push_back(pCaloHit);

    }


    return STATUS_CODE_SUCCESS;

}


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


StatusCode CosmicRaySplittingAlgorithm::SplitCluster(

    const Cluster *const pBranchCluster, const Cluster *const pReplacementCluster, const CaloHitList &caloHitListToMove) const

{

    if (caloHitListToMove.empty())

        return STATUS_CODE_FAILURE;


    for (CaloHitList::const_iterator iter = caloHitListToMove.begin(), iterEnd = caloHitListToMove.end(); iter != iterEnd; ++iter)

    {

        const CaloHit *const pCaloHit = *iter;

        PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::RemoveFromCluster(*this, pBranchCluster, pCaloHit));

        PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToCluster(*this, pReplacementCluster, pCaloHit));

    }


    return STATUS_CODE_SUCCESS;

}


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


StatusCode CosmicRaySplittingAlgorithm::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, "SlidingFitHalfWindow", m_halfWindowLayers));


    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SamplingPitch", m_samplingPitch));


    if (m_samplingPitch < std::numeric_limits<float>::epsilon())

        return STATUS_CODE_INVALID_PARAMETER;


    PANDORA_RETURN_RESULT_IF_AND_IF(

        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxCosSplittingAngle", m_maxCosSplittingAngle));


    PANDORA_RETURN_RESULT_IF_AND_IF(

        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinCosMergingAngle", m_minCosMergingAngle));


    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,

        XmlHelper::ReadValue(xmlHandle, "MaxTransverseDisplacement", m_maxTransverseDisplacement));


    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,

        XmlHelper::ReadValue(xmlHandle, "MaxLongitudinalDisplacement", m_maxLongitudinalDisplacement));

    m_maxLongitudinalDisplacementSquared = m_maxLongitudinalDisplacement * m_maxLongitudinalDisplacement;


    return STATUS_CODE_SUCCESS;

}


} // namespace lar_content

AlgorithmHeaders.h
Grouping of header files for many classes of use in particle flow algorithms.

CosmicRaySplittingAlgorithm.h
Header file for the cosmic ray splitting algorithm class.

LArClusterHelper.h
Header file for the cluster helper class.

LArGeometryHelper.h
Header file for the geometry helper class.

LArPointingClusterHelper.h

PANDORA_RETURN_RESULT_IF_AND_IF
#define PANDORA_RETURN_RESULT_IF_AND_IF(StatusCode1, StatusCode2, Operator, Command)
Definition StatusCodes.h:31

PANDORA_RETURN_RESULT_IF
#define PANDORA_RETURN_RESULT_IF(StatusCode1, Operator, Command)
Definition StatusCodes.h:19

PandoraContentApi::GetCurrentList
static pandora::StatusCode GetCurrentList(const pandora::Algorithm &algorithm, const T *&pT)
Get the current list.
Definition PandoraContentApi.cc:75

PandoraContentApi::MergeAndDeleteClusters
static pandora::StatusCode MergeAndDeleteClusters(const pandora::Algorithm &algorithm, const pandora::Cluster *const pClusterToEnlarge, const pandora::Cluster *const pClusterToDelete)
Merge two clusters in the current list, enlarging one cluster and deleting the second.
Definition PandoraContentApi.cc:301

PandoraContentApi::AddToCluster
static pandora::StatusCode AddToCluster(const pandora::Algorithm &algorithm, const pandora::Cluster *const pCluster, const T *const pT)
Add a calo hit, or a list of calo hits, to a cluster.
Definition PandoraContentApi.cc:215

PandoraContentApi::RemoveFromCluster
static pandora::StatusCode RemoveFromCluster(const pandora::Algorithm &algorithm, const pandora::Cluster *const pCluster, const pandora::CaloHit *const pCaloHit)
Remove a calo hit from a cluster. Note this function will not remove the final calo hit from a cluste...
Definition PandoraContentApi.cc:222

lar_content::CosmicRaySplittingAlgorithm::PerformDoubleSplit
pandora::StatusCode PerformDoubleSplit(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster1, const pandora::Cluster *const pReplacementCluster2, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection1, const pandora::CartesianVector &splitDirection2) const
Split a branch cluster for case of two replacement clusters.
Definition CosmicRaySplittingAlgorithm.cc:371

lar_content::CosmicRaySplittingAlgorithm::FindBestSplitPosition
pandora::StatusCode FindBestSplitPosition(const TwoDSlidingFitResult &slidingFitResult, pandora::CartesianVector &splitPosition, pandora::CartesianVector &splitDirection1, pandora::CartesianVector &splitDirection2) const
Find the position of greatest scatter along a sliding linear fit.
Definition CosmicRaySplittingAlgorithm.cc:207

lar_content::CosmicRaySplittingAlgorithm::m_minCosMergingAngle
float m_minCosMergingAngle
largest relative angle allowed when merging clusters
Definition CosmicRaySplittingAlgorithm.h:165

lar_content::CosmicRaySplittingAlgorithm::ConfirmSplitPosition
pandora::StatusCode ConfirmSplitPosition(const TwoDSlidingFitResult &branchSlidingFitResult, const TwoDSlidingFitResult &replacementSlidingFitResult, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection1, const pandora::CartesianVector &splitDirection2, float &lengthSquared1, float &lengthSquared2) const
Find a second replacement cluster that aligns with the scatter of the first branch cluster.
Definition CosmicRaySplittingAlgorithm.cc:270

lar_content::CosmicRaySplittingAlgorithm::m_maxCosSplittingAngle
float m_maxCosSplittingAngle
smallest scatter angle allowed when splitting cluster
Definition CosmicRaySplittingAlgorithm.h:164

lar_content::CosmicRaySplittingAlgorithm::m_maxLongitudinalDisplacement
float m_maxLongitudinalDisplacement
maximum longitudinal displacement of associated clusters
Definition CosmicRaySplittingAlgorithm.h:167

lar_content::CosmicRaySplittingAlgorithm::ReadSettings
pandora::StatusCode ReadSettings(const pandora::TiXmlHandle xmlHandle)
Read the algorithm settings.
Definition CosmicRaySplittingAlgorithm.cc:531

lar_content::CosmicRaySplittingAlgorithm::m_maxLongitudinalDisplacementSquared
float m_maxLongitudinalDisplacementSquared
longitudinal displacement squared
Definition CosmicRaySplittingAlgorithm.h:168

lar_content::CosmicRaySplittingAlgorithm::GetListOfCleanClusters
void GetListOfCleanClusters(const pandora::ClusterList *const pClusterList, pandora::ClusterVector &clusterVector) const
Populate cluster vector with subset of cluster list, containing clusters judged to be clean.
Definition CosmicRaySplittingAlgorithm.cc:164

lar_content::CosmicRaySplittingAlgorithm::m_maxTransverseDisplacement
float m_maxTransverseDisplacement
maximum transverse displacement of associated clusters
Definition CosmicRaySplittingAlgorithm.h:166

lar_content::CosmicRaySplittingAlgorithm::Run
pandora::StatusCode Run()
Run the algorithm.
Definition CosmicRaySplittingAlgorithm.cc:36

lar_content::CosmicRaySplittingAlgorithm::m_halfWindowLayers
unsigned int m_halfWindowLayers
number of layers to use for half-window of sliding fit
Definition CosmicRaySplittingAlgorithm.h:162

lar_content::CosmicRaySplittingAlgorithm::BuildSlidingFitResultMap
void BuildSlidingFitResultMap(const pandora::ClusterVector &clusterVector, TwoDSlidingFitResultMap &slidingFitResultMap) const
Build the map of sliding fit results.
Definition CosmicRaySplittingAlgorithm.cc:181

lar_content::CosmicRaySplittingAlgorithm::GetCaloHitListToKeep
pandora::StatusCode GetCaloHitListToKeep(const pandora::Cluster *const pBranchCluster, const pandora::CaloHitList &caloHitListToMove, pandora::CaloHitList &caloHitListToKeep) const
Split the branch cluster and add hits to the replacement cluster.
Definition CosmicRaySplittingAlgorithm.cc:491

lar_content::CosmicRaySplittingAlgorithm::GetCaloHitListToMove
void GetCaloHitListToMove(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &forwardDirection, const pandora::CartesianVector &backwardDirection, pandora::CaloHitList &caloHitList) const
Get list of calo hits to move in order to split a branch cluster into two segments for case of one re...
Definition CosmicRaySplittingAlgorithm.cc:397

lar_content::CosmicRaySplittingAlgorithm::IdentifyCrossedTracks
bool IdentifyCrossedTracks(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster1, const pandora::Cluster *const pReplacementCluster2, const pandora::CartesianVector &splitPosition) const
Identify crossed tracks formed from the branch cluster and its replacement cluster.
Definition CosmicRaySplittingAlgorithm.cc:474

lar_content::CosmicRaySplittingAlgorithm::PerformSingleSplit
pandora::StatusCode PerformSingleSplit(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &forwardDirection, const pandora::CartesianVector &backwardDirection) const
Split a branch cluster for case of one replacement cluster.
Definition CosmicRaySplittingAlgorithm.cc:354

lar_content::CosmicRaySplittingAlgorithm::SplitCluster
pandora::StatusCode SplitCluster(const pandora::Cluster *const pBranchCluster, const pandora::Cluster *const pReplacementCluster, const pandora::CaloHitList &caloHitListToMove) const
Split the branch cluster and add hits to the replacement cluster.
Definition CosmicRaySplittingAlgorithm.cc:513

lar_content::CosmicRaySplittingAlgorithm::m_samplingPitch
float m_samplingPitch
sampling pitch for walking along sliding linear fit
Definition CosmicRaySplittingAlgorithm.h:163

lar_content::CosmicRaySplittingAlgorithm::CosmicRaySplittingAlgorithm
CosmicRaySplittingAlgorithm()
Default constructor.
Definition CosmicRaySplittingAlgorithm.cc:22

lar_content::CosmicRaySplittingAlgorithm::m_clusterMinLength
float m_clusterMinLength
minimum length of clusters for this algorithm
Definition CosmicRaySplittingAlgorithm.h:161

lar_content::CosmicRaySplittingAlgorithm::GetCaloHitListsToMove
void GetCaloHitListsToMove(const pandora::Cluster *const pBranchCluster, const pandora::CartesianVector &splitPosition, const pandora::CartesianVector &splitDirection1, const pandora::CartesianVector &splitDirection2, pandora::CaloHitList &caloHitList1, pandora::CaloHitList &caloHitList2) const
Get lists of calo hits to move in order to split a branch cluster into two segments for case of two r...
Definition CosmicRaySplittingAlgorithm.cc:449

lar_content::LArClusterHelper::GetClosestPosition
static pandora::CartesianVector GetClosestPosition(const pandora::CartesianVector &position, const pandora::ClusterList &clusterList)
Get closest position in a list of clusters to a specified input position vector.
Definition LArClusterHelper.cc:221

lar_content::LArClusterHelper::SortByNHits
static bool SortByNHits(const pandora::Cluster *const pLhs, const pandora::Cluster *const pRhs)
Sort clusters by number of hits, then layer span, then inner layer, then position,...
Definition LArClusterHelper.cc:656

lar_content::LArClusterHelper::GetExtremalCoordinates
static void GetExtremalCoordinates(const pandora::ClusterList &clusterList, pandora::CartesianVector &innerCoordinate, pandora::CartesianVector &outerCoordinate)
Get positions of the two most distant calo hits in a list of cluster (ordered by Z)
Definition LArClusterHelper.cc:424

lar_content::LArClusterHelper::GetLengthSquared
static float GetLengthSquared(const pandora::Cluster *const pCluster)
Get length squared of cluster.
Definition LArClusterHelper.cc:65

lar_content::LArClusterHelper::GetClosestDistance
static float GetClosestDistance(const pandora::ClusterList &clusterList1, const pandora::ClusterList &clusterList2)
Get closest distance between clusters in a pair of cluster lists.
Definition LArClusterHelper.cc:146

lar_content::LArGeometryHelper::GetWireZPitch
static float GetWireZPitch(const pandora::Pandora &pandora, const float maxWirePitchDiscrepancy=0.01)
Return the wire pitch.
Definition LArGeometryHelper.h:273

lar_content::LArPointingClusterHelper::GetImpactParameters
static void GetImpactParameters(const LArPointingCluster::Vertex &pointingVertex, const LArPointingCluster::Vertex &targetVertex, float &longitudinal, float &transverse)
Calculate impact parameters between a pair of pointing vertices.
Definition LArPointingClusterHelper.cc:199

lar_content::TwoDSlidingFitResult
TwoDSlidingFitResult class.
Definition LArTwoDSlidingFitResult.h:24

lar_content::TwoDSlidingFitResult::GetGlobalMinLayerPosition
pandora::CartesianVector GetGlobalMinLayerPosition() const
Get global position corresponding to the fit result in minimum fit layer.
Definition LArTwoDSlidingFitResult.cc:239

lar_content::TwoDSlidingFitResult::GetLocalPosition
void GetLocalPosition(const pandora::CartesianVector &position, float &rL, float &rT) const
Get local sliding fit coordinates for a given global position.
Definition LArTwoDSlidingFitResult.cc:197

lar_content::TwoDSlidingFitResult::GetGlobalFitPosition
pandora::StatusCode GetGlobalFitPosition(const float rL, pandora::CartesianVector &position) const
Get global fit position for a given longitudinal coordinate.
Definition LArTwoDSlidingFitResult.cc:347

lar_content::TwoDSlidingFitResult::GetGlobalMinLayerDirection
pandora::CartesianVector GetGlobalMinLayerDirection() const
Get global direction corresponding to the fit result in minimum fit layer.
Definition LArTwoDSlidingFitResult.cc:265

lar_content::TwoDSlidingFitResult::GetGlobalMaxLayerPosition
pandora::CartesianVector GetGlobalMaxLayerPosition() const
Get global position corresponding to the fit result in maximum fit layer.
Definition LArTwoDSlidingFitResult.cc:252

lar_content::TwoDSlidingFitResult::GetLayerFitHalfWindowLength
float GetLayerFitHalfWindowLength() const
Get the layer fit half window length.
Definition LArTwoDSlidingFitResult.cc:153

lar_content::TwoDSlidingFitResult::GetGlobalFitDirection
pandora::StatusCode GetGlobalFitDirection(const float rL, pandora::CartesianVector &direction) const
Get global fit direction for a given longitudinal coordinate.
Definition LArTwoDSlidingFitResult.cc:361

lar_content::TwoDSlidingFitResult::GetCluster
const pandora::Cluster * GetCluster() const
Get the address of the cluster, if originally provided.
Definition LArTwoDSlidingFitResult.cc:143

lar_content::TwoDSlidingFitResult::GetGlobalMaxLayerDirection
pandora::CartesianVector GetGlobalMaxLayerDirection() const
Get global direction corresponding to the fit result in maximum fit layer.
Definition LArTwoDSlidingFitResult.cc:278

pandora::CaloHit
CaloHit class.
Definition CaloHit.h:26

pandora::CaloHit::GetPositionVector
const CartesianVector & GetPositionVector() const
Get the position vector of center of calorimeter cell, units mm.
Definition CaloHit.h:350

pandora::CartesianVector
CartesianVector class.
Definition CartesianVector.h:24

pandora::CartesianVector::GetDotProduct
float GetDotProduct(const CartesianVector &rhs) const
Get the dot product of the cartesian vector with a second cartesian vector.
Definition CartesianVector.h:294

pandora::CartesianVector::GetMagnitude
float GetMagnitude() const
Get the magnitude.
Definition CartesianVector.h:280

pandora::CartesianVector::GetCrossProduct
CartesianVector GetCrossProduct(const CartesianVector &rhs) const
Get the cross product of the cartesian vector with a second cartesian vector.
Definition CartesianVector.h:301

pandora::Cluster
Cluster class.
Definition Cluster.h:31

pandora::Cluster::GetOrderedCaloHitList
const OrderedCaloHitList & GetOrderedCaloHitList() const
Get the ordered calo hit list.
Definition Cluster.h:470

pandora::OrderedCaloHitList::FillCaloHitList
void FillCaloHitList(CaloHitList &caloHitList) const
Fill a provided calo hit list with all the calo hits in the ordered calo hit list.
Definition OrderedCaloHitList.cc:131

pandora::Process::GetPandora
const Pandora & GetPandora() const
Get the associated pandora instance.
Definition Process.h:116

pandora::StatusCodeException
StatusCodeException class.
Definition StatusCodes.h:119

pandora::StatusCodeException::GetStatusCode
StatusCode GetStatusCode() const
Get status code.
Definition StatusCodes.h:191

pandora::TiXmlHandle
Definition tinyxml.h:1646

pandora::XmlHelper::ReadValue
static StatusCode ReadValue(const TiXmlHandle &xmlHandle, const std::string &xmlElementName, T &t)
Read a value from an xml element.
Definition XmlHelper.h:136

lar_content
Definition CheatingBeamParticleIdTool.cc:19

lar_content::TwoDSlidingFitResultMap
std::unordered_map< const pandora::Cluster *, TwoDSlidingFitResult > TwoDSlidingFitResultMap
Definition LArTwoDSlidingFitResult.h:551

pandora
Definition LArContent.h:12

pandora::ClusterVector
std::vector< const Cluster * > ClusterVector
Definition PandoraInternal.h:445

pandora::ClusterList
MANAGED_CONTAINER< const Cluster * > ClusterList
Definition PandoraInternal.h:434

pandora::ClusterSet
std::unordered_set< const Cluster * > ClusterSet
Definition PandoraInternal.h:456

pandora::CaloHitList
MANAGED_CONTAINER< const CaloHit * > CaloHitList
Definition PandoraInternal.h:433

pandora::StatusCode
StatusCode
The StatusCode enum.
Definition StatusCodes.h:100