refactor: Improve Bethe Heitler approximation constructors and logging

Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp

@@ -98,6 +98,9 @@ struct BetheHeitlerApproxSingleCmp {
 /// 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.
 template <int NComponents, int PolyDegree>
 class AtlasBetheHeitlerApprox {
   static_assert(NComponents > 0);
@@ -112,40 +115,47 @@ class AtlasBetheHeitlerApprox {
 
   using Data = std::array<PolyData, NComponents>;
 
-  constexpr static double noChangeLimit = 0.0001;
-  constexpr static double singleGaussianLimit = 0.002;
-  constexpr static double lowerLimit = 0.10;
-  constexpr static double higherLimit = 0.20;
-
  private:
-  Data m_low_data;
-  Data m_high_data;
-  bool m_low_transform;
-  bool m_high_transform;
+  Data m_lowData;
+  Data m_highData;
+  bool m_lowTransform;
+  bool m_highTransform;
+
+  constexpr static double m_noChangeLimit = 0.0001;
+  constexpr static double m_singleGaussianLimit = 0.002;
+  double m_lowLimit = 0.10;
+  double m_highLimit = 0.20;
 
  public:
-  /// Construct the Bethe-Heitler approximation description. Additional to the
-  /// coefficients of the polynomials, the information whether these values need
-  /// to be transformed beforehand must be given (see ATLAS code).
+  /// 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
+  /// to the upper limit of the low-x/x0 data.
   ///
-  /// @param low_data data for the lower x/x0 range
-  /// @param high_data data for the higher x/x0 range
-  /// @param low_transform whether the low data need to be transformed
-  /// @param high_transform whether the high data need to be transformed
-  constexpr AtlasBetheHeitlerApprox(const Data &low_data, const Data &high_data,
-                                    bool low_transform, bool high_transform)
-      : m_low_data(low_data),
-        m_high_data(high_data),
-        m_low_transform(low_transform),
-        m_high_transform(high_transform) {}
+  /// @param lowData data for the lower x/x0 range
+  /// @param highData data for the higher x/x0 range
+  /// @param lowTransform whether the low data need to be transformed
+  /// @param highTransform whether the high data need to be transformed
+  /// @param lowLimit the upper limit for the low data
+  /// @param highLimit the upper limit for the high data
+  constexpr AtlasBetheHeitlerApprox(const Data &lowData, const Data &highData,
+                                    bool lowTransform, bool highTransform,
+                                    double lowLimit = 0.1,
+                                    double highLimit = 0.2)
+      : m_lowData(lowData),
+        m_highData(highData),
+        m_lowTransform(lowTransform),
+        m_highTransform(highTransform),
+        m_lowLimit(lowLimit),
+        m_highLimit(highLimit) {}
 
   /// Returns the number of components the returned mixture will have
   constexpr auto numComponents() const { return NComponents; }
 
   /// Checks if an input is valid for the parameterization
   ///
   /// @param x pathlength in terms of the radiation length
-  constexpr bool validXOverX0(ActsScalar x) const { return x < higherLimit; }
+  constexpr bool validXOverX0(ActsScalar x) const { return x < m_highLimit; }
 
   /// Generates the mixture from the polynomials and reweights them, so
   /// that the sum of all weights is 1
@@ -194,7 +204,7 @@ class AtlasBetheHeitlerApprox {
     };
 
     // Return no change
-    if (x < noChangeLimit) {
+    if (x < m_noChangeLimit) {
       Array ret(1);
 
       ret[0].weight = 1.0;
@@ -204,19 +214,19 @@ class AtlasBetheHeitlerApprox {
       return ret;
     }
     // Return single gaussian approximation
-    if (x < singleGaussianLimit) {
+    if (x < m_singleGaussianLimit) {
       Array ret(1);
       ret[0] = BetheHeitlerApproxSingleCmp::mixture(x)[0];
       return ret;
     }
     // Return a component representation for lower x0
-    if (x < lowerLimit) {
-      return make_mixture(m_low_data, x, m_low_transform);
+    if (x < m_lowLimit) {
+      return make_mixture(m_lowData, x, m_lowTransform);
     }
     // Return a component representation for higher x0
     // Cap the x because beyond the parameterization goes wild
-    const auto high_x = std::min(higherLimit, x);
-    return make_mixture(m_high_data, high_x, m_high_transform);
+    const auto high_x = std::min(m_highLimit, x);
+    return make_mixture(m_highData, high_x, m_highTransform);
   }
 
   /// Loads a parameterization from a file according to the Atlas file
@@ -226,8 +236,11 @@ class AtlasBetheHeitlerApprox {
   /// the parameterization for low x/x0
   /// @param high_parameters_path Path to the foo.par file that stores
   /// the parameterization for high x/x0
+  /// @param lowLimit the upper limit for the low x/x0-data
+  /// @param highLimit the upper limit for the high x/x0-data
   static auto loadFromFiles(const std::string &low_parameters_path,
-                            const std::string &high_parameters_path) {
+                            const std::string &high_parameters_path,
+                            double lowLimit = 0.1, double highLimit = 0.2) {
     auto read_file = [](const std::string &filepath) {
       std::ifstream file(filepath);
 
@@ -267,11 +280,11 @@ class AtlasBetheHeitlerApprox {
       return std::make_tuple(data, transform_code);
     };
 
-    const auto [low_data, low_transform] = read_file(low_parameters_path);
-    const auto [high_data, high_transform] = read_file(high_parameters_path);
+    const auto [lowData, lowTransform] = read_file(low_parameters_path);
+    const auto [highData, highTransform] = read_file(high_parameters_path);
 
-    return AtlasBetheHeitlerApprox(low_data, high_data, low_transform,
-                                   high_transform);
+    return AtlasBetheHeitlerApprox(lowData, highData, lowTransform,
+                                   highTransform, lowLimit, highLimit);
   }
 };
 

Core/include/Acts/TrackFitting/GaussianSumFitter.hpp

@@ -322,6 +322,7 @@ struct GaussianSumFitter {
     ACTS_VERBOSE("- measurement states: " << fwdGsfResult.measurementStates);
 
     std::size_t nInvalidBetheHeitler = fwdGsfResult.nInvalidBetheHeitler;
+    double maxPathXOverX0 = fwdGsfResult.maxPathXOverX0;
 
     //////////////////
     // Backward pass
@@ -408,13 +409,16 @@ struct GaussianSumFitter {
     }
 
     nInvalidBetheHeitler += bwdGsfResult.nInvalidBetheHeitler;
+    maxPathXOverX0 = std::max(maxPathXOverX0, bwdGsfResult.maxPathXOverX0);
 
     if (nInvalidBetheHeitler > 0) {
-      ACTS_WARNING("Encountered "
-                   << nInvalidBetheHeitler
-                   << " cases where the material thickness exceeds the range "
-                      "of the Bethe-Heitler-Approximation. Enable DEBUG output "
-                      "for more information.");
+      ACTS_WARNING("Encountered " << nInvalidBetheHeitler
+                                  << " cases where x/X0 exceeds the range "
+                                     "of the Bethe-Heitler-Approximation. The "
+                                     "maximum x/X0 encountered was "
+                                  << maxPathXOverX0
+                                  << ". Enable DEBUG output "
+                                     "for more information.");
     }
 
     ////////////////////////////////////

Core/include/Acts/TrackFitting/detail/GsfActor.hpp

@@ -56,7 +56,9 @@ struct GsfResult {
   std::vector<const Acts::Surface*> visitedSurfaces;
   std::vector<const Acts::Surface*> surfacesVisitedBwdAgain;
 
+  /// Statistics about encounterings of invalid bethe heitler approximations
   std::size_t nInvalidBetheHeitler = 0;
+  double maxPathXOverX0 = 0.0;
 
   // Propagate potential errors to the outside
   Result<void> result{Result<void>::success()};
@@ -287,7 +289,7 @@ struct GsfActor {
       componentCache.clear();
 
       convoluteComponents(state, stepper, navigator, tmpStates, componentCache,
-                          result.nInvalidBetheHeitler);
+                          result);
 
       if (componentCache.empty()) {
         ACTS_WARNING(
@@ -328,7 +330,7 @@ struct GsfActor {
                            const navigator_t& navigator,
                            const TemporaryStates& tmpStates,
                            std::vector<ComponentCache>& componentCache,
-                           std::size_t& nInvalidBetheHeitler) const {
+                           result_type& result) const {
     auto cmps = stepper.componentIterable(state.stepping);
     for (auto [idx, cmp] : zip(tmpStates.tips, cmps)) {
       auto proxy = tmpStates.traj.getTrackState(idx);
@@ -338,7 +340,7 @@ struct GsfActor {
                                  stepper.particleHypothesis(state.stepping));
 
       applyBetheHeitler(state, navigator, bound, tmpStates.weights.at(idx),
-                        componentCache, nInvalidBetheHeitler);
+                        componentCache, result);
     }
   }
 
@@ -348,7 +350,7 @@ struct GsfActor {
                          const BoundTrackParameters& old_bound,
                          const double old_weight,
                          std::vector<ComponentCache>& componentCaches,
-                         std::size_t& nInvalidBetheHeitler) const {
+                         result_type& result) const {
     const auto& surface = *navigator.currentSurface(state.navigation);
     const auto p_prev = old_bound.absoluteMomentum();
 
@@ -357,22 +359,24 @@ struct GsfActor {
         old_bound.position(state.stepping.geoContext), state.options.direction,
         MaterialUpdateStage::FullUpdate);
 
-    auto pathCorrection = surface.pathCorrection(
+    const auto pathCorrection = surface.pathCorrection(
         state.stepping.geoContext,
         old_bound.position(state.stepping.geoContext), old_bound.direction());
     slab.scaleThickness(pathCorrection);
 
+    const double pathXOverX0 = slab.thicknessInX0();
+
     // Emit a warning if the approximation is not valid for this x/x0
-    if (!m_cfg.bethe_heitler_approx->validXOverX0(slab.thicknessInX0())) {
-      ++nInvalidBetheHeitler;
+    if (!m_cfg.bethe_heitler_approx->validXOverX0(pathXOverX0)) {
+      ++result.nInvalidBetheHeitler;
+      result.maxPathXOverX0 = std::max(result.maxPathXOverX0, pathXOverX0);
       ACTS_DEBUG(
           "Bethe-Heitler approximation encountered invalid value for x/x0="
-          << slab.thicknessInX0() << " at surface " << surface.geometryId());
+          << pathXOverX0 << " at surface " << surface.geometryId());
     }
 
     // Get the mixture
-    const auto mixture =
-        m_cfg.bethe_heitler_approx->mixture(slab.thicknessInX0());
+    const auto mixture = m_cfg.bethe_heitler_approx->mixture(pathXOverX0);
 
     // Create all possible new components
     for (const auto& gaussian : mixture) {

Core/src/TrackFitting/BetheHeitlerApprox.cpp

@@ -51,6 +51,6 @@ Acts::AtlasBetheHeitlerApprox<6, 5> Acts::makeDefaultBetheHeitlerApprox() {
   }};
   // clang-format on
 
-  return AtlasBetheHeitlerApprox<6, 5>(cdf_cmps6_order5_data,
-                                       cdf_cmps6_order5_data, true, true);
+  return AtlasBetheHeitlerApprox<6, 5>(
+      cdf_cmps6_order5_data, cdf_cmps6_order5_data, true, true, 0.2, 0.2);
 }

Examples/Python/src/TrackFitting.cpp

@@ -109,7 +109,8 @@ void addTrackFitting(Context& ctx) {
     py::class_<ActsExamples::BetheHeitlerApprox>(mex, "AtlasBetheHeitlerApprox")
         .def_static("loadFromFiles",
                     &ActsExamples::BetheHeitlerApprox::loadFromFiles,
-                    py::arg("lowParametersPath"), py::arg("lowParametersPath"))
+                    py::arg("lowParametersPath"), py::arg("highParametersPath"),
+                    py::arg("lowLimit") = 0.1, py::arg("highLimit") = 0.2)
         .def_static("makeDefault",
                     []() { return Acts::makeDefaultBetheHeitlerApprox(); });
     mex.def(