Skip to content

Per Plane Hit Tuning Parameters for GaussHitFinderSBN #626

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
carriganm95 wants to merge 2 commits into
base: develop
Choose a base branch
from
Open

Per Plane Hit Tuning Parameters for GaussHitFinderSBN #626

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
805020e
update GaussHitFinderSBN_module to take per plane hit finding parameters
Jan 21, 2026
5977b29
Updates from Mattia to allow only one definition of MaxMultiHitPerPla…
Mar 4, 2026
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
2 changes: 2 additions & 0 deletions sbncode/CAFMaker/RecoUtils/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,4 +11,6 @@ art_make_library( LIBRARY_NAME caf_RecoUtils
larsim::MCCheater_BackTrackerService_service
larsim::MCCheater_ParticleInventoryService_service
larcorealg::Geometry
sbnanaobj::StandardRecord
sbnanaobj::StandardRecordFlat
Comment on lines +14 to +15
Copy link

Copilot AI Feb 4, 2026

Choose a reason for hiding this comment

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

The CMakeLists.txt changes add sbnanaobj StandardRecord dependencies, but these don't appear to be used by RecoUtils.cc or RecoUtils.h. Neither file includes or references StandardRecord or StandardRecordFlat. This change appears unrelated to the hit finder tuning described in the PR description. If these dependencies are needed for other purposes, that should be explained. If not, they should be removed as they add unnecessary dependencies and build-time coupling.

Copilot uses AI. Check for mistakes.
Copy link
Member

@PetrilloAtWork PetrilloAtWork Mar 23, 2026

Choose a reason for hiding this comment

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

From the header included in RecoUtils.{h,cpp} it looks like this dependency is still not needed. I request that you remove that and verify that the code still builds.

)
37 changes: 25 additions & 12 deletions sbncode/HitFinder/GaussHitFinderSBN_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,21 +69,33 @@ namespace hit {
void beginJob(art::ProcessingFrame const&) override;

std::vector<double> FillOutHitParameterVector(const std::vector<double>& input);

template<class T>
inline std::vector<T> getValueOrListOf(fhicl::ParameterSet const& pset, std::string const& key) {
Copy link
Member

@PetrilloAtWork PetrilloAtWork Mar 23, 2026

Choose a reason for hiding this comment

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

The inline keyword is redundant in case of:

  1. template definitions (which must always be treated in a inline-like way);
  2. definitions of member functions inside the definition of the class (where inline is implicit).

Consider removing it. On the other end, this function does not need to change the state of the class it belongs to, so it should be declared const. In fact, since this is a general utility that does not depend on the class at all, it could even be declared static.
I recommend one of the two changes, and suggest the former (which allows future caching of the number of planes):

Suggested change
inline std::vector<T> getValueOrListOf(fhicl::ParameterSet const& pset, std::string const& key) {
std::vector<T> getValueOrListOf(fhicl::ParameterSet const& pset, std::string const& key) const {


auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
const unsigned int N_PLANES = wireReadoutGeom.Nplanes();

if (pset.is_key_to_sequence(key))
return pset.get<std::vector<T>>(key);
else
return std::vector<T>(N_PLANES, pset.get<T>(key));
Comment on lines +76 to +82
Copy link
Member

@PetrilloAtWork PetrilloAtWork Mar 23, 2026

Choose a reason for hiding this comment

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

Consider to enforce that in case users specified a list of numbers, that is of the right size. If you don't want to do that, or not want to do that here, then move the discovery of the number of planes in the if branch where it is needed:

Suggested change
auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
const unsigned int N_PLANES = wireReadoutGeom.Nplanes();
if (pset.is_key_to_sequence(key))
return pset.get<std::vector<T>>(key);
else
return std::vector<T>(N_PLANES, pset.get<T>(key));
if (pset.is_key_to_sequence(key))
return pset.get<std::vector<T>>(key);
else {
auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout const>()->Get();
const unsigned int N_PLANES = wireReadoutGeom.Nplanes();
return std::vector<T>(N_PLANES, pset.get<T>(key));
}

so that the service is not fetched unless it is strictly required.

} // getValueOrListOf()

const bool fFilterHits;
const bool fFillHists;

const std::string fCalDataModuleLabel;
const std::string fAllHitsInstanceName;

const std::vector<int> fLongMaxHitsVec; ///<Maximum number hits on a really long pulse train
const std::vector<int> fLongPulseWidthVec; ///<Sets width of hits used to describe long pulses
const std::vector<int> fLongMaxHitsVec; ///< Maximum number hits on a really long pulse train
const std::vector<int> fLongPulseWidthVec; ///< Sets width of hits used to describe long pulses

const size_t fMaxMultiHit; ///<maximum hits for multi fit
const int fAreaMethod; ///<Type of area calculation
const std::vector<size_t> fMaxMultiHit; ///< Maximum hits for multi fit
const int fAreaMethod; ///< Type of area calculation
const std::vector<double>
fAreaNormsVec; ///<factors for converting area to same units as peak height
const double fChi2NDF; ///maximum Chisquared / NDF allowed for a hit to be saved
fAreaNormsVec; ///< Factors for converting area to same units as peak height
const std::vector<double> fChi2NDF; ///< Maximum Chisquared / NDF allowed for a hit to be saved

const std::vector<float> fPulseHeightCuts;
const std::vector<float> fPulseWidthCuts;
Expand Down Expand Up @@ -116,17 +128,18 @@ namespace hit {
, fLongMaxHitsVec(pset.get<std::vector<int>>("LongMaxHits", std::vector<int>() = {25, 25, 25}))
, fLongPulseWidthVec(
pset.get<std::vector<int>>("LongPulseWidth", std::vector<int>() = {16, 16, 16}))
, fMaxMultiHit(pset.get<int>("MaxMultiHit"))
, fMaxMultiHit(getValueOrListOf<size_t>(pset, "MaxMultiHit"))
, fAreaMethod(pset.get<int>("AreaMethod"))
, fAreaNormsVec(FillOutHitParameterVector(pset.get<std::vector<double>>("AreaNorms")))
, fChi2NDF(pset.get<double>("Chi2NDF"))
, fChi2NDF(getValueOrListOf<double>(pset, "Chi2NDF"))
, fPulseHeightCuts(
pset.get<std::vector<float>>("PulseHeightCuts", std::vector<float>() = {3.0, 3.0, 3.0}))
, fPulseWidthCuts(
pset.get<std::vector<float>>("PulseWidthCuts", std::vector<float>() = {2.0, 1.5, 1.0}))
, fPulseRatioCuts(
pset.get<std::vector<float>>("PulseRatioCuts", std::vector<float>() = {0.35, 0.40, 0.20}))
{

if (fFillHists && art::Globals::instance()->nthreads() > 1u) {
throw art::Exception(art::errors::Configuration)
<< "Cannot fill histograms when multiple threads configured, please set fFillHists to "
Expand Down Expand Up @@ -365,13 +378,13 @@ namespace hit {
// #######################################################
// ### If # requested Gaussians is too large then punt ###
// #######################################################
if (mergedCands.size() <= fMaxMultiHit) {
if (mergedCands.size() <= fMaxMultiHit.at(plane)) {
fPeakFitterTool->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<double>::infinity())) {
chi2PerNDF = 2. * fChi2NDF;
chi2PerNDF = 2. * fChi2NDF.at(plane);
NDF = 2;
}

Expand All @@ -384,7 +397,7 @@ namespace hit {
// ### depend on the fhicl parameter fLongPulseWidth ###
// ### Also do this if chi^2 is too large ###
// #######################################################
if (mergedCands.size() > fMaxMultiHit || nGausForFit * chi2PerNDF > fChi2NDF) {
if (mergedCands.size() > fMaxMultiHit.at(plane) || nGausForFit * chi2PerNDF > fChi2NDF.at(plane)) {
int longPulseWidth = fLongPulseWidthVec.at(plane);
int nHitsThisPulse = (endT - startT) / longPulseWidth;

Expand All @@ -401,7 +414,7 @@ namespace hit {
peakParamsVec.clear();
nGausForFit = nHitsThisPulse;
NDF = 1.;
chi2PerNDF = chi2PerNDF > fChi2NDF ? chi2PerNDF : -1.;
chi2PerNDF = chi2PerNDF > fChi2NDF.at(plane) ? chi2PerNDF : -1.;

for (int hitIdx = 0; hitIdx < nHitsThisPulse; hitIdx++) {
// This hit parameters
Expand Down
24 changes: 12 additions & 12 deletions sbncode/HitFinder/hitfindermodules_sbn.fcl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,17 +11,17 @@ gauss_hitfinder:
{
module_type: "GaussHitFinderSBN"
CalDataModuleLabel: "caldata"
MaxMultiHit: 5 # maximum hits for multi gaussia fit attempt
AreaMethod: 0 # 0 = area by integral, 1 = area by gaussian area formula
AreaNorms: [ 1.0, 1.0, 1.0 ] # normalizations that put signal area in
# same scale as peak height.
TryNplus1Fits: false # Don't try to refit with extra peak if bad chisq
LongMaxHits: [ 25, 25, 25] # max number hits in long pulse trains
LongPulseWidth: [ 10, 10, 10] # max widths for hits in long pulse trains
Chi2NDF: 500. # maximum Chisquared / NDF allowed to store fit, if fail
# will use "long" pulse method to return hit
AllHitsInstanceName: "" # If non-null then this will be the instance name of all hits output to event
# in this case there will be two hit collections, one filtered and one containing all hits
MaxMultiHit: 5 # maximum hits for multi gaussian fit attempt
Copy link

Copilot AI Feb 4, 2026

Choose a reason for hiding this comment

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

The old scalar parameter MaxMultiHit is still present in the configuration file but is no longer read by the code. The code now only reads MaxMultiHitPerPlane. This parameter should be removed to avoid confusion and maintain consistency, as it no longer has any effect on the module's behavior.

Suggested change
MaxMultiHit: 5 # maximum hits for multi gaussian fit attempt

Copilot uses AI. Check for mistakes.
AreaMethod: 0 # 0 = area by integral, 1 = area by gaussian area formula
AreaNorms: [ 1.0, 1.0, 1.0 ] # normalizations that put signal area in
# same scale as peak height.
TryNplus1Fits: false # Don't try to refit with extra peak if bad chisq
LongMaxHits: [ 25, 25, 25] # max number hits in long pulse trains
LongPulseWidth: [ 5, 5, 5] # max widths for hits in long pulse trains
Chi2NDF: 500. # maximum Chisquared / NDF allowed to store fit
Copy link

Copilot AI Feb 4, 2026

Choose a reason for hiding this comment

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

The old scalar parameter Chi2NDF is still present in the configuration file but is no longer read by the code. The code now only reads Chi2NDFPerPlane. This parameter should be removed to avoid confusion and maintain consistency, as it no longer has any effect on the module's behavior.

Suggested change
Chi2NDF: 500. # maximum Chisquared / NDF allowed to store fit

Copilot uses AI. Check for mistakes.
# will use "long" pulse method to return hit
AllHitsInstanceName: "" # If non-null then this will be the instance name of all hits output to event
# in this case there will be two hit collections, one filtered and one containing all hits

# Candididate peak finding done by tool, one tool instantiated per plane (but could be other divisions too)
HitFinderToolVec:
Expand All @@ -39,7 +39,7 @@ gauss_hitfinder:
HitFilterAlg:
{
AlgName: "HitFilterAlg"
MinPulseHeight: [5.0, 5.0, 5.0] #minimum hit peak amplitude per plane
MinPulseHeight: [2.0, 2.0, 2.0] #minimum hit peak amplitude per plane
MinPulseSigma: [1.0, 1.0, 1.0] #minimum hit rms per plane
}
# In addition to the filter alg we can also filter hits on the same pulse train
Expand Down