Skip to content

functional fit option for waveform processing#4131

Merged
pinkenburg merged 5 commits intosPHENIX-Collaboration:masterfrom
Shuonli:function_fit
Jan 21, 2026
Merged

functional fit option for waveform processing#4131
pinkenburg merged 5 commits intosPHENIX-Collaboration:masterfrom
Shuonli:function_fit

Conversation

@Shuonli
Copy link
Copy Markdown
Contributor

@Shuonli Shuonli commented Jan 20, 2026

Types of changes

  • Bug fix (non-breaking change which fixes an issue)
  • New feature (non-breaking change which adds functionality)
  • Breaking change (fix or feature that would cause existing functionality to not work for users)
  • Requiring change in macros repository (Please provide links to the macros pull request in the last section)
  • I am a member of GitHub organization of sPHENIX Collaboration, EIC, or ECCE (contact Chris Pinkenburg to join)

What kind of change does this PR introduce? (Bug fix, feature, ...)

added an option to use functional fit(CaloWaveformProcessing::FUNCFIT) in calotowerbuilder.

The functional form is copied from the testbeam processing code.

TODOs (if applicable)

Links to other PRs in macros and calibration repositories (if applicable)

Summary

Motivation / Context

This PR adds a functional-fit option to the calorimeter waveform processing chain to provide an alternative reconstruction path based on the functional forms used in sPHENIX testbeam processing. The new mode enables per-channel parametric fits in CaloWaveformProcessing/CaloWaveformFitting (and wiring in CaloTowerBuilder), which may improve amplitude and timing estimates for pulse shapes not well-modelled by existing TEMPLATE/NYQUIST/ONNX/FAST methods.

Key changes

  • New processing mode: CaloWaveformProcessing::FUNCFIT (enum added).
  • Two fit models implemented in CaloWaveformFitting:
    • PowerLawExp (5 parameters)
    • PowerLawDoubleExp (7 parameters)
  • New per-channel fitting routine calo_processing_funcfit(...) that:
    • Builds per-channel histograms, handles pedestal estimation, zero-suppression and saturation cases,
    • Runs ROOT TF1 fits for the selected model,
    • Returns per-channel metrics (amplitude, time, pedestal, chi2/ndf).
  • New configuration API and plumbing:
    • CaloWaveformFitting/CaloWaveformProcessing/CaloTowerBuilder setters:
      • set_funcfit_type(...)
      • set_powerlaw_params(power, decay)
      • set_doubleexp_params(power, peaktime1, peaktime2, ratio)
    • CaloTowerBuilder integrates configuration and selects FUNCFIT in InitRun when requested.
  • No on-disk format changes; public interface extended with new setters and enum value (backwards-compatible).

Potential risk areas

  • Reconstruction behavior: FUNCFIT will produce different amplitude/time estimates; calibration and downstream workflows may need retuning and validation.
  • Performance: Per-channel histogram creation and ROOT TF1 fits are substantially more CPU intensive than template/fast methods — may impact throughput and HLT use.
  • Fit robustness: Fits can fail or return unphysical parameters for saturated, noisy, or edge-case waveforms; ensure fallbacks and monitoring.
  • Thread-safety: ROOT histogram/TF1 usage may be problematic in multi-threaded environments unless ROOT is configured for thread-safety — verify in parallel workflows.
  • Defaults & tuning: Hard-coded defaults (from testbeam tuning) may not be optimal for all detectors/conditions; site-specific tuning required.
  • API surface: Public API expanded; source compatibility preserved but clients should be aware of new enum and setters.

Possible future improvements

  • Documentation and recommended per-detector parameter tuning.
  • Monitoring/diagnostics for fit failures, chi2 distributions, and comparisons vs existing methods.
  • Expose fitter diagnostics (parameter errors, covariance) and implement robust fallback strategies for failed fits.
  • Evaluate alternative minimizers/backends or faster approximations to reduce CPU cost.
  • Add regression/validation macros comparing FUNCFIT to existing processing on calibration and physics samples.

Note: This summary was AI-assisted — please verify fit implementation details, default parameter choices, and zero-suppression/saturation handling against the original testbeam code and intended production usage. AI-generated content can contain mistakes; use best judgment when reviewing code and tests.

@coderabbitai
Copy link
Copy Markdown
Contributor

coderabbitai bot commented Jan 20, 2026

📝 Walkthrough

Walkthrough

Adds a functional-fit processing path (FUNCFIT): new PowerLawExp and PowerLawDoubleExp signal shapes, a per-channel fitting routine calo_processing_funcfit(), and public configuration APIs to select fit type and parameters across CaloWaveformFitting, CaloWaveformProcessing, and CaloTowerBuilder.

Changes

Cohort / File(s) Change Summary
Waveform fitting core
offline/packages/CaloReco/CaloWaveformFitting.h, offline/packages/CaloReco/CaloWaveformFitting.cc
Adds FuncFitType enum, SignalShape_PowerLawExp and SignalShape_PowerLawDoubleExp, and calo_processing_funcfit() which builds per-channel histograms and performs parametric fits (PowerLawExp or PowerLawDoubleExp) using ROOT TF1/TH1F. Adds setters and backing members for fit parameters.
Waveform processing integration
offline/packages/CaloReco/CaloWaveformProcessing.h, offline/packages/CaloReco/CaloWaveformProcessing.cc
Adds FUNCFIT = 6 to processing enum. initialize_processing() configures CaloWaveformFitting for FUNCFIT (sets fit type and parameters); process_waveform() routes to calo_processing_funcfit() when FUNCFIT is selected.
Tower builder configuration
offline/packages/CaloReco/CaloTowerBuilder.h, offline/packages/CaloReco/CaloTowerBuilder.cc
Adds public setters set_funcfit_type(int), set_powerlaw_params(double,double), set_doubleexp_params(double,double,double,double) and internal state. InitRun now configures WaveformProcessing for FUNCFIT via these setters. (CaloTowerBuilder.cc also has a minor whitespace change.)
Header / API surface
offline/packages/CaloReco/CaloWaveformFitting.h, offline/packages/CaloReco/CaloWaveformProcessing.h, offline/packages/CaloReco/CaloTowerBuilder.h
New public methods and enum values to configure/select functional-fit behavior; no existing public signatures removed.

Sequence Diagram(s)

sequenceDiagram
    participant Tower as CaloTowerBuilder
    participant Proc as CaloWaveformProcessing
    participant Fit as CaloWaveformFitting
    participant ROOT as ROOT/TF1

    Tower->>Proc: initialize_processing(FUNCFIT) + set fit params
    Proc->>Fit: set_funcfit_type / set_powerlaw_params / set_doubleexp_params
    Tower->>Proc: process_event -> process_waveform()
    Proc->>Fit: calo_processing_funcfit(chnlvector)
    Fit->>Fit: build per-channel histograms, handle ZS/saturation
    Fit->>ROOT: fit histogram (PowerLawExp or PowerLawDoubleExp)
    ROOT-->>Fit: fit results (amplitude, time, pedestal, chi2)
    Fit-->>Proc: per-channel fit vectors
    Proc-->>Tower: processed waveform / tower energies
Loading
✨ Finishing touches
  • 📝 Generate docstrings

Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out.

❤️ Share

Comment @coderabbitai help to get the list of available commands and usage tips.

@Shuonli Shuonli changed the title add functional fit functional fit option for waveform processing Jan 20, 2026
Copy link
Copy Markdown
Contributor

@coderabbitai coderabbitai bot left a comment

Choose a reason for hiding this comment

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

Actionable comments posted: 3

std::string m_zs_fieldname{"zs_threshold"};

// Functional fit parameters
int m_funcfit_type{1}; // 0 = PowerLawExp, 1 = PowerLawDoubleExp
Copy link
Copy Markdown
Contributor

@coderabbitai coderabbitai bot Jan 20, 2026

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

Inconsistent default for m_funcfit_type.

CaloTowerBuilder defaults m_funcfit_type to 1 (POWERLAWDOUBLEEXP), while CaloWaveformFitting and CaloWaveformProcessing default to 0 (POWERLAWEXP). This inconsistency could cause unexpected behavior when users rely on defaults—the fit type used will depend on whether configuration is done through CaloTowerBuilder or directly on the lower-level classes.

Consider aligning all defaults to the same value (likely 0 to match the physics fitting classes).

Suggested fix
   // Functional fit parameters
-  int m_funcfit_type{1};  // 0 = PowerLawExp, 1 = PowerLawDoubleExp
+  int m_funcfit_type{0};  // 0 = PowerLawExp, 1 = PowerLawDoubleExp
📝 Committable suggestion

‼️ IMPORTANT
Carefully review the code before committing. Ensure that it accurately replaces the highlighted code, contains no missing lines, and has no issues with indentation. Thoroughly test & benchmark the code to ensure it meets the requirements.

Suggested change
int m_funcfit_type{1}; // 0 = PowerLawExp, 1 = PowerLawDoubleExp
int m_funcfit_type{0}; // 0 = PowerLawExp, 1 = PowerLawDoubleExp

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

The Rabbit has a point here, the default for the fit type should be identical across classes, otherwise this will confuse everybody. What to use as default is up to you

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

@pinkenburg Thank you for confirming. @Shuonli, could you please decide which default fit type should be used consistently across CaloTowerBuilder, CaloWaveformFitting, and CaloWaveformProcessing? Currently:

  • CaloTowerBuilder.h: defaults to 1 (POWERLAWDOUBLEEXP)
  • CaloWaveformFitting.h and CaloWaveformProcessing.h: default to 0 (POWERLAWEXP)

Once you decide, all three classes should use the same default value to avoid confusion.


🧠 Learnings used
Learnt from: silas-gross
Repo: sPHENIX-Collaboration/coresoftware PR: 4100
File: generators/Herwig/HepMCTrigger/HepMCParticleTrigger.h:93-108
Timestamp: 2026-01-10T07:54:57.021Z
Learning: In HepMCParticleTrigger (generators/Herwig/HepMCTrigger/HepMCParticleTrigger.h), the SetEta* and SetAbsEta* methods intentionally share the same storage variables (_theEtaHigh and _theEtaLow) to emulate the Pythia trigger behavior. This design pattern should be maintained for consistency with existing trigger implementations.

Co-authored-by: coderabbitai[bot] <136622811+coderabbitai[bot]@users.noreply.github.com>
Copy link
Copy Markdown
Contributor

@coderabbitai coderabbitai bot left a comment

Choose a reason for hiding this comment

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

Actionable comments posted: 1

♻️ Duplicate comments (1)
offline/packages/CaloReco/CaloWaveformFitting.cc (1)

717-729: Guard against out-of-bounds access when nsamples < 7.

Line 718 dereferences v.at(6) and Line 724 uses v.at(1) without checking bounds. This can throw for short waveforms (e.g., 3–6 samples). This is the same issue previously raised for the template-fit path and should be handled here too.

✅ Suggested defensive guard
-    if ((_bdosoftwarezerosuppression && v.at(6) - v.at(0) < _nsoftwarezerosuppression) ||
-        (_maxsoftwarezerosuppression && maxheight - pedestal < _nsoftwarezerosuppression))
+    if ((_bdosoftwarezerosuppression && nsamples > 6 && v.at(6) - v.at(0) < _nsoftwarezerosuppression) ||
+        (_maxsoftwarezerosuppression && maxheight - pedestal < _nsoftwarezerosuppression))
     {
-      amp = v.at(6) - v.at(0);
+      amp = (nsamples > 6) ? v.at(6) - v.at(0) : maxheight - pedestal;
       time = std::numeric_limits<float>::quiet_NaN();
       ped = v.at(0);
-      if (v.at(0) != 0 && v.at(1) == 0)
+      if (nsamples > 1 && v.at(0) != 0 && v.at(1) == 0)
       {
         chi2 = 1000000;
       }

Comment on lines +770 to +797
double risetime = m_powerlaw_power / m_powerlaw_decay;
double par[5];
par[0] = maxheight - pedestal; // Amplitude
par[1] = maxbin - risetime; // t0
if (par[1] < 0)
{
par[1] = 0;
}
par[2] = m_powerlaw_power; // Power
par[3] = m_powerlaw_decay; // Decay
par[4] = pedestal; // Pedestal

f.SetParameters(par);
f.SetParLimits(0, (maxheight - pedestal) * 0.5, (maxheight - pedestal) * 10);
f.SetParLimits(1, 0, nsamples);
f.SetParLimits(2, 0, 10.0);
f.SetParLimits(3, 0, 10.0);
f.SetParLimits(4, pedestal - std::abs(maxheight - pedestal), pedestal + std::abs(maxheight - pedestal));

// Perform fit
h.Fit(&f, "QRN0W", "", 0, nsamples);

// Calculate peak amplitude and time from fit parameters
// Peak height is (p0 * Power(p2/p3, p2)) / exp(p2)
fit_amp = (f.GetParameter(0) * pow(f.GetParameter(2) / f.GetParameter(3), f.GetParameter(2))) / exp(f.GetParameter(2));
// Peak time is t0 + power/decay
fit_time = f.GetParameter(1) + f.GetParameter(2) / f.GetParameter(3);
fit_ped = f.GetParameter(4);
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

⚠️ Potential issue | 🟠 Major

Prevent divide-by-zero in POWERLAWEXP risetime and peak calculation.

Line 770 divides by m_powerlaw_decay, and Lines 794–796 divide by the fitted decay. If decay is 0 (or allowed to fit to 0), this yields inf/NaN and contaminates fit outputs. Tighten the lower bound and guard the pre-fit risetime.

✅ Suggested fix
-      double risetime = m_powerlaw_power / m_powerlaw_decay;
+      const double decay0 = std::max(m_powerlaw_decay, 1e-6);
+      double risetime = m_powerlaw_power / decay0;
@@
-      f.SetParLimits(3, 0, 10.0);
+      f.SetParLimits(3, 1e-6, 10.0);
@@
-      fit_amp = (f.GetParameter(0) * pow(f.GetParameter(2) / f.GetParameter(3), f.GetParameter(2))) / exp(f.GetParameter(2));
-      fit_time = f.GetParameter(1) + f.GetParameter(2) / f.GetParameter(3);
+      const double decay = f.GetParameter(3);
+      fit_amp = (f.GetParameter(0) * pow(f.GetParameter(2) / decay, f.GetParameter(2))) / exp(f.GetParameter(2));
+      fit_time = f.GetParameter(1) + f.GetParameter(2) / decay;

@sphenix-jenkins-ci
Copy link
Copy Markdown

Build & test report

Report for commit 2f01ab562a041546fab06579c672cbacb2734568:
Jenkins passed


Automatically generated by sPHENIX Jenkins continuous integration
sPHENIX             jenkins.io

@sphenix-jenkins-ci
Copy link
Copy Markdown

Build & test report

Report for commit 5334735d2a18c44a0d1737216fcf1cdb1d36e4fb:
Jenkins passed


Automatically generated by sPHENIX Jenkins continuous integration
sPHENIX             jenkins.io

@sphenix-jenkins-ci
Copy link
Copy Markdown

Build & test report

Report for commit 0d8b5fc4fd281a852503bb9c20df765a061945b6:
Jenkins passed


Automatically generated by sPHENIX Jenkins continuous integration
sPHENIX             jenkins.io

@sphenix-jenkins-ci
Copy link
Copy Markdown

Build & test report

Report for commit 47d0f386183507d42c2bd37a2751c83b354d73b5:
Jenkins passed


Automatically generated by sPHENIX Jenkins continuous integration
sPHENIX             jenkins.io

@pinkenburg pinkenburg merged commit f704495 into sPHENIX-Collaboration:master Jan 21, 2026
22 checks passed
@coderabbitai coderabbitai bot mentioned this pull request Feb 2, 2026
5 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants