Robust propagation to certain R

shahor02 · shahor02 · commit 074e2ad36235 · 2025-09-28T23:25:04.000+04:00
o2::base::Propagator got method propagateToR(track, r, bool bzOnly, maxSnp ...) which propagates the track (either TrackPar or TrackParCov)
to requested radius (not X!) if it is reachable in any tracking frame. The rotations to avoid track snp goind outside of the allowed range are done automatically.

Moved HelixHelper to ReconstructionDataFormats.
diff --git a/Common/DCAFitter/CMakeLists.txt b/Common/DCAFitter/CMakeLists.txt
@@ -22,8 +22,7 @@ o2_add_library(DCAFitter
                                      O2::DetectorsBase)
 
 o2_target_root_dictionary(DCAFitter
-                          HEADERS include/DCAFitter/HelixHelper.h
-                                  include/DCAFitter/DCAFitterN.h
+                          HEADERS include/DCAFitter/DCAFitterN.h
                                   include/DCAFitter/FwdDCAFitterN.h)
 
 if (OpenMP_CXX_FOUND)
diff --git a/Common/DCAFitter/include/DCAFitter/DCAFitterN.h b/Common/DCAFitter/include/DCAFitter/DCAFitterN.h
@@ -17,7 +17,7 @@
 #ifndef _ALICEO2_DCA_FITTERN_
 #define _ALICEO2_DCA_FITTERN_
 
-#include "DCAFitter/HelixHelper.h"
+#include "ReconstructionDataFormats/HelixHelper.h"
 #include "DetectorsBase/Propagator.h"
 #include "MathUtils/Cartesian.h"
 #include "ReconstructionDataFormats/Track.h"
diff --git a/Common/DCAFitter/include/DCAFitter/FwdDCAFitterN.h b/Common/DCAFitter/include/DCAFitter/FwdDCAFitterN.h
@@ -20,7 +20,7 @@
 #include "MathUtils/Cartesian.h"
 #include "ReconstructionDataFormats/TrackFwd.h"
 #include "ReconstructionDataFormats/Track.h"
-#include "DCAFitter/HelixHelper.h"
+#include "ReconstructionDataFormats/HelixHelper.h"
 #include <TRandom.h>
 #include "DetectorsBase/Propagator.h"
 #include "DetectorsBase/GeometryManager.h"
diff --git a/Common/DCAFitter/src/DCAFitterLinkDef.h b/Common/DCAFitter/src/DCAFitterLinkDef.h
@@ -18,9 +18,6 @@
 #pragma link C++ class o2::vertexing::DCAFitterN < 2, o2::track::TrackParCov> + ;
 #pragma link C++ class o2::vertexing::DCAFitterN < 3, o2::track::TrackParCov> + ;
 
-#pragma link C++ class o2::track::TrackAuxPar + ;
-#pragma link C++ class o2::track::CrossInfo + ;
-
 #pragma link C++ function o2::vertexing::DCAFitter2::process(const o2::track::TrackParCov&, const o2::track::TrackParCov&);
 #pragma link C++ function o2::vertexing::DCAFitter3::process(const o2::track::TrackParCov&, const o2::track::TrackParCov&, const o2::track::TrackParCov&);
 
diff --git a/Common/Field/include/Field/MagFieldFast.h b/Common/Field/include/Field/MagFieldFast.h
@@ -57,7 +57,7 @@ class MagFieldFast
   bool Field(const math_utils::Point3D<double> xyz, double bxyz[3]) const;
   bool GetBcomp(EDim comp, const double xyz[3], double& b) const;
   bool GetBcomp(EDim comp, const float xyz[3], float& b) const;
-  bool GetBcomp(EDim comp, const math_utils::Point3D<float> xyz, double& b) const;
+  bool GetBcomp(EDim comp, const math_utils::Point3D<double> xyz, double& b) const;
   bool GetBcomp(EDim comp, const math_utils::Point3D<float> xyz, float& b) const;
 
   bool GetBx(const double xyz[3], double& bx) const { return GetBcomp(kX, xyz, bx); }
@@ -66,6 +66,8 @@ class MagFieldFast
   bool GetBy(const float xyz[3], float& by) const { return GetBcomp(kY, xyz, by); }
   bool GetBz(const double xyz[3], double& bz) const { return GetBcomp(kZ, xyz, bz); }
   bool GetBz(const float xyz[3], float& bz) const { return GetBcomp(kZ, xyz, bz); }
+  bool GetBz(const math_utils::Point3D<double> xyz, double& bz) const { return GetBcomp(kZ, xyz, bz); }
+  bool GetBz(const math_utils::Point3D<float> xyz, float& bz) const { return GetBcomp(kZ, xyz, bz); }
   void setFactorSol(float v = 1.f) { mFactorSol = v; }
   float getFactorSol() const { return mFactorSol; }
 
diff --git a/Common/Field/src/MagFieldFast.cxx b/Common/Field/src/MagFieldFast.cxx
@@ -145,7 +145,7 @@ bool MagFieldFast::GetBcomp(EDim comp, const double xyz[3], double& b) const
 }
 
 //_______________________________________________________________________
-bool MagFieldFast::GetBcomp(EDim comp, const math_utils::Point3D<float> xyz, double& b) const
+bool MagFieldFast::GetBcomp(EDim comp, const math_utils::Point3D<double> xyz, double& b) const
 {
   // get field
   int zSeg, rSeg, quadrant;
diff --git a/DataFormats/Reconstruction/CMakeLists.txt b/DataFormats/Reconstruction/CMakeLists.txt
@@ -73,6 +73,7 @@ o2_target_root_dictionary(
           include/ReconstructionDataFormats/BCRange.h
           include/ReconstructionDataFormats/TrackHMP.h
           include/ReconstructionDataFormats/MatchInfoHMP.h
+          include/ReconstructionDataFormats/HelixHelper.h
   )
 
 o2_add_test(Vertex
diff --git a/DataFormats/Reconstruction/include/ReconstructionDataFormats/HelixHelper.h b/DataFormats/Reconstruction/include/ReconstructionDataFormats/HelixHelper.h
diff --git a/DataFormats/Reconstruction/src/ReconstructionDataFormatsLinkDef.h b/DataFormats/Reconstruction/src/ReconstructionDataFormatsLinkDef.h
@@ -117,4 +117,7 @@
 #pragma link C++ class o2::dataformats::StrangeTrack + ;
 #pragma link C++ class std::vector < o2::dataformats::StrangeTrack> + ;
 
+#pragma link C++ class o2::track::TrackAuxPar + ;
+#pragma link C++ class o2::track::CrossInfo + ;
+
 #endif
diff --git a/DataFormats/Reconstruction/src/TrackParametrization.cxx b/DataFormats/Reconstruction/src/TrackParametrization.cxx
@@ -651,7 +651,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 auto fy = mP[0], sn = mP[2];
+  const double fy = mP[0], sn = mP[2];
   const value_t kEps = 1.e-6;
   //
   if (gpu::CAMath::Abs(getSnp()) > constants::math::Almost1) {
@@ -670,18 +670,18 @@ 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
     }
-    value_t tR2r0 = 1.f, g = 0.f, tmp = 0.f;
+    double 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);
       tmp = 1.f + g * tR2r0;
     } else {
       tR2r0 = 1.0;
-      g = 0.5f * r * r / (r0 * circle.rC) - 1.f;
-      tmp = 0.5f * r * r / (r0 * r0);
+      g = 0.5 * r * r / (r0 * circle.rC) - 1.;
+      tmp = 0.5 * r * r / (r0 * r0);
     }
-    value_t det = (1.f - g) * (1.f + g);
-    if (det < 0.f) {
+    auto det = (1. - g) * (1. + g);
+    if (det < 0.) {
       return false; // does not reach raduis r
     }
     det = gpu::CAMath::Sqrt(det);
@@ -691,42 +691,44 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
     // where s0 and c0 make direction for the circle center (=circle.xC/r0 and circle.yC/r0)
     //
     x = circle.xC * tmp;
-    value_t y = circle.yC * tmp;
+    auto y = circle.yC * tmp;
     if (gpu::CAMath::Abs(circle.yC) > constants::math::Almost0) { // when circle.yC==0 the x,y is unique
-      value_t dfx = tR2r0 * gpu::CAMath::Abs(circle.yC) * det;
-      value_t dfy = tR2r0 * circle.xC * (circle.yC > 0.f ? det : -det);
+      auto dfx = tR2r0 * gpu::CAMath::Abs(circle.yC) * det;
+      auto dfy = tR2r0 * circle.xC * (circle.yC > 0. ? det : -det);
       if (dir == DirAuto) {                              // chose the one which corresponds to smallest step
-        value_t delta = (x - mX) * dfx - (y - fy) * dfy; // the choice of + in C will lead to smaller step if delta<0
-        x += delta < 0.f ? dfx : -dfx;
+        auto delta = (x - mX) * dfx - (y - fy) * dfy;    // the choice of + in C will lead to smaller step if delta<0
+        x += delta < 0. ? dfx : -dfx;
       } else if (dir == DirOutward) { // along track direction: x must be > mX
         x -= dfx;                     // try the smallest step (dfx is positive)
-        value_t dfeps = mX - x;       // handle special case of very small step
+        auto dfeps = mX - x;          // handle special case of very small step
         if (dfeps < -kEps) {
           return true;
         }
         if (gpu::CAMath::Abs(dfeps) < kEps && gpu::CAMath::Abs(mX * mX + fy * fy - r * r) < kEps) { // are we already in right r?
-          return mX;
+          x = mX;
+          return true;
         }
         x += dfx + dfx;
-        value_t dxm = x - mX;
-        if (dxm > 0.f) {
+        auto dxm = x - mX;
+        if (dxm > 0.) {
           return true;
         } else if (dxm < -kEps) {
           return false;
         }
         x = mX;                 // don't move
       } else {                  // backward: x must be < mX
         x += dfx;               // try the smallest step (dfx is positive)
-        value_t dfeps = x - mX; // handle special case of very small step
+        auto dfeps = x - mX;    // handle special case of very small step
         if (dfeps < -kEps) {
           return true;
         }
         if (gpu::CAMath::Abs(dfeps) < kEps && gpu::CAMath::Abs(mX * mX + fy * fy - r * r) < kEps) { // are we already in right r?
-          return mX;
+          x = mX;
+          return true;
         }
         x -= dfx + dfx;
-        value_t dxm = x - mX;
-        if (dxm < 0.f) {
+        auto dxm = x - mX;
+        if (dxm < 0.) {
           return true;
         }
         if (dxm > kEps) {
@@ -739,11 +741,11 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
         return false;
       }
     }
-    return x;
+    return true;
   }
   // this is a straight track
   if (gpu::CAMath::Abs(sn) >= constants::math::Almost1) { // || to Y axis
-    value_t det = (r - mX) * (r + mX);
+    double det = (r - mX) * (r + mX);
     if (det < 0.f) {
       return false; // does not reach raduis r
     }
@@ -753,7 +755,7 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
     }
     det = gpu::CAMath::Sqrt(det);
     if (dir == DirOutward) { // along the track direction
-      if (sn > 0.f) {
+      if (sn > 0.) {
         if (fy > det) {
           return false; // track is along Y axis and above the circle
         }
@@ -763,7 +765,7 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
         }
       }
     } else if (dir == DirInward) { // against track direction
-      if (sn > 0.f) {
+      if (sn > 0.) {
         if (fy < -det) {
           return false; // track is along Y axis
         }
@@ -772,13 +774,13 @@ 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
-    value_t det = (r - fy) * (r + fy);
-    if (det < 0.f) {
+    double det = (r - fy) * (r + fy);
+    if (det < 0.) {
       return false; // does not reach raduis r
     }
     det = gpu::CAMath::Sqrt(det);
     if (dir == DirAuto) {
-      x = mX > 0.f ? det : -det; // choose the solution requiring the smalest step
+      x = mX > 0. ? det : -det; // choose the solution requiring the smalest step
       return true;
     } else if (dir == DirOutward) { // along the track direction
       if (mX > det) {
@@ -794,17 +796,17 @@ GPUd() bool TrackParametrization<value_T>::getXatLabR(value_t r, value_t& x, val
       }
     }
   } else { // general case of straight line
-    value_t cs = gpu::CAMath::Sqrt((1.f - sn) * (1.f + sn));
-    value_t xsyc = mX * sn - fy * cs;
-    value_t det = (r - xsyc) * (r + xsyc);
-    if (det < 0.f) {
+    auto cs = gpu::CAMath::Sqrt((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);
-    value_t xcys = mX * cs + fy * sn;
-    value_t t = -xcys;
+    auto xcys = mX * cs + fy * sn;
+    auto t = -xcys;
     if (dir == DirAuto) {
-      t += t > 0.f ? -det : det; // chose the solution requiring the smalest step
+      t += t > 0. ? -det : det;  // chose the solution requiring the smalest step
     } else if (dir > 0) {        // go in increasing mX direction. ( t+-det > 0)
       if (t >= -det) {
         t += det; // take minimal step giving t>0
diff --git a/Detectors/Base/CMakeLists.txt b/Detectors/Base/CMakeLists.txt
@@ -8,6 +8,7 @@
 # In applying this license CERN does not waive the privileges and immunities
 # granted to it by virtue of its status as an Intergovernmental Organization
 # or submit itself to any jurisdiction.
+#add_compile_options(-O0 -g -fPIC)
 
 o2_add_library(DetectorsBase
                SOURCES src/Detector.cxx
diff --git a/Detectors/Base/include/DetectorsBase/Propagator.h b/Detectors/Base/include/DetectorsBase/Propagator.h
@@ -92,13 +92,17 @@ class PropagatorImpl
   GPUd() bool propagateTo(track_T& track, value_type x, bool bzOnly = false, value_type maxSnp = MAX_SIN_PHI, value_type maxStep = MAX_STEP,
                           MatCorrType matCorr = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tofInfo = nullptr, int signCorr = 0) const
   {
-    return bzOnly ? propagateToX(track, x, getNominalBz(), maxSnp, maxStep, matCorr, tofInfo, signCorr) : PropagateToXBxByBz(track, x, maxSnp, maxStep, matCorr, tofInfo, signCorr);
+    return bzOnly ? propagateToX(track, x, getBz(track.getXYZGlo()), maxSnp, maxStep, matCorr, tofInfo, signCorr) : PropagateToXBxByBz(track, x, maxSnp, maxStep, matCorr, tofInfo, signCorr);
   }
 
   template <typename track_T>
   GPUd() bool propagateToAlphaX(track_T& track, value_type alpha, value_type x, bool bzOnly = false, value_type maxSnp = MAX_SIN_PHI, value_type maxStep = MAX_STEP, int minSteps = 1,
                                 MatCorrType matCorr = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tofInfo = nullptr, int signCorr = 0) const;
 
+  template <typename track_T>
+  GPUd() bool propagateToR(track_T& track, value_type r, bool bzOnly = false, value_type maxSnp = MAX_SIN_PHI, value_type maxStep = MAX_STEP,
+                           MatCorrType matCorr = MatCorrType::USEMatCorrLUT, track::TrackLTIntegral* tofInfo = nullptr, int signCorr = 0) const;
+
   GPUd() bool propagateToDCA(const o2::dataformats::VertexBase& vtx, o2::track::TrackParametrizationWithError<value_type>& track, value_type bZ,
                              value_type maxStep = MAX_STEP, MatCorrType matCorr = MatCorrType::USEMatCorrLUT,
                              o2::dataformats::DCA* dcaInfo = nullptr, track::TrackLTIntegral* tofInfo = nullptr,
@@ -157,6 +161,10 @@ class PropagatorImpl
 
   GPUd() void getFieldXYZ(const math_utils::Point3D<double> xyz, double* bxyz) const;
 
+  GPUd() float getBz(const math_utils::Point3D<float> xyz) const;
+
+  GPUd() double getBz(const math_utils::Point3D<double> xyz) const;
+
  private:
 #ifndef GPUCA_GPUCODE
   PropagatorImpl(bool uninitialized = false);
@@ -165,6 +173,8 @@ class PropagatorImpl
   static constexpr value_type Epsilon = 0.00001; // precision of propagation to X
   template <typename T>
   GPUd() void getFieldXYZImpl(const math_utils::Point3D<T> xyz, T* bxyz) const;
+  template <typename T>
+  GPUd() T getBzImpl(const math_utils::Point3D<T> xyz) const;
 
   const o2::field::MagFieldFast* mFieldFast = nullptr; ///< External fast field map (barrel only for the moment)
   o2::field::MagneticField* mField = nullptr;          ///< External nominal field map
diff --git a/Detectors/Base/src/Propagator.cxx b/Detectors/Base/src/Propagator.cxx
diff --git a/Detectors/Vertexing/include/DetectorsVertexing/FwdDCAFitterN.h b/Detectors/Vertexing/include/DetectorsVertexing/FwdDCAFitterN.h
diff --git a/Detectors/Vertexing/include/DetectorsVertexing/HelixHelper.h b/Detectors/Vertexing/include/DetectorsVertexing/HelixHelper.h

Original file line number	Diff line number	Diff line change
`@@ -145,7 +145,7 @@ bool MagFieldFast::GetBcomp(EDim comp, const double xyz[3], double& b) const`
`145`	`145`	`}`
`146`	`146`
`147`	`147`	`//_______________________________________________________________________`
`148`		`-bool MagFieldFast::GetBcomp(EDim comp, const math_utils::Point3D<float> xyz, double& b) const`
	`148`	`+bool MagFieldFast::GetBcomp(EDim comp, const math_utils::Point3D<double> xyz, double& b) const`
`149`	`149`	`{`
`150`	`150`	`// get field`
`151`	`151`	`int zSeg, rSeg, quadrant;`
Original file line number	Diff line number	Diff line change
`@@ -73,6 +73,7 @@ o2_target_root_dictionary(`
`73`	`73`	`include/ReconstructionDataFormats/BCRange.h`
`74`	`74`	`include/ReconstructionDataFormats/TrackHMP.h`
`75`	`75`	`include/ReconstructionDataFormats/MatchInfoHMP.h`
	`76`	`+ include/ReconstructionDataFormats/HelixHelper.h`
`76`	`77`	`)`
`77`	`78`
`78`	`79`	`o2_add_test(Vertex`