refactor: Read Bethe-Heitler-Approx from json

Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp

@@ -8,6 +8,8 @@
 
 #pragma once
 
+#include "Acts/Utilities/RangeXD.hpp"
+
 #include <cmath>
 #include <cstddef>
 #include <numbers>
@@ -78,7 +80,7 @@ class BetheHeitlerApprox {
 };
 
 /// This class approximates the Bethe-Heitler with only one component. This is
-/// mainly inside @ref AtlasBetheHeitlerApprox, but can also be used as the
+/// mainly inside @ref PolynomialBetheHeitlerApprox, but can also be used as the
 /// only component approximation (then probably for debugging)
 /// @ingroup material
 class BetheHeitlerApproxSingleCmp final : public BetheHeitlerApprox {
@@ -117,11 +119,8 @@ class BetheHeitlerApproxSingleCmp final : public BetheHeitlerApprox {
 /// mixture. To enable an approximation for continuous input variables, the
 /// weights, means and variances are internally parametrized as a Nth order
 /// polynomial.
-/// @todo This class is rather inflexible: It forces two data representations,
-/// making it a bit awkward to add a single parameterization. It would be good
-/// to generalize this at some point.
 /// @ingroup material
-class AtlasBetheHeitlerApprox : public BetheHeitlerApprox {
+class PolynomialBetheHeitlerApprox : public BetheHeitlerApprox {
  public:
   /// Polynomial coefficient sets for a Gaussian mixture component.
   struct PolyData {
@@ -136,6 +135,19 @@ class AtlasBetheHeitlerApprox : public BetheHeitlerApprox {
   /// Type alias for array of polynomial data for all components
   using Data = std::vector<PolyData>;
 
+  /// Single x/x0 range with its data and transformation flag
+  struct RangeData {
+    Range1D<double> range;
+    Data data;
+    bool transform = false;
+
+    RangeData() = default;
+    RangeData(Range1D<double> r, Data d, bool t)
+        : range(r), data(std::move(d)), transform(t) {}
+    RangeData(double l, double h, Data d, bool t)
+        : range(l, h), data(std::move(d)), transform(t) {}
+  };
+
   /// Loads a parameterization from a file according to the Atlas file
   /// description
   ///
@@ -148,13 +160,30 @@ class AtlasBetheHeitlerApprox : public BetheHeitlerApprox {
   /// @param clampToRange forwarded to constructor
   /// @param noChangeLimit forwarded to constructor
   /// @param singleGaussianLimit forwarded to constructor
-  /// @return AtlasBetheHeitlerApprox instance loaded from parameter files
-  static AtlasBetheHeitlerApprox loadFromFiles(
+  /// @return PolynomialBetheHeitlerApprox instance loaded from parameter files
+  /// @deprecated Use loadBetheHeitlerApproxFromJson instead
+  [[deprecated(
+      "loadFromFiles is deprecated. Use loadBetheHeitlerApproxFromJson "
+      "instead.")]]
+  static PolynomialBetheHeitlerApprox loadFromFiles(
       const std::string &low_parameters_path,
       const std::string &high_parameters_path, double lowLimit,
       double highLimit, bool clampToRange, double noChangeLimit,
       double singleGaussianLimit);
 
+  /// Construct the Bethe-Heitler approximation description with N ranges.
+  /// Each range has its own data and transformation flag.
+  /// The ranges will be sorted by minimum value and validated for
+  /// non-overlapping.
+  ///
+  /// @param ranges Vector of range data
+  /// @param clampToRange whether to clamp the input x/x0 to the allowed range
+  /// @param noChangeLimit limit below which no change is applied
+  /// @param singleGaussianLimit limit below which a single Gaussian is used
+  PolynomialBetheHeitlerApprox(std::vector<RangeData> ranges, bool clampToRange,
+                               double noChangeLimit,
+                               double singleGaussianLimit);
+
   /// Construct the Bethe-Heitler approximation description with two
   /// parameterizations, one for lower ranges, one for higher ranges.
   /// Is it assumed that the lower limit of the high-x/x0 data is equal
@@ -169,24 +198,25 @@ class AtlasBetheHeitlerApprox : public BetheHeitlerApprox {
   /// @param clampToRange whether to clamp the input x/x0 to the allowed range
   /// @param noChangeLimit limit below which no change is applied
   /// @param singleGaussianLimit limit below which a single Gaussian is used
-  AtlasBetheHeitlerApprox(const Data &lowData, const Data &highData,
-                          bool lowTransform, bool highTransform,
-                          double lowLimit, double highLimit, bool clampToRange,
-                          double noChangeLimit, double singleGaussianLimit)
-      : m_lowData(lowData),
-        m_highData(highData),
-        m_lowTransform(lowTransform),
-        m_highTransform(highTransform),
-        m_lowLimit(lowLimit),
-        m_highLimit(highLimit),
-        m_clampToRange(clampToRange),
-        m_noChangeLimit(noChangeLimit),
-        m_singleGaussianLimit(singleGaussianLimit) {}
+  [[deprecated("Use constructor taking std::vector<RangeData> instead")]]
+  PolynomialBetheHeitlerApprox(const Data &lowData, const Data &highData,
+                               bool lowTransform, bool highTransform,
+                               double lowLimit, double highLimit,
+                               bool clampToRange, double noChangeLimit,
+                               double singleGaussianLimit)
+      : PolynomialBetheHeitlerApprox(
+            {{lowLimit, highLimit, lowData, lowTransform},
+             {highLimit, highLimit * 2, highData, highTransform}},
+            clampToRange, noChangeLimit, singleGaussianLimit) {}
 
   /// Returns the number of components the returned mixture will have
   /// @return Number of components in the mixture
   std::size_t maxComponents() const override {
-    return std::max(m_lowData.size(), m_highData.size());
+    std::size_t max = 0;
+    for (const auto &range : m_ranges) {
+      max = std::max(max, range.data.size());
+    }
+    return max;
   }
 
   /// Checks if an input is valid for the parameterization
@@ -197,7 +227,7 @@ class AtlasBetheHeitlerApprox : public BetheHeitlerApprox {
     if (m_clampToRange) {
       return true;
     } else {
-      return xOverX0 < m_highLimit;
+      return xOverX0 < m_ranges.back().range.max();
     }
   }
 
@@ -211,25 +241,22 @@ class AtlasBetheHeitlerApprox : public BetheHeitlerApprox {
       double xOverX0, const std::span<Component> mixture) const override;
 
  private:
-  Data m_lowData;
-  Data m_highData;
-  bool m_lowTransform = false;
-  bool m_highTransform = false;
-  double m_lowLimit = 0;
-  double m_highLimit = 0;
+  std::vector<RangeData> m_ranges;
   bool m_clampToRange = false;
   double m_noChangeLimit = 0;
   double m_singleGaussianLimit = 0;
 };
 
-/// Creates a @ref AtlasBetheHeitlerApprox object based on an ATLAS
-/// configuration, that are stored as static data in the source code.
-/// This may not be an optimal configuration, but should allow to run
-/// the GSF without the need to load files
-/// @param clampToRange Whether to clamp values to the valid range
-/// @return AtlasBetheHeitlerApprox with default ATLAS configuration parameters
-/// @ingroup material
-AtlasBetheHeitlerApprox makeDefaultBetheHeitlerApprox(
+/// @deprecated Use PolynomialBetheHeitlerApprox instead
+using AtlasBetheHeitlerApprox = PolynomialBetheHeitlerApprox;
+
+/// @deprecated Use PolynomialBetheHeitlerApprox instead
+/// @note This function is deprecated for documentation purposes. It still
+///       returns PolynomialBetheHeitlerApprox which is the new default.
+[[deprecated(
+    "AtlasBetheHeitlerApprox is deprecated. Use PolynomialBetheHeitlerApprox "
+    "instead.")]]
+inline PolynomialBetheHeitlerApprox makeDefaultBetheHeitlerApprox(
     bool clampToRange = false);
 
 /// @}

Core/src/TrackFitting/BetheHeitlerApprox.cpp

@@ -8,14 +8,41 @@
 
 #include "Acts/TrackFitting/BetheHeitlerApprox.hpp"
 
+#include "Acts/Utilities/RangeXD.hpp"
+
 #include <algorithm>
 #include <fstream>
 #include <stdexcept>
 #include <tuple>
 
 namespace Acts {
 
-AtlasBetheHeitlerApprox AtlasBetheHeitlerApprox::loadFromFiles(
+PolynomialBetheHeitlerApprox::PolynomialBetheHeitlerApprox(
+    std::vector<RangeData> ranges, bool clampToRange, double noChangeLimit,
+    double singleGaussianLimit)
+    : m_ranges(std::move(ranges)),
+      m_clampToRange(clampToRange),
+      m_noChangeLimit(noChangeLimit),
+      m_singleGaussianLimit(singleGaussianLimit) {
+  if (m_ranges.empty()) {
+    throw std::invalid_argument("At least one range is required");
+  }
+
+  // Sort ranges by minimum value
+  std::sort(m_ranges.begin(), m_ranges.end(), [](const auto &a, const auto &b) {
+    return a.range.min() < b.range.min();
+  });
+
+  // Validate that ranges don't overlap
+  for (std::size_t i = 1; i < m_ranges.size(); ++i) {
+    if (m_ranges[i - 1].range && m_ranges[i].range) {
+      throw std::invalid_argument(
+          "Overlapping ranges detected. Ranges must be non-overlapping.");
+    }
+  }
+}
+
+PolynomialBetheHeitlerApprox PolynomialBetheHeitlerApprox::loadFromFiles(
     const std::string &low_parameters_path,
     const std::string &high_parameters_path, double lowLimit, double highLimit,
     bool clampToRange, double noChangeLimit, double singleGaussianLimit) {
@@ -34,7 +61,12 @@ AtlasBetheHeitlerApprox AtlasBetheHeitlerApprox::loadFromFiles(
 
     Data data;
 
-    for (auto &cmp : data) {
+    for (std::size_t i = 0; i < n_cmps; ++i) {
+      PolyData cmp;
+      cmp.weightCoeffs.resize(degree + 1);
+      cmp.meanCoeffs.resize(degree + 1);
+      cmp.varCoeffs.resize(degree + 1);
+
       for (auto &coeff : cmp.weightCoeffs) {
         file >> coeff;
       }
@@ -44,6 +76,8 @@ AtlasBetheHeitlerApprox AtlasBetheHeitlerApprox::loadFromFiles(
       for (auto &coeff : cmp.varCoeffs) {
         file >> coeff;
       }
+
+      data.push_back(std::move(cmp));
     }
 
     return std::make_tuple(data, transform_code);
@@ -52,14 +86,19 @@ AtlasBetheHeitlerApprox AtlasBetheHeitlerApprox::loadFromFiles(
   const auto [lowData, lowTransform] = read_file(low_parameters_path);
   const auto [highData, highTransform] = read_file(high_parameters_path);
 
-  return {lowData,   highData,     lowTransform,  highTransform,      lowLimit,
-          highLimit, clampToRange, noChangeLimit, singleGaussianLimit};
+  std::vector<PolynomialBetheHeitlerApprox::RangeData> ranges = {
+      {Range1D<double>{lowLimit, highLimit}, lowData, lowTransform},
+      {Range1D<double>{highLimit, highLimit * 2}, highData, highTransform}};
+
+  return PolynomialBetheHeitlerApprox(std::move(ranges), clampToRange,
+                                      noChangeLimit, singleGaussianLimit);
 }
 
-std::span<AtlasBetheHeitlerApprox::Component> AtlasBetheHeitlerApprox::mixture(
+std::span<PolynomialBetheHeitlerApprox::Component>
+PolynomialBetheHeitlerApprox::mixture(
     double xOverX0, const std::span<Component> mixture) const {
   if (m_clampToRange) {
-    xOverX0 = std::clamp(xOverX0, 0.0, m_highLimit);
+    xOverX0 = std::clamp(xOverX0, 0.0, m_ranges.back().range.max());
   }
 
   // Evaluate polynomial at x
@@ -71,10 +110,10 @@ std::span<AtlasBetheHeitlerApprox::Component> AtlasBetheHeitlerApprox::mixture(
     return sum;
   };
 
-  // Lambda which builds the components
+  // Lambda which builds the components for a single range
   const auto make_mixture = [&](const Data &data, double xx,
                                 bool transform) -> std::span<Component> {
-    // Value initialization should garanuee that all is initialized to zero
+    // Value initialization should guarantee that all is initialized to zero
     double weight_sum = 0;
     for (std::size_t i = 0; i < data.size(); ++i) {
       // These transformations must be applied to the data according to ATHENA
@@ -114,24 +153,23 @@ std::span<AtlasBetheHeitlerApprox::Component> AtlasBetheHeitlerApprox::mixture(
     return {mixture.data(), 1};
   }
 
-  // Return a component representation for lower x0
-  if (xOverX0 < m_lowLimit) {
-    return make_mixture(m_lowData, xOverX0, m_lowTransform);
+  // Find the appropriate range and return mixture for that range
+  for (const auto &[range, data, transform] : m_ranges) {
+    if (range.contains(xOverX0)) {
+      return make_mixture(data, xOverX0, transform);
+    }
   }
 
-  // Return a component representation for higher x0
-  // Cap the x because beyond the parameterization goes wild
-  const auto high_x = std::min(m_highLimit, xOverX0);
-  return make_mixture(m_highData, high_x, m_highTransform);
+  // Should not reach here if validXOverX0 is called first
+  // But return the last range's mixture as fallback
+  const auto &[lastRange, lastData, lastTransform] = m_ranges.back();
+  return make_mixture(lastData, lastRange.max(), lastTransform);
 }
 
-}  // namespace Acts
-
-Acts::AtlasBetheHeitlerApprox Acts::makeDefaultBetheHeitlerApprox(
-    bool clampToRange) {
+PolynomialBetheHeitlerApprox makeDefaultBetheHeitlerApprox(bool clampToRange) {
   // Tracking/TrkFitter/TrkGaussianSumFilterUtils/Data/BetheHeitler_cdf_nC6_O5.par
   // clang-format off
-  static AtlasBetheHeitlerApprox::Data cdf_cmps6_order5_data = {{
+  static PolynomialBetheHeitlerApprox::Data cdf_cmps6_order5_data = {{
       // Component #1
       {
           {{3.74397e+004,-1.95241e+004, 3.51047e+003,-2.54377e+002, 1.81080e+001,-3.57643e+000}},
@@ -171,13 +209,11 @@ Acts::AtlasBetheHeitlerApprox Acts::makeDefaultBetheHeitlerApprox(
   }};
   // clang-format on
 
-  return {cdf_cmps6_order5_data,
-          cdf_cmps6_order5_data,
-          true,
-          true,
-          0.2,
-          0.2,
-          clampToRange,
-          0.0001,
-          0.002};
+  std::vector<PolynomialBetheHeitlerApprox::RangeData> ranges = {
+      {Range1D<double>{0.0, 0.2}, cdf_cmps6_order5_data, true}};
+
+  return PolynomialBetheHeitlerApprox(std::move(ranges), clampToRange, 0.0001,
+                                      0.002);
 }
+
+}  // namespace Acts