LArPandoraSliceIdHelper.cxx

Go to the documentation of this file.

#include "larpandora/LArPandoraEventBuilding/LArPandoraSliceIdHelper.h"

#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/PFParticle.h"
#include "nusimdata/SimulationBase/MCTruth.h"

#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/FindManyP.h"

namespace lar_pandora {

  void LArPandoraSliceIdHelper::GetSliceMetadata(const SliceVector& slices,
                                                 const art::Event& evt,
                                                 const std::string& truthLabel,
                                                 const std::string& mcParticleLabel,
                                                 const std::string& hitLabel,
                                                 const std::string& backtrackLabel,
                                                 const std::string& pandoraLabel,
                                                 SliceMetadataVector& sliceMetadata,
                                                 simb::MCNeutrino& mcNeutrino)
  {
    // Find the beam neutrino in MC
    const auto beamNuMCTruth(LArPandoraSliceIdHelper::GetBeamNeutrinoMCTruth(evt, truthLabel));
    mcNeutrino = beamNuMCTruth->GetNeutrino();

    // Collect all MC particles resulting from the beam neutrino
    MCParticleVector mcParticles;
    LArPandoraSliceIdHelper::CollectNeutrinoMCParticles(
      evt, truthLabel, mcParticleLabel, beamNuMCTruth, mcParticles);

    // Get the hits and determine which are neutrino induced
    HitVector hits;
    HitToBoolMap hitToIsNuInducedMap;
    LArPandoraSliceIdHelper::GetHitOrigins(
      evt, hitLabel, backtrackLabel, mcParticles, hits, hitToIsNuInducedMap);
    const unsigned int nNuHits(
      LArPandoraSliceIdHelper::CountNeutrinoHits(hits, hitToIsNuInducedMap));

    // Get the mapping from PFParticle to hits through clusters
    PFParticlesToHits pfParticleToHitsMap;
    LArPandoraSliceIdHelper::GetPFParticleToHitsMap(evt, pandoraLabel, pfParticleToHitsMap);

    // Calculate the metadata for each slice
    LArPandoraSliceIdHelper::GetSliceMetadata(
      slices, pfParticleToHitsMap, hitToIsNuInducedMap, nNuHits, sliceMetadata);
  }

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

  art::Ptr<simb::MCTruth> LArPandoraSliceIdHelper::GetBeamNeutrinoMCTruth(
    const art::Event& evt,
    const std::string& truthLabel)
  {
    // Get the MCTruth handle
    art::Handle<std::vector<simb::MCTruth>> mcTruthHandle;
    evt.getByLabel(truthLabel, mcTruthHandle);

    if (!mcTruthHandle.isValid())
      throw cet::exception("LArPandora")
        << " LArPandoraSliceIdHelper::GetBeamNeutrinoMCTruth - invalid MCTruth handle" << std::endl;

    // Look for the truth block that is from the beam neutrino, and ensure we choose the one with the highest energy if there are multiple
    bool foundNeutrino(false);
    float maxNeutrinoEnergy(-std::numeric_limits<float>::max());
    art::Ptr<simb::MCTruth> beamNuMCTruth;
    for (unsigned int i = 0; i < mcTruthHandle->size(); ++i) {
      const art::Ptr<simb::MCTruth> mcTruth(mcTruthHandle, i);

      if (mcTruth->Origin() != simb::kBeamNeutrino) continue;

      const float nuEnergy(mcTruth->GetNeutrino().Nu().E());
      if (nuEnergy < maxNeutrinoEnergy) continue;

      // Select this as the beam neutrino
      maxNeutrinoEnergy = nuEnergy;
      beamNuMCTruth = mcTruth;
      foundNeutrino = true;
    }

    if (!foundNeutrino)
      throw cet::exception("LArPandora") << " LArPandoraSliceIdHelper::GetBeamNeutrinoMCTruth - "
                                            "found no beam neutrino MCTruth blocks"
                                         << std::endl;

    return beamNuMCTruth;
  }

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

  void LArPandoraSliceIdHelper::CollectNeutrinoMCParticles(
    const art::Event& evt,
    const std::string& truthLabel,
    const std::string& mcParticleLabel,
    const art::Ptr<simb::MCTruth>& beamNuMCTruth,
    MCParticleVector& mcParticles)
  {
    // Get the MCTruth handle
    art::Handle<std::vector<simb::MCTruth>> mcTruthHandle;
    evt.getByLabel(truthLabel, mcTruthHandle);

    if (!mcTruthHandle.isValid())
      throw cet::exception("LArPandora")
        << " LArPandoraSliceIdHelper::CollectNeutrinoMCParticles - invalid MCTruth handle"
        << std::endl;

    // Find MCParticles that are associated to the beam neutrino MCTruth block
    art::FindManyP<simb::MCParticle> truthToMCParticleAssns(mcTruthHandle, evt, mcParticleLabel);
    mcParticles = truthToMCParticleAssns.at(
      beamNuMCTruth
        .key()); // ATTN will throw if association from beamNuMCTruth doesn't exist. We want this!
  }

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

  void LArPandoraSliceIdHelper::GetHitOrigins(const art::Event& evt,
                                              const std::string& hitLabel,
                                              const std::string& backtrackLabel,
                                              const MCParticleVector& mcParticles,
                                              HitVector& hits,
                                              HitToBoolMap& hitToIsNuInducedMap)
  {
    // Collect the hits from the event
    art::Handle<std::vector<recob::Hit>> hitHandle;
    evt.getByLabel(hitLabel, hitHandle);

    if (!hitHandle.isValid())
      throw cet::exception("LArPandora")
        << " LArPandoraSliceIdHelper::GetHitOrigins - invalid hit handle" << std::endl;

    art::FindManyP<simb::MCParticle> hitToMCParticleAssns(hitHandle, evt, backtrackLabel);

    // Find the hits that are associated to a neutrino induced MCParticle using the Hit->MCParticle associations form the backtracker
    for (unsigned int i = 0; i < hitHandle->size(); ++i) {
      const art::Ptr<recob::Hit> hit(hitHandle, i);
      hits.push_back(hit);

      const auto& particles(hitToMCParticleAssns.at(hit.key()));

      bool foundNuParticle(false);
      for (const auto& part : particles) {
        // If the MCParticles isn't in the list of neutrino particles
        if (std::find(mcParticles.begin(), mcParticles.end(), part) == mcParticles.end()) continue;

        foundNuParticle = true;
        break;
      }

      if (!hitToIsNuInducedMap.emplace(hit, foundNuParticle).second)
        throw cet::exception("LArPandora")
          << " LArPandoraSliceIdHelper::GetHitOrigins - repeated hits in input collection"
          << std::endl;
    }
  }

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

  unsigned int LArPandoraSliceIdHelper::CountNeutrinoHits(const HitVector& hits,
                                                          const HitToBoolMap& hitToIsNuInducedMap)
  {
    unsigned int nNuHits(0);
    for (const auto& hit : hits) {
      const auto it(hitToIsNuInducedMap.find(hit));

      if (it == hitToIsNuInducedMap.end())
        throw cet::exception("LArPandora")
          << " LArPandoraSliceIdHelper::CountNeutrinoHits - can't find hit in hitToIsNuInducedMap"
          << std::endl;

      nNuHits += it->second ? 1 : 0;
    }

    return nNuHits;
  }

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

  void LArPandoraSliceIdHelper::GetPFParticleToHitsMap(const art::Event& evt,
                                                       const std::string& pandoraLabel,
                                                       PFParticlesToHits& pfParticleToHitsMap)
  {
    // Get the PFParticles
    art::Handle<std::vector<recob::PFParticle>> pfParticleHandle;
    evt.getByLabel(pandoraLabel, pfParticleHandle);

    if (!pfParticleHandle.isValid())
      throw cet::exception("LArPandora")
        << " LArPandoraSliceIdHelper::GetPFParticleToHitsMap - invalid PFParticle handle"
        << std::endl;

    // Get the Clusters
    art::Handle<std::vector<recob::Cluster>> clusterHandle;
    evt.getByLabel(pandoraLabel, clusterHandle);

    if (!clusterHandle.isValid())
      throw cet::exception("LArPandora")
        << " LArPandoraSliceIdHelper::GetPFParticleToHitsMap - invalid cluster handle" << std::endl;

    // Get the associations between PFParticles -> Clusters -> Hits
    art::FindManyP<recob::Cluster> pfParticleToClusterAssns(pfParticleHandle, evt, pandoraLabel);
    art::FindManyP<recob::Hit> clusterToHitAssns(clusterHandle, evt, pandoraLabel);

    // Get the hits associated to each PFParticles
    for (unsigned int iPart = 0; iPart < pfParticleHandle->size(); ++iPart) {
      const art::Ptr<recob::PFParticle> part(pfParticleHandle, iPart);
      HitVector hits;

      for (const auto& cluster : pfParticleToClusterAssns.at(part.key())) {
        for (const auto& hit : clusterToHitAssns.at(cluster.key())) {
          if (std::find(hits.begin(), hits.end(), hit) != hits.end())
            throw cet::exception("LArPandora")
              << " LArPandoraSliceIdHelper::GetPFParticleToHitsMap - double counted hits!"
              << std::endl;

          hits.push_back(hit);
        }
      }

      if (!pfParticleToHitsMap.emplace(part, hits).second)
        throw cet::exception("LArPandora")
          << " LArPandoraSliceIdHelper::GetPFParticleToHitsMap - repeated input PFParticles"
          << std::endl;
    }
  }

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

  void LArPandoraSliceIdHelper::GetReconstructedHitsInSlice(
    const Slice& slice,
    const PFParticlesToHits& pfParticleToHitsMap,
    HitVector& hits)
  {
    // ATTN here we use the PFParticles from both hypotheses to collect the hits. Hits will not be double counted
    LArPandoraSliceIdHelper::CollectHits(slice.GetTargetHypothesis(), pfParticleToHitsMap, hits);
    LArPandoraSliceIdHelper::CollectHits(slice.GetCosmicRayHypothesis(), pfParticleToHitsMap, hits);
  }

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

  void LArPandoraSliceIdHelper::CollectHits(const PFParticleVector& pfParticles,
                                            const PFParticlesToHits& pfParticleToHitsMap,
                                            HitVector& hits)
  {
    for (const auto& part : pfParticles) {
      const auto it(pfParticleToHitsMap.find(part));
      if (it == pfParticleToHitsMap.end())
        throw cet::exception("LArPandora") << " LArPandoraSliceIdHelper::CollectHits - can't find "
                                              "any hits associated to input PFParticle"
                                           << std::endl;

      for (const auto& hit : it->second) {
        // ATTN here we ensure that we don't double count hits, even if the input PFParticles are from different Pandora instances
        if (std::find(hits.begin(), hits.end(), hit) == hits.end()) hits.push_back(hit);
      }
    }
  }

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

  void LArPandoraSliceIdHelper::GetSliceMetadata(const SliceVector& slices,
                                                 const PFParticlesToHits& pfParticleToHitsMap,
                                                 const HitToBoolMap& hitToIsNuInducedMap,
                                                 const unsigned int nNuHits,
                                                 SliceMetadataVector& sliceMetadata)
  {
    if (!sliceMetadata.empty())
      throw cet::exception("LArPandora")
        << " LArPandoraSliceIdHelper::GetSliceMetadata - non empty input metadata vector"
        << std::endl;

    if (slices.empty()) return;

    unsigned int mostCompleteSliceIndex(0);
    unsigned int maxNuHits(0);

    for (unsigned int sliceIndex = 0; sliceIndex < slices.size(); ++sliceIndex) {
      const Slice& slice(slices.at(sliceIndex));
      HitVector hits;
      LArPandoraSliceIdHelper::GetReconstructedHitsInSlice(slice, pfParticleToHitsMap, hits);

      const unsigned int nHitsInSlice(hits.size());
      const unsigned int nNuHitsInSlice(
        LArPandoraSliceIdHelper::CountNeutrinoHits(hits, hitToIsNuInducedMap));

      if (nNuHitsInSlice > maxNuHits) {
        mostCompleteSliceIndex = sliceIndex;
        maxNuHits = nNuHitsInSlice;
      }

      SliceMetadata metadata;
      metadata.m_nHits = nHitsInSlice;
      metadata.m_purity = ((nHitsInSlice == 0) ?
                             -1.f :
                             static_cast<float>(nNuHitsInSlice) / static_cast<float>(nHitsInSlice));
      metadata.m_completeness =
        ((nNuHits == 0) ? -1.f : static_cast<float>(nNuHitsInSlice) / static_cast<float>(nNuHits));
      metadata.m_isMostComplete = false;

      sliceMetadata.push_back(metadata);
    }

    sliceMetadata.at(mostCompleteSliceIndex).m_isMostComplete = true;
  }

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

  LArPandoraSliceIdHelper::SliceMetadata::SliceMetadata()
    : m_purity(-std::numeric_limits<float>::max())
    , m_completeness(-std::numeric_limits<float>::max())
    , m_nHits(std::numeric_limits<unsigned int>::max())
    , m_isMostComplete(false)
  {}

} // namespace lar_pandora