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 6 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 4 commits
Commits
Show all changes
6 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
c03a60b
address PR comments round 1
Apr 10, 2026
0aaeca7
Merge pull request #2 from carriganm95/update_hitTuning
carriganm95 Apr 10, 2026
90384ac
fix brackets
Apr 10, 2026
cf2f906
Merge pull request #3 from carriganm95/update_hitTuning
carriganm95 Apr 10, 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
43 changes: 31 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,39 @@ namespace hit {
void beginJob(art::ProcessingFrame const&) override;

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

template<class T>
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))
auto vec = pset.get<std::vector<T>>(key);
if (vec.size() != N_PLANES) {
throw cet::exception("Configuration")
<< "Parameter '" << key << "' has size " << vec.size()
Comment thread
carriganm95 marked this conversation as resolved.
Outdated
<< " but expected " << N_PLANES << " (number of planes)";
}
return pset.get<std::vector<T>>(key);
else
return std::vector<T>(N_PLANES, pset.get<T>(key));
Comment on lines +76 to +88
Copy link
Copy Markdown
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.

Copy link
Copy Markdown
Author

@carriganm95 carriganm95 Apr 10, 2026

Choose a reason for hiding this comment

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

see previous comment, number of entries in a list is now enforced to be equal to planes

} // 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 +134,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 +384,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 +403,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 +420,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
Comment thread
carriganm95 marked this conversation as resolved.
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
Comment thread
carriganm95 marked this conversation as resolved.
# 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
Loading