PeakDirectionFinderTool.cc

Go to the documentation of this file.

 
#include "Pandora/AlgorithmHeaders.h"
#include "Pandora/AlgorithmTool.h"
 
#include "larpandoracontent/LArHelpers/LArGeometryHelper.h"
#include "larpandoracontent/LArHelpers/LArPfoHelper.h"
 
#include "larpandoracontent/LArShowerRefinement/PeakDirectionFinderTool.h"
 
using namespace pandora;
 
namespace lar_content
{
 
PeakDirectionFinderTool::PeakDirectionFinderTool() :
    m_pathwaySearchRegion(10.f),
    m_theta0XZBinSize(0.005f),
    m_smoothingWindow(1),
    m_ambiguousParticleMode(false)
{
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode PeakDirectionFinderTool::Run(const ParticleFlowObject *const pShowerPfo, const CartesianVector &nuVertex3D,
    const CaloHitList *const pViewHitList, const HitType hitType, CartesianPointVector &peakDirectionVector)
{
    CaloHitList viewShowerHitList;
    LArPfoHelper::GetCaloHits(pShowerPfo, hitType, viewShowerHitList);
 
    const CartesianVector nuVertex2D(LArGeometryHelper::ProjectPosition(this->GetPandora(), nuVertex3D, hitType));
 
    // Get event hits in region of interest
    CaloHitList viewROIHits;
    this->CollectHitsWithinROI(viewShowerHitList, pViewHitList, nuVertex2D, viewROIHits);
 
    if (viewROIHits.empty())
        return STATUS_CODE_NOT_FOUND;
 
    // Fill angular decomposition map
    AngularDecompositionMap angularDecompositionMap;
    this->FillAngularDecompositionMap(viewROIHits, nuVertex2D, angularDecompositionMap);
 
    if (angularDecompositionMap.empty())
        return STATUS_CODE_NOT_FOUND;
 
    // Smooth angular decomposition map
    this->SmoothAngularDecompositionMap(angularDecompositionMap);
 
    // Get peak directions
    this->RetrievePeakDirections(angularDecompositionMap, peakDirectionVector);
 
    if (peakDirectionVector.empty())
        return STATUS_CODE_NOT_FOUND;
 
    return STATUS_CODE_SUCCESS;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void PeakDirectionFinderTool::CollectHitsWithinROI(const CaloHitList &showerHitList, const CaloHitList *const pViewHitList,
    const CartesianVector &nuVertex2D, CaloHitList &viewROIHits) const
{
    if (m_ambiguousParticleMode)
    {
        for (const CaloHit *const pCaloHit : *pViewHitList)
        {
            if (std::find(showerHitList.begin(), showerHitList.end(), pCaloHit) != showerHitList.end())
                continue;
 
            viewROIHits.push_back(pCaloHit);
        }
    }
    else
    {
        float lowestTheta(std::numeric_limits<float>::max()), highestTheta((-1.f) * std::numeric_limits<float>::max());
 
        this->GetAngularExtrema(showerHitList, nuVertex2D, lowestTheta, highestTheta);
        this->CollectHitsWithinExtrema(pViewHitList, nuVertex2D, lowestTheta, highestTheta, viewROIHits);
    }
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void PeakDirectionFinderTool::GetAngularExtrema(
    const CaloHitList &showerHitList, const CartesianVector &nuVertex2D, float &lowestTheta, float &highestTheta) const
{
    lowestTheta = std::numeric_limits<float>::max();
    highestTheta = (-1.f) * std::numeric_limits<float>::max();
 
    const CartesianVector xAxis(1.f, 0.f, 0.f);
 
    for (const CaloHit *const pCaloHit : showerHitList)
    {
        const CartesianVector &hitPosition(pCaloHit->GetPositionVector());
        const CartesianVector displacementVector(hitPosition - nuVertex2D);
 
        float theta0XZ(xAxis.GetOpeningAngle(displacementVector));
 
        if (displacementVector.GetZ() < 0.f)
            theta0XZ = (2.0 * M_PI) - theta0XZ;
 
        if (theta0XZ < lowestTheta)
            lowestTheta = theta0XZ;
 
        if (theta0XZ > highestTheta)
            highestTheta = theta0XZ;
    }
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void PeakDirectionFinderTool::CollectHitsWithinExtrema(const CaloHitList *const pViewHitList, const CartesianVector &nuVertex2D,
    const float lowestTheta, const float highestTheta, CaloHitList &viewROIHits) const
{
    const CartesianVector xAxis(1.f, 0.f, 0.f);
 
    for (const CaloHit *const pCaloHit : *pViewHitList)
    {
        const CartesianVector &hitPosition(pCaloHit->GetPositionVector());
        const CartesianVector displacementVector(hitPosition - nuVertex2D);
 
        float theta0XZ(xAxis.GetOpeningAngle(displacementVector));
 
        if (displacementVector.GetZ() < 0.f)
            theta0XZ = (2.0 * M_PI) - theta0XZ;
 
        if ((theta0XZ < lowestTheta) || (theta0XZ > highestTheta))
            continue;
 
        viewROIHits.push_back(pCaloHit);
    }
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void PeakDirectionFinderTool::FillAngularDecompositionMap(
    const CaloHitList &viewShowerHitList, const CartesianVector &nuVertex2D, AngularDecompositionMap &angularDecompositionMap) const
{
    const CartesianVector xAxis(1.f, 0.f, 0.f);
 
    for (const CaloHit *const pCaloHit : viewShowerHitList)
    {
        const CartesianVector &hitPosition(pCaloHit->GetPositionVector());
        const CartesianVector displacementVector(hitPosition - nuVertex2D);
 
        if (displacementVector.GetMagnitudeSquared() > (m_pathwaySearchRegion * m_pathwaySearchRegion))
            continue;
 
        float theta0XZ(xAxis.GetOpeningAngle(displacementVector));
 
        if (displacementVector.GetZ() < 0.f)
            theta0XZ *= -1.f;
 
        const int theta0XZFactor(std::floor(theta0XZ / m_theta0XZBinSize));
 
        if (angularDecompositionMap.find(theta0XZFactor) == angularDecompositionMap.end())
            angularDecompositionMap[theta0XZFactor] = 1;
        else
            angularDecompositionMap[theta0XZFactor] += 1;
    }
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void PeakDirectionFinderTool::SmoothAngularDecompositionMap(AngularDecompositionMap &angularDecompositionMap) const
{
    const int loopMin((-1) * m_smoothingWindow);
    const int loopMax(m_smoothingWindow);
 
    const AngularDecompositionMap angularDecompositionMapTemp(angularDecompositionMap);
 
    angularDecompositionMap.clear();
 
    for (const auto &entry : angularDecompositionMapTemp)
    {
        float total(0.f);
        const int currentBin(entry.first);
 
        for (int binOffset = loopMin; binOffset <= loopMax; ++binOffset)
        {
            const int contributingBin = currentBin + binOffset;
            total += (angularDecompositionMapTemp.find(contributingBin) == angularDecompositionMapTemp.end())
                         ? 0.f
                         : angularDecompositionMapTemp.at(contributingBin);
        }
 
        angularDecompositionMap[currentBin] = total / static_cast<float>((2.0 * m_smoothingWindow) + 1);
    }
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
void PeakDirectionFinderTool::RetrievePeakDirections(const AngularDecompositionMap &angularDecompositionMap, CartesianPointVector &peakDirectionVector) const
{
    IntVector orderedBinIndexVector;
 
    for (const auto &entry : angularDecompositionMap)
        orderedBinIndexVector.push_back(entry.first);
 
    // Order peak bin vector from highest to lowest bin height
    // Tie-break: highest index wins
    std::sort(orderedBinIndexVector.begin(), orderedBinIndexVector.end(), [&angularDecompositionMap](const int a, const int b) -> bool {
        const float aWeight(angularDecompositionMap.at(a));
        const float bWeight(angularDecompositionMap.at(b));
 
        if (std::fabs(aWeight - bWeight) < std::numeric_limits<float>::epsilon())
            return a > b;
        else
            return aWeight > bWeight;
    });
 
    for (int binIndex : orderedBinIndexVector)
    {
        const float binWeight(angularDecompositionMap.at(binIndex));
 
        int precedingBin(binIndex - 1);
        bool foundPreceeding(false);
 
        while (!foundPreceeding)
        {
            if ((angularDecompositionMap.find(precedingBin) == angularDecompositionMap.end()) ||
                (std::fabs(angularDecompositionMap.at(precedingBin) - angularDecompositionMap.at(binIndex)) > std::numeric_limits<float>::epsilon()))
            {
                foundPreceeding = true;
                break;
            }
 
            --precedingBin;
        }
 
        int followingBin(binIndex + 1);
        bool foundFollowing(false);
 
        while (!foundFollowing)
        {
            if ((angularDecompositionMap.find(followingBin) == angularDecompositionMap.end()) ||
                (std::fabs(angularDecompositionMap.at(followingBin) - angularDecompositionMap.at(binIndex)) > std::numeric_limits<float>::epsilon()))
            {
                foundFollowing = true;
                break;
            }
 
            ++followingBin;
        }
 
        const float precedingBinWeight(
            angularDecompositionMap.find(precedingBin) == angularDecompositionMap.end() ? 0.f : angularDecompositionMap.at(precedingBin));
        const float followingBinWeight(
            angularDecompositionMap.find(followingBin) == angularDecompositionMap.end() ? 0.f : angularDecompositionMap.at(followingBin));
 
        if ((binWeight < precedingBinWeight) || (binWeight < followingBinWeight))
            continue;
 
        // Convert to a direction
        const float theta0XZ(binIndex * m_theta0XZBinSize);
        const CartesianVector peakDirection(std::cos(theta0XZ), 0.f, std::sin(theta0XZ));
 
        peakDirectionVector.push_back(peakDirection);
    }
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode PeakDirectionFinderTool::ReadSettings(const TiXmlHandle xmlHandle)
{
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "PathwaySearchRegion", m_pathwaySearchRegion));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "Theta0XZBinSize", m_theta0XZBinSize));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "SmoothingWindow", m_smoothingWindow));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "AmbiguousParticleMode", m_ambiguousParticleMode));
 
    return STATUS_CODE_SUCCESS;
}
 
} // namespace lar_content