refactor: Propagator deduce parameter type from stepper (#2413)

The main change is that the `Acts::Propagator` now deduces the parameters type (bound and curvilinear) from the stepper.

This allows the `MultiStepper` to return multi-component-parameters, and makes the whole structure of the GSF more straight.

Related changes to GSF infrastructure:
* Add multi component curvilinear track parameters
* Remove workaround to get final parameters (this leads to changed output, because with the workaround there was one additional transport-to-bound, which causes floating-point-differences)

Core/include/Acts/EventData/MultiComponentBoundTrackParameters.hpp → Core/include/Acts/EventData/MultiComponentTrackParameters.hpp

@@ -10,6 +10,7 @@
 
 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
 #include "Acts/EventData/TrackParameters.hpp"
+#include "Acts/Surfaces/PlaneSurface.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 
 #include <cmath>
@@ -30,9 +31,14 @@ namespace Acts {
 /// covariance does return std::nullopt;
 /// TODO Add constructor from range and projector maybe?
 class MultiComponentBoundTrackParameters {
+ public:
   using Parameters = BoundTrackParameters;
   using ParticleHypothesis = Parameters::ParticleHypothesis;
+  using Scalar = typename Parameters::Scalar;
+  using ParametersVector = typename Parameters::ParametersVector;
+  using CovarianceMatrix = typename Parameters::CovarianceMatrix;
 
+ private:
   std::vector<std::tuple<double, BoundVector, std::optional<BoundSquareMatrix>>>
       m_components;
   std::shared_ptr<const Surface> m_surface;
@@ -62,10 +68,6 @@ class MultiComponentBoundTrackParameters {
   }
 
  public:
-  using Scalar = typename Parameters::Scalar;
-  using ParametersVector = typename Parameters::ParametersVector;
-  using CovarianceMatrix = typename Parameters::CovarianceMatrix;
-
   /// Construct from multiple components
   template <typename covariance_t>
   MultiComponentBoundTrackParameters(
@@ -224,4 +226,78 @@ class MultiComponentBoundTrackParameters {
   }
 };
 
+/// This class mimics the behaviour of the curvilinear parameters for ordinary
+/// track parameters. To adopt this concept, a "common surface" is constructed,
+/// and all parameters are projected onto this surface. The use of this is
+/// questionable, and if the result is reasonable depends largely on the initial
+/// multi component state. However, the propagator infrastructure forces the
+/// existence of this type
+/// @tparam charge_t Helper type to interpret the particle charge/momentum
+class MultiComponentCurvilinearTrackParameters
+    : public MultiComponentBoundTrackParameters {
+  using covariance_t = BoundSquareMatrix;
+
+ public:
+  using ConstructionTuple = std::tuple<double, Acts::Vector4, Acts::Vector3,
+                                       ActsScalar, covariance_t>;
+
+ private:
+  using Base = MultiComponentBoundTrackParameters;
+
+  using BaseConstructionTuple =
+      std::tuple<std::shared_ptr<Acts::Surface>,
+                 std::vector<std::tuple<double, BoundVector, covariance_t>>>;
+
+  /// We need this helper function in order to construct the base class properly
+  static BaseConstructionTuple construct(
+      const std::vector<ConstructionTuple>& curvi) {
+    // TODO where to get a geometry context here
+    Acts::GeometryContext gctx{};
+
+    // Construct and average surface
+    Acts::Vector3 avgPos = Acts::Vector3::Zero();
+    Acts::Vector3 avgDir = Acts::Vector3::Zero();
+    for (const auto& [w, pos4, dir, qop, cov] : curvi) {
+      avgPos += w * pos4.template segment<3>(0);
+      avgDir += w * dir;
+    }
+
+    auto s = Surface::makeShared<PlaneSurface>(avgPos, avgDir);
+
+    std::vector<std::tuple<double, BoundVector, covariance_t>> bound;
+    bound.reserve(curvi.size());
+
+    // Project the position onto the surface, keep everything else as is
+    for (const auto& [w, pos4, dir, qop, cov] : curvi) {
+      Vector3 newPos =
+          s->intersect(gctx, pos4.template segment<3>(eFreePos0), dir, false)
+              .closest()
+              .position();
+
+      BoundVector bv =
+          detail::transformFreeToCurvilinearParameters(pos4[eTime], dir, qop);
+
+      // Because of the projection this should never fail
+      bv.template segment<2>(eBoundLoc0) =
+          *(s->globalToLocal(gctx, newPos, dir));
+
+      bound.emplace_back(w, bv, cov);
+    }
+
+    return {s, bound};
+  }
+
+  /// Private constructor from a tuple
+  MultiComponentCurvilinearTrackParameters(
+      const BaseConstructionTuple& t, ParticleHypothesis particleHypothesis)
+      : Base(std::get<0>(t), std::get<1>(t), particleHypothesis) {}
+
+ public:
+  MultiComponentCurvilinearTrackParameters(
+      const std::vector<ConstructionTuple>& cmps,
+      ParticleHypothesis particleHypothesis)
+      : MultiComponentCurvilinearTrackParameters(construct(cmps),
+                                                 particleHypothesis) {}
+};
+
 }  // namespace Acts

Core/include/Acts/Propagator/MultiEigenStepperLoop.hpp

@@ -15,7 +15,7 @@
 #include "Acts/Definitions/Direction.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/Definitions/Units.hpp"
-#include "Acts/EventData/MultiComponentBoundTrackParameters.hpp"
+#include "Acts/EventData/MultiComponentTrackParameters.hpp"
 #include "Acts/EventData/TrackParameters.hpp"
 #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
@@ -223,10 +223,6 @@ class MultiEigenStepperLoop
   /// surface
   std::size_t m_stepLimitAfterFirstComponentOnSurface = 50;
 
-  /// How to extract a single component state when calling .boundState() or
-  /// .curvilinearState()
-  MixtureReductionMethod m_finalReductionMethod = MixtureReductionMethod::eMean;
-
   /// The logger (used if no logger is provided by caller of methods)
   std::unique_ptr<const Acts::Logger> m_logger;
 
@@ -243,11 +239,17 @@ class MultiEigenStepperLoop
   using SingleState = typename SingleStepper::State;
 
   /// @brief Use the definitions from the Single-stepper
-  using typename SingleStepper::BoundState;
   using typename SingleStepper::Covariance;
-  using typename SingleStepper::CurvilinearState;
   using typename SingleStepper::Jacobian;
 
+  /// @brief Define an own bound state
+  using BoundState =
+      std::tuple<MultiComponentBoundTrackParameters, Jacobian, ActsScalar>;
+
+  /// @brief Define an own curvilinear state
+  using CurvilinearState = std::tuple<MultiComponentCurvilinearTrackParameters,
+                                      Jacobian, ActsScalar>;
+
   /// @brief The reducer type
   using Reducer = component_reducer_t;
 
@@ -256,12 +258,9 @@ class MultiEigenStepperLoop
 
   /// Constructor from a magnetic field and a optionally provided Logger
   MultiEigenStepperLoop(std::shared_ptr<const MagneticFieldProvider> bField,
-                        MixtureReductionMethod finalReductionMethod =
-                            MixtureReductionMethod::eMean,
                         std::unique_ptr<const Logger> logger =
                             getDefaultLogger("GSF", Logging::INFO))
       : EigenStepper<extensionlist_t, auctioneer_t>(std::move(bField)),
-        m_finalReductionMethod(finalReductionMethod),
         m_logger(std::move(logger)) {}
 
   struct State {
@@ -469,7 +468,7 @@ class MultiEigenStepperLoop
           cmp.state, state.navigation, state.options, state.geoContext);
     }
 
-    Result<BoundState> boundState(
+    Result<typename SingleStepper::BoundState> boundState(
         const Surface& surface, bool transportCov,
         const FreeToBoundCorrection& freeToBoundCorrection) {
       return detail::boundState(

Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp

@@ -17,12 +17,8 @@ auto MultiEigenStepperLoop<E, R, A>::boundState(
     -> Result<BoundState> {
   assert(!state.components.empty());
 
-  if (numberComponents(state) == 1) {
-    return SingleStepper::boundState(state.components.front().state, surface,
-                                     transportCov, freeToBoundCorrection);
-  }
-
-  SmallVector<std::tuple<double, BoundVector, BoundSquareMatrix>> states;
+  std::vector<std::tuple<double, BoundVector, Covariance>> cmps;
+  cmps.reserve(numberComponents(state));
   double accumulatedPathLength = 0.0;
 
   for (auto i = 0ul; i < numberComponents(state); ++i) {
@@ -47,28 +43,20 @@ auto MultiEigenStepperLoop<E, R, A>::boundState(
 
     if (bs.ok()) {
       const auto& btp = std::get<BoundTrackParameters>(*bs);
-      states.emplace_back(
+      cmps.emplace_back(
           state.components[i].weight, btp.parameters(),
           btp.covariance().value_or(Acts::BoundSquareMatrix::Zero()));
       accumulatedPathLength +=
           std::get<double>(*bs) * state.components[i].weight;
     }
   }
 
-  if (states.empty()) {
+  if (cmps.empty()) {
     return MultiStepperError::AllComponentsConversionToBoundFailed;
   }
 
-  const auto [finalPars, cov] =
-      Acts::reduceGaussianMixture(states, surface, m_finalReductionMethod);
-
-  std::optional<BoundSquareMatrix> finalCov = std::nullopt;
-  if (cov != BoundSquareMatrix::Zero()) {
-    finalCov = cov;
-  }
-
-  return BoundState{BoundTrackParameters(surface.getSharedPtr(), finalPars,
-                                         finalCov, particleHypothesis(state)),
+  return BoundState{MultiComponentBoundTrackParameters(
+                        surface.getSharedPtr(), cmps, state.particleHypothesis),
                     Jacobian::Zero(), accumulatedPathLength};
 }
 
@@ -78,53 +66,26 @@ auto MultiEigenStepperLoop<E, R, A>::curvilinearState(State& state,
     -> CurvilinearState {
   assert(!state.components.empty());
 
-  if (numberComponents(state) == 1) {
-    return SingleStepper::curvilinearState(state.components.front().state,
-                                           transportCov);
-  } else if (m_finalReductionMethod == MixtureReductionMethod::eMaxWeight) {
-    auto cmpIt = std::max_element(
-        state.components.begin(), state.components.end(),
-        [](const auto& a, const auto& b) { return a.weight < b.weight; });
-
-    return SingleStepper::curvilinearState(cmpIt->state, transportCov);
-  } else {
-    Vector4 pos4 = Vector4::Zero();
-    Vector3 dir = Vector3::Zero();
-    ActsScalar qop = 0.0;
-    BoundSquareMatrix cov = BoundSquareMatrix::Zero();
-    ActsScalar pathLenth = 0.0;
-    ActsScalar sumOfWeights = 0.0;
-
-    for (auto i = 0ul; i < numberComponents(state); ++i) {
-      const auto [cp, jac, pl] = SingleStepper::curvilinearState(
-          state.components[i].state, transportCov);
-
-      pos4 += state.components[i].weight * cp.fourPosition(state.geoContext);
-      dir += state.components[i].weight * cp.direction();
-      qop += state.components[i].weight * (cp.charge() / cp.absoluteMomentum());
-      if (cp.covariance()) {
-        cov += state.components[i].weight * *cp.covariance();
-      }
-      pathLenth += state.components[i].weight * pathLenth;
-      sumOfWeights += state.components[i].weight;
-    }
-
-    pos4 /= sumOfWeights;
-    dir /= sumOfWeights;
-    qop /= sumOfWeights;
-    pathLenth /= sumOfWeights;
-    cov /= sumOfWeights;
-
-    std::optional<BoundSquareMatrix> finalCov = std::nullopt;
-    if (cov != BoundSquareMatrix::Zero()) {
-      finalCov = cov;
-    }
+  std::vector<
+      std::tuple<double, Vector4, Vector3, ActsScalar, BoundSquareMatrix>>
+      cmps;
+  cmps.reserve(numberComponents(state));
+  double accumulatedPathLength = 0.0;
 
-    return CurvilinearState{
-        CurvilinearTrackParameters(pos4, dir, qop, finalCov,
-                                   particleHypothesis(state)),
-        Jacobian::Zero(), pathLenth};
+  for (auto i = 0ul; i < numberComponents(state); ++i) {
+    const auto [cp, jac, pl] = SingleStepper::curvilinearState(
+        state.components[i].state, transportCov);
+
+    cmps.emplace_back(state.components[i].weight,
+                      cp.fourPosition(state.geoContext), cp.direction(),
+                      (cp.charge() / cp.absoluteMomentum()),
+                      cp.covariance().value_or(BoundSquareMatrix::Zero()));
+    accumulatedPathLength += state.components[i].weight * pl;
   }
+
+  return CurvilinearState{
+      MultiComponentCurvilinearTrackParameters(cmps, state.particleHypothesis),
+      Jacobian::Zero(), accumulatedPathLength};
 }
 
 template <typename E, typename R, typename A>

Core/include/Acts/Propagator/Propagator.hpp

@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2019 CERN for the benefit of the Acts project
+// Copyright (C) 2016-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
@@ -16,12 +16,14 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/PdgParticle.hpp"
 #include "Acts/Definitions/Units.hpp"
+#include "Acts/EventData/TrackParametersConcept.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/MagneticField/MagneticFieldContext.hpp"
 #include "Acts/Propagator/AbortList.hpp"
 #include "Acts/Propagator/ActionList.hpp"
 #include "Acts/Propagator/StandardAborters.hpp"
 #include "Acts/Propagator/StepperConcept.hpp"
+#include "Acts/Propagator/detail/ParameterTraits.hpp"
 #include "Acts/Propagator/detail/VoidPropagatorComponents.hpp"
 #include "Acts/Utilities/Logger.hpp"
 #include "Acts/Utilities/Result.hpp"
@@ -199,10 +201,26 @@ struct PropagatorOptions : public PropagatorPlainOptions {
 ///
 template <typename stepper_t, typename navigator_t = detail::VoidNavigator>
 class Propagator final {
+  /// Re-define bound track parameters dependent on the stepper
+  using StepperBoundTrackParameters =
+      detail::stepper_bound_parameters_type_t<stepper_t>;
+  static_assert(
+      Concepts::BoundTrackParametersConcept<StepperBoundTrackParameters>,
+      "Stepper bound track parameters do not fulfill bound "
+      "parameters concept.");
+
+  /// Re-define curvilinear track parameters dependent on the stepper
+  using StepperCurvilinearTrackParameters =
+      detail::stepper_curvilinear_parameters_type_t<stepper_t>;
+  static_assert(
+      Concepts::BoundTrackParametersConcept<StepperCurvilinearTrackParameters>,
+      "Stepper bound track parameters do not fulfill bound "
+      "parameters concept.");
+
   using Jacobian = BoundMatrix;
-  using BoundState = std::tuple<BoundTrackParameters, Jacobian, double>;
+  using BoundState = std::tuple<StepperBoundTrackParameters, Jacobian, double>;
   using CurvilinearState =
-      std::tuple<CurvilinearTrackParameters, Jacobian, double>;
+      std::tuple<StepperCurvilinearTrackParameters, Jacobian, double>;
 
   static_assert(StepperStateConcept<typename stepper_t::State>,
                 "Stepper does not fulfill stepper concept.");
@@ -349,7 +367,7 @@ class Propagator final {
   template <typename parameters_t, typename propagator_options_t,
             typename path_aborter_t = PathLimitReached>
   Result<
-      action_list_t_result_t<CurvilinearTrackParameters,
+      action_list_t_result_t<StepperCurvilinearTrackParameters,
                              typename propagator_options_t::action_list_type>>
   propagate(const parameters_t& start, const propagator_options_t& options,
             bool makeCurvilinear = true) const;
@@ -376,12 +394,12 @@ class Propagator final {
   template <typename parameters_t, typename propagator_options_t,
             typename path_aborter_t = PathLimitReached>
   Result<
-      action_list_t_result_t<CurvilinearTrackParameters,
+      action_list_t_result_t<StepperCurvilinearTrackParameters,
                              typename propagator_options_t::action_list_type>>
   propagate(
       const parameters_t& start, const propagator_options_t& options,
       bool makeCurvilinear,
-      action_list_t_result_t<CurvilinearTrackParameters,
+      action_list_t_result_t<StepperCurvilinearTrackParameters,
                              typename propagator_options_t::action_list_type>&&
           inputResult) const;
 
@@ -406,8 +424,9 @@ class Propagator final {
   template <typename parameters_t, typename propagator_options_t,
             typename target_aborter_t = SurfaceReached,
             typename path_aborter_t = PathLimitReached>
-  Result<action_list_t_result_t<
-      BoundTrackParameters, typename propagator_options_t::action_list_type>>
+  Result<
+      action_list_t_result_t<StepperBoundTrackParameters,
+                             typename propagator_options_t::action_list_type>>
   propagate(const parameters_t& start, const Surface& target,
             const propagator_options_t& options) const;
 
@@ -433,12 +452,13 @@ class Propagator final {
   template <typename parameters_t, typename propagator_options_t,
             typename target_aborter_t = SurfaceReached,
             typename path_aborter_t = PathLimitReached>
-  Result<action_list_t_result_t<
-      BoundTrackParameters, typename propagator_options_t::action_list_type>>
+  Result<
+      action_list_t_result_t<StepperBoundTrackParameters,
+                             typename propagator_options_t::action_list_type>>
   propagate(
       const parameters_t& start, const Surface& target,
       const propagator_options_t& options,
-      action_list_t_result_t<BoundTrackParameters,
+      action_list_t_result_t<StepperBoundTrackParameters,
                              typename propagator_options_t::action_list_type>
           inputResult) const;