fix: FPE in AnnealingUtility by using safeExp (#2406)

CI/physmon/fpe_masks.yml

@@ -4,11 +4,5 @@
   FLTUND: 1
 "Core/include/Acts/TrackFitting/detail/GsfUtils.hpp:197":
   FLTUND: 1  
-"Core/src/Utilities/AnnealingUtility.cpp:43":
-  FLTUND: 1
-"Core/src/Utilities/AnnealingUtility.cpp:40":
-  FLTUND: 1
-"Core/src/Utilities/AnnealingUtility.cpp:38":
-  FLTUND: 1
 "Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.ipp:120":
   FLTUND: 1

Core/include/Acts/Utilities/AlgebraHelpers.hpp

@@ -198,21 +198,35 @@ std::optional<ResultType> safeInverse(const MatrixType& m) noexcept {
   return std::nullopt;
 }
 
+/// Specialization of the exponent limit to be used for safe exponential,
+/// depending on the floating point type.
+/// See https://godbolt.org/z/z53Er6Mzf for reasoning for the concrete numbers.
+template <typename T>
+struct ExpSafeLimit {};
+template <>
+struct ExpSafeLimit<double> {
+  constexpr static double value = 500.0;
+};
+template <>
+struct ExpSafeLimit<float> {
+  constexpr static float value = 50.0;
+};
+
 /// Calculate the exponential function while avoiding FPEs.
-/// @note The boundary values of -50.0 and 50.0 might need to be adapted when
-/// using this function with doubles
 ///
 /// @param val argument for which the exponential function should be evaluated.
 ///
 /// @return 0 in the case of underflow, std::numeric_limits<T>::infinity in the
 /// case of overflow, std::exp(val) else
 template <typename T>
-T safeExp(T val) noexcept {
-  if (val < -50.0) {
+constexpr T safeExp(T val) noexcept {
+  constexpr T maxExponent = ExpSafeLimit<T>::value;
+  constexpr T minExponent = -maxExponent;
+  if (val < minExponent) {
     return 0.0;
   }
 
-  if (val > 50.0) {
+  if (val > maxExponent) {
     return std::numeric_limits<T>::infinity();
   }
 

Core/src/Utilities/AnnealingUtility.cpp

@@ -8,16 +8,20 @@
 
 #include "Acts/Utilities/AnnealingUtility.hpp"
 
+#include "Acts/Utilities/AlgebraHelpers.hpp"
+
+namespace {
 /// @brief Gaussian-like function for weight calculation
 /// Note: Factor 2 in denominator is included in inverse temperature
 ///
 /// @param chi2 Chi^2 value
 /// @param invTemp Denominator 1/(2 * temperature)
 ///
 /// @return exp(-chi2 * invTemp)
-static double computeAnnealingWeight(double chi2, double invTemp) {
-  return std::exp(-chi2 * invTemp);
+double computeAnnealingWeight(double chi2, double invTemp) {
+  return Acts::safeExp(-chi2 * invTemp);
 }
+}  // namespace
 
 Acts::AnnealingUtility::Config::Config() = default;