diff --git a/Core/include/ACTS/EventData/detail/coordinate_transformations.hpp b/Core/include/ACTS/EventData/detail/coordinate_transformations.hpp
index 5d04f0b48b8925d546de61f0f6238906e9687b14..a3c2e5fe5ab0ed1cfecc885d8d5a0d64c4be8223 100644
--- a/Core/include/ACTS/EventData/detail/coordinate_transformations.hpp
+++ b/Core/include/ACTS/EventData/detail/coordinate_transformations.hpp
@@ -12,6 +12,7 @@
 // ACTS includes
 #include "ACTS/Surfaces/Surface.hpp"
 #include "ACTS/Utilities/ParameterDefinitions.hpp"
+#include <cmath>
 
 #ifdef ACTS_COORDINATE_TRANSFORM_PLUGIN
 
@@ -43,7 +44,7 @@ namespace detail {
     parameters2globalMomentum(const ParVector_t& pars)
     {
       ActsVectorD<3> momentum;
-      double         p     = fabs(1. / pars(Acts::eQOP));
+      double         p     = std::abs(1. / pars(Acts::eQOP));
       double         phi   = pars(Acts::ePHI);
       double         theta = pars(Acts::eTHETA);
       momentum << p * sin(theta) * cos(phi), p * sin(theta) * sin(phi),
@@ -59,7 +60,7 @@ namespace detail {
     {
       ParVector_t parameters;
       parameters << 0, 0, mom.phi(), mom.theta(),
-          ((fabs(charge) < 1e-4) ? 1. : charge) / mom.mag();
+          ((std::abs(charge) < 1e-4) ? 1. : charge) / mom.mag();
 
       return parameters;
     }
@@ -74,7 +75,7 @@ namespace detail {
       s.globalToLocal(pos, mom, localPosition);
       ParVector_t result;
       result << localPosition(0), localPosition(1), mom.phi(), mom.theta(),
-          ((fabs(charge) < 1e-4) ? 1. : charge) / mom.mag();
+          ((std::abs(charge) < 1e-4) ? 1. : charge) / mom.mag();
       return result;
     }
 
diff --git a/Core/include/ACTS/Extrapolation/detail/RungeKuttaEngine.ipp b/Core/include/ACTS/Extrapolation/detail/RungeKuttaEngine.ipp
index 4e9c1dd9c866708660aaf797a5760a0df42cebd2..92ed478719606ab29bc49afc7c5eb04f60d6e453 100644
--- a/Core/include/ACTS/Extrapolation/detail/RungeKuttaEngine.ipp
+++ b/Core/include/ACTS/Extrapolation/detail/RungeKuttaEngine.ipp
@@ -13,6 +13,7 @@
 #include "ACTS/Surfaces/PlaneSurface.hpp"
 #include "ACTS/Surfaces/StrawSurface.hpp"
 #include "ACTS/Surfaces/Surface.hpp"
+#include <cmath>
 
 template <class MagneticField>
 template <class T>
@@ -258,7 +259,7 @@ Acts::RungeKuttaEngine<MagneticField>::propagate(
           pCache.pVector[0], pCache.pVector[1], pCache.pVector[2]);
       Acts::Vector3D mom(
           pCache.pVector[3], pCache.pVector[4], pCache.pVector[5]);
-      mom /= fabs(pCache.parameters[4]);
+      mom /= std::abs(pCache.parameters[4]);
       nParameters = std::make_unique<NeutralCurvilinearParameters>(
           std::move(cov), gp, mom);
     }
@@ -470,7 +471,7 @@ Acts::RungeKuttaEngine<MagneticField>::propagateWithJacobian(
   SA[0] = SA[1] = SA[2] = 0.;
   pCache.maxPathLimit   = false;
 
-  if (pCache.mcondition && fabs(pCache.pVector[6]) > 100.) return false;
+  if (pCache.mcondition && std::abs(pCache.pVector[6]) > 100.) return false;
 
   // Step estimation until surface
   bool   Q;
@@ -484,7 +485,7 @@ Acts::RungeKuttaEngine<MagneticField>::propagateWithJacobian(
   }
 
   step > Smax ? S = Smax : step < -Smax ? S = -Smax : S = step;
-  double So                                             = fabs(S);
+  double So                                             = std::abs(S);
   int    iS                                             = 0;
 
   bool InS = false;
@@ -494,7 +495,7 @@ Acts::RungeKuttaEngine<MagneticField>::propagateWithJacobian(
   pCache.newfield = true;
 
   // whie loop over the steps
-  while (fabs(step) > m_cfg.straightStep) {
+  while (std::abs(step) > m_cfg.straightStep) {
     // maximum number of steps
     if (++pCache.niter > 10000) {
       //!< @todo make max number configurable
@@ -541,14 +542,14 @@ Acts::RungeKuttaEngine<MagneticField>::propagateWithJacobian(
     }
 
     // check if the step made sense
-    double aS    = fabs(S);
-    double aStep = fabs(step);
+    double aS    = std::abs(S);
+    double aStep = std::abs(step);
     if (aS > aStep)
       S = step;
     else if (!iS && InS && aS * 2. < aStep)
       S *= 2.;
 
-    if (!dir && fabs(pCache.step) > Wwrong) {
+    if (!dir && std::abs(pCache.step) > Wwrong) {
       EX_MSG_DEBUG(navigationStep,
                    "propagate",
                    "<T> ",
@@ -556,20 +557,20 @@ Acts::RungeKuttaEngine<MagneticField>::propagateWithJacobian(
       return false;
     }
 
-    if (iS > 10 || (iS > 3 && fabs(S) >= So)) {
+    if (iS > 10 || (iS > 3 && std::abs(S) >= So)) {
       if (!kind) break;
       EX_MSG_DEBUG(navigationStep, "propagate", "<T> ", "Abort triggered.");
       return false;
     }
 
-    double dW = pCache.maxPathLength - fabs(pCache.step);
-    if (fabs(S) > dW) {
+    double dW = pCache.maxPathLength - std::abs(pCache.step);
+    if (std::abs(S) > dW) {
       S > 0. ? S = dW : S = -dW;
       step                = S;
       pCache.maxPathLimit = true;
     }
 
-    So = fabs(S);
+    So = std::abs(S);
 
   }  // end of while loop
 
@@ -579,7 +580,7 @@ Acts::RungeKuttaEngine<MagneticField>::propagateWithJacobian(
   // Output track parameteres
   pCache.step += step;
 
-  if (fabs(step) < .001) return true;
+  if (std::abs(step) < .001) return true;
 
   A[0] += (SA[0] * step);
   A[1] += (SA[1] * step);
@@ -630,7 +631,7 @@ Acts::RungeKuttaEngine<MagneticField>::rungeKuttaStep(int navigationStep,
   }
 
   bool Helix                           = false;
-  if (fabs(S) < m_cfg.helixStep) Helix = true;
+  if (std::abs(S) < m_cfg.helixStep) Helix = true;
 
   while (S != 0.) {
     double S3 = (1. / 3.) * S, S4 = .25 * S, PS2 = Pi * S;
@@ -688,8 +689,8 @@ Acts::RungeKuttaEngine<MagneticField>::rungeKuttaStep(int navigationStep,
 
     // Test approximation quality on give step and possible step reduction
     //
-    double EST = fabs((A1 + A6) - (A3 + A4)) + fabs((B1 + B6) - (B3 + B4))
-        + fabs((C1 + C6) - (C3 + C4));
+    double EST = std::abs((A1 + A6) - (A3 + A4)) + std::abs((B1 + B6) - (B3 + B4))
+        + std::abs((C1 + C6) - (C3 + C4));
     if (EST > m_cfg.dlt) {
       S *= .5;
       dltm = 0.;
@@ -899,8 +900,8 @@ Acts::RungeKuttaEngine<MagneticField>::rungeKuttaStepWithGradient(
 
     // Test approximation quality on give step and possible step reduction
     //
-    double EST = fabs((A1 + A6) - (A3 + A4)) + fabs((B1 + B6) - (B3 + B4))
-        + fabs((C1 + C6) - (C3 + C4));
+    double EST = std::abs((A1 + A6) - (A3 + A4)) + std::abs((B1 + B6) - (B3 + B4))
+        + std::abs((C1 + C6) - (C3 + C4));
     if (EST > m_cfg.dlt) {
       S *= .5;
       dltm = 0.;
@@ -1073,7 +1074,7 @@ Acts::RungeKuttaEngine<MagneticField>::newCrossPoint(const CylinderSurface& Su,
   z               = P[2] - T(2, 3);
   RC              = x * Ax[0] + y * Ax[1] + z * Ax[2];
   RS              = x * Ay[0] + y * Ay[1] + z * Ay[2];
-  double dF       = fabs(atan2(RS, RC) - Su.bounds().averagePhi());
+  double dF       = std::abs(atan2(RS, RC) - Su.bounds().averagePhi());
   if (dF > pi) dF = pi2 - pi;
   if (dF <= Su.bounds().halfPhiSector()) return false;
   return true;
@@ -1180,7 +1181,7 @@ Acts::RungeKuttaEngine<MagneticField>::stepEstimatorWithCurvature(
   // Straight step estimation
   double Step = m_rkUtils.stepEstimator(kind, Su, pCache.pVector, Q);
   if (!Q) return 0.;
-  double AStep = fabs(Step);
+  double AStep = std::abs(Step);
   if (kind || AStep < m_cfg.straightStep || !pCache.mcondition) return Step;
 
   const double* SA = &(pCache.pVector[42]);  // Start direction
@@ -1202,6 +1203,6 @@ Acts::RungeKuttaEngine<MagneticField>::stepEstimatorWithCurvature(
     Q = true;
     return Step;
   }
-  if (fabs(StepN) < AStep) return StepN;
+  if (std::abs(StepN) < AStep) return StepN;
   return Step;
 }
diff --git a/Core/include/ACTS/Propagator/AtlasStepper.hpp b/Core/include/ACTS/Propagator/AtlasStepper.hpp
index 93ce90fcc5fa072cd51e82567390b3f51ea331ba..c9286fa2927d20a627638cf2e01974e004e79587 100644
--- a/Core/include/ACTS/Propagator/AtlasStepper.hpp
+++ b/Core/include/ACTS/Propagator/AtlasStepper.hpp
@@ -80,7 +80,7 @@ class AtlasStepper
       pVector[4] = Sf * Se;  // Ay
       pVector[5] = Ce;       // Az
       pVector[6] = Vp[4];    // CM
-      if (fabs(pVector[6]) < .000000000000001) {
+      if (std::abs(pVector[6]) < .000000000000001) {
         pVector[6] < 0. ? pVector[6] = -.000000000000001 : pVector[6]
             = .000000000000001;
       }
@@ -156,7 +156,7 @@ public:
     double         charge = cache.pVector[6] > 0. ? 1. : -1.;
     Acts::Vector3D gp(cache.pVector[0], cache.pVector[1], cache.pVector[2]);
     Acts::Vector3D mom(cache.pVector[3], cache.pVector[4], cache.pVector[5]);
-    mom /= fabs(cache.pVector[6]);
+    mom /= std::abs(cache.pVector[6]);
 
     double P[45];
     for (unsigned int i = 0; i < 45; ++i) P[i] = cache.pVector[i];
@@ -289,7 +289,7 @@ public:
     double         charge = cache.pVector[6] > 0. ? 1. : -1.;
     Acts::Vector3D gp(cache.pVector[0], cache.pVector[1], cache.pVector[2]);
     Acts::Vector3D mom(cache.pVector[3], cache.pVector[4], cache.pVector[5]);
-    mom /= fabs(cache.pVector[6]);
+    mom /= std::abs(cache.pVector[6]);
 
     std::unique_ptr<ActsSymMatrixD<5>> cov = nullptr;
     if (cache.covariance) {
@@ -464,7 +464,7 @@ public:
     }
 
     bool Helix = false;
-    // if (fabs(S) < m_cfg.helixStep) Helix = true;
+    // if (std::abs(S) < m_cfg.helixStep) Helix = true;
 
     while (stepMax != 0.) {
       double S3 = (1. / 3.) * stepMax, S4 = .25 * stepMax, PS2 = Pi * stepMax;
@@ -523,8 +523,8 @@ public:
 
       // Test approximation quality on give step and possible step reduction
       //
-      double EST = fabs((A1 + A6) - (A3 + A4)) + fabs((B1 + B6) - (B3 + B4))
-          + fabs((C1 + C6) - (C3 + C4));
+      double EST = std::abs((A1 + A6) - (A3 + A4)) + std::abs((B1 + B6) - (B3 + B4))
+          + std::abs((C1 + C6) - (C3 + C4));
       if (EST > 0.0002) {
         stepMax *= .5;
         //        dltm = 0.;
diff --git a/Core/include/ACTS/Propagator/EigenStepper.hpp b/Core/include/ACTS/Propagator/EigenStepper.hpp
index f1ff52bc0ab922eb879350eda76184510a8c84df..eccd138889886207eb210cb92bedaed96d39d9b2 100644
--- a/Core/include/ACTS/Propagator/EigenStepper.hpp
+++ b/Core/include/ACTS/Propagator/EigenStepper.hpp
@@ -13,6 +13,7 @@
 #include "ACTS/Surfaces/Surface.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/Units.hpp"
+#include <cmath>
 
 namespace Acts {
 
@@ -134,7 +135,7 @@ public:
       const double theta = c.dir.theta();
 
       ActsMatrixD<5, 7> J = ActsMatrixD<5, 7>::Zero();
-      if (fabs(cos(theta)) < 0.99) {
+      if (std::abs(cos(theta)) < 0.99) {
         J(0, 0) = -sin(phi);
         J(0, 1) = cos(phi);
         J(1, 0) = -cos(phi) * cos(theta);
@@ -179,7 +180,7 @@ public:
     }
 
     return CurvilinearParameters(
-        std::move(cov), c.pos, c.dir / fabs(c.qop), charge);
+        std::move(cov), c.pos, c.dir / std::abs(c.qop), charge);
   }
 
   template <typename S>
@@ -222,7 +223,7 @@ public:
     }
 
     return BoundParameters(
-        std::move(cov), c.pos, c.dir / fabs(c.qop), charge, surface);
+        std::move(cov), c.pos, c.dir / std::abs(c.qop), charge, surface);
   }
 
   static double
diff --git a/Core/include/ACTS/Propagator/Propagator.hpp b/Core/include/ACTS/Propagator/Propagator.hpp
index b90490df3469734e73dd672f5a2b96c8b2f09890..71c313b1b41ec81824fd7f8524bbb65e64c70184 100644
--- a/Core/include/ACTS/Propagator/Propagator.hpp
+++ b/Core/include/ACTS/Propagator/Propagator.hpp
@@ -10,6 +10,7 @@
 #define ACTS_EXTRAPOLATION_PROPAGATOR_H 1
 
 #include <memory>
+#include <cmath>
 #include <type_traits>
 #include "ACTS/Propagator/AbortList.hpp"
 #include "ACTS/Propagator/ObserverList.hpp"
@@ -181,7 +182,7 @@ namespace propagation {
         r.pathLength += m_impl.step(propagation_cache, stepMax);
         return_parameter_type current = m_impl.convert(propagation_cache);
         options.observer_list(current, previous, r);
-        if (fabs(r.pathLength) >= options.max_path_length
+        if (std::abs(r.pathLength) >= options.max_path_length
             || options.stop_conditions(r, current, stepMax)) {
           r.endParameters = std::make_unique<return_parameter_type>(
               m_impl.convert(propagation_cache));
@@ -189,7 +190,7 @@ namespace propagation {
           break;
         }
 
-        if (fabs(stepMax) > fabs(signed_pathLimit - r.pathLength))
+        if (std::abs(stepMax) > std::abs(signed_pathLimit - r.pathLength))
           stepMax = signed_pathLimit - r.pathLength;
 
         previous = current;
@@ -237,15 +238,15 @@ namespace propagation {
       }
 
       double stepMax = options.direction * options.max_step_size;
-      if (fabs(stepMax) > fabs(distance)) stepMax = distance;
+      if (std::abs(stepMax) > std::abs(distance)) stepMax = distance;
 
       for (; r.steps < options.max_steps; ++r.steps) {
         r.pathLength += m_impl.step(cache, stepMax);
         step_parameter_type current = m_impl.convert(cache);
         options.observer_list(current, previous, r);
         distance = m_impl.distance(target, cache.position(), cache.direction());
-        if (fabs(distance) < 1 * units::_um
-            || fabs(r.pathLength) >= options.max_path_length
+        if (std::abs(distance) < 1 * units::_um
+            || std::abs(r.pathLength) >= options.max_path_length
             || options.stop_conditions(r, current, stepMax)) {
           r.endParameters = std::make_unique<return_parameter_type>(
               m_impl.convert(cache, target));
@@ -253,10 +254,10 @@ namespace propagation {
           break;
         }
 
-        if (fabs(stepMax) > fabs(signed_pathLimit - r.pathLength))
+        if (std::abs(stepMax) > std::abs(signed_pathLimit - r.pathLength))
           stepMax = signed_pathLimit - r.pathLength;
 
-        if (fabs(stepMax) > fabs(distance)) stepMax = distance;
+        if (std::abs(stepMax) > std::abs(distance)) stepMax = distance;
 
         previous = current;
       }
diff --git a/Core/include/ACTS/Surfaces/BoundaryCheck.hpp b/Core/include/ACTS/Surfaces/BoundaryCheck.hpp
index 9a8360794f23ea33ca4ab94c8871279f9cc0762b..9dca2be13ba2fc03f565435b4722319102c700b0 100644
--- a/Core/include/ACTS/Surfaces/BoundaryCheck.hpp
+++ b/Core/include/ACTS/Surfaces/BoundaryCheck.hpp
@@ -153,7 +153,7 @@ BoundaryCheck::FastArcTan(double x) const
     x          = 1. / x;  // keep arg between 0 and 1
     complement = true;
   }
-  y = M_PI_4 * x - x * (fabs(x) - 1) * (0.2447 + 0.0663 * fabs(x));
+  y = M_PI_4 * x - x * (std::abs(x) - 1) * (0.2447 + 0.0663 * std::abs(x));
   if (complement) y = M_PI_2 - y;  // correct for 1/x if we did that
   if (sign) y       = -y;          // correct for negative arg
   return y;
diff --git a/Core/include/ACTS/Surfaces/ConeBounds.hpp b/Core/include/ACTS/Surfaces/ConeBounds.hpp
index d93e9ef917df1f207fe6570ea00b51721688b43e..8e14d1cde06f291477893ba725d8249d6eeb9f5b 100644
--- a/Core/include/ACTS/Surfaces/ConeBounds.hpp
+++ b/Core/include/ACTS/Surfaces/ConeBounds.hpp
@@ -266,7 +266,7 @@ ConeBounds::insideLoc1(const Vector2D& lpos, double tol1) const
 inline double
 ConeBounds::r(double z) const
 {
-  return fabs(z * m_tanAlpha);
+  return std::abs(z * m_tanAlpha);
 }
 
 inline double
diff --git a/Core/include/ACTS/Surfaces/CylinderBounds.hpp b/Core/include/ACTS/Surfaces/CylinderBounds.hpp
index 3c22cf082432392fc43db8b70c66f8eea754165e..8234601a63694a007694fce6a1f3bcf9015d7d7b 100644
--- a/Core/include/ACTS/Surfaces/CylinderBounds.hpp
+++ b/Core/include/ACTS/Surfaces/CylinderBounds.hpp
@@ -15,6 +15,7 @@
 
 #include "ACTS/Surfaces/SurfaceBounds.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
+#include <cmath>
 
 namespace Acts {
 
@@ -259,15 +260,15 @@ CylinderBounds::inside(double r,
 inline bool
 CylinderBounds::insideLocZ(double z, double tol1) const
 {
-  return (m_valueStore.at(CylinderBounds::bv_halfZ) + tol1) - fabs(z) > 0.;
+  return (m_valueStore.at(CylinderBounds::bv_halfZ) + tol1) - std::abs(z) > 0.;
 }
 
 inline bool
 CylinderBounds::insideLoc0(const Vector2D& lpos, double tol0) const
 {
   bool insideRphi = false;
-  if (fabs(m_valueStore.at(CylinderBounds::bv_averagePhi)) < 10e-7)
-    insideRphi = (fabs(lpos[Acts::eLOC_RPHI]
+  if (std::abs(m_valueStore.at(CylinderBounds::bv_averagePhi)) < 10e-7)
+    insideRphi = (std::abs(lpos[Acts::eLOC_RPHI]
                        / m_valueStore.at(CylinderBounds::bv_radius))
                   < (m_valueStore.at(CylinderBounds::bv_halfPhiSector) + tol0));
   else {
diff --git a/Core/include/ACTS/Surfaces/CylinderSurface.hpp b/Core/include/ACTS/Surfaces/CylinderSurface.hpp
index 9ff76506e93bddcae157972adbe4cb7ef583494c..0cd644652ab57444548df74324aa5ef16c40c935 100644
--- a/Core/include/ACTS/Surfaces/CylinderSurface.hpp
+++ b/Core/include/ACTS/Surfaces/CylinderSurface.hpp
@@ -16,6 +16,7 @@
 #include "ACTS/Surfaces/CylinderBounds.hpp"
 #include "ACTS/Surfaces/Surface.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
+#include <cmath>
 
 namespace Acts {
 
@@ -270,7 +271,7 @@ CylinderSurface::pathCorrection(const Vector3D& gpos, const Vector3D& mom) const
 {
   Vector3D normalT  = normal(gpos);
   double   cosAlpha = normalT.dot(mom.unit());
-  return fabs(1. / cosAlpha);
+  return std::abs(1. / cosAlpha);
 }
 
 inline const CylinderBounds&
diff --git a/Core/include/ACTS/Surfaces/DiamondBounds.hpp b/Core/include/ACTS/Surfaces/DiamondBounds.hpp
index 74e51dc0a136d8a2746ed59f019c66109ce80655..ffa248fe82b3b1bc67b9f03bc7a8b95553874dbc 100644
--- a/Core/include/ACTS/Surfaces/DiamondBounds.hpp
+++ b/Core/include/ACTS/Surfaces/DiamondBounds.hpp
@@ -13,7 +13,7 @@
 #ifndef ACTS_SURFACES_DIAMONDDBOUNDS_H
 #define ACTS_SURFACES_DIAMONDDBOUNDS_H 1
 
-#include <math.h>
+#include <cmath>
 #include "ACTS/Surfaces/PlanarBounds.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/ParameterDefinitions.hpp"
@@ -241,7 +241,7 @@ DiamondBounds::inside(const Vector2D& lpos, const BoundaryCheck& bcheck) const
       > 2 * m_valueStore.at(DiamondBounds::bv_halfY2) + limit)
     return false;
   // a fast FALSE
-  double fabsX = fabs(lpos[Acts::eLOC_X]);
+  double fabsX = std::abs(lpos[Acts::eLOC_X]);
   if (fabsX > (m_valueStore.at(DiamondBounds::bv_medHalfX) + limit))
     return false;
   // a fast TRUE
@@ -254,7 +254,7 @@ DiamondBounds::inside(const Vector2D& lpos, const BoundaryCheck& bcheck) const
                - limit))
     return true;
   // a fast TRUE
-  if (fabs(lpos[Acts::eLOC_Y])
+  if (std::abs(lpos[Acts::eLOC_Y])
       < (fmin(m_valueStore.at(DiamondBounds::bv_halfY1),
               m_valueStore.at(DiamondBounds::bv_halfY2))
          - limit))
@@ -309,7 +309,7 @@ DiamondBounds::inside(const Vector2D& lpos, const BoundaryCheck& bcheck) const
 inline bool
 DiamondBounds::insideLoc0(const Vector2D& lpos, double tol0) const
 {
-  return (fabs(lpos[Acts::eLOC_X])
+  return (std::abs(lpos[Acts::eLOC_X])
           < m_valueStore.at(DiamondBounds::bv_medHalfX) + tol0);
 }
 
diff --git a/Core/include/ACTS/Surfaces/DiscSurface.hpp b/Core/include/ACTS/Surfaces/DiscSurface.hpp
index b47fe80b3580e02d23a567c8518aa91626533873..552a4f9c62dd2e902a75fe999a05ac7cd57abd81 100644
--- a/Core/include/ACTS/Surfaces/DiscSurface.hpp
+++ b/Core/include/ACTS/Surfaces/DiscSurface.hpp
@@ -338,7 +338,7 @@ inline double
 DiscSurface::pathCorrection(const Vector3D&, const Vector3D& mom) const
 {
   /// we can ignore the global position here
-  return 1. / fabs(normal().dot(mom.unit()));
+  return 1. / std::abs(normal().dot(mom.unit()));
 }
 
 inline Intersection
diff --git a/Core/include/ACTS/Surfaces/DiscTrapezoidalBounds.hpp b/Core/include/ACTS/Surfaces/DiscTrapezoidalBounds.hpp
index 74fbb6aa721283adbd167fcb983aec000b903056..a832a22ad744196fbd4d318d2ba4d7174e103934 100644
--- a/Core/include/ACTS/Surfaces/DiscTrapezoidalBounds.hpp
+++ b/Core/include/ACTS/Surfaces/DiscTrapezoidalBounds.hpp
@@ -14,7 +14,6 @@
 #define ACTS_SURFACES_DISCTRAPEZOIDALBOUNDS_H
 
 #include <cmath>
-#include <math.h>
 #include "ACTS/Surfaces/DiscBounds.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/ParameterDefinitions.hpp"
@@ -184,14 +183,14 @@ DiscTrapezoidalBounds::inside(const Vector2D& lpos, double, double) const
                     Pos[Acts::eLOC_Y] * sin(phi)
                         + Pos[Acts::eLOC_X] * cos(phi));
 
-  bool insideX = (fabs(DeltaPos[Acts::eLOC_X])
+  bool insideX = (std::abs(DeltaPos[Acts::eLOC_X])
                   < m_valueStore.at(DiscTrapezoidalBounds::bv_maxHalfX));
-  bool insideY = (fabs(DeltaPos[Acts::eLOC_Y])
+  bool insideY = (std::abs(DeltaPos[Acts::eLOC_Y])
                   < m_valueStore.at(DiscTrapezoidalBounds::bv_halfY));
 
   if (!insideX || !insideY) return false;
 
-  if (fabs(DeltaPos[Acts::eLOC_X])
+  if (std::abs(DeltaPos[Acts::eLOC_X])
       < m_valueStore.at(DiscTrapezoidalBounds::bv_minHalfX))
     return true;
 
@@ -206,7 +205,7 @@ DiscTrapezoidalBounds::inside(const Vector2D& lpos, double, double) const
          - m_valueStore.at(DiscTrapezoidalBounds::bv_maxHalfX));
 
   bool inside
-      = (DeltaPos[Acts::eLOC_Y] > (m * fabs(DeltaPos[Acts::eLOC_X]) + q));
+      = (DeltaPos[Acts::eLOC_Y] > (m * std::abs(DeltaPos[Acts::eLOC_X]) + q));
 
   return inside;
 }
@@ -219,7 +218,7 @@ DiscTrapezoidalBounds::inside(const Vector2D&      lpos,
       || m_valueStore.at(DiscTrapezoidalBounds::bv_rMin) != 0)
     return DiscTrapezoidalBounds::inside(
         lpos, bcheck.toleranceLoc0, bcheck.toleranceLoc1);
-  double alpha = fabs(lpos[Acts::eLOC_PHI]
+  double alpha = std::abs(lpos[Acts::eLOC_PHI]
                       - m_valueStore.at(DiscTrapezoidalBounds::bv_averagePhi));
   if (alpha > M_PI) alpha = 2 * M_PI - alpha;
   // a fast FALSE
@@ -357,8 +356,8 @@ DiscTrapezoidalBounds::inside(const Vector2D&      lpos,
             double y1,
             double r) const
     {
-      double x = fabs(x1 - x0);
-      double y = fabs(y1 - y0);
+      double x = std::abs(x1 - x0);
+      double y = std::abs(y1 - y0);
       if (x * x + (h - y) * (h - y) <= r * r
           || (w - x) * (w - x) + y * y <= r * r
           || x * h + y * w <= w * h  // collision with (0, h)
@@ -419,7 +418,7 @@ DiscTrapezoidalBounds::inside(const Vector2D&      lpos,
 inline bool
 DiscTrapezoidalBounds::insideLoc0(const Vector2D& lpos, double tol0) const
 {
-  double alpha = fabs(lpos[Acts::eLOC_PHI]
+  double alpha = std::abs(lpos[Acts::eLOC_PHI]
                       - m_valueStore.at(DiscTrapezoidalBounds::bv_averagePhi));
   if (alpha > M_PI) alpha = 2 * M_PI - alpha;
 
@@ -439,7 +438,7 @@ DiscTrapezoidalBounds::insideLoc0(const Vector2D& lpos, double tol0) const
 inline bool
 DiscTrapezoidalBounds::insideLoc1(const Vector2D& lpos, double tol1) const
 {
-  double alpha = fabs(lpos[Acts::eLOC_PHI]
+  double alpha = std::abs(lpos[Acts::eLOC_PHI]
                       - m_valueStore.at(DiscTrapezoidalBounds::bv_averagePhi));
   if (alpha > M_PI) alpha = 2. * M_PI - alpha;
   return (alpha
diff --git a/Core/include/ACTS/Surfaces/EllipseBounds.hpp b/Core/include/ACTS/Surfaces/EllipseBounds.hpp
index 25e0e135a609db862821cf80137a9f2610ee2103..4f6558a77f2509095e48071a74eb7dbd0fe904c2 100644
--- a/Core/include/ACTS/Surfaces/EllipseBounds.hpp
+++ b/Core/include/ACTS/Surfaces/EllipseBounds.hpp
@@ -13,8 +13,8 @@
 #ifndef ACTS_SURFACES_ELLIPSEBOUNDS_H
 #define ACTS_SURFACES_ELLIPSEBOUNDS_H 1
 
-#include <math.h>
-#include <stdlib.h>
+#include <cmath>
+#include <cstdlib>
 #include "ACTS/Surfaces/PlanarBounds.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
 
diff --git a/Core/include/ACTS/Surfaces/LineBounds.hpp b/Core/include/ACTS/Surfaces/LineBounds.hpp
index 6fdd4719c3f7c119d5924a4cbee81cb232d5d658..ffb3f63dc0748d02eaa69d24b0be4884bc5ca821 100644
--- a/Core/include/ACTS/Surfaces/LineBounds.hpp
+++ b/Core/include/ACTS/Surfaces/LineBounds.hpp
@@ -199,13 +199,13 @@ LineBounds::inside(const Vector2D& lpos, double tol0, double tol1) const
 inline bool
 LineBounds::insideLocR(double r, double tol0) const
 {
-  return fabs(m_valueStore.at(LineBounds::bv_radius) - r) < tol0;
+  return std::abs(m_valueStore.at(LineBounds::bv_radius) - r) < tol0;
 }
 
 inline bool
 LineBounds::insideLocZ(double z, double tol1) const
 {
-  return (m_valueStore.at(LineBounds::bv_halfZ) + tol1) - fabs(z) > 0.;
+  return (m_valueStore.at(LineBounds::bv_halfZ) + tol1) - std::abs(z) > 0.;
 }
 
 inline double
diff --git a/Core/include/ACTS/Surfaces/PlaneSurface.hpp b/Core/include/ACTS/Surfaces/PlaneSurface.hpp
index 260b27c12e5645b502d597d3647833c14e0b8877..80ce32fb7eff3df4bb139fc4edfb9b52f183d5d4 100644
--- a/Core/include/ACTS/Surfaces/PlaneSurface.hpp
+++ b/Core/include/ACTS/Surfaces/PlaneSurface.hpp
@@ -247,7 +247,7 @@ inline double
 PlaneSurface::pathCorrection(const Vector3D&, const Vector3D& mom) const
 {
   /// we can ignore the global position here
-  return 1. / fabs(normal().dot(mom.unit()));
+  return 1. / std::abs(normal().dot(mom.unit()));
 }
 
 inline Intersection
diff --git a/Core/include/ACTS/Surfaces/RadialBounds.hpp b/Core/include/ACTS/Surfaces/RadialBounds.hpp
index c34b316a66ad120e2d3c2177d1ab417dfea9c0d1..95d4f553e7b71075e5db1d0d58e1dde79041c9b2 100644
--- a/Core/include/ACTS/Surfaces/RadialBounds.hpp
+++ b/Core/include/ACTS/Surfaces/RadialBounds.hpp
@@ -14,7 +14,6 @@
 #define ACTS_SURFACES_RADIALBOUNDS_H 1
 
 #include <cmath>
-#include <math.h>
 #include "ACTS/Surfaces/DiscBounds.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/ParameterDefinitions.hpp"
@@ -170,7 +169,7 @@ inline bool
 RadialBounds::inside(const Vector2D& lpos, double tol0, double tol1) const
 {
   double alpha
-      = fabs(lpos[Acts::eLOC_PHI] - m_valueStore[RadialBounds::bv_averagePhi]);
+      = std::abs(lpos[Acts::eLOC_PHI] - m_valueStore[RadialBounds::bv_averagePhi]);
   if (alpha > M_PI) alpha = 2 * M_PI - alpha;
   bool insidePhi
       = (alpha <= (m_valueStore[RadialBounds::bv_halfPhiSector] + tol1));
@@ -313,8 +312,8 @@ RadialBounds::inside(const Vector2D& lpos, const BoundaryCheck& bcheck) const
             double y1,
             double r) const
     {
-      double x = fabs(x1 - x0);
-      double y = fabs(y1 - y0);
+      double x = std::abs(x1 - x0);
+      double y = std::abs(y1 - y0);
       if (x * x + (h - y) * (h - y) <= r * r
           || (w - x) * (w - x) + y * y <= r * r
           || x * h + y * w <= w * h  // collision with (0, h)
@@ -383,7 +382,7 @@ inline bool
 RadialBounds::insideLoc1(const Vector2D& lpos, double tol1) const
 {
   double alpha
-      = fabs(lpos[Acts::eLOC_PHI] - m_valueStore[RadialBounds::bv_averagePhi]);
+      = std::abs(lpos[Acts::eLOC_PHI] - m_valueStore[RadialBounds::bv_averagePhi]);
   if (alpha > M_PI) alpha = 2. * M_PI - alpha;
   // alpha -= alpha > M_PI ? 2.*M_PI : 0.;
   // alpha += alpha < -M_PI ? 2.*M_PI : 0.;
diff --git a/Core/include/ACTS/Surfaces/RectangleBounds.hpp b/Core/include/ACTS/Surfaces/RectangleBounds.hpp
index e5f8af1de8898a369a47d5b52f1049db16b9b197..ca3efe4d180182b79308145d6f681d83491fa9c3 100644
--- a/Core/include/ACTS/Surfaces/RectangleBounds.hpp
+++ b/Core/include/ACTS/Surfaces/RectangleBounds.hpp
@@ -210,14 +210,14 @@ RectangleBounds::inside(const Vector2D& lpos, double tol0, double tol1) const
 inline bool
 RectangleBounds::insideLoc0(const Vector2D& lpos, double tol0) const
 {
-  return (fabs(lpos[Acts::eLOC_X])
+  return (std::abs(lpos[Acts::eLOC_X])
           < m_valueStore.at(RectangleBounds::bv_halfX) + tol0);
 }
 
 inline bool
 RectangleBounds::insideLoc1(const Vector2D& lpos, double tol1) const
 {
-  return (fabs(lpos[Acts::eLOC_Y])
+  return (std::abs(lpos[Acts::eLOC_Y])
           < m_valueStore.at(RectangleBounds::bv_halfY) + tol1);
 }
 
diff --git a/Core/include/ACTS/Surfaces/TrapezoidBounds.hpp b/Core/include/ACTS/Surfaces/TrapezoidBounds.hpp
index d46b4670315597d71d0039fe8866a46ed3d13d45..de6fff1fdacf92b9b748c6c3621bd514de9cfef5 100644
--- a/Core/include/ACTS/Surfaces/TrapezoidBounds.hpp
+++ b/Core/include/ACTS/Surfaces/TrapezoidBounds.hpp
@@ -13,7 +13,7 @@
 #ifndef ACTS_SURFACES_TRAPEZOIDBOUNDS_H
 #define ACTS_SURFACES_TRAPEZOIDBOUNDS_H 1
 
-#include <math.h>
+#include <cmath>
 #include "ACTS/Surfaces/PlanarBounds.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/ParameterDefinitions.hpp"
@@ -289,7 +289,7 @@ TrapezoidBounds::inside(const Vector2D& lpos, const BoundaryCheck& bcheck) const
     return inside(lpos, bcheck.toleranceLoc0, bcheck.toleranceLoc1);
 
   // a fast FALSE
-  double fabsY   = fabs(lpos[Acts::eLOC_Y]);
+  double fabsY   = std::abs(lpos[Acts::eLOC_Y]);
   double max_ell = (*bcheck.lCovariance)(0, 0) > (*bcheck.lCovariance)(1, 1)
       ? (*bcheck.lCovariance)(0, 0)
       : (*bcheck.lCovariance)(1, 1);
@@ -297,7 +297,7 @@ TrapezoidBounds::inside(const Vector2D& lpos, const BoundaryCheck& bcheck) const
   if (fabsY > (m_valueStore.at(TrapezoidBounds::bv_halfY) + limit))
     return false;
   // a fast FALSE
-  double fabsX = fabs(lpos[Acts::eLOC_X]);
+  double fabsX = std::abs(lpos[Acts::eLOC_X]);
   if (fabsX > (m_valueStore.at(TrapezoidBounds::bv_maxHalfX) + limit))
     return false;
   // a fast TRUE
@@ -358,14 +358,14 @@ TrapezoidBounds::inside(const Vector2D& lpos, const BoundaryCheck& bcheck) const
 inline bool
 TrapezoidBounds::insideLoc0(const Vector2D& lpos, double tol0) const
 {
-  return (fabs(lpos[Acts::eLOC_X])
+  return (std::abs(lpos[Acts::eLOC_X])
           < m_valueStore.at(TrapezoidBounds::bv_maxHalfX) + tol0);
 }
 
 inline bool
 TrapezoidBounds::insideLoc1(const Vector2D& lpos, double tol1) const
 {
-  return (fabs(lpos[Acts::eLOC_Y])
+  return (std::abs(lpos[Acts::eLOC_Y])
           < m_valueStore.at(TrapezoidBounds::bv_halfY) + tol1);
 }
 
diff --git a/Core/include/ACTS/Surfaces/TriangleBounds.hpp b/Core/include/ACTS/Surfaces/TriangleBounds.hpp
index 74ea3f794d266adefc7508ccc2fffbc20ead1b03..b779b83d7ab7614c17009c98a02d2af09ac65b95 100644
--- a/Core/include/ACTS/Surfaces/TriangleBounds.hpp
+++ b/Core/include/ACTS/Surfaces/TriangleBounds.hpp
@@ -159,7 +159,7 @@ TriangleBounds::inside(const Vector2D& lpos, double tol0, double tol1) const
 
   // special case : lies on base ?
   double db = locB.first * locV.second - locB.second * locV.first;
-  if (fabs(db) < tol0) {
+  if (std::abs(db) < tol0) {
     double a = (locB.first != 0) ? -locV.first / locB.first
                                  : -locV.second / locB.second;
     if (a > -tol1 && a - 1. < tol1) return true;
@@ -168,7 +168,7 @@ TriangleBounds::inside(const Vector2D& lpos, double tol0, double tol1) const
 
   double dn = locB.first * locT.second - locB.second * locT.first;
 
-  if (fabs(dn) > fabs(tol0)) {
+  if (std::abs(dn) > std::abs(tol0)) {
     double t = (locB.first * locV.second - locB.second * locV.first) / dn;
     if (t > 0.) return false;
 
diff --git a/Core/include/ACTS/Utilities/BinningData.hpp b/Core/include/ACTS/Utilities/BinningData.hpp
index 79440b6848cd972343c226296b10eaf180031184..a284c592a4080ecb68e39b3323db964d8dbdc408 100644
--- a/Core/include/ACTS/Utilities/BinningData.hpp
+++ b/Core/include/ACTS/Utilities/BinningData.hpp
@@ -486,7 +486,7 @@ private:
         std::vector<float>::const_iterator mbvalue = m_boundaries.begin();
         for (; mbvalue != m_boundaries.end(); ++mbvalue) {
           // should define numerically stable
-          if (fabs((*mbvalue) - sBinMin) < 10e-10) {
+          if (std::abs((*mbvalue) - sBinMin) < 10e-10) {
             // copy the sub bin boundaries into the vector
             m_totalBoundaries.insert(m_totalBoundaries.begin(),
                                      subBinBoundaries.begin(),
diff --git a/Core/include/ACTS/Utilities/Helpers.hpp b/Core/include/ACTS/Utilities/Helpers.hpp
index 65a475623a813e71162fc697a5d4e0f6ab014196..f2352cf73d8abca774bc042e59f70d5d5de21c26 100644
--- a/Core/include/ACTS/Utilities/Helpers.hpp
+++ b/Core/include/ACTS/Utilities/Helpers.hpp
@@ -43,7 +43,7 @@ namespace Acts {
 inline double
 roundWithPrecision(double val, int precision)
 {
-  if (val < 0 && fabs(val) * std::pow(10, precision) < 1.) return -val;
+  if (val < 0 && std::abs(val) * std::pow(10, precision) < 1.) return -val;
   return val;
 }
 
diff --git a/Core/include/ACTS/Volumes/CuboidVolumeBounds.hpp b/Core/include/ACTS/Volumes/CuboidVolumeBounds.hpp
index f8191d33e2ecc67aa68e2a6ac3ecb53466801586..482d582d6c71fe6aa6b6eefdfcab7bcb9e82ec81 100644
--- a/Core/include/ACTS/Volumes/CuboidVolumeBounds.hpp
+++ b/Core/include/ACTS/Volumes/CuboidVolumeBounds.hpp
@@ -15,6 +15,7 @@
 
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Volumes/VolumeBounds.hpp"
+#include <cmath>
 
 namespace Acts {
 
@@ -149,9 +150,9 @@ CuboidVolumeBounds::clone() const
 inline bool
 CuboidVolumeBounds::inside(const Vector3D& pos, double tol) const
 {
-  return (fabs(pos.x()) <= m_valueStore.at(bv_halfX) + tol
-          && fabs(pos.y()) <= m_valueStore.at(bv_halfY) + tol
-          && fabs(pos.z()) <= m_valueStore.at(bv_halfZ) + tol);
+  return (std::abs(pos.x()) <= m_valueStore.at(bv_halfX) + tol
+          && std::abs(pos.y()) <= m_valueStore.at(bv_halfY) + tol
+          && std::abs(pos.z()) <= m_valueStore.at(bv_halfZ) + tol);
 }
 
 inline double
diff --git a/Core/include/ACTS/Volumes/CylinderVolumeBounds.hpp b/Core/include/ACTS/Volumes/CylinderVolumeBounds.hpp
index 3ff42b4ed0593a8cdc5547f3cb784969064e3da0..3fad79930df27623802bea6fd6b5f8bbe3c3bf64 100644
--- a/Core/include/ACTS/Volumes/CylinderVolumeBounds.hpp
+++ b/Core/include/ACTS/Volumes/CylinderVolumeBounds.hpp
@@ -15,6 +15,7 @@
 
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Volumes/VolumeBounds.hpp"
+#include <cmath>
 
 namespace Acts {
 
@@ -224,7 +225,7 @@ CylinderVolumeBounds::inside(const Vector3D& pos, double tol) const
                               && (ros <= m_valueStore[bv_outerRadius] + tol))
                            : false;
   bool insideZ
-      = insideR ? (fabs(pos.z()) <= m_valueStore[bv_halfZ] + tol) : false;
+      = insideR ? (std::abs(pos.z()) <= m_valueStore[bv_halfZ] + tol) : false;
   return (insideZ && insideR && insidePhi);
 }
 
diff --git a/Core/src/Digitization/RectangleSegmentation.cpp b/Core/src/Digitization/RectangleSegmentation.cpp
index 1fa1a532eda602c037e4307a49159ba9a5cce3f9..3bdd9d186098663845c12488bdbe0ef550aa563d 100644
--- a/Core/src/Digitization/RectangleSegmentation.cpp
+++ b/Core/src/Digitization/RectangleSegmentation.cpp
@@ -14,6 +14,7 @@
 #include "ACTS/Surfaces/PlaneSurface.hpp"
 #include "ACTS/Surfaces/RectangleBounds.hpp"
 #include "ACTS/Utilities/Helpers.hpp"
+#include <cmath>
 
 Acts::RectangularSegmentation::RectangularSegmentation(
     std::shared_ptr<const Acts::RectangleBounds> mBounds,
@@ -91,7 +92,7 @@ Acts::RectangularSegmentation::createSegmentationSurfaces(
   } else {
     // lorentz reduced Bounds
     double lorentzReducedHalfX
-        = m_activeBounds->halflengthX() - fabs(lorentzPlaneShiftX);
+        = m_activeBounds->halflengthX() - std::abs(lorentzPlaneShiftX);
     std::shared_ptr<const PlanarBounds> lorentzReducedBounds(
         new RectangleBounds(lorentzReducedHalfX,
                             m_activeBounds->halflengthY()));
@@ -118,7 +119,7 @@ Acts::RectangularSegmentation::createSegmentationSurfaces(
   std::shared_ptr<const PlanarBounds> xBinBounds(
       new RectangleBounds(m_activeBounds->halflengthY(), halfThickness));
   // now, let's create the shared bounds of all surfaces marking lorentz planes
-  double lorentzPlaneHalfX = fabs(halfThickness / cos(lorentzAngle));
+  double lorentzPlaneHalfX = std::abs(halfThickness / cos(lorentzAngle));
   // teh bounds of the lorentz plane
   std::shared_ptr<const PlanarBounds> lorentzPlaneBounds = (lorentzAngle == 0.)
       ? xBinBounds
diff --git a/Core/src/Digitization/TrapezoidSegmentation.cpp b/Core/src/Digitization/TrapezoidSegmentation.cpp
index 01ec753bc79e33a208d36ef360d8f448861f9e69..b257fcee58934ef6598e2e4cea469008e70beeb9 100644
--- a/Core/src/Digitization/TrapezoidSegmentation.cpp
+++ b/Core/src/Digitization/TrapezoidSegmentation.cpp
@@ -15,6 +15,7 @@
 #include "ACTS/Surfaces/RectangleBounds.hpp"
 #include "ACTS/Surfaces/TrapezoidBounds.hpp"
 #include "ACTS/Utilities/Helpers.hpp"
+#include <cmath>
 
 Acts::TrapezoidSegmentation::TrapezoidSegmentation(
     std::shared_ptr<const TrapezoidBounds> mBounds,
@@ -195,7 +196,7 @@ Acts::TrapezoidSegmentation::digitizationStep(const Vector3D& startStep,
   double lorentzDeltaX = -readoutDirection * stepCenter.z() * tan(lorentzAngle);
   // take the full drift length
   double driftInZ    = (halfThickness - readoutDirection * stepCenter.z());
-  double driftLength = fabs(driftInZ / cos(lorentzAngle));
+  double driftLength = std::abs(driftInZ / cos(lorentzAngle));
   // the projected center
   Vector2D stepCenterProjected(stepCenter.x() + lorentzDeltaX, stepCenter.y());
   // the cell & its center
diff --git a/Core/src/Extrapolation/RungeKuttaUtils.cpp b/Core/src/Extrapolation/RungeKuttaUtils.cpp
index b4772c05a522470b4d0d2487519b73beedd6e67f..9f0a8a9d0bd84bd7faea202039a8e427e99eb0b6 100644
--- a/Core/src/Extrapolation/RungeKuttaUtils.cpp
+++ b/Core/src/Extrapolation/RungeKuttaUtils.cpp
@@ -17,6 +17,7 @@
 #include "ACTS/Surfaces/PerigeeSurface.hpp"
 #include "ACTS/Surfaces/PlaneSurface.hpp"
 #include "ACTS/Surfaces/StrawSurface.hpp"
+#include <cmath>
 
 #ifndef __USE_GNU
 void
@@ -644,7 +645,7 @@ Acts::RungeKuttaUtils::stepEstimatorToCylinder(double*       S,
       Smax += Sq;
     }
   } else {
-    if (fabs(Smax) < .1) {
+    if (std::abs(Smax) < .1) {
       Q = false;
       return 0.;
     }
@@ -673,7 +674,7 @@ Acts::RungeKuttaUtils::stepEstimatorToCylinder(double*       S,
     return Smax;
   }
 
-  // if(fabs(Smin) < .001) {S[8]=-1.; return Smax;}
+  // if(std::abs(Smin) < .001) {S[8]=-1.; return Smax;}
 
   S[8] = 1.;
   return Smin;
@@ -1092,7 +1093,7 @@ Acts::RungeKuttaUtils::transformLocalToGlobal(bool                 useJac,
   P[4] = Sf * Se;  // Ay
   P[5] = Ce;       // Az
   P[6] = p[4];     // CM
-  if (fabs(P[6]) < .000000000000001) {
+  if (std::abs(P[6]) < .000000000000001) {
     P[6] < 0. ? P[6] = -.000000000000001 : P[6] = .000000000000001;
   }
 
diff --git a/Core/src/Surfaces/ConeBounds.cpp b/Core/src/Surfaces/ConeBounds.cpp
index 7e0cfa6d9304e18007773ecdce1a52d543a121de..76ad9825599f36ffa5e217e4258a951fe29f9145 100644
--- a/Core/src/Surfaces/ConeBounds.cpp
+++ b/Core/src/Surfaces/ConeBounds.cpp
@@ -13,7 +13,7 @@
 #include "ACTS/Surfaces/ConeBounds.hpp"
 #include <iomanip>
 #include <iostream>
-#include <math.h>
+#include <cmath>
 
 Acts::ConeBounds::ConeBounds(double alpha,
                              bool   symm,
@@ -84,7 +84,7 @@ Acts::ConeBounds::distanceToBoundary(const Acts::Vector2D& pos) const
   // find the minimum distance along the z direction
   double toMinZ = m_valueStore.at(ConeBounds::bv_minZ) - pos[Acts::eLOC_Z];
   double toMaxZ = pos[Acts::eLOC_Z] - m_valueStore.at(ConeBounds::bv_maxZ);
-  double toZ    = (fabs(toMinZ) < fabs(toMaxZ)) ? toMinZ : toMaxZ;
+  double toZ    = (std::abs(toMinZ) < std::abs(toMaxZ)) ? toMinZ : toMaxZ;
 
   // NB this works only if the localPos is in the same hemisphere as
   // the cone (i.e. if the localPos has z < 0 and the cone only
diff --git a/Core/src/Surfaces/ConeSurface.cpp b/Core/src/Surfaces/ConeSurface.cpp
index 0b4559dc0e5fdec878753e2b6a070a8cb0d94f6a..aad3f03c7ca6c43ed70cffd4cc46d3cb35fd7989 100644
--- a/Core/src/Surfaces/ConeSurface.cpp
+++ b/Core/src/Surfaces/ConeSurface.cpp
@@ -11,9 +11,10 @@
 ///////////////////////////////////////////////////////////////////
 
 #include "ACTS/Surfaces/ConeSurface.hpp"
-#include <assert.h>
+#include <cassert>
 #include <iomanip>
 #include <iostream>
+#include <cmath>
 #include "ACTS/Surfaces/RealQuadraticEquation.hpp"
 
 Acts::ConeSurface::ConeSurface(const ConeSurface& csf)
@@ -130,7 +131,7 @@ Acts::ConeSurface::globalToLocal(const Vector3D& gpos,
   // now decide on the quility of the transformation
   double inttol = r * 0.0001;
   inttol        = (inttol < 0.01) ? 0.01 : 0.01;  // ?
-  return ((fabs(loc3Dframe.perp() - r) > inttol) ? false : true);
+  return ((std::abs(loc3Dframe.perp() - r) > inttol) ? false : true);
 }
 
 Acts::Intersection
@@ -210,5 +211,5 @@ Acts::ConeSurface::pathCorrection(const Vector3D& gpos,
   if (m_transform) normalC = transform() * normalC;
   // back in global frame
   double cAlpha = normalC.dot(mom.unit());
-  return fabs(1. / cAlpha);
+  return std::abs(1. / cAlpha);
 }
diff --git a/Core/src/Surfaces/CylinderBounds.cpp b/Core/src/Surfaces/CylinderBounds.cpp
index 84101dfc21f405fdd35b482492a29e923df946ec..f7192b77ba2b75341b66a623819fa7c31bb578b0 100644
--- a/Core/src/Surfaces/CylinderBounds.cpp
+++ b/Core/src/Surfaces/CylinderBounds.cpp
@@ -13,22 +13,22 @@
 #include "ACTS/Surfaces/CylinderBounds.hpp"
 #include <iomanip>
 #include <iostream>
-#include <math.h>
+#include <cmath>
 
 Acts::CylinderBounds::CylinderBounds(double radius, double halez)
   : SurfaceBounds(CylinderBounds::bv_length), m_checkPhi(false)
 {
-  m_valueStore.at(CylinderBounds::bv_radius)        = fabs(radius);
+  m_valueStore.at(CylinderBounds::bv_radius)        = std::abs(radius);
   m_valueStore.at(CylinderBounds::bv_halfPhiSector) = M_PI;
-  m_valueStore.at(CylinderBounds::bv_halfZ)         = fabs(halez);
+  m_valueStore.at(CylinderBounds::bv_halfZ)         = std::abs(halez);
 }
 
 Acts::CylinderBounds::CylinderBounds(double radius, double haphi, double halez)
   : SurfaceBounds(CylinderBounds::bv_length), m_checkPhi(true)
 {
-  m_valueStore.at(CylinderBounds::bv_radius)        = fabs(radius);
+  m_valueStore.at(CylinderBounds::bv_radius)        = std::abs(radius);
   m_valueStore.at(CylinderBounds::bv_halfPhiSector) = haphi;
-  m_valueStore.at(CylinderBounds::bv_halfZ)         = fabs(halez);
+  m_valueStore.at(CylinderBounds::bv_halfZ)         = std::abs(halez);
 }
 
 Acts::CylinderBounds::CylinderBounds(double radius,
@@ -37,10 +37,10 @@ Acts::CylinderBounds::CylinderBounds(double radius,
                                      double halez)
   : SurfaceBounds(CylinderBounds::bv_length), m_checkPhi(true)
 {
-  m_valueStore.at(CylinderBounds::bv_radius)        = fabs(radius);
+  m_valueStore.at(CylinderBounds::bv_radius)        = std::abs(radius);
   m_valueStore.at(CylinderBounds::bv_averagePhi)    = averagephi;
   m_valueStore.at(CylinderBounds::bv_halfPhiSector) = haphi;
-  m_valueStore.at(CylinderBounds::bv_halfZ)         = fabs(halez);
+  m_valueStore.at(CylinderBounds::bv_halfZ)         = std::abs(halez);
 }
 
 Acts::CylinderBounds::~CylinderBounds()
@@ -63,11 +63,11 @@ Acts::CylinderBounds::distanceToBoundary(const Acts::Vector2D& pos) const
   const double pi2 = 2. * M_PI;
 
   double sZ
-      = fabs(pos[Acts::eLOC_Z]) - m_valueStore.at(CylinderBounds::bv_halfZ);
+      = std::abs(pos[Acts::eLOC_Z]) - m_valueStore.at(CylinderBounds::bv_halfZ);
   double wF = m_valueStore.at(CylinderBounds::bv_halfPhiSector);
   if (wF >= M_PI) return sZ;
   double dF
-      = fabs(pos[Acts::eLOC_RPHI] / m_valueStore.at(CylinderBounds::bv_radius)
+      = std::abs(pos[Acts::eLOC_RPHI] / m_valueStore.at(CylinderBounds::bv_radius)
              - m_valueStore.at(CylinderBounds::bv_averagePhi));
   if (dF > M_PI) dF = pi2 - dF;
   double sF
diff --git a/Core/src/Surfaces/CylinderSurface.cpp b/Core/src/Surfaces/CylinderSurface.cpp
index 3396dad6192cf7bf80d7888a8090499716acabc8..e8431f22ab9294d30f045baed491320ab034f5fb 100644
--- a/Core/src/Surfaces/CylinderSurface.cpp
+++ b/Core/src/Surfaces/CylinderSurface.cpp
@@ -11,9 +11,10 @@
 ///////////////////////////////////////////////////////////////////
 
 #include "ACTS/Surfaces/CylinderSurface.hpp"
-#include <assert.h>
+#include <cassert>
 #include <iomanip>
 #include <iostream>
+#include <cmath>
 #include "ACTS/Surfaces/RealQuadraticEquation.hpp"
 
 Acts::CylinderSurface::CylinderSurface(const CylinderSurface& csf)
@@ -150,7 +151,7 @@ Acts::CylinderSurface::globalToLocal(const Vector3D& gpos,
     radius = gpos.perp();
   }
   // return true or false
-  return ((fabs(radius - bounds().r()) > inttol) ? false : true);
+  return ((std::abs(radius - bounds().r()) > inttol) ? false : true);
 }
 
 bool
diff --git a/Core/src/Surfaces/DiamondBounds.cpp b/Core/src/Surfaces/DiamondBounds.cpp
index 925ae46413c2b0cb939f78a82a1025fa3b1e3975..6edf54e0f57fbeb60717a93ab3bb34587a16dbea 100644
--- a/Core/src/Surfaces/DiamondBounds.cpp
+++ b/Core/src/Surfaces/DiamondBounds.cpp
@@ -13,7 +13,7 @@
 #include "ACTS/Surfaces/DiamondBounds.hpp"
 #include <iomanip>
 #include <iostream>
-#include <math.h>
+#include <cmath>
 
 Acts::DiamondBounds::DiamondBounds(double minhalex,
                                    double medhalex,
@@ -98,11 +98,11 @@ Acts::DiamondBounds::insideFull(const Acts::Vector2D& locpo,
       > 2. * m_valueStore.at(DiamondBounds::bv_halfY2) + tol2)
     return false;
   // (2)
-  if (fabs(locpo[Acts::eLOC_X])
+  if (std::abs(locpo[Acts::eLOC_X])
       > (m_valueStore.at(DiamondBounds::bv_medHalfX) + tol1))
     return false;
   // (3)
-  if (fabs(locpo[Acts::eLOC_X])
+  if (std::abs(locpo[Acts::eLOC_X])
       < (fmin(m_valueStore.at(DiamondBounds::bv_minHalfX),
               m_valueStore.at(DiamondBounds::bv_maxHalfX))
          - tol1))
@@ -115,18 +115,18 @@ Acts::DiamondBounds::insideFull(const Acts::Vector2D& locpo,
            - m_valueStore.at(DiamondBounds::bv_minHalfX))
             / 2 / m_valueStore.at(DiamondBounds::bv_halfY1)
         : 0.;
-    return (fabs(locpo[Acts::eLOC_X])
+    return (std::abs(locpo[Acts::eLOC_X])
             <= m_valueStore.at(DiamondBounds::bv_medHalfX)
-                - k * fabs(locpo[Acts::eLOC_Y]));
+                - k * std::abs(locpo[Acts::eLOC_Y]));
   } else {
     double k = m_valueStore.at(DiamondBounds::bv_halfY2) > 0.
         ? (m_valueStore.at(DiamondBounds::bv_medHalfX)
            - m_valueStore.at(DiamondBounds::bv_maxHalfX))
             / 2 / m_valueStore.at(DiamondBounds::bv_halfY2)
         : 0.;
-    return (fabs(locpo[Acts::eLOC_X])
+    return (std::abs(locpo[Acts::eLOC_X])
             <= m_valueStore.at(DiamondBounds::bv_medHalfX)
-                - k * fabs(locpo[Acts::eLOC_Y]));
+                - k * std::abs(locpo[Acts::eLOC_Y]));
   }
 }
 
diff --git a/Core/src/Surfaces/DiscSurface.cpp b/Core/src/Surfaces/DiscSurface.cpp
index 7c7e82ff7e3fcaccb2be4a4b2ca523f5d64cb773..0a3a49b861829808c458d0a85b03a35d3d33f0e0 100644
--- a/Core/src/Surfaces/DiscSurface.cpp
+++ b/Core/src/Surfaces/DiscSurface.cpp
@@ -13,6 +13,7 @@
 #include "ACTS/Surfaces/DiscSurface.hpp"
 #include <iomanip>
 #include <iostream>
+#include <cmath>
 #include "ACTS/Surfaces/DiscTrapezoidalBounds.hpp"
 #include "ACTS/Surfaces/InfiniteBounds.hpp"
 #include "ACTS/Surfaces/RadialBounds.hpp"
@@ -104,7 +105,7 @@ Acts::DiscSurface::globalToLocal(const Acts::Vector3D& gpos,
   // transport it to the globalframe (very unlikely that this is not needed)
   Vector3D loc3Dframe = (transform().inverse()) * gpos;
   lpos                = Acts::Vector2D(loc3Dframe.perp(), loc3Dframe.phi());
-  return ((fabs(loc3Dframe.z()) > s_onSurfaceTolerance) ? false : true);
+  return ((std::abs(loc3Dframe.z()) > s_onSurfaceTolerance) ? false : true);
 }
 
 const Acts::Vector2D
@@ -149,7 +150,7 @@ Acts::DiscSurface::isOnSurface(const Vector3D&      glopo,
                                const BoundaryCheck& bcheck) const
 {
   Vector3D loc3Dframe = (transform().inverse()) * glopo;
-  if (fabs(loc3Dframe.z()) > (s_onSurfaceTolerance)) return false;
+  if (std::abs(loc3Dframe.z()) > (s_onSurfaceTolerance)) return false;
   return (
       bcheck
           ? bounds().inside(Vector2D(loc3Dframe.perp(), loc3Dframe.phi()), bcheck)
diff --git a/Core/src/Surfaces/DiscTrapezoidalBounds.cpp b/Core/src/Surfaces/DiscTrapezoidalBounds.cpp
index 35bfae91e86a70621548fdb276a2592ce8ffb455..25cd9c39cf53417c71724f299183e5a300e7ff9a 100644
--- a/Core/src/Surfaces/DiscTrapezoidalBounds.cpp
+++ b/Core/src/Surfaces/DiscTrapezoidalBounds.cpp
@@ -13,6 +13,7 @@
 #include "ACTS/Surfaces/DiscTrapezoidalBounds.hpp"
 #include <iomanip>
 #include <iostream>
+#include <cmath>
 
 Acts::DiscTrapezoidalBounds::DiscTrapezoidalBounds(double minhalfx,
                                                    double maxhalfx,
@@ -72,7 +73,7 @@ double
 Acts::DiscTrapezoidalBounds::distanceToBoundary(const Acts::Vector2D& pos) const
 {
   const double pi2        = 2. * M_PI;
-  double       alpha      = fabs(pos[Acts::eLOC_PHI]);
+  double       alpha      = std::abs(pos[Acts::eLOC_PHI]);
   if (alpha > M_PI) alpha = pi2 - alpha;
 
   double r = pos[Acts::eLOC_R];
@@ -96,7 +97,7 @@ Acts::DiscTrapezoidalBounds::distanceToBoundary(const Acts::Vector2D& pos) const
   // check if it is external in phi
   if ((alpha - m_valueStore.at(DiscTrapezoidalBounds::bv_halfPhiSector)) > 0.)
     return r
-        * fabs(sin(alpha
+        * std::abs(sin(alpha
                    - m_valueStore.at(DiscTrapezoidalBounds::bv_halfPhiSector)));
 
   // if here it is inside:
diff --git a/Core/src/Surfaces/EllipseBounds.cpp b/Core/src/Surfaces/EllipseBounds.cpp
index bfdf03e978b10f9e624986222c2edd2c8c469eee..79bba20469e6926b753294d27041143cc6dece7a 100644
--- a/Core/src/Surfaces/EllipseBounds.cpp
+++ b/Core/src/Surfaces/EllipseBounds.cpp
@@ -13,6 +13,7 @@
 #include "ACTS/Surfaces/EllipseBounds.hpp"
 #include <iomanip>
 #include <iostream>
+#include <cmath>
 
 Acts::EllipseBounds::EllipseBounds(double minradX,
                                    double minradY,
@@ -77,7 +78,7 @@ Acts::EllipseBounds::distanceToBoundary(const Vector2D& lpos) const
   if (m_valueStore.at(EllipseBounds::bv_halfPhiSector) < M_PI) {
     dF = atan2(cs, sn) - m_valueStore.at(EllipseBounds::bv_averagePhi);
     dF += (dF > M_PI) ? -pi2 : (dF < -M_PI) ? pi2 : 0;
-    double df = fabs(dF) - m_valueStore.at(EllipseBounds::bv_halfPhiSector);
+    double df = std::abs(dF) - m_valueStore.at(EllipseBounds::bv_halfPhiSector);
     sf        = r * sin(df);
     if (df > 0.) r *= cos(df);
   } else {
diff --git a/Core/src/Surfaces/LineSurface.cpp b/Core/src/Surfaces/LineSurface.cpp
index 9da59975e9053c0e8de573f5a7ee94001961f182..c421bf801230b008f433bf4a57d419b6c83405f3 100644
--- a/Core/src/Surfaces/LineSurface.cpp
+++ b/Core/src/Surfaces/LineSurface.cpp
@@ -13,6 +13,7 @@
 #include "ACTS/Surfaces/LineSurface.hpp"
 #include <iomanip>
 #include <iostream>
+#include <cmath>
 
 Acts::LineSurface::LineSurface(std::shared_ptr<Transform3D> htrans,
                                double                       radius,
@@ -141,7 +142,7 @@ Acts::LineSurface::intersectionEstimate(const Vector3D&      gpos,
   Vector3D mab(mb - ma);
   double   eaTeb = ea.dot(eb);
   double   denom = 1 - eaTeb * eaTeb;
-  if (fabs(denom) > 10e-7) {
+  if (std::abs(denom) > 10e-7) {
     double lambda0 = (mab.dot(ea) - mab.dot(eb) * eaTeb) / denom;
     // evaluate in terms of direction
     bool isValid = forceDir ? (lambda0 > 0.) : true;
diff --git a/Core/src/Surfaces/PlaneSurface.cpp b/Core/src/Surfaces/PlaneSurface.cpp
index 5adaef3b6dad28345c5f63d0068cf7f3e4937490..ddd2c2c295d020da1ec77e347bd95d8fefe7f008 100644
--- a/Core/src/Surfaces/PlaneSurface.cpp
+++ b/Core/src/Surfaces/PlaneSurface.cpp
@@ -13,6 +13,7 @@
 #include "ACTS/Surfaces/PlaneSurface.hpp"
 #include <iomanip>
 #include <iostream>
+#include <cmath>
 #include "ACTS/Surfaces/InfiniteBounds.hpp"
 #include "ACTS/Surfaces/RectangleBounds.hpp"
 #include "ACTS/Utilities/Identifier.hpp"
@@ -37,7 +38,7 @@ Acts::PlaneSurface::PlaneSurface(const Vector3D& center, const Vector3D& normal)
   /// U = Z x T if T not parallel to Z otherwise U = X x T
   /// V = T x U
   Vector3D T = normal.normalized();
-  Vector3D U = fabs(T.dot(Vector3D::UnitZ())) < 0.99
+  Vector3D U = std::abs(T.dot(Vector3D::UnitZ())) < 0.99
       ? Vector3D::UnitZ().cross(T).normalized()
       : Vector3D::UnitX().cross(T).normalized();
   Vector3D         V = T.cross(U);
@@ -111,7 +112,7 @@ Acts::PlaneSurface::isOnSurface(const Vector3D&      glopo,
 {
   /// the chance that there is no transform is almost 0, let's apply it
   Vector3D loc3Dframe = (transform().inverse()) * glopo;
-  if (fabs(loc3Dframe.z()) > s_onSurfaceTolerance) return false;
+  if (std::abs(loc3Dframe.z()) > s_onSurfaceTolerance) return false;
   return (bcheck ? bounds().inside(Vector2D(loc3Dframe.x(), loc3Dframe.y()), bcheck)
                : true);
 }
diff --git a/Core/src/Surfaces/RadialBounds.cpp b/Core/src/Surfaces/RadialBounds.cpp
index 2c2a741d81498fedde185dd9778d595c231bd3b0..b208fdbaae1fdfc7029159e8b71f2bcc79e97693 100644
--- a/Core/src/Surfaces/RadialBounds.cpp
+++ b/Core/src/Surfaces/RadialBounds.cpp
@@ -13,6 +13,7 @@
 #include "ACTS/Surfaces/RadialBounds.hpp"
 #include <iomanip>
 #include <iostream>
+#include <cmath>
 
 Acts::RadialBounds::RadialBounds(double minrad, double maxrad, double hphisec)
   : Acts::DiscBounds(RadialBounds::bv_length)
@@ -65,7 +66,7 @@ Acts::RadialBounds::distanceToBoundary(const Acts::Vector2D& lpos) const
   double dF = 0.;
 
   if (m_valueStore.at(RadialBounds::bv_halfPhiSector) < M_PI) {
-    dF = fabs(lpos[Acts::eLOC_PHI]
+    dF = std::abs(lpos[Acts::eLOC_PHI]
               - m_valueStore.at(RadialBounds::bv_averagePhi));
     if (dF > M_PI) dF = pi2 - dF;
     dF -= m_valueStore.at(RadialBounds::bv_halfPhiSector);
diff --git a/Core/src/Surfaces/RectangleBounds.cpp b/Core/src/Surfaces/RectangleBounds.cpp
index 95b0e533b924913562dac404b65015bf6d5157ba..1b382a3c9307ce25a26b556db2ba7d0daab4d432 100644
--- a/Core/src/Surfaces/RectangleBounds.cpp
+++ b/Core/src/Surfaces/RectangleBounds.cpp
@@ -13,6 +13,7 @@
 #include "ACTS/Surfaces/RectangleBounds.hpp"
 #include <iomanip>
 #include <iostream>
+#include <cmath>
 
 Acts::RectangleBounds::RectangleBounds(double halex, double haley)
   : Acts::PlanarBounds(RectangleBounds::bv_length)
@@ -35,8 +36,8 @@ Acts::RectangleBounds::operator=(const RectangleBounds& recbo)
 double
 Acts::RectangleBounds::distanceToBoundary(const Acts::Vector2D& lpos) const
 {
-  double dx = fabs(lpos[0]) - m_valueStore.at(RectangleBounds::bv_halfX);
-  double dy = fabs(lpos[1]) - m_valueStore.at(RectangleBounds::bv_halfY);
+  double dx = std::abs(lpos[0]) - m_valueStore.at(RectangleBounds::bv_halfX);
+  double dy = std::abs(lpos[1]) - m_valueStore.at(RectangleBounds::bv_halfY);
 
   if (dx <= 0. || dy <= 0.) {
     if (dx > dy)
@@ -44,7 +45,7 @@ Acts::RectangleBounds::distanceToBoundary(const Acts::Vector2D& lpos) const
     else
       return dy;
   }
-  return sqrt(dx * dx + dy * dy);
+  return std::sqrt(dx * dx + dy * dy);
 }
 
 // ostream operator overload
diff --git a/Core/src/Surfaces/TrapezoidBounds.cpp b/Core/src/Surfaces/TrapezoidBounds.cpp
index 60239164e088d1d29686c574777a5827b159c7d9..234a710e0ae70d352fb3dbf7c9741539307cc1dd 100644
--- a/Core/src/Surfaces/TrapezoidBounds.cpp
+++ b/Core/src/Surfaces/TrapezoidBounds.cpp
@@ -13,16 +13,16 @@
 #include "ACTS/Surfaces/TrapezoidBounds.hpp"
 #include <iomanip>
 #include <iostream>
-#include <math.h>
+#include <cmath>
 
 Acts::TrapezoidBounds::TrapezoidBounds(double minhalex,
                                        double maxhalex,
                                        double haley)
   : Acts::PlanarBounds(TrapezoidBounds::bv_length), m_alpha(0.), m_beta(0.)
 {
-  m_valueStore.at(TrapezoidBounds::bv_minHalfX) = fabs(minhalex);
-  m_valueStore.at(TrapezoidBounds::bv_maxHalfX) = fabs(maxhalex);
-  m_valueStore.at(TrapezoidBounds::bv_halfY)    = fabs(haley);
+  m_valueStore.at(TrapezoidBounds::bv_minHalfX) = std::abs(minhalex);
+  m_valueStore.at(TrapezoidBounds::bv_maxHalfX) = std::abs(maxhalex);
+  m_valueStore.at(TrapezoidBounds::bv_halfY)    = std::abs(haley);
   if (m_valueStore.at(TrapezoidBounds::bv_minHalfX)
       > m_valueStore.at(TrapezoidBounds::bv_maxHalfX))
     std::swap(m_valueStore.at(TrapezoidBounds::bv_minHalfX),
@@ -37,10 +37,10 @@ Acts::TrapezoidBounds::TrapezoidBounds(double minhalex,
 {
   double gamma = (alpha > beta) ? (alpha - 0.5 * M_PI) : (beta - 0.5 * M_PI);
   // now fill them
-  m_valueStore.at(TrapezoidBounds::bv_minHalfX) = fabs(minhalex);
+  m_valueStore.at(TrapezoidBounds::bv_minHalfX) = std::abs(minhalex);
   m_valueStore.at(TrapezoidBounds::bv_maxHalfX) = minhalex
       + (2. * m_valueStore.at(TrapezoidBounds::bv_halfY)) * tan(gamma);
-  m_valueStore.at(TrapezoidBounds::bv_halfY) = fabs(haley);
+  m_valueStore.at(TrapezoidBounds::bv_halfY) = std::abs(haley);
 }
 
 Acts::TrapezoidBounds::~TrapezoidBounds()
@@ -74,8 +74,8 @@ Acts::TrapezoidBounds::insideFull(const Acts::Vector2D& lpos,
 {
   // the cases:
   // the cases:
-  double fabsX = fabs(lpos[Acts::eLOC_X]);
-  double fabsY = fabs(lpos[Acts::eLOC_Y]);
+  double fabsX = std::abs(lpos[Acts::eLOC_X]);
+  double fabsY = std::abs(lpos[Acts::eLOC_Y]);
   // (1) a fast FALSE
   if (fabsY > (m_valueStore.at(TrapezoidBounds::bv_halfY) + tol1)) return false;
   // (2) a fast FALSE
@@ -93,7 +93,7 @@ Acts::TrapezoidBounds::insideFull(const Acts::Vector2D& lpos,
       / (m_valueStore.at(TrapezoidBounds::bv_maxHalfX)
          - m_valueStore.at(TrapezoidBounds::bv_minHalfX))
       * ((lpos[Acts::eLOC_X] > 0.) ? 1.0 : -1.0);
-  double d = -fabs(k) * 0.5 * (m_valueStore.at(TrapezoidBounds::bv_maxHalfX)
+  double d = -std::abs(k) * 0.5 * (m_valueStore.at(TrapezoidBounds::bv_maxHalfX)
                                + m_valueStore.at(TrapezoidBounds::bv_minHalfX));
   return (isAbove(lpos, tol0, tol1, k, d));
 }
diff --git a/Core/src/Tools/CylinderVolumeHelper.cpp b/Core/src/Tools/CylinderVolumeHelper.cpp
index b5ec105610305a884b44f341559456a2ad3781a4..22126ea27b67f827734389a357e97bdde66fd167 100644
--- a/Core/src/Tools/CylinderVolumeHelper.cpp
+++ b/Core/src/Tools/CylinderVolumeHelper.cpp
@@ -26,6 +26,7 @@
 #include "ACTS/Volumes/AbstractVolume.hpp"
 #include "ACTS/Volumes/BoundarySurfaceT.hpp"
 #include "ACTS/Volumes/CylinderVolumeBounds.hpp"
+#include <cmath>
 
 Acts::CylinderVolumeHelper::CylinderVolumeHelper(
     const Acts::CylinderVolumeHelper::Config& cvhConfig,
@@ -185,7 +186,7 @@ Acts::CylinderVolumeHelper::createTrackingVolume(
   // create a Transform3D and VolumeBounds out of the zMin/zMax
   double halflengthZ = 0.5 * (zMax - zMin);
   double zPosition   = 0.5 * (zMin + zMax);
-  zPosition          = fabs(zPosition) < 0.1 ? 0. : zPosition;
+  zPosition          = std::abs(zPosition) < 0.1 ? 0. : zPosition;
 
   // now create the cylinder volume bounds
   cBounds = rMin > 0.1 ? new CylinderVolumeBounds(rMin, rMax, halflengthZ)
@@ -275,7 +276,7 @@ Acts::CylinderVolumeHelper::createGapTrackingVolume(
       // create the layer
       layers.push_back(createCylinderLayer(0.5 * (zMinLayer + zMaxLayer),
                                            (*layerPropIter),
-                                           fabs(0.5 * (zMaxLayer - zMinLayer)),
+                                           std::abs(0.5 * (zMaxLayer - zMinLayer)),
                                            m_cfg.passiveLayerThickness,
                                            m_cfg.passiveLayerPhiBins,
                                            m_cfg.passiveLayerRzBins));
@@ -356,7 +357,7 @@ Acts::CylinderVolumeHelper::createContainerTrackingVolume(
   }
   // check whether it is a r-binned case or a z-binned case
   bool rCase
-      = fabs(firstVolumeBounds->innerRadius() - lastVolumeBounds->innerRadius())
+      = std::abs(firstVolumeBounds->innerRadius() - lastVolumeBounds->innerRadius())
       > 0.1;
 
   // fill these ones depending on the rCase though assignment - no parsing at
@@ -388,12 +389,12 @@ Acts::CylinderVolumeHelper::createContainerTrackingVolume(
   double zPos = 0.5 * (zMin + zMax);
   // create the HEP transform from the stuff known so far
   std::shared_ptr<Transform3D> topVolumeTransform
-      = fabs(zPos) > 0.1 ? std::make_shared<Transform3D>() : nullptr;
+      = std::abs(zPos) > 0.1 ? std::make_shared<Transform3D>() : nullptr;
   if (topVolumeTransform) (*topVolumeTransform) = Translation3D(0., 0., zPos);
   // create the bounds from the information gathered so far
-  CylinderVolumeBounds* topVolumeBounds = fabs(rMin) > 0.1
-      ? new CylinderVolumeBounds(rMin, rMax, 0.5 * fabs(zMax - zMin))
-      : new CylinderVolumeBounds(rMax, 0.5 * fabs(zMax - zMin));
+  CylinderVolumeBounds* topVolumeBounds = std::abs(rMin) > 0.1
+      ? new CylinderVolumeBounds(rMin, rMax, 0.5 * std::abs(zMax - zMin))
+      : new CylinderVolumeBounds(rMax, 0.5 * std::abs(zMax - zMin));
 
   // some screen output
   ACTS_VERBOSE("Container voume bounds are " << (*topVolumeBounds));
@@ -521,7 +522,7 @@ Acts::CylinderVolumeHelper::estimateAndCheckDimension(
                "layers + envelope covers");
   // the z from the layers w and w/o envelopes
   double zEstFromLayerEnv    = 0.5 * ((layerZmax) + (layerZmin));
-  double halflengthFromLayer = 0.5 * fabs((layerZmax) - (layerZmin));
+  double halflengthFromLayer = 0.5 * std::abs((layerZmax) - (layerZmin));
 
   bool concentric = (zEstFromLayerEnv * zEstFromLayerEnv < 0.001);
 
@@ -538,7 +539,7 @@ Acts::CylinderVolumeHelper::estimateAndCheckDimension(
     (*transform) = Translation3D(0., 0., zEstFromLayerEnv);
   } else if (transform && !cylinderVolumeBounds) {
     // create the CylinderBounds from parsed layer inputs
-    double halflengthFromLayer = 0.5 * fabs((layerZmax) - (layerZmin));
+    double halflengthFromLayer = 0.5 * std::abs((layerZmax) - (layerZmin));
     cylinderVolumeBounds
         = new CylinderVolumeBounds(layerRmin, layerRmax, halflengthFromLayer);
   }
@@ -864,7 +865,7 @@ Acts::CylinderVolumeHelper::glueTrackingVolumes(
 
     // the transform of the new boundary surface
     std::shared_ptr<Transform3D> transform = nullptr;
-    if (fabs(zMin + zMax) > 0.1) {
+    if (std::abs(zMin + zMax) > 0.1) {
       // it's not a concentric cylinder, so create a transform
       auto pTransform = std::make_shared<Transform3D>();
       (*pTransform)   = Translation3D(Vector3D(0., 0., 0.5 * (zMin + zMax)));
@@ -950,7 +951,7 @@ Acts::CylinderVolumeHelper::createCylinderLayer(double z,
                                                        << r);
   // positioning
   std::shared_ptr<Transform3D> transform = 0;
-  transform = (fabs(z) > 0.1) ? std::make_shared<Transform3D>() : 0;
+  transform = (std::abs(z) > 0.1) ? std::make_shared<Transform3D>() : 0;
   if (transform) (*transform) = Translation3D(0., 0., z);
 
   // z-binning
@@ -998,7 +999,7 @@ Acts::CylinderVolumeHelper::createDiscLayer(double z,
 
   // positioning
   std::shared_ptr<Transform3D> transform
-      = fabs(z) > 0.1 ? std::make_shared<Transform3D>() : 0;
+      = std::abs(z) > 0.1 ? std::make_shared<Transform3D>() : 0;
   if (transform) (*transform) = Translation3D(0., 0., z);
 
   // R is the primary binning for the material
diff --git a/Core/src/Tools/LayerArrayCreator.cpp b/Core/src/Tools/LayerArrayCreator.cpp
index 917fc9b074a5a64ad44ed65b8c2b72fedba41e9b..3751b9eddba8a284973b6e0fce5edda46a0f00d0 100644
--- a/Core/src/Tools/LayerArrayCreator.cpp
+++ b/Core/src/Tools/LayerArrayCreator.cpp
@@ -25,6 +25,7 @@
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/GeometryObjectSorter.hpp"
 #include "ACTS/Utilities/GeometryStatics.hpp"
+#include <cmath>
 
 std::unique_ptr<const Acts::LayerArray>
 Acts::LayerArrayCreator::layerArray(const LayerVector& layersInput,
@@ -93,7 +94,7 @@ Acts::LayerArrayCreator::layerArray(const LayerVector& layersInput,
       if (navigationValue != (layerValue - 0.5 * layerThickness)) {
         // create the navigation layer surface from the layer
         std::unique_ptr<const Surface> navLayerSurface(createNavigationSurface(
-            *layIter, bValue, -fabs(layerValue - navigationValue)));
+            *layIter, bValue, -std::abs(layerValue - navigationValue)));
         ACTS_VERBOSE("arbitrary : creating a  NavigationLayer at "
                      << (navLayerSurface->binningPosition(bValue)).x()
                      << ", "
diff --git a/Core/src/Tools/LayerCreator.cpp b/Core/src/Tools/LayerCreator.cpp
index 6d564af665cf6ba5f293c969376ee6bee7425483..60c16af25a97cf03c8dee91e9a0e4e603a174e84 100644
--- a/Core/src/Tools/LayerCreator.cpp
+++ b/Core/src/Tools/LayerCreator.cpp
@@ -18,6 +18,7 @@
 #include "ACTS/Surfaces/RadialBounds.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/Units.hpp"
+#include <cmath>
 
 Acts::LayerCreator::LayerCreator(const Acts::LayerCreator::Config& lcConfig,
                                  std::unique_ptr<Logger>           logger)
@@ -241,7 +242,7 @@ Acts::LayerCreator::discLayer(const std::vector<const Surface*>&  surfaces,
 {
   // layer parametres
   double layerZ         = 0.5 * (layerZmin + layerZmax);
-  double layerThickness = fabs(layerZmax - layerZmin);
+  double layerThickness = std::abs(layerZmax - layerZmin);
 
   // adjust the layer radius
   ACTS_VERBOSE("Creating a disk Layer:");
@@ -352,7 +353,7 @@ Acts::LayerCreator::radialDistance(const Vector3D& pos1,
   Vector3D mab(mb - ma);
   double   eaTeb = ea.dot(eb);
   double   denom = 1 - eaTeb * eaTeb;
-  if (fabs(denom) > 10e-7) {
+  if (std::abs(denom) > 10e-7) {
     double lambda0 = (mab.dot(ea) - mab.dot(eb) * eaTeb) / denom;
     // evaluate validaty in terms of bounds
     if (lambda0 < 1. && lambda0 > 0.) return (ma + lambda0 * ea).perp();
diff --git a/Core/src/Tools/SurfaceArrayCreator.cpp b/Core/src/Tools/SurfaceArrayCreator.cpp
index 8837af157dae4478c496bcf527dca685f70ac74c..0938f8df8cc33ebbb6877eefcdbcba53b460bab9 100644
--- a/Core/src/Tools/SurfaceArrayCreator.cpp
+++ b/Core/src/Tools/SurfaceArrayCreator.cpp
@@ -17,6 +17,7 @@
 #include "ACTS/Utilities/BinnedArrayXD.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/Units.hpp"
+#include <cmath>
 
 std::unique_ptr<Acts::SurfaceArray>
 Acts::SurfaceArrayCreator::surfaceArrayOnCylinder(
@@ -315,7 +316,7 @@ Acts::SurfaceArrayCreator::createArbitraryBinUtility(
                      end(surf),
                      back_inserter(keys),
                      [](const Acts::Surface* a, const Acts::Surface* b) {
-                       return (fabs(a->center().phi() - b->center().phi())
+                       return (std::abs(a->center().phi() - b->center().phi())
                                < 10e-12);
                      });
     // the phi-center position of the previous surface
@@ -398,7 +399,7 @@ Acts::SurfaceArrayCreator::createArbitraryBinUtility(
                      end(surf),
                      back_inserter(keys),
                      [](const Acts::Surface* a, const Acts::Surface* b) {
-                       return (fabs(a->center().z() - b->center().z())
+                       return (std::abs(a->center().z() - b->center().z())
                                < Acts::units::_um);
                      });
     // the z-center position of the previous surface
@@ -469,7 +470,7 @@ Acts::SurfaceArrayCreator::createArbitraryBinUtility(
                      end(surf),
                      back_inserter(keys),
                      [](const Acts::Surface* a, const Acts::Surface* b) {
-                       return (fabs(a->center().perp() - b->center().perp())
+                       return (std::abs(a->center().perp() - b->center().perp())
                                < Acts::units::_um);
                      });
     // the r-center position of the previous surface
@@ -566,13 +567,13 @@ Acts::SurfaceArrayCreator::createEquidistantBinUtility(
                      end(surf),
                      back_inserter(keys),
                      [](const Acts::Surface* a, const Acts::Surface* b) {
-                       return (fabs(a->center().phi() - b->center().phi())
+                       return (std::abs(a->center().phi() - b->center().phi())
                                < 10e-12);
                      });
     // set minimum and maximum
     double min  = keys.front()->center().phi();
     double max  = keys.back()->center().phi();
-    double step = fabs(max - min) / (keys.size() - 1);
+    double step = std::abs(max - min) / (keys.size() - 1);
     minimum     = min - 0.5 * step;
     maximum     = max + 0.5 * step;
     // phi correction
@@ -598,13 +599,13 @@ Acts::SurfaceArrayCreator::createEquidistantBinUtility(
                      end(surf),
                      back_inserter(keys),
                      [](const Acts::Surface* a, const Acts::Surface* b) {
-                       return (fabs(a->center().z() - b->center().z())
+                       return (std::abs(a->center().z() - b->center().z())
                                < Acts::units::_um);
                      });
     // set minimum and maximum
     double min  = keys.front()->center().z();
     double max  = keys.back()->center().z();
-    double step = fabs(max - min) / (keys.size() - 1);
+    double step = std::abs(max - min) / (keys.size() - 1);
     minimum     = min - 0.5 * step;
     maximum     = max + 0.5 * step;
   } else {
@@ -620,13 +621,13 @@ Acts::SurfaceArrayCreator::createEquidistantBinUtility(
                      end(surf),
                      back_inserter(keys),
                      [](const Acts::Surface* a, const Acts::Surface* b) {
-                       return (fabs(a->center().perp() - b->center().perp())
+                       return (std::abs(a->center().perp() - b->center().perp())
                                < Acts::units::_um);
                      });
     // set the minimum and maximum
     double min  = keys.front()->center().perp();
     double max  = keys.back()->center().perp();
-    double step = fabs(max - min) / (keys.size() - 1);
+    double step = std::abs(max - min) / (keys.size() - 1);
     minimum     = min - 0.5 * step;
     maximum     = max + 0.5 * step;
   }
diff --git a/Core/src/Volumes/CuboidVolumeBounds.cpp b/Core/src/Volumes/CuboidVolumeBounds.cpp
index 39cd7fbebb385c975de6d63eb45d338d6243094a..a0580f155a939eb64f240b2877c63f2a723e8a56 100644
--- a/Core/src/Volumes/CuboidVolumeBounds.cpp
+++ b/Core/src/Volumes/CuboidVolumeBounds.cpp
@@ -12,7 +12,7 @@
 
 #include "ACTS/Volumes/CuboidVolumeBounds.hpp"
 #include <iostream>
-#include <math.h>
+#include <cmath>
 #include "ACTS/Surfaces/PlaneSurface.hpp"
 #include "ACTS/Surfaces/RectangleBounds.hpp"
 #include "ACTS/Surfaces/Surface.hpp"
diff --git a/Core/src/Volumes/CylinderVolumeBounds.cpp b/Core/src/Volumes/CylinderVolumeBounds.cpp
index 7ac3daf3d6d9832226dca6f60236602201b8210b..583e966251bcb45fe0de1190c6b8a712a824e54f 100644
--- a/Core/src/Volumes/CylinderVolumeBounds.cpp
+++ b/Core/src/Volumes/CylinderVolumeBounds.cpp
@@ -12,7 +12,7 @@
 
 #include "ACTS/Volumes/CylinderVolumeBounds.hpp"
 #include <iostream>
-#include <math.h>
+#include <cmath>
 #include "ACTS/Surfaces/CylinderBounds.hpp"
 #include "ACTS/Surfaces/CylinderSurface.hpp"
 #include "ACTS/Surfaces/DiscSurface.hpp"
@@ -31,9 +31,9 @@ Acts::CylinderVolumeBounds::CylinderVolumeBounds(double radius, double halez)
   : VolumeBounds(), m_valueStore(4, 0.)
 {
   m_valueStore.at(bv_innerRadius)   = 0.;
-  m_valueStore.at(bv_outerRadius)   = fabs(radius);
+  m_valueStore.at(bv_outerRadius)   = std::abs(radius);
   m_valueStore.at(bv_halfPhiSector) = M_PI;
-  m_valueStore.at(bv_halfZ)         = fabs(halez);
+  m_valueStore.at(bv_halfZ)         = std::abs(halez);
 }
 
 Acts::CylinderVolumeBounds::CylinderVolumeBounds(double rinner,
@@ -41,10 +41,10 @@ Acts::CylinderVolumeBounds::CylinderVolumeBounds(double rinner,
                                                  double halez)
   : VolumeBounds(), m_valueStore(4, 0.)
 {
-  m_valueStore.at(bv_innerRadius)   = fabs(rinner);
-  m_valueStore.at(bv_outerRadius)   = fabs(router);
+  m_valueStore.at(bv_innerRadius)   = std::abs(rinner);
+  m_valueStore.at(bv_outerRadius)   = std::abs(router);
   m_valueStore.at(bv_halfPhiSector) = M_PI;
-  m_valueStore.at(bv_halfZ)         = fabs(halez);
+  m_valueStore.at(bv_halfZ)         = std::abs(halez);
 }
 
 Acts::CylinderVolumeBounds::CylinderVolumeBounds(double rinner,
@@ -53,10 +53,10 @@ Acts::CylinderVolumeBounds::CylinderVolumeBounds(double rinner,
                                                  double halez)
   : VolumeBounds(), m_valueStore(4, 0.)
 {
-  m_valueStore.at(bv_innerRadius)   = fabs(rinner);
-  m_valueStore.at(bv_outerRadius)   = fabs(router);
-  m_valueStore.at(bv_halfPhiSector) = fabs(haphi);
-  m_valueStore.at(bv_halfZ)         = fabs(halez);
+  m_valueStore.at(bv_innerRadius)   = std::abs(rinner);
+  m_valueStore.at(bv_outerRadius)   = std::abs(router);
+  m_valueStore.at(bv_halfPhiSector) = std::abs(haphi);
+  m_valueStore.at(bv_halfZ)         = std::abs(halez);
 }
 
 Acts::CylinderVolumeBounds::CylinderVolumeBounds(
@@ -111,7 +111,7 @@ Acts::CylinderVolumeBounds::decomposeToSurfaces(
         new CylinderSurface(transformPtr, innerCylinderBounds()));
 
   // the cylinder is sectoral
-  if (fabs(halfPhiSector() - M_PI) > s_numericalStable) {
+  if (std::abs(halfPhiSector() - M_PI) > s_numericalStable) {
     std::shared_ptr<const PlanarBounds> sp12Bounds = sectorPlaneBounds();
     // sectorPlane 1 (negative phi)
     Transform3D* sp1Transform = new Transform3D(
diff --git a/Core/src/Volumes/DoubleTrapezoidVolumeBounds.cpp b/Core/src/Volumes/DoubleTrapezoidVolumeBounds.cpp
index c6886827712d8b028530028fa21bd6bc65a2b9ec..fd2913ca1b2a904ada19c2367ef5a9a2cf3353e5 100644
--- a/Core/src/Volumes/DoubleTrapezoidVolumeBounds.cpp
+++ b/Core/src/Volumes/DoubleTrapezoidVolumeBounds.cpp
@@ -13,7 +13,7 @@
 #include "ACTS/Volumes/DoubleTrapezoidVolumeBounds.hpp"
 #include <iomanip>
 #include <iostream>
-#include <math.h>
+#include <cmath>
 #include "ACTS/Surfaces/DiamondBounds.hpp"
 #include "ACTS/Surfaces/PlaneSurface.hpp"
 #include "ACTS/Surfaces/RectangleBounds.hpp"
@@ -245,7 +245,7 @@ Acts::DoubleTrapezoidVolumeBounds::faceZXRectangleBoundsTop() const
 bool
 Acts::DoubleTrapezoidVolumeBounds::inside(const Vector3D& pos, double tol) const
 {
-  if (fabs(pos.z()) > m_valueStore.at(bv_halfZ) + tol) return false;
+  if (std::abs(pos.z()) > m_valueStore.at(bv_halfZ) + tol) return false;
   if (pos.y() < -2 * m_valueStore.at(bv_halfY1) - tol) return false;
   if (pos.y() > 2 * m_valueStore.at(bv_halfY2) - tol) return false;
   DiamondBounds* faceXYBounds = faceXYDiamondBounds();
diff --git a/Core/src/Volumes/TrapezoidVolumeBounds.cpp b/Core/src/Volumes/TrapezoidVolumeBounds.cpp
index fb0a0afd10c2d85262f9ca88f76c26930e8fb56f..9f9adc4313800c6d430021715db5f0ff1c87c225 100644
--- a/Core/src/Volumes/TrapezoidVolumeBounds.cpp
+++ b/Core/src/Volumes/TrapezoidVolumeBounds.cpp
@@ -13,7 +13,7 @@
 #include "ACTS/Volumes/TrapezoidVolumeBounds.hpp"
 #include <iomanip>
 #include <iostream>
-#include <math.h>
+#include <cmath>
 #include "ACTS/Surfaces/PlaneSurface.hpp"
 #include "ACTS/Surfaces/RectangleBounds.hpp"
 #include "ACTS/Surfaces/TrapezoidBounds.hpp"
@@ -215,8 +215,8 @@ Acts::TrapezoidVolumeBounds::faceZXRectangleBoundsTop() const
 bool
 Acts::TrapezoidVolumeBounds::inside(const Vector3D& pos, double tol) const
 {
-  if (fabs(pos.z()) > m_valueStore.at(bv_halfZ) + tol) return false;
-  if (fabs(pos.y()) > m_valueStore.at(bv_halfY) + tol) return false;
+  if (std::abs(pos.z()) > m_valueStore.at(bv_halfZ) + tol) return false;
+  if (std::abs(pos.y()) > m_valueStore.at(bv_halfY) + tol) return false;
   TrapezoidBounds* faceXYBounds = faceXYTrapezoidBounds();
   Vector2D         locp(pos.x(), pos.y());
   bool inside(faceXYBounds->inside(locp, BoundaryCheck(true, true, tol, tol)));
diff --git a/Examples/bin/main.cpp b/Examples/bin/main.cpp
index 962d061533bc67e06b6097073319b394f1c255f7..cbe05e743b0c15fe22c09820b07a28e390a7e0aa 100644
--- a/Examples/bin/main.cpp
+++ b/Examples/bin/main.cpp
@@ -7,7 +7,7 @@
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #include <iostream>
-#include <math.h>
+#include <cmath>
 #include <memory>
 #include <random>
 #include "ACTS/Detector/TrackingGeometry.hpp"
diff --git a/Examples/src/BuildGenericDetector.cpp b/Examples/src/BuildGenericDetector.cpp
index d44b394b9aef595bd3017f2ed82ad0d81f66e981..bf24ae46dfa33ab8734684a022bd475fd317b4ee 100644
--- a/Examples/src/BuildGenericDetector.cpp
+++ b/Examples/src/BuildGenericDetector.cpp
@@ -10,6 +10,7 @@
 #include <array>
 #include <iostream>
 #include <vector>
+#include <cmath>
 #include "ACTS/Detector/TrackingGeometry.hpp"
 #include "ACTS/Examples/GenericLayerBuilder.hpp"
 #include "ACTS/Material/Material.hpp"
@@ -84,7 +85,7 @@ modulePositionsCylinder(double radius,
   double phiStep = 2 * M_PI / (nPhiBins);
   double minPhi  = -M_PI + 0.5 * phiStep;
   double zStart  = -0.5 * (nZbins - 1) * (2 * moduleHalfLength - lOverlap);
-  double zStep   = 2 * fabs(zStart) / (nZbins - 1);
+  double zStep   = 2 * std::abs(zStart) / (nZbins - 1);
   // loop over the bins
   for (size_t zBin = 0; zBin < size_t(nZbins); ++zBin) {
     // prepare z and r
diff --git a/Tests/EventData/ParameterSetTests.cpp b/Tests/EventData/ParameterSetTests.cpp
index ba06b1072838d9869996eb8e4d8ea41d27fe3d79..837bdfba9f7aab629e76391cbda579922a8cd10f 100644
--- a/Tests/EventData/ParameterSetTests.cpp
+++ b/Tests/EventData/ParameterSetTests.cpp
@@ -210,36 +210,36 @@ namespace Test {
         residual               = parSet_1.residual(parSet_2);
 
         // local parameters are unbound -> check for usual difference
-        if (fabs(residual(0) - delta_loc1) > tol) {
+        if (std::abs(residual(0) - delta_loc1) > tol) {
           BOOST_CHECK(false);
           break;
         }
-        if (fabs(residual(1) - delta_loc2) > tol) {
+        if (std::abs(residual(1) - delta_loc2) > tol) {
           BOOST_CHECK(false);
           break;
         }
         // phi is a cyclic parameter -> check that (unsigned) difference is not
         // larger than half period
         // check that corrected(corrected(phi_2) + residual) == corrected(phi_1)
-        if (fabs(get_cyclic_value(get_cyclic_value(phi_2, phi_min, phi_max)
+        if (std::abs(get_cyclic_value(get_cyclic_value(phi_2, phi_min, phi_max)
                                       + residual(2),
                                   phi_min,
                                   phi_max)
                  - get_cyclic_value(phi_1, phi_min, phi_max))
                 > tol
-            or fabs(residual(2)) > (phi_max - phi_min) / 2) {
+            or std::abs(residual(2)) > (phi_max - phi_min) / 2) {
           BOOST_CHECK(false);
           break;
         }
         // theta is bound -> check that (unsigned) difference is not larger then
         // allowed range, check corrected difference
-        if (fabs(residual(3) - delta_theta) > tol
-            or fabs(residual(3)) > (theta_max - theta_min)) {
+        if (std::abs(residual(3) - delta_theta) > tol
+            or std::abs(residual(3)) > (theta_max - theta_min)) {
           BOOST_CHECK(false);
           break;
         }
         // qop is unbound -> check usual difference
-        if (fabs(residual(4) - delta_qop) > tol) {
+        if (std::abs(residual(4) - delta_qop) > tol) {
           BOOST_CHECK(false);
           break;
         }
@@ -573,21 +573,21 @@ namespace Test {
         = (cyclic.getParameter<ParID_t::ePHI>() - small_number) / (max - min);
     BOOST_CHECK(cyclic.getParameter<ParID_t::ePHI>() >= min);
     BOOST_CHECK(cyclic.getParameter<ParID_t::ePHI>() < max);
-    BOOST_CHECK(fabs(multiple - std::floor(multiple + 0.5)) < tol);
+    BOOST_CHECK(std::abs(multiple - std::floor(multiple + 0.5)) < tol);
 
     cyclic.setParameter<ParID_t::ePHI>(large_number);
     multiple
         = (cyclic.getParameter<ParID_t::ePHI>() - large_number) / (max - min);
     BOOST_CHECK(cyclic.getParameter<ParID_t::ePHI>() >= min);
     BOOST_CHECK(cyclic.getParameter<ParID_t::ePHI>() < max);
-    BOOST_CHECK(fabs(multiple - std::floor(multiple + 0.5)) < tol);
+    BOOST_CHECK(std::abs(multiple - std::floor(multiple + 0.5)) < tol);
 
     cyclic.setParameter<ParID_t::ePHI>(normal_number);
     multiple
         = (cyclic.getParameter<ParID_t::ePHI>() - normal_number) / (max - min);
     BOOST_CHECK(cyclic.getParameter<ParID_t::ePHI>() >= min);
     BOOST_CHECK(cyclic.getParameter<ParID_t::ePHI>() < max);
-    BOOST_CHECK(fabs(multiple - std::floor(multiple + 0.5)) < tol);
+    BOOST_CHECK(std::abs(multiple - std::floor(multiple + 0.5)) < tol);
 
     // check residual calculation
 
diff --git a/Tests/Propagator/propagation.cpp b/Tests/Propagator/propagation.cpp
index 27ef8d98618d79ef42279574980df157c6ef4749..2c4ddbb229ef1805f502469ca3ee1a53d01aa284 100644
--- a/Tests/Propagator/propagation.cpp
+++ b/Tests/Propagator/propagation.cpp
@@ -81,7 +81,7 @@ namespace Test {
     // clang-format on
 
     // calculate bending radius
-    double r = fabs(Nat2SI<units::MOMENTUM>(pT) / (q * Bz));
+    double r = std::abs(Nat2SI<units::MOMENTUM>(pT) / (q * Bz));
     // calculate number of turns of helix
     double turns = options.max_path_length / (2 * M_PI * r) * sin(theta);
     // respect direction of curl
@@ -93,7 +93,7 @@ namespace Test {
     if (exp_phi > M_PI) exp_phi -= 2 * M_PI;
 
     // calculate expected position
-    double exp_z = z + pz / pT * 2 * M_PI * r * fabs(turns);
+    double exp_z = z + pz / pT * 2 * M_PI * r * std::abs(turns);
 
     // calculate center of bending circle in transverse plane
     double xc, yc;
diff --git a/Tests/Utilities/BinningData_tests.cpp b/Tests/Utilities/BinningData_tests.cpp
index e62994a3a8c6b7141f1a8a5fbedf11d25b8f3fa1..869d3910c1ae664bbc26a01290e1625588366af5 100644
--- a/Tests/Utilities/BinningData_tests.cpp
+++ b/Tests/Utilities/BinningData_tests.cpp
@@ -11,6 +11,7 @@
 #include "ACTS/Utilities/BinningData.hpp"
 #include <boost/test/included/unit_test.hpp>
 #include "ACTS/Utilities/BinningType.hpp"
+#include <cmath>
 
 namespace Acts {
 namespace Test {
@@ -139,7 +140,7 @@ namespace Test {
     BOOST_CHECK_EQUAL(phiData_arb_binary.bins(), size_t(nBins_binary));
 
     // h/etaData
-    BOOST_TEST((fabs(etaData_eq.value(eta0Position) - 0.) < 1e-5));
+    BOOST_TEST((std::abs(etaData_eq.value(eta0Position) - 0.) < 1e-5));
   }
 
   // test bin values
@@ -402,7 +403,7 @@ namespace Test {
                                         float(-M_PI + 4.5 * phiStep)};
 
     for (size_t ib = 0; ib < phiCenters_eq.size(); ++ib)
-      BOOST_TEST((fabs(phiData_eq.center(ib) - phiCenters_eq[ib]) < 1e-3));
+      BOOST_TEST((std::abs(phiData_eq.center(ib) - phiCenters_eq[ib]) < 1e-3));
   }
 
   // special test for phi binning