diff --git a/sbncode/CAFMaker/CAFMakerParams.h b/sbncode/CAFMaker/CAFMakerParams.h index dceaa3951..253f0d4c2 100644 --- a/sbncode/CAFMaker/CAFMakerParams.h +++ b/sbncode/CAFMaker/CAFMakerParams.h @@ -350,6 +350,24 @@ namespace caf "" //Empty by default, configured in icaruscode cafmaker_defs }; + Atom NuGraphSliceHitLabel { + Name("NuGraphSliceHitLabel"), + Comment("Label of NuGraph slice hit map."), + "" //Empty by default, please set to e.g. art::InputTag("nuslhits") + }; + + Atom NuGraphFilterLabel { + Name("NuGraphFilterLabel"), + Comment("Label of NuGraph filter."), + "" //Empty by default, please set to e.g. art::InputTag("NuGraph","filter") + }; + + Atom NuGraphSemanticLabel { + Name("NuGraphSemanticLabel"), + Comment("Label of NuGraph semantic."), + "" //Empty by default, please set to e.g. art::InputTag("NuGraph","semantic") + }; + Atom OpFlashLabel { Name("OpFlashLabel"), Comment("Label of PMT flash."), diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index c30f694b2..48da84e8f 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/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> > ng2_slice_hit_map_handle(pandora_tag_suffixes.size()); + std::vector< art::Handle>> > ng2_filter_handle(pandora_tag_suffixes.size()); + std::vector< art::Handle>> > 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>> ng2_filter_vec; + std::vector>> 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 fmOpT0 = FindManyPStrict(sliceList, evt, fParams.OpT0Label() + slice_tag_suff); std::vector> slcOpT0; @@ -1833,6 +1858,25 @@ void CAFMaker::produce(art::Event& evt) noexcept { art::FindManyP fmSpacePointPFPs = FindManyPStrict(slcSpacePoints, evt, fParams.PFParticleLabel() + slice_tag_suff); + art::FindManyP fmPFPClusters = + FindManyPStrict(fmPFPart, evt, fParams.PFParticleLabel() + slice_tag_suff); + + std::vector>> fmPFPartHits; + // make Ptr's to clusters for cluster -> other object associations + if (fmPFPClusters.isValid()) { + for (size_t ipf=0; ipf> pfphits; + std::vector> pfclusters = fmPFPClusters.at(ipf); + art::FindManyP fmCluHits = FindManyPStrict(pfclusters, evt, fParams.PFParticleLabel() + slice_tag_suff); + for (size_t icl=0; icl fmShower = FindManyPStrict(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); diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx index d36e3b719..d81cb5890 100644 --- a/sbncode/CAFMaker/FillReco.cxx +++ b/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> &inputHits, + const std::vector &sliceHitsMap, + const std::vector>> &ngFilterResult, + const std::vector>> &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 &sliceHitsMap, + const std::vector>> &ngFilterResult, + const std::vector>> &ngSemanticResult, + const std::vector> &pfpHits, + caf::SRPFP& srpfp, + bool allowEmpty) + { + + // the nugraph elements are ordered the same as the sliceHitsMap + std::vector 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 ng2sempfpcounts(5,0); + size_t ng2bkgpfpcount = 0; + for (size_t pos : mappedhits) { + auto const& bkgscore = ngFilterResult.at(pos); + if (bkgscore->at(0) ng2semscores; + for (size_t i=0;isize();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, diff --git a/sbncode/CAFMaker/FillReco.h b/sbncode/CAFMaker/FillReco.h index 0e5ec9052..6e8ecf292 100644 --- a/sbncode/CAFMaker/FillReco.h +++ b/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> &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> &inputHits, + const std::vector &sliceHitsMap, + const std::vector>> &ngFilterResult, + const std::vector>> &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 &sliceHitsMap, + const std::vector>> &ngFilterResult, + const std::vector>> &ngSemanticResult, + const std::vector> &pfpHits, + caf::SRPFP& srpfp, + bool allowEmpty= false); + void FillTrackCRTHit(const std::vector> &t0match, const std::vector> &hitmatch, const std::vector> &hitmatchinfo,