Merge branch 'main' into refactor/calibrator-context

Core/include/Acts/Propagator/StandardAborters.hpp

@@ -206,30 +206,4 @@ struct EndOfWorldReached {
   }
 };
 
-/// If the particle stopped (p=0) abort the propagation
-struct ParticleStopped {
-  ParticleStopped() = default;
-
-  /// boolean operator for abort condition without using the result
-  ///
-  /// @tparam propagator_state_t Type of the propagator state
-  /// @tparam stepper_t Type of the stepper
-  /// @tparam navigator_t Type of the navigator
-  ///
-  /// @param [in,out] state The propagation state object
-  /// @param [in] stepper The stepper object
-  /// @param [in] navigator The navigator object
-  template <typename propagator_state_t, typename stepper_t,
-            typename navigator_t>
-  bool operator()(propagator_state_t& state, const stepper_t& stepper,
-                  const navigator_t& navigator,
-                  const Logger& /*logger*/) const {
-    if (stepper.absoluteMomentum(state.stepping) > 0) {
-      return false;
-    }
-    navigator.targetReached(state.navigation, true);
-    return true;
-  }
-};
-
 }  // namespace Acts

Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp

@@ -255,7 +255,7 @@ struct CombinatorialKalmanFilterResult {
   Result<void> result{Result<void>::success()};
 
   // TODO place into options and make them accessible?
-  AbortList<PathLimitReached, EndOfWorldReached, ParticleStopped> abortList;
+  AbortList<PathLimitReached, EndOfWorldReached> abortList;
 };
 
 /// Combinatorial Kalman filter to find tracks.

Core/include/Acts/Vertexing/HelicalTrackLinearizer.hpp

@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019 CERN for the benefit of the Acts project
+// Copyright (C) 2019-2023 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -25,7 +25,7 @@ namespace Acts {
 /// Linearizes the track parameters at the PCA to a user-provided
 /// point (linPoint). The track parameters are written as a function
 /// of the global PCA position and the momentum of the particle at
-/// the PCA. The linearization then reads (see Eq. 5.7 in Ref(1)):
+/// the PCA. The linearization then reads (see Eq. 5.7 in Ref. (1)):
 ///
 /// q = A (r - r_0) + B (p - p_0) + c,
 ///
@@ -37,7 +37,7 @@ namespace Acts {
 ///
 /// This class computes A and B using the analytic formulae of Ref. (1).
 ///
-/// Ref.(1) - CERN-THESIS-2010-027, Giacinto Piacquadio (Freiburg U.)
+/// Ref. (1) - CERN-THESIS-2010-027, Giacinto Piacquadio (Freiburg U.)
 ///
 /// @tparam propagator_t Propagator type
 /// @tparam propagator_options_t Propagator options type
@@ -79,10 +79,12 @@ class HelicalTrackLinearizer {
     // The propagator
     std::shared_ptr<const Propagator_t> propagator;
 
+    /// Tolerance determining how close we need to get to the Perigee surface to
+    /// reach it during propagation
+    ActsScalar targetTolerance = 1e-12;
+
     // Minimum q/p value
     double minQoP = 1e-15;
-    // Maximum curvature value
-    double maxRho = 1e+15;
   };
 
   /// @brief Constructor

Core/include/Acts/Vertexing/HelicalTrackLinearizer.ipp

@@ -1,12 +1,13 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2019 CERN for the benefit of the Acts project
+// Copyright (C) 2019-2023 CERN for the benefit of the Acts project
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #include "Acts/Surfaces/PerigeeSurface.hpp"
+#include "Acts/Vertexing/LinearizerTrackParameters.hpp"
 
 template <typename propagator_t, typename propagator_options_t>
 Acts::Result<Acts::LinearizedTrack> Acts::
@@ -20,16 +21,20 @@ Acts::Result<Acts::LinearizedTrack> Acts::
   const std::shared_ptr<PerigeeSurface> perigeeSurface =
       Surface::makeShared<PerigeeSurface>(linPointPos);
 
+  // Create propagator options
+  propagator_options_t pOptions(gctx, mctx);
+
+  // Length scale at which we consider to be sufficiently close to the Perigee
+  // surface to skip the propagation.
+  pOptions.targetTolerance = m_cfg.targetTolerance;
+
   // Get intersection of the track with the Perigee if the particle would
   // move on a straight line.
   // This allows us to determine whether we need to propagate the track
   // forward or backward to arrive at the PCA.
   auto intersection = perigeeSurface->intersect(gctx, params.position(gctx),
                                                 params.direction(), false);
 
-  // Create propagator options
-  propagator_options_t pOptions(gctx, mctx);
-
   // Setting the propagation direction using the intersection length from
   // above
   // We handle zero path length as forward propagation, but we could actually
@@ -61,132 +66,139 @@ Acts::Result<Acts::LinearizedTrack> Acts::
 
   // Extracting Perigee parameters and compute functions of them for later
   // usage
-  double phi = paramsAtPCA(BoundIndices::eBoundPhi);
-  double sinPhi = std::sin(phi);
-  double cosPhi = std::cos(phi);
+  ActsScalar d0 = paramsAtPCA(BoundIndices::eBoundLoc0);
 
-  double theta = paramsAtPCA(BoundIndices::eBoundTheta);
-  const double sinTheta = std::sin(theta);
-  const double tanTheta = std::tan(theta);
+  ActsScalar phi = paramsAtPCA(BoundIndices::eBoundPhi);
+  ActsScalar sinPhi = std::sin(phi);
+  ActsScalar cosPhi = std::cos(phi);
 
-  // q over p
-  double qOvP = paramsAtPCA(BoundIndices::eBoundQOverP);
+  ActsScalar theta = paramsAtPCA(BoundIndices::eBoundTheta);
+  ActsScalar sinTheta = std::sin(theta);
+  ActsScalar tanTheta = std::tan(theta);
 
-  // Charge of the particle, determines the sign of the helix radius rho
-  double sgnH = (qOvP < 0.) ? -1 : 1;
+  // q over p
+  ActsScalar qOvP = paramsAtPCA(BoundIndices::eBoundQOverP);
+
+  // Get mass hypothesis from propagator options
+  ActsScalar m0 = pOptions.mass;
+  // Assume unit charge
+  // TODO: Get charge hypothesis from propagator options once they are included
+  // there.
+  ActsScalar p = std::abs(1. / qOvP);
+  // Speed in units of c
+  ActsScalar beta = p / std::hypot(p, m0);
+  // Transverse speed (i.e., speed in the x-y plane)
+  ActsScalar betaT = beta * sinTheta;
 
   // Momentu at the PCA
   Vector3 momentumAtPCA(phi, theta, qOvP);
 
-  // Define Jacobians, which will be filled later
-  ActsMatrix<eBoundSize, 4> positionJacobian;
-  positionJacobian.setZero();
-  ActsMatrix<eBoundSize, 3> momentumJacobian;
-  momentumJacobian.setZero();
+  // Complete Jacobian (consists of positionJacobian and momentumJacobian)
+  ActsMatrix<eBoundSize, eLinSize> completeJacobian =
+      ActsMatrix<eBoundSize, eLinSize>::Zero(eBoundSize, eLinSize);
+
   // get the z-component of the B-field at the PCA
   auto field =
       m_cfg.bField->getField(VectorHelpers::position(pca), state.fieldCache);
   if (!field.ok()) {
     return field.error();
   }
-  double Bz = (*field)[eZ];
+  ActsScalar Bz = (*field)[eZ];
 
   // If there is no magnetic field the particle has a straight trajectory
   // If there is a constant magnetic field it has a helical trajectory
   if (Bz == 0. || std::abs(qOvP) < m_cfg.minQoP) {
-    // Fill position Jacobian, i.e., matrix A from Eq. 5.39 in Ref(1)
-    // First row
-    positionJacobian(0, 0) = -sinPhi;
-    positionJacobian(0, 1) = cosPhi;
-
-    // Second row
-    positionJacobian(1, 0) = -cosPhi / tanTheta;
-    positionJacobian(1, 1) = -sinPhi / tanTheta;
-    positionJacobian(1, 2) = 1.;
-
-    // TODO: include timing in track linearization - will be added
-    // in next PR
-    // Sixth row
-    positionJacobian(5, 3) = 1.;
-
-    // Quantities from Eq. 5.41 and 5.42 in Ref (1)
-    double R = (pca(0) - linPointPos.x()) * cosPhi +
-               (pca(1) - linPointPos.y()) * sinPhi;
-    double Q = (pca(0) - linPointPos.x()) * sinPhi -
-               (pca(1) - linPointPos.y()) * cosPhi;
-
-    // Fill momentum Jacobian, i.e., matrix B from Eq. 5.40 in Ref(1)
-    // First row
-    momentumJacobian(0, 0) = -R;
-
-    // Second row
-    momentumJacobian(1, 0) = Q / tanTheta;
-    momentumJacobian(1, 1) = R / (sinTheta * sinTheta);
-
-    // Third row
-    momentumJacobian(2, 0) = 1.;
-
-    // Fourth row
-    momentumJacobian(3, 1) = 1.;
-
-    // Fifth row
-    momentumJacobian(4, 2) = 1.;
-
-    // TODO: Derivatives of time --> Next PR
+    // Derivatives can be found in Eqs. 5.39 and 5.40 of Ref. (1).
+    // Since we propagated to the PCA (point P in Ref. (1)), we evaluate the
+    // Jacobians there. One can show that, in this case, RTilde = 0 and QTilde =
+    // -d0.
+
+    // Derivatives of d0
+    completeJacobian(eBoundLoc0, eLinPos0) = -sinPhi;
+    completeJacobian(eBoundLoc0, eLinPos1) = cosPhi;
+
+    // Derivatives of z0
+    completeJacobian(eBoundLoc1, eLinPos0) = -cosPhi / tanTheta;
+    completeJacobian(eBoundLoc1, eLinPos1) = -sinPhi / tanTheta;
+    completeJacobian(eBoundLoc1, eLinPos2) = 1.;
+    completeJacobian(eBoundLoc1, eLinPhi) = -d0 / tanTheta;
+
+    // Derivatives of phi
+    completeJacobian(eBoundPhi, eLinPhi) = 1.;
+
+    // Derivatives of theta
+    completeJacobian(eBoundTheta, eLinTheta) = 1.;
+
+    // Derivatives of q/p
+    completeJacobian(eBoundQOverP, eLinQOverP) = 1.;
+
+    // Derivatives of time
+    completeJacobian(eBoundTime, eLinPos0) = -cosPhi / betaT;
+    completeJacobian(eBoundTime, eLinPos1) = -sinPhi / betaT;
+    completeJacobian(eBoundTime, eLinTime) = 1.;
+    completeJacobian(eBoundTime, eLinPhi) = -d0 / betaT;
   } else {
     // Helix radius
-    double rho{sinTheta * (1. / qOvP) / Bz};
-
-    // Quantities from Eq. 5.34 in Ref(1) (see .hpp)
-    double X = pca(0) - linPointPos.x() + rho * sinPhi;
-    double Y = pca(1) - linPointPos.y() - rho * cosPhi;
-    const double S2 = (X * X + Y * Y);
+    ActsScalar rho = sinTheta * (1. / qOvP) / Bz;
+    // Sign of helix radius
+    ActsScalar h = (rho < 0.) ? -1 : 1;
+
+    // Quantities from Eq. 5.34 in Ref. (1) (see .hpp)
+    ActsScalar X = pca(0) - linPointPos.x() + rho * sinPhi;
+    ActsScalar Y = pca(1) - linPointPos.y() - rho * cosPhi;
+    ActsScalar S2 = (X * X + Y * Y);
     // S is the 2D distance from the helix center to linPointPos
     // in the x-y plane
-    const double S = std::sqrt(S2);
-
-    // Fill position Jacobian, i.e., matrix A from Eq. 5.36 in Ref(1)
-    // First row
-    positionJacobian(0, 0) = -sgnH * X / S;
-    positionJacobian(0, 1) = -sgnH * Y / S;
-
-    const double S2tanTheta = S2 * tanTheta;
-
-    // Second row
-    positionJacobian(1, 0) = rho * Y / S2tanTheta;
-    positionJacobian(1, 1) = -rho * X / S2tanTheta;
-    positionJacobian(1, 2) = 1.;
-
-    // Third row
-    positionJacobian(2, 0) = -Y / S2;
-    positionJacobian(2, 1) = X / S2;
-
-    // TODO: include timing in track linearization - will be added
-    // in next PR
-    // Sixth row
-    positionJacobian(5, 3) = 1.;
-
-    // Fill momentum Jacobian, i.e., B matrix from Eq. 5.36 in Ref(1).
-    // Since we propagated to the PCA (point P in Ref(1)), the points
-    // P and V coincide and we can choose deltaPhi = 0.
-    // One can show that if deltaPhi = 0 --> R = 0 and Q = sgnH * S.
-    // As a consequence, many terms of the B matrix from Eq. 5.36 vanish.
-    double rhoOverS{sgnH * rho / S};
-
-    // Second row
-    momentumJacobian(1, 0) = (1. - rhoOverS) * rho / tanTheta;
-
-    // Third row
-    momentumJacobian(2, 0) = rhoOverS;
-
-    // Fourth and fifth row
-    momentumJacobian(3, 1) = 1.;
-    momentumJacobian(4, 2) = 1.;
-
-    // TODO: Derivatives of time --> Next PR
+    ActsScalar S = std::sqrt(S2);
+
+    ActsScalar XoverS2 = X / S2;
+    ActsScalar YoverS2 = Y / S2;
+    ActsScalar rhoCotTheta = rho / tanTheta;
+    ActsScalar rhoOverBetaT = rho / betaT;
+    // Absolute value of rho over S
+    ActsScalar absRhoOverS = h * rho / S;
+
+    // Derivatives can be found in Eq. 5.36 in Ref. (1)
+    // Since we propagated to the PCA (point P in Ref. (1)), the points
+    // P and V coincide, and thus deltaPhi = 0.
+    // One can show that if deltaPhi = 0 --> R = 0 and Q = h * S.
+    // As a consequence, many terms of the B matrix vanish.
+
+    // Derivatives of d0
+    completeJacobian(eBoundLoc0, eLinPos0) = -h * X / S;
+    completeJacobian(eBoundLoc0, eLinPos1) = -h * Y / S;
+
+    // Derivatives of z0
+    completeJacobian(eBoundLoc1, eLinPos0) = rhoCotTheta * YoverS2;
+    completeJacobian(eBoundLoc1, eLinPos1) = -rhoCotTheta * XoverS2;
+    completeJacobian(eBoundLoc1, eLinPos2) = 1.;
+    completeJacobian(eBoundLoc1, eLinPhi) = rhoCotTheta * (1. - absRhoOverS);
+
+    // Derivatives of phi
+    completeJacobian(eBoundPhi, eLinPos0) = -YoverS2;
+    completeJacobian(eBoundPhi, eLinPos1) = XoverS2;
+    completeJacobian(eBoundPhi, eLinPhi) = absRhoOverS;
+
+    // Derivatives of theta
+    completeJacobian(eBoundTheta, eLinTheta) = 1.;
+
+    // Derivatives of q/p
+    completeJacobian(eBoundQOverP, eLinQOverP) = 1.;
+
+    // Derivatives of time
+    completeJacobian(eBoundTime, eLinPos0) = rhoOverBetaT * YoverS2;
+    completeJacobian(eBoundTime, eLinPos1) = -rhoOverBetaT * XoverS2;
+    completeJacobian(eBoundTime, eLinTime) = 1.;
+    completeJacobian(eBoundTime, eLinPhi) = rhoOverBetaT * (1. - absRhoOverS);
   }
 
-  // const term in Talyor expansion from Eq. 5.38 in Ref(1)
+  // Extracting positionJacobian and momentumJacobian from the complete Jacobian
+  ActsMatrix<eBoundSize, eLinPosSize> positionJacobian =
+      completeJacobian.block<eBoundSize, eLinPosSize>(0, 0);
+  ActsMatrix<eBoundSize, eLinMomSize> momentumJacobian =
+      completeJacobian.block<eBoundSize, eLinMomSize>(0, eLinPosSize);
+
+  // const term in Taylor expansion from Eq. 5.38 in Ref. (1)
   BoundVector constTerm =
       paramsAtPCA - positionJacobian * pca - momentumJacobian * momentumAtPCA;
 

Core/include/Acts/Vertexing/NumericalTrackLinearizer.ipp

@@ -46,6 +46,9 @@ Acts::NumericalTrackLinearizer<propagator_t, propagator_options_t>::
 
   // Propagate to the PCA of linPointPos
   auto result = m_cfg.propagator->propagate(params, *perigeeSurface, pOptions);
+  if (not result.ok()) {
+    return result.error();
+  }
 
   // Extracting the Perigee representation of the track wrt linPointPos
   auto endParams = *result->endParameters;
@@ -123,6 +126,9 @@ Acts::NumericalTrackLinearizer<propagator_t, propagator_options_t>::
     // Propagate to the new PCA and extract Perigee parameters
     auto newResult = m_cfg.propagator->propagate(wiggledCurvilinearParams,
                                                  *perigeeSurface, pOptions);
+    if (not newResult.ok()) {
+      return newResult.error();
+    }
     newPerigeeParams = (*newResult->endParameters).parameters();
 
     // Computing the numerical derivatives and filling the Jacobian