Add GPUdoubleCalc and use it for the double-precision track math

The track propagation computes its Jacobian and covariance intermediates in
double even when the track is float, because the terms cancel:
jj = dx * (dy2dx - f2 * r2inv) is a difference of nearly equal quantities. MSL
has no double, and plain float is not an option there.

GPUdoubleCalc is compensated two-float arithmetic: the value is mHi + mLo, so
the rounding error of each operation is carried explicitly. Over 500k samples of
the isolated jj expression, float intermediates reach 3.6e-2 relative error with
23 samples past 1e-3, while the two-float form lands at 5.9e-8, half a float ulp,
with none past 1e-3. End to end over 524288 tracks with strongly correlated
covariances it keeps the covariance within 1.8e-6 of sqrt(C_ii C_jj) of the CPU
double result, where plain float is at 1.2e-5 and the CPU double is itself
1.1e-6 away from a cancellation-free reference.

Off Metal it is a plain double, so nothing else moves: substituting the alias
back reproduces the previous source exactly, and the files compile unchanged
with their real build flags.

Cost on an M1 Max, over the whole propagateTo kernel with fast math off: 2.85x
the plain float path. Metal has no fp64 at all, so the comparison is against not
compiling.

The type is defined for every backend so the arithmetic can be tested on the
host, where GPUCA_FORCE_DOUBLECALC and GPUCA_FORCE_FLOATCALC select the
representation explicitly.

Author: ktf

DataFormats/Reconstruction/include/ReconstructionDataFormats/HelixHelper.h

@@ -17,6 +17,7 @@
 #define _ALICEO2_HELIX_HELPER_
 
 #include "CommonConstants/MathConstants.h"
+#include "GPUCommonDouble.h"
 #include "MathUtils/Utils.h"
 #include "MathUtils/Primitive2D.h"
 
@@ -247,8 +248,8 @@ struct CrossInfo {
     auto tgp = trcL.getSnp() * cspi;
     float kx = traxL.c - traxL.s * tgp;
     float ky = traxL.s + traxL.c * tgp;
-    double dk = dx * kx + dy * ky;
-    double det = dk * dk - cspi2 * (dx * dx + dy * dy - traxH.rC * traxH.rC);
+    o2::gpu::GPUdoubleCalc dk = dx * kx + dy * ky;
+    o2::gpu::GPUdoubleCalc det = dk * dk - cspi2 * (dx * dx + dy * dy - traxH.rC * traxH.rC);
     if (det > 0) { // 2 crossings
       det = o2::gpu::GPUCommonMath::Sqrt(det);
       float t0 = (-dk + det) * cspi2;

DataFormats/Reconstruction/src/TrackParametrization.cxx

@@ -15,6 +15,7 @@
 /// @brief
 
 #include "ReconstructionDataFormats/TrackParametrization.h"
+#include "GPUCommonDouble.h"
 #include "ReconstructionDataFormats/Vertex.h"
 #include "ReconstructionDataFormats/DCA.h"
 #include <MathUtils/Cartesian.h>
@@ -477,7 +478,7 @@ GPUd() bool TrackParametrization<value_T>::getYZAt(value_t xk, value_t b, value_
   if (gpu::CAMath::Abs(r2) < constants::math::Almost0) {
     return false;
   }
-  double dy2dx = (f1 + f2) / (r1 + r2);
+  GPUdoubleCalc dy2dx = (f1 + f2) / (r1 + r2);
   y += dx * dy2dx;
   if (gpu::CAMath::Abs(x2r) < 0.05f) {
     z += dx * (r2 + f2 * dy2dx) * getTgl();
@@ -692,7 +693,7 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
   // DirOutward (==1) - go along the track (increasing mX)
   // DirInward (==-1) - go backward (decreasing mX)
   //
-  const double fy = mP[0], sn = mP[2];
+  const GPUdoubleCalc fy = mP[0], sn = mP[2];
   const value_t kEps = 1.e-6;
   //
   if (gpu::CAMath::Abs(getSnp()) > constants::math::Almost1) {
@@ -711,7 +712,7 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
     if (r0 <= constants::math::Almost0) {
       return false; // the track is concentric to circle
     }
-    double tR2r0 = 1., g = 0., tmp = 0.;
+    GPUdoubleCalc tR2r0 = 1., g = 0., tmp = 0.;
     if (gpu::CAMath::Abs(circle.rC - r0) > kEps) {
       tR2r0 = circle.rC / r0;
       g = 0.5f * (r * r / (r0 * circle.rC) - tR2r0 - 1.f / tR2r0);
@@ -786,7 +787,7 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
   }
   // this is a straight track
   if (gpu::CAMath::Abs(sn) >= constants::math::Almost1) { // || to Y axis
-    double det = (r - mX) * (r + mX);
+    GPUdoubleCalc det = (r - mX) * (r + mX);
     if (det < 0.f) {
       return false; // does not reach raduis r
     }
@@ -815,7 +816,7 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
       }
     }
   } else if (gpu::CAMath::Abs(sn) <= constants::math::Almost0) { // || to X axis
-    double det = (r - fy) * (r + fy);
+    GPUdoubleCalc det = (r - fy) * (r + fy);
     if (det < 0.) {
       return false; // does not reach raduis r
     }