PFParticleCosmicAna_module.cc

Go to the documentation of this file.

 
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
 
#include "larpandora/LArPandoraInterface/LArPandoraHelper.h"
 
#include "TTree.h"
 
#include <string>
 
//------------------------------------------------------------------------------------------------------------------------------------------
 
namespace lar_pandora {
 
  class PFParticleCosmicAna : public art::EDAnalyzer {
  public:
    PFParticleCosmicAna(fhicl::ParameterSet const& pset);
 
    virtual ~PFParticleCosmicAna();
 
    void beginJob();
    void endJob();
    void analyze(const art::Event& evt);
    void reconfigure(fhicl::ParameterSet const& pset);
 
  private:
    void FillRecoTree(const PFParticlesToHits& recoParticlesToHits,
                      const PFParticlesToTracks& recoParticlesToTracks,
                      const TracksToCosmicTags& recoTracksToCosmicTags);
 
    void FillTrueTree(const HitVector& hitVector,
                      const HitsToMCParticles& trueHitsToParticles,
                      const HitsToPFParticles& recoHitsToParticles,
                      const MCParticlesToMCTruth& particlesToTruth,
                      const PFParticlesToTracks& particlesToTracks,
                      const TracksToCosmicTags& tracksToCosmicTags);
 
    float GetCosmicScore(const art::Ptr<recob::PFParticle> particle,
                         const PFParticlesToTracks& recoParticlesToTracks,
                         const TracksToCosmicTags& recoTracksToCosmicTags) const;
 
    TTree* m_pRecoTree; 
    TTree* m_pTrueTree; 
 
    int m_run;   
    int m_event; 
    int m_index; 
 
    int m_self;          
    int m_pdgCode;       
    int m_isTrackLike;   
    int m_isPrimary;     
    float m_cosmicScore; 
    int m_nTracks;       
    int m_nHits;         
 
    float m_trackVtxX;       
    float m_trackVtxY;       
    float m_trackVtxZ;       
    float m_trackEndX;       
    float m_trackEndY;       
    float m_trackEndZ;       
    float m_trackVtxDirX;    
    float m_trackVtxDirY;    
    float m_trackVtxDirZ;    
    float m_trackEndDirX;    
    float m_trackEndDirY;    
    float m_trackEndDirZ;    
    float m_trackLength;     
    float m_trackWidthX;     
    float m_trackWidthY;     
    float m_trackWidthZ;     
    float m_trackVtxDeltaYZ; 
    float m_trackEndDeltaYZ; 
 
    int m_trackVtxContained; 
    int m_trackEndContained; 
 
    int m_nNeutrinoHits;
    int m_nNeutrinoHitsFullyTagged;
    int m_nNeutrinoHitsSemiTagged;
    int m_nNeutrinoHitsNotTagged;
    int m_nNeutrinoHitsNotReconstructed;
    int m_nNeutrinoHitsReconstructed;
 
    int m_nCosmicHits;
    int m_nCosmicHitsFullyTagged;
    int m_nCosmicHitsSemiTagged;
    int m_nCosmicHitsNotTagged;
    int m_nCosmicHitsNotReconstructed;
    int m_nCosmicHitsReconstructed;
 
    std::string m_hitfinderLabel;   
    std::string m_trackfitLabel;    
    std::string m_particleLabel;    
    std::string m_cosmicLabel;      
    std::string m_geantModuleLabel; 
 
    bool m_useDaughterPFParticles; 
    bool m_useDaughterMCParticles; 
 
    double m_cosmicContainmentCut; 
  };
 
  DEFINE_ART_MODULE(PFParticleCosmicAna)
 
} // namespace lar_pandora
 
//------------------------------------------------------------------------------------------------------------------------------------------
// implementation follows
 
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art_root_io/TFileDirectory.h"
#include "art_root_io/TFileService.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
 
#include "larcore/Geometry/Geometry.h"
#include "lardataobj/AnalysisBase/CosmicTag.h"
#include "lardataobj/RecoBase/PFParticle.h"
#include "lardataobj/RecoBase/Track.h"
#include "nusimdata/SimulationBase/MCTruth.h"
 
#include <iostream>
 
namespace lar_pandora {
 
  PFParticleCosmicAna::PFParticleCosmicAna(fhicl::ParameterSet const& pset) : art::EDAnalyzer(pset)
  {
    this->reconfigure(pset);
  }
 
  //------------------------------------------------------------------------------------------------------------------------------------------
 
  PFParticleCosmicAna::~PFParticleCosmicAna() {}
 
  //------------------------------------------------------------------------------------------------------------------------------------------
 
  void PFParticleCosmicAna::reconfigure(fhicl::ParameterSet const& pset)
  {
    m_cosmicLabel = pset.get<std::string>("CosmicTagModule", "cosmictagger");
    m_particleLabel = pset.get<std::string>("PFParticleModule", "pandora");
    m_trackfitLabel = pset.get<std::string>("TrackFitModule", "trackfit");
    m_hitfinderLabel = pset.get<std::string>("HitFinderModule", "gaushit");
    m_geantModuleLabel = pset.get<std::string>("GeantModule", "largeant");
 
    m_useDaughterPFParticles = pset.get<bool>("UseDaughterPFParticles", true);
    m_useDaughterMCParticles = pset.get<bool>("UseDaughterMCParticles", true);
 
    m_cosmicContainmentCut = pset.get<double>("CosmicContainmentCut", 5.0);
  }
 
  //------------------------------------------------------------------------------------------------------------------------------------------
 
  void PFParticleCosmicAna::beginJob()
  {
    mf::LogDebug("LArPandora") << " *** PFParticleCosmicAna::beginJob() *** " << std::endl;
 
    //
    art::ServiceHandle<art::TFileService const> tfs;
 
    m_pRecoTree = tfs->make<TTree>("recoTree", "LAr Cosmic Reco Tree");
    m_pRecoTree->Branch("run", &m_run, "run/I");
    m_pRecoTree->Branch("event", &m_event, "event/I");
    m_pRecoTree->Branch("index", &m_index, "index/I");
    m_pRecoTree->Branch("self", &m_self, "self/I");
    m_pRecoTree->Branch("pdgCode", &m_pdgCode, "pdgCode/I");
    m_pRecoTree->Branch("isTrackLike", &m_isTrackLike, "isTrackLike/I");
    m_pRecoTree->Branch("isPrimary", &m_isPrimary, "isPrimary/I");
    m_pRecoTree->Branch("cosmicScore", &m_cosmicScore, "cosmicScore/F");
    m_pRecoTree->Branch("trackVtxX", &m_trackVtxX, "trackVtxX/F");
    m_pRecoTree->Branch("trackVtxY", &m_trackVtxY, "trackVtxY/F");
    m_pRecoTree->Branch("trackVtxZ", &m_trackVtxZ, "trackVtxZ/F");
    m_pRecoTree->Branch("trackEndX", &m_trackEndX, "trackEndX/F");
    m_pRecoTree->Branch("trackEndY", &m_trackEndY, "trackEndY/F");
    m_pRecoTree->Branch("trackEndZ", &m_trackEndZ, "trackEndZ/F");
    m_pRecoTree->Branch("trackVtxDirX", &m_trackVtxDirX, "trackVtxDirX/F");
    m_pRecoTree->Branch("trackVtxDirY", &m_trackVtxDirY, "trackVtxDirY/F");
    m_pRecoTree->Branch("trackVtxDirZ", &m_trackVtxDirZ, "trackVtxDirZ/F");
    m_pRecoTree->Branch("trackEndDirX", &m_trackEndDirX, "trackEndDirX/F");
    m_pRecoTree->Branch("trackEndDirY", &m_trackEndDirY, "trackEndDirY/F");
    m_pRecoTree->Branch("trackEndDirZ", &m_trackEndDirZ, "trackEndDirZ/F");
    m_pRecoTree->Branch("trackLength", &m_trackLength, "trackLength/F");
    m_pRecoTree->Branch("trackWidthX", &m_trackWidthX, "trackWidthX/F");
    m_pRecoTree->Branch("trackWidthY", &m_trackWidthY, "trackWidthY/F");
    m_pRecoTree->Branch("trackWidthZ", &m_trackWidthZ, "trackWidthZ/F");
    m_pRecoTree->Branch("trackVtxDeltaYZ", &m_trackVtxDeltaYZ, "trackVtxDeltaYZ/F");
    m_pRecoTree->Branch("trackEndDeltaYZ", &m_trackEndDeltaYZ, "trackEndDeltaYZ/F");
    m_pRecoTree->Branch("trackVtxContained", &m_trackVtxContained, "trackVtxContained/I");
    m_pRecoTree->Branch("trackEndContained", &m_trackEndContained, "trackEndContained/I");
    m_pRecoTree->Branch("nTracks", &m_nTracks, "nTracks/I");
    m_pRecoTree->Branch("nHits", &m_nHits, "nHits/I");
 
    m_pTrueTree = tfs->make<TTree>("trueTree", "LAr Cosmic True Tree");
    m_pTrueTree->Branch("run", &m_run, "run/I");
    m_pTrueTree->Branch("event", &m_event, "event/I");
    m_pTrueTree->Branch("nHits", &m_nHits, "nHits/I");
    m_pTrueTree->Branch("nNeutrinoHits", &m_nNeutrinoHits, "nNeutrinoHits/I");
    m_pTrueTree->Branch(
      "nNeutrinoHitsFullyTagged", &m_nNeutrinoHitsFullyTagged, "nNeutrinoHitsFullyTagged/I");
    m_pTrueTree->Branch(
      "nNeutrinoHitsSemiTagged", &m_nNeutrinoHitsSemiTagged, "nNeutrinoHitsSemiTagged/I");
    m_pTrueTree->Branch(
      "nNeutrinoHitsNotTagged", &m_nNeutrinoHitsNotTagged, "nNeutrinoHitsNotTagged/I");
    m_pTrueTree->Branch("nNeutrinoHitsNotReconstructed",
                        &m_nNeutrinoHitsNotReconstructed,
                        "nNeutrinoHitsNotReconstructed/I");
    m_pTrueTree->Branch(
      "nNeutrinoHitsReconstructed", &m_nNeutrinoHitsReconstructed, "nNeutrinoHitsReconstructed/I");
    m_pTrueTree->Branch("nCosmicHits", &m_nCosmicHits, "nCosmicHits/I");
    m_pTrueTree->Branch(
      "nCosmicHitsFullyTagged", &m_nCosmicHitsFullyTagged, "nCosmicHitsFullyTagged/I");
    m_pTrueTree->Branch(
      "nCosmicHitsSemiTagged", &m_nCosmicHitsSemiTagged, "nCosmicHitsSemiTagged/I");
    m_pTrueTree->Branch("nCosmicHitsNotTagged", &m_nCosmicHitsNotTagged, "nCosmicHitsNotTagged/I");
    m_pTrueTree->Branch("nCosmicHitsNotReconstructed",
                        &m_nCosmicHitsNotReconstructed,
                        "nCosmicHitsNotReconstructed/I");
    m_pTrueTree->Branch(
      "nCosmicHitsReconstructed", &m_nCosmicHitsReconstructed, "nCosmicHitsReconstructed/I");
  }
 
  //------------------------------------------------------------------------------------------------------------------------------------------
 
  void PFParticleCosmicAna::endJob() {}
 
  //------------------------------------------------------------------------------------------------------------------------------------------
 
  void PFParticleCosmicAna::analyze(const art::Event& evt)
  {
    std::cout << " *** PFParticleCosmicAna::analyze(...) *** " << std::endl;
 
    //
    // Note: I've made this is MicroBooNE-only module
    //
 
    m_run = evt.run();
    m_event = evt.id().event();
 
    std::cout << "  Run: " << m_run << std::endl;
    std::cout << "  Event: " << m_event << std::endl;
 
    // Collect True Particles
    // ======================
    HitVector hitVector;
    MCTruthToMCParticles truthToParticles;
    MCParticlesToMCTruth particlesToTruth;
    MCParticlesToHits trueParticlesToHits;
    HitsToMCParticles trueHitsToParticles;
 
    LArPandoraHelper::CollectHits(evt, m_hitfinderLabel, hitVector);
    LArPandoraHelper::CollectMCParticles(
      evt, m_geantModuleLabel, truthToParticles, particlesToTruth);
    LArPandoraHelper::BuildMCParticleHitMaps(evt,
                                             m_geantModuleLabel,
                                             hitVector,
                                             trueParticlesToHits,
                                             trueHitsToParticles,
                                             (m_useDaughterMCParticles ?
                                                LArPandoraHelper::kAddDaughters :
                                                LArPandoraHelper::kIgnoreDaughters));
 
    // Collect Reco Particles
    // ======================
    PFParticleVector recoParticleVector;
    PFParticlesToHits recoParticlesToHits;
    HitsToPFParticles recoHitsToParticles;
 
    LArPandoraHelper::CollectPFParticles(evt, m_particleLabel, recoParticleVector);
    LArPandoraHelper::BuildPFParticleHitMaps(evt,
                                             m_particleLabel,
                                             m_particleLabel,
                                             recoParticlesToHits,
                                             recoHitsToParticles,
                                             (m_useDaughterPFParticles ?
                                                LArPandoraHelper::kAddDaughters :
                                                LArPandoraHelper::kIgnoreDaughters));
 
    std::cout << "  PFParticles: " << recoParticleVector.size() << std::endl;
 
    // Collect Reco Tracks
    // ===================
    TrackVector recoTrackVector;
    PFParticlesToTracks recoParticlesToTracks;
    LArPandoraHelper::CollectTracks(evt, m_trackfitLabel, recoTrackVector, recoParticlesToTracks);
 
    // Collect Cosmic Tags
    // =====================
    CosmicTagVector recoCosmicTagVector;
    TracksToCosmicTags recoTracksToCosmicTags;
    LArPandoraHelper::CollectCosmicTags(
      evt, m_cosmicLabel, recoCosmicTagVector, recoTracksToCosmicTags);
 
    // Analyse Reconstructed Particles
    // ===============================
    this->FillRecoTree(recoParticlesToHits, recoParticlesToTracks, recoTracksToCosmicTags);
 
    // Analyse True Hits
    // =================
    this->FillTrueTree(hitVector,
                       trueHitsToParticles,
                       recoHitsToParticles,
                       particlesToTruth,
                       recoParticlesToTracks,
                       recoTracksToCosmicTags);
  }
 
  //------------------------------------------------------------------------------------------------------------------------------------------
 
  void PFParticleCosmicAna::FillRecoTree(const PFParticlesToHits& recoParticlesToHits,
                                         const PFParticlesToTracks& recoParticlesToTracks,
                                         const TracksToCosmicTags& recoTracksToCosmicTags)
  {
    // Set up Geometry Service
    // =======================
    art::ServiceHandle<geo::Geometry const> theGeometry;
 
    const double xmin(0.0);
    const double xmax(2.0 * theGeometry->DetHalfWidth());
    const double ymin(-theGeometry->DetHalfHeight());
    const double ymax(+theGeometry->DetHalfHeight());
    const double zmin(0.0);
    const double zmax(theGeometry->DetLength());
    const double xyzCut(m_cosmicContainmentCut);
 
    m_index = 0;
 
    m_self = 0;
    m_pdgCode = 0;
    m_isTrackLike = 0;
    m_isPrimary = 0;
    m_cosmicScore = 0.f;
 
    m_trackVtxX = 0.f;
    m_trackVtxY = 0.f;
    m_trackVtxZ = 0.f;
    m_trackEndX = 0.f;
    m_trackEndY = 0.f;
    m_trackEndZ = 0.f;
    m_trackVtxDirX = 0.f;
    m_trackVtxDirY = 0.f;
    m_trackVtxDirZ = 0.f;
    m_trackEndDirX = 0.f;
    m_trackEndDirY = 0.f;
    m_trackEndDirZ = 0.f;
    m_trackLength = 0.f;
    m_trackWidthX = 0.f;
    m_trackWidthY = 0.f;
    m_trackWidthZ = 0.f;
    m_trackVtxDeltaYZ = 0.f;
    m_trackEndDeltaYZ = 0.f;
 
    m_trackVtxContained = 0;
    m_trackEndContained = 0;
 
    m_nTracks = 0;
    m_nHits = 0;
 
    // Loop over Reco Particles
    // ========================
    for (PFParticlesToHits::const_iterator iter1 = recoParticlesToHits.begin(),
                                           iterEnd1 = recoParticlesToHits.end();
         iter1 != iterEnd1;
         ++iter1) {
      const art::Ptr<recob::PFParticle> recoParticle = iter1->first;
 
      const HitVector& hitVector = iter1->second;
      if (hitVector.empty()) continue;
 
      PFParticlesToTracks::const_iterator iter2 = recoParticlesToTracks.find(recoParticle);
      if (recoParticlesToTracks.end() == iter2) continue;
 
      const TrackVector& trackVector = iter2->second;
      if (trackVector.empty()) continue;
 
      m_nHits = hitVector.size();
      m_nTracks = trackVector.size();
 
      m_self = recoParticle->Self();
      m_pdgCode = recoParticle->PdgCode();
      m_isPrimary = recoParticle->IsPrimary();
      m_isTrackLike = LArPandoraHelper::IsTrack(recoParticle);
      m_cosmicScore =
        this->GetCosmicScore(recoParticle, recoParticlesToTracks, recoTracksToCosmicTags);
 
      m_trackVtxX = 0.f;
      m_trackVtxY = 0.f;
      m_trackVtxZ = 0.f;
      m_trackEndX = 0.f;
      m_trackEndY = 0.f;
      m_trackEndZ = 0.f;
      m_trackVtxDirX = 0.f;
      m_trackVtxDirY = 0.f;
      m_trackVtxDirZ = 0.f;
      m_trackEndDirX = 0.f;
      m_trackEndDirY = 0.f;
      m_trackEndDirZ = 0.f;
      m_trackLength = 0.f;
      m_trackWidthX = 0.f;
      m_trackWidthY = 0.f;
      m_trackWidthZ = 0.f;
      m_trackVtxDeltaYZ = 0.f;
      m_trackEndDeltaYZ = 0.f;
 
      m_trackVtxContained = 0;
      m_trackEndContained = 0;
 
      for (TrackVector::const_iterator iter3 = trackVector.begin(), iterEnd3 = trackVector.end();
           iter3 != iterEnd3;
           ++iter3) {
        const art::Ptr<recob::Track> track = *iter3;
        const float trackLength(track->Length());
 
        if (trackLength < m_trackLength) continue;
 
        m_trackLength = trackLength;
 
        const auto& trackVtxPosition = track->Vertex();
        const auto& trackVtxDirection = track->VertexDirection();
        const auto& trackEndPosition = track->End();
        const auto& trackEndDirection = track->EndDirection();
 
        m_trackVtxX = trackVtxPosition.x();
        m_trackVtxY = trackVtxPosition.y();
        m_trackVtxZ = trackVtxPosition.z();
        m_trackVtxDirX = trackVtxDirection.x();
        m_trackVtxDirY = trackVtxDirection.y();
        m_trackVtxDirZ = trackVtxDirection.z();
        m_trackEndX = trackEndPosition.x();
        m_trackEndY = trackEndPosition.y();
        m_trackEndZ = trackEndPosition.z();
        m_trackEndDirX = trackEndDirection.x();
        m_trackEndDirY = trackEndDirection.y();
        m_trackEndDirZ = trackEndDirection.z();
 
        m_trackWidthX = std::fabs(m_trackEndX - m_trackVtxX);
        m_trackWidthY = std::fabs(m_trackEndY - m_trackVtxY);
        m_trackWidthZ = std::fabs(m_trackEndZ - m_trackVtxZ);
 
        m_trackVtxDeltaYZ =
          std::min((ymax - m_trackVtxY), std::min((m_trackVtxZ - zmin), (zmax - m_trackVtxZ)));
        m_trackEndDeltaYZ =
          std::min((m_trackEndY - ymin), std::min((m_trackEndZ - zmin), (zmax - m_trackEndZ)));
 
        m_trackVtxContained = ((m_trackVtxX > xmin + xyzCut && m_trackVtxX < xmax - xyzCut) &&
                               (m_trackVtxY > ymin + xyzCut && m_trackVtxY < ymax - xyzCut) &&
                               (m_trackVtxZ > zmin + xyzCut && m_trackVtxZ < zmax - xyzCut));
        m_trackEndContained = ((m_trackEndX > xmin + xyzCut && m_trackEndX < xmax - xyzCut) &&
                               (m_trackEndY > ymin + xyzCut && m_trackEndY < ymax - xyzCut) &&
                               (m_trackEndZ > zmin + xyzCut && m_trackEndZ < zmax - xyzCut));
      }
 
      std::cout << "   PFParticle: [" << m_index << "] nHits=" << m_nHits
                << ", nTracks=" << m_nTracks << ", cosmicScore=" << m_cosmicScore << std::endl;
 
      m_pRecoTree->Fill();
      ++m_index;
    }
  }
 
  //------------------------------------------------------------------------------------------------------------------------------------------
 
  void PFParticleCosmicAna::FillTrueTree(const HitVector& hitVector,
                                         const HitsToMCParticles& trueHitsToParticles,
                                         const HitsToPFParticles& recoHitsToParticles,
                                         const MCParticlesToMCTruth& particlesToTruth,
                                         const PFParticlesToTracks& particlesToTracks,
                                         const TracksToCosmicTags& tracksToCosmicTags)
  {
    m_nHits = 0;
 
    m_nNeutrinoHits = 0;
    m_nNeutrinoHitsFullyTagged = 0;
    m_nNeutrinoHitsSemiTagged = 0;
    m_nNeutrinoHitsNotTagged = 0;
    m_nNeutrinoHitsNotReconstructed = 0;
    m_nNeutrinoHitsReconstructed = 0;
 
    m_nCosmicHits = 0;
    m_nCosmicHitsFullyTagged = 0;
    m_nCosmicHitsSemiTagged = 0;
    m_nCosmicHitsNotTagged = 0;
    m_nCosmicHitsNotReconstructed = 0;
    m_nCosmicHitsReconstructed = 0;
 
    for (HitVector::const_iterator iter2 = hitVector.begin(), iterEnd2 = hitVector.end();
         iter2 != iterEnd2;
         ++iter2) {
      const art::Ptr<recob::Hit> hit = *iter2;
 
      HitsToMCParticles::const_iterator iter3 = trueHitsToParticles.find(hit);
      if (trueHitsToParticles.end() == iter3) continue;
 
      const art::Ptr<simb::MCParticle> trueParticle = iter3->second;
 
      MCParticlesToMCTruth::const_iterator iter4 = particlesToTruth.find(trueParticle);
      if (particlesToTruth.end() == iter4)
        throw cet::exception("LArPandora") << " PFParticleCosmicAna::analyze --- Found a true "
                                              "particle without any ancestry information ";
 
      const art::Ptr<simb::MCTruth> truth = iter4->second;
 
      float cosmicScore(-0.2);
 
      HitsToPFParticles::const_iterator iter5 = recoHitsToParticles.find(hit);
      if (recoHitsToParticles.end() != iter5) {
        const art::Ptr<recob::PFParticle> particle = iter5->second;
        cosmicScore = this->GetCosmicScore(particle, particlesToTracks, tracksToCosmicTags);
      }
 
      ++m_nHits;
 
      if (truth->NeutrinoSet()) {
        ++m_nNeutrinoHits;
 
        if (cosmicScore >= 0)
          ++m_nNeutrinoHitsReconstructed;
        else
          ++m_nNeutrinoHitsNotReconstructed;
 
        if (cosmicScore > 0.51)
          ++m_nNeutrinoHitsFullyTagged;
        else if (cosmicScore > 0.39)
          ++m_nNeutrinoHitsSemiTagged;
        else
          ++m_nNeutrinoHitsNotTagged;
      }
      else {
        ++m_nCosmicHits;
 
        if (cosmicScore >= 0)
          ++m_nCosmicHitsReconstructed;
        else
          ++m_nCosmicHitsNotReconstructed;
 
        if (cosmicScore > 0.51)
          ++m_nCosmicHitsFullyTagged;
        else if (cosmicScore > 0.39)
          ++m_nCosmicHitsSemiTagged;
        else
          ++m_nCosmicHitsNotTagged;
      }
    }
 
    m_pTrueTree->Fill();
  }
 
  //------------------------------------------------------------------------------------------------------------------------------------------
 
  float PFParticleCosmicAna::GetCosmicScore(const art::Ptr<recob::PFParticle> particle,
                                            const PFParticlesToTracks& recoParticlesToTracks,
                                            const TracksToCosmicTags& recoTracksToCosmicTags) const
  {
    float cosmicScore(0.f);
 
    // Get cosmic tags associated with this particle
    PFParticlesToTracks::const_iterator iter2 = recoParticlesToTracks.find(particle);
    if (recoParticlesToTracks.end() != iter2) {
      for (TrackVector::const_iterator iter3 = iter2->second.begin(),
                                       iterEnd3 = iter2->second.end();
           iter3 != iterEnd3;
           ++iter3) {
        const art::Ptr<recob::Track> track = *iter3;
 
        TracksToCosmicTags::const_iterator iter4 = recoTracksToCosmicTags.find(track);
        if (recoTracksToCosmicTags.end() != iter4) {
          for (CosmicTagVector::const_iterator iter5 = iter4->second.begin(),
                                               iterEnd5 = iter4->second.end();
               iter5 != iterEnd5;
               ++iter5) {
            const art::Ptr<anab::CosmicTag> cosmicTag = *iter5;
            if (cosmicTag->CosmicScore() > cosmicScore) cosmicScore = cosmicTag->CosmicScore();
          }
        }
      }
    }
 
    return cosmicScore;
  }
 
} //namespace lar_pandora