BranchSplittingAlgorithm.cc

Go to the documentation of this file.

#include "Pandora/AlgorithmHeaders.h"

#include "larpandoracontent/LArHelpers/LArPointingClusterHelper.h"

#include "larpandoracontent/LArTwoDReco/LArClusterSplitting/BranchSplittingAlgorithm.h"

using namespace pandora;

namespace lar_content
{

BranchSplittingAlgorithm::BranchSplittingAlgorithm() :
    m_maxTransverseDisplacement(1.5f),
    m_maxLongitudinalDisplacement(10.f),
    m_minLongitudinalExtension(3.f),
    m_minCosRelativeAngle(0.966f),
    m_projectionAngularAllowance(20.f)
{
}

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

void BranchSplittingAlgorithm::FindBestSplitPosition(const TwoDSlidingFitResult &branchSlidingFit, const TwoDSlidingFitResult &principalSlidingFit, 
    CartesianVector &principalStartPosition, CartesianVector &branchSplitPosition, CartesianVector &branchSplitDirection) const
{
    // Conventions:
    // (1) Delta ray is split from the branch cluster
    // (2) Delta ray occurs where the vertex of the principal cluster meets the vertex of the branch cluster
    //     Method loops over the inner and outer positions of the principal and branch clusters, trying all
    //     possible assignments of vertex and end position until a split is found
    for (unsigned int principalForward = 0; principalForward < 2; ++principalForward)
    {
        const CartesianVector principalVertexPosition(1==principalForward ? principalSlidingFit.GetGlobalMinLayerPosition() : principalSlidingFit.GetGlobalMaxLayerPosition());
        const CartesianVector principalEndPosition(1!=principalForward ? principalSlidingFit.GetGlobalMinLayerPosition() : principalSlidingFit.GetGlobalMaxLayerPosition());
        const CartesianVector principalVertexDirection(1==principalForward ? principalSlidingFit.GetGlobalMinLayerDirection() : principalSlidingFit.GetGlobalMaxLayerDirection() * -1.f);

        CartesianVector projectedBranchPosition(0.f,0.f,0.f);
        bool projectedPositionFound(false), projectedPositionFail(false);

        for (unsigned int branchForward = 0; branchForward < 2; ++branchForward)
        {
            const CartesianVector branchVertexPosition(1==branchForward ? branchSlidingFit.GetGlobalMinLayerPosition() : branchSlidingFit.GetGlobalMaxLayerPosition());
            const CartesianVector branchEndPosition(1!=branchForward ? branchSlidingFit.GetGlobalMinLayerPosition() : branchSlidingFit.GetGlobalMaxLayerPosition());
            const CartesianVector branchEndDirection(1!=branchForward ? branchSlidingFit.GetGlobalMinLayerDirection() : branchSlidingFit.GetGlobalMaxLayerDirection() * -1.f);

            if (principalVertexDirection.GetDotProduct(branchEndDirection) < 0.5f)
                continue;

            if ((principalEndPosition - branchEndPosition).GetMagnitudeSquared() < (principalVertexPosition - branchVertexPosition).GetMagnitudeSquared())
                continue;

            // Project the principal vertex onto the branch cluster
            try
            {
                if (!projectedPositionFound && !projectedPositionFail)
                {
                    projectedBranchPosition = LArPointingClusterHelper::GetProjectedPosition(principalVertexPosition, principalVertexDirection, branchSlidingFit.GetCluster(), m_projectionAngularAllowance);
                    projectedPositionFound = true;
                }
            }
            catch (StatusCodeException &)
            {
                projectedPositionFail = true;
            }

            if (!projectedPositionFound || projectedPositionFail)
                continue;

            const float projectedDistanceSquared((projectedBranchPosition - principalVertexPosition).GetMagnitudeSquared());

            if (projectedDistanceSquared > m_maxLongitudinalDisplacement * m_maxLongitudinalDisplacement)
                continue;

            const float commonDistanceSquared((projectedBranchPosition - branchEndPosition).GetMagnitudeSquared());

            if (projectedDistanceSquared > commonDistanceSquared)
                continue;

            const float replacementDistanceSquared((projectedBranchPosition - principalEndPosition).GetMagnitudeSquared());

            if (replacementDistanceSquared < m_minLongitudinalExtension * m_minLongitudinalExtension)
                continue;

            const float branchDistanceSquared((projectedBranchPosition - branchVertexPosition).GetMagnitudeSquared());

            if (branchDistanceSquared > 4.f * replacementDistanceSquared)
                continue;

            // Require that principal vertex and branch projection have good (and improved) pointing
            bool foundSplit(false);

            const float halfWindowLength(branchSlidingFit.GetLayerFitHalfWindowLength());
            const float deltaL(1==branchForward ? +halfWindowLength : -halfWindowLength);

            float localL(0.f), localT(0.f);
            CartesianVector forwardDirection(0.f,0.f,0.f);
            branchSlidingFit.GetLocalPosition(projectedBranchPosition, localL, localT);

            if (STATUS_CODE_SUCCESS != branchSlidingFit.GetGlobalFitDirection(localL + deltaL, forwardDirection))
                continue;

            CartesianVector projectedBranchDirection(1==branchForward ? forwardDirection : forwardDirection * -1.f);
            const float cosTheta(-projectedBranchDirection.GetDotProduct(principalVertexDirection));

            try
            {
                const float currentCosTheta(branchSlidingFit.GetCosScatteringAngle(localL)); 

                if (cosTheta < currentCosTheta)
                    continue;
            }
            catch (StatusCodeException &)
            {
            }

            float rT1(0.f), rL1(0.f), rT2(0.f), rL2(0.f);
            LArPointingClusterHelper::GetImpactParameters(projectedBranchPosition, projectedBranchDirection, principalVertexPosition, rL1, rT1);
            LArPointingClusterHelper::GetImpactParameters(principalVertexPosition, principalVertexDirection, projectedBranchPosition, rL2, rT2);

            if ((cosTheta > m_minCosRelativeAngle) && (rT1 < m_maxTransverseDisplacement) && (rT2 < m_maxTransverseDisplacement))
            {
                foundSplit = true;
                principalStartPosition = principalVertexPosition;
                branchSplitPosition = projectedBranchPosition;
                branchSplitDirection = projectedBranchDirection * -1.f;
            }

            if (foundSplit)
                return;
        }
    }

    throw StatusCodeException(STATUS_CODE_NOT_FOUND);
}

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

StatusCode BranchSplittingAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
{
    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));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
        "MinLongitudinalExtension", m_minLongitudinalExtension));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
        "MinCosRelativeAngle", m_minCosRelativeAngle));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
        "ProjectionAngularAllowance", m_projectionAngularAllowance));

    return TwoDSlidingFitSplittingAndSplicingAlgorithm::ReadSettings(xmlHandle);
}

} // namespace lar_content