fix: Compute complete covariance for max weight reduction (#1927)

In case of returning the parameters of the component with the largest weight, it is of course not correct to return also the covariance matrix of this component, since it does not reflect the uncertainty of the whole distribution. Rather, return the the covariance computed from all components.
acts-project · Mar 13, 2023 · 7f0429f · 7f0429f
1 parent 250ea03
commit 7f0429f
Show file tree

Hide file tree

Showing 3 changed files with 40 additions and 44 deletions.
diff --git a/CI/physmon/reference/performance_gsf.root b/CI/physmon/reference/performance_gsf.root
diff --git a/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp b/Core/include/Acts/Propagator/MultiEigenStepperLoop.ipp
@@ -20,55 +20,51 @@ auto MultiEigenStepperLoop<E, R, A>::boundState(
   if (numberComponents(state) == 1) {
     return SingleStepper::boundState(state.components.front().state, surface,
                                      transportCov, freeToBoundCorrection);
-  } else if (m_finalReductionMethod == FinalReductionMethod::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::boundState(cmpIt->state, surface, transportCov,
-                                     freeToBoundCorrection);
-  } else {
-    SmallVector<std::tuple<double, BoundVector, BoundSymMatrix>> states;
-    double accumulatedPathLength = 0.0;
-    int failedBoundTransforms = 0;
-
-    for (auto i = 0ul; i < numberComponents(state); ++i) {
-      auto bs = SingleStepper::boundState(state.components[i].state, surface,
-                                          transportCov, freeToBoundCorrection);
-
-      if (bs.ok()) {
-        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 {
-        failedBoundTransforms++;
-      }
-    }
+  }
 
-    if (states.empty()) {
-      return MultiStepperError::AllComponentsConversionToBoundFailed;
-    }
+  SmallVector<std::tuple<double, BoundVector, BoundSymMatrix>> states;
+  double accumulatedPathLength = 0.0;
 
-    if (failedBoundTransforms > 0) {
-      return MultiStepperError::SomeComponentsConversionToBoundFailed;
+  for (auto i = 0ul; i < numberComponents(state); ++i) {
+    auto bs = SingleStepper::boundState(state.components[i].state, surface,
+                                        transportCov, freeToBoundCorrection);
+
+    if (bs.ok()) {
+      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;
     }
+  }
 
-    auto [params, cov] =
-        detail::angleDescriptionSwitch(surface, [&](const auto& desc) {
-          return detail::combineGaussianMixture(states, Acts::Identity{}, desc);
-        });
+  if (states.empty()) {
+    return MultiStepperError::AllComponentsConversionToBoundFailed;
+  }
 
-    std::optional<BoundSymMatrix> finalCov = std::nullopt;
-    if (cov != BoundSymMatrix::Zero()) {
-      finalCov = cov;
-    }
+  const auto [mean, cov] =
+      detail::angleDescriptionSwitch(surface, [&](const auto& desc) {
+        return detail::combineGaussianMixture(states, Acts::Identity{}, desc);
+      });
 
-    return BoundState{BoundTrackParameters(surface.getSharedPtr(), params, cov),
-                      Jacobian::Zero(), accumulatedPathLength};
+  const auto finalPars =
+      (m_finalReductionMethod == FinalReductionMethod::eMaxWeight)
+          ? std::get<BoundVector>(*std::max_element(
+                states.begin(), states.end(),
+                [](const auto& a, const auto& b) {
+                  return std::get<double>(a) < std::get<double>(b);
+                }))
+          : mean;
+
+  std::optional<BoundSymMatrix> finalCov = std::nullopt;
+  if (cov != BoundSymMatrix::Zero()) {
+    finalCov = cov;
   }
+
+  return BoundState{
+      BoundTrackParameters(surface.getSharedPtr(), finalPars, finalCov),
+      Jacobian::Zero(), accumulatedPathLength};
 }
 
 template <typename E, typename R, typename A>

diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt
@@ -35,9 +35,9 @@ test_truth_tracking_kalman[odd-1000.0]__trackstates_fitter.root: c76b77c50bd0545
 test_truth_tracking_kalman[odd-1000.0]__tracksummary_fitter.root: 0f26a2e7c67407e0c7092b2d861753d658529147ff2dc4ea968a9c9f28a19bf9
 test_truth_tracking_kalman[odd-1000.0]__performance_track_finder.root: dacf32346bca4c7010ca6a21f4868e2a9cdbc764d0918cb498cefdf322567716
 test_truth_tracking_gsf[generic]__trackstates_gsf.root: d69377c836d1825864aeca41bc7b880ce86f38f3692bd186667a3a36377366cd
-test_truth_tracking_gsf[generic]__tracksummary_gsf.root: d01a5841520fdbd2d2c8483468ffc730c8e1dc459eb2179eb91808e23e5039cf
+test_truth_tracking_gsf[generic]__tracksummary_gsf.root: 6fb6f8980f6407bfc612f28187ebd504f3a75b54e427427313bdf2f60e4fa7af
 test_truth_tracking_gsf[odd]__trackstates_gsf.root: 24259acbb2d31acfb88ee5704cf7174eb5cec312ab7e32d1f3102005e33581dd
-test_truth_tracking_gsf[odd]__tracksummary_gsf.root: e89871a8f2607b7da7176e5f36cde861231d2e24885bdb4a7041d2f73f3517c9
+test_truth_tracking_gsf[odd]__tracksummary_gsf.root: b3560f2b3bce7346d2da4bcd5f9a6bcffb4cd66e1c5ac62a0dcb59ae4dd728cf
 test_particle_gun__particles.root: 8549ba6e20338004ab8ba299fc65e1ee5071985b46df8f77f887cb6fef56a8ec
 test_material_mapping__material-map_tracks.root: 90a7deab073691ef431a3711e9cba40e902ce4de92a9b57162d84fbef5f97039
 test_material_mapping__propagation-material.root: c60b17ddc4921847bb3e9ebf58de10322d0d3888596ee1106f160055647d8445