From 98bcd1514c08c0e76c62decb9c40cb673664a622 Mon Sep 17 00:00:00 2001
From: Matthias Kleiner <matthias.kleiner@cern.ch>
Date: Tue, 20 Oct 2026 19:15:45 +0200
Subject: [PATCH] TrackParametrization: use double math functions for TrackParD
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Follow-up to #15986. TrackParametrization, TrackParametrizationWithError and
TrackUtils call the float-only gpu::CAMath functions (Sqrt, ASin, ATan2, ...),
so TrackParD / PropagatorD is effectively float. Compared with a std:: double
reference, propagateParamTo(x, b[3]) from x = 85 to 250 cm in 100 steps is off
by 0.4 to 74 µm, depending on pt and the field (toy fields with a transverse
component, no material), and finite-difference derivatives through it are
noisy. This affects the barrel alignment (Detectors/Align) and the ITS3
alignment; no reconstruction workflow uses the double types.

Call math_utils::detail::{sqrt, asin, ...}<value_t> instead, as Propagator.cxx
already does: for double these dispatch to the double functions, for float to
the same functions as before. Add GPUCommonMath::Expd and
math_utils::detail::exp for the remaining Exp calls. This restores the
math_utils::detail calls of cd9835aa55 (switched back to CAMath in bf65d3c754)
and extends them to TrackParametrizationWithError and TrackUtils.

Float results are bit-identical; TrackParD now agrees exactly with the double
reference.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 .../include/MathUtils/detail/basicMath.h      | 11 +++
 .../TrackParametrization.h                    | 20 ++--
 .../TrackParametrizationWithError.h           |  2 +-
 .../ReconstructionDataFormats/TrackUtils.h    | 12 +--
 .../src/TrackParametrization.cxx              | 84 ++++++++--------
 .../src/TrackParametrizationWithError.cxx     | 98 +++++++++----------
 GPU/Common/GPUCommonMath.h                    |  2 +
 7 files changed, 121 insertions(+), 108 deletions(-)

diff --git a/Common/MathUtils/include/MathUtils/detail/basicMath.h b/Common/MathUtils/include/MathUtils/detail/basicMath.h
index cd9ea213dd510..c5d8a91288136 100644
--- a/Common/MathUtils/include/MathUtils/detail/basicMath.h
+++ b/Common/MathUtils/include/MathUtils/detail/basicMath.h
@@ -91,6 +91,12 @@ GPUdi() T log(T x)
   return o2::gpu::GPUCommonMath::Log(x);
 };
 
+template <class T>
+GPUdi() T exp(T x)
+{
+  return o2::gpu::GPUCommonMath::Exp(x);
+};
+
 #if !defined(GPUCA_GPUCODE_DEVICE) || defined(__CUDACC__) || defined(__HIPCC__) // Double specialisations, not for OpenCL and Metal
 template <>
 GPUhdi() double copysign(double x, double y)
@@ -107,6 +113,11 @@ GPUdi() double log(double x)
 {
   return o2::gpu::GPUCommonMath::Logd(x);
 };
+template <>
+GPUdi() double exp(double x)
+{
+  return o2::gpu::GPUCommonMath::Expd(x);
+};
 #endif
 
 #ifndef GPUCA_GPUCODE_DEVICE // Double overloads using std namespace not visible on GPU
diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h
index ee9b3c10e05b7..83680bbe0332e 100644
--- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h
+++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrization.h
@@ -402,7 +402,7 @@ GPUdi() auto TrackParametrization<value_T>::getR2() const -> value_t
 template <typename value_T>
 GPUdi() auto TrackParametrization<value_T>::getR() const -> value_t
 {
-  return gpu::CAMath::Sqrt(getR2());
+  return math_utils::detail::sqrt<value_t>(getR2());
 }
 
 //____________________________________________________________
@@ -441,7 +441,7 @@ GPUdi() auto TrackParametrization<value_T>::getCsp2() const -> value_t
 template <typename value_T>
 GPUdi() auto TrackParametrization<value_T>::getCsp() const -> value_t
 {
-  return gpu::CAMath::Sqrt(getCsp2());
+  return math_utils::detail::sqrt<value_t>(getCsp2());
 }
 
 //____________________________________________________________
@@ -518,7 +518,7 @@ GPUdi() void TrackParametrization<value_T>::getCircleParamsLoc(value_t bz, o2::m
   constexpr value_t MinSagitta = 0.01f, TPCMidR = 160.f, MinCurv = 8 * MinSagitta / (TPCMidR * TPCMidR);
   if (gpu::CAMath::Abs(c.rC) > MinCurv) {
     c.rC = 1.f / getCurvature(bz);
-    value_t sn = getSnp(), cs = gpu::CAMath::Sqrt((1.f - sn) * (1.f + sn));
+    value_t sn = getSnp(), cs = math_utils::detail::sqrt<value_t>((1.f - sn) * (1.f + sn));
     c.xC = getX() - sn * c.rC; // center in tracking
     c.yC = getY() + cs * c.rC; // frame. Note: r is signed!!!
     c.rC = gpu::CAMath::Abs(c.rC);
@@ -546,7 +546,7 @@ GPUdi() void TrackParametrization<value_T>::getLineParams(o2::math_utils::Interv
   // get line parameterization as { x = x0 + xSlp*t, y = y0 + ySlp*t }
   o2::math_utils::detail::sincos(getAlpha(), sna, csa);
   o2::math_utils::detail::rotateZ<value_t>(getX(), getY(), ln.getX0(), ln.getY0(), sna, csa); // reference point in global frame
-  value_t snp = getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp));
+  value_t snp = getSnp(), csp = math_utils::detail::sqrt<value_t>((1.f - snp) * (1.f + snp));
   ln.setDX(csp * csa - snp * sna);
   ln.setDY(snp * csa + csp * sna);
 }
@@ -577,7 +577,7 @@ template <typename value_T>
 GPUdi() auto TrackParametrization<value_T>::getPhi() const -> value_t
 {
   // track pt direction phi (in 0:2pi range)
-  value_t phi = gpu::CAMath::ASin(getSnp()) + getAlpha();
+  value_t phi = math_utils::detail::asin<value_t>(getSnp()) + getAlpha();
   math_utils::detail::bringTo02Pi<value_t>(phi);
   return phi;
 }
@@ -587,7 +587,7 @@ template <typename value_T>
 GPUdi() auto TrackParametrization<value_T>::getPhiPos() const -> value_t
 {
   // angle of track position (in -pi:pi range)
-  value_t phi = gpu::CAMath::ATan2(getY(), getX()) + getAlpha();
+  value_t phi = math_utils::detail::atan2<value_t>(getY(), getX()) + getAlpha();
   math_utils::detail::bringTo02Pi<value_t>(phi);
   return phi;
 }
@@ -638,7 +638,7 @@ template <typename value_T>
 GPUdi() auto TrackParametrization<value_T>::getPInv() const -> value_t
 {
   // return the inverted track momentum
-  return getPtInv() / gpu::CAMath::Sqrt(1.f + getTgl() * getTgl()); // getPtInv() is protected against being 0, full charge accounted
+  return getPtInv() / math_utils::detail::sqrt<value_t>(1.f + getTgl() * getTgl()); // getPtInv() is protected against being 0, full charge accounted
 }
 
 //____________________________________________________________
@@ -662,7 +662,7 @@ template <typename value_T>
 GPUdi() auto TrackParametrization<value_T>::getE() const -> value_t
 {
   // return the track energy
-  return gpu::CAMath::Sqrt(getE2());
+  return math_utils::detail::sqrt<value_t>(getE2());
 }
 
 //____________________________________________________________
@@ -677,14 +677,14 @@ GPUdi() auto TrackParametrization<value_T>::getPt() const -> value_t
 template <typename value_T>
 GPUdi() auto TrackParametrization<value_T>::getTheta() const -> value_t
 {
-  return constants::math::PIHalf - gpu::CAMath::ATan(mP[3]);
+  return constants::math::PIHalf - math_utils::detail::atan<value_t>(mP[3]);
 }
 
 //____________________________________________________________
 template <typename value_T>
 GPUdi() auto TrackParametrization<value_T>::getEta() const -> value_t
 {
-  return -gpu::CAMath::Log(gpu::CAMath::Tan(0.5f * getTheta()));
+  return -math_utils::detail::log<value_t>(math_utils::detail::tan<value_t>(0.5f * getTheta()));
 }
 
 //_______________________________________________________
diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h
index 81280d090be71..435a6d8daeba4 100644
--- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h
+++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackParametrizationWithError.h
@@ -449,7 +449,7 @@ GPUdi() void TrackParametrizationWithError<value_T>::updateCovCorr(const value_t
     oldDiag[i] = mC[diagI];
     mC[diagI] += delta2[i];
     for (int j = 0; j < i; j++) {
-      mC[CovarMap[i][j]] *= gpu::CAMath::Sqrt(mC[diagI] * mC[DiagMap[j]] / (oldDiag[i] * oldDiag[j]));
+      mC[CovarMap[i][j]] *= math_utils::detail::sqrt<value_t>(mC[diagI] * mC[DiagMap[j]] / (oldDiag[i] * oldDiag[j]));
     }
   }
 #pragma GCC diagnostic pop
diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackUtils.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackUtils.h
index 6befcd0dfc898..b17cc2c7ab405 100644
--- a/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackUtils.h
+++ b/DataFormats/Reconstruction/include/ReconstructionDataFormats/TrackUtils.h
@@ -76,10 +76,10 @@ GPUd() void g3helx3(value_T qfield, value_T step, std::array<value_T, 7>& vect)
 
   value_T tsint, sintt, sint, cos1t;
   if (gpu::CAMath::Abs(tet) > value_T(0.03)) {
-    sint = gpu::CAMath::Sin(tet);
+    sint = math_utils::detail::sin<value_T>(tet);
     sintt = sint / tet;
     tsint = (tet - sint) / tet;
-    value_T t = gpu::CAMath::Sin(value_T(0.5) * tet);
+    value_T t = math_utils::detail::sin<value_T>(value_T(0.5) * tet);
     cos1t = value_T(2) * t * t / tet;
   } else {
     tsint = tet * tet / value_T(6);
@@ -133,15 +133,15 @@ GPUd() value_T BetheBlochSolid(value_T bg, value_T rho, value_T kp1, value_T kp2
 
   //*** Density effect
   value_T d2 = value_T(0);
-  const value_T x = gpu::CAMath::Log(bg);
-  const value_T lhwI = gpu::CAMath::Log(value_T(28.816e-9) * gpu::CAMath::Sqrt(rho * meanZA) / meanI);
+  const value_T x = math_utils::detail::log<value_T>(bg);
+  const value_T lhwI = math_utils::detail::log<value_T>(value_T(28.816e-9) * math_utils::detail::sqrt<value_T>(rho * meanZA) / meanI);
   if (x > kp2) {
     d2 = lhwI + x - value_T(0.5);
   } else if (x > kp1) {
     double r = (kp2 - x) / (kp2 - kp1);
     d2 = lhwI + x - value_T(0.5) + (value_T(0.5) - lhwI - kp1) * r * r * r;
   }
-  auto dedx = mK * meanZA / beta2 * (value_T(0.5) * gpu::CAMath::Log(value_T(2) * me * bg2 * maxT / (meanI * meanI)) - beta2 - d2);
+  auto dedx = mK * meanZA / beta2 * (value_T(0.5) * math_utils::detail::log<value_T>(value_T(2) * me * bg2 * maxT / (meanI * meanI)) - beta2 - d2);
   return dedx > value_T(0) ? dedx : value_T(0);
 }
 
@@ -177,7 +177,7 @@ GPUd() value_T BetheBlochSolidOpt(value_T bg)
 
   //*** Density effect
   value_T d2 = value_T(0);
-  const value_T x = gpu::CAMath::Log(bg);
+  const value_T x = math_utils::detail::log<value_T>(bg);
   if (x > kp2) {
     d2 = lhwI - value_T(0.5) + x;
   } else if (x > kp1) {
diff --git a/DataFormats/Reconstruction/src/TrackParametrization.cxx b/DataFormats/Reconstruction/src/TrackParametrization.cxx
index fb398bbcf07cf..beafe79cb54ce 100644
--- a/DataFormats/Reconstruction/src/TrackParametrization.cxx
+++ b/DataFormats/Reconstruction/src/TrackParametrization.cxx
@@ -49,9 +49,9 @@ GPUd() TrackParametrization<value_T>::TrackParametrization(const dim3_t& xyz, co
   value_t radPos2 = xyz[0] * xyz[0] + xyz[1] * xyz[1];
   value_t alp = 0;
   if (sectorAlpha || radPos2 < 1) {
-    alp = gpu::CAMath::ATan2(pxpypz[1], pxpypz[0]);
+    alp = math_utils::detail::atan2<value_t>(pxpypz[1], pxpypz[0]);
   } else {
-    alp = gpu::CAMath::ATan2(xyz[1], xyz[0]);
+    alp = math_utils::detail::atan2<value_t>(xyz[1], xyz[0]);
   }
   if (sectorAlpha) {
     alp = math_utils::detail::angle2Alpha<value_t>(alp);
@@ -62,7 +62,7 @@ GPUd() TrackParametrization<value_T>::TrackParametrization(const dim3_t& xyz, co
   // protection against cosp<0
   if (cs * pxpypz[0] + sn * pxpypz[1] < 0) {
     LOG(debug) << "alpha from phiPos() will invalidate this track parameters, overriding to alpha from phi()";
-    alp = gpu::CAMath::ATan2(pxpypz[1], pxpypz[0]);
+    alp = math_utils::detail::atan2<value_t>(pxpypz[1], pxpypz[0]);
     if (sectorAlpha) {
       alp = math_utils::detail::angle2Alpha<value_t>(alp);
     }
@@ -93,7 +93,7 @@ GPUd() TrackParametrization<value_T>::TrackParametrization(const dim3_t& xyz, co
   math_utils::detail::rotateZ<value_t>(ver, -alp);
   math_utils::detail::rotateZ<value_t>(mom, -alp);
   //
-  value_t ptI = 1.f / gpu::CAMath::Sqrt(mom[0] * mom[0] + mom[1] * mom[1]);
+  value_t ptI = 1.f / math_utils::detail::sqrt<value_t>(mom[0] * mom[0] + mom[1] * mom[1]);
   mX = ver[0];
   mAlpha = alp;
   mP[kY] = ver[1];
@@ -121,7 +121,7 @@ GPUd() bool TrackParametrization<value_T>::getPxPyPzGlo(dim3_t& pxyz) const
     return false;
   }
   value_t cs, sn, pt = getPt();
-  value_t r = gpu::CAMath::Sqrt((1.f - getSnp()) * (1.f + getSnp()));
+  value_t r = math_utils::detail::sqrt<value_t>((1.f - getSnp()) * (1.f + getSnp()));
   math_utils::detail::sincos(getAlpha(), sn, cs);
   pxyz[0] = pt * (r * cs - getSnp() * sn);
   pxyz[1] = pt * (getSnp() * cs + r * sn);
@@ -140,8 +140,8 @@ GPUd() bool TrackParametrization<value_T>::getPosDirGlo(std::array<value_t, 9>&
     return false;
   }
   value_t &sn = posdirp[7], &cs = posdirp[8];
-  value_t csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp));
-  value_t cstht = gpu::CAMath::Sqrt(1.f + getTgl() * getTgl());
+  value_t csp = math_utils::detail::sqrt<value_t>((1.f - snp) * (1.f + snp));
+  value_t cstht = math_utils::detail::sqrt<value_t>(1.f + getTgl() * getTgl());
   value_t csthti = 1.f / cstht;
   math_utils::detail::sincos(getAlpha(), sn, cs);
   posdirp[0] = getX() * cs - getY() * sn;
@@ -168,7 +168,7 @@ GPUd() bool TrackParametrization<value_T>::rotateParam(value_t alpha)
   //
   value_t ca = 0, sa = 0;
   math_utils::detail::sincos(alpha - getAlpha(), sa, ca);
-  value_t snp = getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp)); // Improve precision
+  value_t snp = getSnp(), csp = math_utils::detail::sqrt<value_t>((1.f - snp) * (1.f + snp)); // Improve precision
   // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
   // direction in local frame is along the X axis
   if ((csp * ca + snp * sa) < 0) {
@@ -201,7 +201,7 @@ GPUd() bool TrackParametrization<value_T>::rotateParam(value_t& alpha, value_t&
   //
   math_utils::detail::bringToPMPi<value_t>(alpha);
   math_utils::detail::sincos(alpha - getAlpha(), sa, ca);
-  value_t snp = getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp)); // Improve precision
+  value_t snp = getSnp(), csp = math_utils::detail::sqrt<value_t>((1.f - snp) * (1.f + snp)); // Improve precision
   // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle direction in local frame is along the X axis
   if ((csp * ca + snp * sa) < 0) {
     // LOGF(warning,"Rotation failed: local cos(phi) would become {:.2f}", csp * ca + snp * sa);
@@ -250,18 +250,18 @@ GPUd() bool TrackParametrization<value_T>::propagateParamTo(value_t xk, const di
   if (gpu::CAMath::Abs(f1) > constants::math::Almost1 || gpu::CAMath::Abs(f2) > constants::math::Almost1) {
     return false;
   }
-  value_t r1 = gpu::CAMath::Sqrt((1.f - f1) * (1.f + f1));
+  value_t r1 = math_utils::detail::sqrt<value_t>((1.f - f1) * (1.f + f1));
   if (gpu::CAMath::Abs(r1) < constants::math::Almost0) {
     return false;
   }
-  value_t r2 = gpu::CAMath::Sqrt((1.f - f2) * (1.f + f2));
+  value_t r2 = math_utils::detail::sqrt<value_t>((1.f - f2) * (1.f + f2));
   if (gpu::CAMath::Abs(r2) < constants::math::Almost0) {
     return false;
   }
   value_t dy2dx = (f1 + f2) / (r1 + r2);
-  value_t step = (gpu::CAMath::Abs(x2r) < 0.05f) ? dx * gpu::CAMath::Abs(r2 + f2 * dy2dx)                                              // chord
-                                                 : 2.f * CAMath::ASin(0.5f * dx * gpu::CAMath::Sqrt(1.f + dy2dx * dy2dx) * crv) / crv; // arc
-  step *= gpu::CAMath::Sqrt(1.f + getTgl() * getTgl());
+  value_t step = (gpu::CAMath::Abs(x2r) < 0.05f) ? dx * gpu::CAMath::Abs(r2 + f2 * dy2dx)                                                                                   // chord
+                                                 : 2.f * math_utils::detail::asin<value_t>(0.5f * dx * math_utils::detail::sqrt<value_t>(1.f + dy2dx * dy2dx) * crv) / crv; // arc
+  step *= math_utils::detail::sqrt<value_t>(1.f + getTgl() * getTgl());
   //
   // get the track x,y,z,px/p,py/p,pz/p,p,sinAlpha,cosAlpha in the Global System
   std::array<value_t, 9> vecLab{0.f};
@@ -271,13 +271,13 @@ GPUd() bool TrackParametrization<value_T>::propagateParamTo(value_t xk, const di
 
   // rotate to the system where Bx=By=0.
   value_t bxy2 = b[0] * b[0] + b[1] * b[1];
-  value_t bt = gpu::CAMath::Sqrt(bxy2);
+  value_t bt = math_utils::detail::sqrt<value_t>(bxy2);
   value_t cosphi = 1.f, sinphi = 0.f;
   if (bt > constants::math::Almost0) {
     cosphi = b[0] / bt;
     sinphi = b[1] / bt;
   }
-  value_t bb = gpu::CAMath::Sqrt(bxy2 + b[2] * b[2]);
+  value_t bb = math_utils::detail::sqrt<value_t>(bxy2 + b[2] * b[2]);
   value_t costet = 1.f, sintet = 0.f;
   if (bb > constants::math::Almost0) {
     costet = b[2] / bb;
@@ -326,7 +326,7 @@ GPUd() bool TrackParametrization<value_T>::propagateParamTo(value_t xk, const di
   }
 
   // Calculate the track parameters
-  t = 1.f / gpu::CAMath::Sqrt(vecLab[3] * vecLab[3] + vecLab[4] * vecLab[4]);
+  t = 1.f / math_utils::detail::sqrt<value_t>(vecLab[3] * vecLab[3] + vecLab[4] * vecLab[4]);
   mX = xk;
   mP[kY] = y;
   mP[kZ] = z;
@@ -356,11 +356,11 @@ GPUd() bool TrackParametrization<value_T>::propagateParamTo(value_t xk, value_t
   if ((gpu::CAMath::Abs(f1) > constants::math::Almost1) || (gpu::CAMath::Abs(f2) > constants::math::Almost1)) {
     return false;
   }
-  value_t r1 = gpu::CAMath::Sqrt((1.f - f1) * (1.f + f1));
+  value_t r1 = math_utils::detail::sqrt<value_t>((1.f - f1) * (1.f + f1));
   if (gpu::CAMath::Abs(r1) < constants::math::Almost0) {
     return false;
   }
-  value_t r2 = gpu::CAMath::Sqrt((1.f - f2) * (1.f + f2));
+  value_t r2 = math_utils::detail::sqrt<value_t>((1.f - f2) * (1.f + f2));
   if (gpu::CAMath::Abs(r2) < constants::math::Almost0) {
     return false;
   }
@@ -379,7 +379,7 @@ GPUd() bool TrackParametrization<value_T>::propagateParamTo(value_t xk, value_t
     if (gpu::CAMath::Abs(arg) > constants::math::Almost1) {
       return false;
     }
-    value_t rot = CAMath::ASin(arg);                // more economic version from Yura.
+    value_t rot = math_utils::detail::asin<value_t>(arg); // more economic version from Yura.
     if (f1 * f1 + f2 * f2 > 1.f && f1 * f2 < 0.f) { // special cases of large rotations or large abs angles
       if (f2 > 0.f) {
         rot = constants::math::PI - rot; //
@@ -404,7 +404,7 @@ GPUd() bool TrackParametrization<value_T>::propagateParamToDCA(const math_utils:
   // propagate track to DCA to the vertex
   value_t sn, cs, alp = getAlpha();
   math_utils::detail::sincos(alp, sn, cs);
-  value_t x = getX(), y = getY(), snp = getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp));
+  value_t x = getX(), y = getY(), snp = getSnp(), csp = math_utils::detail::sqrt<value_t>((1.f - snp) * (1.f + snp));
   value_t xv = vtx.X() * cs + vtx.Y() * sn, yv = -vtx.X() * sn + vtx.Y() * cs, zv = vtx.Z();
   x -= xv;
   y -= yv;
@@ -419,8 +419,8 @@ GPUd() bool TrackParametrization<value_T>::propagateParamToDCA(const math_utils:
   }
   value_t crv = getCurvature(b);
   value_t tgfv = -(crv * x - snp) / (crv * y + csp);
-  sn = tgfv / gpu::CAMath::Sqrt(1.f + tgfv * tgfv);
-  cs = gpu::CAMath::Sqrt((1.f - sn) * (1.f + sn));
+  sn = tgfv / math_utils::detail::sqrt<value_t>(1.f + tgfv * tgfv);
+  cs = math_utils::detail::sqrt<value_t>((1.f - sn) * (1.f + sn));
   cs = (gpu::CAMath::Abs(tgfv) > constants::math::Almost0) ? sn / tgfv : constants::math::Almost1;
 
   x = xv * cs + yv * sn;
@@ -428,7 +428,7 @@ GPUd() bool TrackParametrization<value_T>::propagateParamToDCA(const math_utils:
   xv = x;
 
   auto tmpT(*this); // operate on the copy to recover after the failure
-  alp += gpu::CAMath::ASin(sn);
+  alp += math_utils::detail::asin<value_t>(sn);
   if (!tmpT.rotateParam(alp) || !tmpT.propagateParamTo(xv, b)) {
 #ifndef GPUCA_ALIGPUCODE
     LOG(debug) << "failed to propagate to alpha=" << alp << " X=" << xv << " for vertex "
@@ -469,11 +469,11 @@ GPUd() bool TrackParametrization<value_T>::getYZAt(value_t xk, value_t b, value_
   if ((gpu::CAMath::Abs(f1) > constants::math::Almost1) || (gpu::CAMath::Abs(f2) > constants::math::Almost1)) {
     return false;
   }
-  value_t r1 = gpu::CAMath::Sqrt((1.f - f1) * (1.f + f1));
+  value_t r1 = math_utils::detail::sqrt<value_t>((1.f - f1) * (1.f + f1));
   if (gpu::CAMath::Abs(r1) < constants::math::Almost0) {
     return false;
   }
-  value_t r2 = gpu::CAMath::Sqrt((1.f - f2) * (1.f + f2));
+  value_t r2 = math_utils::detail::sqrt<value_t>((1.f - f2) * (1.f + f2));
   if (gpu::CAMath::Abs(r2) < constants::math::Almost0) {
     return false;
   }
@@ -490,7 +490,7 @@ GPUd() bool TrackParametrization<value_T>::getYZAt(value_t xk, value_t b, value_
     //    double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
     //    track1 += rot/crv*track3;
     //
-    value_t rot = CAMath::ASin(r1 * f2 - r2 * f1);  // more economic version from Yura.
+    value_t rot = math_utils::detail::asin<value_t>(r1 * f2 - r2 * f1); // more economic version from Yura.
     if (f1 * f1 + f2 * f2 > 1.f && f1 * f2 < 0.f) { // special cases of large rotations or large abs angles
       if (f2 > 0.f) {
         rot = constants::math::PI - rot; //
@@ -564,7 +564,7 @@ GPUd() typename TrackParametrization<value_T>::value_t TrackParametrization<valu
   value_t snp = mP[kSnp] + x2r;
   value_t phi = 999.;
   if (gpu::CAMath::Abs(snp) < constants::math::Almost1) {
-    phi = gpu::CAMath::ASin(snp) + getAlpha();
+    phi = math_utils::detail::asin<value_t>(snp) + getAlpha();
     math_utils::detail::bringTo02Pi<value_t>(phi);
   }
   return phi;
@@ -578,7 +578,7 @@ GPUd() typename TrackParametrization<value_T>::value_t TrackParametrization<valu
   value_t phi = 999.;
   auto y = getYAt(xk, b);
   if (y > -9998.) {
-    phi = gpu::CAMath::ATan2(y, xk) + getAlpha();
+    phi = math_utils::detail::atan2<value_t>(y, xk) + getAlpha();
     math_utils::detail::bringTo02Pi<value_t>(phi);
   }
   return phi;
@@ -592,7 +592,7 @@ GPUd() typename TrackParametrization<value_T>::value_t TrackParametrization<valu
   math_utils::detail::bringToPMPi<value_t>(alpha);
   value_t ca = 0, sa = 0;
   math_utils::detail::sincos(alpha - getAlpha(), sa, ca);
-  value_t snp = getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp)); // Improve precision
+  value_t snp = getSnp(), csp = math_utils::detail::sqrt<value_t>((1.f - snp) * (1.f + snp)); // Improve precision
   // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle direction in local frame is along the X axis
   if ((csp * ca + snp * sa) < 0.) {
     // LOGF(warning,"Rotation failed: local cos(phi) would become {:.2f}", csp * ca + snp * sa);
@@ -707,7 +707,7 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
       crv = 0.;
       break;
     }
-    value_t r0 = gpu::CAMath::Sqrt(circle.getCenterD2());
+    value_t r0 = math_utils::detail::sqrt<value_t>(circle.getCenterD2());
     if (r0 <= constants::math::Almost0) {
       return false; // the track is concentric to circle
     }
@@ -725,7 +725,7 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
     if (det < 0.) {
       return false; // does not reach raduis r
     }
-    det = gpu::CAMath::Sqrt(det);
+    det = math_utils::detail::sqrt<value_t>(det);
     //
     // the intersection happens in 2 points: {circle.xC+tR*C,circle.yC+tR*S}
     // with C=f*c0+-|s0|*det and S=f*s0-+c0 sign(s0)*det
@@ -794,7 +794,7 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
     if (dir == DirAuto) {
       return true;
     }
-    det = gpu::CAMath::Sqrt(det);
+    det = math_utils::detail::sqrt<value_t>(det);
     if (dir == DirOutward) { // along the track direction
       if (sn > 0.) {
         if (fy > det) {
@@ -819,7 +819,7 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
     if (det < 0.) {
       return false; // does not reach raduis r
     }
-    det = gpu::CAMath::Sqrt(det);
+    det = math_utils::detail::sqrt<value_t>(det);
     if (dir == DirAuto) {
       x = mX > 0. ? det : -det; // choose the solution requiring the smalest step
       return true;
@@ -837,13 +837,13 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
       }
     }
   } else { // general case of straight line
-    auto cs = gpu::CAMath::Sqrt((1. - sn) * (1. + sn));
+    auto cs = math_utils::detail::sqrt<value_t>((1. - sn) * (1. + sn));
     auto xsyc = mX * sn - fy * cs;
     auto det = (r - xsyc) * (r + xsyc);
     if (det < 0.) {
       return false; // does not reach raduis r
     }
-    det = gpu::CAMath::Sqrt(det);
+    det = math_utils::detail::sqrt<value_t>(det);
     auto xcys = mX * cs + fy * sn;
     auto t = -xcys;
     if (dir == DirAuto) {
@@ -879,10 +879,10 @@ GPUd() int TrackParametrization<value_T>::getELossSteps(value_t xrho, bool angle
   if (anglecorr) {
     const value_t csp2 = (1.f - getSnp()) * (1.f + getSnp()); // cos(phi)^2
     const value_t cst2I = (1.f + getTgl() * getTgl());        // 1/cos(lambda)^2
-    xrho *= gpu::CAMath::Sqrt(cst2I / csp2);
+    xrho *= math_utils::detail::sqrt<value_t>(cst2I / csp2);
   }
   const value_t p = getP(), massInv = 1.f / m;
-  const value_t e = gpu::CAMath::Sqrt(p * p + getPID().getMass2()), ekin = e - m;
+  const value_t e = math_utils::detail::sqrt<value_t>(p * p + getPID().getMass2()), ekin = e - m;
   value_t dedx = getdEdxBBOpt(p * massInv);
   const int charge2 = getAbsCharge() * getAbsCharge();
   if (charge2 != 1) {
@@ -911,12 +911,12 @@ GPUd() bool TrackParametrization<value_T>::correctForELoss(value_t xrho, bool an
     if (anglecorr) {
       value_t csp2 = (1.f - getSnp()) * (1.f + getSnp()); // cos(phi)^2
       value_t cst2I = (1.f + getTgl() * getTgl());        // 1/cos(lambda)^2
-      value_t angle = gpu::CAMath::Sqrt(cst2I / (csp2));
+      value_t angle = math_utils::detail::sqrt<value_t>(cst2I / (csp2));
       xrho *= angle;
     }
     int charge2 = getAbsCharge() * getAbsCharge();
     value_t p = getP(), p0 = p, p2 = p * p, e2 = p2 + getPID().getMass2(), massInv = 1.f / m, bg = p * massInv;
-    value_t e = gpu::CAMath::Sqrt(e2), ekin = e - m, dedx = getdEdxBBOpt(bg);
+    value_t e = math_utils::detail::sqrt<value_t>(e2), ekin = e - m, dedx = getdEdxBBOpt(bg);
 #ifdef _BB_NONCONST_CORR_
     value_t dedxDer = 0., dedx1 = dedx;
 #endif
@@ -941,13 +941,13 @@ GPUd() bool TrackParametrization<value_T>::correctForELoss(value_t xrho, bool an
         if (xrho < 0) {
           dedxDer = -dedxDer; // E.loss ( -> positive derivative)
         }
-        auto corrC = (gpu::CAMath::Exp(dedxDer) - 1.) / dedxDer;
+        auto corrC = (math_utils::detail::exp<value_t>(dedxDer) - 1.) / dedxDer;
         dE *= corrC;
       }
 #endif
       e += dE;
       if (e > m) { // stopped
-        p = gpu::CAMath::Sqrt(e * e - getPID().getMass2());
+        p = math_utils::detail::sqrt<value_t>(e * e - getPID().getMass2());
       } else {
         return false;
       }
diff --git a/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx b/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx
index 94091a007f94a..4f787d9d7647f 100644
--- a/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx
+++ b/DataFormats/Reconstruction/src/TrackParametrizationWithError.cxx
@@ -59,11 +59,11 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, valu
   if ((gpu::CAMath::Abs(f1) > constants::math::Almost1) || (gpu::CAMath::Abs(f2) > constants::math::Almost1)) {
     return false;
   }
-  value_t r1 = gpu::CAMath::Sqrt((1.f - f1) * (1.f + f1));
+  value_t r1 = math_utils::detail::sqrt<value_t>((1.f - f1) * (1.f + f1));
   if (gpu::CAMath::Abs(r1) < constants::math::Almost0) {
     return false;
   }
-  value_t r2 = gpu::CAMath::Sqrt((1.f - f2) * (1.f + f2));
+  value_t r2 = math_utils::detail::sqrt<value_t>((1.f - f2) * (1.f + f2));
   if (gpu::CAMath::Abs(r2) < constants::math::Almost0) {
     return false;
   }
@@ -85,7 +85,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, valu
     if (gpu::CAMath::Abs(arg) > constants::math::Almost1) {
       return false;
     }
-    value_t rot = gpu::CAMath::ASin(arg);           // more economic version from Yura.
+    value_t rot = math_utils::detail::asin<value_t>(arg); // more economic version from Yura.
     if (f1 * f1 + f2 * f2 > 1.f && f1 * f2 < 0.f) { // special cases of large rotations or large abs angles
       if (f2 > 0.f) {
         rot = constants::math::PI - rot; //
@@ -180,8 +180,8 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, Trac
   }
   value_t kb = bz * constants::math::B2C;
   // evaluate in double prec.
-  double snpRef0 = linRef0.getSnp(), cspRef0 = gpu::CAMath::Sqrt((1 - snpRef0) * (1 + snpRef0));
-  double snpRef1 = linRef1.getSnp(), cspRef1 = gpu::CAMath::Sqrt((1 - snpRef1) * (1 + snpRef1));
+  double snpRef0 = linRef0.getSnp(), cspRef0 = math_utils::detail::sqrt<value_t>((1 - snpRef0) * (1 + snpRef0));
+  double snpRef1 = linRef1.getSnp(), cspRef1 = math_utils::detail::sqrt<value_t>((1 - snpRef1) * (1 + snpRef1));
   double cspRef0Inv = 1 / cspRef0, cspRef1Inv = 1 / cspRef1, cc = cspRef0 + cspRef1, ccInv = 1 / cc, dy2dx = (snpRef0 + snpRef1) * ccInv;
   double dxccInv = dx * ccInv, hh = dxccInv * cspRef1Inv * (1 + cspRef0 * cspRef1 + snpRef0 * snpRef1), jj = dx * (dy2dx - snpRef1 * cspRef1Inv);
 
@@ -272,7 +272,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::rotate(value_t alpha)
   //
   value_t ca = 0, sa = 0;
   math_utils::detail::sincos(alpha - this->getAlpha(), sa, ca);
-  value_t snp = this->getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp)); // Improve precision
+  value_t snp = this->getSnp(), csp = math_utils::detail::sqrt<value_t>((1.f - snp) * (1.f + snp)); // Improve precision
   // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
   // direction in local frame is along the X axis
   if ((csp * ca + snp * sa) < 0) {
@@ -339,7 +339,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::rotate(value_t alpha, TrackP
   }
 
   // now rotate the current track
-  value_t snp = this->getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp)), updSnp = snp * ca - csp * sa;
+  value_t snp = this->getSnp(), csp = math_utils::detail::sqrt<value_t>((1.f - snp) * (1.f + snp)), updSnp = snp * ca - csp * sa;
   if ((csp * ca + snp * sa) < 0 || gpu::CAMath::Abs(updSnp) > constants::math::Almost1) {
     // LOGP(warning,"Rotation failed: local cos(phi) would become {:.2f}", csp * ca + snp * sa);
     return false;
@@ -350,7 +350,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::rotate(value_t alpha, TrackP
   this->setAlpha(alpha);
 
   // rotate covariance, accounting for the extra error from the rotated X
-  value_t snpRef0 = linRef0.getSnp(), cspRef0 = gpu::CAMath::Sqrt((value_t(1) - snpRef0) * (value_t(1) + snpRef0)); // original reference
+  value_t snpRef0 = linRef0.getSnp(), cspRef0 = math_utils::detail::sqrt<value_t>((value_t(1) - snpRef0) * (value_t(1) + snpRef0)); // original reference
   value_t snpRef1 = linRef1.getSnp(), cspRef1 = ca * cspRef0 + sa * snpRef0;                                        // rotated reference
   value_t rr = cspRef1 / cspRef0;                                                                                   // cos1_ref / cos0_ref
 
@@ -415,7 +415,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateToDCA(const o2::dat
   // propagate track to DCA to the vertex
   value_t sn, cs, alp = this->getAlpha();
   o2::math_utils::detail::sincos(alp, sn, cs);
-  value_t x = this->getX(), y = this->getY(), snp = this->getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp));
+  value_t x = this->getX(), y = this->getY(), snp = this->getSnp(), csp = math_utils::detail::sqrt<value_t>((1.f - snp) * (1.f + snp));
   value_t xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
   x -= xv;
   y -= yv;
@@ -430,8 +430,8 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateToDCA(const o2::dat
   }
   value_t crv = this->getCurvature(b);
   value_t tgfv = -(crv * x - snp) / (crv * y + csp);
-  sn = tgfv / gpu::CAMath::Sqrt(1.f + tgfv * tgfv);
-  cs = gpu::CAMath::Sqrt((1.f - sn) * (1.f + sn));
+  sn = tgfv / math_utils::detail::sqrt<value_t>(1.f + tgfv * tgfv);
+  cs = math_utils::detail::sqrt<value_t>((1.f - sn) * (1.f + sn));
   cs = (gpu::CAMath::Abs(tgfv) > constants::math::Almost0) ? sn / tgfv : constants::math::Almost1;
 
   x = xv * cs + yv * sn;
@@ -439,7 +439,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateToDCA(const o2::dat
   xv = x;
 
   auto tmpT(*this); // operate on the copy to recover after the failure
-  alp += gpu::CAMath::ASin(sn);
+  alp += math_utils::detail::asin<value_t>(sn);
   if (!tmpT.rotate(alp) || !tmpT.propagateTo(xv, b)) {
 #if !defined(GPUCA_ALIGPUCODE)
     LOG(debug) << "failed to propagate to alpha=" << alp << " X=" << xv << vtx << " | Track is: " << tmpT.asString();
@@ -485,9 +485,9 @@ GPUd() void TrackParametrizationWithError<value_T>::set(const dim3_t& xyz, const
   value_t radPos2 = xyz[0] * xyz[0] + xyz[1] * xyz[1];
   value_t alp = 0;
   if (sectorAlpha || radPos2 < 1) {
-    alp = gpu::CAMath::ATan2(pxpypz[1], pxpypz[0]);
+    alp = math_utils::detail::atan2<value_t>(pxpypz[1], pxpypz[0]);
   } else {
-    alp = gpu::CAMath::ATan2(xyz[1], xyz[0]);
+    alp = math_utils::detail::atan2<value_t>(xyz[1], xyz[0]);
   }
   if (sectorAlpha) {
     alp = math_utils::detail::angle2Alpha<value_t>(alp);
@@ -498,7 +498,7 @@ GPUd() void TrackParametrizationWithError<value_T>::set(const dim3_t& xyz, const
   // protection against cosp<0
   if (cs * pxpypz[0] + sn * pxpypz[1] < 0) {
     LOG(debug) << "alpha from phiPos() will invalidate this track parameters, overriding to alpha from phi()";
-    alp = gpu::CAMath::ATan2(pxpypz[1], pxpypz[0]);
+    alp = math_utils::detail::atan2<value_t>(pxpypz[1], pxpypz[0]);
     if (sectorAlpha) {
       alp = math_utils::detail::angle2Alpha<value_t>(alp);
     }
@@ -529,7 +529,7 @@ GPUd() void TrackParametrizationWithError<value_T>::set(const dim3_t& xyz, const
   math_utils::detail::rotateZ<value_t>(mom, -alp);
   //
   const value_t pt2 = mom[0] * mom[0] + mom[1] * mom[1];
-  const value_t pt = gpu::CAMath::Sqrt(pt2);
+  const value_t pt = math_utils::detail::sqrt<value_t>(pt2);
   const value_t ptI = 1.f / pt;
   this->setX(ver[0]);
   this->setAlpha(alp);
@@ -628,19 +628,19 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, cons
   if ((gpu::CAMath::Abs(f1) > constants::math::Almost1) || (gpu::CAMath::Abs(f2) > constants::math::Almost1)) {
     return false;
   }
-  value_t r1 = gpu::CAMath::Sqrt((1.f - f1) * (1.f + f1));
+  value_t r1 = math_utils::detail::sqrt<value_t>((1.f - f1) * (1.f + f1));
   if (gpu::CAMath::Abs(r1) < constants::math::Almost0) {
     return false;
   }
-  value_t r2 = gpu::CAMath::Sqrt((1.f - f2) * (1.f + f2));
+  value_t r2 = math_utils::detail::sqrt<value_t>((1.f - f2) * (1.f + f2));
   if (gpu::CAMath::Abs(r2) < constants::math::Almost0) {
     return false;
   }
   double r1pr2Inv = 1. / (r1 + r2), r2inv = 1. / r2, r1inv = 1. / r1;
   double dy2dx = (f1 + f2) * r1pr2Inv, dx2r1pr2 = dx * r1pr2Inv;
-  value_t step = (gpu::CAMath::Abs(x2r) < 0.05f) ? value_t(dx * gpu::CAMath::Abs(r2 + f2 * dy2dx))                                          // chord
-                                                 : 2.f * gpu::CAMath::ASin(0.5f * dx * gpu::CAMath::Sqrt(1.f + dy2dx * dy2dx) * crv) / crv; // arc
-  step *= gpu::CAMath::Sqrt(1.f + this->getTgl() * this->getTgl());
+  value_t step = (gpu::CAMath::Abs(x2r) < 0.05f) ? value_t(dx * gpu::CAMath::Abs(r2 + f2 * dy2dx))                                                                          // chord
+                                                 : 2.f * math_utils::detail::asin<value_t>(0.5f * dx * math_utils::detail::sqrt<value_t>(1.f + dy2dx * dy2dx) * crv) / crv; // arc
+  step *= math_utils::detail::sqrt<value_t>(1.f + this->getTgl() * this->getTgl());
   //
   // get the track x,y,z,px/p,py/p,pz/p,p,sinAlpha,cosAlpha in the Global System
   std::array<value_t, 9> vecLab{0.f};
@@ -699,13 +699,13 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, cons
 
   // Rotate to the system where Bx=By=0.
   value_t bxy2 = b[0] * b[0] + b[1] * b[1];
-  value_t bt = gpu::CAMath::Sqrt(bxy2);
+  value_t bt = math_utils::detail::sqrt<value_t>(bxy2);
   value_t cosphi = 1.f, sinphi = 0.f;
   if (bt > constants::math::Almost0) {
     cosphi = b[0] / bt;
     sinphi = b[1] / bt;
   }
-  value_t bb = gpu::CAMath::Sqrt(bxy2 + b[2] * b[2]);
+  value_t bb = math_utils::detail::sqrt<value_t>(bxy2 + b[2] * b[2]);
   value_t costet = 1., sintet = 0.;
   if (bb > constants::math::Almost0) {
     costet = b[2] / bb;
@@ -754,7 +754,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, cons
   }
 
   // Calculate the track parameters
-  t = 1.f / gpu::CAMath::Sqrt(vecLab[3] * vecLab[3] + vecLab[4] * vecLab[4]);
+  t = 1.f / math_utils::detail::sqrt<value_t>(vecLab[3] * vecLab[3] + vecLab[4] * vecLab[4]);
   this->setX(xk);
   this->setY(y);
   this->setZ(z);
@@ -802,13 +802,13 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, Trac
   if ((gpu::CAMath::Abs(snpRef0) > constants::math::Almost1) || (gpu::CAMath::Abs(snpRef1) > constants::math::Almost1)) {
     return false;
   }
-  value_t cspRef0 = gpu::CAMath::Sqrt((1 - snpRef0) * (1 + snpRef0)), cspRef1 = gpu::CAMath::Sqrt((1 - snpRef1) * (1 + snpRef1));
+  value_t cspRef0 = math_utils::detail::sqrt<value_t>((1 - snpRef0) * (1 + snpRef0)), cspRef1 = math_utils::detail::sqrt<value_t>((1 - snpRef1) * (1 + snpRef1));
   if (gpu::CAMath::Abs(cspRef0) < constants::math::Almost0 || gpu::CAMath::Abs(cspRef1) < constants::math::Almost0) {
     return false;
   }
   value_t cspRef0Inv = value_t(1) / cspRef0, cspRef1Inv = value_t(1) / cspRef1, cc = cspRef0 + cspRef1, ccInv = value_t(1) / cc, dy2dx = (snpRef0 + snpRef1) * ccInv;
-  value_t step = (gpu::CAMath::Abs(crv * dx) < 0.05f) ? dx * (cspRef1 + snpRef1 * dy2dx) : 2. * gpu::CAMath::ASin(0.5 * dx * gpu::CAMath::Sqrt(1.f + dy2dx * dy2dx) * crv) / crv; // arc
-  step *= gpu::CAMath::Sqrt(1.f + linRef0.getTgl() * linRef0.getTgl());
+  value_t step = (gpu::CAMath::Abs(crv * dx) < 0.05f) ? dx * (cspRef1 + snpRef1 * dy2dx) : 2. * math_utils::detail::asin<value_t>(0.5 * dx * math_utils::detail::sqrt<value_t>(1.f + dy2dx * dy2dx) * crv) / crv; // arc
+  step *= math_utils::detail::sqrt<value_t>(1.f + linRef0.getTgl() * linRef0.getTgl());
 
   //
   // get the track x,y,z,px/p,py/p,pz/p,p,sinAlpha,cosAlpha in the Global System
@@ -819,13 +819,13 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, Trac
   //
   // Rotate to the system where Bx=By=0.
   value_t bxy2 = b[0] * b[0] + b[1] * b[1];
-  value_t bt = gpu::CAMath::Sqrt(bxy2);
+  value_t bt = math_utils::detail::sqrt<value_t>(bxy2);
   value_t cosphi = 1.f, sinphi = 0.f;
   if (bt > constants::math::Almost0) {
     cosphi = b[0] / bt;
     sinphi = b[1] / bt;
   }
-  value_t bb = gpu::CAMath::Sqrt(bxy2 + b[2] * b[2]);
+  value_t bb = math_utils::detail::sqrt<value_t>(bxy2 + b[2] * b[2]);
   value_t costet = 1., sintet = 0.;
   if (bb > constants::math::Almost0) {
     costet = b[2] / bb;
@@ -875,7 +875,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, Trac
 
   // Calculate the track parameters
   auto linRef1 = linRef0;
-  t = 1.f / gpu::CAMath::Sqrt(vecLab[3] * vecLab[3] + vecLab[4] * vecLab[4]);
+  t = 1.f / math_utils::detail::sqrt<value_t>(vecLab[3] * vecLab[3] + vecLab[4] * vecLab[4]);
   linRef1.setX(xk);
   linRef1.setY(y);
   linRef1.setZ(z);
@@ -884,7 +884,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, Trac
   linRef1.setQ2Pt(q * t / vecLab[6]);
 
   // recalculate parameters of the transported ref track needed for transport of this:
-  cspRef1 = gpu::CAMath::Sqrt((1 - snpRef1) * (1 + snpRef1));
+  cspRef1 = math_utils::detail::sqrt<value_t>((1 - snpRef1) * (1 + snpRef1));
   cspRef1Inv = value_t(1) / cspRef1;
   cc = cspRef0 + cspRef1;
   ccInv = value_t(1) / cc;
@@ -968,7 +968,7 @@ GPUd() void TrackParametrizationWithError<value_T>::checkCorrelations()
       auto sig2 = mC[DiagMap[i]] * mC[DiagMap[j]];
       auto& cov = mC[CovarMap[i][j]];
       if (cov * cov >= MaxCorr * sig2) { // constrain correlation
-        cov = gpu::CAMath::Sqrt(sig2) * (cov > 0. ? MaxCorr : -MaxCorr);
+        cov = math_utils::detail::sqrt<value_t>(sig2) * (cov > 0. ? MaxCorr : -MaxCorr);
       }
     }
   }
@@ -984,7 +984,7 @@ GPUd() void TrackParametrizationWithError<value_T>::checkCovariance()
 
   mC[kSigY2] = gpu::CAMath::Abs(mC[kSigY2]);
   if (mC[kSigY2] > kCY2max) {
-    value_t scl = gpu::CAMath::Sqrt(kCY2max / mC[kSigY2]);
+    value_t scl = math_utils::detail::sqrt<value_t>(kCY2max / mC[kSigY2]);
     mC[kSigY2] = kCY2max;
     mC[kSigZY] *= scl;
     mC[kSigSnpY] *= scl;
@@ -993,7 +993,7 @@ GPUd() void TrackParametrizationWithError<value_T>::checkCovariance()
   }
   mC[kSigZ2] = gpu::CAMath::Abs(mC[kSigZ2]);
   if (mC[kSigZ2] > kCZ2max) {
-    value_t scl = gpu::CAMath::Sqrt(kCZ2max / mC[kSigZ2]);
+    value_t scl = math_utils::detail::sqrt<value_t>(kCZ2max / mC[kSigZ2]);
     mC[kSigZ2] = kCZ2max;
     mC[kSigZY] *= scl;
     mC[kSigSnpZ] *= scl;
@@ -1002,7 +1002,7 @@ GPUd() void TrackParametrizationWithError<value_T>::checkCovariance()
   }
   mC[kSigSnp2] = gpu::CAMath::Abs(mC[kSigSnp2]);
   if (mC[kSigSnp2] > kCSnp2max) {
-    value_t scl = gpu::CAMath::Sqrt(kCSnp2max / mC[kSigSnp2]);
+    value_t scl = math_utils::detail::sqrt<value_t>(kCSnp2max / mC[kSigSnp2]);
     mC[kSigSnp2] = kCSnp2max;
     mC[kSigSnpY] *= scl;
     mC[kSigSnpZ] *= scl;
@@ -1011,7 +1011,7 @@ GPUd() void TrackParametrizationWithError<value_T>::checkCovariance()
   }
   mC[kSigTgl2] = gpu::CAMath::Abs(mC[kSigTgl2]);
   if (mC[kSigTgl2] > kCTgl2max) {
-    value_t scl = gpu::CAMath::Sqrt(kCTgl2max / mC[kSigTgl2]);
+    value_t scl = math_utils::detail::sqrt<value_t>(kCTgl2max / mC[kSigTgl2]);
     mC[kSigTgl2] = kCTgl2max;
     mC[kSigTglY] *= scl;
     mC[kSigTglZ] *= scl;
@@ -1020,7 +1020,7 @@ GPUd() void TrackParametrizationWithError<value_T>::checkCovariance()
   }
   mC[kSigQ2Pt2] = gpu::CAMath::Abs(mC[kSigQ2Pt2]);
   if (mC[kSigQ2Pt2] > kC1Pt2max) {
-    value_t scl = gpu::CAMath::Sqrt(kC1Pt2max / mC[kSigQ2Pt2]);
+    value_t scl = math_utils::detail::sqrt<value_t>(kC1Pt2max / mC[kSigQ2Pt2]);
     mC[kSigQ2Pt2] = kC1Pt2max;
     mC[kSigQ2PtY] *= scl;
     mC[kSigQ2PtZ] *= scl;
@@ -1421,14 +1421,14 @@ GPUd() bool TrackParametrizationWithError<value_T>::correctForMaterial(value_t x
   value_t csp2 = (1.f - this->getSnp()) * (1.f + this->getSnp()); // cos(phi)^2
   value_t cst2I = (1.f + this->getTgl() * this->getTgl());        // 1/cos(lambda)^2
   if (anglecorr) {                                                // Apply angle correction, if requested
-    value_t angle = gpu::CAMath::Sqrt(cst2I / csp2);
+    value_t angle = math_utils::detail::sqrt<value_t>(cst2I / csp2);
     x2x0 *= angle;
     xrho *= angle;
   }
   auto m = this->getPID().getMass();
   int charge2 = this->getAbsCharge() * this->getAbsCharge();
   value_t p = this->getP(), p0 = p, p02 = p * p, e2 = p02 + this->getPID().getMass2(), massInv = 1.f / m, bg = p * massInv, dETot = 0.f;
-  value_t e = gpu::CAMath::Sqrt(e2), e0 = e;
+  value_t e = math_utils::detail::sqrt<value_t>(e2), e0 = e;
   if (m > 0 && xrho != 0.f) {
     value_t ekin = e - m, dedx = this->getdEdxBBOpt(bg);
 #ifdef _BB_NONCONST_CORR_
@@ -1455,13 +1455,13 @@ GPUd() bool TrackParametrizationWithError<value_T>::correctForMaterial(value_t x
         if (xrho < 0) {
           dedxDer = -dedxDer; // E.loss ( -> positive derivative)
         }
-        auto corrC = (gpu::CAMath::Exp(dedxDer) - 1.f) / dedxDer;
+        auto corrC = (math_utils::detail::exp<value_t>(dedxDer) - 1.f) / dedxDer;
         dE *= corrC;
       }
 #endif
       e += dE;
       if (e > m) { // stopped
-        p = gpu::CAMath::Sqrt(e * e - this->getPID().getMass2());
+        p = math_utils::detail::sqrt<value_t>(e * e - this->getPID().getMass2());
       } else {
         return false;
       }
@@ -1518,7 +1518,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::correctForMaterial(value_t x
 
   // the energy loss correction contribution to cov.matrix: approximate energy loss fluctuation (M.Ivanov)
   constexpr value_t knst = 0.0007f; // To be tuned.
-  value_t sigmadE = knst * gpu::CAMath::Sqrt(gpu::CAMath::Abs(dETot)) * e0 / p02 * this->getCharge2Pt();
+  value_t sigmadE = knst * math_utils::detail::sqrt<value_t>(gpu::CAMath::Abs(dETot)) * e0 / p02 * this->getCharge2Pt();
   cC44 += sigmadE * sigmadE;
 
   // Applying the corrections*****************************
@@ -1552,7 +1552,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::correctForMaterial(TrackPara
   value_t csp2 = (1.f - linRef.getSnp()) * (1.f + linRef.getSnp()); // cos(phi)^2
   value_t cst2I = (1.f + linRef.getTgl() * linRef.getTgl());        // 1/cos(lambda)^2
   if (anglecorr) {                                                  // Apply angle correction, if requested
-    value_t angle = gpu::CAMath::Sqrt(cst2I / csp2);
+    value_t angle = math_utils::detail::sqrt<value_t>(cst2I / csp2);
     x2x0 *= angle;
     xrho *= angle;
   }
@@ -1560,7 +1560,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::correctForMaterial(TrackPara
   auto m = pid.getMass();
   int charge2 = linRef.getAbsCharge() * linRef.getAbsCharge();
   value_t p = linRef.getP(), p0 = p, p02 = p * p, e2 = p02 + pid.getMass2(), massInv = 1.f / m, bg = p * massInv, dETot = 0.f;
-  value_t e = gpu::CAMath::Sqrt(e2), e0 = e;
+  value_t e = math_utils::detail::sqrt<value_t>(e2), e0 = e;
   if (m > 0 && xrho != 0.f) {
     value_t ekin = e - m, dedx = this->getdEdxBBOpt(bg);
 #ifdef _BB_NONCONST_CORR_
@@ -1587,13 +1587,13 @@ GPUd() bool TrackParametrizationWithError<value_T>::correctForMaterial(TrackPara
         if (xrho < 0) {
           dedxDer = -dedxDer; // E.loss ( -> positive derivative)
         }
-        auto corrC = (gpu::CAMath::Exp(dedxDer) - 1.f) / dedxDer;
+        auto corrC = (math_utils::detail::exp<value_t>(dedxDer) - 1.f) / dedxDer;
         dE *= corrC;
       }
 #endif
       e += dE;
       if (e > m) { // stopped
-        p = gpu::CAMath::Sqrt(e * e - pid.getMass2());
+        p = math_utils::detail::sqrt<value_t>(e * e - pid.getMass2());
       } else {
         return false;
       }
@@ -1650,7 +1650,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::correctForMaterial(TrackPara
 
   // the energy loss correction contribution to cov.matrix: approximate energy loss fluctuation (M.Ivanov)
   constexpr value_t knst = 0.0007f; // To be tuned.
-  value_t sigmadE = knst * gpu::CAMath::Sqrt(gpu::CAMath::Abs(dETot)) * e0 / p02 * linRef.getCharge2Pt();
+  value_t sigmadE = knst * math_utils::detail::sqrt<value_t>(gpu::CAMath::Abs(dETot)) * e0 / p02 * linRef.getCharge2Pt();
   cC44 += sigmadE * sigmadE;
 
   // Applying the corrections*****************************
@@ -1695,7 +1695,7 @@ GPUd() bool TrackParametrizationWithError<value_T>::getCovXYZPxPyPzGlo(std::arra
   value_t sn = 0.f, cs = 0.f;
   o2::math_utils::detail::sincos(this->getAlpha(), sn, cs);
   const value_t snp = this->getSnp();
-  const value_t csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp));
+  const value_t csp = math_utils::detail::sqrt<value_t>((1.f - snp) * (1.f + snp));
   const value_t pXLoc = pt * csp;
   const value_t pYLoc = pt * snp;
   const value_t pZ = pt * this->getTgl();
@@ -1818,7 +1818,7 @@ bool TrackParametrizationWithError<value_T>::toFwdTrackParCov(TrackParCovFwd& t)
       gpu::CAMath::Abs(this->getTgl()) <= o2::constants::math::Almost0) {
     return false;
   }
-  value_T csa, sna, csP, snP, csp = gpu::CAMath::Sqrt((1. - this->getSnp()) * (1. + this->getSnp()));
+  value_T csa, sna, csP, snP, csp = math_utils::detail::sqrt<value_t>((1. - this->getSnp()) * (1. + this->getSnp()));
   math_utils::detail::sincos(value_T(this->getAlpha()), sna, csa);
   math_utils::detail::sincos(value_T(t.getPhi()), snP, csP);
   /*
diff --git a/GPU/Common/GPUCommonMath.h b/GPU/Common/GPUCommonMath.h
index 54b6fcfabcb9e..f32cf61a057ae 100644
--- a/GPU/Common/GPUCommonMath.h
+++ b/GPU/Common/GPUCommonMath.h
@@ -101,6 +101,7 @@ class GPUCommonMath
   GPUd() constexpr static float Log(float x);
   GPUd() constexpr static double Logd(double x);
   GPUd() constexpr static float Exp(float x);
+  GPUd() constexpr static double Expd(double x);
   GPUhdni() constexpr static float Copysign(float x, float y) { return GPUCA_CHOICE(std::copysignf(x, y), copysignf(x, y), copysign(x, y)); }
   GPUd() constexpr static double Copysignd(double x, double y);
   GPUd() constexpr static float TwoPi() { return 6.2831853f; }
@@ -377,6 +378,7 @@ GPUdi() constexpr double GPUCommonMath::Sind(double x) { return GPUCA_CHOICE(std
 GPUdi() constexpr double GPUCommonMath::Cosd(double x) { return GPUCA_CHOICE(std::cos(x), cos(x), cos(x)); }
 GPUdi() constexpr double GPUCommonMath::Tand(double x) { return GPUCA_CHOICE(std::tan(x), tan(x), tan(x)); }
 GPUdi() constexpr double GPUCommonMath::Logd(double x) { return GPUCA_CHOICE(std::log(x), log(x), log(x)); }
+GPUdi() constexpr double GPUCommonMath::Expd(double x) { return GPUCA_CHOICE(std::exp(x), exp(x), exp(x)); }
 GPUdi() constexpr double GPUCommonMath::Copysignd(double x, double y) { return GPUCA_CHOICE(std::copysign(x, y), copysign(x, y), copysign(x, y)); }
 #endif