CheatingNeutrinoCreationAlgorithm.cc

Go to the documentation of this file.

#include "Pandora/AlgorithmHeaders.h"

#include "larpandoracontent/LArHelpers/LArMCParticleHelper.h"
#include "larpandoracontent/LArHelpers/LArPfoHelper.h"

#include "larpandoracontent/LArCheating/CheatingNeutrinoCreationAlgorithm.h"

using namespace pandora;

namespace lar_content
{

CheatingNeutrinoCreationAlgorithm::CheatingNeutrinoCreationAlgorithm() :
    m_collapseToPrimaryMCParticles(false),
    m_vertexTolerance(0.5f)
{
}

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

StatusCode CheatingNeutrinoCreationAlgorithm::Run()
{
    MCParticleVector mcNeutrinoVector;
    this->GetMCNeutrinoVector(mcNeutrinoVector);

    for (const MCParticle *const pMCNeutrino : mcNeutrinoVector)
    {
        const ParticleFlowObject *pNeutrinoPfo(nullptr);
        this->CreateAndSaveNeutrinoPfo(pMCNeutrino, pNeutrinoPfo);

        if (!pNeutrinoPfo)
            return STATUS_CODE_FAILURE;

        this->AddNeutrinoVertex(pMCNeutrino, pNeutrinoPfo);

        MCParticleToPfoMap mcParticleToPfoMap;
        this->GetMCParticleToDaughterPfoMap(mcParticleToPfoMap);
        this->CreatePfoHierarchy(pMCNeutrino, pNeutrinoPfo, mcParticleToPfoMap);
    }

    return STATUS_CODE_SUCCESS;
}

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

void CheatingNeutrinoCreationAlgorithm::GetMCNeutrinoVector(MCParticleVector &mcNeutrinoVector) const
{
    const MCParticleList *pMCParticleList(nullptr);
    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetList(*this, m_mcParticleListName, pMCParticleList));

    MCParticleVector allMCNeutrinoVector;
    LArMCParticleHelper::GetTrueNeutrinos(pMCParticleList, allMCNeutrinoVector);

    for (const MCParticle *const pMCNeutrino : allMCNeutrinoVector)
    {
        if (LArMCParticleHelper::IsNeutrino(pMCNeutrino) && (MC_3D == pMCNeutrino->GetMCParticleType()))
            mcNeutrinoVector.push_back(pMCNeutrino);
    }
}

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

void CheatingNeutrinoCreationAlgorithm::CreateAndSaveNeutrinoPfo(const MCParticle *const pMCNeutrino, const ParticleFlowObject *&pNeutrinoPfo) const
{
    pNeutrinoPfo = nullptr;

    PandoraContentApi::ParticleFlowObject::Parameters pfoParameters;
    pfoParameters.m_particleId = pMCNeutrino->GetParticleId();
    pfoParameters.m_charge = PdgTable::GetParticleCharge(pfoParameters.m_particleId.Get());
    pfoParameters.m_mass = PdgTable::GetParticleMass(pfoParameters.m_particleId.Get());
    pfoParameters.m_energy = pMCNeutrino->GetEnergy();
    pfoParameters.m_momentum = pMCNeutrino->GetMomentum();

    std::string neutrinoPfoListName;
    const PfoList *pNeutrinoPfoList(nullptr);
    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::CreateTemporaryListAndSetCurrent(*this, pNeutrinoPfoList, neutrinoPfoListName));
    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::Create(*this, pfoParameters, pNeutrinoPfo));

    if (!pNeutrinoPfoList || pNeutrinoPfoList->empty())
        throw StatusCodeException(STATUS_CODE_FAILURE);

    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SaveList<Pfo>(*this, m_neutrinoPfoListName));
    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ReplaceCurrentList<Pfo>(*this, m_neutrinoPfoListName));
}

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

void CheatingNeutrinoCreationAlgorithm::AddNeutrinoVertex(const MCParticle *const pMCNeutrino, const ParticleFlowObject *const pNeutrinoPfo) const
{
    const VertexList *pVertexList(nullptr);
    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::GetCurrentList(*this, pVertexList));

    const Vertex *pNeutrinoVertex(nullptr);
    float closestVertexDistance(std::numeric_limits<float>::max());

    for (const Vertex *const pVertex : *pVertexList)
    {
        const float distance((pVertex->GetPosition() - pMCNeutrino->GetEndpoint()).GetMagnitude());

        if (distance < closestVertexDistance)
        {
            pNeutrinoVertex = pVertex;
            closestVertexDistance = distance;
        }
    }

    if (!pNeutrinoVertex || (VERTEX_3D != pNeutrinoVertex->GetVertexType()) || ((pNeutrinoVertex->GetPosition() - pMCNeutrino->GetEndpoint()).GetMagnitude() > m_vertexTolerance))
        throw StatusCodeException(STATUS_CODE_INVALID_PARAMETER);

    PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::AddToPfo(*this, pNeutrinoPfo, pNeutrinoVertex));
}

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

void CheatingNeutrinoCreationAlgorithm::GetMCParticleToDaughterPfoMap(MCParticleToPfoMap &mcParticleToPfoMap) const
{
    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);
    }

    for (const std::string &daughterPfoListName : m_daughterPfoListNames)
    {
        const PfoList *pDaughterPfoList(nullptr);

        if (STATUS_CODE_SUCCESS != PandoraContentApi::GetList(*this, daughterPfoListName, pDaughterPfoList))
        {
            if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
                std::cout << "CheatingNeutrinoCreationAlgorithm: pfo list " << daughterPfoListName << " unavailable." << std::endl;

            continue;
        }

        for (const ParticleFlowObject *const pDaughterPfo : *pDaughterPfoList)
        {
            const MCParticle *pMCParticle(LArMCParticleHelper::GetMainMCParticle(pDaughterPfo));

            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 (!mcParticleToPfoMap.insert(MCParticleToPfoMap::value_type(pMCParticle, pDaughterPfo)).second)
                throw StatusCodeException(STATUS_CODE_OUT_OF_RANGE);
        }
    }
}

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

void CheatingNeutrinoCreationAlgorithm::CreatePfoHierarchy(const MCParticle *const pParentMCParticle, const ParticleFlowObject *const pParentPfo,
    const MCParticleToPfoMap &mcParticleToPfoMap) const
{
    for (const MCParticle *const pDaughterMCParticle : pParentMCParticle->GetDaughterList())
    {
        MCParticleToPfoMap::const_iterator mapIter = mcParticleToPfoMap.find(pDaughterMCParticle);

        if (mcParticleToPfoMap.end() == mapIter)
        {
            this->CreatePfoHierarchy(pDaughterMCParticle, pParentPfo, mcParticleToPfoMap);
            continue;
        }

        const ParticleFlowObject *const pDaughterPfo(mapIter->second);
        PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::SetPfoParentDaughterRelationship(*this, pParentPfo, pDaughterPfo));
        this->CreatePfoHierarchy(pDaughterMCParticle, pDaughterPfo, mcParticleToPfoMap);
    }
}

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

StatusCode CheatingNeutrinoCreationAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
{
    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
        "CollapseToPrimaryMCParticles", m_collapseToPrimaryMCParticles));

    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
        "MCParticleListName", m_mcParticleListName));

    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
        "NeutrinoPfoListName", m_neutrinoPfoListName));

    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
        "VertexListName", m_vertexListName));

    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadVectorOfValues(xmlHandle,
        "DaughterPfoListNames", m_daughterPfoListNames));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
        "VertexTolerance", m_vertexTolerance));

    return STATUS_CODE_SUCCESS;
}

} // namespace lar_content