From 997650e8c8fbf47eddf36246bccf98ac3b1111d4 Mon Sep 17 00:00:00 2001 From: "W. David Dagenhart" Date: Tue, 21 Jul 2026 15:24:11 -0500 Subject: [PATCH] Version of GausHitFinder with 2 unfolds --- migration/gauss_hit_finder/CMakeLists.txt | 55 ++ migration/gauss_hit_finder/README.md | 23 +- .../find_hits_with_gaussians_design1.cpp | 5 +- .../find_hits_with_gaussians_design2.cpp | 504 ++++++++++++++++++ .../find_hits_with_gaussians_design2.hpp | 137 +++++ ...ister_find_hits_with_gaussians_design2.cpp | 129 +++++ ...t_find_hits_with_gaussians_design2.jsonnet | 73 +++ 7 files changed, 920 insertions(+), 6 deletions(-) create mode 100644 migration/gauss_hit_finder/find_hits_with_gaussians_design2.cpp create mode 100644 migration/gauss_hit_finder/find_hits_with_gaussians_design2.hpp create mode 100644 migration/gauss_hit_finder/register_find_hits_with_gaussians_design2.cpp create mode 100644 migration/gauss_hit_finder/test_find_hits_with_gaussians_design2.jsonnet diff --git a/migration/gauss_hit_finder/CMakeLists.txt b/migration/gauss_hit_finder/CMakeLists.txt index db41623..19aba77 100644 --- a/migration/gauss_hit_finder/CMakeLists.txt +++ b/migration/gauss_hit_finder/CMakeLists.txt @@ -3,6 +3,7 @@ add_library(find_hits_with_gaussians SHARED find_hits_with_gaussians.cpp find_hits_with_gaussians_design1.cpp + find_hits_with_gaussians_design2.cpp copied_from_larsoft_minor_edits/Hit.cxx copied_from_larsoft_minor_edits/geo_types.cxx copied_from_larsoft_minor_edits/CandHitStandard.cxx @@ -40,6 +41,11 @@ add_library(find_hits_with_gaussians_design1_hof MODULE target_link_libraries(find_hits_with_gaussians_design1_hof PRIVATE phlex::module find_hits_with_gaussians) +add_library(find_hits_with_gaussians_design2_hof MODULE + register_find_hits_with_gaussians_design2.cpp) +target_link_libraries(find_hits_with_gaussians_design2_hof PRIVATE + phlex::module find_hits_with_gaussians) + file(GLOB test_files ${CMAKE_CURRENT_SOURCE_DIR}/*.dat) file(COPY ${test_files} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) @@ -142,6 +148,54 @@ else() endforeach() endif() +# Repeat for design2 version with unfold-transform-fold. +add_test(NAME GausHitFinderTestDesign2 + COMMAND phlex -c ${CMAKE_CURRENT_SOURCE_DIR}/test_find_hits_with_gaussians_design2.jsonnet + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +# The tests must not be run concurrently because they read and write to the same file. +set_tests_properties(GausHitFinderTestDesign2 PROPERTIES DEPENDS + "GausHitFinderCompareWithUnfold_0;GausHitFinderCompareWithUnfold_1;GausHitFinderCompareWithUnfold_2;GausHitFinderCompareWithUnfold_3;GausHitFinderCompareWithUnfold_4") + +# The reference files were generated on a Linux machine +# so we expect identical output on Linux +if(CMAKE_SYSTEM_NAME STREQUAL "Linux") + foreach(n RANGE 4) + add_test( + NAME GausHitFinderCompareDesign2_${n} + COMMAND diff -q -s + ${CMAKE_CURRENT_SOURCE_DIR}/art_hits_${n}.txt + ${CMAKE_CURRENT_BINARY_DIR}/hits_${n}.txt + ) + set_tests_properties(GausHitFinderCompareDesign2_${n} PROPERTIES DEPENDS GausHitFinderTestDesign2) + endforeach() +# On a Mac, we expect the output to be slightly +# different because the reference files were +# generated on a Linux machine. +elseif(CMAKE_SYSTEM_NAME STREQUAL "Darwin") + foreach(n RANGE 4) + add_test( + NAME GausHitFinderCompareDesign2_${n} + COMMAND ${Python3_EXECUTABLE} + ${CMAKE_CURRENT_SOURCE_DIR}/compare_hits.py + ${CMAKE_CURRENT_SOURCE_DIR}/art_hits_${n}.txt + ${CMAKE_CURRENT_BINARY_DIR}/hits_${n}.txt + ) + set_tests_properties(GausHitFinderCompareDesign2_${n} PROPERTIES DEPENDS GausHitFinderTestDesign2) + endforeach() +# Always fail on other platforms +# This will warn us to add a new comparison +# method if we ever run on a new platform. +else() + foreach(n RANGE 4) + add_test( + NAME GausHitFinderCompareDesign2_${n} + COMMAND ${CMAKE_COMMAND} -E false + ) + set_tests_properties(GausHitFinderCompareDesign2_${n} PROPERTIES DEPENDS GausHitFinderTestDesign2) + endforeach() +endif() + # We need to include the library directory for 'generate_layers' on the PHLEX_PLUGIN_PATH. # Unfortunately, phlex does not have a target that easily does this (plugins are not intended # to be accessible downstream), so we have to work around that. Right now, the library directory @@ -151,6 +205,7 @@ endif() set_tests_properties( GausHitFinderTest GausHitFinderTestWithUnfold + GausHitFinderTestDesign2 PROPERTIES ENVIRONMENT_MODIFICATION "PHLEX_PLUGIN_PATH=path_list_prepend:${CMAKE_CURRENT_BINARY_DIR}:$" diff --git a/migration/gauss_hit_finder/README.md b/migration/gauss_hit_finder/README.md index c9ecdb6..7dcc1a1 100644 --- a/migration/gauss_hit_finder/README.md +++ b/migration/gauss_hit_finder/README.md @@ -62,6 +62,18 @@ iteration is replaced by that sequence of algorithms. 3. register_find_hits_with_gaussians_design1.cpp 4. test_find_hits_with_gaussians_design1.jsonnet +These files extend design1 by also replacing the inner +parallel_for (over ROIs) with a second unfold-transform-fold. +The resulting layer hierarchy is: +spill -> wire -> roi. A `wire_roi_data` struct bundles +a single ROI (`datarange_t`) with the wire-level context +(channel, view) needed by the transform. + +1. find_hits_with_gaussians_design2.hpp +2. find_hits_with_gaussians_design2.cpp +3. register_find_hits_with_gaussians_design2.cpp +4. test_find_hits_with_gaussians_design2.jsonnet + We plan to implement more prototype migrations of GausHitFinder in the future to explore the possibilities and execute tests. @@ -114,10 +126,11 @@ Getting the single transform to work was a good first step. The next thing we implemented was replacing the outer `parallel_for` with an `unfold`, `transform`, and `fold` -sequence of algorithms. In the future, we plan to create other -versions pushing the division into separate algorithms to lower -levels. Then we plan to run tests and compare the different -versions. +sequence of algorithms (design1). After that, we created +design2 which also replaces the inner `parallel_for` (over ROIs) +with a second `unfold`, `transform`, and `fold` sequence, +resulting in a three-layer hierarchy (spill -> wire -> roi). +We plan to run tests and compare the different versions. `Phlex` will not support the `Tools` feature that existed in `art`. One possibility is that algorithm nodes scheduled @@ -223,7 +236,7 @@ Two output data members of `Hit` were ignored because they depend on the `Geomet geo::WireID fWireID; ///< WireID for the hit (Cryostat, TPC, Plane, Wire) ``` -This was done for both `phlex` versions of GausHitFinder. The output of both are identical with each other and with the `art` version. +This was done for all `phlex` versions of GausHitFinder. The output of all versions are identical with each other and with the `art` version. The `phlex` process was run with multithreading and that causes the order of `Hit` objects to vary from one execution to the next and also the order of events to vary. In the comparison the `Hit` objects were sorted and we had to be careful to compare matching events. diff --git a/migration/gauss_hit_finder/find_hits_with_gaussians_design1.cpp b/migration/gauss_hit_finder/find_hits_with_gaussians_design1.cpp index f26c397..9598b09 100644 --- a/migration/gauss_hit_finder/find_hits_with_gaussians_design1.cpp +++ b/migration/gauss_hit_finder/find_hits_with_gaussians_design1.cpp @@ -45,7 +45,7 @@ namespace examples { { // Probably eventually delete the following line // (or convert to logging utility) - std::cout << "Finding hits with Gaussians." << std::endl; + std::cout << "Finding hits with Gaussians (design 1)." << std::endl; } unfold_wire_vector_design1::const_iterator unfold_wire_vector_design1::initial_value() const @@ -59,6 +59,9 @@ namespace examples { unfold_wire_vector_design1::unfold(const_iterator current) const { recob::Wire const& wire = *current; + // Note this copies the Wire object. + // We discussed this in May and did it intentionally. + // We may revisit this decision later. return std::make_pair(++current, wire); }; diff --git a/migration/gauss_hit_finder/find_hits_with_gaussians_design2.cpp b/migration/gauss_hit_finder/find_hits_with_gaussians_design2.cpp new file mode 100644 index 0000000..c4210db --- /dev/null +++ b/migration/gauss_hit_finder/find_hits_with_gaussians_design2.cpp @@ -0,0 +1,504 @@ +// See REAMDME.md for some general comments about this example. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "copied_from_larsoft_minor_edits/geo_types.h" // geo::View_t, geo::SignalType, geo::WireID +#include "copied_from_larsoft_minor_edits/ICandidateHitFinder.h" +#include "copied_from_larsoft_minor_edits/RawTypes.h" // raw::ChannelID_t +#include "find_hits_with_gaussians_design2.hpp" + +namespace { + // This is an edited copy of the TMath::Gaus function from ROOT, since we + // don't want to depend on ROOT in this example. + double Gaus(double x, double mean, double sigma, bool norm) + { + if (sigma == 0) + return 1.e30; + double arg = (x - mean) / sigma; + // for |arg| > 39 result is zero in double precision + if (arg < -39.0 || arg > 39.0) + return 0.0; + double res = std::exp(-0.5 * arg * arg); + if (!norm) + return res; + return res / (2.50662827463100024 * sigma); //sqrt(2*Pi)=2.50662827463100024 + } +} + +namespace examples { + + // --------------------------------------------------------------- + // First unfold: vector -> individual Wire objects + // --------------------------------------------------------------- + + unfold_wire_vector_design2::unfold_wire_vector_design2( + std::vector const& wires) : + begin_{wires.begin()}, end_{wires.end()} + { + // Probably eventually delete the following line + // (or convert to logging utility) + std::cout << "Finding hits with Gaussians (design 2)." << std::endl; + } + + unfold_wire_vector_design2::const_iterator unfold_wire_vector_design2::initial_value() const + { + return begin_; + } + + bool unfold_wire_vector_design2::predicate(const_iterator current) const { return current != end_; } + + std::pair + unfold_wire_vector_design2::unfold(const_iterator current) const + { + recob::Wire const& wire = *current; + // Note this copies the Wire object. + // We discussed this in May and did it intentionally. + // We may revisit this decision later. + return std::make_pair(++current, wire); + }; + + // --------------------------------------------------------------- + // Second unfold: Wire -> individual wire_roi_data objects + // --------------------------------------------------------------- + + unfold_wire_design2::unfold_wire_design2(recob::Wire const& wire) : + wire_{wire}, n_ranges_{wire.SignalROI().n_ranges()} + {} + + unfold_wire_design2::state_type unfold_wire_design2::initial_value() const + { + return 0; + } + + bool unfold_wire_design2::predicate(state_type current) const + { + return current < n_ranges_; + } + + std::pair + unfold_wire_design2::unfold(state_type current) const + { + recob::Wire::RegionsOfInterest_t const& signalROI = wire_.SignalROI(); + wire_roi_data data{signalROI.range(current), wire_.Channel(), static_cast(wire_.View())}; + return std::make_pair(current + 1, std::move(data)); + } + + // --------------------------------------------------------------- + // Transform: processes a single ROI and returns the hits found + // --------------------------------------------------------------- + + std::vector find_hits_with_gaussians_design2( + find_hits_with_gaussians_design2_cfg const& cfg, + wire_roi_data const& roi_data, + std::vector> const& cand_hit_standard, + PeakFitterMrqdt const& peak_fitter_mrqdt, + HitFilterAlg const& hit_filter_alg) + { + std::vector hitstruct_vec; + std::vector filthitstruct_vec; + + //################################################# + //### Set the charge determination method ### + //### Default is to compute the normalized area ### + //################################################# + std::function chargeFunc = + [](double /* peakMean */, + double peakAmp, + double peakWidth, + double areaNorm, + int /* low */, + int /* hi */) { return std::sqrt(2 * std::numbers::pi) * peakAmp * peakWidth / areaNorm; }; + + //############################################## + //### Alternative is to integrate over pulse ### + //############################################## + if (cfg.area_method == 0) + chargeFunc = [](double peakMean, + double peakAmp, + double peakWidth, + double /* areaNorm */, + int low, + int hi) { + double charge(0); + for (int sigPos = low; sigPos < hi; sigPos++) + charge += peakAmp * Gaus(sigPos, peakMean, peakWidth, false); + return charge; + }; + + // --- Setting Channel Number and Signal type --- + + raw::ChannelID_t channel = roi_data.channel; + + // TODO. The plane number should come from the Geometry. + // That is not implemented yet... Below is the commented out + // version of the code that worked with the art framework. + geo::PlaneID::PlaneID_t plane = 0; + + // get the WireID for this hit + // std::vector wids = wireReadoutGeom.ChannelToWire(channel); + // for now, just take the first option returned from ChannelToWire + // geo::WireID wid = wids[0]; + // We need to know the plane to look up parameters + // geo::PlaneID::PlaneID_t plane = wid.Plane; + + // ---------------------------------------------------------- + // -- Setting the appropriate signal widths and thresholds -- + // -- for the right plane. -- + // ---------------------------------------------------------- + + auto const& range = roi_data.range; + + // ROI start time + raw::TDCtick_t roiFirstBinTick = range.begin_index(); + + // ########################################################### + // ### Scan the waveform and find candidate peaks + merge ### + // ########################################################### + + examples::ICandidateHitFinder::HitCandidateVec hitCandidateVec; + examples::ICandidateHitFinder::MergeHitCandidateVec mergedCandidateHitVec; + + cand_hit_standard.at(plane)->findHitCandidates(range, 0, channel, hitCandidateVec); + cand_hit_standard.at(plane)->MergeHitCandidates( + range, hitCandidateVec, mergedCandidateHitVec); + + // ####################################################### + // ### Lets loop over the pulses we found on this wire ### + // ####################################################### + + for (auto& mergedCands : mergedCandidateHitVec) { + int startT = mergedCands.front().startTick; + int endT = mergedCands.back().stopTick; + + // ### Putting in a protection in case things went wrong ### + // ### In the end, this primarily catches the case where ### + // ### a fake pulse is at the start of the ROI ### + if (endT - startT < 5) + continue; + + // ####################################################### + // ### Clearing the parameter vector for the new Pulse ### + // ####################################################### + + // === Setting The Number Of Gaussians to try === + int nGausForFit = mergedCands.size(); + + // ################################################## + // ### Calling the function for fitting Gaussians ### + // ################################################## + double chi2PerNDF(0.); + int NDF(1); + /*stand alone + reco_tool::IPeakFitter::PeakParamsVec peakParamsVec(nGausForFit); + */ + examples::IPeakFitter::PeakParamsVec peakParamsVec; + + // ####################################################### + // ### If # requested Gaussians is too large then punt ### + // ####################################################### + if (mergedCands.size() <= cfg.max_multi_hit) { + peak_fitter_mrqdt.findPeakParameters( + range.data(), mergedCands, peakParamsVec, chi2PerNDF, NDF); + + // If the chi2 is infinite then there is a real problem so we bail + if (!(chi2PerNDF < std::numeric_limits::infinity())) { + chi2PerNDF = 2. * cfg.chi2_ndf; + NDF = 2; + } + } + + // ####################################################### + // ### If too large then force alternate solution ### + // ### - Make n hits from pulse train where n will ### + // ### depend on the fhicl parameter fLongPulseWidth ### + // ### Also do this if chi^2 is too large ### + // ####################################################### + if (mergedCands.size() > cfg.max_multi_hit || nGausForFit * chi2PerNDF > cfg.chi2_ndf) { + int longPulseWidth = cfg.long_pulse_width_vec.at(plane); + int nHitsThisPulse = (endT - startT) / longPulseWidth; + + if (nHitsThisPulse > cfg.long_max_hits_vec.at(plane)) { + nHitsThisPulse = cfg.long_max_hits_vec.at(plane); + longPulseWidth = (endT - startT) / nHitsThisPulse; + } + + if (nHitsThisPulse * longPulseWidth < endT - startT) + nHitsThisPulse++; + + int firstTick = startT; + int lastTick = std::min(firstTick + longPulseWidth, endT); + + peakParamsVec.clear(); + nGausForFit = nHitsThisPulse; + NDF = 1.; + chi2PerNDF = chi2PerNDF > cfg.chi2_ndf ? chi2PerNDF : -1.; + + for (int hitIdx = 0; hitIdx < nHitsThisPulse; hitIdx++) { + // This hit parameters + double ROIsumADC = + std::accumulate(range.begin() + firstTick, range.begin() + lastTick, 0.); + double peakSigma = (lastTick - firstTick) / 3.; // Set the width... + double peakAmp = 0.3989 * ROIsumADC / peakSigma; // Use gaussian formulation + double peakMean = (firstTick + lastTick) / 2.; + + // Store hit params + examples::IPeakFitter::PeakFitParams_t peakParams; + + peakParams.peakCenter = peakMean; + peakParams.peakCenterError = 0.1 * peakMean; + peakParams.peakSigma = peakSigma; + peakParams.peakSigmaError = 0.1 * peakSigma; + peakParams.peakAmplitude = peakAmp; + peakParams.peakAmplitudeError = 0.1 * peakAmp; + + peakParamsVec.push_back(peakParams); + + // set for next loop + firstTick = lastTick; + lastTick = std::min(lastTick + longPulseWidth, endT); + } + } + + // ####################################################### + // ### Loop through returned peaks and make recob hits ### + // ####################################################### + + int numHits(0); + + // Make a container for what will be the filtered collection + std::vector filteredHitVec; + + float nextpeak(0); + float prevpeak(0); + float nextpeakSig(0); + float prevpeakSig(0); + float nsigmaADC(2.0); + float newright(0); + float newleft(0); + for (auto const& peakParams : peakParamsVec) { + // Extract values for this hit + float peakAmp = peakParams.peakAmplitude; + float peakMean = peakParams.peakCenter; + float peakWidth = peakParams.peakSigma; + + //std::cout<<" ans hits "< 0) { + prevpeak = (peakParamsVec.at(numHits - 1)).peakCenter; + prevpeakSig = (peakParamsVec.at(numHits - 1)).peakSigma; + //std::cout<<" ans size "<::const_iterator sumStartItr = range.begin() + startT; + std::vector::const_iterator sumEndItr = range.begin() + endT; + + //### limits for the sum of the Hit based on the gaussian peak and sigma + std::vector::const_iterator HitsumStartItr = + range.begin() + peakMean - nsigmaADC * peakWidth; + std::vector::const_iterator HitsumEndItr = + range.begin() + peakMean + nsigmaADC * peakWidth; + + if (nGausForFit > 1) { + if (numHits > 0) { + if ((peakMean - nsigmaADC * peakWidth) < (prevpeak + nsigmaADC * prevpeakSig)) { + float difPeak = peakMean - prevpeak; + float weightpeak = prevpeakSig / (prevpeakSig + peakWidth); + HitsumStartItr = range.begin() + prevpeak + difPeak * weightpeak; + newleft = prevpeak + difPeak * weightpeak; + } + } + + if (numHits < nGausForFit - 1) { + if ((peakMean + nsigmaADC * peakWidth) > (nextpeak - nsigmaADC * nextpeakSig)) { + float difPeak = nextpeak - peakMean; + float weightpeak = peakWidth / (nextpeakSig + peakWidth); + HitsumEndItr = range.begin() + peakMean + difPeak * weightpeak; + newright = peakMean + difPeak * weightpeak; + } + } + } + + //protection to avoid negative ranges + if (newright - newleft < 0) + continue; + + //avoid ranges out of ROI if it happens + if (HitsumStartItr < sumStartItr) + HitsumStartItr = sumStartItr; + + if (HitsumEndItr > sumEndItr) + HitsumEndItr = sumEndItr; + + if (HitsumStartItr > HitsumEndItr) + continue; + + // ### Sum of ADC counts + double ROIsumADC = std::accumulate(sumStartItr, sumEndItr, 0.); + double HitsumADC = std::accumulate(HitsumStartItr, HitsumEndItr, 0.); + + recob::Hit hit( + roi_data.channel, + startT + roiFirstBinTick, + endT + roiFirstBinTick, + peakMean + roiFirstBinTick, + peakMeanErr, + peakWidth, + peakAmp, + peakAmpErr, + ROIsumADC, + HitsumADC, + charge, + chargeErr, + nGausForFit, + numHits, + chi2PerNDF, + NDF, + static_cast(roi_data.view), + // Geometry system for Phlex is not yet implemented, + // so we just set signal type to 0 (kInduction) for now. + geo::kInduction, + // art::ServiceHandle()->Get().SignalType(wire.Channel()), + // wid also comes from the Geometry, which is not yet implemented, so we just set + // it to a default value for now. + geo::WireID()); + // wid); + + if (cfg.filter_hits) + filteredHitVec.push_back(hit); + + // This loop will store ALL hits + hitstruct_vec.push_back(std::move(hit)); + + numHits++; + } // <---End loop over gaussians + + // Should we filter hits? + if (cfg.filter_hits && !filteredHitVec.empty()) { + // ####################################################################### + // Is all this sorting really necessary? Would it be faster to just loop + // through the hits and perform simple cuts on amplitude and width on a + // hit-by-hit basis, either here in the module (using fPulseHeightCuts and + // fPulseWidthCuts) or in HitFilterAlg? + // ####################################################################### + + // Sort in ascending peak height + // (I believe the preceding comment is incorrect. The sort below + // is in descending order of peak height, not ascending.) + std::sort(filteredHitVec.begin(), + filteredHitVec.end(), + [](auto const& left, auto const& right) { + return left.PeakAmplitude() > right.PeakAmplitude(); + }); + + // Reject if the first hit fails the PH/wid cuts + if (filteredHitVec.front().PeakAmplitude() < cfg.pulse_height_cuts.at(plane) || + filteredHitVec.front().RMS() < cfg.pulse_width_cuts.at(plane)) + filteredHitVec.clear(); + + // Now check other hits in the snippet + if (filteredHitVec.size() > 1) { + // The largest pulse height will now be at the front... + float largestPH = filteredHitVec.front().PeakAmplitude(); + + // Find where the pulse heights drop below threshold + float threshold(cfg.pulse_ratio_cuts.at(plane)); + + std::vector::iterator smallHitItr = std::find_if( + filteredHitVec.begin(), + filteredHitVec.end(), + [largestPH, threshold](auto const& hit) { + return hit.PeakAmplitude() < 8. && hit.PeakAmplitude() / largestPH < threshold; + }); + + // Shrink to fit + if (smallHitItr != filteredHitVec.end()) + filteredHitVec.resize(std::distance(filteredHitVec.begin(), smallHitItr)); + + // Resort in time order + std::sort(filteredHitVec.begin(), + filteredHitVec.end(), + [](auto const& left, auto const& right) { + return left.PeakTime() < right.PeakTime(); + }); + } + + // Copy the hits we want to keep to the filtered hit collection + for (auto const& filteredHit : filteredHitVec) { + if (!cfg.filter_hits || hit_filter_alg.IsGoodHit(filteredHit)) { + filthitstruct_vec.push_back(std::move(filteredHit)); + } + } + } + } //<---End loop over merged candidate hits + + if (cfg.filter_hits) { + return filthitstruct_vec; + } else { + return hitstruct_vec; + } + } + + // --------------------------------------------------------------- + // Inner fold: collects hits from individual ROIs into a + // per-wire vector (roi layer -> wire layer) + // --------------------------------------------------------------- + void fold_roi_hits_design2(std::vector& hits, + std::vector const& hits_from_roi) + { + hits.insert(hits.end(), hits_from_roi.begin(), hits_from_roi.end()); + } + + // --------------------------------------------------------------- + // Outer fold: collects per-wire hit vectors into the final + // output vector (wire layer -> spill layer) + // --------------------------------------------------------------- + void fold_hits_into_vector_design2(std::vector& hits, + std::vector const& hits_from_wire) + { + hits.insert(hits.end(), hits_from_wire.begin(), hits_from_wire.end()); + } + +} // end of examples namespace diff --git a/migration/gauss_hit_finder/find_hits_with_gaussians_design2.hpp b/migration/gauss_hit_finder/find_hits_with_gaussians_design2.hpp new file mode 100644 index 0000000..cfbdac3 --- /dev/null +++ b/migration/gauss_hit_finder/find_hits_with_gaussians_design2.hpp @@ -0,0 +1,137 @@ +#ifndef PHLEX_EXAMPLES_FIND_HITS_WITH_GAUSSIANS_DESIGN2_HPP +#define PHLEX_EXAMPLES_FIND_HITS_WITH_GAUSSIANS_DESIGN2_HPP + +// Design2 extends design1 by also replacing the inner parallel_for +// (over ROIs) with a second unfold-transform-fold. + +// See REAMDME.md for some general comments about this example. + +//////////////////////////////////////////////////////////////////////// +// +// GaussHitFinder class +// +// jaasaadi@syr.edu +// +// This algorithm is designed to find hits on wires after deconvolution. +// ----------------------------------- +// This algorithm is based on the FFTHitFinder written by Brian Page, +// Michigan State University, for the ArgoNeuT experiment. +// +// +// The algorithm walks along the wire and looks for pulses above threshold +// The algorithm then attempts to fit n-gaussians to these pulses where n +// is set by the number of peaks found in the pulse +// If the Chi2/NDF returned is "bad" it attempts to fit n+1 gaussians to +// the pulse. If this is a better fit it then uses the parameters of the +// Gaussian fit to characterize the "hit" object +// +// To use this simply include the following in your producers: +// gaushit: @local::microboone_gaushitfinder +// gaushit: @local::argoneut_gaushitfinder +//////////////////////////////////////////////////////////////////////// + +#include +#include +#include +#include + +#include "copied_from_larsoft_minor_edits/CandHitStandard.h" +#include "copied_from_larsoft_minor_edits/Hit.h" +#include "copied_from_larsoft_minor_edits/HitFilterAlg.h" +#include "copied_from_larsoft_minor_edits/PeakFitterMrqdt.h" +#include "copied_from_larsoft_minor_edits/Wire.h" + +namespace examples { + + // --------------------------------------------------------------- + // A lightweight struct that bundles a single ROI (datarange_t) + // together with the wire-level context needed by the transform. + // --------------------------------------------------------------- + struct wire_roi_data { + recob::Wire::RegionsOfInterest_t::datarange_t range; + raw::ChannelID_t channel; + int view; + }; + + // --------------------------------------------------------------- + // First unfold: vector -> individual Wire objects + // (unchanged from design1) + // --------------------------------------------------------------- + class unfold_wire_vector_design2 { + public: + explicit unfold_wire_vector_design2(std::vector const& wires); + + using const_iterator = std::vector::const_iterator; + + const_iterator initial_value() const; + bool predicate(const_iterator current) const; + std::pair unfold(const_iterator current) const; + + private: + const_iterator begin_; + const_iterator end_; + }; + + // --------------------------------------------------------------- + // Second unfold: Wire -> individual wire_roi_data objects + // Replaces the inner tbb::parallel_for over ROI ranges. + // --------------------------------------------------------------- + class unfold_wire_design2 { + public: + explicit unfold_wire_design2(recob::Wire const& wire); + + using state_type = std::size_t; // index into signalROI ranges + + state_type initial_value() const; + bool predicate(state_type current) const; + std::pair unfold(state_type current) const; + + private: + recob::Wire const& wire_; + std::size_t n_ranges_; + }; + + // --------------------------------------------------------------- + // Configuration for the transform + // --------------------------------------------------------------- + struct find_hits_with_gaussians_design2_cfg { + bool filter_hits; + + std::vector long_max_hits_vec; /// long_pulse_width_vec; /// + area_norms_vec; /// pulse_height_cuts; + std::vector pulse_width_cuts; + std::vector pulse_ratio_cuts; + }; + + // --------------------------------------------------------------- + // Transform: processes a single ROI and returns the hits found + // --------------------------------------------------------------- + std::vector find_hits_with_gaussians_design2( + find_hits_with_gaussians_design2_cfg const& cfg, + wire_roi_data const& roi_data, + std::vector> const& cand_hit_standard, + PeakFitterMrqdt const& peak_fitter_mrqdt, + HitFilterAlg const& hit_filter_alg); + + // --------------------------------------------------------------- + // Inner fold: collects hits from individual ROIs into a + // per-wire vector (roi layer -> wire layer) + // --------------------------------------------------------------- + void fold_roi_hits_design2(std::vector& hits, + std::vector const& hits_from_roi); + + // --------------------------------------------------------------- + // Outer fold: collects per-wire hit vectors into the final + // output vector (wire layer -> spill layer) + // --------------------------------------------------------------- + void fold_hits_into_vector_design2(std::vector& hits, + std::vector const& hits_from_wire); +} +#endif // PHLEX_EXAMPLES_FIND_HITS_WITH_GAUSSIANS_DESIGN2_HPP diff --git a/migration/gauss_hit_finder/register_find_hits_with_gaussians_design2.cpp b/migration/gauss_hit_finder/register_find_hits_with_gaussians_design2.cpp new file mode 100644 index 0000000..aeaa6ab --- /dev/null +++ b/migration/gauss_hit_finder/register_find_hits_with_gaussians_design2.cpp @@ -0,0 +1,129 @@ +// See REAMDME.md for some general comments about this example. + +#include +#include +#include +#include +#include + +#include "phlex/concurrency.hpp" +#include "phlex/configuration.hpp" +#include "phlex/core/product_selector.hpp" +#include "phlex/module.hpp" + +#include "copied_from_larsoft_minor_edits/CandHitStandard.h" +#include "copied_from_larsoft_minor_edits/HitFilterAlg.h" +#include "copied_from_larsoft_minor_edits/PeakFitterMrqdt.h" +#include "copied_from_larsoft_minor_edits/Wire.h" +#include "find_hits_with_gaussians_design2.hpp" + +using namespace phlex; + +PHLEX_REGISTER_ALGORITHMS(m, config) +{ + auto const layer_vector_of_wires = config.get("layer_vector_of_wires"); + auto const layer_wire = config.get("layer_wire"); + auto const layer_roi = config.get("layer_roi"); + + // --------------------------------------------------------------- + // Outer unfold: spill -> wire + // --------------------------------------------------------------- + m.unfold("unfold_wire_vector_design2", + &examples::unfold_wire_vector_design2::predicate, + &examples::unfold_wire_vector_design2::unfold, + layer_wire, + concurrency::unlimited) + .input_family(product_selector{.creator = "wires", .layer = layer_vector_of_wires, .suffix = ""}); + + // --------------------------------------------------------------- + // Inner unfold: wire -> roi + // --------------------------------------------------------------- + m.unfold("unfold_wire_design2", + &examples::unfold_wire_design2::predicate, + &examples::unfold_wire_design2::unfold, + layer_roi, + concurrency::unlimited) + .input_family(product_selector{.creator = "unfold_wire_vector_design2", .layer = layer_wire}); + + // --------------------------------------------------------------- + // Configuration + // --------------------------------------------------------------- + examples::find_hits_with_gaussians_design2_cfg cfg = { + .filter_hits = config.get("filter_hits"), + .long_max_hits_vec = config.get>("long_max_hits_vec"), + .long_pulse_width_vec = config.get>("long_pulse_width_vec"), + .max_multi_hit = config.get("max_multi_hit"), + .area_method = config.get("area_method"), + .area_norms_vec = config.get>("area_norms_vec"), + .chi2_ndf = config.get("chi2_ndf"), + .pulse_height_cuts = config.get>("pulse_height_cuts"), + .pulse_width_cuts = config.get>("pulse_width_cuts"), + .pulse_ratio_cuts = config.get>("pulse_ratio_cuts") + }; + + std::vector> cand_hit_standard_vec; + + auto finder_configs = config.get("cand_hit_standard_configs"); + cand_hit_standard_vec.resize(finder_configs.keys().size()); + for (auto const& key : finder_configs.keys()) { + auto finder_config = finder_configs.get(key); + auto hit_finder = std::make_shared( + examples::CandHitStandardCfg{ + .fRoiThreshold = finder_config.get("roiThreshold") + }); + unsigned int plane = finder_config.get("Plane"); + if (plane >= cand_hit_standard_vec.size()) { + std::cerr << "Error: plane number " << plane << " is out of range for cand_hit_standard_vec of size " << cand_hit_standard_vec.size() << std::endl; + throw std::runtime_error("Invalid plane number"); + } + cand_hit_standard_vec[plane] = std::move(hit_finder); + } + + std::shared_ptr peak_fitter_mrqdt; + auto fitter_config = config.get("peak_fitter_mrqdt_config"); + peak_fitter_mrqdt = std::make_shared( + examples::PeakFitterMrqdtCfg{ + .fMinWidth = fitter_config.get("min_width"), + .fMaxWidthMult = fitter_config.get("max_width_mult"), + .fPeakRange = fitter_config.get("peak_range_fact"), + .fAmpRange = fitter_config.get("peak_amp_range") + }); + + std::shared_ptr hit_filter_alg; + auto filter_config = config.get("hit_filter_alg_config"); + hit_filter_alg = std::make_shared( + examples::HitFilterAlgCfg{ + .fMinPulseHeight = filter_config.get>("min_pulse_height"), + .fMinPulseSigma = filter_config.get>("min_pulse_sigma") + }); + + // --------------------------------------------------------------- + // Transform: processes a single ROI + // --------------------------------------------------------------- + m.transform("find_hits_with_gaussians_design2", + [cfg = std::move(cfg), + cand_hit_standard_vec = std::move(cand_hit_standard_vec), + peak_fitter_mrqdt = std::move(peak_fitter_mrqdt), + hit_filter_alg = std::move(hit_filter_alg)] + (examples::wire_roi_data const& roi_data) { + return examples::find_hits_with_gaussians_design2(cfg, + roi_data, + cand_hit_standard_vec, + *peak_fitter_mrqdt, + *hit_filter_alg); + }, + concurrency::unlimited) + .input_family(product_selector{.creator = "unfold_wire_design2", .layer = layer_roi}); + + // --------------------------------------------------------------- + // Inner fold: roi -> wire (collects hits from ROIs of one wire) + // --------------------------------------------------------------- + m.fold("fold_roi_hits_design2", examples::fold_roi_hits_design2, concurrency::serial, layer_wire) + .input_family(product_selector{.creator = "find_hits_with_gaussians_design2", .layer = layer_roi}); + + // --------------------------------------------------------------- + // Outer fold: wire -> spill (collects hits from all wires) + // --------------------------------------------------------------- + m.fold("fold_hits_into_vector_design2", examples::fold_hits_into_vector_design2, concurrency::serial, layer_vector_of_wires) + .input_family(product_selector{.creator = "fold_roi_hits_design2", .layer = layer_wire}); +} diff --git a/migration/gauss_hit_finder/test_find_hits_with_gaussians_design2.jsonnet b/migration/gauss_hit_finder/test_find_hits_with_gaussians_design2.jsonnet new file mode 100644 index 0000000..054a4d1 --- /dev/null +++ b/migration/gauss_hit_finder/test_find_hits_with_gaussians_design2.jsonnet @@ -0,0 +1,73 @@ +{ + driver: { + cpp: 'generate_layers', + layers: { + spill: { parent: 'job', total: 5, starting_number: 0 }, + }, + }, + sources: { + wires_source: { + cpp: 'wires_source', + layer: 'spill', + }, + cell_id_source: { + cpp: 'find_hits_with_gaussians_cell_id_hof', + layer: 'spill', + }, + }, + modules: { + find_hits_with_gaussians_design2_cpp: { + cpp: 'find_hits_with_gaussians_design2_hof', + layer_vector_of_wires: 'spill', + layer_wire: 'wire', + layer_roi: 'roi', + // Copied in the configuration parameters that + // I obtained from "Far Detector Simulation/Reconstruction conveners": + // Dominic Brailsford and Laura Paulucci (in cc). See Dom and Kyle + // 2/17/2026 on SLACK and also email around the same time. These were + // copied out of the expanded fcl configuration. Note that the ones + // that have different values for different planes are not used when + // filter_hits is false. + filter_hits: false, + long_max_hits_vec: [1, 1, 1], + long_pulse_width_vec: [16, 16, 16], + max_multi_hit: 4, + area_method: 0, + area_norms_vec: [13.25, 13.25, 13.25], + chi2_ndf: 50.0, + pulse_height_cuts: [3.0, 3.0, 3.0], + pulse_width_cuts: [2.0, 1.5, 1.0], + pulse_ratio_cuts: [0.35, 0.40, 0.20], + cand_hit_standard_configs: { + CandidateHitsPlane0: { + Plane: 0, + roiThreshold: 6.0, + }, + CandidateHitsPlane1: { + Plane: 1, + roiThreshold: 6.0, + }, + CandidateHitsPlane2: { + Plane: 2, + roiThreshold: 6.0, + }, + }, + peak_fitter_mrqdt_config: { + min_width: 0.5, + max_width_mult: 3.0, + peak_range_fact: 2.0, + peak_amp_range: 2.0, + }, + hit_filter_alg_config: { + min_pulse_height: [5.0, 5.0, 5.0], + min_pulse_sigma: [1.0, 1.0, 1.0], + }, + }, + print_hits_to_file_cpp: { + cpp: 'print_hits_to_file_hof', + creator: 'fold_hits_into_vector_design2', + layer: 'spill', + }, + + }, +}