refactor!: Distinct MultiComponentTrackParameters (#5150)

Decouples the `MultiComponentTrackParameters` interface from the usual bound track parameters. This is motivated by their distinct nature. A `MultiComponentTrackParameters` cannot be uniquely represented as a single bound track parameter but requires a user specified merging procedure.

Additionally, refactors the component storage to use SoA over AoS to model optional covariance more efficiently.

Author: andiwand

Core/include/Acts/EventData/MultiComponentTrackParameters.hpp

@@ -13,6 +13,7 @@
 #include "Acts/Geometry/GeometryContext.hpp"
 
 #include <optional>
+#include <tuple>
 
 namespace Acts {
 

Core/include/Acts/EventData/TrackParametersConcept.hpp

@@ -0,0 +1,71 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/TrackParametrization.hpp"
+
+#include <tuple>
+
+namespace Acts::detail {
+
+template <typename cmp_t>
+concept ComponentWithoutCovarianceConcept =
+    requires { typename std::tuple_size<std::remove_cvref_t<cmp_t>>::type; } &&
+    (std::tuple_size_v<std::remove_cvref_t<cmp_t>> == 2) &&
+    requires(const cmp_t& cmp) {
+      { std::get<0>(cmp) } -> std::convertible_to<const double&>;
+      { std::get<1>(cmp) } -> std::convertible_to<const BoundVector&>;
+    };
+
+template <typename cmp_t>
+concept ComponentWithCovarianceConcept =
+    requires { typename std::tuple_size<std::remove_cvref_t<cmp_t>>::type; } &&
+    (std::tuple_size_v<std::remove_cvref_t<cmp_t>> == 3) &&
+    requires(const cmp_t& cmp) {
+      { std::get<0>(cmp) } -> std::convertible_to<const double&>;
+      { std::get<1>(cmp) } -> std::convertible_to<const BoundVector&>;
+      { std::get<2>(cmp) } -> std::convertible_to<const BoundMatrix&>;
+    };
+
+template <typename cmp_t>
+concept ComponentConcept = ComponentWithoutCovarianceConcept<cmp_t> ||
+                           ComponentWithCovarianceConcept<cmp_t>;
+
+template <typename projector_t, typename cmp_t>
+concept ComponentWithoutCovarianceProjectorConcept =
+    requires(const projector_t& proj, const cmp_t& cmp) {
+      { proj(cmp) } -> ComponentWithoutCovarianceConcept;
+    };
+
+template <typename projector_t, typename cmp_t>
+concept ComponentWithCovarianceProjectorConcept =
+    requires(const projector_t& proj, const cmp_t& cmp) {
+      { proj(cmp) } -> ComponentWithCovarianceConcept;
+    };
+
+template <typename projector_t, typename cmp_t>
+concept ComponentProjectorConcept =
+    ComponentWithoutCovarianceProjectorConcept<projector_t, cmp_t> ||
+    ComponentWithCovarianceProjectorConcept<projector_t, cmp_t>;
+
+template <typename component_range_t, typename projector_t>
+concept ComponentRangeAndProjectorWithoutCovarianceConcept =
+    std::ranges::range<component_range_t> &&
+    requires(const component_range_t& cmps, const projector_t& proj) {
+      { proj(*cmps.begin()) } -> ComponentWithoutCovarianceConcept;
+    };
+
+template <typename component_range_t, typename projector_t>
+concept ComponentRangeAndProjectorWithCovarianceConcept =
+    std::ranges::range<component_range_t> &&
+    requires(const component_range_t& cmps, const projector_t& proj) {
+      { proj(*cmps.begin()) } -> ComponentWithCovarianceConcept;
+    };
+
+}  // namespace Acts::detail

Core/include/Acts/EventData/detail/MultiComponentTrackParametersConcept.hpp

@@ -31,7 +31,6 @@
 #include <cstddef>
 #include <limits>
 #include <sstream>
-#include <vector>
 
 #include <boost/container/small_vector.hpp>
 
@@ -274,27 +273,19 @@ class MultiStepperLoop final {
   /// @param par The multi-component bound track parameters
   void initialize(State& state,
                   const MultiComponentBoundTrackParameters& par) const {
-    if (par.components().empty()) {
-      throw std::invalid_argument(
-          "Cannot construct MultiEigenStepperLoop::State with empty "
-          "multi-component parameters");
-    }
-
     state.particleHypothesis = par.particleHypothesis();
 
     const auto surface = par.referenceSurface().getSharedPtr();
 
-    for (auto i = 0ul; i < par.components().size(); ++i) {
+    for (std::size_t i = 0; i < par.size(); ++i) {
       const auto& [weight, singlePars] = par[i];
       auto& cmp = state.components.emplace_back(
           m_singleStepper.makeState(state.options), weight,
           IntersectionStatus::onSurface);
       m_singleStepper.initialize(cmp.state, singlePars);
     }
 
-    if (std::get<2>(par.components().front())) {
-      state.covTransport = true;
-    }
+    state.covTransport = par.hasCovariance();
   }
 
   /// A proxy struct which allows access to a single component of the

Core/include/Acts/Propagator/MultiStepperLoop.hpp

@@ -21,8 +21,10 @@ auto MultiStepperLoop<S, R>::boundState(
     -> Result<BoundState> {
   assert(!state.components.empty());
 
-  std::vector<std::tuple<double, BoundVector, Covariance>> cmps;
-  cmps.reserve(numberComponents(state));
+  MultiComponentBoundTrackParameters params(
+      surface.getSharedPtr(), transportCov, state.particleHypothesis);
+  params.reserve(numberComponents(state));
+
   double accumulatedPathLength = 0.0;
 
   for (auto i = 0ul; i < numberComponents(state); ++i) {
@@ -48,20 +50,20 @@ auto MultiStepperLoop<S, R>::boundState(
 
     if (bs.ok()) {
       const auto& btp = std::get<BoundTrackParameters>(*bs);
-      cmps.emplace_back(state.components[i].weight, btp.parameters(),
-                        btp.covariance().value_or(Acts::BoundMatrix::Zero()));
+      params.pushComponent(state.components[i].weight, btp.parameters(),
+                           btp.covariance());
       accumulatedPathLength +=
           std::get<double>(*bs) * state.components[i].weight;
     }
   }
 
-  if (cmps.empty()) {
+  if (params.empty()) {
     return MultiStepperError::AllComponentsConversionToBoundFailed;
   }
 
-  return BoundState{MultiComponentBoundTrackParameters(
-                        surface.getSharedPtr(), cmps, state.particleHypothesis),
-                    Jacobian::Zero(), accumulatedPathLength};
+  params.normalizeWeights();
+
+  return BoundState{params, Jacobian::Zero(), accumulatedPathLength};
 }
 
 template <Concepts::SingleStepper S, typename R>

Core/include/Acts/Propagator/MultiStepperLoop.ipp

@@ -93,9 +93,6 @@ class Propagator final
   /// Re-define bound track parameters dependent on the stepper
   using StepperBoundTrackParameters =
       detail::stepper_bound_parameters_type_t<stepper_t>;
-  static_assert(BoundTrackParametersConcept<StepperBoundTrackParameters>,
-                "Stepper bound track parameters do not fulfill bound "
-                "parameters concept.");
   static_assert(std::copy_constructible<StepperBoundTrackParameters>,
                 "return track parameter type must be copy-constructible");
 

Core/include/Acts/Propagator/Propagator.hpp

@@ -311,9 +311,6 @@ template <typename propagator_state_t, typename parameters_t,
           typename path_aborter_t>
 Result<void> Propagator<S, N>::initialize(propagator_state_t& state,
                                           const parameters_t& start) const {
-  static_assert(BoundTrackParametersConcept<parameters_t>,
-                "Parameters do not fulfill bound parameters concept.");
-
   m_stepper.initialize(state.stepping, start);
 
   state.position = m_stepper.position(state.stepping);

Core/include/Acts/Propagator/Propagator.ipp

@@ -305,7 +305,7 @@ struct GaussianSumFitter {
       if constexpr (!IsMultiParameters::value) {
         params = MultiComponentBoundTrackParameters(
             sParameters.referenceSurface().getSharedPtr(),
-            sParameters.parameters(), *sParameters.covariance(),
+            sParameters.parameters(), sParameters.covariance(),
             sParameters.particleHypothesis());
       } else {
         params = sParameters;
@@ -377,18 +377,19 @@ struct GaussianSumFitter {
                                   ? *options.referenceSurface
                                   : sParameters.referenceSurface();
 
-      std::vector<std::tuple<double, BoundVector, std::optional<BoundMatrix>>>
-          inflatedParamVector;
       assert(!fwdGsfResult.lastMeasurementComponents.empty());
       assert(fwdGsfResult.lastMeasurementSurface != nullptr);
-      for (auto& [w, p, cov] : fwdGsfResult.lastMeasurementComponents) {
-        inflatedParamVector.emplace_back(
-            w, p, cov * options.reverseFilteringCovarianceScaling);
-      }
 
       MultiComponentBoundTrackParameters inflatedParams(
           fwdGsfResult.lastMeasurementSurface->getSharedPtr(),
-          std::move(inflatedParamVector), sParameters.particleHypothesis());
+          fwdGsfResult.lastMeasurementComponents,
+          [&options](const auto& cmp)
+              -> std::tuple<double, const BoundVector&, BoundMatrix> {
+            return {
+                std::get<0>(cmp), std::get<1>(cmp),
+                std::get<2>(cmp) * options.reverseFilteringCovarianceScaling};
+          },
+          sParameters.particleHypothesis());
 
       auto state = m_propagator.template makeState<decltype(bwdPropOptions),
                                                    MultiStepperSurfaceReached>(
@@ -504,17 +505,10 @@ struct GaussianSumFitter {
 
     if (options.referenceSurface) {
       const auto& params = *bwdResult->endParameters;
+      const auto singleParams = params.merge(options.componentMergeMethod);
 
-      const auto [finalPars, finalCov] = detail::Gsf::mergeGaussianMixture(
-          params.components(),
-          [](const auto& cmp) {
-            auto&& [weight_l, pars_l, opt_cov_l] = cmp;
-            return std::tie(weight_l, pars_l, *opt_cov_l);
-          },
-          params.referenceSurface(), options.componentMergeMethod);
-
-      track.parameters() = finalPars;
-      track.covariance() = finalCov;
+      track.parameters() = singleParams.parameters();
+      track.covariance() = singleParams.covariance().value();
 
       track.setReferenceSurface(params.referenceSurface().getSharedPtr());
 

Core/include/Acts/TrackFitting/GaussianSumFitter.hpp

@@ -11,32 +11,113 @@
 #include "Acts/EventData/TrackParameters.hpp"
 #include "Acts/TrackFitting/detail/GsfComponentMerging.hpp"
 
+#include <cstddef>
+#include <numeric>
+#include <ranges>
+
 namespace Acts {
 
 BoundTrackParameters MultiComponentBoundTrackParameters::merge(
     const ComponentMergeMethod method) const {
-  const bool hasCov = std::get<2>(m_components.front()).has_value();
+  if (size() == 0) {
+    throw std::logic_error(
+        "Cannot merge MultiComponentBoundTrackParameters with zero components");
+  }
+
+  if (size() == 1) {
+    return BoundTrackParameters(
+        m_surface, m_parameters[0],
+        hasCovariance() ? std::optional(m_covariances[0]) : std::nullopt,
+        m_particleHypothesis);
+  }
 
-  if (!hasCov) {
+  if (!hasCovariance()) {
     const BoundVector singleParams = detail::Gsf::mergeGaussianMixtureParams(
-        m_components,
-        [](const auto& cmp) -> std::tuple<double, const BoundVector&> {
-          return {std::get<0>(cmp), std::get<1>(cmp)};
+        std::views::iota(std::size_t{0}, size()),
+        [this](const std::size_t i) -> std::tuple<double, const BoundVector&> {
+          return {m_weights[i], m_parameters[i]};
         },
         *m_surface, method);
     return BoundTrackParameters(m_surface, singleParams, std::nullopt,
                                 m_particleHypothesis);
   }
 
   const auto [singleParams, singleCov] = detail::Gsf::mergeGaussianMixture(
-      m_components,
-      [](const auto& cmp)
+      std::views::iota(std::size_t{0}, size()),
+      [this](const std::size_t i)
           -> std::tuple<double, const BoundVector&, const BoundMatrix&> {
-        return {std::get<0>(cmp), std::get<1>(cmp), std::get<2>(cmp).value()};
+        return {m_weights[i], m_parameters[i], m_covariances[i]};
       },
       *m_surface, method);
   return BoundTrackParameters(m_surface, singleParams, singleCov,
                               m_particleHypothesis);
 }
 
+void MultiComponentBoundTrackParameters::reserve(const std::size_t n) {
+  m_weights.reserve(n);
+  m_parameters.reserve(n);
+  if (hasCovariance()) {
+    m_covariances.reserve(n);
+  }
+}
+
+void MultiComponentBoundTrackParameters::clear() {
+  m_weights.clear();
+  m_parameters.clear();
+  m_covariances.clear();
+}
+
+void MultiComponentBoundTrackParameters::pushComponent(
+    const double weight, const ParametersVector& params) {
+  if (hasCovariance()) {
+    throw std::logic_error(
+        "Cannot push component without covariance to "
+        "MultiComponentBoundTrackParameters with covariance");
+  }
+
+  m_weights.push_back(weight);
+  m_parameters.push_back(params);
+}
+
+void MultiComponentBoundTrackParameters::pushComponent(
+    const double weight, const ParametersVector& params,
+    const CovarianceMatrix& cov) {
+  if (!hasCovariance()) {
+    throw std::logic_error(
+        "Cannot push component with covariance to "
+        "MultiComponentBoundTrackParameters without covariance");
+  }
+
+  m_weights.push_back(weight);
+  m_parameters.push_back(params);
+  m_covariances.push_back(cov);
+}
+
+void MultiComponentBoundTrackParameters::pushComponent(
+    const double weight, const ParametersVector& params,
+    const std::optional<CovarianceMatrix>& cov) {
+  if (hasCovariance() != cov.has_value()) {
+  }
+
+  m_weights.push_back(weight);
+  m_parameters.push_back(params);
+  if (hasCovariance()) {
+    m_covariances.push_back(*cov);
+  }
+}
+
+void MultiComponentBoundTrackParameters::normalizeWeights() {
+  const double sumWeights =
+      std::accumulate(m_weights.begin(), m_weights.end(), 0.0);
+  if (sumWeights <= 0.0) {
+    throw std::logic_error(
+        "Cannot normalize weights of MultiComponentBoundTrackParameters: sum "
+        "of weights is not strictly positive: " +
+        std::to_string(sumWeights));
+  }
+  for (double& weight : m_weights) {
+    weight /= sumWeights;
+  }
+}
+
 }  // namespace Acts