refactor: use mathematical constants from std::numbers (#3781)

blocked:
- #3783

Asked co-pilot to summarise the PR:

Summary:

The pull request aims to refactor the code to use mathematical constants from the C++20 <numbers> header, replacing existing instances of M_PI with std::numbers::pi. The changes span multiple files, primarily focusing on improving the consistency and modernity of mathematical constant usage across the codebase.

Core/include/Acts/AmbiguityResolution/ScoreBasedAmbiguityResolution.hpp

@@ -18,6 +18,7 @@
 #include <cstddef>
 #include <map>
 #include <memory>
+#include <numbers>
 #include <string>
 #include <tuple>
 #include <vector>
@@ -103,8 +104,8 @@ class ScoreBasedAmbiguityResolution {
     double pTMin = 0 * UnitConstants::GeV;
     double pTMax = 1e5 * UnitConstants::GeV;
 
-    double phiMin = -M_PI * UnitConstants::rad;
-    double phiMax = M_PI * UnitConstants::rad;
+    double phiMin = -std::numbers::pi * UnitConstants::rad;
+    double phiMax = std::numbers::pi * UnitConstants::rad;
 
     double etaMin = -5;
     double etaMax = 5;

Core/include/Acts/Definitions/Units.hpp

@@ -8,6 +8,8 @@
 
 #pragma once
 
+#include <numbers>
+
 namespace Acts {
 
 /// @verbatim embed:rst:leading-slashes
@@ -170,7 +172,7 @@ constexpr double h = 3600.0 * s;
 // Angles, native unit radian
 constexpr double mrad = 1e-3;
 constexpr double rad = 1.0;
-constexpr double degree = 0.017453292519943295;  // = M_PI / 180.0 * rad;
+constexpr double degree = std::numbers::pi / 180. / rad;
 // Energy/mass/momentum, native unit GeV
 constexpr double GeV = 1.0;
 constexpr double eV = 1e-9 * GeV;

Core/include/Acts/EventData/TransformationHelpers.hpp

@@ -15,6 +15,8 @@
 #include "Acts/Utilities/Result.hpp"
 #include "Acts/Utilities/detail/periodic.hpp"
 
+#include <numbers>
+
 namespace Acts {
 
 class Surface;
@@ -25,8 +27,9 @@ class Surface;
 /// @return Reflected bound track parameters vector
 inline BoundVector reflectBoundParameters(const BoundVector& boundParams) {
   BoundVector reflected = boundParams;
-  auto [phi, theta] = detail::normalizePhiTheta(
-      boundParams[eBoundPhi] - M_PI, M_PI - boundParams[eBoundTheta]);
+  auto [phi, theta] =
+      detail::normalizePhiTheta(boundParams[eBoundPhi] - std::numbers::pi,
+                                std::numbers::pi - boundParams[eBoundTheta]);
   reflected[eBoundPhi] = phi;
   reflected[eBoundTheta] = theta;
   reflected[eBoundQOverP] = -boundParams[eBoundQOverP];

Core/include/Acts/EventData/detail/ParameterTraits.hpp

@@ -13,6 +13,7 @@
 #include <algorithm>
 #include <cmath>
 #include <cstddef>
+#include <numbers>
 
 namespace Acts::detail {
 
@@ -106,12 +107,12 @@ struct CyclicParameterTraits {
 //
 // The functions names are chosen to be consistent w/ std::numeric_limits
 struct PhiBoundParameterLimits {
-  static constexpr double lowest() { return -M_PI; }
-  static constexpr double max() { return M_PI; }
+  static constexpr double lowest() { return -std::numbers::pi; }
+  static constexpr double max() { return std::numbers::pi; }
 };
 struct ThetaBoundParameterLimits {
   static constexpr double lowest() { return 0; }
-  static constexpr double max() { return M_PI; }
+  static constexpr double max() { return std::numbers::pi; }
 };
 
 // Traits implementation structs for single parameters.

Core/include/Acts/Geometry/CylinderVolumeBounds.hpp

@@ -18,6 +18,7 @@
 #include <initializer_list>
 #include <iosfwd>
 #include <memory>
+#include <numbers>
 #include <ostream>
 #include <vector>
 
@@ -105,8 +106,9 @@ class CylinderVolumeBounds : public VolumeBounds {
   /// @param bevelMinZ The bevel angle, in radians, for the negative side
   /// @param bevelMaxZ The bevel angle, in radians, for the positive side
   CylinderVolumeBounds(ActsScalar rmin, ActsScalar rmax, ActsScalar halfz,
-                       ActsScalar halfphi = M_PI, ActsScalar avgphi = 0.,
-                       ActsScalar bevelMinZ = 0., ActsScalar bevelMaxZ = 0.);
+                       ActsScalar halfphi = std::numbers::pi_v<ActsScalar>,
+                       ActsScalar avgphi = 0., ActsScalar bevelMinZ = 0.,
+                       ActsScalar bevelMaxZ = 0.);
 
   /// Constructor - from a fixed size array
   ///

Core/include/Acts/Geometry/SurfaceArrayCreator.hpp

@@ -24,6 +24,7 @@
 #include <functional>
 #include <iterator>
 #include <memory>
+#include <numbers>
 #include <optional>
 #include <ostream>
 #include <tuple>
@@ -257,7 +258,7 @@ class SurfaceArrayCreator {
 
       // ...so by injecting them into atan2, we get the angle between them
       auto dPhi = std::atan2(sin_dPhi_n2, cos_dPhi_n2);
-      return std::abs(dPhi) < M_PI / 180.;
+      return std::abs(dPhi) < std::numbers::pi / 180.;
     }
 
     if (bValue == Acts::BinningValue::binZ) {

Core/include/Acts/Geometry/VolumeBounds.hpp

@@ -17,6 +17,7 @@
 #include <cmath>
 #include <iostream>
 #include <memory>
+#include <numbers>
 #include <utility>
 #include <vector>
 
@@ -33,12 +34,14 @@ struct OrientedSurface {
 
 // Planar definitions to help construct the boundary surfaces
 static const Transform3 s_planeXY = Transform3::Identity();
-static const Transform3 s_planeYZ = AngleAxis3(0.5 * M_PI, Vector3::UnitY()) *
-                                    AngleAxis3(0.5 * M_PI, Vector3::UnitZ()) *
-                                    Transform3::Identity();
-static const Transform3 s_planeZX = AngleAxis3(-0.5 * M_PI, Vector3::UnitX()) *
-                                    AngleAxis3(-0.5 * M_PI, Vector3::UnitZ()) *
-                                    Transform3::Identity();
+static const Transform3 s_planeYZ =
+    AngleAxis3(std::numbers::pi / 2., Vector3::UnitY()) *
+    AngleAxis3(std::numbers::pi / 2., Vector3::UnitZ()) *
+    Transform3::Identity();
+static const Transform3 s_planeZX =
+    AngleAxis3(-std::numbers::pi / 2., Vector3::UnitX()) *
+    AngleAxis3(-std::numbers::pi / 2., Vector3::UnitZ()) *
+    Transform3::Identity();
 
 /// @class VolumeBounds
 ///

Core/include/Acts/Propagator/RiddersPropagator.ipp

@@ -8,6 +8,8 @@
 
 #include "Acts/Definitions/TrackParametrization.hpp"
 
+#include <numbers>
+
 template <typename propagator_t>
 template <typename parameters_t, typename propagator_options_t>
 auto Acts::RiddersPropagator<propagator_t>::propagate(
@@ -150,8 +152,8 @@ bool Acts::RiddersPropagator<propagator_t>::inconsistentDerivativesOnDisc(
     for (unsigned int j = 0; j < derivatives.size(); j++) {
       // If there is at least one with a similar angle then it seems to work
       // properly
-      if (i != j &&
-          std::abs(derivatives[i](1) - derivatives[j](1)) < 0.5 * M_PI) {
+      if (i != j && std::abs(derivatives[i](1) - derivatives[j](1)) <
+                        std::numbers::pi / 2.) {
         jumpedAngle = false;
         break;
       }
@@ -179,8 +181,8 @@ Acts::RiddersPropagator<propagator_t>::wiggleParameter(
     // Treatment for theta
     if (param == eBoundTheta) {
       const double current_theta = start.template get<eBoundTheta>();
-      if (current_theta + h > M_PI) {
-        h = M_PI - current_theta;
+      if (current_theta + h > std::numbers::pi) {
+        h = std::numbers::pi - current_theta;
       }
       if (current_theta + h < 0) {
         h = -current_theta;
@@ -202,10 +204,14 @@ Acts::RiddersPropagator<propagator_t>::wiggleParameter(
     if (param == eBoundPhi) {
       double phi0 = nominal(Acts::eBoundPhi);
       double phi1 = r.endParameters->parameters()(Acts::eBoundPhi);
-      if (std::abs(phi1 + 2. * M_PI - phi0) < std::abs(phi1 - phi0)) {
-        derivatives.back()[Acts::eBoundPhi] = (phi1 + 2. * M_PI - phi0) / h;
-      } else if (std::abs(phi1 - 2. * M_PI - phi0) < std::abs(phi1 - phi0)) {
-        derivatives.back()[Acts::eBoundPhi] = (phi1 - 2. * M_PI - phi0) / h;
+      if (std::abs(phi1 + 2. * std::numbers::pi - phi0) <
+          std::abs(phi1 - phi0)) {
+        derivatives.back()[Acts::eBoundPhi] =
+            (phi1 + 2. * std::numbers::pi - phi0) / h;
+      } else if (std::abs(phi1 - 2. * std::numbers::pi - phi0) <
+                 std::abs(phi1 - phi0)) {
+        derivatives.back()[Acts::eBoundPhi] =
+            (phi1 - 2. * std::numbers::pi - phi0) / h;
       }
     }
   }

Core/include/Acts/Propagator/detail/LoopProtection.hpp

@@ -11,6 +11,8 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Utilities/Logger.hpp"
 
+#include <numbers>
+
 namespace Acts::detail {
 
 /// Estimate the loop protection limit
@@ -46,7 +48,8 @@ void setupLoopProtection(propagator_state_t& state, const stepper_t& stepper,
   // Transverse component at start is taken for the loop protection
   const double p = stepper.absoluteMomentum(state.stepping);
   // Calculate the full helix path
-  const double helixPath = state.options.direction * 2 * M_PI * p / B;
+  const double helixPath =
+      state.options.direction * 2 * std::numbers::pi * p / B;
   // And set it as the loop limit if it overwrites the internal limit
   const double loopLimit = state.options.loopFraction * helixPath;
   const double previousLimit = pathAborter.internalLimit;

Core/include/Acts/Seeding/GbtsDataStorage.hpp

@@ -14,6 +14,7 @@
 
 #include <algorithm>
 #include <map>
+#include <numbers>
 #include <vector>
 
 namespace Acts {
@@ -117,12 +118,12 @@ class GbtsEtaBin {
       // float phi = pN->m_sp.phi();
       // float phi = (std::atan(pN->m_sp.x() / pN->m_sp.y()));
       float phi = pN->m_spGbts.phi();
-      if (phi <= M_PI - dphi) {
+      if (phi <= std::numbers::pi_v<float> - dphi) {
         continue;
       }
 
-      m_vPhiNodes.push_back(
-          std::pair<float, unsigned int>(phi - 2 * M_PI, nIdx));
+      m_vPhiNodes.push_back(std::pair<float, unsigned int>(
+          phi - static_cast<float>(2. * std::numbers::pi), nIdx));
     }
 
     for (unsigned int nIdx = 0; nIdx < m_vn.size(); nIdx++) {
@@ -134,11 +135,11 @@ class GbtsEtaBin {
     for (unsigned int nIdx = 0; nIdx < m_vn.size(); nIdx++) {
       GbtsNode<space_point_t> *pN = m_vn.at(nIdx);
       float phi = pN->m_spGbts.phi();
-      if (phi >= -M_PI + dphi) {
+      if (phi >= -std::numbers::pi_v<float> + dphi) {
         break;
       }
-      m_vPhiNodes.push_back(
-          std::pair<float, unsigned int>(phi + 2 * M_PI, nIdx));
+      m_vPhiNodes.push_back(std::pair<float, unsigned int>(
+          phi + static_cast<float>(2. * std::numbers::pi), nIdx));
     }
   }
 

Core/include/Acts/Seeding/SeedFinderConfig.hpp

@@ -15,6 +15,7 @@
 
 #include <limits>
 #include <memory>
+#include <numbers>
 #include <vector>
 
 namespace Acts {
@@ -33,8 +34,8 @@ struct SeedFinderConfig {
   /// Geometry Settings + Detector ROI
   /// (r, z, phi) range for limiting location of all measurements and grid
   /// creation
-  float phiMin = -M_PI;
-  float phiMax = M_PI;
+  float phiMin = -std::numbers::pi_v<float>;
+  float phiMax = std::numbers::pi_v<float>;
   float zMin = -2800 * Acts::UnitConstants::mm;
   float zMax = 2800 * Acts::UnitConstants::mm;
   float rMax = 600 * Acts::UnitConstants::mm;

Core/include/Acts/Seeding/SeedFinderGbts.ipp

@@ -9,7 +9,6 @@
 // SeedFinderGbts.ipp
 // TODO: update to C++17 style
 
-#include "Acts/Definitions/Algebra.hpp"  //for M_PI
 #include "Acts/Geometry/Extent.hpp"
 #include "Acts/Seeding/SeedFilter.hpp"
 #include "Acts/Seeding/SeedFinder.hpp"
@@ -23,6 +22,7 @@
 #include <functional>
 #include <iostream>
 #include <list>
+#include <numbers>
 #include <numeric>
 #include <type_traits>
 #include <vector>
@@ -54,7 +54,7 @@ void SeedFinderGbts<external_spacepoint_t>::loadSpacePoints(
     m_storage->addSpacePoint(gbtssp, (m_config.m_useClusterWidth > 0));
   }
 
-  m_config.m_phiSliceWidth = 2 * M_PI / m_config.m_nMaxPhiSlice;
+  m_config.m_phiSliceWidth = 2 * std::numbers::pi / m_config.m_nMaxPhiSlice;
 
   m_storage->sortByPhi();
 
@@ -390,10 +390,10 @@ void SeedFinderGbts<external_spacepoint_t>::runGbts_TrackFinder(
 
         float dPhi = pNS->m_p[3] - Phi1;
 
-        if (dPhi < -M_PI) {
-          dPhi += 2 * M_PI;
-        } else if (dPhi > M_PI) {
-          dPhi -= 2 * M_PI;
+        if (dPhi < -std::numbers::pi_v<float>) {
+          dPhi += static_cast<float>(2 * std::numbers::pi);
+        } else if (dPhi > std::numbers::pi_v<float>) {
+          dPhi -= static_cast<float>(2 * std::numbers::pi);
         }
 
         if (dPhi < -m_config.cut_dphi_max || dPhi > m_config.cut_dphi_max) {

Core/include/Acts/Seeding/SeedFinderOrthogonalConfig.hpp

@@ -14,6 +14,7 @@
 #include "Acts/Utilities/Delegate.hpp"
 
 #include <memory>
+#include <numbers>
 
 namespace Acts {
 // forward declaration to avoid cyclic dependence
@@ -28,8 +29,8 @@ struct SeedFinderOrthogonalConfig {
   /// Seeding parameters for geometry settings and detector ROI
 
   // Limiting location of all measurements
-  float phiMin = -M_PI;
-  float phiMax = M_PI;
+  float phiMin = -std::numbers::pi_v<float>;
+  float phiMax = std::numbers::pi_v<float>;
   /// limiting location of measurements
   float zMin = -2800 * Acts::UnitConstants::mm;
   float zMax = 2800 * Acts::UnitConstants::mm;

Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp

@@ -71,7 +71,7 @@ struct CylindricalSpacePointGridConfig {
   // maximum phi value for phiAxis construction
   float phiMax = std::numbers::pi_v<float>;
   // Multiplicator for the number of phi-bins. The minimum number of phi-bins
-  // depends on min_pt, magnetic field: 2*M_PI/(minPT particle phi-deflection).
+  // depends on min_pt, magnetic field: 2*pi/(minPT particle phi-deflection).
   // phiBinDeflectionCoverage is a multiplier for this number. If
   // numPhiNeighbors (in the configuration of the BinFinders) is configured to
   // return 1 neighbor on either side of the current phi-bin (and you want to

Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp

@@ -7,6 +7,7 @@
 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
 #include <concepts>
+#include <numbers>
 
 template <typename external_spacepoint_t>
 Acts::CylindricalSpacePointGrid<external_spacepoint_t>
@@ -91,16 +92,14 @@ Acts::CylindricalSpacePointGridCreator::createGrid(
 
     // divide 2pi by angle delta to get number of phi-bins
     // size is always 2pi even for regions of interest
-    phiBins = static_cast<int>(std::ceil(2 * M_PI / deltaPhi));
+    phiBins = static_cast<int>(std::ceil(2 * std::numbers::pi / deltaPhi));
     // need to scale the number of phi bins accordingly to the number of
     // consecutive phi bins in the seed making step.
     // Each individual bin should be approximately a fraction (depending on this
     // number) of the maximum expected azimutal deflection.
 
     // set protection for large number of bins, by default it is large
-    if (phiBins > config.maxPhiBins) {
-      phiBins = config.maxPhiBins;
-    }
+    phiBins = std::min(phiBins, config.maxPhiBins);
   }
 
   Acts::Axis<AxisType::Equidistant, AxisBoundaryType::Closed> phiAxis(

Core/include/Acts/Surfaces/AnnulusBounds.hpp

@@ -20,6 +20,7 @@
 #include <cmath>
 #include <exception>
 #include <iosfwd>
+#include <numbers>
 #include <stdexcept>
 #include <vector>
 
@@ -231,7 +232,7 @@ inline double AnnulusBounds::phiMax() const {
 }
 
 inline bool AnnulusBounds::coversFullAzimuth() const {
-  return (std::abs((get(eMinPhiRel) - get(eMaxPhiRel)) - M_PI) <
+  return (std::abs((get(eMinPhiRel) - get(eMaxPhiRel)) - std::numbers::pi) <
           s_onSurfaceTolerance);
 }