NeutrinoPropertiesAlgorithm.cc

Go to the documentation of this file.

#include "Pandora/AlgorithmHeaders.h"

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

#include "larpandoracontent/LArThreeDReco/LArEventBuilding/NeutrinoPropertiesAlgorithm.h"

using namespace pandora;

namespace lar_content
{

NeutrinoPropertiesAlgorithm::NeutrinoPropertiesAlgorithm() :
    m_includeIsolatedHits(false)
{
}

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

StatusCode NeutrinoPropertiesAlgorithm::Run()
{
    const PfoList *pPfoList(nullptr);
    PANDORA_THROW_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_INITIALIZED, !=, PandoraContentApi::GetList(*this, m_neutrinoPfoListName, pPfoList));

    if (!pPfoList || pPfoList->empty())
    {
        if (PandoraContentApi::GetSettings(*this)->ShouldDisplayAlgorithmInfo())
            std::cout << "NeutrinoPropertiesAlgorithm: unable to find pfo list " << m_neutrinoPfoListName << std::endl;

        return STATUS_CODE_SUCCESS;
    }

    // ATTN Enforces that only one pfo, of neutrino-type, be in the specified input list
    const ParticleFlowObject *const pNeutrinoPfo((1 == pPfoList->size()) ? *(pPfoList->begin()) : nullptr);

    if (!pNeutrinoPfo || !LArPfoHelper::IsNeutrino(pNeutrinoPfo))
        return STATUS_CODE_FAILURE;

    // ATTN At this (maybe unconventional) stage, remove any lone, placeholder neutrinos, with no daughter particles
    if (pNeutrinoPfo->GetDaughterPfoList().empty())
    {
        PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::Delete(*this, pNeutrinoPfo, m_neutrinoPfoListName));
    }
    else
    {
        this->SetNeutrinoId(pNeutrinoPfo);
    }

    return STATUS_CODE_SUCCESS;
}

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

void NeutrinoPropertiesAlgorithm::SetNeutrinoId(const ParticleFlowObject *const pNeutrinoPfo) const
{
    unsigned int nPrimaryTwoDHits(0);
    const ParticleFlowObject *pPrimaryDaughter(nullptr);

    PfoVector daughterPfoVector(pNeutrinoPfo->GetDaughterPfoList().begin(), pNeutrinoPfo->GetDaughterPfoList().end());
    std::sort(daughterPfoVector.begin(), daughterPfoVector.end(), LArPfoHelper::SortByNHits);

    for (const ParticleFlowObject *const pDaughterPfo : daughterPfoVector)
    {
        const unsigned int nTwoDHits(this->GetNTwoDHitsInPfoChain(pDaughterPfo));

        if (!pPrimaryDaughter || (nTwoDHits > nPrimaryTwoDHits))
        {
            nPrimaryTwoDHits = nTwoDHits;
            pPrimaryDaughter = pDaughterPfo;
        }
    }

    if (!pPrimaryDaughter)
        throw StatusCodeException(STATUS_CODE_FAILURE);

    PandoraContentApi::ParticleFlowObject::Metadata metadata;

    if (E_MINUS == std::abs(pPrimaryDaughter->GetParticleId()))
    {
        metadata.m_particleId = NU_E;
    }
    else if (MU_MINUS == std::abs(pPrimaryDaughter->GetParticleId()))
    {
        metadata.m_particleId = NU_MU;
    }

    if (metadata.m_particleId.IsInitialized())
    {
        metadata.m_charge = PdgTable::GetParticleCharge(metadata.m_particleId.Get());
        metadata.m_mass = PdgTable::GetParticleMass(metadata.m_particleId.Get());
        PANDORA_THROW_RESULT_IF(STATUS_CODE_SUCCESS, !=, PandoraContentApi::ParticleFlowObject::AlterMetadata(*this, pNeutrinoPfo, metadata));
    }
}

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

unsigned int NeutrinoPropertiesAlgorithm::GetNTwoDHitsInPfoChain(const ParticleFlowObject *const pPfo) const
{
    unsigned int nTwoDHits(0);

    for (const Cluster *const pCluster : pPfo->GetClusterList())
    {
        const HitType hitType(LArClusterHelper::GetClusterHitType(pCluster));

        if ((TPC_VIEW_U == hitType) || (TPC_VIEW_V == hitType) || (TPC_VIEW_W == hitType))
        {
            nTwoDHits += pCluster->GetNCaloHits();

            if (m_includeIsolatedHits)
                nTwoDHits += pCluster->GetNIsolatedCaloHits();
        }
    }

    for (const ParticleFlowObject *const pDaughterPfo : pPfo->GetDaughterPfoList())
    {
        nTwoDHits += this->GetNTwoDHitsInPfoChain(pDaughterPfo);
    }

    return nTwoDHits;
}

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

StatusCode NeutrinoPropertiesAlgorithm::ReadSettings(const TiXmlHandle xmlHandle)
{
    PANDORA_RETURN_RESULT_IF(STATUS_CODE_SUCCESS, !=, XmlHelper::ReadValue(xmlHandle,
        "NeutrinoPfoListName", m_neutrinoPfoListName));

    PANDORA_RETURN_RESULT_IF_AND_IF(STATUS_CODE_SUCCESS, STATUS_CODE_NOT_FOUND, !=, XmlHelper::ReadValue(xmlHandle,
        "IncludeIsolatedHits", m_includeIsolatedHits));

    return STATUS_CODE_SUCCESS;
}

} // namespace lar_content