CheatingPfoCreationAlgorithm.cc

Go to the documentation of this file.

#include "Pandora/AlgorithmHeaders.h"

#include "larpandoracontent/LArCheating/CheatingPfoCreationAlgorithm.h"

#include "larpandoracontent/LArHelpers/LArClusterHelper.h"

using namespace pandora;

namespace lar_content
{

CheatingPfoCreationAlgorithm::CheatingPfoCreationAlgorithm() :
    m_collapseToPrimaryMCParticles(false),
    m_useOnlyAvailableClusters(true),
    m_addVertices(true),
    m_replaceCurrentVertexList(false),
    m_minGoodHitTypes(0),
    m_nHitsForGoodHitType(10)
{
}

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

StatusCode CheatingPfoCreationAlgorithm::Run()
{
    LArMCParticleHelper::MCRelationMap mcPrimaryMap;

    if (m_collapseToPrimaryMCParticles)
    {
        const MCParticleList *pMCParticleList(nullptr);
        PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetList(*this, m_mcParticleListName, pMCParticleList));

        LArMCParticleHelper::GetMCPrimaryMap(pMCParticleList, mcPrimaryMap);
    }

    MCParticleToClusterListMap mcParticleToClusterListMap;

    for (const std::string &clusterListName : m_inputClusterListNames)
    {
        const ClusterList *pClusterList(nullptr);

        if (STATUS_CODE_SUCCESS != PandoraContentApi::GetList(*this, clusterListName, pClusterList))
        {
            if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
                std::cout << "CheatingPfoCreationAlgorithm - Could not access cluster list with name " << clusterListName << std::endl;

            continue;
        }

        this->GetMCParticleToClusterListMap(pClusterList, mcPrimaryMap, mcParticleToClusterListMap);
    }

    this->CreatePfos(mcParticleToClusterListMap);
    return STATUS_CODE_SUCCESS;
}

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

void CheatingPfoCreationAlgorithm::GetMCParticleToClusterListMap(const ClusterList *const pClusterList,
    const LArMCParticleHelper::MCRelationMap &mcPrimaryMap, MCParticleToClusterListMap &mcParticleToClusterListMap) const
{
    for (const Cluster *const pCluster : *pClusterList)
    {
        try
        {
            if (m_useOnlyAvailableClusters && !PandoraContentApi::IsAvailable(*this, pCluster))
                continue;

            const MCParticle *pMCParticle(MCParticleHelper::GetMainMCParticle(pCluster));

            if (m_collapseToPrimaryMCParticles)
            {
                LArMCParticleHelper::MCRelationMap::const_iterator primaryIter = mcPrimaryMap.find(pMCParticle);

                if (mcPrimaryMap.end() == primaryIter)
                    throw StatusCodeException(STATUS_CODE_NOT_FOUND);

                pMCParticle = primaryIter->second;
            }

            if (!m_particleIdList.empty() && !m_particleIdList.count(pMCParticle->GetParticleId()))
                continue;

            mcParticleToClusterListMap[pMCParticle].push_back(pCluster);
        }
        catch (const StatusCodeException &)
        {
        }
    }
}

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

void CheatingPfoCreationAlgorithm::CreatePfos(const MCParticleToClusterListMap &mcParticleToClusterListMap) const
{
    if (mcParticleToClusterListMap.empty())
        return;

    const PfoList *pPfoList(nullptr);
    std::string pfoListName;
    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*this, pPfoList, pfoListName));

    const VertexList *pVertexList(nullptr);
    std::string vertexListName;
    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*this, pVertexList, vertexListName));

    MCParticleVector mcParticleVector;
    for (const auto &mapEntry : mcParticleToClusterListMap)
        mcParticleVector.push_back(mapEntry.first);
    std::sort(mcParticleVector.begin(), mcParticleVector.end(), LArMCParticleHelper::SortByMomentum);

    for (const MCParticle *const pMCParticle : mcParticleVector)
    {
        const ClusterList &clusterList(mcParticleToClusterListMap.at(pMCParticle));

        if (clusterList.empty())
            continue;

        if (this->GetNHitTypesAboveThreshold(clusterList, m_nHitsForGoodHitType) < m_minGoodHitTypes)
            continue;

        try
        {
            PandoraContentApi::ParticleFlowObject::Parameters pfoParameters;
            pfoParameters.m_particleId = (this->IsShower(pMCParticle) ? E_MINUS : MU_MINUS);
            pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
            pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
            pfoParameters.m_energy = pMCParticle->GetEnergy();
            pfoParameters.m_momentum = pMCParticle->GetMomentum();
            pfoParameters.m_clusterList.insert(pfoParameters.m_clusterList.end(), clusterList.begin(), clusterList.end());

            const ParticleFlowObject *pPfo(nullptr);
            PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*this, pfoParameters, pPfo));

            if (m_addVertices)
            {
                PandoraContentApi::Vertex::Parameters parameters;
                parameters.m_position = pMCParticle->GetVertex();
                parameters.m_vertexLabel = VERTEX_START;
                parameters.m_vertexType = VERTEX_3D;

                const Vertex *pVertex(nullptr);
                PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Vertex::Create(*this, parameters, pVertex));
                PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToPfo<Vertex>(*this, pPfo, pVertex));
            }
        }
        catch (const StatusCodeException &)
        {
            std::cout << "CheatingPfoCreationAlgorithm: Could not create PFO for MCParticle with pdg code " << pMCParticle->GetParticleId()
                      << std::endl;
        }
    }

    if (!pPfoList->empty())
    {
        PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*this, m_outputPfoListName));
        PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ReplaceCurrentList<Pfo>(*this, m_outputPfoListName));
    }

    if (!pVertexList->empty())
    {
        PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Vertex>(*this, m_outputVertexListName));

        if (m_replaceCurrentVertexList)
            PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ReplaceCurrentList<Vertex>(*this, m_outputVertexListName));
    }
}

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

unsigned int CheatingPfoCreationAlgorithm::GetNHitTypesAboveThreshold(const ClusterList &clusterList, const unsigned int nHitsThreshold) const
{
    HitTypeMap hitTypeMap;

    for (const Cluster *const pCluster : clusterList)
    {
        hitTypeMap[LArClusterHelper::GetClusterHitType(pCluster)] += pCluster->GetNCaloHits();
    }

    unsigned int nGoodViews(0);

    for (const HitTypeMap::value_type &mapEntry : hitTypeMap)
    {
        if (mapEntry.second > nHitsThreshold)
            ++nGoodViews;
    }

    return nGoodViews;
}

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

bool CheatingPfoCreationAlgorithm::IsTrack(const MCParticle *const pMCParticle) const
{
    const int pdg(pMCParticle->GetParticleId());

    return ((MU_MINUS == std::abs(pdg)) || (PI_PLUS == std::abs(pdg)) || (PROTON == std::abs(pdg)) || (K_PLUS == std::abs(pdg)));
}

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

bool CheatingPfoCreationAlgorithm::IsShower(const MCParticle *const pMCParticle) const
{
    const int pdg(pMCParticle->GetParticleId());

    return ((E_MINUS == std::abs(pdg)) || (PHOTON == std::abs(pdg)));
}

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

StatusCode CheatingPfoCreationAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
{
    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadVectorOfValues(xmlHandle, "InputClusterListNames", m_inputClusterListNames));

    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "OutputPfoListName", m_outputPfoListName));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
        XmlHelper::ReadValue(xmlHandle, "CollapseToPrimaryMCParticles", m_collapseToPrimaryMCParticles));

    if (m_collapseToPrimaryMCParticles)
    {
        PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "MCParticleListName", m_mcParticleListName));
    }

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
        XmlHelper::ReadValue(xmlHandle, "UseOnlyAvailableClusters", m_useOnlyAvailableClusters));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "AddVertices", m_addVertices));

    if (m_addVertices)
    {
        PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle, "OutputVertexListName", m_outputVertexListName));

        PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=,
            XmlHelper::ReadValue(xmlHandle, "ReplaceCurrentVertexList", m_replaceCurrentVertexList));
    }

    IntVector particleIdVector;
    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadVectorOfValues(xmlHandle, "ParticleIdList", particleIdVector));

    m_particleIdList.insert(particleIdVector.begin(), particleIdVector.end());

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "MinGoodHitTypes", m_minGoodHitTypes));

    PANDORA_RETURN_RESULT_IF_AND_IF(
        STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle, "NHitsForGoodHitType", m_nHitsForGoodHitType));

    return STATUS_CODE_SUCCESS;
}

} // namespace lar_content