PFParticleValidation_module.cc

Go to the documentation of this file.

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

#include "larpandora/LArPandoraInterface/LArPandoraHelper.h"

#include <string>

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

namespace lar_pandora
{

class PFParticleValidation : public art::EDAnalyzer
{
public:
     PFParticleValidation(fhicl::ParameterSet const &pset);

     virtual ~PFParticleValidation();

     void beginJob();
     void endJob();
     void analyze(const art::Event &evt);
     void reconfigure(fhicl::ParameterSet const &pset);

private:
    class SimpleMCPrimary
    {
    public:
        SimpleMCPrimary();

        bool operator<(const SimpleMCPrimary &rhs) const;

        int                                 m_id;                       
        int                                 m_pdgCode;                  
        int                                 m_nMCHitsTotal;             
        int                                 m_nMCHitsU;                 
        int                                 m_nMCHitsV;                 
        int                                 m_nMCHitsW;                 
        float                               m_energy;                   
        int                                 m_nMatchedPfos;             
        const simb::MCParticle             *m_pAddress;                 
    };

    typedef std::vector<SimpleMCPrimary> SimpleMCPrimaryList;

    class SimpleMatchedPfo
    {
    public:
        SimpleMatchedPfo();

        int                                 m_id;                       
        int                                 m_parentId;                 
        int                                 m_pdgCode;                  
        int                                 m_nPfoHitsTotal;            
        int                                 m_nPfoHitsU;                
        int                                 m_nPfoHitsV;                
        int                                 m_nPfoHitsW;                
        int                                 m_nMatchedHitsTotal;        
        int                                 m_nMatchedHitsU;            
        int                                 m_nMatchedHitsV;            
        int                                 m_nMatchedHitsW;            
        const recob::PFParticle            *m_pAddress;                 
    };

    typedef std::vector<SimpleMatchedPfo> SimpleMatchedPfoList;

    class MatchingDetails
    {
    public:
        MatchingDetails();

        int                                 m_matchedPrimaryId;         
        int                                 m_nMatchedHits;             
        float                               m_completeness;             
    };

    typedef std::map<int, MatchingDetails> MatchingDetailsMap;
    typedef std::map<SimpleMCPrimary, SimpleMatchedPfoList> MCPrimaryMatchingMap; // SimpleMCPrimary has a defined operator<

    typedef std::map< art::Ptr<recob::PFParticle>, HitVector > PFParticleToMatchedHits;
    typedef std::map< art::Ptr<simb::MCParticle>,  PFParticleToMatchedHits > MCParticleMatchingMap;

    void GetMCParticleMatchingMap(const PFParticlesToHits &recoParticlesToHits, const MCParticlesToHits &trueParticlesToHits,
        const HitsToMCParticles &hitsToTrueParticles, MCParticleMatchingMap &mcParticleMatchingMap) const;

    void GetSimpleMCPrimaryList(const art::Event &evt, const MCParticlesToHits &mcParticlesToHits, const HitsToMCParticles &hitsToMCParticles,
        const MCParticleMatchingMap &mcParticleMatchingMap, SimpleMCPrimaryList &simpleMCPrimaryList) const;

    void GetMCPrimaryMatchingMap(const SimpleMCPrimaryList &simpleMCPrimaryList, const MCParticleMatchingMap &mcParticleMatchingMap,
        const PFParticlesToHits &pfParticlesToHits, MCPrimaryMatchingMap &mcPrimaryMatchingMap) const;

    bool IsNeutrinoInduced(const art::Ptr<simb::MCParticle> pMCParticle, const MCParticlesToMCTruth &artMCParticlesToMCTruth) const;

    void GetMCTruth(const art::Event &evt, MCTruthVector &mcTruthVector) const;

    void GetRecoNeutrinos(const art::Event &evt, PFParticleVector &recoNeutrinoVector) const;

    void PrintAllOutput(const MCTruthVector &mcTruthVector, const PFParticleVector &recoNeutrinoVector, const MCPrimaryMatchingMap &mcPrimaryMatchingMap) const;

    void PerformMatching(const MCPrimaryMatchingMap &mcPrimaryMatchingMap, MatchingDetailsMap &matchingDetailsMap) const;

    typedef std::set<int> IntSet;

    bool GetStrongestPfoMatch(const MCPrimaryMatchingMap &mcPrimaryMatchingMap, IntSet &usedMCIds, IntSet &usedPfoIds, MatchingDetailsMap &matchingDetailsMap) const;

    void GetRemainingPfoMatches(const MCPrimaryMatchingMap &mcPrimaryMatchingMap, const IntSet &usedPfoIds, MatchingDetailsMap &matchingDetailsMap) const;

    void PrintMatchingOutput(const MCPrimaryMatchingMap &mcPrimaryMatchingMap, const MatchingDetailsMap &matchingDetailsMap) const;

    bool IsGoodMCPrimary(const SimpleMCPrimary &simpleMCPrimary) const;

    bool HasMatch(const SimpleMCPrimary &simpleMCPrimary, const SimpleMatchedPfoList &simpleMatchedPfoList, const MatchingDetailsMap &matchingDetailsMap) const;

    bool IsGoodMatch(const SimpleMCPrimary &simpleMCPrimary, const SimpleMatchedPfo &simpleMatchedPfo) const;

    unsigned int CountHitsByType(const geo::View_t view, const HitVector &hitVector) const;

    static bool SortSimpleMCPrimaries(const SimpleMCPrimary &lhs, const SimpleMCPrimary &rhs);

    static bool SortSimpleMatchedPfos(const SimpleMatchedPfo &lhs, const SimpleMatchedPfo &rhs);

    std::string         m_hitfinderLabel;               
    std::string         m_particleLabel;                
    std::string         m_geantModuleLabel;             
    std::string         m_backtrackerLabel;             

    bool                m_printAllToScreen;             
    bool                m_printMatchingToScreen;        

    bool                m_neutrinoInducedOnly;          

    int                 m_matchingMinPrimaryHits;       
    int                 m_matchingMinHitsForGoodView;   
    int                 m_matchingMinPrimaryGoodViews;  

    bool                m_useSmallPrimaries;            
    int                 m_matchingMinSharedHits;        
    float               m_matchingMinCompleteness;      
    float               m_matchingMinPurity;            
};

DEFINE_ART_MODULE(PFParticleValidation)

} // namespace lar_pandora

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


#include "art/Framework/Principal/Event.h"

#include "fhiclcpp/ParameterSet.h"

#include "nusimdata/SimulationBase/MCParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"

#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/PFParticle.h"

#include <algorithm>
#include <iostream>

namespace lar_pandora
{

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

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

PFParticleValidation::~PFParticleValidation()
{
}

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

void PFParticleValidation::reconfigure(fhicl::ParameterSet const &pset)
{
    m_particleLabel = pset.get<std::string>("PFParticleModule", "pandoraNu");
    m_hitfinderLabel = pset.get<std::string>("HitFinderModule", "gaushit");
    m_geantModuleLabel = pset.get<std::string>("GeantModule","largeant");
    m_backtrackerLabel = pset.get<std::string>("BackTrackerModule","gaushitTruthMatch");
    m_printAllToScreen = pset.get<bool>("PrintAllToScreen", true);
    m_printMatchingToScreen = pset.get<bool>("PrintMatchingToScreen", true);
    m_neutrinoInducedOnly = pset.get<bool>("NeutrinoInducedOnly", true);
    m_matchingMinPrimaryHits = pset.get<int>("MatchingMinPrimaryHits", 15);
    m_matchingMinHitsForGoodView = pset.get<int>("MatchingMinHitsForGoodView", 5);
    m_matchingMinPrimaryGoodViews = pset.get<int>("MatchingMinPrimaryGoodViews", 2);
    m_useSmallPrimaries = pset.get<bool>("UseSmallPrimaries", true);
    m_matchingMinSharedHits = pset.get<int>("MatchingMinSharedHits", 5);
    m_matchingMinCompleteness = pset.get<float>("MatchingMinCompleteness", 0.1f);
    m_matchingMinPurity = pset.get<float>("MatchingMinPurity", 0.5f);
}

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

void PFParticleValidation::beginJob()
{
}

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

void PFParticleValidation::endJob()
{
}

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

void PFParticleValidation::analyze(const art::Event &evt)
{
    HitVector hitVector;
    LArPandoraHelper::CollectHits(evt, m_hitfinderLabel, hitVector);

    PFParticlesToHits pfParticlesToHits;
    HitsToPFParticles hitsToPfParticles;
    LArPandoraHelper::BuildPFParticleHitMaps(evt, m_particleLabel, pfParticlesToHits, hitsToPfParticles, LArPandoraHelper::kAddDaughters);

    MCParticlesToHits mcParticlesToHits;
    HitsToMCParticles hitsToMCParticles;

    LArPandoraHelper::BuildMCParticleHitMaps(evt, m_geantModuleLabel, hitVector,
        mcParticlesToHits, hitsToMCParticles, LArPandoraHelper::kAddDaughters);

    if (hitsToMCParticles.empty())
    {
        if (m_backtrackerLabel.empty())
            throw cet::exception("LArPandora") << " PFParticleValidation::analyze - no sim channels found, backtracker module must be set in FHiCL " << std::endl;

        LArPandoraHelper::BuildMCParticleHitMaps(evt, m_geantModuleLabel, m_hitfinderLabel, m_backtrackerLabel,
            mcParticlesToHits, hitsToMCParticles, LArPandoraHelper::kAddDaughters);
    }

    MCParticleMatchingMap mcParticleMatchingMap;
    this->GetMCParticleMatchingMap(pfParticlesToHits, mcParticlesToHits, hitsToMCParticles, mcParticleMatchingMap);

    SimpleMCPrimaryList simpleMCPrimaryList;
    this->GetSimpleMCPrimaryList(evt, mcParticlesToHits, hitsToMCParticles, mcParticleMatchingMap, simpleMCPrimaryList);

    MCPrimaryMatchingMap mcPrimaryMatchingMap;
    this->GetMCPrimaryMatchingMap(simpleMCPrimaryList, mcParticleMatchingMap, pfParticlesToHits, mcPrimaryMatchingMap);

    MCTruthVector mcTruthVector;
    this->GetMCTruth(evt, mcTruthVector);

    PFParticleVector recoNeutrinoVector;
    this->GetRecoNeutrinos(evt, recoNeutrinoVector);

    if (m_printAllToScreen)
        this->PrintAllOutput(mcTruthVector, recoNeutrinoVector, mcPrimaryMatchingMap);

    if (m_printMatchingToScreen)
    {
        MatchingDetailsMap matchingDetailsMap;
        this->PerformMatching(mcPrimaryMatchingMap, matchingDetailsMap);
        this->PrintMatchingOutput(mcPrimaryMatchingMap, matchingDetailsMap);
    }
}

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

void PFParticleValidation::GetMCParticleMatchingMap(const PFParticlesToHits &pfParticlesToHits, const MCParticlesToHits &mcParticlesToHits,
    const HitsToMCParticles &hitsToMCParticles, MCParticleMatchingMap &mcParticleMatchingMap) const
{
    // Create a placeholder entry for all mc particles with >0 hits
    for (const MCParticlesToHits::value_type &mcParticleToHitsEntry : mcParticlesToHits)
    {
        if (!mcParticleToHitsEntry.second.empty())
            (void) mcParticleMatchingMap.insert(MCParticleMatchingMap::value_type(mcParticleToHitsEntry.first, PFParticleToMatchedHits()));
    }

    // Store true to reco matching details
    for (const PFParticlesToHits::value_type &recoParticleToHits : pfParticlesToHits)
    {
        const art::Ptr<recob::PFParticle> pRecoParticle(recoParticleToHits.first);
        const HitVector &hitVector(recoParticleToHits.second);

        for (const art::Ptr<recob::Hit> pHit : hitVector)
        {
            HitsToMCParticles::const_iterator mcParticleIter = hitsToMCParticles.find(pHit);

            if (hitsToMCParticles.end() == mcParticleIter)
                continue;

            const art::Ptr<simb::MCParticle> pTrueParticle = mcParticleIter->second;
            mcParticleMatchingMap[pTrueParticle][pRecoParticle].push_back(pHit);
        }
    }
}

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

void PFParticleValidation::GetSimpleMCPrimaryList(const art::Event &evt, const MCParticlesToHits &mcParticlesToHits,
    const HitsToMCParticles &hitsToMCParticles, const MCParticleMatchingMap &mcParticleMatchingMap, SimpleMCPrimaryList &simpleMCPrimaryList) const
{
    MCTruthToMCParticles artMCTruthToMCParticles;
    MCParticlesToMCTruth artMCParticlesToMCTruth;
    LArPandoraHelper::CollectMCParticles(evt, m_geantModuleLabel, artMCTruthToMCParticles, artMCParticlesToMCTruth);

    for (const MCParticlesToHits::value_type &mapEntry : mcParticlesToHits)
    {
        const art::Ptr<simb::MCParticle> pMCPrimary(mapEntry.first);

        if (m_neutrinoInducedOnly && !this->IsNeutrinoInduced(pMCPrimary, artMCParticlesToMCTruth))
            continue;

        SimpleMCPrimary simpleMCPrimary;
        // ATTN simpleMCPrimary.m_id assigned later, after sorting
        simpleMCPrimary.m_pAddress = pMCPrimary.get();
        simpleMCPrimary.m_pdgCode = pMCPrimary->PdgCode();
        simpleMCPrimary.m_energy = pMCPrimary->E();

        MCParticlesToHits::const_iterator trueHitsIter = mcParticlesToHits.find(pMCPrimary);

        if (mcParticlesToHits.end() != trueHitsIter)
        {
            const HitVector &hitVector(trueHitsIter->second);
            simpleMCPrimary.m_nMCHitsTotal = hitVector.size();
            simpleMCPrimary.m_nMCHitsU = this->CountHitsByType(geo::kU, hitVector);
            simpleMCPrimary.m_nMCHitsV = this->CountHitsByType(geo::kV, hitVector);
            simpleMCPrimary.m_nMCHitsW = this->CountHitsByType(geo::kW, hitVector);
        }

        MCParticleMatchingMap::const_iterator matchedPfoIter = mcParticleMatchingMap.find(pMCPrimary);

        if (mcParticleMatchingMap.end() != matchedPfoIter)
            simpleMCPrimary.m_nMatchedPfos = matchedPfoIter->second.size();

        simpleMCPrimaryList.push_back(simpleMCPrimary);
    }

    std::sort(simpleMCPrimaryList.begin(), simpleMCPrimaryList.end(), PFParticleValidation::SortSimpleMCPrimaries);

    int mcPrimaryId(0);
    for (SimpleMCPrimary &simpleMCPrimary : simpleMCPrimaryList)
        simpleMCPrimary.m_id = mcPrimaryId++;
}

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

bool PFParticleValidation::IsNeutrinoInduced(const art::Ptr<simb::MCParticle> pMCParticle, const MCParticlesToMCTruth &artMCParticlesToMCTruth) const
{
    MCParticlesToMCTruth::const_iterator iter = artMCParticlesToMCTruth.find(pMCParticle);

    if (artMCParticlesToMCTruth.end() == iter)
        return false;

    const art::Ptr<simb::MCTruth> pMCTruth = iter->second;
    return (simb::kBeamNeutrino == pMCTruth->Origin());
}

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

void PFParticleValidation::GetMCPrimaryMatchingMap(const SimpleMCPrimaryList &simpleMCPrimaryList, const MCParticleMatchingMap &mcParticleMatchingMap,
    const PFParticlesToHits &pfParticlesToHits, MCPrimaryMatchingMap &mcPrimaryMatchingMap) const
{
    for (const SimpleMCPrimary &simpleMCPrimary : simpleMCPrimaryList)
    {
        SimpleMatchedPfoList simpleMatchedPfoList;
        MCParticleMatchingMap::const_iterator matchedPfoIter = mcParticleMatchingMap.end();

        // ATTN Nasty workaround I
        for (MCParticleMatchingMap::const_iterator iter = mcParticleMatchingMap.begin(), iterEnd = mcParticleMatchingMap.end(); iter != iterEnd; ++iter)
        {
            if (simpleMCPrimary.m_pAddress == iter->first.get())
            {
                matchedPfoIter = iter;
                break;
            };
        }

        if (mcParticleMatchingMap.end() != matchedPfoIter)
        {
            for (const PFParticleToMatchedHits::value_type &contribution : matchedPfoIter->second)
            {
                const art::Ptr<recob::PFParticle> pMatchedPfo(contribution.first);
                const HitVector &matchedHitVector(contribution.second);

                SimpleMatchedPfo simpleMatchedPfo;
                simpleMatchedPfo.m_pAddress = pMatchedPfo.get();
                simpleMatchedPfo.m_id = pMatchedPfo->Self();

                // ATTN Assume pfos have either zero or one parents. Ignore parent neutrino.
                PFParticlesToHits::const_iterator parentPfoIter = pfParticlesToHits.end();

                // ATTN Nasty workaround II, bad place for another loop.
                for (PFParticlesToHits::const_iterator iter = pfParticlesToHits.begin(), iterEnd = pfParticlesToHits.end(); iter != iterEnd; ++iter)
                {
                    if (pMatchedPfo->Parent() == iter->first->Self())
                    {
                        parentPfoIter = iter;
                        break;
                    };
                }

                if ((pfParticlesToHits.end() != parentPfoIter) && !LArPandoraHelper::IsNeutrino(parentPfoIter->first))
                    simpleMatchedPfo.m_parentId = parentPfoIter->first->Self();

                simpleMatchedPfo.m_pdgCode = pMatchedPfo->PdgCode();
                simpleMatchedPfo.m_nMatchedHitsTotal = matchedHitVector.size();
                simpleMatchedPfo.m_nMatchedHitsU = this->CountHitsByType(geo::kU, matchedHitVector);
                simpleMatchedPfo.m_nMatchedHitsV = this->CountHitsByType(geo::kV, matchedHitVector);
                simpleMatchedPfo.m_nMatchedHitsW = this->CountHitsByType(geo::kW, matchedHitVector);

                PFParticlesToHits::const_iterator pfoHitsIter = pfParticlesToHits.find(pMatchedPfo);

                if (pfParticlesToHits.end() == pfoHitsIter)
                    throw cet::exception("LArPandora") << " PFParticleValidation::analyze --- Presence of PFParticle in map mandatory.";

                const HitVector &pfoHitVector(pfoHitsIter->second);

                simpleMatchedPfo.m_nPfoHitsTotal = pfoHitVector.size();
                simpleMatchedPfo.m_nPfoHitsU = this->CountHitsByType(geo::kU, pfoHitVector);
                simpleMatchedPfo.m_nPfoHitsV = this->CountHitsByType(geo::kV, pfoHitVector);
                simpleMatchedPfo.m_nPfoHitsW = this->CountHitsByType(geo::kW, pfoHitVector);

                simpleMatchedPfoList.push_back(simpleMatchedPfo);
            }
        }

        // Store the ordered vectors of matched pfo details
        std::sort(simpleMatchedPfoList.begin(), simpleMatchedPfoList.end(), PFParticleValidation::SortSimpleMatchedPfos);

        if (!mcPrimaryMatchingMap.insert(MCPrimaryMatchingMap::value_type(simpleMCPrimary, simpleMatchedPfoList)).second)
            throw cet::exception("LArPandora") << " PFParticleValidation::analyze --- Double-counting MC primaries.";
    }
}

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

void PFParticleValidation::GetMCTruth(const art::Event &evt, MCTruthVector &mcTruthVector) const
{
    MCTruthToMCParticles artMCTruthToMCParticles;
    MCParticlesToMCTruth artMCParticlesToMCTruth;
    LArPandoraHelper::CollectMCParticles(evt, m_geantModuleLabel, artMCTruthToMCParticles, artMCParticlesToMCTruth);

    for (const auto &mapEntry : artMCTruthToMCParticles)
    {
        const art::Ptr<simb::MCTruth> truth = mapEntry.first;

        if (!truth->NeutrinoSet())
            continue;

        if (mcTruthVector.end() == std::find(mcTruthVector.begin(), mcTruthVector.end(), truth))
            mcTruthVector.push_back(truth);
    }
}

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

void PFParticleValidation::GetRecoNeutrinos(const art::Event &evt, PFParticleVector &recoNeutrinoVector) const
{
    PFParticleVector allPFParticles;
    LArPandoraHelper::CollectPFParticles(evt, m_particleLabel, allPFParticles);
    LArPandoraHelper::SelectNeutrinoPFParticles(allPFParticles, recoNeutrinoVector);
}

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

void PFParticleValidation::PrintAllOutput(const MCTruthVector &mcTruthVector, const PFParticleVector &recoNeutrinoVector,
    const MCPrimaryMatchingMap &mcPrimaryMatchingMap) const
{
    std::cout << "---RAW-MATCHING-OUTPUT--------------------------------------------------------------------------" << std::endl;

    for (const art::Ptr<simb::MCTruth> pMCTruth : mcTruthVector)
    {
        std::cout << "MCNeutrino, PDG " << pMCTruth->GetNeutrino().Nu().PdgCode() << ", InteractionType " << pMCTruth->GetNeutrino().InteractionType() << std::endl;
    }

    for (const art::Ptr<recob::PFParticle> pPfo : recoNeutrinoVector)
    {
        std::cout << "RecoNeutrino, PDG " << pPfo->PdgCode() << ", nDaughters " << pPfo->NumDaughters() << std::endl;
    }

    for (const MCPrimaryMatchingMap::value_type &mapValue : mcPrimaryMatchingMap)
    {
        const SimpleMCPrimary &simpleMCPrimary(mapValue.first);

        std::cout << std::endl << "Primary " << simpleMCPrimary.m_id << ", PDG " << simpleMCPrimary.m_pdgCode << ", nMCHits " << simpleMCPrimary.m_nMCHitsTotal
            << " (" << simpleMCPrimary.m_nMCHitsU << ", " << simpleMCPrimary.m_nMCHitsV << ", " << simpleMCPrimary.m_nMCHitsW << ")" << std::endl;

        for (const SimpleMatchedPfo &simpleMatchedPfo : mapValue.second)
        {
            std::cout << "-MatchedPfo " << simpleMatchedPfo.m_id;

            if (simpleMatchedPfo.m_parentId >= 0)
                std::cout << ", ParentPfo " << simpleMatchedPfo.m_parentId;

            std::cout << ", PDG " << simpleMatchedPfo.m_pdgCode << ", nMatchedHits " << simpleMatchedPfo.m_nMatchedHitsTotal
                << " (" << simpleMatchedPfo.m_nMatchedHitsU << ", " << simpleMatchedPfo.m_nMatchedHitsV << ", " << simpleMatchedPfo.m_nMatchedHitsW << ")"
                << ", nPfoHits " << simpleMatchedPfo.m_nPfoHitsTotal << " (" << simpleMatchedPfo.m_nPfoHitsU << ", " << simpleMatchedPfo.m_nPfoHitsV << ", "
                << simpleMatchedPfo.m_nPfoHitsW << ")" << std::endl;
        }
    }

    std::cout << "------------------------------------------------------------------------------------------------" << std::endl;
}

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

void PFParticleValidation::PerformMatching(const MCPrimaryMatchingMap &mcPrimaryMatchingMap, MatchingDetailsMap &matchingDetailsMap) const
{
    // Get best matches, one-by-one, until no more strong matches possible
    IntSet usedMCIds, usedPfoIds;
    while (GetStrongestPfoMatch(mcPrimaryMatchingMap, usedMCIds, usedPfoIds, matchingDetailsMap)) {}

    // Assign any remaining pfos to primaries, based on number of matched hits
    GetRemainingPfoMatches(mcPrimaryMatchingMap, usedPfoIds, matchingDetailsMap);
}

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

bool PFParticleValidation::GetStrongestPfoMatch(const MCPrimaryMatchingMap &mcPrimaryMatchingMap, IntSet &usedMCIds, IntSet &usedPfoIds,
    MatchingDetailsMap &matchingDetailsMap) const
{
    int bestPfoMatchId(-1);
    MatchingDetails bestMatchingDetails;

    for (const MCPrimaryMatchingMap::value_type &mapValue : mcPrimaryMatchingMap)
    {
        const SimpleMCPrimary &simpleMCPrimary(mapValue.first);

        if (!m_useSmallPrimaries && !this->IsGoodMCPrimary(simpleMCPrimary))
            continue;

        if (usedMCIds.count(simpleMCPrimary.m_id))
            continue;

        for (const SimpleMatchedPfo &simpleMatchedPfo : mapValue.second)
        {
            if (usedPfoIds.count(simpleMatchedPfo.m_id))
                continue;

            if (!this->IsGoodMatch(simpleMCPrimary, simpleMatchedPfo))
                continue;

            if (simpleMatchedPfo.m_nMatchedHitsTotal > bestMatchingDetails.m_nMatchedHits)
            {
                bestPfoMatchId = simpleMatchedPfo.m_id;
                bestMatchingDetails.m_matchedPrimaryId = simpleMCPrimary.m_id;
                bestMatchingDetails.m_nMatchedHits = simpleMatchedPfo.m_nMatchedHitsTotal;
                bestMatchingDetails.m_completeness = static_cast<float>(simpleMatchedPfo.m_nMatchedHitsTotal) / static_cast<float>(simpleMCPrimary.m_nMCHitsTotal);
            }
        }
    }

    if (bestPfoMatchId > -1)
    {
        matchingDetailsMap[bestPfoMatchId] = bestMatchingDetails;
        usedMCIds.insert(bestMatchingDetails.m_matchedPrimaryId);
        usedPfoIds.insert(bestPfoMatchId);
        return true;
    }

    return false;
}

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

void PFParticleValidation::GetRemainingPfoMatches(const MCPrimaryMatchingMap &mcPrimaryMatchingMap, const IntSet &usedPfoIds,
    MatchingDetailsMap &matchingDetailsMap) const
{
    for (const MCPrimaryMatchingMap::value_type &mapValue : mcPrimaryMatchingMap)
    {
        const SimpleMCPrimary &simpleMCPrimary(mapValue.first);

        if (!m_useSmallPrimaries && !this->IsGoodMCPrimary(simpleMCPrimary))
            continue;

        for (const SimpleMatchedPfo &simpleMatchedPfo : mapValue.second)
        {
            if (usedPfoIds.count(simpleMatchedPfo.m_id))
                continue;

            MatchingDetails &matchingDetails(matchingDetailsMap[simpleMatchedPfo.m_id]);

            if (simpleMatchedPfo.m_nMatchedHitsTotal > matchingDetails.m_nMatchedHits)
            {
                matchingDetails.m_matchedPrimaryId = simpleMCPrimary.m_id;
                matchingDetails.m_nMatchedHits = simpleMatchedPfo.m_nMatchedHitsTotal;
                matchingDetails.m_completeness = static_cast<float>(simpleMatchedPfo.m_nMatchedHitsTotal) / static_cast<float>(simpleMCPrimary.m_nMCHitsTotal);
            }
        }
    }
}

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

void PFParticleValidation::PrintMatchingOutput(const MCPrimaryMatchingMap &mcPrimaryMatchingMap, const MatchingDetailsMap &matchingDetailsMap) const
{
    std::cout << "---PROCESSED-MATCHING-OUTPUT--------------------------------------------------------------------" << std::endl;
    std::cout << "MinPrimaryGoodHits " << m_matchingMinPrimaryHits << ", MinHitsForGoodView " << m_matchingMinHitsForGoodView << ", MinPrimaryGoodViews " << m_matchingMinPrimaryGoodViews << std::endl;
    std::cout << "UseSmallPrimaries " << m_useSmallPrimaries << ", MinSharedHits " << m_matchingMinSharedHits << ", MinCompleteness " << m_matchingMinCompleteness << ", MinPurity " << m_matchingMinPurity << std::endl;

    bool isCorrect(true), isCalculable(false);

    for (const MCPrimaryMatchingMap::value_type &mapValue : mcPrimaryMatchingMap)
    {
        const SimpleMCPrimary &simpleMCPrimary(mapValue.first);
        const bool hasMatch(this->HasMatch(simpleMCPrimary, mapValue.second, matchingDetailsMap));
        const bool isTargetPrimary(this->IsGoodMCPrimary(simpleMCPrimary) && (2112 != simpleMCPrimary.m_pdgCode));

        if (!hasMatch && !isTargetPrimary)
            continue;

        std::cout << std::endl << (!isTargetPrimary ? "(Non target) " : "")
                  << "Primary " << simpleMCPrimary.m_id << ", PDG " << simpleMCPrimary.m_pdgCode << ", nMCHits " << simpleMCPrimary.m_nMCHitsTotal
                  << " (" << simpleMCPrimary.m_nMCHitsU << ", " << simpleMCPrimary.m_nMCHitsV << ", " << simpleMCPrimary.m_nMCHitsW << ")" << std::endl;

        if (2112 != simpleMCPrimary.m_pdgCode)
            isCalculable = true;

        unsigned int nMatches(0);

        for (const SimpleMatchedPfo &simpleMatchedPfo : mapValue.second)
        {
            if (matchingDetailsMap.count(simpleMatchedPfo.m_id) && (simpleMCPrimary.m_id == matchingDetailsMap.at(simpleMatchedPfo.m_id).m_matchedPrimaryId))
            {
                const bool isGoodMatch(this->IsGoodMatch(simpleMCPrimary, simpleMatchedPfo));

                if (isGoodMatch) ++nMatches;
                std::cout << "-" << (!isGoodMatch ? "(Below threshold) " : "") << "MatchedPfo " << simpleMatchedPfo.m_id;

                if (simpleMatchedPfo.m_parentId >= 0) std::cout << ", ParentPfo " << simpleMatchedPfo.m_parentId;

                std::cout << ", PDG " << simpleMatchedPfo.m_pdgCode << ", nMatchedHits " << simpleMatchedPfo.m_nMatchedHitsTotal
                          << " (" << simpleMatchedPfo.m_nMatchedHitsU << ", " << simpleMatchedPfo.m_nMatchedHitsV << ", " << simpleMatchedPfo.m_nMatchedHitsW << ")"
                          << ", nPfoHits " << simpleMatchedPfo.m_nPfoHitsTotal
                          << " (" << simpleMatchedPfo.m_nPfoHitsU << ", " << simpleMatchedPfo.m_nPfoHitsV << ", " << simpleMatchedPfo.m_nPfoHitsW << ")" << std::endl;
            }
        }

        if (isTargetPrimary && (1 != nMatches))
            isCorrect = false;
    }

    std::cout << std::endl << "Is correct? " << (isCorrect && isCalculable) << std::endl;
    std::cout << "------------------------------------------------------------------------------------------------" << std::endl;
}

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

bool PFParticleValidation::IsGoodMCPrimary(const SimpleMCPrimary &simpleMCPrimary) const
{
    if (simpleMCPrimary.m_nMCHitsTotal < m_matchingMinPrimaryHits)
        return false;

    int nGoodViews(0);
    if (simpleMCPrimary.m_nMCHitsU >= m_matchingMinHitsForGoodView) ++nGoodViews;
    if (simpleMCPrimary.m_nMCHitsV >= m_matchingMinHitsForGoodView) ++nGoodViews;
    if (simpleMCPrimary.m_nMCHitsW >= m_matchingMinHitsForGoodView) ++nGoodViews;

    if (nGoodViews < m_matchingMinPrimaryGoodViews)
        return false;

    return true;
}

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

bool PFParticleValidation::HasMatch(const SimpleMCPrimary &simpleMCPrimary, const SimpleMatchedPfoList &simpleMatchedPfoList,
    const MatchingDetailsMap &matchingDetailsMap) const
{
    for (const SimpleMatchedPfo &simpleMatchedPfo : simpleMatchedPfoList)
    {
        if (matchingDetailsMap.count(simpleMatchedPfo.m_id) && (simpleMCPrimary.m_id == matchingDetailsMap.at(simpleMatchedPfo.m_id).m_matchedPrimaryId))
            return true;
    }

    return false;
}

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

bool PFParticleValidation::IsGoodMatch(const SimpleMCPrimary &simpleMCPrimary, const SimpleMatchedPfo &simpleMatchedPfo) const
{
    const float purity((simpleMatchedPfo.m_nPfoHitsTotal > 0) ? static_cast<float>(simpleMatchedPfo.m_nMatchedHitsTotal) / static_cast<float>(simpleMatchedPfo.m_nPfoHitsTotal) : 0.f);
    const float completeness((simpleMCPrimary.m_nMCHitsTotal > 0) ? static_cast<float>(simpleMatchedPfo.m_nMatchedHitsTotal) / static_cast<float>(simpleMCPrimary.m_nMCHitsTotal) : 0.f);

    return ((simpleMatchedPfo.m_nMatchedHitsTotal >= m_matchingMinSharedHits) && (purity >= m_matchingMinPurity) && (completeness >= m_matchingMinCompleteness));
}

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

unsigned int PFParticleValidation::CountHitsByType(const geo::View_t view, const HitVector &hitVector) const
{
    unsigned int nHitsOfSpecifiedType(0);

    for (const art::Ptr<recob::Hit> pHit : hitVector)
    {
        if (view == pHit->View())
            ++nHitsOfSpecifiedType;
    }

    return nHitsOfSpecifiedType;
}

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

bool PFParticleValidation::SortSimpleMCPrimaries(const SimpleMCPrimary &lhs, const SimpleMCPrimary &rhs)
{
    if (lhs.m_nMCHitsTotal != rhs.m_nMCHitsTotal)
        return (lhs.m_nMCHitsTotal > rhs.m_nMCHitsTotal);

    return (lhs.m_energy > rhs.m_energy);
}

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

bool PFParticleValidation::SortSimpleMatchedPfos(const SimpleMatchedPfo &lhs, const SimpleMatchedPfo &rhs)
{
    if (lhs.m_nMatchedHitsTotal != rhs.m_nMatchedHitsTotal)
        return (lhs.m_nMatchedHitsTotal > rhs.m_nMatchedHitsTotal);

    if (lhs.m_nPfoHitsTotal != rhs.m_nPfoHitsTotal)
        return (lhs.m_nPfoHitsTotal > rhs.m_nPfoHitsTotal);

    return (lhs.m_id < rhs.m_id);
}

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

PFParticleValidation::SimpleMCPrimary::SimpleMCPrimary() :
    m_id(-1),
    m_pdgCode(0),
    m_nMCHitsTotal(0),
    m_nMCHitsU(0),
    m_nMCHitsV(0),
    m_nMCHitsW(0),
    m_energy(0.f),
    m_nMatchedPfos(0),
    m_pAddress(nullptr)
{
}

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

bool PFParticleValidation::SimpleMCPrimary::operator<(const SimpleMCPrimary &rhs) const
{
    if (this == &rhs)
        return false;

    return (m_id < rhs.m_id);
}

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

PFParticleValidation::SimpleMatchedPfo::SimpleMatchedPfo() :
    m_id(-1),
    m_parentId(-1),
    m_pdgCode(0),
    m_nPfoHitsTotal(0),
    m_nPfoHitsU(0),
    m_nPfoHitsV(0),
    m_nPfoHitsW(0),
    m_nMatchedHitsTotal(0),
    m_nMatchedHitsU(0),
    m_nMatchedHitsV(0),
    m_nMatchedHitsW(0),
    m_pAddress(nullptr)
{
}

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

PFParticleValidation::MatchingDetails::MatchingDetails() :
    m_matchedPrimaryId(-1),
    m_nMatchedHits(0),
    m_completeness(0.f)
{
}

} //namespace lar_pandora