LArConnectionPathwayHelper.cc

Go to the documentation of this file.

 
#include "Pandora/Algorithm.h"
#include "Pandora/Pandora.h"
#include "Pandora/PandoraInternal.h"
 
#include "larpandoracontent/LArShowerRefinement/LArProtoShower.h"
 
#include "larpandoracontent/LArHelpers/LArClusterHelper.h"
#include "larpandoracontent/LArHelpers/LArConnectionPathwayHelper.h"
#include "larpandoracontent/LArHelpers/LArGeometryHelper.h"
#include "larpandoracontent/LArHelpers/LArPfoHelper.h"
 
#include <limits>
 
using namespace pandora;
 
namespace lar_content
{
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool LArConnectionPathwayHelper::SortByDistanceToPoint::operator()(const CartesianVector &lhs, const CartesianVector &rhs)
{
    return (m_referencePoint.GetDistanceSquared(lhs) < m_referencePoint.GetDistanceSquared(rhs));
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool LArConnectionPathwayHelper::SortByDistanceToPoint::operator()(const CaloHit *const lhs, const CaloHit *const rhs)
{
    return (m_referencePoint.GetDistanceSquared(lhs->GetPositionVector()) < m_referencePoint.GetDistanceSquared(rhs->GetPositionVector()));
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
// Finding 3D Shower Start Position
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool LArConnectionPathwayHelper::FindShowerStarts3D(const Algorithm *const pAlgorithm, const ParticleFlowObject *const pShowerPfo,
    const ProtoShowerMatch &protoShowerMatch, const CartesianVector &nuVertexPosition, const float maxSeparationFromHit,
    const float maxProjectionSeparation, const float maxXSeparation, CartesianPointVector &showerStarts3D)
{
    const ProtoShower &protoShowerU(protoShowerMatch.GetProtoShowerU());
    const ProtoShower &protoShowerV(protoShowerMatch.GetProtoShowerV());
    const ProtoShower &protoShowerW(protoShowerMatch.GetProtoShowerW());
    const Consistency consistency(protoShowerMatch.GetConsistencyType());
 
    bool uFound(false), vFound(false), wFound(false);
    CartesianVector uShowerStart3D(0.f, 0.f, 0.f), vShowerStart3D(0.f, 0.f, 0.f), wShowerStart3D(0.f, 0.f, 0.f);
 
    CaloHitList caloHitList3D;
    LArPfoHelper::GetCaloHits(pShowerPfo, TPC_3D, caloHitList3D);
 
    if (consistency == Consistency::POSITION)
    {
        LArConnectionPathwayHelper::FindShowerStartFromPosition(pAlgorithm, protoShowerU, protoShowerV, protoShowerW, uShowerStart3D);
        vShowerStart3D = uShowerStart3D;
        wShowerStart3D = uShowerStart3D;
 
        if (LArClusterHelper::GetClosestDistance(uShowerStart3D, caloHitList3D) < maxSeparationFromHit)
        {
            uFound = true;
            vFound = true;
            wFound = true;
        }
    }
    else if (consistency == Consistency::DIRECTION)
    {
        if (LArConnectionPathwayHelper::FindShowerStartFromDirection(
                pAlgorithm, protoShowerU, protoShowerV, protoShowerW, uShowerStart3D, vShowerStart3D, wShowerStart3D))
        {
            if (LArClusterHelper::GetClosestDistance(uShowerStart3D, caloHitList3D) < maxSeparationFromHit)
                uFound = true;
 
            if (LArClusterHelper::GetClosestDistance(vShowerStart3D, caloHitList3D) < maxSeparationFromHit)
                vFound = true;
 
            if (LArClusterHelper::GetClosestDistance(wShowerStart3D, caloHitList3D) < maxSeparationFromHit)
                wFound = true;
        }
    }
 
    if (!uFound || (consistency == Consistency::X_PROJECTION))
    {
        if (LArConnectionPathwayHelper::FindShowerStartFromXProjection(
                pAlgorithm, protoShowerU, protoShowerV, protoShowerW, maxProjectionSeparation, maxXSeparation, uShowerStart3D) ||
            LArConnectionPathwayHelper::FindShowerStartFromXProjectionRelaxed(
                pAlgorithm, protoShowerU, protoShowerV, protoShowerW, maxProjectionSeparation, maxXSeparation, uShowerStart3D) ||
            LArConnectionPathwayHelper::FindShowerStartFromXProjectionRelaxed(
                pAlgorithm, protoShowerU, protoShowerW, protoShowerV, maxProjectionSeparation, maxXSeparation, uShowerStart3D))
        {
            uFound = true;
        }
 
        if (LArConnectionPathwayHelper::FindShowerStartFromXProjection(
                pAlgorithm, protoShowerV, protoShowerW, protoShowerU, maxProjectionSeparation, maxXSeparation, vShowerStart3D) ||
            LArConnectionPathwayHelper::FindShowerStartFromXProjectionRelaxed(
                pAlgorithm, protoShowerV, protoShowerW, protoShowerU, maxProjectionSeparation, maxXSeparation, vShowerStart3D) ||
            LArConnectionPathwayHelper::FindShowerStartFromXProjectionRelaxed(
                pAlgorithm, protoShowerV, protoShowerU, protoShowerW, maxProjectionSeparation, maxXSeparation, vShowerStart3D))
        {
            vFound = true;
        }
 
        if (LArConnectionPathwayHelper::FindShowerStartFromXProjection(
                pAlgorithm, protoShowerW, protoShowerU, protoShowerV, maxProjectionSeparation, maxXSeparation, wShowerStart3D) ||
            LArConnectionPathwayHelper::FindShowerStartFromXProjectionRelaxed(
                pAlgorithm, protoShowerW, protoShowerU, protoShowerV, maxProjectionSeparation, maxXSeparation, wShowerStart3D) ||
            LArConnectionPathwayHelper::FindShowerStartFromXProjectionRelaxed(
                pAlgorithm, protoShowerW, protoShowerV, protoShowerU, maxProjectionSeparation, maxXSeparation, wShowerStart3D))
        {
            wFound = true;
        }
    }
 
    CartesianPointVector foundShowerStarts3D;
 
    if (uFound)
        foundShowerStarts3D.push_back(uShowerStart3D);
 
    if (vFound)
        foundShowerStarts3D.push_back(vShowerStart3D);
 
    if (wFound)
        foundShowerStarts3D.push_back(wShowerStart3D);
 
    if (foundShowerStarts3D.empty())
        return false;
 
    std::sort(foundShowerStarts3D.begin(), foundShowerStarts3D.end(), LArConnectionPathwayHelper::SortByDistanceToPoint(nuVertexPosition));
 
    // ATTN: We need there to be three
    showerStarts3D.push_back(foundShowerStarts3D.front());
    showerStarts3D.push_back((foundShowerStarts3D.size() == 3) ? foundShowerStarts3D[1] : foundShowerStarts3D[0]);
    showerStarts3D.push_back(foundShowerStarts3D.back());
 
    return true;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void LArConnectionPathwayHelper::FindShowerStartFromPosition(const Algorithm *const pAlgorithm, const ProtoShower &protoShowerU,
    const ProtoShower &protoShowerV, const ProtoShower &protoShowerW, CartesianVector &showerStart3D)
{
    const CartesianVector showerStartU(protoShowerU.GetShowerCore().GetStartPosition());
    const CartesianVector showerStartV(protoShowerV.GetShowerCore().GetStartPosition());
    const CartesianVector showerStartW(protoShowerW.GetShowerCore().GetStartPosition());
 
    float chi2(0.0);
    LArGeometryHelper::MergeThreePositions3D(
        pAlgorithm->GetPandora(), TPC_VIEW_U, TPC_VIEW_V, TPC_VIEW_W, showerStartU, showerStartV, showerStartW, showerStart3D, chi2);
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool LArConnectionPathwayHelper::FindShowerStartFromDirection(const Algorithm *const pAlgorithm, const ProtoShower &protoShowerU,
    const ProtoShower &protoShowerV, const ProtoShower &protoShowerW, CartesianVector &uShowerStart3D, CartesianVector &vShowerStart3D,
    CartesianVector &wShowerStart3D)
{
 
    if (!LArConnectionPathwayHelper::FindShowerStartFromDirection(pAlgorithm, protoShowerU, protoShowerV, protoShowerW, uShowerStart3D))
        return false;
 
    if (!LArConnectionPathwayHelper::FindShowerStartFromDirection(pAlgorithm, protoShowerV, protoShowerW, protoShowerU, vShowerStart3D))
        return false;
 
    if (!LArConnectionPathwayHelper::FindShowerStartFromDirection(pAlgorithm, protoShowerW, protoShowerU, protoShowerV, wShowerStart3D))
        return false;
 
    return true;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool LArConnectionPathwayHelper::FindShowerStartFromDirection(const Algorithm *const pAlgorithm, const ProtoShower &protoShower,
    const ProtoShower &protoShowerA, const ProtoShower &protoShowerB, CartesianVector &showerStart3D)
{
    const CartesianVector &showerStart(protoShower.GetShowerCore().GetStartPosition());
    const float x(showerStart.GetX());
    CartesianVector showerStartA(0.f, 0.f, 0.f), showerStartB(0.f, 0.f, 0.f);
 
    if (!LArConnectionPathwayHelper::ProjectShowerStartByDirection(protoShowerA, x, showerStartA))
        return false;
 
    if (!LArConnectionPathwayHelper::ProjectShowerStartByDirection(protoShowerB, x, showerStartB))
        return false;
 
    float chi2(0.0);
    const HitType hitType(protoShower.GetSpineHitList().front()->GetHitType());
    const HitType hitTypeA(protoShowerA.GetSpineHitList().front()->GetHitType());
    const HitType hitTypeB(protoShowerB.GetSpineHitList().front()->GetHitType());
 
    LArGeometryHelper::MergeThreePositions3D(
        pAlgorithm->GetPandora(), hitType, hitTypeA, hitTypeB, showerStart, showerStartA, showerStartB, showerStart3D, chi2);
 
    return false;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool LArConnectionPathwayHelper::ProjectShowerStartByDirection(const ProtoShower &protoShower, const float x, CartesianVector &showerStart2D)
{
    const CartesianVector &viewShowerStart(protoShower.GetShowerCore().GetStartPosition());
    const CartesianVector &viewNuVertex(protoShower.GetConnectionPathway().GetStartPosition());
    const CartesianVector &viewPeakDirection(protoShower.GetConnectionPathway().GetStartDirection());
    const CartesianVector &displacement(viewShowerStart - viewNuVertex);
    const float transverse(viewPeakDirection.GetCrossProduct(displacement).GetMagnitude());
 
    if (transverse > 1.f)
        return false;
 
    const CartesianVector xAxis(1.f, 0.f, 0.f);
    const float cosTheta(std::fabs(viewPeakDirection.GetCosOpeningAngle(xAxis)));
 
    showerStart2D = viewNuVertex + (viewPeakDirection * (std::fabs(x - viewNuVertex.GetX()) * cosTheta));
 
    return true;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool LArConnectionPathwayHelper::FindShowerStartFromXProjection(const Algorithm *const pAlgorithm, const ProtoShower &protoShower,
    const ProtoShower &protoShower1, const ProtoShower &protoShower2, const float maxSeparation, const float maxXSeparation, CartesianVector &showerStart3D)
{
    const CartesianVector showerStart(protoShower.GetShowerCore().GetStartPosition());
    CartesianVector showerStart1(0.f, 0.f, 0.f), showerStart2(0.f, 0.f, 0.f);
 
    const HitType hitType(protoShower.GetSpineHitList().front()->GetHitType());
    const HitType hitType1(protoShower1.GetSpineHitList().front()->GetHitType());
    const HitType hitType2(protoShower2.GetSpineHitList().front()->GetHitType());
 
    if (!LArConnectionPathwayHelper::FindClosestSpinePosition(protoShower1, showerStart, maxXSeparation, showerStart1))
        return false;
 
    if (!LArConnectionPathwayHelper::FindClosestSpinePosition(protoShower2, showerStart, maxXSeparation, showerStart2))
        return false;
 
    float chi2(0.f);
    CartesianVector projection(0.f, 0.f, 0.f), projection1(0.f, 0.f, 0.f), projection2(0.f, 0.f, 0.f);
 
    LArGeometryHelper::MergeTwoPositions(pAlgorithm->GetPandora(), hitType1, hitType2, showerStart1, showerStart2, projection, chi2);
    const float separationSquared((projection - showerStart).GetMagnitudeSquared());
 
    LArGeometryHelper::MergeTwoPositions(pAlgorithm->GetPandora(), hitType, hitType2, showerStart, showerStart2, projection1, chi2);
    const float separationSquared1((projection1 - showerStart1).GetMagnitudeSquared());
 
    LArGeometryHelper::MergeTwoPositions(pAlgorithm->GetPandora(), hitType, hitType1, showerStart, showerStart1, projection2, chi2);
    const float separationSquared2((projection2 - showerStart2).GetMagnitudeSquared());
 
    if ((separationSquared > maxSeparation * maxSeparation) || (separationSquared1 > maxSeparation * maxSeparation) ||
        (separationSquared2 > maxSeparation * maxSeparation))
    {
        return false;
    }
 
    LArGeometryHelper::MergeThreePositions3D(
        pAlgorithm->GetPandora(), hitType, hitType1, hitType2, showerStart, showerStart1, showerStart2, showerStart3D, chi2);
 
    return true;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool LArConnectionPathwayHelper::FindClosestSpinePosition(
    const ProtoShower &protoShower, const CartesianVector &showerStart3D, const float maxXSeparation, CartesianVector &foundShowerStart)
{
    bool found(false);
    float lowestL(std::numeric_limits<float>::max());
 
    for (const CaloHit *const pCaloHit : protoShower.GetSpineHitList())
    {
        if (std::fabs(pCaloHit->GetPositionVector().GetX() - showerStart3D.GetX()) > maxXSeparation)
            continue;
 
        const float lVertex(protoShower.GetConnectionPathway().GetStartDirection().GetDotProduct(
            pCaloHit->GetPositionVector() - protoShower.GetConnectionPathway().GetStartPosition()));
 
        if ((lVertex > 0.f) && (lVertex < lowestL))
        {
            lowestL = lVertex;
            foundShowerStart = pCaloHit->GetPositionVector();
            found = true;
        }
    }
 
    return found;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool LArConnectionPathwayHelper::FindShowerStartFromXProjectionRelaxed(const Algorithm *const pAlgorithm, const ProtoShower &protoShower,
    const ProtoShower &protoShower1, const ProtoShower &protoShower2, const float maxSeparation, const float maxXSeparation, CartesianVector &showerStart3D)
{
    const CartesianVector showerStart(protoShower.GetShowerCore().GetStartPosition());
    CartesianVector showerStart1(0.f, 0.f, 0.f);
 
    const HitType hitType(protoShower.GetSpineHitList().front()->GetHitType());
    const HitType hitType1(protoShower1.GetSpineHitList().front()->GetHitType());
 
    if (!LArConnectionPathwayHelper::FindClosestSpinePosition(protoShower1, showerStart, maxXSeparation, showerStart1))
        return false;
 
    float chi2(0.f);
    CartesianVector projection2(0.f, 0.f, 0.f);
 
    LArGeometryHelper::MergeTwoPositions(pAlgorithm->GetPandora(), hitType, hitType1, showerStart, showerStart1, projection2, chi2);
 
    // Now make sure that they agree...
    const CaloHitList &caloHitList2(protoShower2.GetSpineHitList());
    const float separation(LArClusterHelper::GetClosestDistance(projection2, caloHitList2));
 
    if (separation > maxSeparation)
        return false;
 
    LArGeometryHelper::MergeTwoPositions3D(pAlgorithm->GetPandora(), hitType, hitType1, showerStart, showerStart1, showerStart3D, chi2);
 
    return true;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
// Ordeing Functions
//------------------------------------------------------------------------------------------------------------------------------------------
 
void LArConnectionPathwayHelper::GetMinMiddleMax(
    const float value1, const float value2, const float value3, float &minValue, float &middleValue, float &maxValue)
{
    FloatVector values({value1, value2, value3});
 
    std::sort(values.begin(), values.end());
 
    minValue = values.at(0);
    middleValue = values.at(1);
    maxValue = values.at(2);
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
} // namespace lar_content