Skip to content

Commit d5d32a0

Browse files
rtriozzirtriozzi
committed
fill nugraph hip tagger
1 parent 78a0667 commit d5d32a0

3 files changed

Lines changed: 132 additions & 48 deletions

File tree

‎sbncode/CAFMaker/CAFMaker_module.cc‎

Lines changed: 92 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -1794,38 +1794,38 @@ void CAFMaker::produce(art::Event& evt) noexcept {
17941794
}
17951795
}
17961796

1797-
if (std::find(nuGraphSlices.begin(), nuGraphSlices.end(), slice) != nuGraphSlices.end()) {
1798-
std::vector<art::Ptr<anab::FeatureVector<1>>> ng2_filter_vec;
1799-
std::vector<art::Ptr<anab::FeatureVector<5>>> ng2_semantic_vec;
1800-
art::FindOneP<anab::FeatureVector<1>> findOneFilter(slcHits, evt, fParams.NuGraphFilterLabel().label() + slice_tag_suff + ":" + fParams.NuGraphFilterLabel().instance());
1801-
art::FindOneP<anab::FeatureVector<5>> findOneSemantic(slcHits, evt, fParams.NuGraphSemanticLabel().label() + slice_tag_suff + ":" + fParams.NuGraphSemanticLabel().instance());
1802-
// iteration is the only way to get something out of a FindOne or FindMany...
1803-
if (findOneFilter.isValid()) {
1804-
for (size_t hitIdx = 0; hitIdx < slcHits.size(); ++hitIdx) {
1805-
if (findOneFilter.at(hitIdx).isNull()) {
1806-
slcHits.erase(slcHits.begin()+hitIdx);
1807-
hitIdx--;
1808-
continue;
1809-
}
1810-
ng2_filter_vec.emplace_back(findOneFilter.at(hitIdx));
1811-
}
1812-
}
1813-
1814-
if (findOneSemantic.isValid()) {
1815-
for (size_t hitIdx = 0; hitIdx < slcHits.size(); ++hitIdx) {
1816-
if (findOneSemantic.at(hitIdx).isNull()) {
1817-
slcHits.erase(slcHits.begin()+hitIdx);
1818-
hitIdx--;
1819-
continue;
1820-
}
1821-
ng2_semantic_vec.emplace_back(findOneSemantic.at(hitIdx));
1822-
}
1823-
}
1824-
1825-
if (ng2_filter_vec.size() > 0 || ng2_semantic_vec.size() > 0) {
1826-
FillSliceNuGraph(slcHits, ng2_filter_vec, ng2_semantic_vec, recslc);
1827-
}
1828-
}
1797+
// if (std::find(nuGraphSlices.begin(), nuGraphSlices.end(), slice) != nuGraphSlices.end()) {
1798+
// std::vector<art::Ptr<anab::FeatureVector<1>>> ng2_filter_vec;
1799+
// std::vector<art::Ptr<anab::FeatureVector<5>>> ng2_semantic_vec;
1800+
// art::FindOneP<anab::FeatureVector<1>> findOneFilter(slcHits, evt, fParams.NuGraphFilterLabel().label() + slice_tag_suff + ":" + fParams.NuGraphFilterLabel().instance());
1801+
// art::FindOneP<anab::FeatureVector<5>> findOneSemantic(slcHits, evt, fParams.NuGraphSemanticLabel().label() + slice_tag_suff + ":" + fParams.NuGraphSemanticLabel().instance());
1802+
1803+
// if (findOneFilter.isValid()) {
1804+
// for (size_t hitIdx = 0; hitIdx < slcHits.size(); ++hitIdx) {
1805+
// if (findOneFilter.at(hitIdx).isNull()) {
1806+
// slcHits.erase(slcHits.begin()+hitIdx);
1807+
// hitIdx--;
1808+
// continue;
1809+
// }
1810+
// ng2_filter_vec.emplace_back(findOneFilter.at(hitIdx));
1811+
// }
1812+
// }
1813+
1814+
// if (findOneSemantic.isValid()) {
1815+
// for (size_t hitIdx = 0; hitIdx < slcHits.size(); ++hitIdx) {
1816+
// if (findOneSemantic.at(hitIdx).isNull()) {
1817+
// slcHits.erase(slcHits.begin()+hitIdx);
1818+
// hitIdx--;
1819+
// continue;
1820+
// }
1821+
// ng2_semantic_vec.emplace_back(findOneSemantic.at(hitIdx));
1822+
// }
1823+
// }
1824+
1825+
// if (ng2_filter_vec.size() > 0 || ng2_semantic_vec.size() > 0) {
1826+
// FillSliceNuGraph(slcHits, ng2_filter_vec, ng2_semantic_vec, recslc);
1827+
// }
1828+
// }
18291829

18301830
art::FindManyP<sbn::OpT0Finder> fmOpT0 =
18311831
FindManyPStrict<sbn::OpT0Finder>(sliceList, evt, fParams.OpT0Label() + slice_tag_suff);
@@ -1880,15 +1880,15 @@ void CAFMaker::produce(art::Event& evt) noexcept {
18801880
// make Ptr's to clusters for cluster -> other object associations
18811881
if (fmPFPClusters.isValid()) {
18821882
for (size_t ipf=0; ipf<fmPFPart.size();++ipf) {
1883-
std::vector<art::Ptr<recob::Hit>> pfphits;
1884-
std::vector<art::Ptr<recob::Cluster>> pfclusters = fmPFPClusters.at(ipf);
1885-
art::FindManyP<recob::Hit> fmCluHits = FindManyPStrict<recob::Hit>(pfclusters, evt, fParams.PFParticleLabel() + slice_tag_suff);
1886-
for (size_t icl=0; icl<fmCluHits.size();icl++) {
1887-
for (auto hit : fmCluHits.at(icl)) {
1888-
pfphits.push_back(hit);
1889-
}
1890-
}
1891-
fmPFPartHits.push_back(pfphits);
1883+
std::vector<art::Ptr<recob::Hit>> pfphits;
1884+
std::vector<art::Ptr<recob::Cluster>> pfclusters = fmPFPClusters.at(ipf);
1885+
art::FindManyP<recob::Hit> fmCluHits = FindManyPStrict<recob::Hit>(pfclusters, evt, fParams.PFParticleLabel() + slice_tag_suff);
1886+
for (size_t icl=0; icl<fmCluHits.size();icl++) {
1887+
for (auto hit : fmCluHits.at(icl)) {
1888+
pfphits.push_back(hit);
1889+
}
1890+
}
1891+
fmPFPartHits.push_back(pfphits);
18921892
}
18931893
}
18941894

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

2062+
// get the flash match
20632063
std::map<std::string, const sbn::SimpleFlashMatch*> fmatch_map;
20642064
std::map<std::string, art::FindManyP<sbn::SimpleFlashMatch>>::iterator fmatch_it;
20652065
for(fmatch_it = fmatch_assn_map.begin();fmatch_it != fmatch_assn_map.end();fmatch_it++) {
@@ -2075,12 +2075,62 @@ void CAFMaker::produce(art::Event& evt) noexcept {
20752075
}
20762076
}
20772077
}
2078+
20782079
// get the primary vertex
20792080
const recob::Vertex *vertex = (iPart == fmPFPart.size() || !fmVertex.at(iPart).size()) ? NULL : fmVertex.at(iPart).at(0).get();
20802081

20812082
//#######################################################
20822083
// Add slice info.
20832084
//#######################################################
2085+
if (std::find(nuGraphSlices.begin(), nuGraphSlices.end(), slice) != nuGraphSlices.end()) {
2086+
std::vector<art::Ptr<anab::FeatureVector<1>>> ng2_filter_vec;
2087+
std::vector<art::Ptr<anab::FeatureVector<5>>> ng2_semantic_vec;
2088+
art::FindOneP<anab::FeatureVector<1>> findOneFilter(slcHits, evt, fParams.NuGraphFilterLabel().label() + slice_tag_suff + ":" + fParams.NuGraphFilterLabel().instance());
2089+
art::FindOneP<anab::FeatureVector<5>> findOneSemantic(slcHits, evt, fParams.NuGraphSemanticLabel().label() + slice_tag_suff + ":" + fParams.NuGraphSemanticLabel().instance());
2090+
2091+
// filter
2092+
if (findOneFilter.isValid()) {
2093+
for (size_t hitIdx = 0; hitIdx < slcHits.size(); ++hitIdx) {
2094+
if (findOneFilter.at(hitIdx).isNull()) {
2095+
slcHits.erase(slcHits.begin()+hitIdx);
2096+
hitIdx--;
2097+
continue;
2098+
}
2099+
ng2_filter_vec.emplace_back(findOneFilter.at(hitIdx));
2100+
}
2101+
}
2102+
2103+
// semantic tagging
2104+
if (findOneSemantic.isValid()) {
2105+
for (size_t hitIdx = 0; hitIdx < slcHits.size(); ++hitIdx) {
2106+
if (findOneSemantic.at(hitIdx).isNull()) {
2107+
slcHits.erase(slcHits.begin()+hitIdx);
2108+
hitIdx--;
2109+
continue;
2110+
}
2111+
ng2_semantic_vec.emplace_back(findOneSemantic.at(hitIdx));
2112+
}
2113+
}
2114+
2115+
// vertex projection onto the three wire planes
2116+
float vtx_wire[3];
2117+
float vtx_tick[3];
2118+
2119+
if (vertex != NULL) {
2120+
auto const& tpcID = geom->FindTPCAtPosition(vertex->position());
2121+
for (geo::PlaneID const& p : wireReadout.Iterate<geo::PlaneID>()) {
2122+
auto const& planeID = geo::PlaneID{tpcID, p.Plane};
2123+
const geo::PlaneGeo& planeGeo = wireReadout.Plane(planeID);
2124+
vtx_wire[p.Plane] = planeGeo.WireCoordinate(vertex->position()); ///< wire projection
2125+
vtx_tick[p.Plane] = dprop.ConvertXToTicks(vertex->position().X(), planeID); ///< drift projection
2126+
}
2127+
}
2128+
2129+
if (ng2_filter_vec.size() > 0 || ng2_semantic_vec.size() > 0) {
2130+
FillSliceNuGraph(slcHits, ng2_filter_vec, ng2_semantic_vec, vtx_wire, vtx_tick, recslc);
2131+
}
2132+
}
2133+
20842134
FillSliceVars(*slice, primary, producer, recslc);
20852135
FillSliceMetadata(primary_meta, recslc);
20862136
FillSliceFlashMatch(fmatch_map["fmatch"], recslc.fmatch);

‎sbncode/CAFMaker/FillReco.cxx‎

Lines changed: 36 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -583,20 +583,51 @@ namespace caf
583583
void FillSliceNuGraph(const std::vector<art::Ptr<recob::Hit>> &inputHits,
584584
const std::vector<art::Ptr<anab::FeatureVector<1>>> &ngFilterResult,
585585
const std::vector<art::Ptr<anab::FeatureVector<5>>> &ngSemanticResult,
586+
const float vtx_wire[3],
587+
const float vtx_tick[3],
586588
caf::SRSlice &slice)
587589
{
588590

589-
//need to double check that the slice processed by NuGraph is the same under consideration
590-
// std::cout << "sizes=" << inputHits.size() << " " << " " << ngFilterResult.size() << " " << ngSemanticResult.size() << std::endl;
591+
// nugraph filter fraction
591592
unsigned int nHits = inputHits.size();
592-
// implementing this same-slice check with the updated architecture would be a major headache and I believe is not required given how inputHits is retrieved.
593-
// if (nHits==0 || nHits!=sliceHitsMap.size() || inputHits[0].key()!=sliceHitsMap[0]) return;//not the same slice!
594-
595593
unsigned int npass = 0;
594+
595+
std::cout << "number of hits: " << nHits << std::endl;
596+
596597
for ( unsigned int i = 0; i < nHits; i++ ) {
597598
if (ngFilterResult.at(i)->at(0)>=ng_filter_cut) npass++;
598599
}
599600
slice.ng_filt_pass_frac = float(npass)/nHits;
601+
602+
// nugraph HIP vertex tagging
603+
for (unsigned int plane = 0; plane < 3; ++plane) {
604+
605+
unsigned int nHIPHits = 0;
606+
607+
for ( unsigned int i = 0; i < nHits; i++ ) {
608+
const recob::Hit& hit = *inputHits.at(i);
609+
if (hit.WireID().Plane != plane) continue;
610+
611+
auto const& sem = *ngSemanticResult.at(i);
612+
std::vector<float> semVec;
613+
for (size_t k = 0; k < sem.size(); ++k) semVec.push_back(sem.at(k));
614+
auto highestScoreIdx = std::distance(
615+
semVec.begin(),
616+
std::max_element(semVec.begin(), semVec.end())
617+
);
618+
if (highestScoreIdx != 1) continue; // look for HIP hits
619+
620+
float dwire = std::abs(float(hit.WireID().Wire) - vtx_wire[plane]);
621+
if (dwire > 10) continue;
622+
623+
float dtick = std::abs(hit.PeakTime() - vtx_tick[plane]);
624+
if (dtick > 50) continue;
625+
626+
nHIPHits += 1;
627+
}
628+
629+
slice.ng_plane[plane].ng_vtx_hip_hits = nHIPHits;
630+
}
600631
}
601632

602633
//......................................................................

‎sbncode/CAFMaker/FillReco.h‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -109,16 +109,19 @@ namespace caf
109109
/**
110110
* @brief Fills the results from NuGraph at slice level
111111
* @param inputHits (pointers to) the hits associated to the slice
112-
* @param sliceHitsMap maps position of hits in collection input to NuGraph (slice only) to the one input to Pandora (all gaus hits)
113112
* @param ngFilterResult NuGraph filter result, for each hit
114113
* @param ngSemanticResult NuGraph semnatic result, for each hit (MIP track, HIP, shower, Michel electron, diffuse activity)
114+
* @param vtx_wire vertex coordinates projected onto wires, per plane
115+
* @param vtx_tick vertex coordinates projected onto ticks, per plane
115116
* @param[out] slice the destination slice object
116117
*
117118
* Hits with filter value (`ngFilterResult`) lower than `ng_filter_cut` are counted as background.
118119
*/
119120
void FillSliceNuGraph(const std::vector<art::Ptr<recob::Hit>> &inputHits,
120121
const std::vector<art::Ptr<anab::FeatureVector<1>>> &ngFilterResult,
121122
const std::vector<art::Ptr<anab::FeatureVector<5>>> &ngSemanticResult,
123+
const float vtx_wire[3],
124+
const float vtx_tick[3],
122125
caf::SRSlice &slice);
123126

124127
bool SelectSlice(const caf::SRSlice &slice, bool cut_clear_cosmic);

0 commit comments

Comments
 (0)