TrackParametrization: use double math functions for TrackParD (#15990)

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 cd9835a (switched back to CAMath in bf65d3c)
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>

Author: matthias-kleiner

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

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()));
 }
 
 //_______________________________________________________

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

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) {