fix: Change the component reducer for the GSF (to fix navigation issues)

Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp

@@ -48,6 +48,10 @@ using namespace Acts::UnitLiterals;
 
 /// @brief Reducer struct for the Loop MultiEigenStepper which reduces the
 /// multicomponent state to simply by summing the weighted values
+///
+/// @note Usage is not encouraged, since it can lead to navigation failures
+/// as the global position might not be on surface, even if all components
+/// are on surface
 struct WeightedComponentReducerLoop {
   template <typename component_range_t>
   static Vector3 toVector3(const component_range_t& comps,
@@ -140,74 +144,96 @@ struct WeightedComponentReducerLoop {
   }
 };
 
-struct MaxMomentumReducerLoop {
+namespace detail {
+
+struct MaxMomentumComponent {
+  template <typename component_range_t>
+  auto operator()(const component_range_t& cmps) const {
+    return std::max_element(cmps.begin(), cmps.end(),
+                            [&](const auto& a, const auto& b) {
+                              return std::abs(a.state.pars[eFreeQOverP]) >
+                                     std::abs(b.state.pars[eFreeQOverP]);
+                            });
+  }
+};
+
+struct MaxWeightComponent {
   template <typename component_range_t>
-  static const auto& maxAbsoluteMomentumIt(const component_range_t& cmps) {
-    return *std::max_element(cmps.begin(), cmps.end(),
-                             [&](const auto& a, const auto& b) {
-                               return std::abs(a.state.pars[eFreeQOverP]) >
-                                      std::abs(b.state.pars[eFreeQOverP]);
-                             });
+  auto operator()(const component_range_t& cmps) {
+    return std::max_element(
+        cmps.begin(), cmps.end(),
+        [&](const auto& a, const auto& b) { return a.weight < b.weight; });
   }
+};
 
+template <typename component_chooser_t>
+struct SingleComponentReducer {
   template <typename stepper_state_t>
   static Vector3 position(const stepper_state_t& s) {
-    return maxAbsoluteMomentumIt(s.components)
-        .state.pars.template segment<3>(eFreePos0);
+    return component_chooser_t{}(s.components)
+        ->state.pars.template segment<3>(eFreePos0);
   }
 
   template <typename stepper_state_t>
   static Vector3 direction(const stepper_state_t& s) {
-    return maxAbsoluteMomentumIt(s.components)
-        .state.pars.template segment<3>(eFreeDir0);
+    return component_chooser_t{}(s.components)
+        ->state.pars.template segment<3>(eFreeDir0);
   }
 
   template <typename stepper_state_t>
   static ActsScalar qOverP(const stepper_state_t& s) {
-    const auto& cmp = maxAbsoluteMomentumIt(s.components);
-    return cmp.state.pars[eFreeQOverP];
+    const auto cmp = component_chooser_t{}(s.components);
+    return cmp->state.pars[eFreeQOverP];
   }
 
   template <typename stepper_state_t>
   static ActsScalar absoluteMomentum(const stepper_state_t& s) {
-    const auto& cmp = maxAbsoluteMomentumIt(s.components);
-    return std::abs(cmp.state.absCharge / cmp.state.pars[eFreeQOverP]);
+    const auto cmp = component_chooser_t{}(s.components);
+    return s.particleHypothesis.extractMomentum(cmp->state.pars[eFreeQOverP]);
   }
 
   template <typename stepper_state_t>
   static Vector3 momentum(const stepper_state_t& s) {
-    const auto& cmp = maxAbsoluteMomentumIt(s.components);
-    return std::abs(cmp.state.absCharge / cmp.state.pars[eFreeQOverP]) *
-           cmp.state.pars.template segment<3>(eFreeDir0);
+    const auto cmp = component_chooser_t{}(s.components);
+    return s.particleHypothesis.extractMomentum(cmp->state.pars[eFreeQOverP]) *
+           cmp->state.pars.template segment<3>(eFreeDir0);
   }
 
   template <typename stepper_state_t>
   static ActsScalar charge(const stepper_state_t& s) {
-    return maxAbsoluteMomentumIt(s.components).state.absCharge;
+    const auto cmp = component_chooser_t{}(s.components);
+    return s.particleHypothesis.extractCharge(cmp->state.pars[eFreeQOverP]);
   }
 
   template <typename stepper_state_t>
   static ActsScalar time(const stepper_state_t& s) {
-    return maxAbsoluteMomentumIt(s.components).state.pars[eFreeTime];
+    return component_chooser_t{}(s.components)->state.pars[eFreeTime];
   }
 
   template <typename stepper_state_t>
   static FreeVector pars(const stepper_state_t& s) {
-    return maxAbsoluteMomentumIt(s.components).state.pars;
+    return component_chooser_t{}(s.components)->state.pars;
   }
 
   template <typename stepper_state_t>
   static FreeVector cov(const stepper_state_t& s) {
-    return maxAbsoluteMomentumIt(s.components).state.cov;
+    return component_chooser_t{}(s.components)->state.cov;
   }
 };
 
+}  // namespace detail
+
+using MaxMomentumReducerLoop =
+    detail::SingleComponentReducer<detail::MaxMomentumComponent>;
+using MaxWeightReducerLoop =
+    detail::SingleComponentReducer<detail::MaxWeightComponent>;
+
 /// @brief Stepper based on the EigenStepper, but handles Multi-Component Tracks
 /// (e.g., for the GSF). Internally, this only manages a vector of
 /// EigenStepper::States. This simplifies implementation, but has several
 /// drawbacks:
-/// * There are certain redundancies between the global State and the component
-/// states
+/// * There are certain redundancies between the global State and the
+/// component states
 /// * The components do not share a single magnetic-field-cache
 /// @tparam extensionlist_t See EigenStepper for details
 /// @tparam component_reducer_t How to map the multi-component state to a single

Examples/Algorithms/TrackFitting/src/GsfFitterFunction.cpp

@@ -19,6 +19,7 @@
 #include "Acts/Propagator/MultiEigenStepperLoop.hpp"
 #include "Acts/Propagator/Navigator.hpp"
 #include "Acts/Propagator/Propagator.hpp"
+#include "Acts/Propagator/StepperExtensionList.hpp"
 #include "Acts/TrackFitting/GainMatrixUpdater.hpp"
 #include "Acts/TrackFitting/GaussianSumFitter.hpp"
 #include "Acts/TrackFitting/GsfMixtureReduction.hpp"
@@ -57,7 +58,9 @@ using namespace ActsExamples;
 
 namespace {
 
-using MultiStepper = Acts::MultiEigenStepperLoop<>;
+using MultiStepper = Acts::MultiEigenStepperLoop<
+    Acts::StepperExtensionList<Acts::DefaultExtension>,
+    Acts::MaxWeightReducerLoop>;
 using Propagator = Acts::Propagator<MultiStepper, Acts::Navigator>;
 using DirectPropagator = Acts::Propagator<MultiStepper, Acts::DirectNavigator>;
 

Examples/Python/tests/root_file_hashes.txt

@@ -19,8 +19,8 @@ test_propagation__propagation_summary.root: 85d877c3a48a6235d3d5a398b47130d5d17d
 test_material_recording__geant4_material_tracks.root: c022b9362249b29f57a07926b20644e3ab4ab8ebcf03f773fbf46c446fc1a0a1
 test_truth_tracking_gsf[generic]__trackstates_gsf.root: 1466671bee7b7cfffe90459d9aef316dcee104e6c32a2df4f2f11ea9f12ddc14
 test_truth_tracking_gsf[generic]__tracksummary_gsf.root: b698e3d21eacc34fc8b0ce1d3fbe07405a4b8b549e07f0160573e64c3b401f04
-test_truth_tracking_gsf[odd]__trackstates_gsf.root: 194198c9c9657b02192e573ba9969230663f179104e32e180101ffcd677b4745
-test_truth_tracking_gsf[odd]__tracksummary_gsf.root: 80e60784d4d6eabe734f6fc1ffe04aefa0276cf8798430c590a06f40a07a484f
+test_truth_tracking_gsf[odd]__trackstates_gsf.root: b0d05586724ad2b08e40f81535929375ae55a687d5d3436e9ab2ac4661876e9a
+test_truth_tracking_gsf[odd]__tracksummary_gsf.root: df8da514590ed5ff7236adbcd47030132e1b34873cf3ab264dc5c5f1319bdf98
 test_particle_gun__particles.root: 5fe7dda2933ee6b9615b064d192322fe07831133cd998e5ed99a3b992b713a10
 test_material_mapping__material-map_tracks.root: 938b1a855369e9304401cb10d2751df3fd7acf32a2573f2057eb1691cd94edf3
 test_material_mapping__propagation-material.root: 5eacd0cb804d381171c8fb65d7b2d36e1d57db9f4cb2b58c0f24703479218cbf

Tests/UnitTests/Core/Propagator/MultiStepperTests.cpp

@@ -59,7 +59,8 @@ const MagneticFieldContext magCtx;
 const GeometryContext geoCtx;
 
 using MultiStepperLoop =
-    MultiEigenStepperLoop<StepperExtensionList<DefaultExtension>>;
+    MultiEigenStepperLoop<StepperExtensionList<DefaultExtension>,
+                          MaxWeightReducerLoop>;
 using SingleStepper = EigenStepper<StepperExtensionList<DefaultExtension>>;
 
 const double defaultStepSize = 123.;
@@ -134,6 +135,86 @@ auto makeDefaultBoundPars(bool cov = true, std::size_t n = 4,
   return MultiComponentBoundTrackParameters(surface, cmps, particleHypothesis);
 }
 
+//////////////////////
+/// Test the reducers
+//////////////////////
+BOOST_AUTO_TEST_CASE(test_weighted_reducer) {
+  // Can use this multistepper since we only care about the state which is
+  // invariant
+  using MultiState = typename MultiStepperLoop::State;
+
+  constexpr std::size_t N = 4;
+  const auto multi_pars = makeDefaultBoundPars(false, N);
+
+  MultiState state(geoCtx, magCtx, defaultBField, multi_pars, defaultStepSize);
+  SingleStepper singleStepper(defaultBField);
+
+  WeightedComponentReducerLoop reducer{};
+
+  Acts::Vector3 pos = Acts::Vector3::Zero();
+  Acts::Vector3 dir = Acts::Vector3::Zero();
+  for (const auto &[sstate, weight, _] : state.components) {
+    pos += weight * singleStepper.position(sstate);
+    dir += weight * singleStepper.direction(sstate);
+  }
+  dir.normalize();
+
+  BOOST_CHECK_EQUAL(reducer.position(state), pos);
+  BOOST_CHECK_EQUAL(reducer.direction(state), dir);
+}
+
+BOOST_AUTO_TEST_CASE(test_max_weight_reducer) {
+  // Can use this multistepper since we only care about the state which is
+  // invariant
+  using MultiState = typename MultiStepperLoop::State;
+
+  constexpr std::size_t N = 4;
+  const auto multi_pars = makeDefaultBoundPars(false, N);
+  MultiState state(geoCtx, magCtx, defaultBField, multi_pars, defaultStepSize);
+  SingleStepper singleStepper(defaultBField);
+
+  double w = 0.1;
+  double wSum = 0.0;
+  for (auto &[sstate, weight, _] : state.components) {
+    weight = w;
+    wSum += w;
+    w += 0.1;
+  }
+  BOOST_CHECK_EQUAL(wSum, 1.0);
+  BOOST_CHECK_EQUAL(state.components.back().weight, 0.4);
+
+  MaxWeightReducerLoop reducer{};
+  BOOST_CHECK_EQUAL(reducer.position(state),
+                    singleStepper.position(state.components.back().state));
+  BOOST_CHECK_EQUAL(reducer.direction(state),
+                    singleStepper.direction(state.components.back().state));
+}
+
+BOOST_AUTO_TEST_CASE(test_max_momentum_reducer) {
+  // Can use this multistepper since we only care about the state which is
+  // invariant
+  using MultiState = typename MultiStepperLoop::State;
+
+  constexpr std::size_t N = 4;
+  const auto multi_pars = makeDefaultBoundPars(false, N);
+  MultiState state(geoCtx, magCtx, defaultBField, multi_pars, defaultStepSize);
+  SingleStepper singleStepper(defaultBField);
+
+  double p = 1.0;
+  double q = 1.0;
+  for (auto &[sstate, weight, _] : state.components) {
+    sstate.pars[eFreeQOverP] = q / p;
+    p *= 2.0;
+  }
+  BOOST_CHECK_EQUAL(state.components.back().state.pars[eFreeQOverP], q / 8.0);
+
+  MaxMomentumReducerLoop reducer{};
+  BOOST_CHECK_EQUAL(reducer.position(state),
+                    singleStepper.position(state.components.back().state));
+  BOOST_CHECK_EQUAL(reducer.direction(state),
+                    singleStepper.direction(state.components.back().state));
+}
+
 //////////////////////////////////////////////////////
 /// Test the construction of the MultiStepper::State
 //////////////////////////////////////////////////////