Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
32 changes: 28 additions & 4 deletions sbncode/CAFMaker/CAFMakerParams.h
Comment thread
rtriozzi marked this conversation as resolved.
Original file line number Diff line number Diff line change
Expand Up @@ -350,10 +350,10 @@ 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> NuGraphSlicesLabel {
Name("NuGraphSlicesLabel"),
Comment("Label of slices that have NuGraph inference."),
"" //Empty by default, please set to e.g. art::InputTag("NCCSlices")
};

Atom<art::InputTag> NuGraphFilterLabel {
Expand All @@ -368,6 +368,30 @@ namespace caf
"" //Empty by default, please set to e.g. art::InputTag("NuGraph","semantic")
};

Atom<bool> UsePandoraAfterNuGraph {
Name("UsePandoraAfterNuGraph"),
Comment("Whether to use the second pass Pandora outputs for NuGraph reco."),
false
};

Atom<float> NuGraphFilterCut {
Name("NuGraphFilterCut"),
Comment("Cut on the NuGraph2 filter score to define hit as signal or noise."),
0.5
};

Atom<float> NuGraphHIPTagWireDist {
Name("NuGraphHIPTagWireDist"),
Comment("TPC wire distance from the vertex used to count NuGraph2–tagged HIP hits."),
10
};

Atom<float> NuGraphHIPTagTickDist {
Name("NuGraphHIPTagTickDist"),
Comment("TPC tick distance from the vertex used to count NuGraph-2–tagged HIP hits."),
50
};

Atom<string> OpFlashLabel {
Name("OpFlashLabel"),
Comment("Label of PMT flash."),
Expand Down
141 changes: 105 additions & 36 deletions sbncode/CAFMaker/CAFMaker_module.cc
Comment thread
rtriozzi marked this conversation as resolved.
Original file line number Diff line number Diff line change
Expand Up @@ -1717,33 +1717,30 @@ void CAFMaker::produce(art::Event& evt) noexcept {

// collect the TPC slices
std::vector<art::Ptr<recob::Slice>> slices;
std::vector<art::Ptr<recob::Slice>> nuGraphSlices;
std::vector<std::string> slice_tag_suffixes;
std::vector<unsigned> slice_tag_indices;
for (unsigned i_tag = 0; i_tag < pandora_tag_suffixes.size(); i_tag++) {
const std::string &pandora_tag_suffix = pandora_tag_suffixes[i_tag];

// Get a handle on the slices
art::Handle<std::vector<recob::Slice>> thisSlices;
GetByLabelStrict(evt, fParams.PFParticleLabel() + pandora_tag_suffix, thisSlices);

if (thisSlices.isValid()) {
art::fill_ptr_vector(slices, thisSlices);
if (fParams.UsePandoraAfterNuGraph()) {
nuGraphSlices = slices;
} else {
nuGraphSlices = evt.getProduct<std::vector<art::Ptr<recob::Slice>>>(fParams.NuGraphSlicesLabel().label() + pandora_tag_suffix);
}

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't understand this logic.
If fParams.UsePandoraAfterNuGraph() is true, nugraphSlices will be the same as slices (so that case could easily live outside the loop); in the other case, nuGraphSlices is overwritten with the content of each tag until the last one, hence holding only slices from the data product with the last tag.
I suspect what was meant to add slices, something like

      if (!fParams.UsePandoraAfterNuGraph()) {
        auto const& moreSlices = evt.getProduct<std::vector<art::Ptr<recob::Slice>>>(fParams.NuGraphSlicesLabel().label() + pandora_tag_suffix);
        nuGraphSlices.insert(nuGraphSlices.end(), moreSlices.begin(), moreSlices.end());
      }

for (unsigned i = 0; i < thisSlices->size(); i++) {
slice_tag_suffixes.push_back(pandora_tag_suffix);
slice_tag_indices.push_back(i_tag);
}
}
}

// 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;
Expand Down Expand Up @@ -1799,19 +1796,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
fmatch_assn_map.emplace(std::make_pair(fname_opdet, sfm_assn));
}
}

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;
Expand Down Expand Up @@ -1865,15 +1850,15 @@ void CAFMaker::produce(art::Event& evt) noexcept {
// 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);
std::vector<art::Ptr<recob::Hit>> pfphits;
std::vector<art::Ptr<recob::Cluster>> pfclusters = fmPFPClusters.at(ipf);

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Avoid copies:

Suggested change
std::vector<art::Ptr<recob::Cluster>> pfclusters = fmPFPClusters.at(ipf);
std::vector<art::Ptr<recob::Cluster>> const& 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)) {

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Avoid copies:

Suggested change
for (auto hit : fmCluHits.at(icl)) {
for (auto const& hit : fmCluHits.at(icl)) {

pfphits.push_back(hit);
}
}
fmPFPartHits.push_back(pfphits);
}
}

Expand Down Expand Up @@ -2043,8 +2028,8 @@ void CAFMaker::produce(art::Event& evt) noexcept {
// primary particle and meta-data
const recob::PFParticle *primary = (iPart == fmPFPart.size()) ? NULL : fmPFPart[iPart].get();
const larpandoraobj::PFParticleMetadata *primary_meta = (iPart == fmPFPart.size()) ? NULL : fmPFPMeta.at(iPart).at(0).get();
// get the flash match

// get the flash match
std::map<std::string, const sbn::SimpleFlashMatch*> fmatch_map;
std::map<std::string, art::FindManyP<sbn::SimpleFlashMatch>>::iterator fmatch_it;
for(fmatch_it = fmatch_assn_map.begin();fmatch_it != fmatch_assn_map.end();fmatch_it++) {
Expand All @@ -2060,12 +2045,72 @@ void CAFMaker::produce(art::Event& evt) noexcept {
}
}
}

// get the primary vertex
const recob::Vertex *vertex = (iPart == fmPFPart.size() || !fmVertex.at(iPart).size()) ? NULL : fmVertex.at(iPart).at(0).get();

//#######################################################
// Add slice info.
//#######################################################
if (std::find(nuGraphSlices.begin(), nuGraphSlices.end(), slice) != nuGraphSlices.end()) {
std::vector<art::Ptr<anab::FeatureVector<1>>> ng2_filter_vec;
std::vector<art::Ptr<anab::FeatureVector<5>>> ng2_semantic_vec;
art::FindOneP<anab::FeatureVector<1>> findOneFilter(slcHits, evt, fParams.NuGraphFilterLabel().label() + slice_tag_suff + ":" + fParams.NuGraphFilterLabel().instance());
art::FindOneP<anab::FeatureVector<5>> findOneSemantic(slcHits, evt, fParams.NuGraphSemanticLabel().label() + slice_tag_suff + ":" + fParams.NuGraphSemanticLabel().instance());
Comment on lines +2058 to +2059

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It's hard to see whether the labels are the same. Define the tag first:

Suggested change
art::FindOneP<anab::FeatureVector<1>> findOneFilter(slcHits, evt, fParams.NuGraphFilterLabel().label() + slice_tag_suff + ":" + fParams.NuGraphFilterLabel().instance());
art::FindOneP<anab::FeatureVector<5>> findOneSemantic(slcHits, evt, fParams.NuGraphSemanticLabel().label() + slice_tag_suff + ":" + fParams.NuGraphSemanticLabel().instance());
art::InputTag const filtTag{
fParams.NuGraphFilterLabel().label() + slice_tag_suff,
fParams.NuGraphFilterLabel().instance()
};
art::FindOneP<anab::FeatureVector<1>> findOneFilter(slcHits, filtTag);
art::FindOneP<anab::FeatureVector<5>> findOneSemantic(slcHits, filtTag);


// filter
if (findOneFilter.isValid()) {
ng2_filter_vec.reserve(slcHits.size());
for (size_t hitIdx = 0; hitIdx < slcHits.size(); ++hitIdx) {
ng2_filter_vec.emplace_back(findOneFilter.at(hitIdx));
}
}

// semantic tagging
if (findOneSemantic.isValid()) {
ng2_semantic_vec.reserve(slcHits.size());
for (size_t hitIdx = 0; hitIdx < slcHits.size(); ++hitIdx) {
ng2_semantic_vec.emplace_back(findOneSemantic.at(hitIdx));
}
}

// vertex projection onto the three wire planes
float vtx_wire[3];
float vtx_tick[3];

if (vertex != NULL) {

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What kind of garbage is going to end into the tree if the vertex is null? The arrays need to be initialised somehow.
This is a request, as we don't want uninitialised values to end into the trees as it compromises tree comparisons.
Also, NULL is C, and we need to keep a façade here.

Suggested change
if (vertex != NULL) {
if (vertex) {

@rtriozzi rtriozzi Jul 9, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

going through the comments right now! about this: every slice has to have a vertex, sort of by definition in Pandora (e.g., even if the reconstruction is sure that the slice is a cosmic, the vertex will just be the highest-y hit)

I think there are also protections for failed vertexing (which should be very rare: we really just need a hit to define a vertex...) within Pandora, so that we always get something when we grab the vertex in the CAF maker...

anyway, I agree it's bad practice from the CAF maker stand point (+ the NULL check was therefore strange), I'll fix this! thanks!

auto const& tpcID = geom->FindTPCAtPosition(vertex->position());
if (tpcID.isValid) {
for (geo::PlaneID const& p : wireReadout.Iterate<geo::PlaneID>()) {
auto const& planeID = geo::PlaneID{tpcID, p.Plane};
const geo::PlaneGeo& planeGeo = wireReadout.Plane(planeID);
Comment on lines +2084 to +2086

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should work:

Suggested change
for (geo::PlaneID const& p : wireReadout.Iterate<geo::PlaneID>()) {
auto const& planeID = geo::PlaneID{tpcID, p.Plane};
const geo::PlaneGeo& planeGeo = wireReadout.Plane(planeID);
for (auto const& planeGeo : wireReadout.Iterate<geo::PlaneGeo>(tpcID)) {

(plane number is planeGeo.ID().Plane).

vtx_wire[p.Plane] = planeGeo.WireCoordinate(vertex->position()); ///< wire projection
vtx_tick[p.Plane] = dprop.ConvertXToTicks(vertex->position().X(), planeID); ///< drift projection
}
}
}

if (ng2_filter_vec.size() > 0 || ng2_semantic_vec.size() > 0) {

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Semantics:

Suggested change
if (ng2_filter_vec.size() > 0 || ng2_semantic_vec.size() > 0) {
if (!ng2_filter_vec.empty() || !ng2_semantic_vec.empty()) {

FillSliceNuGraph(slcHits, ng2_filter_vec, ng2_semantic_vec, fmPFPartHits,
vtx_wire, vtx_tick, fParams.NuGraphHIPTagWireDist(), fParams.NuGraphHIPTagTickDist(),
fParams.NuGraphFilterCut(), recslc);
}

// overwrite the `ng_filt_pass_frac` variable with post-NuGraph2-filtering
if (fParams.UsePandoraAfterNuGraph()) {
art::Handle<std::vector<anab::FeatureVector<1>>> ngFilterVecHandle;
evt.getByLabel(art::InputTag("NGMultiSlice" + slice_tag_suff, "filter"), ngFilterVecHandle);
Comment on lines +2101 to +2102

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Smoother (less repetitions, less maintenance):

Suggested change
art::Handle<std::vector<anab::FeatureVector<1>>> ngFilterVecHandle;
evt.getByLabel(art::InputTag("NGMultiSlice" + slice_tag_suff, "filter"), ngFilterVecHandle);
auto ngFilterVecHandle = evt.getHandle<std::vector<anab::FeatureVector<1>>>
(art::InputTag("NGMultiSlice" + slice_tag_suff, "filter"));

art::Handle<std::vector<recob::Hit>> ngFilteredHitsHandle;
GetByLabelStrict(evt, "ngfilteredhits" + slice_tag_suff, ngFilteredHitsHandle);

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this "strict" and the previous is not?

if (ngFilterVecHandle.isValid() &&
ngFilteredHitsHandle.isValid() &&
ngFilterVecHandle->size() > 0) {
recslc.ng_filt_pass_frac = float(ngFilteredHitsHandle->size()) / float(ngFilterVecHandle->size());

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Single precision floats have only <7 digits, it may lose information when there are millions of hits. Probably not a problem.
Anyway, since the hardware likely converts to intermediate double anyway...

Suggested change
recslc.ng_filt_pass_frac = float(ngFilteredHitsHandle->size()) / float(ngFilterVecHandle->size());
recslc.ng_filt_pass_frac = double(ngFilteredHitsHandle->size()) / double(ngFilterVecHandle->size());

}
}

}

FillSliceVars(*slice, primary, producer, recslc);
FillSliceMetadata(primary_meta, recslc);
FillSliceFlashMatch(fmatch_map["fmatch"], recslc.fmatch);
Expand Down Expand Up @@ -2189,8 +2234,32 @@ 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 (std::find(nuGraphSlices.begin(), nuGraphSlices.end(), slice) != nuGraphSlices.end()) {
std::vector<art::Ptr<recob::Hit>>& PFPHits = fmPFPartHits.at(iPart);
art::FindOneP<anab::FeatureVector<1>> findOneFilter(PFPHits, evt, fParams.NuGraphFilterLabel().label() + slice_tag_suff + ":" + fParams.NuGraphFilterLabel().instance());
art::FindOneP<anab::FeatureVector<5>> findOneSemantic(PFPHits, evt, fParams.NuGraphSemanticLabel().label() + slice_tag_suff + ":" + fParams.NuGraphSemanticLabel().instance());
std::vector<art::Ptr<anab::FeatureVector<1>>> ng2_filter_vec;
std::vector<art::Ptr<anab::FeatureVector<5>>> ng2_semantic_vec;

// filter
if (findOneFilter.isValid()) {
ng2_filter_vec.reserve(PFPHits.size());
for (size_t hitIdx = 0; hitIdx < PFPHits.size(); ++hitIdx) {
ng2_filter_vec.emplace_back(findOneFilter.at(hitIdx));
}
}

// semantic tagging
if (findOneSemantic.isValid()) {
ng2_semantic_vec.reserve(PFPHits.size());
for (size_t hitIdx = 0; hitIdx < PFPHits.size(); ++hitIdx) {
ng2_semantic_vec.emplace_back(findOneSemantic.at(hitIdx));
}
}
Comment on lines +2238 to +2258

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This feels quite familiar... same comments as for twin code above.
Also, consider moving into a member function

std::tuple<
  std::vector<art::Ptr<anab::FeatureVector<1>>>,
  std::vector<art::Ptr<anab::FeatureVector<5>>>
  >
CAFMaker::readNuGraphHitFeatures(
  std::vector<art::Ptr<recob::Hit>> const& hits,
  art::Event const& event, std::string const& suffix
) const;


if (ng2_filter_vec.size() > 0 || ng2_semantic_vec.size() > 0) {
FillPFPNuGraph(PFPHits, ng2_filter_vec, ng2_semantic_vec, fParams.NuGraphFilterCut(), pfp);
}
}

if (!thisTrack.empty()) { // it has a track!
Expand Down
Loading