PulseFitter.cxx

#include "Tools/PulseFitter.h"
 
#include <cmath>
#include <limits>
 
namespace ldmx {
 
PulseFitter::PulseFitter(const PulseShape& shape, double t0_offset_ns,
                         double sampling_interval_ns, double pedestal_adc,
                         double noise_sigma_adc, double t_scan_min_ns,
                         double t_scan_max_ns, double t_scan_step_ns)
    : shape_(&shape),
      t0_offset_ns_(t0_offset_ns),
      sampling_interval_ns_(sampling_interval_ns),
      pedestal_adc_(pedestal_adc),
      inv_sigma2_(1.0 / (noise_sigma_adc * noise_sigma_adc)),
      t_scan_min_ns_(t_scan_min_ns),
      t_scan_max_ns_(t_scan_max_ns),
      t_scan_step_ns_(t_scan_step_ns) {}
 
// ---------------------------------------------------------------------------
// Private helper: chi2 and amplitude for a fixed hit time T
// ---------------------------------------------------------------------------
 
std::pair<double, double> PulseFitter::evalAtT(
    const std::vector<short>& samples, double T) const {
  const int n = static_cast<int>(samples.size());
 
  // Evaluate pulse shape at each sample time.
  // t_i = t0_offset + i * dt  →  argument to shape = t_i - T
  double sum_sf = 0.0;  // Σ (S_i - ped) * f_i
  double sum_f2 = 0.0;  // Σ f_i²
 
  for (int i = 0; i < n; ++i) {
    const double t_arg = t0_offset_ns_ + i * sampling_interval_ns_ - T;
    const double fi = shape_->eval(t_arg);
    const double si = static_cast<double>(samples[i]) - pedestal_adc_;
    sum_sf += si * fi;
    sum_f2 += fi * fi;
  }
 
  // Degenerate: pulse shape is zero everywhere in the window.
  if (sum_f2 < 1.0e-12) {
    return {std::numeric_limits<double>::max(), 0.0};
  }
 
  const double amplitude = sum_sf / sum_f2;
 
  // Compute chi2 for this (T, amplitude) pair.
  double chi2 = 0.0;
  for (int i = 0; i < n; ++i) {
    const double t_arg = t0_offset_ns_ + i * sampling_interval_ns_ - T;
    const double fi = shape_->eval(t_arg);
    const double si = static_cast<double>(samples[i]) - pedestal_adc_;
    const double resid = si - amplitude * fi;
    chi2 += resid * resid * inv_sigma2_;
  }
 
  return {chi2, amplitude};
}
 
// ---------------------------------------------------------------------------
// Public: fit
// ---------------------------------------------------------------------------
 
PulseFitter::FitResult PulseFitter::fit(
    const std::vector<short>& samples) const {
  FitResult result;
  result.ndf = static_cast<int>(samples.size()) - 2;
 
  // Quick sanity: check that at least one sample is above pedestal.
  bool any_signal = false;
  for (auto s : samples) {
    if (static_cast<double>(s) > pedestal_adc_) {
      any_signal = true;
      break;
    }
  }
  if (!any_signal) {
    return result;  // converged = false
  }
 
  // -------------------------------------------------------------------
  // Coarse scan over T to find the approximate minimum.
  // -------------------------------------------------------------------
  double best_chi2 = std::numeric_limits<double>::max();
  double best_t = t_scan_min_ns_;
  double best_amp = 0.0;
 
  // Keep track of the three-point neighbourhood for parabolic refinement.
  double t_lo = t_scan_min_ns_, chi2_lo = std::numeric_limits<double>::max();
  double t_hi = t_scan_min_ns_, chi2_hi = std::numeric_limits<double>::max();
 
  double t_prev = t_scan_min_ns_;
  double chi2_prev = std::numeric_limits<double>::max();
 
  for (double t = t_scan_min_ns_; t <= t_scan_max_ns_; t += t_scan_step_ns_) {
    auto [chi2, amp] = evalAtT(samples, t);
 
    if (chi2 < best_chi2) {
      // Update neighbourhood: the previous point becomes T_lo.
      t_lo = t_prev;
      chi2_lo = chi2_prev;
      best_chi2 = chi2;
      best_t = t;
      best_amp = amp;
      // T_hi is the next point, forget one found for an earlier minimum
      chi2_hi = std::numeric_limits<double>::max();
      t_hi = t;
    } else if (t_prev == best_t) {
      // First point to the right of the minimum.
      t_hi = t;
      chi2_hi = chi2;
    }
 
    t_prev = t;
    chi2_prev = chi2;
  }
 
  // -------------------------------------------------------------------
  // Parabolic refinement around the minimum if neighbours are available.
  // -------------------------------------------------------------------
  if (t_lo < best_t && t_hi > best_t && chi2_lo > best_chi2 &&
      chi2_hi > best_chi2) {
    // Fit a parabola through (T_lo, chi2_lo), (best_T, best_chi2), (T_hi,
    // chi2_hi).
    const double d1 = best_t - t_lo;
    const double d2 = t_hi - best_t;
    const double dc1 = chi2_lo - best_chi2;
    const double dc2 = chi2_hi - best_chi2;
    // Vertex of parabola: ΔT = (d2²·dc1 - d1²·dc2) / (2·(d1·dc2 + d2·dc1))
    const double denom = 2.0 * (d1 * dc2 + d2 * dc1);
    if (std::abs(denom) > 1.0e-12) {
      const double delta = (d2 * d2 * dc1 - d1 * d1 * dc2) / denom;
      // Only accept the refinement if it stays within the bracket.
      if (std::abs(delta) < std::max(d1, d2)) {
        const double t_refined = best_t + delta;
        auto [chi2_ref, amp_ref] = evalAtT(samples, t_refined);
        if (chi2_ref < best_chi2) {
          best_t = t_refined;
          best_amp = amp_ref;
          best_chi2 = chi2_ref;
        }
      }
    }
  }
 
  result.amplitude = best_amp;
  result.t0 = best_t;
  result.chi2 = best_chi2;
  result.converged = true;
  return result;
}
 
// ---------------------------------------------------------------------------
// Public: fitFixed
// ---------------------------------------------------------------------------
 
PulseFitter::FitResult PulseFitter::fitFixed(const std::vector<short>& samples,
                                             double T) const {
  FitResult result;
  result.ndf = static_cast<int>(samples.size()) - 1;
  auto [chi2, amp] = evalAtT(samples, T);
  result.amplitude = amp;
  result.t0 = T;
  result.chi2 = chi2;
  result.converged = chi2 < std::numeric_limits<double>::max();
  return result;
}
 
}  // namespace ldmx