diff --git a/.github/workflows/builds.yml b/.github/workflows/builds.yml
index 06312a6b339..b50f3369f4d 100644
--- a/.github/workflows/builds.yml
+++ b/.github/workflows/builds.yml
@@ -303,6 +303,62 @@ jobs:
       - name: Downstream run
         run: ./build-downstream/bin/ShowActsVersion
 
+  linux_float:
+    # This job tests the correct usage of the type ActsScalar by testing it with float instead of double
+    runs-on: ubuntu-latest
+    container: ghcr.io/acts-project/ubuntu2404:58
+    env:
+      INSTALL_DIR: ${{ github.workspace }}/install
+      ACTS_LOG_FAILURE_THRESHOLD: WARNING
+    steps:
+      - name: Install git lfs
+        run: apt-get update && apt-get install -y git-lfs
+
+      - uses: actions/checkout@v4
+        with:
+          submodules: true
+          lfs: true
+
+      - name: Restore ccache
+        uses: actions/cache/restore@v4
+        id: ccache-restore
+        with:
+          path: ${{ github.workspace }}/ccache
+          key: ${{ runner.os  }}-ccache-linux_float_${{ env.CCACHE_KEY_SUFFIX }}_${{ github.sha }}
+          restore-keys: |
+            ${{ runner.os  }}-ccache-linux_float_${{ env.CCACHE_KEY_SUFFIX }}_
+
+      - name: Configure
+        # setting CMAKE_CXX_STANDARD=20 is a workaround for a bug in the
+        # dd4hep CMake configuration that gets triggered on recent CMake
+        # versions
+        # Need to set git user & email for patching to work (GeoModel plugin)
+        run: >
+          git config --global user.name 'CI' &&
+          git config --global user.email '<>' &&
+          ccache -z &&
+          cmake -B build -S .
+          -DCMAKE_CXX_FLAGS="-Wno-narrowing"
+          -DACTS_CUSTOM_SCALARTYPE=float
+
+      - name: Build
+        run: cmake --build build
+
+      - name: ccache stats
+        run: ccache -s
+
+      - name: Save ccache
+        uses: actions/cache/save@v4
+        if: always()
+        with:
+          path: ${{ github.workspace }}/ccache
+          key: ${{ steps.ccache-restore.outputs.cache-primary-key }}
+
+      - uses: actions/upload-artifact@v4
+        with:
+          name: acts-linux-ubuntu
+          path: build.tar.gz
+
   macos:
     runs-on: macos-14
     env:
diff --git a/Core/include/Acts/Definitions/Algebra.hpp b/Core/include/Acts/Definitions/Algebra.hpp
index 843958a3282..07a393d30e4 100644
--- a/Core/include/Acts/Definitions/Algebra.hpp
+++ b/Core/include/Acts/Definitions/Algebra.hpp
@@ -8,6 +8,8 @@
 
 #pragma once
 
+#include <string>
+
 // for GNU: ignore this specific warning, otherwise just include Eigen/Dense
 #if defined(__GNUC__) && !defined(__clang__) && !defined(__INTEL_COMPILER)
 #pragma GCC diagnostic push
@@ -102,4 +104,12 @@ using Transform3 = Eigen::Transform<ActsScalar, 3, Eigen::Affine>;
 
 constexpr ActsScalar s_transformEquivalentTolerance = 1e-9;
 
+// User-defined literals for creating ActsScalar. Use as 3_scalar
+constexpr ActsScalar operator""_scalar(long double n) {
+  return static_cast<ActsScalar>(n);
+}
+constexpr ActsScalar operator""_scalar(unsigned long long n) {
+  return static_cast<ActsScalar>(n);
+}
+
 }  // namespace Acts
diff --git a/Core/include/Acts/EventData/ChargeConcept.hpp b/Core/include/Acts/EventData/ChargeConcept.hpp
index 448722bd4b3..288cf573c1e 100644
--- a/Core/include/Acts/EventData/ChargeConcept.hpp
+++ b/Core/include/Acts/EventData/ChargeConcept.hpp
@@ -20,16 +20,16 @@
 namespace Acts {
 
 template <typename C>
-concept ChargeConcept = requires(C c, C c2, float f, double d) {
+concept ChargeConcept = requires(C c, C c2, float f, ActsScalar a) {
   { C{f} };
 
   { c == c2 } -> std::same_as<bool>;
   { c != c2 } -> std::same_as<bool>;
 
   { c.absQ() } -> std::same_as<float>;
-  { c.extractCharge(d) } -> std::same_as<float>;
-  { c.extractMomentum(d) } -> std::same_as<double>;
-  { c.qOverP(d, d) } -> std::same_as<double>;
+  { c.extractCharge(a) } -> std::same_as<float>;
+  { c.extractMomentum(a) } -> std::same_as<ActsScalar>;
+  { c.qOverP(a, f) } -> std::same_as<ActsScalar>;
 };
 
 }  // namespace Acts
diff --git a/Core/include/Acts/EventData/TrackParameterHelpers.hpp b/Core/include/Acts/EventData/TrackParameterHelpers.hpp
index cd68b5ae6a4..1ae4f2ce55f 100644
--- a/Core/include/Acts/EventData/TrackParameterHelpers.hpp
+++ b/Core/include/Acts/EventData/TrackParameterHelpers.hpp
@@ -12,6 +12,8 @@
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/Utilities/detail/periodic.hpp"
 
+#include <numbers>
+
 namespace Acts {
 
 /// Normalize the bound parameter angles
@@ -36,8 +38,8 @@ inline BoundVector normalizeBoundParameters(const BoundVector& boundParams) {
 inline BoundVector subtractBoundParameters(const BoundVector& lhs,
                                            const BoundVector& rhs) {
   BoundVector result = lhs - rhs;
-  result[eBoundPhi] =
-      detail::difference_periodic(lhs[eBoundPhi], rhs[eBoundPhi], 2 * M_PI);
+  result[eBoundPhi] = detail::difference_periodic(
+      lhs[eBoundPhi], rhs[eBoundPhi], ActsScalar{2 * std::numbers::pi});
   return result;
 }
 
diff --git a/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp b/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp
index fb9d6ca15a7..e150db4a5d9 100644
--- a/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp
@@ -56,8 +56,8 @@ class CutoutCylinderVolumeBounds : public VolumeBounds {
   /// @param rmax The outer radius of the overall shape
   /// @param hlZ The longer halflength of the shape
   /// @param hlZc The cutout halflength of the shape
-  CutoutCylinderVolumeBounds(double rmin, double rmed, double rmax, double hlZ,
-                             double hlZc) noexcept(false)
+  CutoutCylinderVolumeBounds(ActsScalar rmin, ActsScalar rmed, ActsScalar rmax,
+                             ActsScalar hlZ, ActsScalar hlZc) noexcept(false)
       : m_values({rmin, rmed, rmax, hlZ, hlZc}) {
     checkConsistency();
     buildSurfaceBounds();
@@ -66,8 +66,8 @@ class CutoutCylinderVolumeBounds : public VolumeBounds {
   /// Constructor - from a fixed size array
   ///
   /// @param values The bound values
-  CutoutCylinderVolumeBounds(const std::array<double, eSize>& values) noexcept(
-      false)
+  CutoutCylinderVolumeBounds(
+      const std::array<ActsScalar, eSize>& values) noexcept(false)
       : m_values(values) {
     checkConsistency();
     buildSurfaceBounds();
@@ -82,14 +82,14 @@ class CutoutCylinderVolumeBounds : public VolumeBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// Inside method to test whether a point is inside the shape
   ///
   /// @param gpos The point to test
   /// @param tol The tolerance to test with
   /// @return Whether the point is inside or not.
-  bool inside(const Vector3& gpos, double tol = 0) const override;
+  bool inside(const Vector3& gpos, ActsScalar tol = 0) const override;
 
   /// Oriented surfaces, i.e. the decomposed boundary surfaces and the
   /// according navigation direction into the volume given the normal
@@ -134,7 +134,7 @@ class CutoutCylinderVolumeBounds : public VolumeBounds {
   double get(BoundValues bValue) const { return m_values[bValue]; }
 
  private:
-  std::array<double, eSize> m_values;
+  std::array<ActsScalar, eSize> m_values;
 
   // The surface bound objects
   std::shared_ptr<const CylinderBounds> m_innerCylinderBounds{nullptr};
@@ -151,8 +151,8 @@ class CutoutCylinderVolumeBounds : public VolumeBounds {
   void checkConsistency() noexcept(false);
 };
 
-inline std::vector<double> CutoutCylinderVolumeBounds::values() const {
-  std::vector<double> valvector;
+inline std::vector<ActsScalar> CutoutCylinderVolumeBounds::values() const {
+  std::vector<ActsScalar> valvector;
   valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
   return valvector;
 }
diff --git a/Core/include/Acts/Geometry/GenericCuboidVolumeBounds.hpp b/Core/include/Acts/Geometry/GenericCuboidVolumeBounds.hpp
index 5f063a62842..14735c3b29d 100644
--- a/Core/include/Acts/Geometry/GenericCuboidVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/GenericCuboidVolumeBounds.hpp
@@ -46,7 +46,7 @@ class GenericCuboidVolumeBounds : public VolumeBounds {
   ///
   /// @param values The input values
   GenericCuboidVolumeBounds(
-      const std::array<double, BoundValues::eSize>& values) noexcept(false);
+      const std::array<ActsScalar, BoundValues::eSize>& values) noexcept(false);
 
   ~GenericCuboidVolumeBounds() override = default;
 
@@ -57,7 +57,7 @@ class GenericCuboidVolumeBounds : public VolumeBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// Checking if position given in volume frame is inside
   ///
@@ -65,7 +65,7 @@ class GenericCuboidVolumeBounds : public VolumeBounds {
   /// @param tol is the tolerance applied for the inside check
   ///
   /// @return boolean indicating if the position is inside
-  bool inside(const Vector3& gpos, double tol = 0.) const override;
+  bool inside(const Vector3& gpos, ActsScalar tol = 0.) const override;
 
   /// Oriented surfaces, i.e. the decomposed boundary surfaces and the
   /// according navigation direction into the volume given the normal
diff --git a/Core/include/Acts/Geometry/VolumeBounds.hpp b/Core/include/Acts/Geometry/VolumeBounds.hpp
index 3b5810ff4c8..ac29462dde4 100644
--- a/Core/include/Acts/Geometry/VolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/VolumeBounds.hpp
@@ -86,7 +86,7 @@ class VolumeBounds {
   /// vector containing the parameters needed to describe these bounds
   ///
   /// @return of the stored values for this SurfaceBounds object
-  virtual std::vector<double> values() const = 0;
+  virtual std::vector<ActsScalar> values() const = 0;
 
   /// Checking if position given in volume frame is inside
   ///
@@ -94,7 +94,7 @@ class VolumeBounds {
   /// @param tol is the tolerance applied for the inside check
   ///
   /// @return boolean indicating if the position is inside
-  virtual bool inside(const Vector3& gpos, double tol = 0.) const = 0;
+  virtual bool inside(const Vector3& gpos, ActsScalar tol = 0.) const = 0;
 
   /// Oriented surfaces, i.e. the decomposed boundary surfaces and the
   /// according navigation direction into the volume given the normal
@@ -143,7 +143,7 @@ class VolumeBounds {
   /// @param bValue is the binning schema used
   ///
   /// @return float offset to be used for the binning
-  virtual double binningBorder(BinningValue bValue) const;
+  virtual ActsScalar binningBorder(BinningValue bValue) const;
 
   /// Output Method for std::ostream, to be overloaded by child classes
   ///
@@ -157,8 +157,8 @@ inline Vector3 VolumeBounds::binningOffset(
   return Vector3(0., 0., 0.);
 }
 
-inline double VolumeBounds::binningBorder(BinningValue /*bValue*/) const {
-  return 0.;
+inline ActsScalar VolumeBounds::binningBorder(BinningValue /*bValue*/) const {
+  return 0_scalar;
 }
 
 /// Overload of << operator for std::ostream for debug output
diff --git a/Core/include/Acts/MagneticField/InterpolatedBFieldMap.hpp b/Core/include/Acts/MagneticField/InterpolatedBFieldMap.hpp
index d65f8113d9e..20414138e44 100644
--- a/Core/include/Acts/MagneticField/InterpolatedBFieldMap.hpp
+++ b/Core/include/Acts/MagneticField/InterpolatedBFieldMap.hpp
@@ -32,12 +32,12 @@ class InterpolatedMagneticField : public MagneticFieldProvider {
   /// @brief get the minimum value of all axes of the field map
   ///
   /// @return vector returning the minima of all field map axes
-  virtual std::vector<double> getMin() const = 0;
+  virtual std::vector<ActsScalar> getMin() const = 0;
 
   /// @brief get the maximum value of all axes of the field map
   ///
   /// @return vector returning the maxima of all field map axes
-  virtual std::vector<double> getMax() const = 0;
+  virtual std::vector<ActsScalar> getMax() const = 0;
 
   /// @brief check whether given 3D position is inside look-up domain
   ///
@@ -98,8 +98,8 @@ class InterpolatedBFieldMap : public InterpolatedMagneticField {
     ///                         each Dimension)
     /// @param [in] fieldValues field values at the hyper box corners sorted in
     ///                         the canonical order defined in Acts::interpolate
-    FieldCell(std::array<double, DIM_POS> lowerLeft,
-              std::array<double, DIM_POS> upperRight,
+    FieldCell(std::array<ActsScalar, DIM_POS> lowerLeft,
+              std::array<ActsScalar, DIM_POS> upperRight,
               std::array<Vector3, N> fieldValues)
         : m_lowerLeft(std::move(lowerLeft)),
           m_upperRight(std::move(upperRight)),
@@ -132,10 +132,10 @@ class InterpolatedBFieldMap : public InterpolatedMagneticField {
 
    private:
     /// generalized lower-left corner of the confining hyper-box
-    std::array<double, DIM_POS> m_lowerLeft;
+    std::array<ActsScalar, DIM_POS> m_lowerLeft;
 
     /// generalized upper-right corner of the confining hyper-box
-    std::array<double, DIM_POS> m_upperRight;
+    std::array<ActsScalar, DIM_POS> m_upperRight;
 
     /// @brief magnetic field vectors at the hyper-box corners
     ///
@@ -225,15 +225,15 @@ class InterpolatedBFieldMap : public InterpolatedMagneticField {
   /// @brief get the minimum value of all axes of the field map
   ///
   /// @return vector returning the minima of all field map axes
-  std::vector<double> getMin() const final {
-    return std::vector<double>(m_lowerLeft.begin(), m_lowerLeft.end());
+  std::vector<ActsScalar> getMin() const final {
+    return std::vector<ActsScalar>(m_lowerLeft.begin(), m_lowerLeft.end());
   }
 
   /// @brief get the maximum value of all axes of the field map
   ///
   /// @return vector returning the maxima of all field map axes
-  std::vector<double> getMax() const final {
-    return std::vector<double>(m_upperRight.begin(), m_upperRight.end());
+  std::vector<ActsScalar> getMax() const final {
+    return std::vector<ActsScalar>(m_upperRight.begin(), m_upperRight.end());
   }
 
   /// @brief check whether given 3D position is inside look-up domain
diff --git a/Core/include/Acts/Surfaces/AnnulusBounds.hpp b/Core/include/Acts/Surfaces/AnnulusBounds.hpp
index f1283ec63b8..7c0d5840d2a 100644
--- a/Core/include/Acts/Surfaces/AnnulusBounds.hpp
+++ b/Core/include/Acts/Surfaces/AnnulusBounds.hpp
@@ -78,7 +78,7 @@ class AnnulusBounds : public DiscBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// Inside check for the bounds object driven by the boundary check directive
   /// Each Bounds has a method inside, which checks if a LocalPosition is inside
@@ -248,8 +248,8 @@ inline double AnnulusBounds::binningValuePhi() const {
   return get(eAveragePhi);
 }
 
-inline std::vector<double> AnnulusBounds::values() const {
-  std::vector<double> valvector;
+inline std::vector<ActsScalar> AnnulusBounds::values() const {
+  std::vector<ActsScalar> valvector;
   valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
   return valvector;
 }
diff --git a/Core/include/Acts/Surfaces/ConeBounds.hpp b/Core/include/Acts/Surfaces/ConeBounds.hpp
index 55b2203686c..ce17c6a89de 100644
--- a/Core/include/Acts/Surfaces/ConeBounds.hpp
+++ b/Core/include/Acts/Surfaces/ConeBounds.hpp
@@ -82,7 +82,7 @@ class ConeBounds : public SurfaceBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// inside method for local position
   ///
@@ -134,8 +134,8 @@ inline double ConeBounds::tanAlpha() const {
   return m_tanAlpha;
 }
 
-inline std::vector<double> ConeBounds::values() const {
-  std::vector<double> valvector;
+inline std::vector<ActsScalar> ConeBounds::values() const {
+  std::vector<ActsScalar> valvector;
   valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
   return valvector;
 }
diff --git a/Core/include/Acts/Surfaces/ConeSurface.hpp b/Core/include/Acts/Surfaces/ConeSurface.hpp
index a2b2572bba8..ee34d3a67a8 100644
--- a/Core/include/Acts/Surfaces/ConeSurface.hpp
+++ b/Core/include/Acts/Surfaces/ConeSurface.hpp
@@ -186,7 +186,7 @@ class ConeSurface : public RegularSurface {
       const Vector3& direction,
       const BoundaryTolerance& boundaryTolerance =
           BoundaryTolerance::Infinite(),
-      double tolerance = s_onSurfaceTolerance) const final;
+      ActsScalar tolerance = s_onSurfaceTolerance) const final;
 
   /// The pathCorrection for derived classes with thickness
   ///
diff --git a/Core/include/Acts/Surfaces/ConvexPolygonBounds.hpp b/Core/include/Acts/Surfaces/ConvexPolygonBounds.hpp
index e802056e7fd..b3de40bd350 100644
--- a/Core/include/Acts/Surfaces/ConvexPolygonBounds.hpp
+++ b/Core/include/Acts/Surfaces/ConvexPolygonBounds.hpp
@@ -40,7 +40,7 @@ class ConvexPolygonBoundsBase : public PlanarBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
  protected:
   /// Return a rectangle bounds instance that encloses a set of vertices.
diff --git a/Core/include/Acts/Surfaces/CylinderBounds.hpp b/Core/include/Acts/Surfaces/CylinderBounds.hpp
index 51508055e06..aa8fbed2b48 100644
--- a/Core/include/Acts/Surfaces/CylinderBounds.hpp
+++ b/Core/include/Acts/Surfaces/CylinderBounds.hpp
@@ -69,9 +69,9 @@ class CylinderBounds : public SurfaceBounds {
   /// @param avgPhi (optional) The phi value from which the opening angle spans
   /// @param bevelMinZ (optional) The bevel on the negative z side
   /// @param bevelMaxZ (optional) The bevel on the positive z sid The bevel on the positive z side
-  CylinderBounds(double r, double halfZ, double halfPhi = M_PI,
-                 double avgPhi = 0., double bevelMinZ = 0.,
-                 double bevelMaxZ = 0.) noexcept(false)
+  CylinderBounds(ActsScalar r, ActsScalar halfZ, ActsScalar halfPhi = M_PI,
+                 ActsScalar avgPhi = 0., ActsScalar bevelMinZ = 0.,
+                 ActsScalar bevelMaxZ = 0.) noexcept(false)
       : m_values({r, halfZ, halfPhi, avgPhi, bevelMinZ, bevelMaxZ}),
         m_closed(std::abs(halfPhi - M_PI) < s_epsilon) {
     checkConsistency();
@@ -80,7 +80,7 @@ class CylinderBounds : public SurfaceBounds {
   /// Constructor - from fixed size array
   ///
   /// @param values The parameter values
-  CylinderBounds(const std::array<double, eSize>& values) noexcept(false)
+  CylinderBounds(const std::array<ActsScalar, eSize>& values) noexcept(false)
       : m_values(values),
         m_closed(std::abs(values[eHalfPhiSector] - M_PI) < s_epsilon) {
     checkConsistency();
@@ -93,7 +93,7 @@ class CylinderBounds : public SurfaceBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// Inside check for the bounds object driven by the boundary check directive
   /// Each Bounds has a method inside, which checks if a LocalPosition is inside
@@ -107,7 +107,7 @@ class CylinderBounds : public SurfaceBounds {
 
   /// Access to the bound values
   /// @param bValue the class nested enum for the array access
-  double get(BoundValues bValue) const { return m_values[bValue]; }
+  ActsScalar get(BoundValues bValue) const { return m_values[bValue]; }
 
   /// Returns true for full phi coverage
   bool coversFullAzimuth() const;
@@ -130,7 +130,7 @@ class CylinderBounds : public SurfaceBounds {
 
  private:
   /// The bound radius, half Z, half phi and average phi
-  std::array<double, eSize> m_values;
+  std::array<ActsScalar, eSize> m_values;
   /// Indicator if the bounds are closed
   bool m_closed{false};
 
@@ -146,8 +146,8 @@ class CylinderBounds : public SurfaceBounds {
   ActsMatrix<2, 2> jacobian() const;
 };
 
-inline std::vector<double> CylinderBounds::values() const {
-  std::vector<double> valvector;
+inline std::vector<ActsScalar> CylinderBounds::values() const {
+  std::vector<ActsScalar> valvector;
   valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
   return valvector;
 }
diff --git a/Core/include/Acts/Surfaces/DiamondBounds.hpp b/Core/include/Acts/Surfaces/DiamondBounds.hpp
index 7497bd22309..4283493d7fe 100644
--- a/Core/include/Acts/Surfaces/DiamondBounds.hpp
+++ b/Core/include/Acts/Surfaces/DiamondBounds.hpp
@@ -78,7 +78,7 @@ class DiamondBounds : public PlanarBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// Inside check for the bounds object driven by the boundary check directive
   /// Each Bounds has a method inside, which checks if a LocalPosition is inside
@@ -119,8 +119,8 @@ class DiamondBounds : public PlanarBounds {
   void checkConsistency() noexcept(false);
 };
 
-inline std::vector<double> DiamondBounds::values() const {
-  std::vector<double> valvector;
+inline std::vector<ActsScalar> DiamondBounds::values() const {
+  std::vector<ActsScalar> valvector;
   valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
   return valvector;
 }
diff --git a/Core/include/Acts/Surfaces/DiscTrapezoidBounds.hpp b/Core/include/Acts/Surfaces/DiscTrapezoidBounds.hpp
index a3e4e5ea28f..c28987a9471 100644
--- a/Core/include/Acts/Surfaces/DiscTrapezoidBounds.hpp
+++ b/Core/include/Acts/Surfaces/DiscTrapezoidBounds.hpp
@@ -70,7 +70,7 @@ class DiscTrapezoidBounds : public DiscBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   ///  This method checks if the radius given in the LocalPosition is inside
   ///  [rMin,rMax]
@@ -207,8 +207,8 @@ inline double DiscTrapezoidBounds::binningValuePhi() const {
   return get(eAveragePhi);
 }
 
-inline std::vector<double> DiscTrapezoidBounds::values() const {
-  std::vector<double> valvector;
+inline std::vector<ActsScalar> DiscTrapezoidBounds::values() const {
+  std::vector<ActsScalar> valvector;
   valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
   return valvector;
 }
diff --git a/Core/include/Acts/Surfaces/EllipseBounds.hpp b/Core/include/Acts/Surfaces/EllipseBounds.hpp
index aec08d10381..dcb5560285b 100644
--- a/Core/include/Acts/Surfaces/EllipseBounds.hpp
+++ b/Core/include/Acts/Surfaces/EllipseBounds.hpp
@@ -78,7 +78,7 @@ class EllipseBounds : public PlanarBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// This method checks if the point given in the local coordinates is between
   /// two ellipsoids if only tol0 is given and additional in the phi sector is
@@ -119,8 +119,8 @@ class EllipseBounds : public PlanarBounds {
   void checkConsistency() noexcept(false);
 };
 
-inline std::vector<double> EllipseBounds::values() const {
-  std::vector<double> valvector;
+inline std::vector<ActsScalar> EllipseBounds::values() const {
+  std::vector<ActsScalar> valvector;
   valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
   return valvector;
 }
diff --git a/Core/include/Acts/Surfaces/InfiniteBounds.hpp b/Core/include/Acts/Surfaces/InfiniteBounds.hpp
index 261f40e4729..a41ce2d25ea 100644
--- a/Core/include/Acts/Surfaces/InfiniteBounds.hpp
+++ b/Core/include/Acts/Surfaces/InfiniteBounds.hpp
@@ -28,7 +28,7 @@ class InfiniteBounds : public SurfaceBounds {
     return SurfaceBounds::eBoundless;
   }
 
-  std::vector<double> values() const final { return {}; }
+  std::vector<ActsScalar> values() const final { return {}; }
 
   /// Method inside() returns true for any case
   ///
diff --git a/Core/include/Acts/Surfaces/LineBounds.hpp b/Core/include/Acts/Surfaces/LineBounds.hpp
index e3169d8a845..cf0135a2dac 100644
--- a/Core/include/Acts/Surfaces/LineBounds.hpp
+++ b/Core/include/Acts/Surfaces/LineBounds.hpp
@@ -51,7 +51,7 @@ class LineBounds : public SurfaceBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// Inside check for the bounds object driven by the boundary check directive
   /// Each Bounds has a method inside, which checks if a LocalPosition is inside
@@ -81,8 +81,8 @@ class LineBounds : public SurfaceBounds {
   void checkConsistency() noexcept(false);
 };
 
-inline std::vector<double> LineBounds::values() const {
-  std::vector<double> valvector;
+inline std::vector<ActsScalar> LineBounds::values() const {
+  std::vector<ActsScalar> valvector;
   valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
   return valvector;
 }
diff --git a/Core/include/Acts/Surfaces/RadialBounds.hpp b/Core/include/Acts/Surfaces/RadialBounds.hpp
index 0c7aea81400..5e83146a8b1 100644
--- a/Core/include/Acts/Surfaces/RadialBounds.hpp
+++ b/Core/include/Acts/Surfaces/RadialBounds.hpp
@@ -47,8 +47,8 @@ class RadialBounds : public DiscBounds {
   /// @param maxR The outer radius
   /// @param halfPhi The half opening angle (Pi for full angular coverage)
   /// @param avgPhi The average phi for the disc/ring sector
-  RadialBounds(double minR, double maxR, double halfPhi = M_PI,
-               double avgPhi = 0.) noexcept(false)
+  RadialBounds(ActsScalar minR, ActsScalar maxR, ActsScalar halfPhi = M_PI,
+               ActsScalar avgPhi = 0.) noexcept(false)
       : m_values({minR, maxR, halfPhi, avgPhi}) {
     checkConsistency();
   }
@@ -56,7 +56,7 @@ class RadialBounds : public DiscBounds {
   /// Constructor from array values
   ///
   /// @param values The bound values
-  RadialBounds(const std::array<double, eSize>& values) noexcept(false)
+  RadialBounds(const std::array<ActsScalar, eSize>& values) noexcept(false)
       : m_values(values) {
     checkConsistency();
   }
@@ -68,7 +68,7 @@ class RadialBounds : public DiscBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// For disc surfaces the local position in (r,phi) is checked
   ///
@@ -108,7 +108,7 @@ class RadialBounds : public DiscBounds {
   double binningValuePhi() const final;
 
  private:
-  std::array<double, eSize> m_values;
+  std::array<ActsScalar, eSize> m_values;
 
   /// Check the input values for consistency, will throw a logic_exception
   /// if consistency is not given
@@ -157,8 +157,8 @@ inline double RadialBounds::binningValuePhi() const {
   return get(eAveragePhi);
 }
 
-inline std::vector<double> RadialBounds::values() const {
-  std::vector<double> valvector;
+inline std::vector<ActsScalar> RadialBounds::values() const {
+  std::vector<ActsScalar> valvector;
   valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
   return valvector;
 }
diff --git a/Core/include/Acts/Surfaces/RectangleBounds.hpp b/Core/include/Acts/Surfaces/RectangleBounds.hpp
index aba1743fa74..f63a3fa189b 100644
--- a/Core/include/Acts/Surfaces/RectangleBounds.hpp
+++ b/Core/include/Acts/Surfaces/RectangleBounds.hpp
@@ -45,7 +45,7 @@ class RectangleBounds : public PlanarBounds {
   ///
   /// @param halfX halflength in X
   /// @param halfY halflength in Y
-  RectangleBounds(double halfX, double halfY) noexcept(false)
+  RectangleBounds(ActsScalar halfX, ActsScalar halfY) noexcept(false)
       : m_min({-halfX, -halfY}), m_max({halfX, halfY}) {
     checkConsistency();
   }
@@ -53,7 +53,7 @@ class RectangleBounds : public PlanarBounds {
   /// Constructor - from fixed size array - generic
   ///
   /// @param values The parameter values
-  RectangleBounds(const std::array<double, eSize>& values) noexcept(false)
+  RectangleBounds(const std::array<ActsScalar, eSize>& values) noexcept(false)
       : m_min({values[eMinX], values[eMinY]}),
         m_max({values[eMaxX], values[eMaxY]}) {
     checkConsistency();
@@ -72,7 +72,7 @@ class RectangleBounds : public PlanarBounds {
 
   BoundsType type() const final;
 
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// Inside check for the bounds object driven by the boundary check directive
   /// Each Bounds has a method inside, which checks if a LocalPosition is inside
@@ -148,7 +148,7 @@ inline double RectangleBounds::halfLengthY() const {
   return 0.5 * (m_max.y() - m_min.y());
 }
 
-inline std::vector<double> RectangleBounds::values() const {
+inline std::vector<ActsScalar> RectangleBounds::values() const {
   return {m_min.x(), m_min.y(), m_max.x(), m_max.y()};
 }
 
diff --git a/Core/include/Acts/Surfaces/SurfaceBounds.hpp b/Core/include/Acts/Surfaces/SurfaceBounds.hpp
index 62216f75c29..2f8d19ddf1c 100644
--- a/Core/include/Acts/Surfaces/SurfaceBounds.hpp
+++ b/Core/include/Acts/Surfaces/SurfaceBounds.hpp
@@ -56,7 +56,7 @@ class SurfaceBounds {
   /// vector containing the parameters needed to describe these bounds
   ///
   /// @return of the stored values for this SurfaceBounds object
-  virtual std::vector<double> values() const = 0;
+  virtual std::vector<ActsScalar> values() const = 0;
 
   /// Inside check for the bounds object driven by the boundary check directive
   /// Each Bounds has a method inside, which checks if a LocalPosition is inside
diff --git a/Core/include/Acts/Surfaces/TrapezoidBounds.hpp b/Core/include/Acts/Surfaces/TrapezoidBounds.hpp
index 8ccd3a5c928..172c41f6f64 100644
--- a/Core/include/Acts/Surfaces/TrapezoidBounds.hpp
+++ b/Core/include/Acts/Surfaces/TrapezoidBounds.hpp
@@ -48,19 +48,19 @@ class TrapezoidBounds : public PlanarBounds {
   /// @param halfXposY maximal half length X, definition at positive Y
   /// @param halfY half length Y - defined at x=0
   /// @param rotAngle: rotation angle of the bounds w.r.t coordinate axes
-  TrapezoidBounds(double halfXnegY, double halfXposY, double halfY,
-                  double rotAngle = 0.) noexcept(false);
+  TrapezoidBounds(ActsScalar halfXnegY, ActsScalar halfXposY, ActsScalar halfY,
+                  ActsScalar rotAngle = 0.) noexcept(false);
 
   /// Constructor for symmetric Trapezoid - from fixed size array
   ///
   /// @param values the values to be stream in
-  TrapezoidBounds(const std::array<double, eSize>& values) noexcept(false);
+  TrapezoidBounds(const std::array<ActsScalar, eSize>& values) noexcept(false);
 
   ~TrapezoidBounds() override;
 
   BoundsType type() const final;
 
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// The orientation of the Trapezoid is according to the figure above,
   /// in words: the shorter of the two parallel sides of the trapezoid
@@ -126,10 +126,10 @@ class TrapezoidBounds : public PlanarBounds {
 
   /// Access to the bound values
   /// @param bValue the class nested enum for the array access
-  double get(BoundValues bValue) const { return m_values[bValue]; }
+  ActsScalar get(BoundValues bValue) const { return m_values[bValue]; }
 
  private:
-  std::array<double, eSize> m_values;
+  std::array<ActsScalar, eSize> m_values;
   RectangleBounds m_boundingBox;
 
   void rotateBoundingBox() noexcept(false);
diff --git a/Core/include/Acts/TrackFinding/MeasurementSelector.hpp b/Core/include/Acts/TrackFinding/MeasurementSelector.hpp
index fd348d44075..d4092daa7c1 100644
--- a/Core/include/Acts/TrackFinding/MeasurementSelector.hpp
+++ b/Core/include/Acts/TrackFinding/MeasurementSelector.hpp
@@ -116,7 +116,8 @@ class MeasurementSelector {
   Result<Cuts> getCuts(const GeometryIdentifier& geoID, double theta) const;
 
   double calculateChi2(
-      const double* fullCalibrated, const double* fullCalibratedCovariance,
+      const ActsScalar* fullCalibrated,
+      const ActsScalar* fullCalibratedCovariance,
       TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
                        false>::Parameters predicted,
       TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
diff --git a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp
index 666ebf6aa23..41c3ea9424e 100644
--- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp
+++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp
@@ -240,7 +240,7 @@ class AtlasBetheHeitlerApprox {
     }
     // Return a component representation for higher x0
     // Cap the x because beyond the parameterization goes wild
-    const auto high_x = std::min(m_highLimit, x);
+    const auto high_x = std::min(m_highLimit, static_cast<double>(x));
     return make_mixture(m_highData, high_x, m_highTransform);
   }
 
diff --git a/Core/include/Acts/TrackFitting/GainMatrixUpdater.hpp b/Core/include/Acts/TrackFitting/GainMatrixUpdater.hpp
index 0ac45a79b1d..388d84a3554 100644
--- a/Core/include/Acts/TrackFitting/GainMatrixUpdater.hpp
+++ b/Core/include/Acts/TrackFitting/GainMatrixUpdater.hpp
@@ -27,8 +27,8 @@ class GainMatrixUpdater {
   struct InternalTrackState {
     unsigned int calibratedSize;
     // This is used to build a covariance matrix view in the .cpp file
-    const double* calibrated;
-    const double* calibratedCovariance;
+    const ActsScalar* calibrated;
+    const ActsScalar* calibratedCovariance;
     BoundSubspaceIndices projector;
 
     TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
diff --git a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp
index 92a6a4f76ba..423cbc0942f 100644
--- a/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp
+++ b/Core/include/Acts/TrackFitting/GlobalChiSquareFitter.hpp
@@ -298,11 +298,11 @@ struct ScatteringProperties {
 /// @param trackState The track state to analyse
 /// @param logger A logger instance
 template <std::size_t kMeasDim, typename track_state_t>
-void addMeasurementToGx2fSums(Eigen::MatrixXd& aMatrixExtended,
-                              Eigen::VectorXd& bVectorExtended, double& chi2sum,
-                              const std::vector<BoundMatrix>& jacobianFromStart,
-                              const track_state_t& trackState,
-                              const Logger& logger) {
+void addMeasurementToGx2fSums(
+    Eigen::Matrix<ActsScalar, Eigen::Dynamic, Eigen::Dynamic>& aMatrixExtended,
+    Eigen::Matrix<ActsScalar, Eigen::Dynamic, 1>& bVectorExtended,
+    double& chi2sum, const std::vector<BoundMatrix>& jacobianFromStart,
+    const track_state_t& trackState, const Logger& logger) {
   // First we get back the covariance and try to invert it. If the inversion
   // fails, we can already abort.
   const ActsSquareMatrix<kMeasDim> covarianceMeasurement =
@@ -321,8 +321,9 @@ void addMeasurementToGx2fSums(Eigen::MatrixXd& aMatrixExtended,
   const std::size_t dimsExtendedParams = aMatrixExtended.rows();
 
   // We create an empty Jacobian and fill it in the next steps
-  Eigen::MatrixXd extendedJacobian =
-      Eigen::MatrixXd::Zero(eBoundSize, dimsExtendedParams);
+  Eigen::Matrix<ActsScalar, Eigen::Dynamic, Eigen::Dynamic> extendedJacobian =
+      Eigen::Matrix<ActsScalar, Eigen::Dynamic, Eigen::Dynamic>::Zero(
+          eBoundSize, dimsExtendedParams);
 
   // This part of the Jacobian comes from the material-less propagation
   extendedJacobian.topLeftCorner<eBoundSize, eBoundSize>() =
@@ -352,7 +353,8 @@ void addMeasurementToGx2fSums(Eigen::MatrixXd& aMatrixExtended,
   const ActsMatrix<kMeasDim, eBoundSize> projector =
       trackState.projector().template topLeftCorner<kMeasDim, eBoundSize>();
 
-  const Eigen::MatrixXd projJacobian = projector * extendedJacobian;
+  const Eigen::Matrix<ActsScalar, Eigen::Dynamic, Eigen::Dynamic> projJacobian =
+      projector * extendedJacobian;
 
   const ActsMatrix<kMeasDim, 1> projPredicted = projector * predicted;
 
@@ -418,7 +420,8 @@ void addMeasurementToGx2fSums(Eigen::MatrixXd& aMatrixExtended,
 /// @param logger A logger instance
 template <typename track_state_t>
 void addMaterialToGx2fSums(
-    Eigen::MatrixXd& aMatrixExtended, Eigen::VectorXd& bVectorExtended,
+    Eigen::Matrix<ActsScalar, Eigen::Dynamic, Eigen::Dynamic>& aMatrixExtended,
+    Eigen::Matrix<ActsScalar, Eigen::Dynamic, 1>& bVectorExtended,
     double& chi2sum, const std::size_t nMaterialsHandled,
     const std::unordered_map<GeometryIdentifier, ScatteringProperties>&
         scatteringMap,
@@ -498,9 +501,10 @@ void addMaterialToGx2fSums(
 /// @param ndfSystem The number of degrees of freedom, determining the size of meaning full block
 ///
 /// @return deltaParams The calculated delta parameters.
-void updateGx2fCovarianceParams(BoundMatrix& fullCovariancePredicted,
-                                Eigen::MatrixXd& aMatrixExtended,
-                                const std::size_t ndfSystem);
+void updateGx2fCovarianceParams(
+    BoundMatrix& fullCovariancePredicted,
+    Eigen::Matrix<ActsScalar, Eigen::Dynamic, Eigen::Dynamic>& aMatrixExtended,
+    const std::size_t ndfSystem);
 
 /// Global Chi Square fitter (GX2F) implementation.
 ///
@@ -1180,10 +1184,13 @@ class Gx2Fitter {
 
       // Set to zero before filling
       chi2sum = 0;
-      Eigen::MatrixXd aMatrixExtended =
-          Eigen::MatrixXd::Zero(dimsExtendedParams, dimsExtendedParams);
-      Eigen::VectorXd bVectorExtended =
-          Eigen::VectorXd::Zero(dimsExtendedParams);
+      Eigen::Matrix<ActsScalar, Eigen::Dynamic, Eigen::Dynamic>
+          aMatrixExtended =
+              Eigen::Matrix<ActsScalar, Eigen::Dynamic, Eigen::Dynamic>::Zero(
+                  dimsExtendedParams, dimsExtendedParams);
+      Eigen::Matrix<ActsScalar, Eigen::Dynamic, 1> bVectorExtended =
+          Eigen::Matrix<ActsScalar, Eigen::Dynamic, 1>::Zero(
+              dimsExtendedParams);
 
       std::vector<BoundMatrix> jacobianFromStart;
       jacobianFromStart.emplace_back(BoundMatrix::Identity());
@@ -1300,7 +1307,7 @@ class Gx2Fitter {
       bVector = bVectorExtended.topLeftCorner<eBoundSize, 1>().eval();
 
       // calculate delta params [a] * delta = b
-      Eigen::VectorXd deltaParamsExtended =
+      Eigen::Matrix<ActsScalar, Eigen::Dynamic, 1> deltaParamsExtended =
           aMatrixExtended.colPivHouseholderQr().solve(bVectorExtended);
 
       deltaParams = deltaParamsExtended.topLeftCorner<eBoundSize, 1>().eval();
diff --git a/Core/include/Acts/TrackFitting/MbfSmoother.hpp b/Core/include/Acts/TrackFitting/MbfSmoother.hpp
index 02066a642e7..03dfaf10588 100644
--- a/Core/include/Acts/TrackFitting/MbfSmoother.hpp
+++ b/Core/include/Acts/TrackFitting/MbfSmoother.hpp
@@ -103,8 +103,8 @@ class MbfSmoother {
     struct Measurement final {
       unsigned int calibratedSize{0};
       // This is used to build a covariance matrix view in the .cpp file
-      const double* calibrated{nullptr};
-      const double* calibratedCovariance{nullptr};
+      const ActsScalar* calibrated{nullptr};
+      const ActsScalar* calibratedCovariance{nullptr};
       Projector projector;
 
       template <typename TrackStateProxy>
diff --git a/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp b/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp
index 501208b7e09..ba6ba5d64eb 100644
--- a/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp
+++ b/Core/include/Acts/TrackFitting/detail/GsfUtils.hpp
@@ -151,7 +151,7 @@ class ScopedGsfInfoPrinterAndChecker {
 };
 
 ActsScalar calculateDeterminant(
-    const double *fullCalibratedCovariance,
+    const ActsScalar *fullCalibratedCovariance,
     TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
                      true>::Covariance predictedCovariance,
     TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax, true>::Projector
diff --git a/Core/src/Geometry/CompositePortalLink.cpp b/Core/src/Geometry/CompositePortalLink.cpp
index d8010bbbce2..ea95d2b5a2f 100644
--- a/Core/src/Geometry/CompositePortalLink.cpp
+++ b/Core/src/Geometry/CompositePortalLink.cpp
@@ -212,7 +212,7 @@ std::unique_ptr<GridPortalLink> CompositePortalLink::makeGrid(
       throw std::runtime_error{"Unsupported binning direction"};
     }
 
-    std::vector<double> edges;
+    std::vector<ActsScalar> edges;
     edges.reserve(m_children.size() + 1);
 
     const Transform3& groupTransform = m_surface->transform(gctx);
@@ -256,7 +256,7 @@ std::unique_ptr<GridPortalLink> CompositePortalLink::makeGrid(
       throw std::runtime_error{"Unsupported binning direction"};
     }
 
-    std::vector<double> edges;
+    std::vector<ActsScalar> edges;
     edges.reserve(m_children.size() + 1);
 
     std::ranges::sort(trivialLinks, [](const auto& a, const auto& b) {
diff --git a/Core/src/Geometry/CuboidVolumeBuilder.cpp b/Core/src/Geometry/CuboidVolumeBuilder.cpp
index 614ce28a67f..75d03496484 100644
--- a/Core/src/Geometry/CuboidVolumeBuilder.cpp
+++ b/Core/src/Geometry/CuboidVolumeBuilder.cpp
@@ -136,10 +136,12 @@ std::pair<double, double> Acts::CuboidVolumeBuilder::binningRange(
   }
 
   // Use the volume boundaries as limits for the binning
-  minMax.first = std::min(
-      minMax.first, minVolumeBoundaries(toUnderlying(cfg.binningDimension)));
-  minMax.second = std::max(
-      minMax.second, maxVolumeBoundaries(toUnderlying(cfg.binningDimension)));
+  minMax.first =
+      std::min(minMax.first, static_cast<double>(minVolumeBoundaries(
+                                 toUnderlying(cfg.binningDimension))));
+  minMax.second =
+      std::max(minMax.second, static_cast<double>(maxVolumeBoundaries(
+                                  toUnderlying(cfg.binningDimension))));
 
   return minMax;
 }
diff --git a/Core/src/Geometry/CutoutCylinderVolumeBounds.cpp b/Core/src/Geometry/CutoutCylinderVolumeBounds.cpp
index 9d70f49b4b9..9b6368d63b3 100644
--- a/Core/src/Geometry/CutoutCylinderVolumeBounds.cpp
+++ b/Core/src/Geometry/CutoutCylinderVolumeBounds.cpp
@@ -27,7 +27,7 @@
 #include <utility>
 
 bool Acts::CutoutCylinderVolumeBounds::inside(const Acts::Vector3& gpos,
-                                              double tol) const {
+                                              Acts::ActsScalar tol) const {
   // first check whether we are in the outer envelope at all (ignore r_med)
   using VectorHelpers::perp;
   using VectorHelpers::phi;
diff --git a/Core/src/Geometry/CylinderVolumeBounds.cpp b/Core/src/Geometry/CylinderVolumeBounds.cpp
index d9af620ea45..7a6ca13cd4e 100644
--- a/Core/src/Geometry/CylinderVolumeBounds.cpp
+++ b/Core/src/Geometry/CylinderVolumeBounds.cpp
@@ -98,16 +98,15 @@ std::vector<OrientedSurface> CylinderVolumeBounds::orientedSurfaces(
   if (bevelMinZ != 0.) {
     ActsScalar sy = 1 - 1 / std::cos(bevelMinZ);
     transMinZ = transform * vMinZ *
-                Eigen::AngleAxisd(-bevelMinZ, Eigen::Vector3d(1., 0., 0.)) *
-                Eigen::Scaling(1., 1. + sy, 1.);
+                AngleAxis3(-bevelMinZ, Vector3(1., 0., 0.)) *
+                Eigen::DiagonalMatrix<ActsScalar, 3>(1., 1. + sy, 1.);
   } else {
     transMinZ = transform * vMinZ;
   }
   if (bevelMaxZ != 0.) {
     ActsScalar sy = 1 - 1 / std::cos(bevelMaxZ);
-    transMaxZ = transform * vMaxZ *
-                Eigen::AngleAxisd(bevelMaxZ, Eigen::Vector3d(1., 0., 0.)) *
-                Eigen::Scaling(1., 1. + sy, 1.);
+    transMaxZ = transform * vMaxZ * AngleAxis3(bevelMaxZ, Vector3(1., 0., 0.)) *
+                Eigen::DiagonalMatrix<ActsScalar, 3>(1., 1. + sy, 1.);
   } else {
     transMaxZ = transform * vMaxZ;
   }
diff --git a/Core/src/Geometry/CylinderVolumeBuilder.cpp b/Core/src/Geometry/CylinderVolumeBuilder.cpp
index 47de379d9be..1a7713c2706 100644
--- a/Core/src/Geometry/CylinderVolumeBuilder.cpp
+++ b/Core/src/Geometry/CylinderVolumeBuilder.cpp
@@ -614,14 +614,18 @@ Acts::VolumeConfig Acts::CylinderVolumeBuilder::analyzeContent(
       const CylinderVolumeBounds* cvBounds =
           dynamic_cast<const CylinderVolumeBounds*>(&volume->volumeBounds());
       if (cvBounds != nullptr) {
-        lConfig.rMin =
-            std::min(lConfig.rMin, cvBounds->get(CylinderVolumeBounds::eMinR));
-        lConfig.rMax =
-            std::max(lConfig.rMax, cvBounds->get(CylinderVolumeBounds::eMaxR));
-        lConfig.zMin = std::min(
-            lConfig.zMin, -cvBounds->get(CylinderVolumeBounds::eHalfLengthZ));
-        lConfig.zMax = std::max(
-            lConfig.zMax, cvBounds->get(CylinderVolumeBounds::eHalfLengthZ));
+        lConfig.rMin = std::min(
+            lConfig.rMin,
+            static_cast<double>(cvBounds->get(CylinderVolumeBounds::eMinR)));
+        lConfig.rMax = std::max(
+            lConfig.rMax,
+            static_cast<double>(cvBounds->get(CylinderVolumeBounds::eMaxR)));
+        lConfig.zMin =
+            std::min(lConfig.zMin, static_cast<double>(-cvBounds->get(
+                                       CylinderVolumeBounds::eHalfLengthZ)));
+        lConfig.zMax =
+            std::max(lConfig.zMax, static_cast<double>(cvBounds->get(
+                                       CylinderVolumeBounds::eHalfLengthZ)));
       }
     }
   }
diff --git a/Core/src/Geometry/Extent.cpp b/Core/src/Geometry/Extent.cpp
index 2f07f6358bc..3f04fffa734 100644
--- a/Core/src/Geometry/Extent.cpp
+++ b/Core/src/Geometry/Extent.cpp
@@ -43,7 +43,7 @@ void Acts::Extent::extend(const Vector3& vtx,
     ActsScalar mValue = cValue - lEnv;
     // Special protection for radial value
     if (bValue == BinningValue::binR && mValue < 0.) {
-      mValue = std::max(mValue, 0.);
+      mValue = std::max(mValue, 0_scalar);
     }
     if (constrains(bValue)) {
       m_range[toUnderlying(bValue)].expand(mValue, cValue + hEnv);
diff --git a/Core/src/Geometry/GenericCuboidVolumeBounds.cpp b/Core/src/Geometry/GenericCuboidVolumeBounds.cpp
index bdf8f960dc2..1ee64f2ee62 100644
--- a/Core/src/Geometry/GenericCuboidVolumeBounds.cpp
+++ b/Core/src/Geometry/GenericCuboidVolumeBounds.cpp
@@ -34,7 +34,8 @@ Acts::GenericCuboidVolumeBounds::GenericCuboidVolumeBounds(
 }
 
 Acts::GenericCuboidVolumeBounds::GenericCuboidVolumeBounds(
-    const std::array<double, GenericCuboidVolumeBounds::BoundValues::eSize>&
+    const std::array<Acts::ActsScalar,
+                     GenericCuboidVolumeBounds::BoundValues::eSize>&
         values) noexcept(false)
     : m_vertices() {
   for (std::size_t iv = 0; iv < 8; ++iv) {
@@ -45,12 +46,12 @@ Acts::GenericCuboidVolumeBounds::GenericCuboidVolumeBounds(
 }
 
 bool Acts::GenericCuboidVolumeBounds::inside(const Acts::Vector3& gpos,
-                                             double tol) const {
+                                             Acts::ActsScalar tol) const {
   constexpr std::array<std::size_t, 6> vtxs = {0, 4, 0, 1, 2, 1};
   // needs to be on same side, get ref
   bool ref = std::signbit((gpos - m_vertices[vtxs[0]]).dot(m_normals[0]));
   for (std::size_t i = 1; i < 6; i++) {
-    double dot = (gpos - m_vertices[vtxs[i]]).dot(m_normals[i]);
+    Acts::ActsScalar dot = (gpos - m_vertices[vtxs[i]]).dot(m_normals[i]);
     if (std::signbit(dot) != ref) {
       // technically outside, but how far?
       if (std::abs(dot) > tol) {
@@ -187,8 +188,8 @@ void Acts::GenericCuboidVolumeBounds::construct() noexcept(false) {
   handle_face(m_vertices[1], m_vertices[0], m_vertices[4], m_vertices[5]);
 }
 
-std::vector<double> Acts::GenericCuboidVolumeBounds::values() const {
-  std::vector<double> rvalues;
+std::vector<Acts::ActsScalar> Acts::GenericCuboidVolumeBounds::values() const {
+  std::vector<Acts::ActsScalar> rvalues;
   rvalues.reserve(BoundValues::eSize);
   for (std::size_t iv = 0; iv < 8; ++iv) {
     for (std::size_t ic = 0; ic < 3; ++ic) {
diff --git a/Core/src/Geometry/Portal.cpp b/Core/src/Geometry/Portal.cpp
index 69cab13b8e0..675a3148ef6 100644
--- a/Core/src/Geometry/Portal.cpp
+++ b/Core/src/Geometry/Portal.cpp
@@ -324,8 +324,8 @@ bool Portal::isSameSurface(const GeometryContext& gctx, const Surface& a,
     return false;
   }
 
-  std::vector<double> aValues = a.bounds().values();
-  std::vector<double> bValues = b.bounds().values();
+  std::vector<ActsScalar> aValues = a.bounds().values();
+  std::vector<ActsScalar> bValues = b.bounds().values();
   bool different = false;
   for (auto [aVal, bVal] : zip(aValues, bValues)) {
     if (std::abs(aVal - bVal) > s_onSurfaceTolerance) {
diff --git a/Core/src/Propagator/SympyStepper.cpp b/Core/src/Propagator/SympyStepper.cpp
index a33e1f567bb..62365c8fff6 100644
--- a/Core/src/Propagator/SympyStepper.cpp
+++ b/Core/src/Propagator/SympyStepper.cpp
@@ -110,12 +110,12 @@ Result<double> SympyStepper::stepImpl(
     double stepSizeCutOff, std::size_t maxRungeKuttaStepTrials) const {
   auto pos = position(state);
   auto dir = direction(state);
-  double t = time(state);
-  double qop = qOverP(state);
-  double m = particleHypothesis(state).mass();
-  double p_abs = absoluteMomentum(state);
+  ActsScalar t = time(state);
+  ActsScalar qop = qOverP(state);
+  ActsScalar m = particleHypothesis(state).mass();
+  ActsScalar p_abs = absoluteMomentum(state);
 
-  auto getB = [&](const double* p) -> Result<Vector3> {
+  auto getB = [&](const ActsScalar* p) -> Result<Vector3> {
     return getField(state, {p[0], p[1], p[2]});
   };
 
@@ -138,10 +138,10 @@ Result<double> SympyStepper::stepImpl(
     return std::clamp(x, lower, upper);
   };
 
-  double h = state.stepSize.value() * stepDirection;
+  ActsScalar h = state.stepSize.value() * stepDirection;
   double initialH = h;
   std::size_t nStepTrials = 0;
-  double errorEstimate = 0.;
+  ActsScalar errorEstimate = 0.;
 
   while (true) {
     nStepTrials++;
@@ -149,7 +149,8 @@ Result<double> SympyStepper::stepImpl(
     // For details about the factor 4 see ATL-SOFT-PUB-2009-001
     Result<bool> res =
         rk4(pos.data(), dir.data(), t, h, qop, m, p_abs, getB, &errorEstimate,
-            4 * stepTolerance, state.pars.template segment<3>(eFreePos0).data(),
+            ActsScalar{4 * stepTolerance},
+            state.pars.template segment<3>(eFreePos0).data(),
             state.pars.template segment<3>(eFreeDir0).data(),
             state.pars.template segment<1>(eFreeTime).data(),
             state.derivative.data(),
@@ -158,13 +159,13 @@ Result<double> SympyStepper::stepImpl(
       return res.error();
     }
     // Protect against division by zero
-    errorEstimate = std::max(1e-20, errorEstimate);
+    errorEstimate = std::max(1e-20_scalar, errorEstimate);
 
     if (*res) {
       break;
     }
 
-    const double stepSizeScaling = calcStepSizeScaling(errorEstimate);
+    const ActsScalar stepSizeScaling = calcStepSizeScaling(errorEstimate);
     h *= stepSizeScaling;
 
     // If step size becomes too small the particle remains at the initial
@@ -186,10 +187,10 @@ Result<double> SympyStepper::stepImpl(
   ++state.nSteps;
   state.nStepTrials += nStepTrials;
 
-  const double stepSizeScaling = calcStepSizeScaling(errorEstimate);
-  const double nextAccuracy = std::abs(h * stepSizeScaling);
-  const double previousAccuracy = std::abs(state.stepSize.accuracy());
-  const double initialStepLength = std::abs(initialH);
+  const ActsScalar stepSizeScaling = calcStepSizeScaling(errorEstimate);
+  const ActsScalar nextAccuracy = std::abs(h * stepSizeScaling);
+  const ActsScalar previousAccuracy = std::abs(state.stepSize.accuracy());
+  const ActsScalar initialStepLength = std::abs(initialH);
   if (nextAccuracy < initialStepLength || nextAccuracy > previousAccuracy) {
     state.stepSize.setAccuracy(nextAccuracy);
   }
diff --git a/Core/src/Surfaces/ConvexPolygonBounds.cpp b/Core/src/Surfaces/ConvexPolygonBounds.cpp
index f6b08624308..79fd129edbb 100644
--- a/Core/src/Surfaces/ConvexPolygonBounds.cpp
+++ b/Core/src/Surfaces/ConvexPolygonBounds.cpp
@@ -30,8 +30,8 @@ std::ostream& Acts::ConvexPolygonBoundsBase::toStream(std::ostream& sl) const {
   return sl;
 }
 
-std::vector<double> Acts::ConvexPolygonBoundsBase::values() const {
-  std::vector<double> values;
+std::vector<Acts::ActsScalar> Acts::ConvexPolygonBoundsBase::values() const {
+  std::vector<Acts::ActsScalar> values;
   for (const auto& vtx : vertices()) {
     values.push_back(vtx.x());
     values.push_back(vtx.y());
diff --git a/Core/src/Surfaces/TrapezoidBounds.cpp b/Core/src/Surfaces/TrapezoidBounds.cpp
index 3450e8cea97..15645ec4a61 100644
--- a/Core/src/Surfaces/TrapezoidBounds.cpp
+++ b/Core/src/Surfaces/TrapezoidBounds.cpp
@@ -22,9 +22,9 @@
 /// @param halfXposY maximal half length X, definition at positive Y
 /// @param halfY half length Y - defined at x=0
 /// @param rotAngle: rotation angle of the bounds w.r.t coordinate axes
-Acts::TrapezoidBounds::TrapezoidBounds(double halfXnegY, double halfXposY,
-                                       double halfY,
-                                       double rotAngle) noexcept(false)
+Acts::TrapezoidBounds::TrapezoidBounds(
+    Acts::ActsScalar halfXnegY, Acts::ActsScalar halfXposY,
+    Acts::ActsScalar halfY, Acts::ActsScalar rotAngle) noexcept(false)
     : m_values({halfXnegY, halfXposY, halfY, rotAngle}),
       m_boundingBox(std::max(halfXnegY, halfXposY), halfY) {
   rotateBoundingBox();
@@ -35,7 +35,7 @@ Acts::TrapezoidBounds::TrapezoidBounds(double halfXnegY, double halfXposY,
 ///
 /// @param values the values to be stream in
 Acts::TrapezoidBounds::TrapezoidBounds(
-    const std::array<double, eSize>& values) noexcept(false)
+    const std::array<Acts::ActsScalar, eSize>& values) noexcept(false)
     : m_values(values),
       m_boundingBox(
           std::max(values[eHalfLengthXnegY], values[eHalfLengthXposY]),
@@ -57,19 +57,20 @@ bool Acts::TrapezoidBounds::inside(
     return true;
   }
 
-  const double hlXnY = get(TrapezoidBounds::eHalfLengthXnegY);
-  const double hlXpY = get(TrapezoidBounds::eHalfLengthXposY);
-  const double hlY = get(TrapezoidBounds::eHalfLengthY);
-  const double rotAngle = get(TrapezoidBounds::eRotationAngle);
+  const Acts::ActsScalar hlXnY = get(TrapezoidBounds::eHalfLengthXnegY);
+  const Acts::ActsScalar hlXpY = get(TrapezoidBounds::eHalfLengthXposY);
+  const Acts::ActsScalar hlY = get(TrapezoidBounds::eHalfLengthY);
+  const Acts::ActsScalar rotAngle = get(TrapezoidBounds::eRotationAngle);
 
-  const Acts::Vector2 extPosition = Eigen::Rotation2Dd(rotAngle) * lposition;
-  const double x = extPosition[0];
-  const double y = extPosition[1];
+  const Acts::Vector2 extPosition =
+      Eigen::Rotation2D<Acts::ActsScalar>(rotAngle) * lposition;
+  const Acts::ActsScalar x = extPosition[0];
+  const Acts::ActsScalar y = extPosition[1];
 
   if (auto absoluteBound = boundaryTolerance.asAbsoluteBoundOpt(true);
       absoluteBound.has_value()) {
-    double tolX = absoluteBound->tolerance0;
-    double tolY = absoluteBound->tolerance1;
+    Acts::ActsScalar tolX = absoluteBound->tolerance0;
+    Acts::ActsScalar tolY = absoluteBound->tolerance1;
 
     if (std::abs(y) - hlY > tolY) {
       // outside y range
@@ -97,15 +98,15 @@ bool Acts::TrapezoidBounds::inside(
 
 std::vector<Acts::Vector2> Acts::TrapezoidBounds::vertices(
     unsigned int /*ignoredSegments*/) const {
-  const double hlXnY = get(TrapezoidBounds::eHalfLengthXnegY);
-  const double hlXpY = get(TrapezoidBounds::eHalfLengthXposY);
-  const double hlY = get(TrapezoidBounds::eHalfLengthY);
-  const double rotAngle = get(TrapezoidBounds::eRotationAngle);
+  const Acts::ActsScalar hlXnY = get(TrapezoidBounds::eHalfLengthXnegY);
+  const Acts::ActsScalar hlXpY = get(TrapezoidBounds::eHalfLengthXposY);
+  const Acts::ActsScalar hlY = get(TrapezoidBounds::eHalfLengthY);
+  const Acts::ActsScalar rotAngle = get(TrapezoidBounds::eRotationAngle);
 
   std::vector<Acts::Vector2> vertices = {
       {-hlXnY, -hlY}, {hlXnY, -hlY}, {hlXpY, hlY}, {-hlXpY, hlY}};
   for (auto& v : vertices) {
-    v = Eigen::Rotation2Dd(-rotAngle) * v;
+    v = Eigen::Rotation2D<Acts::ActsScalar>(-rotAngle) * v;
   }
   return vertices;
 }
@@ -124,14 +125,14 @@ std::ostream& Acts::TrapezoidBounds::toStream(std::ostream& sl) const {
   return sl;
 }
 
-std::vector<double> Acts::TrapezoidBounds::values() const {
-  std::vector<double> valvector;
+std::vector<Acts::ActsScalar> Acts::TrapezoidBounds::values() const {
+  std::vector<Acts::ActsScalar> valvector;
   valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
   return valvector;
 }
 
 void Acts::TrapezoidBounds::rotateBoundingBox() noexcept(false) {
-  const double rotAngle = get(eRotationAngle);
+  const Acts::ActsScalar rotAngle = get(eRotationAngle);
 
   if (rotAngle != 0.) {
     m_boundingBox = ConvexPolygonBounds<4>(vertices()).boundingBox();
diff --git a/Core/src/Surfaces/detail/AnnulusBoundsHelper.cpp b/Core/src/Surfaces/detail/AnnulusBoundsHelper.cpp
index 68e22e4c72b..c90782d2e6e 100644
--- a/Core/src/Surfaces/detail/AnnulusBoundsHelper.cpp
+++ b/Core/src/Surfaces/detail/AnnulusBoundsHelper.cpp
@@ -17,7 +17,7 @@ std::tuple<std::shared_ptr<Acts::AnnulusBounds>, Acts::Transform3>
 Acts::detail::AnnulusBoundsHelper::create(const Transform3& transform,
                                           ActsScalar rMin, ActsScalar rMax,
                                           std::vector<Vector2> vertices) {
-  using Line2D = Eigen::Hyperplane<double, 2>;
+  using Line2D = Eigen::Hyperplane<Acts::ActsScalar, 2>;
 
   // Construct the bound lines
   std::vector<std::pair<Vector2, Vector2>> boundLines;
@@ -45,7 +45,7 @@ Acts::detail::AnnulusBoundsHelper::create(const Transform3& transform,
   Line2D lB = Line2D::Through(boundLines[1].first, boundLines[1].second);
   Vector2 ix = lA.intersection(lB);
 
-  const Eigen::Translation3d originTranslation(ix.x(), ix.y(), 0.);
+  const Translation3 originTranslation(ix.x(), ix.y(), 0);
   const Vector2 originShift = -ix;
 
   // Update transform by prepending the origin shift translation
@@ -53,9 +53,9 @@ Acts::detail::AnnulusBoundsHelper::create(const Transform3& transform,
   // Transform phi line point to new origin and get phi
   double phi1 = VectorHelpers::phi(boundLines[0].second - boundLines[0].first);
   double phi2 = VectorHelpers::phi(boundLines[1].second - boundLines[1].first);
-  double phiMax = std::max(phi1, phi2);
-  double phiMin = std::min(phi1, phi2);
-  double phiShift = 0.;
+  ActsScalar phiMax = std::max(phi1, phi2);
+  ActsScalar phiMin = std::min(phi1, phi2);
+  ActsScalar phiShift{0.};
 
   // Create the bounds
   auto annulusBounds = std::make_shared<AnnulusBounds>(
diff --git a/Core/src/Surfaces/detail/MergeHelper.cpp b/Core/src/Surfaces/detail/MergeHelper.cpp
index 8d0b1fa0b5a..c4eca991a3a 100644
--- a/Core/src/Surfaces/detail/MergeHelper.cpp
+++ b/Core/src/Surfaces/detail/MergeHelper.cpp
@@ -8,6 +8,8 @@
 
 #include "Acts/Surfaces/detail/MergeHelper.hpp"
 
+#include <numbers>
+
 namespace Acts::detail {
 
 std::tuple<ActsScalar, ActsScalar, bool> mergedPhiSector(
@@ -15,23 +17,25 @@ std::tuple<ActsScalar, ActsScalar, bool> mergedPhiSector(
     ActsScalar avgPhi2, const Logger& logger, ActsScalar tolerance) {
   using namespace Acts::UnitLiterals;
 
-  if (std::abs(hlPhi1 - M_PI / 2.0) < tolerance &&
-      std::abs(hlPhi2 - M_PI / 2.0) < tolerance) {
+  if (std::abs(hlPhi1 - ActsScalar{std::numbers::pi / 2.}) < tolerance &&
+      std::abs(hlPhi2 - ActsScalar{std::numbers::pi / 2.}) < tolerance) {
     ACTS_VERBOSE("Both phi sectors cover a half circle");
 
     ACTS_VERBOSE("-> distance between sectors: "
-                 << detail::difference_periodic(avgPhi1, avgPhi2, 2 * M_PI) /
+                 << detail::difference_periodic(
+                        avgPhi1, avgPhi2, ActsScalar{2. * std::numbers::pi}) /
                         1_degree);
 
-    if (std::abs(
-            std::abs(detail::difference_periodic(avgPhi1, avgPhi2, 2 * M_PI)) -
-            M_PI) > tolerance) {
+    if (std::abs(std::abs(detail::difference_periodic(
+                     avgPhi1, avgPhi2, ActsScalar{2 * std::numbers::pi})) -
+                 std::numbers::pi_v<ActsScalar>) > tolerance) {
       throw std::invalid_argument(
           "Phi sectors cover half a circle but are not opposite");
     }
 
-    ActsScalar newAvgPhi = detail::radian_sym(avgPhi1 + M_PI / 2.0);
-    ActsScalar newHlPhi = M_PI;
+    ActsScalar newAvgPhi =
+        detail::radian_sym(avgPhi1 + ActsScalar{std::numbers::pi / 2.});
+    ActsScalar newHlPhi = std::numbers::pi_v<ActsScalar>;
     ACTS_VERBOSE("merged: ["
                  << detail::radian_sym(newAvgPhi - newHlPhi) / 1_degree << ", "
                  << detail::radian_sym(newAvgPhi + newHlPhi) / 1_degree
diff --git a/Core/src/TrackFinding/MeasurementSelector.cpp b/Core/src/TrackFinding/MeasurementSelector.cpp
index 130c8381536..1ab7be9f289 100644
--- a/Core/src/TrackFinding/MeasurementSelector.cpp
+++ b/Core/src/TrackFinding/MeasurementSelector.cpp
@@ -73,7 +73,8 @@ MeasurementSelector::InternalCutBins MeasurementSelector::convertCutBins(
 }
 
 double MeasurementSelector::calculateChi2(
-    const double* fullCalibrated, const double* fullCalibratedCovariance,
+    const ActsScalar* fullCalibrated,
+    const ActsScalar* fullCalibratedCovariance,
     TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
                      false>::Parameters predicted,
     TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax,
diff --git a/Core/src/TrackFitting/GlobalChiSquareFitter.cpp b/Core/src/TrackFitting/GlobalChiSquareFitter.cpp
index 0784322308b..3d012f7ed30 100644
--- a/Core/src/TrackFitting/GlobalChiSquareFitter.cpp
+++ b/Core/src/TrackFitting/GlobalChiSquareFitter.cpp
@@ -11,7 +11,8 @@
 #include "Acts/Definitions/TrackParametrization.hpp"
 
 void Acts::Experimental::updateGx2fCovarianceParams(
-    BoundMatrix& fullCovariancePredicted, Eigen::MatrixXd& aMatrixExtended,
+    BoundMatrix& fullCovariancePredicted,
+    Eigen::Matrix<ActsScalar, Eigen::Dynamic, Eigen::Dynamic>& aMatrixExtended,
     const std::size_t ndfSystem) {
   // make invertible
   for (int i = 0; i < aMatrixExtended.rows(); ++i) {
diff --git a/Core/src/TrackFitting/GsfUtils.cpp b/Core/src/TrackFitting/GsfUtils.cpp
index ae5b47ca700..a4a2870deaa 100644
--- a/Core/src/TrackFitting/GsfUtils.cpp
+++ b/Core/src/TrackFitting/GsfUtils.cpp
@@ -18,7 +18,7 @@ using TrackStateTraits =
     TrackStateTraits<MultiTrajectoryTraits::MeasurementSizeMax, true>;
 
 ActsScalar calculateDeterminant(
-    const double* fullCalibratedCovariance,
+    const ActsScalar* fullCalibratedCovariance,
     TrackStateTraits::Covariance predictedCovariance,
     TrackStateTraits::Projector projector, unsigned int calibratedSize) {
   return visit_measurement(calibratedSize, [&](auto N) {
diff --git a/Core/src/Vertexing/NumericalTrackLinearizer.cpp b/Core/src/Vertexing/NumericalTrackLinearizer.cpp
index 035e42ab3c8..60839bb5cb2 100644
--- a/Core/src/Vertexing/NumericalTrackLinearizer.cpp
+++ b/Core/src/Vertexing/NumericalTrackLinearizer.cpp
@@ -13,6 +13,8 @@
 #include "Acts/Utilities/UnitVectors.hpp"
 #include "Acts/Vertexing/LinearizerTrackParameters.hpp"
 
+#include <numbers>
+
 Acts::Result<Acts::LinearizedTrack>
 Acts::NumericalTrackLinearizer::linearizeTrack(
     const BoundTrackParameters& params, double linPointTime,
@@ -141,7 +143,8 @@ Acts::NumericalTrackLinearizer::linearizeTrack(
     // previously computed value for better readability.
     completeJacobian(eLinPhi, i) =
         Acts::detail::difference_periodic(newPerigeeParams(eLinPhi),
-                                          perigeeParams(eLinPhi), 2 * M_PI) /
+                                          perigeeParams(eLinPhi),
+                                          ActsScalar{2. * std::numbers::pi}) /
         m_cfg.delta;
   }
 
diff --git a/Core/src/Visualization/ObjVisualization3D.cpp b/Core/src/Visualization/ObjVisualization3D.cpp
index 777e883b3d3..67f387cc354 100644
--- a/Core/src/Visualization/ObjVisualization3D.cpp
+++ b/Core/src/Visualization/ObjVisualization3D.cpp
@@ -60,7 +60,9 @@ void ObjVisualization3D::faces(const std::vector<Vector3>& vtxs,
     if (color != Color{0, 0, 0}) {
       o.vertexColors[o.vertices.size()] = color;
     }
-    o.vertices.insert(o.vertices.end(), vtxs.begin(), vtxs.end());
+    std::ranges::transform(
+        vtxs, std::back_inserter(o.vertices),
+        [](const Vector3& v) { return v.template cast<double>(); });
     for (const auto& face : faces) {
       if (face.size() == 2) {
         o.lines.push_back({face[0] + vtxoffs, face[2] + vtxoffs});
diff --git a/docs/core/definitions/algebra.rst b/docs/core/definitions/algebra.rst
index 9d35ef34992..7beb81f7ad7 100644
--- a/docs/core/definitions/algebra.rst
+++ b/docs/core/definitions/algebra.rst
@@ -2,12 +2,12 @@ Algebra definitions
 ===================
 
 The main algebra classes for ACTS are defined in the `Acts/Definitions/Algebra.hpp` header file.
-The basic scalar type can be defined via this file and is set per default to `double`, however, if `ACTS_CUSTOM_SCALAR` is set it will be used instead.
+The basic scalar type can be defined via this file and is set per default to `double`, however, if `ACTS_CUSTOM_SCALARTYPE` is set it will be used instead.
 
 .. code-block:: cpp
 
-   #ifdef ACTS_CUSTOM_SCALAR
-   using ActsScalar = ACTS_CUSTOM_SCALAR;
+   #ifdef ACTS_CUSTOM_SCALARTYPE
+   using ActsScalar = ACTS_CUSTOM_SCALARTYPE;
    #else
    using ActsScalar = double;
    #endif