diff --git a/Core/CMakeLists.txt b/Core/CMakeLists.txt
index eeb6bc8b92a..554656f514c 100644
--- a/Core/CMakeLists.txt
+++ b/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)
diff --git a/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp b/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp
index 27b411a99cc..f624cd88aee 100644
--- a/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp
+++ b/Core/include/Acts/EventData/GenericBoundTrackParameters.hpp
@@ -250,6 +250,17 @@ class GenericBoundTrackParameters {
     return m_surface->referenceFrame(geoCtx, position(geoCtx), momentum());
   }
 
+  /// Reflect the parameters in place.
+  void reflectInPlace() { m_params = reflectBoundParameters(m_params); }
+
+  /// Reflect the parameters.
+  /// @return Reflected parameters.
+  GenericBoundTrackParameters<ParticleHypothesis> reflect() const {
+    GenericBoundTrackParameters<ParticleHypothesis> reflected = *this;
+    reflected.reflectInPlace();
+    return reflected;
+  }
+
  private:
   BoundVector m_params;
   std::optional<BoundSquareMatrix> m_cov;
diff --git a/Core/include/Acts/EventData/GenericCurvilinearTrackParameters.hpp b/Core/include/Acts/EventData/GenericCurvilinearTrackParameters.hpp
index a42ac2f116f..6e4ed9bae7a 100644
--- a/Core/include/Acts/EventData/GenericCurvilinearTrackParameters.hpp
+++ b/Core/include/Acts/EventData/GenericCurvilinearTrackParameters.hpp
@@ -111,6 +111,14 @@ class GenericCurvilinearTrackParameters
   Vector3 position() const {
     return GenericBoundTrackParameters<ParticleHypothesis>::position({});
   }
+
+  /// Reflect the parameters.
+  /// @return Reflected parameters.
+  GenericCurvilinearTrackParameters<ParticleHypothesis> reflect() const {
+    GenericCurvilinearTrackParameters<ParticleHypothesis> reflected = *this;
+    reflected.reflectInPlace();
+    return reflected;
+  }
 };
 
 }  // namespace Acts
diff --git a/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp b/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp
index 1e7846318ef..214aa2e1551 100644
--- a/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp
+++ b/Core/include/Acts/EventData/GenericFreeTrackParameters.hpp
@@ -12,9 +12,11 @@
 #include "Acts/Definitions/Common.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/EventData/TrackParametersConcept.hpp"
+#include "Acts/EventData/TransformationHelpers.hpp"
 #include "Acts/EventData/detail/PrintParameters.hpp"
 #include "Acts/Utilities/MathHelpers.hpp"
 #include "Acts/Utilities/UnitVectors.hpp"
+#include "Acts/Utilities/VectorHelpers.hpp"
 
 #include <cassert>
 #include <cmath>
@@ -55,6 +57,29 @@ class GenericFreeTrackParameters {
         m_cov(std::move(cov)),
         m_particleHypothesis(std::move(particleHypothesis)) {}
 
+  /// Construct from four-position, direction, absolute momentum, and charge.
+  ///
+  /// @param pos4 Track position/time four-vector
+  /// @param dir Track direction three-vector; normalization is ignored.
+  /// @param qOverP Charge over momentum
+  /// @param cov Free parameters covariance matrix
+  /// @param particleHypothesis Particle hypothesis
+  GenericFreeTrackParameters(const Vector4& pos4, const Vector3& dir,
+                             Scalar qOverP, std::optional<CovarianceMatrix> cov,
+                             ParticleHypothesis particleHypothesis)
+      : m_params(FreeVector::Zero()),
+        m_cov(std::move(cov)),
+        m_particleHypothesis(std::move(particleHypothesis)) {
+    m_params[eFreePos0] = pos4[ePos0];
+    m_params[eFreePos1] = pos4[ePos1];
+    m_params[eFreePos2] = pos4[ePos2];
+    m_params[eFreeTime] = pos4[eTime];
+    m_params[eFreeDir0] = dir[eMom0];
+    m_params[eFreeDir1] = dir[eMom1];
+    m_params[eFreeDir2] = dir[eMom2];
+    m_params[eFreeQOverP] = qOverP;
+  }
+
   /// Construct from four-position, angles, absolute momentum, and charge.
   ///
   /// @param pos4 Track position/time four-vector
@@ -135,9 +160,9 @@ class GenericFreeTrackParameters {
   Scalar time() const { return m_params[eFreeTime]; }
 
   /// Phi direction.
-  Scalar phi() const { return phi(direction()); }
+  Scalar phi() const { return VectorHelpers::phi(direction()); }
   /// Theta direction.
-  Scalar theta() const { return theta(direction()); }
+  Scalar theta() const { return VectorHelpers::theta(direction()); }
   /// Charge over momentum.
   Scalar qOverP() const { return m_params[eFreeQOverP]; }
 
@@ -175,6 +200,17 @@ class GenericFreeTrackParameters {
     return m_particleHypothesis;
   }
 
+  /// Reflect the parameters in place.
+  void reflectInPlace() { m_params = reflectFreeParameters(m_params); }
+
+  /// Reflect the parameters.
+  /// @return Reflected parameters.
+  GenericFreeTrackParameters<ParticleHypothesis> reflect() const {
+    GenericFreeTrackParameters<ParticleHypothesis> reflected = *this;
+    reflected.reflectInPlace();
+    return reflected;
+  }
+
  private:
   FreeVector m_params;
   std::optional<FreeSquareMatrix> m_cov;
diff --git a/Core/include/Acts/EventData/TrackProxy.hpp b/Core/include/Acts/EventData/TrackProxy.hpp
index ce674e9cfa4..5939022e2e8 100644
--- a/Core/include/Acts/EventData/TrackProxy.hpp
+++ b/Core/include/Acts/EventData/TrackProxy.hpp
@@ -324,7 +324,8 @@ class TrackProxy {
     return std::distance(tsRange.begin(), tsRange.end());
   }
 
-  /// Return the number of measurements for the track. Const version
+  /// Return a mutable reference to the number of measurements for the track.
+  /// Mutable version
   /// @note Only available if the track proxy is not read-only
   /// @return The number of measurements
   unsigned int& nMeasurements()
@@ -333,8 +334,7 @@ class TrackProxy {
     return component<unsigned int, hashString("nMeasurements")>();
   }
 
-  /// Return a mutable reference to the number of measurements for the track.
-  /// Mutable version
+  /// Return the number of measurements for the track. Const version
   /// @return The number of measurements
   unsigned int nMeasurements() const {
     return component<unsigned int, hashString("nMeasurements")>();
diff --git a/Core/include/Acts/EventData/TransformationHelpers.hpp b/Core/include/Acts/EventData/TransformationHelpers.hpp
index 7fa297f4abc..39240bf469e 100644
--- a/Core/include/Acts/EventData/TransformationHelpers.hpp
+++ b/Core/include/Acts/EventData/TransformationHelpers.hpp
@@ -13,11 +13,39 @@
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Utilities/Result.hpp"
+#include "Acts/Utilities/detail/periodic.hpp"
 
 namespace Acts {
 
 class Surface;
 
+/// Reflect bound track parameters.
+///
+/// @param boundParams Bound track parameters vector
+/// @return Reflected bound track parameters vector
+inline BoundVector reflectBoundParameters(const BoundVector& boundParams) {
+  BoundVector reflected = boundParams;
+  auto [phi, theta] = detail::normalizePhiTheta(
+      boundParams[eBoundPhi] - M_PI, M_PI - boundParams[eBoundTheta]);
+  reflected[eBoundPhi] = phi;
+  reflected[eBoundTheta] = theta;
+  reflected[eBoundQOverP] = -boundParams[eBoundQOverP];
+  return reflected;
+}
+
+/// Reflect free track parameters.
+///
+/// @param freeParams Free track parameters vector
+/// @return Reflected free track parameters vector
+inline FreeVector reflectFreeParameters(const FreeVector& freeParams) {
+  FreeVector reflected = freeParams;
+  reflected[eFreeDir0] = -freeParams[eFreeDir0];
+  reflected[eFreeDir1] = -freeParams[eFreeDir1];
+  reflected[eFreeDir2] = -freeParams[eFreeDir2];
+  reflected[eFreeQOverP] = -freeParams[eFreeQOverP];
+  return reflected;
+}
+
 /// Transform bound track parameters into equivalent free track parameters.
 ///
 /// @param surface Surface onto which the input parameters are bound
diff --git a/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp b/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp
index 7a5c597f626..f3850e102ce 100644
--- a/Core/include/Acts/Seeding/EstimateTrackParamsFromSeed.hpp
+++ b/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
diff --git a/Core/include/Acts/TrackFinding/TrackSelector.hpp b/Core/include/Acts/TrackFinding/TrackSelector.hpp
index d0e8526dc35..e782b7feb9d 100644
--- a/Core/include/Acts/TrackFinding/TrackSelector.hpp
+++ b/Core/include/Acts/TrackFinding/TrackSelector.hpp
@@ -142,7 +142,7 @@ class TrackSelector {
     std::vector<Config> cutSets = {};
 
     /// Eta bin edges for varying cuts by eta
-    std::vector<double> absEtaEdges = {};
+    std::vector<double> absEtaEdges = {0, inf};
 
     /// Get the number of eta bins
     /// @return Number of eta bins
@@ -150,7 +150,7 @@ class TrackSelector {
 
     /// Construct an empty (accepts everything) configuration.
     /// Results in a single cut set and one abs eta bin from 0 to infinity.
-    EtaBinnedConfig() : cutSets{{}}, absEtaEdges{{0, inf}} {};
+    EtaBinnedConfig() : cutSets{{}} {};
 
     /// Constructor to create a config object that is not upper-bounded.
     /// This is useful to use the "fluent" API to populate the configuration.
@@ -163,13 +163,12 @@ class TrackSelector {
     /// @param absEtaEdgesIn is the vector of eta bin edges
     EtaBinnedConfig(std::vector<double> absEtaEdgesIn)
         : absEtaEdges{std::move(absEtaEdgesIn)} {
-      cutSets.resize(absEtaEdges.size() - 1);
+      cutSets.resize(nEtaBins());
     }
 
     /// Auto-converting constructor from a single cut configuration.
     /// Results in a single absolute eta bin from 0 to infinity.
-    EtaBinnedConfig(Config cutSet)
-        : cutSets{std::move(cutSet)}, absEtaEdges{{0, inf}} {}
+    EtaBinnedConfig(Config cutSet) : cutSets{std::move(cutSet)} {}
 
     /// Add a new eta bin with the given upper bound.
     /// @param etaMax Upper bound of the new eta bin
@@ -195,11 +194,17 @@ class TrackSelector {
     /// @return True if the configuration has a bin for the given eta
     bool hasCuts(double eta) const;
 
-    /// Get the index of the eta bin for a given eta
+    /// Get the index of the eta bin for a given eta.
+    /// throws an exception if Eta is outside the abs eta bin edges.
     /// @param eta Eta value
     /// @return Index of the eta bin
     std::size_t binIndex(double eta) const;
 
+    /// Get the index of the eta bin for a given eta
+    /// @param eta Eta value
+    /// @return Index of the eta bin, or >= nEtaBins() if Eta is outside the abs eta bin edges.
+    std::size_t binIndexNoCheck(double eta) const;
+
     /// Get the cuts for a given eta
     /// @param eta Eta value
     /// @return Cuts for the given eta
@@ -237,8 +242,7 @@ class TrackSelector {
 
  private:
   EtaBinnedConfig m_cfg;
-  bool m_isUnbinned;
-  bool m_noEtaCuts;
+  bool m_isUnbinned = false;
 };
 
 inline TrackSelector::Config& TrackSelector::Config::loc0(double min,
@@ -350,14 +354,22 @@ inline bool TrackSelector::EtaBinnedConfig::hasCuts(double eta) const {
 }
 
 inline std::size_t TrackSelector::EtaBinnedConfig::binIndex(double eta) const {
-  if (!hasCuts(eta)) {
+  std::size_t index = binIndexNoCheck(eta);
+  if (!(index < nEtaBins())) {
     throw std::invalid_argument{"Eta is outside the abs eta bin edges"};
   }
+  return index;
+}
 
+inline std::size_t TrackSelector::EtaBinnedConfig::binIndexNoCheck(
+    double eta) const {
   auto binIt =
       std::upper_bound(absEtaEdges.begin(), absEtaEdges.end(), std::abs(eta));
-  std::size_t index = std::distance(absEtaEdges.begin(), binIt) - 1;
-  return index;
+  std::size_t index = std::distance(absEtaEdges.begin(), binIt);
+  if (index == 0) {
+    index = absEtaEdges.size() + 1;  // positive value to check for underflow
+  }
+  return index - 1;
 }
 
 inline const TrackSelector::Config& TrackSelector::EtaBinnedConfig::getCuts(
@@ -428,8 +440,8 @@ bool TrackSelector::isValidTrack(const track_proxy_t& track) const {
 
   return track.hasReferenceSurface() &&
          within(track.transverseMomentum(), cuts.ptMin, cuts.ptMax) &&
-         (m_noEtaCuts || (within(absEta(), cuts.absEtaMin, cuts.absEtaMax) &&
-                          within(_eta, cuts.etaMin, cuts.etaMax))) &&
+         (!m_isUnbinned || (within(absEta(), cuts.absEtaMin, cuts.absEtaMax) &&
+                            within(_eta, cuts.etaMin, cuts.etaMax))) &&
          within(track.phi(), cuts.phiMin, cuts.phiMax) &&
          within(track.loc0(), cuts.loc0Min, cuts.loc0Max) &&
          within(track.loc1(), cuts.loc1Min, cuts.loc1Max) &&
@@ -452,26 +464,19 @@ inline TrackSelector::TrackSelector(
         "TrackSelector cut / eta bin configuration is inconsistent"};
   }
 
-  m_isUnbinned = false;
   if (m_cfg.nEtaBins() == 1) {
     static const std::vector<double> infVec = {0, inf};
-    bool limitEta = m_cfg.absEtaEdges != infVec;
-    m_isUnbinned = !limitEta;  // single bin, no eta edges given
-
-    const Config& cuts = m_cfg.cutSets[0];
-
-    if (limitEta && (cuts.etaMin != -inf || cuts.etaMax != inf ||
-                     cuts.absEtaMin != 0.0 || cuts.absEtaMax != inf)) {
-      throw std::invalid_argument{
-          "Explicit eta cuts are only valid for single eta bin"};
-    }
+    m_isUnbinned =
+        m_cfg.absEtaEdges == infVec;  // single bin, no eta edges given
   }
 
-  m_noEtaCuts = m_isUnbinned;
-  for (const auto& cuts : m_cfg.cutSets) {
-    if (cuts.etaMin != -inf || cuts.etaMax != inf || cuts.absEtaMin != 0.0 ||
-        cuts.absEtaMax != inf) {
-      m_noEtaCuts = false;
+  if (!m_isUnbinned) {
+    for (const auto& cuts : m_cfg.cutSets) {
+      if (cuts.etaMin != -inf || cuts.etaMax != inf || cuts.absEtaMin != 0.0 ||
+          cuts.absEtaMax != inf) {
+        throw std::invalid_argument{
+            "Explicit eta cuts are only valid for single eta bin"};
+      }
     }
   }
 }
diff --git a/Core/src/Seeding/CMakeLists.txt b/Core/src/Seeding/CMakeLists.txt
new file mode 100644
index 00000000000..770037b1dfd
--- /dev/null
+++ b/Core/src/Seeding/CMakeLists.txt
@@ -0,0 +1 @@
+target_sources(ActsCore PRIVATE EstimateTrackParamsFromSeed.cpp)
diff --git a/Core/src/Seeding/EstimateTrackParamsFromSeed.cpp b/Core/src/Seeding/EstimateTrackParamsFromSeed.cpp
new file mode 100644
index 00000000000..83d950fc3a1
--- /dev/null
+++ b/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;
+}
diff --git a/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp
index a54b10d2e7d..c8a30c7147d 100644
--- a/Examples/Algorithms/TrackFinding/src/TrackParamsEstimationAlgorithm.cpp
+++ b/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);
diff --git a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp b/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp
index 937a5bc1863..cfc62ab9da9 100644
--- a/Examples/Algorithms/TruthTracking/ActsExamples/TruthTracking/ParticleSmearing.cpp
+++ b/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
 
diff --git a/Examples/Framework/src/Framework/Sequencer.cpp b/Examples/Framework/src/Framework/Sequencer.cpp
index 78444c58ac4..f6bb69d44c0 100644
--- a/Examples/Framework/src/Framework/Sequencer.cpp
+++ b/Examples/Framework/src/Framework/Sequencer.cpp
@@ -611,7 +611,7 @@ void Sequencer::fpeReport() const {
     auto merged = std::accumulate(
         fpe.begin(), fpe.end(), Acts::FpeMonitor::Result{},
         [](const auto& lhs, const auto& rhs) { return lhs.merged(rhs); });
-    if (!merged) {
+    if (!merged.hasStackTraces()) {
       // no FPEs to report
       continue;
     }
diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt
index 5d673eb3678..54c81ea9ba1 100644
--- a/Examples/Python/tests/root_file_hashes.txt
+++ b/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
diff --git a/Fatras/include/ActsFatras/Kernel/Simulation.hpp b/Fatras/include/ActsFatras/Kernel/Simulation.hpp
index cbc1115225a..2d8280737f0 100644
--- a/Fatras/include/ActsFatras/Kernel/Simulation.hpp
+++ b/Fatras/include/ActsFatras/Kernel/Simulation.hpp
@@ -245,6 +245,9 @@ struct Simulation {
           continue;
         }
 
+        assert(result->particle.particleId() == initialParticle.particleId() &&
+               "Particle id must not change during simulation");
+
         copyOutputs(result.value(), simulatedParticlesInitial,
                     simulatedParticlesFinal, hits);
         // since physics processes are independent, there can be particle id
@@ -256,6 +259,10 @@ struct Simulation {
       }
     }
 
+    assert(
+        (simulatedParticlesInitial.size() == simulatedParticlesFinal.size()) &&
+        "Inconsistent final sizes of the simulated particle containers");
+
     // the overall function call succeeded, i.e. no fatal errors occurred.
     // yet, there might have been some particle for which the propagation
     // failed. thus, the successful result contains a list of failed particles.
@@ -284,12 +291,13 @@ struct Simulation {
     // initial particle state was already pushed to the container before
     // store final particle state at the end of the simulation
     particlesFinal.push_back(result.particle);
+    std::copy(result.hits.begin(), result.hits.end(), std::back_inserter(hits));
+
     // move generated secondaries that should be simulated to the output
     std::copy_if(
         result.generatedParticles.begin(), result.generatedParticles.end(),
         std::back_inserter(particlesInitial),
         [this](const Particle &particle) { return selectParticle(particle); });
-    std::copy(result.hits.begin(), result.hits.end(), std::back_inserter(hits));
   }
 
   /// Renumber particle ids in the tail of the container.
diff --git a/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp b/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp
index cc6d7200868..acbaaaa713d 100644
--- a/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp
+++ b/Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp
@@ -69,7 +69,17 @@ struct SimulationActor {
   void act(propagator_state_t &state, stepper_t &stepper,
            navigator_t &navigator, result_type &result,
            const Acts::Logger &logger) const {
-    assert(generator && "The generator pointer must be valid");
+    assert(generator != nullptr && "The generator pointer must be valid");
+
+    if (state.stage == Acts::PropagatorStage::prePropagation) {
+      // first step is special: there is no previous state and we need to arm
+      // the decay simulation for all future steps.
+      result.particle =
+          makeParticle(initialParticle, state, stepper, navigator);
+      result.properTimeLimit =
+          decay.generateProperTimeLimit(*generator, initialParticle);
+      return;
+    }
 
     // actors are called once more after the propagation terminated
     if (!result.isAlive) {
@@ -82,28 +92,11 @@ struct SimulationActor {
       return;
     }
 
-    // check if we are still on the start surface and skip if so
-    if ((navigator.startSurface(state.navigation) != nullptr) &&
-        (navigator.startSurface(state.navigation) ==
-         navigator.currentSurface(state.navigation))) {
-      return;
-    }
-
     // update the particle state first. this also computes the proper time which
     // needs the particle state from the previous step for reference. that means
     // this must happen for every step (not just on surface) and before
     // everything, e.g. any interactions that could modify the state.
-    if (std::isnan(result.properTimeLimit)) {
-      // first step is special: there is no previous state and we need to arm
-      // the decay simulation for all future steps.
-      result.particle =
-          makeParticle(initialParticle, state, stepper, navigator);
-      result.properTimeLimit =
-          decay.generateProperTimeLimit(*generator, initialParticle);
-    } else {
-      result.particle =
-          makeParticle(result.particle, state, stepper, navigator);
-    }
+    result.particle = makeParticle(result.particle, state, stepper, navigator);
 
     // decay check. needs to happen at every step, not just on surfaces.
     if (std::isfinite(result.properTimeLimit) &&
diff --git a/Plugins/FpeMonitoring/include/Acts/Plugins/FpeMonitoring/FpeMonitor.hpp b/Plugins/FpeMonitoring/include/Acts/Plugins/FpeMonitoring/FpeMonitor.hpp
index d1b9fabec64..022a07a3426 100644
--- a/Plugins/FpeMonitoring/include/Acts/Plugins/FpeMonitoring/FpeMonitor.hpp
+++ b/Plugins/FpeMonitoring/include/Acts/Plugins/FpeMonitoring/FpeMonitor.hpp
@@ -40,7 +40,7 @@ std::ostream &operator<<(std::ostream &os, FpeType type);
 class FpeMonitor {
  public:
   struct Buffer {
-    Buffer(std::size_t bufferSize)
+    explicit Buffer(std::size_t bufferSize)
         : m_data{std::make_unique<std::byte[]>(bufferSize)},
           m_size{bufferSize} {}
 
@@ -105,12 +105,12 @@ class FpeMonitor {
 
     Result() = default;
 
-    operator bool() const { return !m_stracktraces.empty(); }
+    bool hasStackTraces() const { return !m_stackTraces.empty(); }
 
     void add(Acts::FpeType type, void *stackPtr, std::size_t bufferSize);
 
    private:
-    std::vector<FpeInfo> m_stracktraces;
+    std::vector<FpeInfo> m_stackTraces;
     std::array<unsigned int, 32> m_counts{};
 
     friend FpeMonitor;
diff --git a/Plugins/FpeMonitoring/src/FpeMonitor.cpp b/Plugins/FpeMonitoring/src/FpeMonitor.cpp
index db41dc1f8b2..09dce1f4b0a 100644
--- a/Plugins/FpeMonitoring/src/FpeMonitor.cpp
+++ b/Plugins/FpeMonitoring/src/FpeMonitor.cpp
@@ -61,10 +61,10 @@ FpeMonitor::Result FpeMonitor::Result::merged(const Result &with) const {
     result.m_counts[i] = m_counts[i] + with.m_counts[i];
   }
 
-  std::copy(with.m_stracktraces.begin(), with.m_stracktraces.end(),
-            std::back_inserter(result.m_stracktraces));
-  std::copy(m_stracktraces.begin(), m_stracktraces.end(),
-            std::back_inserter(result.m_stracktraces));
+  std::copy(with.m_stackTraces.begin(), with.m_stackTraces.end(),
+            std::back_inserter(result.m_stackTraces));
+  std::copy(m_stackTraces.begin(), m_stackTraces.end(),
+            std::back_inserter(result.m_stackTraces));
 
   result.deduplicate();
 
@@ -76,8 +76,8 @@ void FpeMonitor::Result::merge(const Result &with) {
     m_counts[i] = m_counts[i] + with.m_counts[i];
   }
 
-  std::copy(with.m_stracktraces.begin(), with.m_stracktraces.end(),
-            std::back_inserter(m_stracktraces));
+  std::copy(with.m_stackTraces.begin(), with.m_stackTraces.end(),
+            std::back_inserter(m_stackTraces));
 
   deduplicate();
 }
@@ -87,20 +87,20 @@ void FpeMonitor::Result::add(FpeType type, void *stackPtr,
   auto st = std::make_unique<boost::stacktrace::stacktrace>(
       boost::stacktrace::stacktrace::from_dump(stackPtr, bufferSize));
 
-  auto it = std::ranges::find_if(m_stracktraces, [&](const FpeInfo &el) {
+  auto it = std::ranges::find_if(m_stackTraces, [&](const FpeInfo &el) {
     return areFpesEquivalent({el.type, *el.st}, {type, *st});
   });
 
-  if (it != m_stracktraces.end()) {
+  if (it != m_stackTraces.end()) {
     it->count += 1;
   } else {
-    m_stracktraces.push_back({1, type, std::move(st)});
+    m_stackTraces.push_back({1, type, std::move(st)});
   }
 }
 
 bool FpeMonitor::Result::contains(
     FpeType type, const boost::stacktrace::stacktrace &st) const {
-  return std::ranges::any_of(m_stracktraces, [&](const FpeInfo &el) {
+  return std::ranges::any_of(m_stackTraces, [&](const FpeInfo &el) {
     return areFpesEquivalent({el.type, *el.st}, {type, st});
   });
 }
@@ -128,12 +128,12 @@ unsigned int FpeMonitor::Result::count(FpeType type) const {
 }
 
 unsigned int FpeMonitor::Result::numStackTraces() const {
-  return m_stracktraces.size();
+  return m_stackTraces.size();
 }
 
 const std::vector<FpeMonitor::Result::FpeInfo> &
 FpeMonitor::Result::stackTraces() const {
-  return m_stracktraces;
+  return m_stackTraces;
 }
 
 bool FpeMonitor::Result::encountered(FpeType type) const {
@@ -161,18 +161,18 @@ void FpeMonitor::Result::summary(std::ostream &os, std::size_t depth) const {
 
 void FpeMonitor::Result::deduplicate() {
   std::vector<FpeInfo> copy{};
-  copy = std::move(m_stracktraces);
-  m_stracktraces.clear();
+  copy = std::move(m_stackTraces);
+  m_stackTraces.clear();
 
   for (auto &info : copy) {
-    auto it = std::ranges::find_if(m_stracktraces, [&info](const FpeInfo &el) {
+    auto it = std::ranges::find_if(m_stackTraces, [&info](const FpeInfo &el) {
       return areFpesEquivalent({el.type, *el.st}, {info.type, *info.st});
     });
-    if (it != m_stracktraces.end()) {
+    if (it != m_stackTraces.end()) {
       it->count += info.count;
       continue;
     }
-    m_stracktraces.push_back({info.count, info.type, std::move(info.st)});
+    m_stackTraces.push_back({info.count, info.type, std::move(info.st)});
   }
 }
 
diff --git a/Plugins/Json/include/Acts/Plugins/Json/GridJsonConverter.hpp b/Plugins/Json/include/Acts/Plugins/Json/GridJsonConverter.hpp
index 20065d488bb..3215bd66eae 100644
--- a/Plugins/Json/include/Acts/Plugins/Json/GridJsonConverter.hpp
+++ b/Plugins/Json/include/Acts/Plugins/Json/GridJsonConverter.hpp
@@ -9,6 +9,7 @@
 #pragma once
 
 #include "Acts/Plugins/Json/ActsJson.hpp"
+#include "Acts/Plugins/Json/TrackParametersJsonConverter.hpp"
 #include "Acts/Utilities/AxisFwd.hpp"
 #include "Acts/Utilities/GridAccessHelpers.hpp"
 #include "Acts/Utilities/IAxis.hpp"
@@ -266,15 +267,17 @@ auto fromJson(const nlohmann::json& jGrid,
   if constexpr (GridType::DIM == 1u) {
     for (const auto& jd : jData) {
       std::array<std::size_t, 1u> lbin = jd[0u];
-      value_type values = jd[1u];
-      grid.atLocalBins(lbin) = values;
+      if (!jd[1u].is_null()) {
+        grid.atLocalBins(lbin) = jd[1u].get<value_type>();
+      }
     }
   }
   if constexpr (GridType::DIM == 2u) {
     for (const auto& jd : jData) {
       std::array<std::size_t, 2u> lbin = jd[0u];
-      value_type values = jd[1u];
-      grid.atLocalBins(lbin) = values;
+      if (!jd[1u].is_null()) {
+        grid.atLocalBins(lbin) = jd[1u].get<value_type>();
+      }
     }
   }
   return grid;
diff --git a/Plugins/Json/include/Acts/Plugins/Json/TrackParametersJsonConverter.hpp b/Plugins/Json/include/Acts/Plugins/Json/TrackParametersJsonConverter.hpp
new file mode 100644
index 00000000000..ebf7d5c6054
--- /dev/null
+++ b/Plugins/Json/include/Acts/Plugins/Json/TrackParametersJsonConverter.hpp
@@ -0,0 +1,221 @@
+// 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/EventData/TrackParameters.hpp"
+#include "Acts/Plugins/Json/ActsJson.hpp"
+#include "Acts/Plugins/Json/SurfaceJsonConverter.hpp"
+
+#include <nlohmann/json.hpp>
+
+namespace {
+
+// Alias to bound adl_serializer specialization
+// only to track parameters
+template <class parameters_t>
+concept TrackParameters = Acts::FreeTrackParametersConcept<parameters_t> ||
+                          Acts::BoundTrackParametersConcept<parameters_t>;
+
+// Shorthand for bound track parameters
+template <class parameters_t>
+concept IsGenericBound =
+    std::same_as<parameters_t, Acts::GenericBoundTrackParameters<
+                                   typename parameters_t::ParticleHypothesis>>;
+
+}  // namespace
+
+namespace Acts {
+NLOHMANN_JSON_SERIALIZE_ENUM(Acts::PdgParticle,
+
+                             {{Acts::PdgParticle::eInvalid, "Invalid"},
+                              {Acts::PdgParticle::eElectron, "Electron"},
+                              {Acts::PdgParticle::eAntiElectron,
+                               "AntiElectron"},
+                              {Acts::PdgParticle::ePositron, "Positron"},
+                              {Acts::PdgParticle::eMuon, "Muon"},
+                              {Acts::PdgParticle::eAntiMuon, "AntiMuon"},
+                              {Acts::PdgParticle::eTau, "Tau"},
+                              {Acts::PdgParticle::eAntiTau, "AntiTau"},
+                              {Acts::PdgParticle::eGamma, "Gamma"},
+                              {Acts::PdgParticle::ePionZero, "PionZero"},
+                              {Acts::PdgParticle::ePionPlus, "PionPlus"},
+                              {Acts::PdgParticle::ePionMinus, "PionMinus"},
+                              {Acts::PdgParticle::eKaonPlus, "KaonPlus"},
+                              {Acts::PdgParticle::eKaonMinus, "KaonMinus"},
+                              {Acts::PdgParticle::eNeutron, "Neutron"},
+                              {Acts::PdgParticle::eAntiNeutron, "AntiNeutron"},
+                              {Acts::PdgParticle::eProton, "Proton"},
+                              {Acts::PdgParticle::eAntiProton, "AntiProton"},
+                              {Acts::PdgParticle::eLead, "Lead"}}
+
+)
+}
+
+namespace nlohmann {
+
+/// @brief Serialize a track parameters object to json
+///
+/// nlohmann::json serializer specialized for track parameters
+/// as they are not default constructible. Is able to serialize
+/// either bound or free track parameters given that the constructor
+/// convention is followed.
+///
+/// @tparam parameters_t The track parameters type
+template <TrackParameters parameters_t>
+struct adl_serializer<parameters_t> {
+  /// Covariance matrix type attached to the parameters
+  using CovarianceMatrix = typename parameters_t::CovarianceMatrix;
+
+  /// @brief Serialize track parameters object to json
+  ///
+  /// @param j Json object to write to
+  /// @param t Track parameters object to serialize
+  static void to_json(nlohmann::json& j, const parameters_t& t) {
+    // Serialize parameters
+    // common to all track parameters
+    j["direction"] = t.direction();
+    j["qOverP"] = t.qOverP();
+    j["particleHypothesis"] = t.particleHypothesis().absolutePdg();
+
+    // Covariance is optional
+    j["covariance"];
+    if (t.covariance().has_value()) {
+      // Extract covariance matrix
+      // parameters and serialize
+      auto cov = t.covariance().value();
+      constexpr unsigned int size = cov.rows();
+      std::array<Acts::ActsScalar, size * size> covData{};
+      for (std::size_t n = 0; n < size; ++n) {
+        for (std::size_t m = 0; m < size; ++m) {
+          covData[n * size + m] = cov(n, m);
+        }
+      }
+      j["covariance"] = covData;
+    }
+    // Bound track parameters have
+    // reference surface attached
+    // and position takes a geometry context
+    if constexpr (IsGenericBound<parameters_t>) {
+      Acts::GeometryContext gctx;
+      j["position"] = t.fourPosition(gctx);
+
+      j["referenceSurface"] =
+          Acts::SurfaceJsonConverter::toJson(gctx, t.referenceSurface());
+    } else {
+      j["position"] = t.fourPosition();
+    }
+  }
+
+  /// @brief Deserialize track parameters object from json
+  ///
+  /// @param j Json object to read from
+  /// @return Track parameters object
+  static parameters_t from_json(const nlohmann::json& j) {
+    // Extract common parameters
+    std::array<Acts::ActsScalar, 4> posData = j.at("position");
+    Acts::Vector4 position(posData[0], posData[1], posData[2], posData[3]);
+
+    std::array<Acts::ActsScalar, 3> dirData = j.at("direction");
+    Acts::Vector3 direction(dirData[0], dirData[1], dirData[2]);
+
+    Acts::ActsScalar qOverP = j.at("qOverP");
+    Acts::PdgParticle absPdg = j.at("particleHypothesis");
+
+    // Covariance is optional
+    std::optional<CovarianceMatrix> cov;
+    if (j.at("covariance").is_null()) {
+      cov = std::nullopt;
+    } else {
+      // Extract covariance matrix
+      // parameters and deserialize
+      CovarianceMatrix mat;
+      constexpr unsigned int size = mat.rows();
+      std::array<Acts::ActsScalar, size * size> covData = j.at("covariance");
+      for (std::size_t n = 0; n < size; ++n) {
+        for (std::size_t m = 0; m < size; ++m) {
+          mat(n, m) = covData[n * size + m];
+        }
+      }
+      cov.emplace(std::move(mat));
+    }
+
+    // Create particle hypothesis
+    typename parameters_t::ParticleHypothesis particle(absPdg);
+
+    // Bound track parameters have
+    // reference surface attached
+    // and constructor is hidden
+    // behind a factory method
+    if constexpr (IsGenericBound<parameters_t>) {
+      Acts::GeometryContext gctx;
+      auto referenceSurface =
+          Acts::SurfaceJsonConverter::fromJson(j.at("referenceSurface"));
+
+      auto res = parameters_t::create(referenceSurface, gctx, position,
+                                      direction, qOverP, cov, particle);
+
+      if (!res.ok()) {
+        throw std::invalid_argument("Invalid bound track parameters");
+      }
+      return res.value();
+    } else {
+      return parameters_t(position, direction, qOverP, cov, particle);
+    }
+  }
+};
+
+/// @brief Serialize a shared pointer to track parameters object to json
+///
+/// nlohmann::json serializer specialized for shared pointers to track
+/// parameters as they are not default constructible. Is able to serialize
+/// either bound or free track parameters given that the constructor
+/// convention is followed.
+///
+/// @tparam parameters_t The track parameters type
+template <TrackParameters parameters_t>
+struct adl_serializer<std::shared_ptr<parameters_t>> {
+  using CovarianceMatrix = typename parameters_t::CovarianceMatrix;
+  static void to_json(nlohmann::json& j,
+                      const std::shared_ptr<parameters_t>& t) {
+    if (t == nullptr) {
+      return;
+    }
+    j = *t;
+  }
+
+  static std::shared_ptr<parameters_t> from_json(const nlohmann::json& j) {
+    return std::make_shared<parameters_t>(j.get<parameters_t>());
+  }
+};
+
+/// @brief Serialize a unique pointer to track parameters object to json
+///
+/// nlohmann::json serializer specialized for unique pointers to track
+/// parameters as they are not default constructible. Is able to serialize
+/// either bound or free track parameters given that the constructor
+/// convention is followed.
+///
+/// @tparam parameters_t The track parameters type
+template <TrackParameters parameters_t>
+struct adl_serializer<std::unique_ptr<parameters_t>> {
+  using CovarianceMatrix = typename parameters_t::CovarianceMatrix;
+  static void to_json(nlohmann::json& j,
+                      const std::unique_ptr<parameters_t>& t) {
+    if (t == nullptr) {
+      return;
+    }
+    j = *t;
+  }
+
+  static std::unique_ptr<parameters_t> from_json(const nlohmann::json& j) {
+    return std::make_unique<parameters_t>(j.get<parameters_t>());
+  }
+};
+
+}  // namespace nlohmann
diff --git a/Tests/UnitTests/Core/EventData/BoundTrackParametersTests.cpp b/Tests/UnitTests/Core/EventData/BoundTrackParametersTests.cpp
index a05d273f7a1..9daf24c8a79 100644
--- a/Tests/UnitTests/Core/EventData/BoundTrackParametersTests.cpp
+++ b/Tests/UnitTests/Core/EventData/BoundTrackParametersTests.cpp
@@ -76,6 +76,14 @@ void checkParameters(const BoundTrackParameters& params, double l0, double l1,
                        eps);
   CHECK_CLOSE_OR_SMALL(params.momentum(), p * unitDir, eps, eps);
   BOOST_CHECK_EQUAL(params.charge(), q);
+
+  // reflection
+  BoundTrackParameters reflectedParams = params;
+  reflectedParams.reflectInPlace();
+  CHECK_CLOSE_OR_SMALL(params.reflect().parameters(),
+                       reflectedParams.parameters(), eps, eps);
+  CHECK_CLOSE_OR_SMALL(reflectedParams.reflect().parameters(),
+                       params.parameters(), eps, eps);
 }
 
 void runTest(const std::shared_ptr<const Surface>& surface, double l0,
diff --git a/Tests/UnitTests/Core/EventData/CurvilinearTrackParametersTests.cpp b/Tests/UnitTests/Core/EventData/CurvilinearTrackParametersTests.cpp
index 122d4b075d0..143515148d4 100644
--- a/Tests/UnitTests/Core/EventData/CurvilinearTrackParametersTests.cpp
+++ b/Tests/UnitTests/Core/EventData/CurvilinearTrackParametersTests.cpp
@@ -72,6 +72,15 @@ void checkParameters(const CurvilinearTrackParameters& params, double phi,
   // curvilinear reference surface
   CHECK_CLOSE_OR_SMALL(referenceSurface->center(geoCtx), pos, eps, eps);
   CHECK_CLOSE_OR_SMALL(referenceSurface->normal(geoCtx), unitDir, eps, eps);
+
+  // reflection
+  CurvilinearTrackParameters reflectedParams = params;
+  reflectedParams.reflectInPlace();
+  CHECK_CLOSE_OR_SMALL(params.reflect().parameters(),
+                       reflectedParams.parameters(), eps, eps);
+  CHECK_CLOSE_OR_SMALL(reflectedParams.reflect().parameters(),
+                       params.parameters(), eps, eps);
+
   // TODO verify reference frame
 }
 
diff --git a/Tests/UnitTests/Core/EventData/FreeTrackParametersTests.cpp b/Tests/UnitTests/Core/EventData/FreeTrackParametersTests.cpp
index 71b0dc0a20e..fb7a11f3f68 100644
--- a/Tests/UnitTests/Core/EventData/FreeTrackParametersTests.cpp
+++ b/Tests/UnitTests/Core/EventData/FreeTrackParametersTests.cpp
@@ -67,6 +67,14 @@ void checkParameters(const FreeTrackParameters& params, const Vector4& pos4,
                        eps);
   CHECK_CLOSE_OR_SMALL(params.time(), params.template get<eFreeTime>(), eps,
                        eps);
+
+  // reflection
+  FreeTrackParameters reflectedParams = params;
+  reflectedParams.reflectInPlace();
+  CHECK_CLOSE_OR_SMALL(params.reflect().parameters(),
+                       reflectedParams.parameters(), eps, eps);
+  CHECK_CLOSE_OR_SMALL(reflectedParams.reflect().parameters(),
+                       params.parameters(), eps, eps);
 }
 
 }  // namespace
diff --git a/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp b/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp
index e52d586c30b..3306c5702d6 100644
--- a/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp
+++ b/Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp
@@ -244,7 +244,13 @@ BOOST_AUTO_TEST_CASE(HitsOnEmptySurface) {
   BOOST_CHECK_EQUAL(f.actor.initialParticle.absoluteMomentum(), f.p);
   BOOST_CHECK_EQUAL(f.actor.initialParticle.energy(), f.e);
 
+  // call.actor: pre propagation
+  f.state.stage = Acts::PropagatorStage::prePropagation;
+  f.actor.act(f.state, f.stepper, f.navigator, f.result,
+              Acts::getDummyLogger());
+
   // call.actor: surface selection -> one hit, no material -> no secondary
+  f.state.stage = Acts::PropagatorStage::postStep;
   f.actor.act(f.state, f.stepper, f.navigator, f.result,
               Acts::getDummyLogger());
   BOOST_CHECK(f.result.isAlive);
@@ -271,6 +277,7 @@ BOOST_AUTO_TEST_CASE(HitsOnEmptySurface) {
   BOOST_CHECK_EQUAL(f.state.stepping.p, f.result.particle.absoluteMomentum());
 
   // call.actor again: one more hit, still no secondary
+  f.state.stage = Acts::PropagatorStage::postStep;
   f.actor.act(f.state, f.stepper, f.navigator, f.result,
               Acts::getDummyLogger());
   BOOST_CHECK(f.result.isAlive);
@@ -317,7 +324,13 @@ BOOST_AUTO_TEST_CASE(HitsOnMaterialSurface) {
   BOOST_CHECK_EQUAL(f.actor.initialParticle.absoluteMomentum(), f.p);
   BOOST_CHECK_EQUAL(f.actor.initialParticle.energy(), f.e);
 
+  // call.actor: pre propagation
+  f.state.stage = Acts::PropagatorStage::prePropagation;
+  f.actor.act(f.state, f.stepper, f.navigator, f.result,
+              Acts::getDummyLogger());
+
   // call.actor: surface selection -> one hit, material -> one secondary
+  f.state.stage = Acts::PropagatorStage::postStep;
   f.actor.act(f.state, f.stepper, f.navigator, f.result,
               Acts::getDummyLogger());
   BOOST_CHECK(f.result.isAlive);
@@ -346,6 +359,7 @@ BOOST_AUTO_TEST_CASE(HitsOnMaterialSurface) {
                   tol);
 
   // call.actor again: one more hit, one more secondary
+  f.state.stage = Acts::PropagatorStage::postStep;
   f.actor.act(f.state, f.stepper, f.navigator, f.result,
               Acts::getDummyLogger());
   BOOST_CHECK(f.result.isAlive);
@@ -392,7 +406,13 @@ BOOST_AUTO_TEST_CASE(NoHitsEmptySurface) {
   BOOST_CHECK_EQUAL(f.actor.initialParticle.absoluteMomentum(), f.p);
   BOOST_CHECK_EQUAL(f.actor.initialParticle.energy(), f.e);
 
+  // call.actor: pre propagation
+  f.state.stage = Acts::PropagatorStage::prePropagation;
+  f.actor.act(f.state, f.stepper, f.navigator, f.result,
+              Acts::getDummyLogger());
+
   // call.actor: no surface sel. -> no hit, no material -> no secondary
+  f.state.stage = Acts::PropagatorStage::postStep;
   f.actor.act(f.state, f.stepper, f.navigator, f.result,
               Acts::getDummyLogger());
   BOOST_CHECK(f.result.isAlive);
@@ -419,6 +439,7 @@ BOOST_AUTO_TEST_CASE(NoHitsEmptySurface) {
   BOOST_CHECK_EQUAL(f.state.stepping.p, f.result.particle.absoluteMomentum());
 
   // call.actor again: no hit, still no secondary
+  f.state.stage = Acts::PropagatorStage::postStep;
   f.actor.act(f.state, f.stepper, f.navigator, f.result,
               Acts::getDummyLogger());
   BOOST_CHECK(f.result.isAlive);
@@ -455,7 +476,13 @@ BOOST_AUTO_TEST_CASE(NoHitsEmptySurface) {
 BOOST_AUTO_TEST_CASE(NoHitsMaterialSurface) {
   Fixture<NoSurface> f(125_MeV, makeMaterialSurface());
 
+  // call.actor: pre propagation
+  f.state.stage = Acts::PropagatorStage::prePropagation;
+  f.actor.act(f.state, f.stepper, f.navigator, f.result,
+              Acts::getDummyLogger());
+
   // call.actor: no surface sel. -> no hit, material -> one secondary
+  f.state.stage = Acts::PropagatorStage::postStep;
   f.actor.act(f.state, f.stepper, f.navigator, f.result,
               Acts::getDummyLogger());
   BOOST_CHECK(f.result.isAlive);
@@ -483,6 +510,7 @@ BOOST_AUTO_TEST_CASE(NoHitsMaterialSurface) {
                   tol);
 
   // call.actor again: still no hit, one more secondary
+  f.state.stage = Acts::PropagatorStage::postStep;
   f.actor.act(f.state, f.stepper, f.navigator, f.result,
               Acts::getDummyLogger());
   BOOST_CHECK(f.result.isAlive);
@@ -523,7 +551,13 @@ BOOST_AUTO_TEST_CASE(Decay) {
   // inverse Lorentz factor for proper time dilation: 1/gamma = m/E
   const auto gammaInv = f.m / f.e;
 
+  // call.actor: pre propagation
+  f.state.stage = Acts::PropagatorStage::prePropagation;
+  f.actor.act(f.state, f.stepper, f.navigator, f.result,
+              Acts::getDummyLogger());
+
   // first step w/ defaults leaves particle alive
+  f.state.stage = Acts::PropagatorStage::postStep;
   f.actor.act(f.state, f.stepper, f.navigator, f.result,
               Acts::getDummyLogger());
   BOOST_CHECK(f.result.isAlive);
@@ -536,6 +570,7 @@ BOOST_AUTO_TEST_CASE(Decay) {
   BOOST_CHECK_EQUAL(f.result.particle.properTime(), 0_ns);
 
   // second step w/ defaults increases proper time
+  f.state.stage = Acts::PropagatorStage::postStep;
   f.state.stepping.time += 1_ns;
   f.actor.act(f.state, f.stepper, f.navigator, f.result,
               Acts::getDummyLogger());
@@ -549,6 +584,7 @@ BOOST_AUTO_TEST_CASE(Decay) {
   CHECK_CLOSE_REL(f.result.particle.properTime(), gammaInv * 1_ns, tol);
 
   // third step w/ proper time limit decays the particle
+  f.state.stage = Acts::PropagatorStage::postStep;
   f.state.stepping.time += 1_ns;
   f.result.properTimeLimit = f.result.particle.properTime() + gammaInv * 0.5_ns;
   f.actor.act(f.state, f.stepper, f.navigator, f.result,
diff --git a/Tests/UnitTests/Plugins/Json/CMakeLists.txt b/Tests/UnitTests/Plugins/Json/CMakeLists.txt
index 1f8573fdf50..141f86fb3fe 100644
--- a/Tests/UnitTests/Plugins/Json/CMakeLists.txt
+++ b/Tests/UnitTests/Plugins/Json/CMakeLists.txt
@@ -15,3 +15,4 @@ add_unittest(UtilitiesJsonConverter UtilitiesJsonConverterTests.cpp)
 add_unittest(SurfaceBoundsJsonConverter SurfaceBoundsJsonConverterTests.cpp)
 add_unittest(SurfaceJsonConverter SurfaceJsonConverterTests.cpp)
 add_unittest(VolumeBoundsJsonConverter VolumeBoundsJsonConverterTests.cpp)
+add_unittest(TrackParametersJsonConverter TrackParametersJsonConverterTests.cpp)
diff --git a/Tests/UnitTests/Plugins/Json/TrackParametersJsonConverterTests.cpp b/Tests/UnitTests/Plugins/Json/TrackParametersJsonConverterTests.cpp
new file mode 100644
index 00000000000..566464a1889
--- /dev/null
+++ b/Tests/UnitTests/Plugins/Json/TrackParametersJsonConverterTests.cpp
@@ -0,0 +1,97 @@
+// 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 <boost/test/unit_test.hpp>
+
+#include "Acts/Plugins/Json/TrackParametersJsonConverter.hpp"
+#include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Surfaces/RectangleBounds.hpp"
+#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+
+#include <fstream>
+#include <memory>
+
+#include <nlohmann/json.hpp>
+
+BOOST_AUTO_TEST_SUITE(TrackParametersJsonIO)
+
+BOOST_AUTO_TEST_CASE(TrackParametersJsonIO) {
+  Acts::GeometryContext gctx;
+
+  // Track parameters
+  Acts::Vector4 position(1., 2., 3., 4.);
+  Acts::ActsScalar phi = 0.1;
+  Acts::ActsScalar theta = 0.2;
+  Acts::ActsScalar qOverP = 3.0;
+  Acts::ParticleHypothesis particle = Acts::ParticleHypothesis::electron();
+  Acts::FreeMatrix freeCov = Acts::FreeMatrix::Identity();
+  Acts::BoundMatrix boundCov = Acts::BoundMatrix::Identity();
+
+  auto surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
+      Acts::Transform3::Identity(),
+      std::make_shared<Acts::RectangleBounds>(10., 10.));
+  surface->assignGeometryId(Acts::GeometryIdentifier(1u));
+
+  // Free track parameters conversion
+  Acts::FreeTrackParameters ftp(position, phi, theta, qOverP, freeCov,
+                                particle);
+
+  nlohmann::json ftpJson = ftp;
+
+  Acts::FreeTrackParameters ftpRead = ftpJson;
+
+  BOOST_CHECK_EQUAL(ftp.position(), ftpRead.position());
+  BOOST_CHECK_EQUAL(ftp.direction(), ftpRead.direction());
+  BOOST_CHECK_EQUAL(ftp.qOverP(), ftpRead.qOverP());
+  BOOST_CHECK_EQUAL(ftp.covariance().value(), ftpRead.covariance().value());
+  BOOST_CHECK_EQUAL(ftp.particleHypothesis(), ftpRead.particleHypothesis());
+
+  // Curvilinear track parameters conversion
+  Acts::CurvilinearTrackParameters ctp(position, phi, theta, qOverP, boundCov,
+                                       particle);
+
+  nlohmann::json ctpJson = ctp;
+
+  Acts::CurvilinearTrackParameters ctpRead = ctpJson;
+
+  BOOST_CHECK_EQUAL(ctp.position(), ctpRead.position());
+  BOOST_CHECK_EQUAL(ctp.direction(), ctpRead.direction());
+  BOOST_CHECK_EQUAL(ctp.qOverP(), ctpRead.qOverP());
+  BOOST_CHECK_EQUAL(ctp.covariance().value(), ctpRead.covariance().value());
+  BOOST_CHECK_EQUAL(ctp.particleHypothesis(), ctpRead.particleHypothesis());
+
+  BOOST_CHECK(ctp.referenceSurface().transform(gctx).isApprox(
+      ctpRead.referenceSurface().transform(gctx)));
+  BOOST_CHECK_EQUAL(ctp.referenceSurface().geometryId(),
+                    ctpRead.referenceSurface().geometryId());
+  BOOST_CHECK_EQUAL(ctp.referenceSurface().bounds(),
+                    ctpRead.referenceSurface().bounds());
+
+  // Bound track parameters conversion
+  Acts::BoundVector boundPosition{1., 2., 3., 4., 5., 6.};
+  Acts::BoundTrackParameters btp(surface, boundPosition, boundCov, particle);
+
+  nlohmann::json btpJson = btp;
+
+  Acts::BoundTrackParameters btpRead = btpJson;
+
+  BOOST_CHECK_EQUAL(btp.position(gctx), btpRead.position(gctx));
+  BOOST_CHECK_EQUAL(btp.direction(), btpRead.direction());
+  BOOST_CHECK_EQUAL(btp.qOverP(), btpRead.qOverP());
+  BOOST_CHECK_EQUAL(btp.covariance().value(), btpRead.covariance().value());
+  BOOST_CHECK_EQUAL(btp.particleHypothesis(), btpRead.particleHypothesis());
+
+  BOOST_CHECK(btp.referenceSurface().transform(gctx).isApprox(
+      btpRead.referenceSurface().transform(gctx)));
+  BOOST_CHECK_EQUAL(btp.referenceSurface().geometryId(),
+                    btpRead.referenceSurface().geometryId());
+  BOOST_CHECK_EQUAL(btp.referenceSurface().bounds(),
+                    btpRead.referenceSurface().bounds());
+}
+
+BOOST_AUTO_TEST_SUITE_END()