ShowerSpineFinderTool.cc

Go to the documentation of this file.

 
#include "Pandora/AlgorithmHeaders.h"
#include "Pandora/AlgorithmTool.h"
 
#include "larpandoracontent/LArHelpers/LArConnectionPathwayHelper.h"
#include "larpandoracontent/LArHelpers/LArGeometryHelper.h"
#include "larpandoracontent/LArHelpers/LArHitWidthHelper.h"
 
#include "larpandoracontent/LArObjects/LArTwoDSlidingFitResult.h"
 
#include "larpandoracontent/LArShowerRefinement/ShowerSpineFinderTool.h"
 
using namespace pandora;
 
namespace lar_content
{
 
ShowerSpineFinderTool::ShowerSpineFinderTool() :
    m_hitThresholdForSpine(7),
    m_growingFitInitialLength(10.f),
    m_initialFitDistanceToLine(0.5f),
    m_minInitialHitsFound(7),
    m_maxFittingHits(15),
    m_localSlidingFitWindow(10),
    m_growingFitSegmentLength(5.f),
    m_highResolutionSlidingFitWindow(5),
    m_distanceToLine(0.75f),
    m_hitConnectionDistance(0.75f)
{
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode ShowerSpineFinderTool::Run(const CartesianVector &nuVertex3D, const CaloHitList *const pViewHitList, const HitType hitType,
    const CartesianVector &peakDirection, CaloHitList &unavailableHitList, CaloHitList &showerSpineHitList)
{
    const CartesianVector nuVertex2D(LArGeometryHelper::ProjectPosition(this->GetPandora(), nuVertex3D, hitType));
 
    this->FindShowerSpine(pViewHitList, nuVertex2D, peakDirection, unavailableHitList, showerSpineHitList);
 
    // Demand that spine is significant, be lenient here as some have small stubs and a gap
    if (showerSpineHitList.size() < m_hitThresholdForSpine)
        return STATUS_CODE_NOT_FOUND;
 
    return STATUS_CODE_SUCCESS;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void ShowerSpineFinderTool::FindShowerSpine(const CaloHitList *const pViewHitList, const CartesianVector &nuVertex2D,
    const CartesianVector &initialDirection, CaloHitList &unavailableHitList, CaloHitList &showerSpineHitList) const
{
    // Use initial direction to find seed hits for a starting fit
    float highestL(0.f);
    CartesianPointVector runningFitPositionVector;
 
    for (const CaloHit *const pCaloHit : *pViewHitList)
    {
        const CartesianVector &hitPosition(pCaloHit->GetPositionVector());
        const CartesianVector &displacementVector(hitPosition - nuVertex2D);
        const float l(initialDirection.GetDotProduct(displacementVector));
        const float t(initialDirection.GetCrossProduct(displacementVector).GetMagnitude());
 
        if ((l < m_growingFitInitialLength) && (l > 0.f) && (t < m_initialFitDistanceToLine))
        {
            if (l > highestL)
                highestL = l;
 
            if (std::find(unavailableHitList.begin(), unavailableHitList.end(), pCaloHit) == unavailableHitList.end())
            {
                showerSpineHitList.push_back(pCaloHit);
            }
 
            runningFitPositionVector.push_back(hitPosition);
        }
    }
 
    // Require significant number of initial hits
    if (runningFitPositionVector.size() < m_minInitialHitsFound)
    {
        showerSpineHitList.clear();
        return;
    }
 
    // Perform a running fit to collect a pathway of hits
    unsigned int count(0);
    bool hitsCollected(true);
    const bool isEndDownstream(initialDirection.GetZ() > 0.f);
    const float slidingFitPitch(LArGeometryHelper::GetWirePitch(this->GetPandora(), showerSpineHitList.front()->GetHitType()));
    CartesianVector extrapolatedStartPosition(nuVertex2D), extrapolatedDirection(initialDirection);
    CartesianVector extrapolatedEndPosition(extrapolatedStartPosition + (extrapolatedDirection * highestL));
 
    while (hitsCollected)
    {
        ++count;
 
        try
        {
            const int excessHitsInFit(runningFitPositionVector.size() - m_maxFittingHits);
 
            if (excessHitsInFit > 0)
            {
                // Remove furthest away hits
                std::sort(runningFitPositionVector.begin(), runningFitPositionVector.end(),
                    LArConnectionPathwayHelper::SortByDistanceToPoint(extrapolatedEndPosition));
 
                for (int i = 0; i < excessHitsInFit; ++i)
                    runningFitPositionVector.erase(std::prev(runningFitPositionVector.end()));
            }
 
            const TwoDSlidingFitResult extrapolatedFit(&runningFitPositionVector, m_localSlidingFitWindow, slidingFitPitch);
 
            extrapolatedStartPosition =
                count == 1 ? extrapolatedEndPosition
                           : isEndDownstream ? extrapolatedFit.GetGlobalMaxLayerPosition() : extrapolatedFit.GetGlobalMinLayerPosition();
            extrapolatedDirection =
                isEndDownstream ? extrapolatedFit.GetGlobalMaxLayerDirection() : extrapolatedFit.GetGlobalMinLayerDirection() * (-1.f);
            extrapolatedEndPosition = extrapolatedStartPosition + (extrapolatedDirection * m_growingFitSegmentLength);
 
            hitsCollected = this->CollectSubsectionHits(extrapolatedFit, extrapolatedStartPosition, extrapolatedEndPosition,
                extrapolatedDirection, isEndDownstream, pViewHitList, runningFitPositionVector, unavailableHitList, showerSpineHitList);
 
            // If no hits found, as a final effort, reduce the sliding fit window
            if (!hitsCollected)
            {
                const TwoDSlidingFitResult microExtrapolatedFit(&runningFitPositionVector, m_highResolutionSlidingFitWindow, slidingFitPitch);
 
                extrapolatedStartPosition = count == 1 ? extrapolatedStartPosition
                                                       : isEndDownstream ? microExtrapolatedFit.GetGlobalMaxLayerPosition()
                                                                         : microExtrapolatedFit.GetGlobalMinLayerPosition();
                extrapolatedDirection = isEndDownstream ? microExtrapolatedFit.GetGlobalMaxLayerDirection()
                                                        : microExtrapolatedFit.GetGlobalMinLayerDirection() * (-1.f);
                extrapolatedEndPosition = extrapolatedStartPosition + (extrapolatedDirection * m_growingFitSegmentLength);
 
                hitsCollected = this->CollectSubsectionHits(microExtrapolatedFit, extrapolatedStartPosition, extrapolatedEndPosition,
                    extrapolatedDirection, isEndDownstream, pViewHitList, runningFitPositionVector, unavailableHitList, showerSpineHitList);
            }
        }
        catch (const StatusCodeException &)
        {
            return;
        }
    }
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool ShowerSpineFinderTool::CollectSubsectionHits(const TwoDSlidingFitResult &extrapolatedFit,
    const CartesianVector &extrapolatedStartPosition, const CartesianVector &extrapolatedEndPosition,
    const CartesianVector &extrapolatedDirection, const bool isEndDownstream, const CaloHitList *const pViewHitList,
    CartesianPointVector &runningFitPositionVector, CaloHitList &unavailableHitList, CaloHitList &showerSpineHitList) const
{
    float extrapolatedStartL(0.f), extrapolatedStartT(0.f);
    extrapolatedFit.GetLocalPosition(extrapolatedStartPosition, extrapolatedStartL, extrapolatedStartT);
 
    float extrapolatedEndL(0.f), extrapolatedEndT(0.f);
    extrapolatedFit.GetLocalPosition(extrapolatedEndPosition, extrapolatedEndL, extrapolatedEndT);
 
    CaloHitList collectedHits;
 
    for (const CaloHit *const pCaloHit : *pViewHitList)
    {
        if (std::find(showerSpineHitList.begin(), showerSpineHitList.end(), pCaloHit) != showerSpineHitList.end())
            continue;
 
        if (std::find(unavailableHitList.begin(), unavailableHitList.end(), pCaloHit) != unavailableHitList.end())
            continue;
 
        const CartesianVector &hitPosition(pCaloHit->GetPositionVector());
 
        float hitL(0.f), hitT(0.f);
        extrapolatedFit.GetLocalPosition(hitPosition, hitL, hitT);
 
        // Assess whether hit is within section boundaries
        if (isEndDownstream && ((hitL < extrapolatedStartL) || (hitL > extrapolatedEndL)))
            continue;
 
        if (!isEndDownstream && ((hitL > extrapolatedStartL) || (hitL < extrapolatedEndL)))
            continue;
 
        // Assess whether hit is close to connecting line - taking account hit width if necessary
        if (this->IsCloseToLine(hitPosition, extrapolatedStartPosition, extrapolatedDirection, m_distanceToLine))
        {
            collectedHits.push_back(pCaloHit);
        }
        else
        {
            const CartesianVector closestPointInHit(
                LArHitWidthHelper::GetClosestPointToLine2D(extrapolatedStartPosition, extrapolatedDirection, pCaloHit));
 
            if (this->IsCloseToLine(closestPointInHit, extrapolatedStartPosition, extrapolatedDirection, m_distanceToLine))
                collectedHits.push_back(pCaloHit);
        }
    }
 
    const int nInitialHits(showerSpineHitList.size());
 
    // Now find a continuous path of collected hits
    this->CollectConnectedHits(collectedHits, extrapolatedStartPosition, extrapolatedDirection, runningFitPositionVector, showerSpineHitList);
 
    const int nFinalHits(showerSpineHitList.size());
 
    return (nFinalHits != nInitialHits);
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
bool ShowerSpineFinderTool::IsCloseToLine(const CartesianVector &hitPosition, const CartesianVector &lineStart,
    const CartesianVector &lineDirection, const float distanceToLine) const
{
    const float transverseDistanceFromLine(lineDirection.GetCrossProduct(hitPosition - lineStart).GetMagnitude());
 
    if (transverseDistanceFromLine > distanceToLine)
        return false;
 
    return true;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void ShowerSpineFinderTool::CollectConnectedHits(const CaloHitList &collectedHits, const CartesianVector &extrapolatedStartPosition,
    const CartesianVector &extrapolatedDirection, CartesianPointVector &runningFitPositionVector, CaloHitList &showerSpineHitList) const
{
    bool found(true);
 
    while (found)
    {
        found = false;
 
        for (const CaloHit *const pCaloHit : collectedHits)
        {
            if (std::find(showerSpineHitList.begin(), showerSpineHitList.end(), pCaloHit) != showerSpineHitList.end())
                continue;
 
            CartesianVector hitPosition(pCaloHit->GetPositionVector());
 
            if (this->GetClosestDistance(hitPosition, runningFitPositionVector) > m_hitConnectionDistance)
            {
                if (LArHitWidthHelper::GetClosestDistance(pCaloHit, showerSpineHitList) > m_hitConnectionDistance)
                    continue;
 
                hitPosition = LArHitWidthHelper::GetClosestPointToLine2D(extrapolatedStartPosition, extrapolatedDirection, pCaloHit);
            }
 
            found = true;
 
            runningFitPositionVector.push_back(hitPosition);
            showerSpineHitList.push_back(pCaloHit);
        }
    }
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
float ShowerSpineFinderTool::GetClosestDistance(const CartesianVector &position, const CartesianPointVector &testPositions) const
{
    float closestDistanceSqaured(std::numeric_limits<float>::max());
 
    for (const CartesianVector &testPosition : testPositions)
    {
        const float separationSquared((testPosition - position).GetMagnitudeSquared());
 
        if (separationSquared < closestDistanceSqaured)
            closestDistanceSqaured = separationSquared;
    }
 
    return std::sqrt(closestDistanceSqaured);
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode ShowerSpineFinderTool::ReadSettings(const TiXmlHandle xmlHandle)
{
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "HitThresholdForSpine", m_hitThresholdForSpine));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
        XmlHelper::ReadValue(xmlHandle, "GrowingFitInitialLength", m_growingFitInitialLength));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
        XmlHelper::ReadValue(xmlHandle, "InitialFitDistanceToLine", m_initialFitDistanceToLine));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinInitialHitsFound", m_minInitialHitsFound));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxFittingHits", m_maxFittingHits));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "LocalSlidingFitWindow", m_localSlidingFitWindow));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
        XmlHelper::ReadValue(xmlHandle, "GrowingFitSegmentLength", m_growingFitSegmentLength));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
        XmlHelper::ReadValue(xmlHandle, "HighResolutionSlidingFitWindow", m_highResolutionSlidingFitWindow));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "DistanceToLine", m_distanceToLine));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "HitConnectionDistance", m_hitConnectionDistance));
 
    return STATUS_CODE_SUCCESS;
}
 
} // namespace lar_content