Merge pull request #532 from cerati/feature/cerati_ng2caf

Feature/cerati ng2caf

Author: kjplows

sbncode/CAFMaker/CAFMakerParams.h

@@ -350,6 +350,24 @@ namespace caf
       "" //Empty by default, configured in icaruscode cafmaker_defs
     };
 
+    Atom<art::InputTag> NuGraphSliceHitLabel {
+      Name("NuGraphSliceHitLabel"),
+      Comment("Label of NuGraph slice hit map."),
+      "" //Empty by default, please set to e.g. art::InputTag("nuslhits")
+    };
+
+    Atom<art::InputTag> NuGraphFilterLabel {
+      Name("NuGraphFilterLabel"),
+      Comment("Label of NuGraph filter."),
+      "" //Empty by default, please set to e.g. art::InputTag("NuGraph","filter")
+    };
+
+    Atom<art::InputTag> NuGraphSemanticLabel {
+      Name("NuGraphSemanticLabel"),
+      Comment("Label of NuGraph semantic."),
+      "" //Empty by default, please set to e.g. art::InputTag("NuGraph","semantic")
+    };
+
     Atom<string> OpFlashLabel {
       Name("OpFlashLabel"),
       Comment("Label of PMT flash."),

sbncode/CAFMaker/CAFMaker_module.cc

@@ -95,6 +95,8 @@
 #include "lardataobj/RecoBase/Vertex.h"
 #include "lardataobj/RecoBase/Shower.h"
 #include "lardataobj/RecoBase/MCSFitResult.h"
+#include "lardataobj/RecoBase/Cluster.h"
+#include "lardataobj/AnalysisBase/MVAOutput.h"
 
 #include "nusimdata/SimulationBase/MCFlux.h"
 #include "nusimdata/SimulationBase/MCTruth.h"
@@ -1731,6 +1733,17 @@ void CAFMaker::produce(art::Event& evt) noexcept {
     }
   }
 
+  // nu graph
+  std::vector< art::Handle<std::vector<unsigned int>> > ng2_slice_hit_map_handle(pandora_tag_suffixes.size());
+  std::vector< art::Handle<std::vector<anab::FeatureVector<1>>> > ng2_filter_handle(pandora_tag_suffixes.size());
+  std::vector< art::Handle<std::vector<anab::FeatureVector<5>>> > ng2_semantic_handle(pandora_tag_suffixes.size());
+  for (unsigned i_tag = 0; i_tag < pandora_tag_suffixes.size(); i_tag++) {
+    const std::string &pandora_tag_suffix = pandora_tag_suffixes[i_tag];
+    GetByLabelIfExists(evt, fParams.NuGraphSliceHitLabel().encode() + pandora_tag_suffix, ng2_slice_hit_map_handle[i_tag]);
+    GetByLabelIfExists(evt, fParams.NuGraphFilterLabel().label() + pandora_tag_suffix + ":" + fParams.NuGraphFilterLabel().instance(), ng2_filter_handle[i_tag]);
+    GetByLabelIfExists(evt, fParams.NuGraphSemanticLabel().label() + pandora_tag_suffix + ":" + fParams.NuGraphSemanticLabel().instance(), ng2_semantic_handle[i_tag]);
+  }
+
   // The Standard Record
   // Branch entry definition -- contains list of slices, CRT information, and truth information
   StandardRecord rec;
@@ -1787,6 +1800,18 @@ void CAFMaker::produce(art::Event& evt) noexcept {
       }
     }
 
+    std::vector<art::Ptr<anab::FeatureVector<1>>> ng2_filter_vec;
+    std::vector<art::Ptr<anab::FeatureVector<5>>> ng2_semantic_vec;
+    if (ng2_filter_handle[producer].isValid()) {
+      art::fill_ptr_vector(ng2_filter_vec,ng2_filter_handle[producer]);
+    }
+    if (ng2_semantic_handle[producer].isValid()) {
+      art::fill_ptr_vector(ng2_semantic_vec,ng2_semantic_handle[producer]);
+    }
+    if (ng2_slice_hit_map_handle[producer].isValid()) {
+      FillSliceNuGraph(slcHits,*ng2_slice_hit_map_handle[producer],ng2_filter_vec,ng2_semantic_vec,recslc);
+    }
+
     art::FindManyP<sbn::OpT0Finder> fmOpT0 =
       FindManyPStrict<sbn::OpT0Finder>(sliceList, evt, fParams.OpT0Label() + slice_tag_suff);
     std::vector<art::Ptr<sbn::OpT0Finder>> slcOpT0;
@@ -1833,6 +1858,25 @@ void CAFMaker::produce(art::Event& evt) noexcept {
     art::FindManyP<recob::PFParticle> fmSpacePointPFPs =
       FindManyPStrict<recob::PFParticle>(slcSpacePoints, evt, fParams.PFParticleLabel() + slice_tag_suff);
 
+    art::FindManyP<recob::Cluster> fmPFPClusters =
+      FindManyPStrict<recob::Cluster>(fmPFPart, evt, fParams.PFParticleLabel() + slice_tag_suff);
+
+    std::vector<std::vector<art::Ptr<recob::Hit>>> fmPFPartHits;
+    // make Ptr's to clusters for cluster -> other object associations
+    if (fmPFPClusters.isValid()) {
+      for (size_t ipf=0; ipf<fmPFPart.size();++ipf) {
+	std::vector<art::Ptr<recob::Hit>> pfphits;
+	std::vector<art::Ptr<recob::Cluster>> pfclusters = fmPFPClusters.at(ipf);
+	art::FindManyP<recob::Hit> fmCluHits = FindManyPStrict<recob::Hit>(pfclusters, evt, fParams.PFParticleLabel() + slice_tag_suff);
+	for (size_t icl=0; icl<fmCluHits.size();icl++) {
+	  for (auto hit : fmCluHits.at(icl)) {
+	    pfphits.push_back(hit);
+	  }
+	}
+	fmPFPartHits.push_back(pfphits);
+      }
+    }
+
     art::FindManyP<recob::Shower> fmShower =
       FindManyPStrict<recob::Shower>(fmPFPart, evt, fParams.RecoShowerLabel() + slice_tag_suff);
 
@@ -2145,6 +2189,10 @@ void CAFMaker::produce(art::Event& evt) noexcept {
         FillCNNScores(thisParticle, cnnScores, pfp);
       }
 
+      if (ng2_slice_hit_map_handle[producer].isValid()) {
+	FillPFPNuGraph(*ng2_slice_hit_map_handle[producer], ng2_filter_vec, ng2_semantic_vec, fmPFPartHits.at(iPart), pfp);
+      }
+
       if (!thisTrack.empty())  { // it has a track!
         assert(thisTrack.size() == 1);
 

sbncode/CAFMaker/FillReco.cxx

@@ -12,6 +12,7 @@
 
 namespace caf
 {
+  const float ng_filter_cut = 0.5;
 
   //......................................................................
   bool SelectSlice(const caf::SRSlice &slice, bool cut_clear_cosmic) {
@@ -579,6 +580,25 @@ namespace caf
     }
   }
 
+  void FillSliceNuGraph(const std::vector<art::Ptr<recob::Hit>> &inputHits,
+			const std::vector<unsigned int> &sliceHitsMap,
+			const std::vector<art::Ptr<anab::FeatureVector<1>>> &ngFilterResult,
+			const std::vector<art::Ptr<anab::FeatureVector<5>>> &ngSemanticResult,
+			caf::SRSlice &slice)
+  {
+
+    //need to double check that the slice processed by NuGraph is the same under consideration
+    //std::cout << "sizes=" << inputHits.size() << " " << sliceHitsMap.size() << " " << ngFilterResult.size() << " " << ngSemanticResult.size() << std::endl;
+    unsigned int nHits = inputHits.size();
+    if (nHits==0 || nHits!=sliceHitsMap.size() || inputHits[0].key()!=sliceHitsMap[0]) return;//not the same slice!
+
+    unsigned int npass = 0;
+    for ( unsigned int i = 0; i < nHits; i++ ) {
+      if (ngFilterResult.at(i)->at(0)>=ng_filter_cut) npass++;
+    }
+    slice.ng_filt_pass_frac = float(npass)/nHits;
+  }
+
   //......................................................................
 
 
@@ -1004,6 +1024,59 @@ namespace caf
     srpfp.cnnscore.nclusters = cnnscore->nClusters;
   }
 
+  void FillPFPNuGraph(const std::vector<unsigned int> &sliceHitsMap,
+		      const std::vector<art::Ptr<anab::FeatureVector<1>>> &ngFilterResult,
+		      const std::vector<art::Ptr<anab::FeatureVector<5>>> &ngSemanticResult,
+		      const std::vector<art::Ptr<recob::Hit>> &pfpHits,
+		      caf::SRPFP& srpfp,
+		      bool allowEmpty)
+  {
+
+    // the nugraph elements are ordered the same as the sliceHitsMap
+    std::vector<size_t> mappedhits;
+    for (auto& hit : pfpHits) {
+      auto it = std::find(sliceHitsMap.begin(), sliceHitsMap.end(), hit.key());
+      if (it != sliceHitsMap.end()) {
+	size_t index = std::distance(sliceHitsMap.begin(), it);
+	mappedhits.push_back(index);
+      }
+    }
+
+    if (mappedhits.size()>0) {
+      std::vector<float> ng2sempfpcounts(5,0);
+      size_t ng2bkgpfpcount = 0;
+      for (size_t pos : mappedhits) {
+	auto const& bkgscore = ngFilterResult.at(pos);
+	if (bkgscore->at(0)<ng_filter_cut) {
+	  ng2bkgpfpcount++;
+	} else {
+	  auto const& scores = ngSemanticResult.at(pos);
+	  std::vector<float> ng2semscores;
+	  for (size_t i=0;i<scores->size();i++) ng2semscores.push_back(scores->at(i));
+	  size_t sem_label = std::distance(ng2semscores.begin(), std::max_element(ng2semscores.begin(), ng2semscores.end()));//arg_max(ng2semscores);
+	  ng2sempfpcounts[sem_label]++;
+	}
+      }
+      srpfp.ngscore.sem_cat = SRNuGraphScore::NuGraphCategory(std::distance(ng2sempfpcounts.begin(), std::max_element(ng2sempfpcounts.begin(), ng2sempfpcounts.end())));//arg_max(ng2sempfpcounts);
+      size_t nonBkgHits = (pfpHits.size() > ng2bkgpfpcount ? pfpHits.size()-ng2bkgpfpcount : 0);
+      srpfp.ngscore.mip_frac = (nonBkgHits>0 ? float(ng2sempfpcounts[0])/nonBkgHits : -1.);
+      srpfp.ngscore.hip_frac = (nonBkgHits>0 ? float(ng2sempfpcounts[1])/nonBkgHits : -1.);
+      srpfp.ngscore.shr_frac = (nonBkgHits>0 ? float(ng2sempfpcounts[2])/nonBkgHits : -1.);
+      srpfp.ngscore.mhl_frac = (nonBkgHits>0 ? float(ng2sempfpcounts[3])/nonBkgHits : -1.);
+      srpfp.ngscore.dif_frac = (nonBkgHits>0 ? float(ng2sempfpcounts[4])/nonBkgHits : -1.);
+      srpfp.ngscore.bkg_frac = float(ng2bkgpfpcount)/pfpHits.size();
+    } else {
+      srpfp.ngscore.sem_cat = SRNuGraphScore::NuGraphCategory::Unset;
+      srpfp.ngscore.mip_frac = -1.;
+      srpfp.ngscore.hip_frac = -1.;
+      srpfp.ngscore.shr_frac = -1.;
+      srpfp.ngscore.mhl_frac = -1.;
+      srpfp.ngscore.dif_frac = -1.;
+      srpfp.ngscore.bkg_frac = -1.;
+    }
+
+  }
+
   void FillHitVars(const recob::Hit& hit,
                    unsigned producer,
                    const recob::SpacePoint& spacepoint,

sbncode/CAFMaker/FillReco.h

@@ -20,6 +20,7 @@
 #include "lardataobj/RecoBase/OpFlash.h"
 #include "lardataobj/RecoBase/OpHit.h"
 #include "lardataobj/AnalysisBase/Calorimetry.h"
+#include "lardataobj/AnalysisBase/MVAOutput.h"
 #include "lardataobj/AnalysisBase/ParticleID.h"
 #include "lardataobj/AnalysisBase/T0.h"
 #include "lardataobj/RecoBase/PFParticleMetadata.h"
@@ -105,6 +106,22 @@ namespace caf
                            const std::vector<art::Ptr<recob::SpacePoint>> &inputPoints,
                            caf::SRSlice &slice);
 
+  /**
+   * @brief Fills the results from NuGraph at slice level
+   * @param inputHits (pointers to) the hits associated to the slice
+   * @param sliceHitsMap maps position of hits in collection input to NuGraph (slice only) to the one input to Pandora (all gaus hits)
+   * @param ngFilterResult NuGraph filter result, for each hit
+   * @param ngSemanticResult NuGraph semnatic result, for each hit (MIP track, HIP, shower, Michel electron, diffuse activity)
+   * @param[out] slice the destination slice object
+   *
+   * Hits with filter value (`ngFilterResult`) lower than `ng_filter_cut` are counted as background.
+   */
+  void FillSliceNuGraph(const std::vector<art::Ptr<recob::Hit>> &inputHits,
+			const std::vector<unsigned int> &sliceHitsMap,
+			const std::vector<art::Ptr<anab::FeatureVector<1>>> &ngFilterResult,
+			const std::vector<art::Ptr<anab::FeatureVector<5>>> &ngSemanticResult,
+			caf::SRSlice &slice);
+
   bool SelectSlice(const caf::SRSlice &slice, bool cut_clear_cosmic);
 
   void FillTrackVars(const recob::Track& track,
@@ -131,6 +148,23 @@ namespace caf
                      caf::SRPFP& srpfp,
                      bool allowEmpty = false);
 
+  /**
+   * @brief Fills the results from NuGraph at slice level
+   * @param sliceHitsMap maps position of hits in collection input to NuGraph (slice only) to the one input to Pandora (all gaus hits)
+   * @param ngFilterResult NuGraph filter result, for each hit
+   * @param ngSemanticResult NuGraph semnatic result, for each hit (MIP track, HIP, shower, Michel electron, diffuse activity)
+   * @param pfpHits Vector of hits associated to the PFParticle
+   * @param[out] srpfp the destination PFParticle object
+   *
+   * Hits with filter value (`ngFilterResult`) lower than `ng_filter_cut` are counted as background.
+   */
+  void FillPFPNuGraph(const std::vector<unsigned int> &sliceHitsMap,
+		      const std::vector<art::Ptr<anab::FeatureVector<1>>> &ngFilterResult,
+		      const std::vector<art::Ptr<anab::FeatureVector<5>>> &ngSemanticResult,
+		      const std::vector<art::Ptr<recob::Hit>> &pfpHits,
+		      caf::SRPFP& srpfp,
+		      bool allowEmpty= false);
+
   void FillTrackCRTHit(const std::vector<art::Ptr<anab::T0>> &t0match,
                        const std::vector<art::Ptr<sbn::crt::CRTHit>> &hitmatch,
                        const std::vector<art::Ptr<sbn::crt::CRTHitT0TaggingInfo>> &hitmatchinfo,