fix: Gaussian histogram fit implementation in Core

Core/include/Acts/Utilities/Histogram.hpp

@@ -10,11 +10,17 @@
 
 #include "Acts/Utilities/RangeXD.hpp"
 
+#include <algorithm>
 #include <array>
+#include <cassert>
+#include <cmath>
+#include <concepts>
+#include <ranges>
 #include <string>
 #include <tuple>
 
 #include <boost/histogram.hpp>
+#include <boost/histogram/accumulators/weighted_sum.hpp>
 
 namespace Acts::Experimental {
 
@@ -36,18 +42,31 @@ using AxisVariant =
     boost::histogram::axis::variant<BoostVariableAxis, BoostRegularAxis,
                                     BoostLogAxis>;
 
-/// @brief Underlying Boost type for histograms
-using BoostHist = decltype(boost::histogram::make_histogram(
-    std::declval<std::vector<AxisVariant>>()));
-
 /// @brief Underlying Boost type for ProfileHistogram
 using BoostProfileHist = decltype(boost::histogram::make_profile(
     std::declval<std::vector<AxisVariant>>()));
 
+/// @brief Underlying Boost type for @c Histogram
+///
+/// Uses a weighted-sum accumulator so that every bin carries both a content
+/// and a variance. An unweighted @c fill() accumulates the count into both,
+/// giving the usual @f$ \sqrt{N} @f$ error; @c setBin() can instead write an
+/// arbitrary value/error pair, e.g. the result of a Gaussian fit.
+using BoostHist = decltype(boost::histogram::make_histogram_with(
+    boost::histogram::dense_storage<
+        boost::histogram::accumulators::weighted_sum<double>>(),
+    std::declval<std::vector<AxisVariant>>()));
+
 /// @brief Multi-dimensional histogram wrapper using boost::histogram for data collection
 ///
 /// This class wraps boost::histogram to provide a ROOT-independent histogram
-/// implementation with compile-time dimensionality.
+/// implementation with compile-time dimensionality. Every bin carries both a
+/// content and an error: @c fill() accumulates the usual @f$ \sqrt{N} @f$
+/// error, while @c setBin() can write an arbitrary value/error pair -- for
+/// example the mean and width extracted from a Gaussian fit to each slice of
+/// a residual histogram. This makes @c Histogram the ROOT-independent
+/// equivalent of both a @c TH1 filled through @c Fill() and one populated
+/// through @c SetBinContent / @c SetBinError.
 ///
 /// @tparam Dim Number of dimensions
 template <std::size_t Dim>
@@ -62,7 +81,10 @@ class Histogram {
             const std::array<AxisVariant, Dim>& axes)
       : m_name(std::move(name)),
         m_title(std::move(title)),
-        m_hist(boost::histogram::make_histogram(axes.begin(), axes.end())) {}
+        m_hist(boost::histogram::make_histogram_with(
+            boost::histogram::dense_storage<
+                boost::histogram::accumulators::weighted_sum<double>>(),
+            std::vector<AxisVariant>(axes.begin(), axes.end()))) {}
 
   /// Copy constructor
   /// @param other The other histogram to copy from
@@ -89,6 +111,62 @@ class Histogram {
     std::apply([this](auto... v) { m_hist(v...); }, std::tuple_cat(values));
   }
 
+  /// Set the content of a single bin, discarding whatever was there before
+  ///
+  /// @param indices Zero-based bin index per axis, excluding under-/overflow
+  /// @param content The content to store
+  /// @remark Indices must be in `[0, axis.size())` for every axis
+  /// @remark The bin's error is set to @f$ \sqrt{content} @f$, matching the
+  ///         error @c fill() would give a bin with that many entries. Use
+  ///         @c setBin to set an arbitrary value/error pair instead.
+  void setBinContent(const std::array<int, Dim>& indices, double content) {
+    std::apply(
+        [&](auto... i) {
+          m_hist.at(i...) =
+              boost::histogram::accumulators::weighted_sum<double>(content);
+        },
+        indices);
+  }
+
+  /// Set the content and error of a single bin, discarding whatever was
+  /// there before
+  ///
+  /// @param indices Zero-based bin index per axis, excluding under-/overflow
+  /// @param content The content to store
+  /// @param error The uncertainty on @p content
+  /// @remark Indices must be in `[0, axis.size())` for every axis
+  void setBin(const std::array<int, Dim>& indices, double content,
+              double error) {
+    std::apply(
+        [&](auto... i) {
+          m_hist.at(i...) =
+              boost::histogram::accumulators::weighted_sum<double>(
+                  content, error * error);
+        },
+        indices);
+  }
+
+  /// Get the content of a single bin
+  ///
+  /// @param indices Zero-based bin index per axis, excluding under-/overflow
+  /// @return The bin content
+  /// @remark Indices must be in `[0, axis.size())` for every axis
+  double binContent(const std::array<int, Dim>& indices) const {
+    return std::apply([&](auto... i) { return m_hist.at(i...).value(); },
+                      indices);
+  }
+
+  /// Get the error of a single bin
+  ///
+  /// @param indices Zero-based bin index per axis, excluding under-/overflow
+  /// @return The uncertainty on the bin content
+  /// @remark Indices must be in `[0, axis.size())` for every axis
+  double binError(const std::array<int, Dim>& indices) const {
+    return std::apply(
+        [&](auto... i) { return std::sqrt(m_hist.at(i...).variance()); },
+        indices);
+  }
+
   /// Get histogram name
   /// @return The histogram name
   const std::string& name() const { return m_name; }
@@ -211,8 +289,14 @@ class Efficiency {
              const std::array<AxisVariant, Dim>& axes)
       : m_name(std::move(name)),
         m_title(std::move(title)),
-        m_accepted(boost::histogram::make_histogram(axes.begin(), axes.end())),
-        m_total(boost::histogram::make_histogram(axes.begin(), axes.end())) {}
+        m_accepted(boost::histogram::make_histogram_with(
+            boost::histogram::dense_storage<
+                boost::histogram::accumulators::weighted_sum<double>>(),
+            std::vector<AxisVariant>(axes.begin(), axes.end()))),
+        m_total(boost::histogram::make_histogram_with(
+            boost::histogram::dense_storage<
+                boost::histogram::accumulators::weighted_sum<double>>(),
+            std::vector<AxisVariant>(axes.begin(), axes.end()))) {}
 
   /// Fill efficiency histogram
   ///
@@ -266,14 +350,90 @@ using Efficiency2 = Efficiency<2>;
 ///
 /// @param hist2d The 2D histogram to project
 /// @return A 1D histogram containing the projection
+/// @note Unlike ROOT's `TH2::ProjectionX`, the sum runs over the under- and
+///       overflow bins of the Y axis as well.
 Histogram1 projectionX(const Histogram2& hist2d);
 
 /// Project a 2D histogram onto the Y axis (axis 1)
 ///
 /// @param hist2d The 2D histogram to project
 /// @return A 1D histogram containing the projection
+/// @note Unlike ROOT's `TH2::ProjectionY`, the sum runs over the under- and
+///       overflow bins of the X axis as well.
 Histogram1 projectionY(const Histogram2& hist2d);
 
+namespace detail {
+
+/// Core of @c sliceLastAxis, taking the outer bin indices as an array rather
+/// than a parameter pack so it can also be called from generic code that
+/// already has them packed (see @c extractMeanWidthProfiles)
+template <std::size_t Dim>
+Histogram1 sliceLastAxis(const Histogram<Dim>& hist,
+                         const std::array<int, Dim - 1>& outerBins) {
+  const auto& lastAxis = hist.histogram().axis(Dim - 1);
+
+  assert(std::ranges::all_of(std::views::iota(std::size_t{0}, Dim - 1),
+                             [&](std::size_t d) {
+                               return outerBins[d] >= 0 &&
+                                      outerBins[d] <
+                                          hist.histogram().axis(d).size();
+                             }) &&
+         "outer bin index out of range");
+
+  std::string sliceName = hist.name() + "_slice";
+  for (const int bin : outerBins) {
+    sliceName += "_" + std::to_string(bin);
+  }
+
+  std::array<AxisVariant, 1> axes = {lastAxis};
+  Histogram1 slice(std::move(sliceName), hist.title(), axes);
+
+  for (int k = 0; k < lastAxis.size(); ++k) {
+    std::array<int, Dim> indices{};
+    std::ranges::copy(outerBins, indices.begin());
+    indices[Dim - 1] = k;
+    slice.setBin({k}, hist.binContent(indices), hist.binError(indices));
+  }
+
+  return slice;
+}
+
+}  // namespace detail
+
+/// Extract the distribution along the last axis at fixed bins of the others
+///
+/// Equivalent to ROOT's `TH2::ProjectionY(name, xBin + 1, xBin + 1)` for a 2D
+/// histogram, or `TH3::ProjectionZ(name, xBin + 1, xBin + 1, yBin + 1, yBin +
+/// 1)` for a 3D one.
+///
+/// @param hist The histogram to slice
+/// @param outerBins Zero-based bin index for every axis but the last, in
+///                  axis order
+/// @return A 1D histogram over the last axis
+/// @remark Every entry of @p outerBins must be in range for its axis
+template <std::size_t Dim, std::integral... Ints>
+  requires(sizeof...(Ints) + 1 == Dim)
+Histogram1 sliceLastAxis(const Histogram<Dim>& hist, Ints... outerBins) {
+  return detail::sliceLastAxis<Dim>(
+      hist, std::array<int, Dim - 1>{static_cast<int>(outerBins)...});
+}
+
+/// Total content of a histogram's in-range bins
+///
+/// The ROOT-independent equivalent of `TH1::GetEntries()` on a filled
+/// histogram.
+///
+/// @param hist The histogram to sum
+/// @return The sum of all in-range bin contents
+template <std::size_t Dim>
+double totalContent(const Histogram<Dim>& hist) {
+  double total = 0;
+  for (auto&& bin : boost::histogram::indexed(hist.histogram())) {
+    total += (*bin).value();
+  }
+  return total;
+}
+
 /// Extract bin edges from an AxisVariant
 ///
 /// Works with all axis types (regular, variable, log) in the variant by

Core/include/Acts/Utilities/HistogramFit.hpp

@@ -0,0 +1,40 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 CERN for the benefit of the ACTS project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Utilities/Histogram.hpp"
+
+#include <functional>
+#include <optional>
+#include <tuple>
+#include <utility>
+
+namespace Acts::Experimental {
+
+/// @brief Outcome of a Gaussian fit to a 1D histogram: `(mean, sigma,
+///        meanError, sigmaError)`
+///
+/// A plain @c std::tuple rather than a named struct so that a Python fit
+/// backend can return an ordinary 4-tuple and have it convert automatically
+/// via `pybind11/stl.h`, with no dedicated binding required.
+using HistogramFitResult = std::tuple<double, double, double, double>;
+
+/// @brief Fit range `[xMin, xMax]`, closed, selected by bin centre
+using HistogramFitRange = std::pair<double, double>;
+
+/// A single Gaussian fit to a 1D histogram, optionally restricted to a range
+///
+/// Any backend -- @c ActsExamples::gaussianHistogramFit, a callable
+/// `ActsPlugins::RootHistogramFit`, or a Python callable -- can be adapted to
+/// this signature. Living in Core rather than Examples or a plugin lets both
+/// sides share the same vocabulary types without depending on each other.
+using HistogramFitFunction = std::function<std::optional<HistogramFitResult>(
+    const Histogram1&, std::optional<HistogramFitRange>)>;
+
+}  // namespace Acts::Experimental

Core/src/Utilities/Histogram.cpp

@@ -8,32 +8,47 @@
 
 #include "Acts/Utilities/Histogram.hpp"
 
+#include <array>
 #include <cassert>
+#include <string>
 #include <vector>
 
 namespace Acts::Experimental {
 
+namespace {
+
+/// Wrap an already projected boost histogram into a Histogram1, carrying over
+/// both the axis (including its metadata) and the bin contents.
+Histogram1 wrapProjection(const BoostHist& projected, std::string name,
+                          std::string title) {
+  std::array<AxisVariant, 1> axes = {projected.axis(0)};
+  Histogram1 result(std::move(name), std::move(title), axes);
+
+  for (int i = 0; i < projected.axis(0).size(); ++i) {
+    const auto& bin = projected.at(i);
+    result.setBin({i}, bin.value(), std::sqrt(bin.variance()));
+  }
+
+  return result;
+}
+
+}  // namespace
+
 // Projection free functions
 Histogram1 projectionX(const Histogram2& hist2d) {
-  auto projectedHist = boost::histogram::algorithm::project(
+  const BoostHist projectedHist = boost::histogram::algorithm::project(
       hist2d.histogram(), std::integral_constant<unsigned, 0>{});
 
-  // Extract single axis from projected histogram
-  std::array<AxisVariant, 1> axes = {projectedHist.axis(0)};
-
-  return Histogram1(hist2d.name() + "_projX", hist2d.title() + " projection X",
-                    axes);
+  return wrapProjection(projectedHist, hist2d.name() + "_projX",
+                        hist2d.title() + " projection X");
 }
 
 Histogram1 projectionY(const Histogram2& hist2d) {
-  auto projectedHist = boost::histogram::algorithm::project(
+  const BoostHist projectedHist = boost::histogram::algorithm::project(
       hist2d.histogram(), std::integral_constant<unsigned, 1>{});
 
-  // Extract single axis from projected histogram
-  std::array<AxisVariant, 1> axes = {projectedHist.axis(0)};
-
-  return Histogram1(hist2d.name() + "_projY", hist2d.title() + " projection Y",
-                    axes);
+  return wrapProjection(projectedHist, hist2d.name() + "_projY",
+                        hist2d.title() + " projection Y");
 }
 
 std::vector<double> extractBinEdges(const AxisVariant& axis) {

Examples/Framework/CMakeLists.txt

@@ -24,6 +24,8 @@ acts_add_library(
     src/Validation/DuplicationPlotTool.cpp
     src/Validation/EffPlotTool.cpp
     src/Validation/FakePlotTool.cpp
+    src/Validation/GaussianHistogramFit.cpp
+    src/Validation/HistogramFit.cpp
     src/Validation/ResPlotTool.cpp
     src/Validation/TrackClassification.cpp
     src/Validation/PatternRecognitionPerformanceCollector.cpp

Examples/Framework/include/ActsExamples/Validation/GaussianHistogramFit.hpp

@@ -0,0 +1,41 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 CERN for the benefit of the ACTS project
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Utilities/Histogram.hpp"
+#include "ActsExamples/Validation/HistogramFit.hpp"
+
+#include <optional>
+
+namespace ActsExamples {
+
+/// Fit a Gaussian to a 1D histogram by chi-square minimisation
+///
+/// The model is @f$ \mu_i = A \exp(-(x_i - m)^2 / (2 s^2)) @f$, @f$ x_i @f$
+/// the bin centre, compared directly against the bin content: no bin-width
+/// factor and no integration of the model over the bin. Minimises
+/// @f$ \sum_i (n_i - \mu_i)^2 / n_i @f$ over bins with @f$ n_i > 0 @f$ --
+/// Neyman's chi-square with Poisson counting errors, matching ROOT's
+/// predefined `"gaus"` under `TH1::Fit(..., "SQ0")`, which gives zero-content
+/// bins zero error and drops them from the sum. See
+/// @c ActsPlugins::RootHistogramFit for a ROOT-backed fit with the same
+/// interface.
+///
+/// @param hist The histogram to fit
+/// @param range If set, only bins whose centre lies in `[range->first,
+///              range->second]` enter the fit; matches how ROOT restricts a
+///              fit range. If unset, the fit uses every bin.
+/// @return `(mean, sigma, meanError, sigmaError)`, or @c std::nullopt if the
+///         fit could not be performed
+/// @note Under- and overflow bins are always ignored.
+std::optional<HistogramFitResult> gaussianHistogramFit(
+    const Acts::Experimental::Histogram1& hist,
+    std::optional<HistogramFitRange> range = std::nullopt);
+
+}  // namespace ActsExamples