refactor: no optional covariance in GSF (#1826)

The GSF in principle worked with `std::optional<BoundSymMatrix>` for the covariance in almost all cases. However, it does not make any sense to run the GSF without covariance, and removing this does also make the code simpler.
Actually, it was assumed to have a value in almost all cases, so this should also make the code a bit saver.

Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp

@@ -28,7 +28,7 @@ auto MultiEigenStepperLoop<E, R, A>::boundState(
     return SingleStepper::boundState(cmpIt->state, surface, transportCov,
                                      freeToBoundCorrection);
   } else {
-    SmallVector<std::pair<double, BoundTrackParameters>> states;
+    SmallVector<std::tuple<double, BoundVector, BoundSymMatrix>> states;
     double accumulatedPathLength = 0.0;
     int failedBoundTransforms = 0;
 
@@ -37,8 +37,10 @@ auto MultiEigenStepperLoop<E, R, A>::boundState(
                                           transportCov, freeToBoundCorrection);
 
       if (bs.ok()) {
-        states.push_back(
-            {state.components[i].weight, std::get<BoundTrackParameters>(*bs)});
+        const auto& btp = std::get<BoundTrackParameters>(*bs);
+        states.emplace_back(
+            state.components[i].weight, btp.parameters(),
+            btp.covariance().value_or(Acts::BoundSymMatrix::Zero()));
         accumulatedPathLength +=
             std::get<double>(*bs) * state.components[i].weight;
       } else {
@@ -54,16 +56,16 @@ auto MultiEigenStepperLoop<E, R, A>::boundState(
       return MultiStepperError::SomeComponentsConversionToBoundFailed;
     }
 
-    const auto proj = [&](const auto& wbs) {
-      return std::tie(wbs.first, wbs.second.parameters(),
-                      wbs.second.covariance());
-    };
-
     auto [params, cov] =
         detail::angleDescriptionSwitch(surface, [&](const auto& desc) {
-          return detail::combineGaussianMixture(states, proj, desc);
+          return detail::combineGaussianMixture(states, Acts::Identity{}, desc);
         });
 
+    std::optional<BoundSymMatrix> finalCov = std::nullopt;
+    if (cov != BoundSymMatrix::Zero()) {
+      finalCov = cov;
+    }
+
     return BoundState{BoundTrackParameters(surface.getSharedPtr(), params, cov),
                       Jacobian::Zero(), accumulatedPathLength};
   }
@@ -90,6 +92,7 @@ auto MultiEigenStepperLoop<E, R, A>::curvilinearState(State& state,
     ActsScalar qop = 0.0;
     BoundSymMatrix cov = BoundSymMatrix::Zero();
     ActsScalar pathLenth = 0.0;
+    ActsScalar sumOfWeights = 0.0;
 
     for (auto i = 0ul; i < numberComponents(state); ++i) {
       const auto [cp, jac, pl] = SingleStepper::curvilinearState(
@@ -102,10 +105,23 @@ auto MultiEigenStepperLoop<E, R, A>::curvilinearState(State& state,
         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<BoundSymMatrix> finalCov = std::nullopt;
+    if (cov != BoundSymMatrix::Zero()) {
+      finalCov = cov;
     }
 
-    return CurvilinearState{CurvilinearTrackParameters(pos4, dir, qop, cov),
-                            Jacobian::Zero(), pathLenth};
+    return CurvilinearState{
+        CurvilinearTrackParameters(pos4, dir, qop, finalCov), Jacobian::Zero(),
+        pathLenth};
   }
 }
 

Core/include/Acts/TrackFitting/detail/GsfActor.hpp

@@ -133,7 +133,7 @@ struct GsfActor {
   struct ParameterCache {
     ActsScalar weight = 0;
     BoundVector boundPars;
-    std::optional<BoundSymMatrix> boundCov;
+    BoundSymMatrix boundCov;
   };
 
   struct TemporaryStates {
@@ -422,25 +422,22 @@ struct GsfActor {
       new_pars[eBoundQOverP] = old_bound.charge() / (p_prev + delta_p);
 
       // compute inverse variance of p from mixture and update covariance
-      auto new_cov = old_bound.covariance();
-
-      if (new_cov.has_value()) {
-        const auto varInvP = [&]() {
-          if (state.stepping.navDir == NavigationDirection::Forward) {
-            const auto f = 1. / (p_prev * gaussian.mean);
-            return f * f * gaussian.var;
-          } else {
-            return gaussian.var / (p_prev * p_prev);
-          }
-        }();
-
-        (*new_cov)(eBoundQOverP, eBoundQOverP) += varInvP;
-        throw_assert(
-            std::isfinite((*new_cov)(eBoundQOverP, eBoundQOverP)),
-            "cov not finite, varInvP=" << varInvP << ", p_prev=" << p_prev
-                                       << ", gaussian.mean=" << gaussian.mean
-                                       << ", gaussian.var=" << gaussian.var);
-      }
+      auto new_cov = old_bound.covariance().value();
+
+      const auto varInvP = [&]() {
+        if (state.stepping.navDir == NavigationDirection::Forward) {
+          const auto f = 1. / (p_prev * gaussian.mean);
+          return f * f * gaussian.var;
+        } else {
+          return gaussian.var / (p_prev * p_prev);
+        }
+      }();
+
+      new_cov(eBoundQOverP, eBoundQOverP) += varInvP;
+      throw_assert(std::isfinite(new_cov(eBoundQOverP, eBoundQOverP)),
+                   "cov not finite, varInvP="
+                       << varInvP << ", p_prev=" << p_prev << ", gaussian.mean="
+                       << gaussian.mean << ", gaussian.var=" << gaussian.var);
 
       // Set the remaining things and push to vector
       componentCaches.push_back(
@@ -612,7 +609,7 @@ struct GsfActor {
     for (const auto& idx : tmpStates.tips) {
       const auto [w, p, c] = proj(idx);
       if (w > 0.0) {
-        v.push_back({w, p, *c});
+        v.push_back({w, p, c});
       }
     }
 
@@ -746,9 +743,9 @@ struct GsfActor {
       // We set predicted & filtered the same so that the fields are not
       // uninitialized when not finding this state in the reverse pass.
       proxy.predicted() = filtMean;
-      proxy.predictedCovariance() = filtCov.value();
+      proxy.predictedCovariance() = filtCov;
       proxy.filtered() = filtMean;
-      proxy.filteredCovariance() = filtCov.value();
+      proxy.filteredCovariance() = filtCov;
     } else {
       assert((result.currentTip != MultiTrajectoryTraits::kInvalid &&
               "tip not valid"));
@@ -764,7 +761,7 @@ struct GsfActor {
                   });
 
               trackState.filtered() = filtMean;
-              trackState.filteredCovariance() = filtCov.value();
+              trackState.filteredCovariance() = filtCov;
               return false;
             }
             return true;

Core/include/Acts/TrackFitting/detail/GsfUtils.hpp

@@ -225,13 +225,13 @@ struct MultiTrajectoryProjector {
     switch (type) {
       case StatesType::ePredicted:
         return std::make_tuple(weights.at(idx), proxy.predicted(),
-                               std::optional{proxy.predictedCovariance()});
+                               proxy.predictedCovariance());
       case StatesType::eFiltered:
         return std::make_tuple(weights.at(idx), proxy.filtered(),
-                               std::optional{proxy.filteredCovariance()});
+                               proxy.filteredCovariance());
       case StatesType::eSmoothed:
         return std::make_tuple(weights.at(idx), proxy.smoothed(),
-                               std::optional{proxy.smoothedCovariance()});
+                               proxy.smoothedCovariance());
     }
   }
 };

Core/include/Acts/TrackFitting/detail/KLMixtureReduction.hpp

@@ -26,8 +26,8 @@ auto computeSymmetricKlDivergence(const component_t &a, const component_t &b,
                                   const component_projector_t &proj) {
   const auto parsA = proj(a).boundPars[eBoundQOverP];
   const auto parsB = proj(b).boundPars[eBoundQOverP];
-  const auto covA = (*proj(a).boundCov)(eBoundQOverP, eBoundQOverP);
-  const auto covB = (*proj(b).boundCov)(eBoundQOverP, eBoundQOverP);
+  const auto covA = proj(a).boundCov(eBoundQOverP, eBoundQOverP);
+  const auto covB = proj(b).boundCov(eBoundQOverP, eBoundQOverP);
 
   throw_assert(covA != 0.0, "");
   throw_assert(std::isfinite(covA), "");

Core/include/Acts/Utilities/detail/gaussian_mixture_helpers.hpp

@@ -84,7 +84,7 @@ auto combineCov(const components_t components, const ActsVector<D> &mean,
   for (const auto &cmp : components) {
     const auto &[weight_l, pars_l, cov_l] = projector(cmp);
 
-    cov += weight_l * *cov_l;
+    cov += weight_l * cov_l;
 
     ActsVector<D> diff = pars_l - mean;
 
@@ -133,7 +133,7 @@ auto combineGaussianMixture(const components_t components,
 
   // Assert type properties
   using ParsType = std::decay_t<decltype(beginPars)>;
-  using CovType = std::decay_t<decltype(*beginCov)>;
+  using CovType = std::decay_t<decltype(beginCov)>;
   using WeightType = std::decay_t<decltype(beginWeight)>;
 
   constexpr int D = ParsType::RowsAtCompileTime;
@@ -152,11 +152,11 @@ auto combineGaussianMixture(const components_t components,
 #endif
 
   // Define the return type
-  using RetType = std::tuple<ActsVector<D>, std::optional<ActsSymMatrix<D>>>;
+  using RetType = std::tuple<ActsVector<D>, ActsSymMatrix<D>>;
 
   // Early return in case of range with length 1
   if (components.size() == 1) {
-    return RetType{beginPars / beginWeight, *beginCov / beginWeight};
+    return RetType{beginPars / beginWeight, beginCov / beginWeight};
   }
 
   // Zero initialized values for aggregation
@@ -194,14 +194,10 @@ auto combineGaussianMixture(const components_t components,
 
   std::apply([&](auto... dsc) { (wrap(dsc), ...); }, angleDesc);
 
-  if (not beginCov) {
-    return RetType{mean, std::nullopt};
-  } else {
-    const auto cov =
-        combineCov(components, mean, sumOfWeights, projector, angleDesc);
+  const auto cov =
+      combineCov(components, mean, sumOfWeights, projector, angleDesc);
 
-    return RetType{mean, cov};
-  }
+  return RetType{mean, cov};
 }
 
 }  // namespace detail

Tests/UnitTests/Core/TrackFitting/GsfComponentMergingTests.cpp

@@ -32,7 +32,7 @@ template <int D>
 struct DummyComponent {
   Acts::ActsScalar weight = 0;
   Acts::ActsVector<D> boundPars;
-  std::optional<Acts::ActsSymMatrix<D>> boundCov;
+  Acts::ActsSymMatrix<D> boundCov;
 };
 
 // A Multivariate distribution object working in the same way as the
@@ -72,7 +72,7 @@ auto sampleFromMultivariate(const std::vector<DummyComponent<D>> &cmps,
   std::vector<MultiNormal> dists;
   std::vector<double> weights;
   for (const auto &cmp : cmps) {
-    dists.push_back(MultiNormal(cmp.boundPars, *cmp.boundCov));
+    dists.push_back(MultiNormal(cmp.boundPars, cmp.boundCov));
     weights.push_back(cmp.weight);
   }
 
@@ -205,6 +205,9 @@ void test_surface(const Surface &surface, const angle_description_t &desc,
               detail::wrap_periodic(phi + dphi, -M_PI, 2 * M_PI);
           a.boundPars[eBoundTheta] = theta + dtheta;
 
+          // We don't look at covariance in this test
+          a.boundCov = BoundSymMatrix::Zero();
+
           cmps.push_back(a);
           ++p_it;
         }
@@ -213,9 +216,6 @@ void test_surface(const Surface &surface, const angle_description_t &desc,
       const auto [mean_approx, cov_approx] =
           detail::combineGaussianMixture(cmps, proj, desc);
 
-      // We don't have a boundCovariance in this test
-      BOOST_CHECK(not cov_approx);
-
       const auto mean_ref = meanFromFree(cmps, surface);
 
       CHECK_CLOSE_MATRIX(mean_approx, mean_ref, expectedError);
@@ -228,13 +228,11 @@ BOOST_AUTO_TEST_CASE(test_with_data) {
   std::vector<DummyComponent<2>> cmps(2);
 
   cmps[0].boundPars << 1.0, 1.0;
-  cmps[0].boundCov = ActsSymMatrix<2>::Zero();
-  *cmps[0].boundCov << 1.0, 0.0, 0.0, 1.0;
+  cmps[0].boundCov << 1.0, 0.0, 0.0, 1.0;
   cmps[0].weight = 0.5;
 
   cmps[1].boundPars << -2.0, -2.0;
-  cmps[1].boundCov = ActsSymMatrix<2>::Zero();
-  *cmps[1].boundCov << 1.0, 1.0, 1.0, 2.0;
+  cmps[1].boundCov << 1.0, 1.0, 1.0, 2.0;
   cmps[1].weight = 0.5;
 
   const auto samples = sampleFromMultivariate(cmps, 10000, gen);
@@ -245,21 +243,19 @@ BOOST_AUTO_TEST_CASE(test_with_data) {
       detail::combineGaussianMixture(cmps, Identity{}, std::tuple<>{});
 
   CHECK_CLOSE_MATRIX(mean_data, mean_test, 1.e-1);
-  CHECK_CLOSE_MATRIX(boundCov_data, *boundCov_test, 1.e-1);
+  CHECK_CLOSE_MATRIX(boundCov_data, boundCov_test, 1.e-1);
 }
 
 BOOST_AUTO_TEST_CASE(test_with_data_circular) {
   std::mt19937 gen(42);
   std::vector<DummyComponent<2>> cmps(2);
 
   cmps[0].boundPars << 175_degree, 5_degree;
-  cmps[0].boundCov = ActsSymMatrix<2>::Zero();
-  *cmps[0].boundCov << 20_degree, 0.0, 0.0, 20_degree;
+  cmps[0].boundCov << 20_degree, 0.0, 0.0, 20_degree;
   cmps[0].weight = 0.5;
 
   cmps[1].boundPars << -175_degree, -5_degree;
-  cmps[1].boundCov = ActsSymMatrix<2>::Zero();
-  *cmps[1].boundCov << 20_degree, 20_degree, 20_degree, 40_degree;
+  cmps[1].boundCov << 20_degree, 20_degree, 20_degree, 40_degree;
   cmps[1].weight = 0.5;
 
   const auto samples = sampleFromMultivariate(cmps, 10000, gen);
@@ -281,7 +277,7 @@ BOOST_AUTO_TEST_CASE(test_with_data_circular) {
                                                    2 * M_PI)) < 1.e-1);
   BOOST_CHECK(std::abs(detail::difference_periodic(mean_data[1], mean_test[1],
                                                    2 * M_PI)) < 1.e-1);
-  CHECK_CLOSE_MATRIX(boundCov_data, *boundCov_test, 1.e-1);
+  CHECK_CLOSE_MATRIX(boundCov_data, boundCov_test, 1.e-1);
 }
 
 BOOST_AUTO_TEST_CASE(test_plane_surface) {

Tests/UnitTests/Core/TrackFitting/GsfMixtureReductionTests.cpp

@@ -16,7 +16,7 @@ using namespace Acts::UnitLiterals;
 struct DummyComponent {
   double weight = 0.0;
   BoundVector boundPars = BoundVector::Zero();
-  std::optional<BoundSymMatrix> boundCov = BoundSymMatrix::Identity();
+  BoundSymMatrix boundCov = BoundSymMatrix::Identity();
 };
 
 BOOST_AUTO_TEST_CASE(test_distance_matrix_min_distance) {