From a88b33a5c5cfc7eb255c5ef40eae189f61f95ee4 Mon Sep 17 00:00:00 2001 From: Mattia Sotgia Date: Mon, 11 May 2026 12:58:31 -0500 Subject: [PATCH 1/2] Added checks to avoid segfault failures on hit creation / develop branch --- sbncode/HitFinder/GaussHitFinderSBN_module.cc | 69 +++++++++++++------ sbncode/HitFinder/hitfindermodules_sbn.fcl | 26 +++---- 2 files changed, 61 insertions(+), 34 deletions(-) diff --git a/sbncode/HitFinder/GaussHitFinderSBN_module.cc b/sbncode/HitFinder/GaussHitFinderSBN_module.cc index c0a91cc91..e6a3220c2 100644 --- a/sbncode/HitFinder/GaussHitFinderSBN_module.cc +++ b/sbncode/HitFinder/GaussHitFinderSBN_module.cc @@ -69,6 +69,18 @@ namespace hit { void beginJob(art::ProcessingFrame const&) override; std::vector FillOutHitParameterVector(const std::vector& input); + + template + inline std::vector getValueOrListOf(fhicl::ParameterSet const& pset, std::string const& key) { + + auto const& wireReadoutGeom = art::ServiceHandle()->Get(); + const unsigned int N_PLANES = wireReadoutGeom.Nplanes(); + + if (pset.is_key_to_sequence(key)) + return pset.get>(key); + else + return std::vector(N_PLANES, pset.get(key)); + } // getValueOrListOf() const bool fFilterHits; const bool fFillHists; @@ -76,14 +88,14 @@ namespace hit { const std::string fCalDataModuleLabel; const std::string fAllHitsInstanceName; - const std::vector fLongMaxHitsVec; /// fLongPulseWidthVec; /// fLongMaxHitsVec; ///< Maximum number hits on a really long pulse train + const std::vector fLongPulseWidthVec; ///< Sets width of hits used to describe long pulses - const size_t fMaxMultiHit; /// fMaxMultiHit; ///< Maximum hits for multi fit + const int fAreaMethod; ///< Type of area calculation const std::vector - fAreaNormsVec; /// fChi2NDF; ///< Maximum Chisquared / NDF allowed for a hit to be saved const std::vector fPulseHeightCuts; const std::vector fPulseWidthCuts; @@ -116,10 +128,10 @@ namespace hit { , fLongMaxHitsVec(pset.get>("LongMaxHits", std::vector() = {25, 25, 25})) , fLongPulseWidthVec( pset.get>("LongPulseWidth", std::vector() = {16, 16, 16})) - , fMaxMultiHit(pset.get("MaxMultiHit")) + , fMaxMultiHit(getValueOrListOf(pset, "MaxMultiHit")) , fAreaMethod(pset.get("AreaMethod")) , fAreaNormsVec(FillOutHitParameterVector(pset.get>("AreaNorms"))) - , fChi2NDF(pset.get("Chi2NDF")) + , fChi2NDF(getValueOrListOf(pset, "Chi2NDF")) , fPulseHeightCuts( pset.get>("PulseHeightCuts", std::vector() = {3.0, 3.0, 3.0})) , fPulseWidthCuts( @@ -365,13 +377,19 @@ 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::infinity())) { - chi2PerNDF = 2. * fChi2NDF; + chi2PerNDF = 2. * fChi2NDF.at(plane); + NDF = 2; + } + + // If the chi2 is also low (low-er than -inf) we bail + if (!(chi2PerNDF > -std::numeric_limits::infinity())) { + chi2PerNDF = 2. * fChi2NDF.at(plane); NDF = 2; } @@ -384,7 +402,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; @@ -401,7 +419,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 @@ -445,14 +463,13 @@ namespace hit { float nsigmaADC(2.0); float newright(0); float newleft(0); + for (const auto& 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 "<= endT)) continue; + // Extract errors float peakAmpErr = peakParams.peakAmplitudeError; float peakMeanErr = peakParams.peakCenterError; @@ -531,6 +550,9 @@ namespace hit { if (HitsumEndItr > sumEndItr) HitsumEndItr = sumEndItr; + // This prevents L577 and L579 to make any possible boundaty flip + if (HitsumStartItr > HitsumEndItr) continue; + // ### Sum of ADC counts float ROIsumADC = std::accumulate(sumStartItr, sumEndItr, 0.); float HitsumADC = std::accumulate(HitsumStartItr, HitsumEndItr, 0.); @@ -569,6 +591,11 @@ namespace hit { numHits++; } // <---End loop over gaussians + // THIS IS NOT USED + // THIS IS NOT USED + // THIS IS NOT USED + // THIS IS NOT USED + // Should we filter hits? if (filteredHitCol && !filteredHitVec.empty()) { // ####################################################################### @@ -627,11 +654,11 @@ namespace hit { // if (fFillHists) fChi2->Fill(chi2PerNDF); } - } //<---End loop over merged candidate hits - } //<---End looping over ROI's - ); //end tbb parallel for - } //<---End looping over all the wires - ); //end tbb parallel for + } //< End loop over merged candidate hits + } //< End looping over ROI's + ); //< End tbb::parallel_for(ROI) + } //< End looping over all the wires + ); //< End tbb::parallel_for(channelROI) for (size_t i = 0; i < hitstruct_vec.size(); i++) { allHitCol.emplace_back(hitstruct_vec[i].hit_tbb, hitstruct_vec[i].channelROI_tbb); diff --git a/sbncode/HitFinder/hitfindermodules_sbn.fcl b/sbncode/HitFinder/hitfindermodules_sbn.fcl index 868e5a2b8..a9f579d72 100644 --- a/sbncode/HitFinder/hitfindermodules_sbn.fcl +++ b/sbncode/HitFinder/hitfindermodules_sbn.fcl @@ -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 + 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 + #< 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: @@ -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 @@ -51,4 +51,4 @@ gauss_hitfinder: # Define sbn version of gaushit finder gausshit_sbn: @local::gauss_hitfinder -END_PROLOG +END_PROLOG \ No newline at end of file From 9766a63f56c3e0401a26bbaae1efeb4a0dc413c8 Mon Sep 17 00:00:00 2001 From: Mattia Sotgia Date: Fri, 3 Jul 2026 10:15:58 -0500 Subject: [PATCH 2/2] Replaced +/-std::infinity() checks with std::isfinite() --- sbncode/HitFinder/GaussHitFinderSBN_module.cc | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/sbncode/HitFinder/GaussHitFinderSBN_module.cc b/sbncode/HitFinder/GaussHitFinderSBN_module.cc index 16ee2f3c4..e02ee448a 100644 --- a/sbncode/HitFinder/GaussHitFinderSBN_module.cc +++ b/sbncode/HitFinder/GaussHitFinderSBN_module.cc @@ -388,14 +388,8 @@ namespace hit { 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::infinity())) { - chi2PerNDF = 2. * fChi2NDF.at(plane); - NDF = 2; - } - - // If the chi2 is also low (low-er than -inf) we bail - if (!(chi2PerNDF > -std::numeric_limits::infinity())) { + // If the chi2 is not finite then there is a real problem so we bail + if (!std::isfinite(chi2PerNDF)) { chi2PerNDF = 2. * fChi2NDF.at(plane); NDF = 2; }