LayerSplittingAlgorithm.cc

Go to the documentation of this file.

 
#include "Pandora/AlgorithmHeaders.h"
 
#include "larpandoracontent/LArTwoDReco/LArClusterSplitting/LayerSplittingAlgorithm.h"
 
using namespace pandora;
 
namespace lar_content
{
 
LayerSplittingAlgorithm::LayerSplittingAlgorithm() :
    m_minClusterLayers(20),
    m_layerWindow(10),
    m_maxScatterRms(0.35f),
    m_maxScatterCosTheta(0.5f),
    m_maxSlidingCosTheta(0.866f)
{
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode LayerSplittingAlgorithm::DivideCaloHits(const Cluster *const pCluster, CaloHitList &firstHitList, CaloHitList &secondHitList) const
{
    unsigned int splitLayer(0);
 
    if (STATUS_CODE_SUCCESS == this->FindBestSplitLayer(pCluster, splitLayer))
        return this->DivideCaloHits(pCluster, splitLayer, firstHitList, secondHitList);
 
    return STATUS_CODE_NOT_FOUND;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode LayerSplittingAlgorithm::FindBestSplitLayer(const Cluster *const pCluster, unsigned int &splitLayer) const
{
    const OrderedCaloHitList &orderedCaloHitList(pCluster->GetOrderedCaloHitList());
 
    if (orderedCaloHitList.size() < m_minClusterLayers)
        return STATUS_CODE_NOT_FOUND;
 
    bool foundSplit(false);
 
    float bestCosTheta(1.f);
    CartesianVector bestPosition(0.f, 0.f, 0.f);
 
    for (unsigned int iLayer = pCluster->GetInnerPseudoLayer() + 4; iLayer + 4 <= pCluster->GetOuterPseudoLayer(); ++iLayer)
    {
        if (orderedCaloHitList.find(iLayer) == orderedCaloHitList.end())
            continue;
 
        unsigned int innerLayer((pCluster->GetInnerPseudoLayer() + m_layerWindow > iLayer) ? pCluster->GetInnerPseudoLayer() : iLayer - m_layerWindow);
        unsigned int outerLayer((iLayer + m_layerWindow > pCluster->GetOuterPseudoLayer()) ? pCluster->GetOuterPseudoLayer() : iLayer + m_layerWindow);
 
        for (; innerLayer >= pCluster->GetInnerPseudoLayer(); --innerLayer)
        {
            if (orderedCaloHitList.find(innerLayer) != orderedCaloHitList.end())
                break;
        }
 
        for (; outerLayer <= pCluster->GetOuterPseudoLayer(); ++outerLayer)
        {
            if (orderedCaloHitList.find(outerLayer) != orderedCaloHitList.end())
                break;
        }
 
        const CartesianVector splitPosition(pCluster->GetCentroid(iLayer));
        const CartesianVector innerPosition(pCluster->GetCentroid(innerLayer));
        const CartesianVector outerPosition(pCluster->GetCentroid(outerLayer));
 
        const CartesianVector r1(innerPosition - splitPosition);
        const CartesianVector r2(outerPosition - splitPosition);
        const CartesianVector p1(r1.GetUnitVector());
        const CartesianVector p2(r2.GetUnitVector());
 
        const float cosTheta(-p1.GetDotProduct(p2));
        const float rms1(this->CalculateRms(pCluster, innerLayer, iLayer));
        const float rms2(this->CalculateRms(pCluster, outerLayer, iLayer));
        const float rms(std::max(rms1, rms2));
 
        float rmsCut(std::numeric_limits<float>::max());
 
        if (cosTheta > 0.f)
        {
            rmsCut = m_maxScatterRms;
 
            if (cosTheta > m_maxScatterCosTheta)
            {
                rmsCut *= ((m_maxSlidingCosTheta > cosTheta) ? (m_maxSlidingCosTheta - cosTheta) / (m_maxSlidingCosTheta - m_maxScatterCosTheta) : 0.f);
            }
        }
 
        if (rms < rmsCut && cosTheta < bestCosTheta)
        {
            bestCosTheta = cosTheta;
            bestPosition = splitPosition;
 
            splitLayer = iLayer;
            foundSplit = true;
        }
    }
 
    if (!foundSplit)
        return STATUS_CODE_NOT_FOUND;
 
    return STATUS_CODE_SUCCESS;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
float LayerSplittingAlgorithm::CalculateRms(const Cluster *const pCluster, const unsigned int &firstLayer, const unsigned int &secondLayer) const
{
    const OrderedCaloHitList &orderedCaloHitList(pCluster->GetOrderedCaloHitList());
 
    const unsigned int innerLayer(std::min(firstLayer, secondLayer));
    const unsigned int outerLayer(std::max(firstLayer, secondLayer));
 
    const CartesianVector innerPosition(pCluster->GetCentroid(innerLayer));
    const CartesianVector outerPosition(pCluster->GetCentroid(outerLayer));
    const CartesianVector predictedDirection((outerPosition - innerPosition).GetUnitVector());
 
    float totalChi2(0.f);
    float totalLayers(0.f);
 
    for (unsigned int iLayer = innerLayer + 1; iLayer + 1 < outerLayer; ++iLayer)
    {
        if (orderedCaloHitList.find(iLayer) == orderedCaloHitList.end())
            continue;
 
        const CartesianVector hitPosition(pCluster->GetCentroid(iLayer));
        const CartesianVector predictedPosition(innerPosition + predictedDirection * predictedDirection.GetDotProduct(hitPosition - innerPosition));
 
        totalChi2 += (predictedPosition - hitPosition).GetMagnitudeSquared();
        totalLayers += 1.f;
    }
 
    if (totalLayers > 0.f)
        return std::sqrt(totalChi2 / totalLayers);
 
    return 0.f;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode LayerSplittingAlgorithm::DivideCaloHits(
    const Cluster *const pCluster, const unsigned int &splitLayer, CaloHitList &firstHitList, CaloHitList &secondHitList) const
{
    const OrderedCaloHitList &orderedCaloHitList(pCluster->GetOrderedCaloHitList());
 
    for (OrderedCaloHitList::const_iterator iter = orderedCaloHitList.begin(); iter != orderedCaloHitList.end(); ++iter)
    {
        const unsigned int thisLayer(iter->first);
 
        for (CaloHitList::const_iterator hitIter = iter->second->begin(), hitIterEnd = iter->second->end(); hitIter != hitIterEnd; ++hitIter)
        {
            const CaloHit *const pCaloHit = *hitIter;
 
            if (thisLayer < splitLayer)
            {
                firstHitList.push_back(pCaloHit);
            }
            else
            {
                secondHitList.push_back(pCaloHit);
            }
        }
    }
 
    if (firstHitList.empty() || secondHitList.empty())
        return STATUS_CODE_NOT_FOUND;
 
    return STATUS_CODE_SUCCESS;
}
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
StatusCode LayerSplittingAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
{
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinClusterLayers", m_minClusterLayers));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "LayerWindow", m_layerWindow));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxScatterRms", m_maxScatterRms));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxScatterCosTheta", m_maxScatterCosTheta));
 
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxSlidingCosTheta", m_maxSlidingCosTheta));
 
    return ClusterSplittingAlgorithm::ReadSettings(xmlHandle);
}
 
} // namespace lar_content