feat: Add estimateTrackParamCovariance to Core

Core/CMakeLists.txt

@@ -110,6 +110,7 @@ add_subdirectory(src/MagneticField)
 add_subdirectory(src/Material)
 add_subdirectory(src/Navigation)
 add_subdirectory(src/Propagator)
+add_subdirectory(src/Seeding)
 add_subdirectory(src/Surfaces)
 add_subdirectory(src/TrackFinding)
 add_subdirectory(src/TrackFitting)

Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp

@@ -287,4 +287,41 @@ std::optional<BoundVector> estimateTrackParamsFromSeed(
   return params;
 }
 
+/// Configuration for the estimation of the covariance matrix of the track
+/// parameters with `estimateTrackParamCovariance`.
+struct EstimateTrackParamCovarianceConfig {
+  /// The initial sigmas for the track parameters
+  BoundVector initialSigmas = {1. * UnitConstants::mm,
+                               1. * UnitConstants::mm,
+                               1. * UnitConstants::degree,
+                               1. * UnitConstants::degree,
+                               1. * UnitConstants::e / UnitConstants::GeV,
+                               1. * UnitConstants::ns};
+
+  /// The initial relative uncertainty of the q/pt
+  double initialSigmaPtRel = 0.1;
+
+  /// The inflation factors for the variances of the track parameters
+  BoundVector initialVarInflation = {1., 1., 1., 1., 1., 1.};
+  /// The inflation factor for time uncertainty if the time parameter was not
+  /// estimated
+  double noTimeVarInflation = 100.;
+};
+
+/// Estimate the covariance matrix of the given track parameters based on the
+/// provided configuration. The assumption is that we can model the uncertainty
+/// of the track parameters as a diagonal matrix with the provided initial
+/// sigmas. The inflation factors are used to inflate the initial variances
+/// based on the provided configuration. The uncertainty of q/p is estimated
+/// based on the relative uncertainty of the q/pt and the theta uncertainty.
+///
+/// @param config is the configuration for the estimation
+/// @param params is the track parameters
+/// @param hasTime is true if the track parameters have time
+///
+/// @return the covariance matrix of the track parameters
+BoundMatrix estimateTrackParamCovariance(
+    const EstimateTrackParamCovarianceConfig& config, const BoundVector& params,
+    bool hasTime);
+
 }  // namespace Acts

Core/src/Seeding/CMakeLists.txt

@@ -0,0 +1 @@
+target_sources(ActsCore PRIVATE EstimateTrackParamsFromSeed.cpp)

Core/src/Seeding/EstimateTrackParamsFromSeed.cpp

@@ -0,0 +1,50 @@
+// 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/.
+
+#include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
+
+#include "Acts/Definitions/TrackParametrization.hpp"
+
+Acts::BoundMatrix Acts::estimateTrackParamCovariance(
+    const EstimateTrackParamCovarianceConfig& config, const BoundVector& params,
+    bool hasTime) {
+  assert((params[eBoundTheta] > 0 && params[eBoundTheta] < M_PI) &&
+         "Theta must be in the range (0, pi)");
+
+  BoundSquareMatrix result = BoundSquareMatrix::Zero();
+
+  for (std::size_t i = eBoundLoc0; i < eBoundSize; ++i) {
+    double sigma = config.initialSigmas[i];
+    double variance = sigma * sigma;
+
+    if (i == eBoundQOverP) {
+      // note that we rely on the fact that sigma theta is already computed
+      double varianceTheta = result(eBoundTheta, eBoundTheta);
+
+      // transverse momentum contribution
+      variance += std::pow(config.initialSigmaPtRel * params[eBoundQOverP], 2);
+
+      // theta contribution
+      variance +=
+          varianceTheta *
+          std::pow(params[eBoundQOverP] / std::tan(params[eBoundTheta]), 2);
+    }
+
+    if (i == eBoundTime && !hasTime) {
+      // Inflate the time uncertainty if no time measurement is available
+      variance *= config.noTimeVarInflation;
+    }
+
+    // Inflate the initial covariance
+    variance *= config.initialVarInflation[i];
+
+    result(i, i) = variance;
+  }
+
+  return result;
+}

Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp

@@ -32,48 +32,6 @@
 
 namespace ActsExamples {
 
-namespace {
-
-Acts::BoundSquareMatrix makeInitialCovariance(
-    const TrackParamsEstimationAlgorithm::Config& config,
-    const Acts::BoundVector& params, const SimSpacePoint& sp) {
-  Acts::BoundSquareMatrix result = Acts::BoundSquareMatrix::Zero();
-
-  for (std::size_t i = Acts::eBoundLoc0; i < Acts::eBoundSize; ++i) {
-    double sigma = config.initialSigmas[i];
-    double variance = sigma * sigma;
-
-    if (i == Acts::eBoundQOverP) {
-      // note that we rely on the fact that sigma theta is already computed
-      double varianceTheta = result(Acts::eBoundTheta, Acts::eBoundTheta);
-
-      // transverse momentum contribution
-      variance +=
-          std::pow(config.initialSigmaPtRel * params[Acts::eBoundQOverP], 2);
-
-      // theta contribution
-      variance +=
-          varianceTheta * std::pow(params[Acts::eBoundQOverP] /
-                                       std::tan(params[Acts::eBoundTheta]),
-                                   2);
-    }
-
-    // Inflate the time uncertainty if no time measurement is available
-    if (i == Acts::eBoundTime && !sp.t().has_value()) {
-      variance *= config.noTimeVarInflation;
-    }
-
-    // Inflate the initial covariance
-    variance *= config.initialVarInflation[i];
-
-    result(i, i) = variance;
-  }
-
-  return result;
-}
-
-}  // namespace
-
 TrackParamsEstimationAlgorithm::TrackParamsEstimationAlgorithm(
     TrackParamsEstimationAlgorithm::Config cfg, Acts::Logging::Level lvl)
     : IAlgorithm("TrackParamsEstimationAlgorithm", lvl), m_cfg(std::move(cfg)) {
@@ -166,8 +124,15 @@ ProcessCode TrackParamsEstimationAlgorithm::execute(
 
     const auto& params = optParams.value();
 
-    Acts::BoundSquareMatrix cov =
-        makeInitialCovariance(m_cfg, params, *bottomSP);
+    Acts::EstimateTrackParamCovarianceConfig config{
+        .initialSigmas =
+            Eigen::Map<const Acts::BoundVector>{m_cfg.initialSigmas.data()},
+        .initialSigmaPtRel = m_cfg.initialSigmaPtRel,
+        .initialVarInflation = Eigen::Map<const Acts::BoundVector>{
+            m_cfg.initialVarInflation.data()}};
+
+    Acts::BoundSquareMatrix cov = Acts::estimateTrackParamCovariance(
+        config, params, bottomSP->t().has_value());
 
     trackParameters.emplace_back(surface->getSharedPtr(), params, cov,
                                  m_cfg.particleHypothesis);

Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp

@@ -11,6 +11,7 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/EventData/ParticleHypothesis.hpp"
+#include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
 #include "Acts/Surfaces/PerigeeSurface.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Helpers.hpp"
@@ -128,31 +129,16 @@ ActsExamples::ProcessCode ActsExamples::ParticleSmearing::execute(
       Acts::BoundSquareMatrix cov = Acts::BoundSquareMatrix::Zero();
       if (m_cfg.initialSigmas) {
         // use the initial sigmas if set
-        for (std::size_t i = Acts::eBoundLoc0; i < Acts::eBoundSize; ++i) {
-          double sigma = (*m_cfg.initialSigmas)[i];
-          double variance = sigma * sigma;
-
-          if (i == Acts::eBoundQOverP) {
-            // note that we rely on the fact that sigma theta is already
-            // computed
-            double varianceTheta = cov(Acts::eBoundTheta, Acts::eBoundTheta);
 
-            // transverse momentum contribution
-            variance += std::pow(
-                m_cfg.initialSigmaPtRel * params[Acts::eBoundQOverP], 2);
+        Acts::EstimateTrackParamCovarianceConfig config{
+            .initialSigmas =
+                Eigen::Map<const Acts::BoundVector>{
+                    m_cfg.initialSigmas->data()},
+            .initialSigmaPtRel = m_cfg.initialSigmaPtRel,
+            .initialVarInflation = Eigen::Map<const Acts::BoundVector>{
+                m_cfg.initialVarInflation.data()}};
 
-            // theta contribution
-            variance += varianceTheta *
-                        std::pow(params[Acts::eBoundQOverP] /
-                                     std::tan(params[Acts::eBoundTheta]),
-                                 2);
-          }
-
-          // Inflate the initial covariance
-          variance *= m_cfg.initialVarInflation[i];
-
-          cov(i, i) = variance;
-        }
+        cov = Acts::estimateTrackParamCovariance(config, params, false);
       } else {
         // otherwise use the smearing sigmas
 

Examples/Python/tests/root_file_hashes.txt

@@ -39,16 +39,16 @@ test_ckf_tracks_example[generic-full_seeding]__performance_seeding_trees.root: 0
 test_ckf_tracks_example[generic-truth_estimated]__trackstates_ckf.root: a8c5c6f6c1e6303b887d47b509b7f71a2ffa5f38638fe46ce5bce76fd20d64ca
 test_ckf_tracks_example[generic-truth_estimated]__tracksummary_ckf.root: 417f7326e1e1bb4519f1378145ac733bdda6653eb9871fd69e455e0269d996a6
 test_ckf_tracks_example[generic-truth_estimated]__performance_seeding.root: 1facb05c066221f6361b61f015cdf0918e94d9f3fce2269ec7b6a4dffeb2bc7e
-test_ckf_tracks_example[generic-truth_smeared]__trackstates_ckf.root: de11c0868a70ade0dcc80465d4e6dcf1dd7fcf8149603b47ee7d87d862a6534a
-test_ckf_tracks_example[generic-truth_smeared]__tracksummary_ckf.root: f18e9ecce6d9585fd150c5aafc9ac225a5bab342aaab50a28283ba879691af1f
+test_ckf_tracks_example[generic-truth_smeared]__trackstates_ckf.root: edf0b06ce9ee0e4fcb153e41859af7b5153271de18f49a6842a23ad2d66b7e09
+test_ckf_tracks_example[generic-truth_smeared]__tracksummary_ckf.root: 06d6ae1d05cb611b19df3c59531997c9b0108f5ef6027d76c4827bd2d9edb921
 test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: 463d6aaed4d869652b5b184940e789cde0fb441bdd135813b85462a515e6480a
 test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: a8ad83a07b48d4cfcf70d0e6fdc3c8997eb03c1f8c2a7be27ea888b099000d79
 test_ckf_tracks_example[odd-full_seeding]__performance_seeding_trees.root: 43c58577aafe07645e5660c4f43904efadf91d8cda45c5c04c248bbe0f59814f
 test_ckf_tracks_example[odd-truth_estimated]__trackstates_ckf.root: 247dd581cc177625c0286718261c004e2149536d70c8281dfaf697879a84d76d
 test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: 1b08a80e73aedf5cf38a3a407794b82297bec37f556ad4efcda3489a1b17d4d2
 test_ckf_tracks_example[odd-truth_estimated]__performance_seeding.root: 1a36b7017e59f1c08602ef3c2cb0483c51df248f112e3780c66594110719c575
-test_ckf_tracks_example[odd-truth_smeared]__trackstates_ckf.root: 7adfc2bf5ee35a126b713187dd8b11f4497cf864a4a83e57a40885688974413e
-test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: 7a9de8a8bd1c09f7b4d1c547f824af6c8123afb044dd429180b0d13e47d6f975
+test_ckf_tracks_example[odd-truth_smeared]__trackstates_ckf.root: a9621b535ea2912d172142394f51f68e4e7dc255b32d479d6305fa599152b420
+test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: af1a6bb16a070db7ed8043e2188d56f0034843099fc3c332731c4cf86ba39c57
 test_vertex_fitting_reading[Truth-False-100]__performance_vertexing.root: 76ef6084d758dfdfc0151ddec2170e12d73394424e3dac4ffe46f0f339ec8293
 test_vertex_fitting_reading[Iterative-False-100]__performance_vertexing.root: 60372210c830a04f95ceb78c6c68a9b0de217746ff59e8e73053750c837b57eb
 test_vertex_fitting_reading[Iterative-True-100]__performance_vertexing.root: e34f217d524a5051dbb04a811d3407df3ebe2cc4bb7f54f6bda0847dbd7b52c3