PFParticleValidation_module.cc

Go to the documentation of this file.

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

#include "larcoreobj/SimpleTypesAndConstants/geo_types.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