ConsolidatedPFParticleAnalysisTemplate_module.cc

Go to the documentation of this file.

#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Core/EDAnalyzer.h"

#include "lardataobj/RecoBase/PFParticle.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/Shower.h"

#include <string>

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

namespace lar_pandora
{

class ConsolidatedPFParticleAnalysisTemplate : public art::EDAnalyzer
{
public:
    typedef art::Handle< std::vector<recob::PFParticle> > PFParticleHandle;
    typedef std::map< size_t, art::Ptr<recob::PFParticle> > PFParticleIdMap;
    typedef std::vector< art::Ptr<recob::PFParticle> > PFParticleVector;
    typedef std::vector< art::Ptr<recob::Track> > TrackVector;
    typedef std::vector< art::Ptr<recob::Shower> > ShowerVector;

    ConsolidatedPFParticleAnalysisTemplate(fhicl::ParameterSet const &pset);
    
    void reconfigure(fhicl::ParameterSet const &pset);

    void analyze(const art::Event &evt);

private:
    void GetPFParticleIdMap(const PFParticleHandle &pfParticleHandle, PFParticleIdMap &pfParticleMap);

    void PrintOutScores(const art::Event &evt, const PFParticleHandle &pfParticleHandle) const;

    void GetFinalStatePFParticleVectors(const PFParticleIdMap &pfParticleMap, PFParticleVector &crParticles, PFParticleVector &nuParticles);

    void CollectTracksAndShowers(const PFParticleVector &particles, const PFParticleHandle &pfParticleHandle, const art::Event &evt, TrackVector &tracks, ShowerVector &showers);

    std::string m_pandoraLabel;         
    std::string m_trackLabel;           
    std::string m_showerLabel;          
    bool        m_printOutScores;       
};

DEFINE_ART_MODULE(ConsolidatedPFParticleAnalysisTemplate)

} // namespace lar_pandora

//------------------------------------------------------------------------------------------------------------------------------------------
// implementation follows

#include "art/Framework/Principal/Event.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "lardataobj/RecoBase/PFParticleMetadata.h"

#include "Pandora/PdgTable.h"

#include <iostream>

namespace lar_pandora
{

ConsolidatedPFParticleAnalysisTemplate::ConsolidatedPFParticleAnalysisTemplate(fhicl::ParameterSet const &pset) : art::EDAnalyzer(pset)
{
    this->reconfigure(pset);
}

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

void ConsolidatedPFParticleAnalysisTemplate::reconfigure(fhicl::ParameterSet const &pset)
{
    m_pandoraLabel = pset.get<std::string>("PandoraLabel");
    m_trackLabel = pset.get<std::string>("TrackLabel");
    m_showerLabel = pset.get<std::string>("ShowerLabel");
    m_printOutScores = pset.get<bool>("PrintOutScores",true);
}

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

void ConsolidatedPFParticleAnalysisTemplate::analyze(const art::Event &evt)
{
    // Collect the PFParticles from the event
    PFParticleHandle pfParticleHandle;
    evt.getByLabel(m_pandoraLabel, pfParticleHandle);
   
    if (!pfParticleHandle.isValid())
    {
        mf::LogDebug("ConsolidatedPFParticleAnalysisTemplate") << "  Failed to find the PFParticles." << std::endl;
        return;
    }

    // Produce a map of the PFParticle IDs for fast navigation through the hierarchy
    PFParticleIdMap pfParticleMap;
    this->GetPFParticleIdMap(pfParticleHandle, pfParticleMap);
    
    if (m_printOutScores)
        this->PrintOutScores(evt, pfParticleHandle);

    // Produce two PFParticle vectors containing final-state particles:
    // 1. Particles identified as cosmic-rays - recontructed under cosmic-hypothesis
    // 2. Daughters of the neutrino PFParticle - reconstructed under the neutrino hypothesis
    std::vector< art::Ptr<recob::PFParticle> > crParticles;
    std::vector< art::Ptr<recob::PFParticle> > nuParticles;
    this->GetFinalStatePFParticleVectors(pfParticleMap, crParticles, nuParticles);

    // Use as required!
    // -----------------------------
    //   What follows is an example showing how one might access the reconstructed neutrino final-state tracks and showers
    
    // These are the vectors to hold the tracks and showers for the final-states of the reconstructed neutrino
    std::vector< art::Ptr<recob::Track> > tracks;
    std::vector< art::Ptr<recob::Shower> > showers;
    this->CollectTracksAndShowers(nuParticles, pfParticleHandle, evt, tracks, showers);

    // Print a summary of the consolidated event
    std::cout << "Consolidated event summary:" << std::endl;
    std::cout << "  - Number of primary cosmic-ray PFParticles   : " << crParticles.size() << std::endl;
    std::cout << "  - Number of neutrino final-state PFParticles : " << nuParticles.size() << std::endl;
    std::cout << "    ... of which are track-like   : " << tracks.size() << std::endl;
    std::cout << "    ... of which are showers-like : " << showers.size() << std::endl;
}

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

void ConsolidatedPFParticleAnalysisTemplate::GetPFParticleIdMap(const PFParticleHandle &pfParticleHandle, PFParticleIdMap &pfParticleMap)
{
    for (unsigned int i = 0; i < pfParticleHandle->size(); ++i)
    {
        const art::Ptr<recob::PFParticle> pParticle(pfParticleHandle, i);
        if (!pfParticleMap.insert(PFParticleIdMap::value_type(pParticle->Self(), pParticle)).second)
        {
            throw cet::exception("ConsolidatedPFParticleAnalysisTemplate") << "  Unable to get PFParticle ID map, the input PFParticle collection has repeat IDs!";
        }
    }
}

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

void ConsolidatedPFParticleAnalysisTemplate::PrintOutScores(const art::Event &evt, const PFParticleHandle &pfParticleHandle) const
{
    // Get the associations between PFParticles and larpandoraobj::PFParticleMetadata
    art::FindManyP< larpandoraobj::PFParticleMetadata > pfPartToMetadataAssoc(pfParticleHandle, evt, m_pandoraLabel);

    for (unsigned int i = 0; i < pfParticleHandle->size(); ++i)
    {
        const std::vector< art::Ptr<larpandoraobj::PFParticleMetadata> > &pfParticleMetadataList(pfPartToMetadataAssoc.at(i));
        if (!pfParticleMetadataList.empty())
        {
            const art::Ptr<recob::PFParticle> pParticle(pfParticleHandle, i);
            for (unsigned int j=0; j<pfParticleMetadataList.size(); ++j)
            {
                const art::Ptr<larpandoraobj::PFParticleMetadata> &pfParticleMetadata(pfParticleMetadataList.at(j));
                const larpandoraobj::PFParticleMetadata::PropertiesMap &pfParticlePropertiesMap(pfParticleMetadata->GetPropertiesMap());
                if (!pfParticlePropertiesMap.empty())
                    std::cout << " Found PFParticle " << pParticle->Self() << " with: " << std::endl;
                for (larpandoraobj::PFParticleMetadata::PropertiesMap::const_iterator it = pfParticlePropertiesMap.begin(); it != pfParticlePropertiesMap.end(); ++it)
                    std::cout << "  - " << it->first << " = " << it->second << std::endl;
            }
        }
    }
}

//------------------------------------------------------------------------------------------------------------------------------------------
    
void ConsolidatedPFParticleAnalysisTemplate::GetFinalStatePFParticleVectors(const PFParticleIdMap &pfParticleMap, PFParticleVector &crParticles, PFParticleVector &nuParticles)
{
    for (PFParticleIdMap::const_iterator it = pfParticleMap.begin(); it != pfParticleMap.end(); ++it)
    {
        const art::Ptr<recob::PFParticle> pParticle(it->second);

        // Only look for primary particles
        if (!pParticle->IsPrimary()) continue;

        // Check if this particle is identified as the neutrino
        const int pdg(pParticle->PdgCode());
        const bool isNeutrino(std::abs(pdg) == pandora::NU_E || std::abs(pdg) == pandora::NU_MU || std::abs(pdg) == pandora::NU_TAU);

        // All non-neutrino primary particles are reconstructed under the cosmic hypothesis
        if (!isNeutrino)
        {
            crParticles.push_back(pParticle);
            continue;
        }

        // ATTN. We are filling nuParticles under the assumption that there is only one reconstructed neutrino identified per event.
        //       If this is not the case please handle accordingly
        if (!nuParticles.empty())
        {
            throw cet::exception("ConsolidatedPFParticleAnalysisTemplate") << "  This event contains multiple reconstructed neutrinos!";
        }

        // Add the daughters of the neutrino PFParticle to the nuPFParticles vector
        for (const size_t daughterId : pParticle->Daughters())
        {
            if (pfParticleMap.find(daughterId) == pfParticleMap.end())
                throw cet::exception("ConsolidatedPFParticleAnalysisTemplate") << "  Invalid PFParticle collection!";

            nuParticles.push_back(pfParticleMap.at(daughterId));
        }
    }
}

//------------------------------------------------------------------------------------------------------------------------------------------
    
void ConsolidatedPFParticleAnalysisTemplate::CollectTracksAndShowers(const PFParticleVector &particles, const PFParticleHandle &pfParticleHandle, const art::Event &evt, TrackVector &tracks, ShowerVector &showers)
{
    // Get the associations between PFParticles and tracks/showers from the event
    art::FindManyP< recob::Track > pfPartToTrackAssoc(pfParticleHandle, evt, m_trackLabel);
    art::FindManyP< recob::Shower > pfPartToShowerAssoc(pfParticleHandle, evt, m_showerLabel);
   
    for (const art::Ptr<recob::PFParticle> &pParticle : particles)
    {
        const std::vector< art::Ptr<recob::Track> > associatedTracks(pfPartToTrackAssoc.at(pParticle.key()));
        const std::vector< art::Ptr<recob::Shower> > associatedShowers(pfPartToShowerAssoc.at(pParticle.key()));
        const unsigned int nTracks(associatedTracks.size());
        const unsigned int nShowers(associatedShowers.size());

        // Check if the PFParticle has no associated tracks or showers
        if (nTracks == 0 && nShowers == 0)
        {
            mf::LogDebug("ConsolidatedPFParticleAnalysisTemplate") << "  No tracks or showers were associated to PFParticle " << pParticle->Self() << std::endl;
            continue;
        }

        // Check if there is an associated track
        if (nTracks == 1 && nShowers == 0)
        {
            tracks.push_back(associatedTracks.front());
            continue;
        }

        // Check if there is an associated shower
        if (nTracks == 0 && nShowers == 1)
        {
            showers.push_back(associatedShowers.front());
            continue;
        }

        throw cet::exception("ConsolidatedPFParticleAnalysisTemplate") << "  There were " << nTracks << " tracks and " << nShowers << " showers associated with PFParticle " << pParticle->Self();
    }
}

} //namespace lar_pandora