CheatingSliceSelectionTool.cc

Go to the documentation of this file.

#include "Pandora/AlgorithmHeaders.h"

#include "larpandoracontent/LArCheating/CheatingSliceSelectionTool.h"

#include "larpandoracontent/LArHelpers/LArMCParticleHelper.h"

using namespace pandora;

namespace lar_content
{

CheatingSliceSelectionTool::CheatingSliceSelectionTool() :
    m_maxSlices{1},
    m_threshold{-1.f},
    m_cutVariable{"completeness"}
{
}

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

void CheatingSliceSelectionTool::SelectSlices(const pandora::Algorithm *const /*pAlgorithm*/, const SliceVector &inputSliceVector, SliceVector &outputSliceVector)
{
    // ATTN Ensure this only runs if slicing enabled
    unsigned int sliceCounter{0};
    MetricSliceIndexMap reducedSliceVarIndexMap;

    FloatVector targetWeights(inputSliceVector.size());
    FloatVector otherWeights(inputSliceVector.size());

    // Calculate target and total weight for each slice
    for (const CaloHitList &sliceHits : inputSliceVector)
    {
        CaloHitList localHitList{};
        // ATTN Must ensure we copy the hit actually owned by master instance; access differs with/without slicing enabled
        for (const CaloHit *const pSliceCaloHit : sliceHits)
            localHitList.push_back(static_cast<const CaloHit *>(pSliceCaloHit->GetParentAddress()));

        for (const CaloHit *const pCaloHit : localHitList)
        {
            float thisTargetParticleWeight{0.f}, thisTotalWeight{0.f};
            const MCParticleWeightMap &hitMCParticleWeightMap(pCaloHit->GetMCParticleWeightMap());

            if (hitMCParticleWeightMap.empty())
                continue;

            MCParticleList mcParticleList;
            for (const auto &mapEntry : hitMCParticleWeightMap)
                mcParticleList.push_back(mapEntry.first);
            mcParticleList.sort(LArMCParticleHelper::SortByMomentum);

            for (const MCParticle *const pMCParticle : mcParticleList)
            {
                const float weight{hitMCParticleWeightMap.at(pMCParticle)};
                const MCParticle *const pParentMCParticle{LArMCParticleHelper::GetParentMCParticle(pMCParticle)};

                if (this->IsTarget(pParentMCParticle))
                    thisTargetParticleWeight += weight;

                thisTotalWeight += weight;
            }

            // ATTN normalise arbitrary input weights at this point
            if (thisTotalWeight > std::numeric_limits<float>::epsilon())
            {
                thisTargetParticleWeight *= 1.f / thisTotalWeight;
                thisTotalWeight = 1.f;
            }
            else
            {
                thisTargetParticleWeight = 0.f;
                thisTotalWeight = 0.f;
            }
            targetWeights[sliceCounter] += thisTargetParticleWeight;
            otherWeights[sliceCounter] += (1. - thisTargetParticleWeight);
        }

        ++sliceCounter;
    }

    float totalTargetWeight{0.f};
    for (const float weight : targetWeights)
        totalTargetWeight += weight;

    // Add slices passing cut threshold to map
    for (unsigned int i = 0; i < targetWeights.size(); ++i)
    {
        const float sliceWeight = targetWeights[i] + otherWeights[i];
        const float completeness{totalTargetWeight > std::numeric_limits<float>::epsilon() ? targetWeights[i] / totalTargetWeight : 0.f};
        const float purity{sliceWeight > std::numeric_limits<float>::epsilon() ? targetWeights[i] / sliceWeight : 0.f};
        // Already checked correctness of variable in ReadSettings
        if (m_cutVariable == "completeness")
        {
            if (completeness >= m_threshold)
                reducedSliceVarIndexMap.emplace(completeness, i);
        }
        else if (m_cutVariable == "purity")
        {
            if (purity >= m_threshold)
                reducedSliceVarIndexMap.emplace(purity, i);
        }
    }
    // Select the best m_maxSlices slices - prefix increment ensures all slices retained if m_maxSlices == 0
    std::vector<int> reducedSliceIndices{};
    int i = 0;
    for (const auto [cutVariable, index] : reducedSliceVarIndexMap)
    {                      // ATTN: Map is sorted on cut variable from max to min
        (void)cutVariable; // GCC 7 support, versions 8+ do not need this
        reducedSliceIndices.push_back(index);
        outputSliceVector.push_back(inputSliceVector[index]);
        if (++i == m_maxSlices)
            break;
    }
}

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

StatusCode CheatingSliceSelectionTool::ReadSettings(const TiXmlHandle xmlHandle)
{
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MaxSlices", m_maxSlices));
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "Threshold", m_threshold));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "CutVariable", m_cutVariable));
    std::transform(m_cutVariable.begin(), m_cutVariable.end(), m_cutVariable.begin(), [](unsigned char c) { return std::tolower(c); });
    if (m_cutVariable != "completeness" && m_cutVariable != "purity")
    {
        std::cout << "CheatingSliceSelectionTool::ReadSettings: Unknown cut variable \'" << m_cutVariable << "\'" << std::endl;
        return STATUS_CODE_INVALID_PARAMETER;
    }

    return STATUS_CODE_SUCCESS;
}

} // namespace lar_content