CosmicRayShowerMatchingAlgorithm.cc

Go to the documentation of this file.

#include "Pandora/AlgorithmHeaders.h"

#include "larpandoracontent/LArThreeDReco/LArCosmicRay/CosmicRayShowerMatchingAlgorithm.h"

#include "larpandoracontent/LArHelpers/LArClusterHelper.h"
#include "larpandoracontent/LArHelpers/LArGeometryHelper.h"

using namespace pandora;

namespace lar_content
{

CosmicRayShowerMatchingAlgorithm::CosmicRayShowerMatchingAlgorithm() :
    m_minCaloHitsPerCluster(10),
    m_minXOverlap(1.f),
    m_minXOverlapFraction(0.5f),
    m_pseudoChi2Cut(5.f)
{
}

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

void CosmicRayShowerMatchingAlgorithm::SelectCleanClusters(const ClusterVector &inputVector, ClusterVector &outputVector) const
{
    for (const Cluster *const pCluster : inputVector)
    {
        if (pCluster->GetNCaloHits() < m_minCaloHitsPerCluster)
            continue;

        outputVector.push_back(pCluster);
    }
}

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

bool CosmicRayShowerMatchingAlgorithm::MatchClusters(const Cluster *const pCluster1, const Cluster *const pCluster2) const
{
    float xMin1(0.f), xMax1(0.f), xMin2(0.f), xMax2(0.f);
    pCluster1->GetClusterSpanX(xMin1, xMax1);
    pCluster2->GetClusterSpanX(xMin2, xMax2);

    const float xOverlap(std::min(xMax1, xMax2) - std::max(xMin1, xMin2));
    const float xSpan(std::max(xMax1, xMax2) - std::min(xMin1, xMin2));

    if (xSpan < std::numeric_limits<float>::epsilon())
        return false;

    if (xOverlap > m_minXOverlap && xOverlap / xSpan > m_minXOverlapFraction)
        return true;

    return false;
}

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

bool CosmicRayShowerMatchingAlgorithm::CheckMatchedClusters3D(
    const Cluster *const pCluster1, const Cluster *const pCluster2, const Cluster *const pCluster3) const
{
    // Check that three clusters have a consistent 3D position
    const HitType hitType1(LArClusterHelper::GetClusterHitType(pCluster1));
    const HitType hitType2(LArClusterHelper::GetClusterHitType(pCluster2));
    const HitType hitType3(LArClusterHelper::GetClusterHitType(pCluster3));

    if (hitType1 == hitType2 || hitType2 == hitType3 || hitType3 == hitType1)
        throw StatusCodeException(STATUS_CODE_FAILURE);

    // Requirements on X matching
    float xMin1(0.f), xMin2(0.f), xMin3(0.f), xMax1(0.f), xMax2(0.f), xMax3(0.f);
    pCluster1->GetClusterSpanX(xMin1, xMax1);
    pCluster2->GetClusterSpanX(xMin2, xMax2);
    pCluster3->GetClusterSpanX(xMin3, xMax3);

    const float xMin(std::max(xMin1, std::max(xMin2, xMin3)));
    const float xMax(std::min(xMax1, std::min(xMax2, xMax3)));
    const float xOverlap(xMax - xMin);

    if (xOverlap < std::numeric_limits<float>::epsilon())
        return false;

    float p1(std::numeric_limits<float>::max()), p2(std::numeric_limits<float>::max()), p3(std::numeric_limits<float>::max());

    if ((STATUS_CODE_SUCCESS != LArClusterHelper::GetAverageZ(pCluster1, xMin, xMax, p1)) ||
        (STATUS_CODE_SUCCESS != LArClusterHelper::GetAverageZ(pCluster2, xMin, xMax, p2)) ||
        (STATUS_CODE_SUCCESS != LArClusterHelper::GetAverageZ(pCluster3, xMin, xMax, p3)))
    {
        return false;
    }

    const float q3(LArGeometryHelper::MergeTwoPositions(this->GetPandora(), hitType1, hitType2, p1, p2));
    const float q1(LArGeometryHelper::MergeTwoPositions(this->GetPandora(), hitType2, hitType3, p2, p3));
    const float q2(LArGeometryHelper::MergeTwoPositions(this->GetPandora(), hitType3, hitType1, p3, p1));

    const float pseudoChi2(((q1 - p1) * (q1 - p1) + (q2 - p2) * (q2 - p2) + (q3 - p3) * (q3 - p3)) / 3.f);

    if (pseudoChi2 > m_pseudoChi2Cut)
        return false;

    return true;
}

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

void CosmicRayShowerMatchingAlgorithm::SetPfoParameters(const Particle &particle, PandoraContentApi::ParticleFlowObject::Parameters &pfoParameters) const
{
    ClusterList clusterList;
    if (particle.m_pClusterU)
        clusterList.push_back(particle.m_pClusterU);
    if (particle.m_pClusterV)
        clusterList.push_back(particle.m_pClusterV);
    if (particle.m_pClusterW)
        clusterList.push_back(particle.m_pClusterW);

    // TODO Correct these placeholder parameters
    pfoParameters.m_particleId = E_MINUS; // SHOWER
    pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
    pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
    pfoParameters.m_energy = 0.f;
    pfoParameters.m_momentum = CartesianVector(0.f, 0.f, 0.f);
    pfoParameters.m_clusterList = clusterList;
}

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

StatusCode CosmicRayShowerMatchingAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
{
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinCaloHitsPerCluster", m_minCaloHitsPerCluster));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinXOverlap", m_minXOverlap));

    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinXOverlapFraction", m_minXOverlapFraction));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "PseudoChi2Cut", m_pseudoChi2Cut));

    return CosmicRayBaseMatchingAlgorithm::ReadSettings(xmlHandle);
}

} // namespace lar_content