fix: Fix stepper charge hypothesis (#2140)

I found at least one case where we rely on the sign of the charge hypothesis in the stepper. The charge might flip during KF but the one inside the stepper would stay the same causing inconsistencies and in my case a terrible fit.

In this PR I remove the sign of the charge hypothesis, refactored the `momentum` and `charge` accessor methods and added a `qop` accessor method that should be preferred over explicit `charge / momentum` IMO

Core/include/Acts/Navigation/NavigationState.hpp

@@ -60,8 +60,8 @@ struct NavigationState {
   /// The current absolute momentum
   ActsScalar absMomentum = 0.;
 
-  /// The current charge
-  ActsScalar charge = 0.;
+  /// The current absolute charge
+  ActsScalar absCharge = 0.;
 
   /// The current magnetic field
   Vector3 magneticField = Vector3(0., 0., 0.);

Core/include/Acts/Navigation/NextNavigator.hpp

@@ -414,7 +414,6 @@ class NextNavigator {
     nState.position = stepper.position(state.stepping);
     nState.direction = stepper.direction(state.stepping);
     nState.absMomentum = stepper.momentum(state.stepping);
-    nState.charge = stepper.charge(state.stepping);
     auto fieldResult = stepper.getField(state.stepping, nState.position);
     if (!fieldResult.ok()) {
       ACTS_ERROR(volInfo(state) << posInfo(state, stepper)

Core/include/Acts/Propagator/AtlasStepper.hpp

@@ -353,14 +353,14 @@ class AtlasStepper {
     return Vector3(state.pVector[4], state.pVector[5], state.pVector[6]);
   }
 
+  double qop(const State& state) const { return state.pVector[7]; }
+
   double momentum(const State& state) const {
-    return 1. / std::abs(state.pVector[7]);
+    return 1. / std::abs(qop(state));
   }
 
   /// Charge access
-  double charge(const State& state) const {
-    return state.pVector[7] > 0. ? 1. : -1.;
-  }
+  double charge(const State& state) const { return qop(state) > 0. ? 1. : -1.; }
 
   /// Overstep limit
   double overstepLimit(const State& /*state*/) const { return m_overstepLimit; }
@@ -705,10 +705,10 @@ class AtlasStepper {
   /// @param state The state object
   /// @param uposition the updated position
   /// @param udirection the updated direction
-  /// @param up the updated momentum value
+  /// @param qop the updated momentum value
   /// @param time the update time
   void update(State& state, const Vector3& uposition, const Vector3& udirection,
-              double up, double time) const {
+              double qop, double time) const {
     // update the vector
     state.pVector[0] = uposition[0];
     state.pVector[1] = uposition[1];
@@ -717,7 +717,7 @@ class AtlasStepper {
     state.pVector[4] = udirection[0];
     state.pVector[5] = udirection[1];
     state.pVector[6] = udirection[2];
-    state.pVector[7] = charge(state) / up;
+    state.pVector[7] = qop;
   }
 
   /// Method for on-demand transport of the covariance

Core/include/Acts/Propagator/EigenStepper.hpp

@@ -78,7 +78,7 @@ class EigenStepper {
                    Direction ndir = Direction::Forward,
                    double ssize = std::numeric_limits<double>::max(),
                    double stolerance = s_onSurfaceTolerance)
-        : q(par.charge()),
+        : absCharge(std::abs(par.charge())),
           navDir(ndir),
           stepSize(ndir * std::abs(ssize)),
           tolerance(stolerance),
@@ -103,8 +103,8 @@ class EigenStepper {
     /// Internal free vector parameters
     FreeVector pars = FreeVector::Zero();
 
-    /// The charge as the free vector can be 1/p or q/p
-    double q = 1.;
+    /// The absolute charge as the free vector can be 1/p or q/p
+    double absCharge = UnitConstants::e;
 
     /// Covariance matrix (and indicator)
     /// associated with the initial error on track parameters
@@ -213,17 +213,25 @@ class EigenStepper {
     return state.pars.template segment<3>(eFreeDir0);
   }
 
+  /// QoP direction accessor
+  ///
+  /// @param state [in] The stepping state (thread-local cache)
+  double qop(const State& state) const { return state.pars[eFreeQOverP]; }
+
   /// Absolute momentum accessor
   ///
   /// @param state [in] The stepping state (thread-local cache)
   double momentum(const State& state) const {
-    return std::abs((state.q == 0. ? 1. : state.q) / state.pars[eFreeQOverP]);
+    auto q = charge(state);
+    return std::abs((q == 0 ? 1 : q) / qop(state));
   }
 
   /// Charge access
   ///
   /// @param state [in] The stepping state (thread-local cache)
-  double charge(const State& state) const { return state.q; }
+  double charge(const State& state) const {
+    return std::copysign(state.absCharge, qop(state));
+  }
 
   /// Time access
   ///
@@ -350,15 +358,15 @@ class EigenStepper {
               const BoundVector& boundParams, const Covariance& covariance,
               const Surface& surface) const;
 
-  /// Method to update momentum, direction and p
+  /// Method to update the stepper state
   ///
   /// @param [in,out] state State object that will be updated
   /// @param [in] uposition the updated position
   /// @param [in] udirection the updated direction
-  /// @param [in] up the updated momentum value
+  /// @param [in] qop the updated qop value
   /// @param [in] time the updated time value
   void update(State& state, const Vector3& uposition, const Vector3& udirection,
-              double up, double time) const;
+              double qop, double time) const;
 
   /// Method for on-demand transport of the covariance
   /// to a new curvilinear frame at current  position,

Core/include/Acts/Propagator/EigenStepper.ipp

@@ -84,12 +84,12 @@ void Acts::EigenStepper<E, A>::update(State& state,
 
 template <typename E, typename A>
 void Acts::EigenStepper<E, A>::update(State& state, const Vector3& uposition,
-                                      const Vector3& udirection, double up,
+                                      const Vector3& udirection, double qop,
                                       double time) const {
   state.pars.template segment<3>(eFreePos0) = uposition;
   state.pars.template segment<3>(eFreeDir0) = udirection;
   state.pars[eFreeTime] = time;
-  state.pars[eFreeQOverP] = (state.q != 0. ? state.q / up : 1. / up);
+  state.pars[eFreeQOverP] = qop;
 }
 
 template <typename E, typename A>

Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp

@@ -63,23 +63,35 @@ struct WeightedComponentReducerLoop {
     return toVector3(s.components, eFreeDir0).normalized();
   }
 
+  // TODO: Maybe we can cache this value and only update it when the parameters
+  // change
+  template <typename stepper_state_t>
+  static ActsScalar qop(const stepper_state_t& s) {
+    return std::accumulate(
+        s.components.begin(), s.components.end(), ActsScalar{0.},
+        [](const auto& sum, const auto& cmp) -> ActsScalar {
+          return sum + cmp.weight * cmp.state.pars[eFreeQOverP];
+        });
+  }
+
   template <typename stepper_state_t>
   static ActsScalar momentum(const stepper_state_t& s) {
     return std::accumulate(
         s.components.begin(), s.components.end(), ActsScalar{0.},
         [](const auto& sum, const auto& cmp) -> ActsScalar {
-          return sum +
-                 cmp.weight * (1 / (cmp.state.pars[eFreeQOverP] / cmp.state.q));
+          return sum + cmp.weight * std::abs(cmp.state.absCharge /
+                                             cmp.state.pars[eFreeQOverP]);
         });
   }
 
   template <typename stepper_state_t>
   static ActsScalar charge(const stepper_state_t& s) {
-    return std::accumulate(s.components.begin(), s.components.end(),
-                           ActsScalar{0.},
-                           [](const auto& sum, const auto& cmp) -> ActsScalar {
-                             return sum + cmp.weight * cmp.state.q;
-                           });
+    return std::accumulate(
+        s.components.begin(), s.components.end(), ActsScalar{0.},
+        [](const auto& sum, const auto& cmp) -> ActsScalar {
+          return sum + cmp.weight * std::copysign(cmp.state.absCharge,
+                                                  cmp.state.pars[eFreeQOverP]);
+        });
   }
 
   template <typename stepper_state_t>
@@ -135,12 +147,12 @@ struct MaxMomentumReducerLoop {
   template <typename stepper_state_t>
   static ActsScalar momentum(const stepper_state_t& s) {
     const auto& cmp = maxMomenutmIt(s.components);
-    return 1.0 / (cmp.state.pars[eFreeQOverP] / cmp.state.q);
+    return std::abs(cmp.state.absCharge / cmp.state.pars[eFreeQOverP]);
   }
 
   template <typename stepper_state_t>
   static ActsScalar charge(const stepper_state_t& s) {
-    return maxMomenutmIt(s.components).state.q;
+    return maxMomenutmIt(s.components).state.absCharge;
   }
 
   template <typename stepper_state_t>
@@ -355,7 +367,6 @@ class MultiEigenStepperLoop
     // ComponentProxy and the ConstComponentProxy
     auto status() const { return cmp.status; }
     auto weight() const { return cmp.weight; }
-    auto charge() const { return cmp.state.q; }
     auto pathAccumulated() const { return cmp.state.pathAccumulated; }
     const auto& pars() const { return cmp.state.pars; }
     const auto& derivative() const { return cmp.state.derivative; }
@@ -395,7 +406,6 @@ class MultiEigenStepperLoop
     using Base = ComponentProxyBase<typename State::Component>;
 
     // Import the const accessors from ComponentProxyBase
-    using Base::charge;
     using Base::cmp;
     using Base::cov;
     using Base::derivative;
@@ -419,7 +429,6 @@ class MultiEigenStepperLoop
     // ComponentProxyBase
     auto& status() { return cmp.status; }
     auto& weight() { return cmp.weight; }
-    auto& charge() { return cmp.state.q; }
     auto& pathAccumulated() { return cmp.state.pathAccumulated; }
     auto& pars() { return cmp.state.pars; }
     auto& derivative() { return cmp.state.derivative; }

Core/include/Acts/Propagator/Propagator.hpp

@@ -63,7 +63,7 @@ struct PropagatorPlainOptions {
   /// pion default
   int absPdgCode = PdgParticle::ePionPlus;
 
-  /// The mass for the particle for (eventual) material integration -
+  /// The mass of the particle for (eventual) material integration -
   /// pion default
   double mass = 139.57018 * UnitConstants::MeV;
 

Core/include/Acts/Propagator/StraightLineStepper.hpp

@@ -66,7 +66,7 @@ class StraightLineStepper {
                    Direction ndir = Direction::Forward,
                    double ssize = std::numeric_limits<double>::max(),
                    double stolerance = s_onSurfaceTolerance)
-        : q(par.charge()),
+        : absCharge(std::abs(par.charge())),
           navDir(ndir),
           stepSize(ndir * std::abs(ssize)),
           tolerance(stolerance),
@@ -101,8 +101,8 @@ class StraightLineStepper {
     /// Internal free vector parameters
     FreeVector pars = FreeVector::Zero();
 
-    /// The charge as the free vector can be 1/p or q/p
-    double q = 1.;
+    /// The absolute charge as the free vector can be 1/p or q/p
+    double absCharge = UnitConstants::e;
 
     /// Boolean to indiciate if you need covariance transport
     bool covTransport = false;
@@ -182,17 +182,25 @@ class StraightLineStepper {
     return state.pars.template segment<3>(eFreeDir0);
   }
 
+  /// QoP direction accessor
+  ///
+  /// @param state [in] The stepping state (thread-local cache)
+  double qop(const State& state) const { return state.pars[eFreeQOverP]; }
+
   /// Absolute momentum accessor
   ///
   /// @param state [in] The stepping state (thread-local cache)
   double momentum(const State& state) const {
-    return std::abs((state.q == 0. ? 1. : state.q) / state.pars[eFreeQOverP]);
+    auto q = charge(state);
+    return std::abs((q == 0 ? 1 : q) / qop(state));
   }
 
   /// Charge access
   ///
   /// @param state [in] The stepping state (thread-local cache)
-  double charge(const State& state) const { return state.q; }
+  double charge(const State& state) const {
+    return std::copysign(state.absCharge, qop(state));
+  }
 
   /// Time access
   ///
@@ -319,15 +327,15 @@ class StraightLineStepper {
               const BoundVector& boundParams, const Covariance& covariance,
               const Surface& surface) const;
 
-  /// Method to update momentum, direction and p
+  /// Method to update the stepper state
   ///
   /// @param [in,out] state State object that will be updated
   /// @param [in] uposition the updated position
   /// @param [in] udirection the updated direction
-  /// @param [in] up the updated momentum value
+  /// @param [in] qop the updated qop value
   /// @param [in] time the updated time value
   void update(State& state, const Vector3& uposition, const Vector3& udirection,
-              double up, double time) const;
+              double qop, double time) const;
 
   /// Method for on-demand transport of the covariance
   /// to a new curvilinear frame at current  position,

Core/include/Acts/Propagator/detail/GenericDefaultExtension.hpp

@@ -62,8 +62,7 @@ struct GenericDefaultExtension {
          const navigator_t& /*navigator*/, ThisVector3& knew,
          const Vector3& bField, std::array<Scalar, 4>& kQoP, const int i = 0,
          const double h = 0., const ThisVector3& kprev = ThisVector3::Zero()) {
-    auto qop =
-        stepper.charge(state.stepping) / stepper.momentum(state.stepping);
+    auto qop = stepper.qop(state.stepping);
     // First step does not rely on previous data
     if (i == 0) {
       knew = qop * stepper.direction(state.stepping).cross(bField);
@@ -183,8 +182,7 @@ struct GenericDefaultExtension {
 
     auto& sd = state.stepping.stepData;
     auto dir = stepper.direction(state.stepping);
-    auto qop =
-        stepper.charge(state.stepping) / stepper.momentum(state.stepping);
+    auto qop = stepper.qop(state.stepping);
 
     D = FreeMatrix::Identity();
 
@@ -245,9 +243,8 @@ struct GenericDefaultExtension {
 
     D(3, 7) =
         h * state.options.mass * state.options.mass *
-        stepper.charge(state.stepping) /
-        (stepper.momentum(state.stepping) *
-         std::hypot(1., state.options.mass / stepper.momentum(state.stepping)));
+        stepper.qop(state.stepping) /
+        std::hypot(1., state.options.mass / stepper.momentum(state.stepping));
     return true;
   }
 };

Core/include/Acts/Propagator/detail/GenericDenseEnvironmentExtension.hpp

@@ -115,6 +115,8 @@ struct GenericDenseEnvironmentExtension {
          const navigator_t& navigator, ThisVector3& knew, const Vector3& bField,
          std::array<Scalar, 4>& kQoP, const int i = 0, const double h = 0.,
          const ThisVector3& kprev = ThisVector3::Zero()) {
+    auto q = stepper.charge(state.stepping);
+
     // i = 0 is used for setup and evaluation of k
     if (i == 0) {
       // Set up container for energy loss
@@ -123,14 +125,12 @@ struct GenericDenseEnvironmentExtension {
       material = (volumeMaterial->material(position.template cast<double>()));
       initialMomentum = stepper.momentum(state.stepping);
       currentMomentum = initialMomentum;
-      qop[0] = stepper.charge(state.stepping) / initialMomentum;
+      qop[0] = q / initialMomentum;
       initializeEnergyLoss(state);
       // Evaluate k
       knew = qop[0] * stepper.direction(state.stepping).cross(bField);
       // Evaluate k for the time propagation
-      Lambdappi[0] =
-          -qop[0] * qop[0] * qop[0] * g * energy[0] /
-          (stepper.charge(state.stepping) * stepper.charge(state.stepping));
+      Lambdappi[0] = -qop[0] * qop[0] * qop[0] * g * energy[0] / (q * q);
       //~ tKi[0] = std::hypot(1, state.options.mass / initialMomentum);
       tKi[0] = hypot(1, state.options.mass * qop[0]);
       kQoP[0] = Lambdappi[0];
@@ -145,10 +145,8 @@ struct GenericDenseEnvironmentExtension {
              (stepper.direction(state.stepping) + h * kprev).cross(bField);
       // Evaluate k_i for the time propagation
       auto qopNew = qop[0] + h * Lambdappi[i - 1];
-      Lambdappi[i] =
-          -qopNew * qopNew * qopNew * g * energy[i] /
-          (stepper.charge(state.stepping) * stepper.charge(state.stepping) *
-           UnitConstants::C * UnitConstants::C);
+      Lambdappi[i] = -qopNew * qopNew * qopNew * g * energy[i] /
+                     (q * q * UnitConstants::C * UnitConstants::C);
       tKi[i] = hypot(1, state.options.mass * qopNew);
       kQoP[i] = Lambdappi[i];
     }