diff --git a/CMakePresets.json b/CMakePresets.json
index bee5ab514df..2aef0f3fe94 100644
--- a/CMakePresets.json
+++ b/CMakePresets.json
@@ -46,7 +46,7 @@
             "inherits": "common",
             "cacheVariables": {
                 "CMAKE_BUILD_TYPE": "Release",
-                "CMAKE_CXX_FLAGS": "-Werror",
+                "CMAKE_COMPILE_WARNING_AS_ERROR": "ON",
                 "ACTS_FORCE_ASSERTIONS": "ON",
                 "ACTS_ENABLE_LOG_FAILURE_THRESHOLD": "ON",
                 "ACTS_BUILD_BENCHMARKS": "ON",
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/Geometry/CylinderVolumeBounds.hpp b/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
index 1951d96b2da..002af04cd59 100644
--- a/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
@@ -128,7 +128,7 @@ class CylinderVolumeBounds : public VolumeBounds {
   /// Copy Constructor
   ///
   /// @param cylbo is the source cylinder volume bounds for the copy
-  CylinderVolumeBounds(const CylinderVolumeBounds& cylbo) = default;
+  CylinderVolumeBounds(const CylinderVolumeBounds& cylbo);
 
   ~CylinderVolumeBounds() override = default;
   CylinderVolumeBounds& operator=(const CylinderVolumeBounds& cylbo) = default;
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/Seeding/HoughTransformUtils.hpp b/Core/include/Acts/Seeding/HoughTransformUtils.hpp
index a7a84d8908c..4fcece48034 100644
--- a/Core/include/Acts/Seeding/HoughTransformUtils.hpp
+++ b/Core/include/Acts/Seeding/HoughTransformUtils.hpp
@@ -14,6 +14,7 @@
 #include "Acts/Utilities/Logger.hpp"
 #include "Acts/Utilities/Result.hpp"
 
+#include <algorithm>
 #include <array>
 #include <map>
 #include <optional>
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/include/Acts/TrackFitting/BetheHeitlerApprox.hpp b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp
index 653bef7e899..666ebf6aa23 100644
--- a/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp
+++ b/Core/include/Acts/TrackFitting/BetheHeitlerApprox.hpp
@@ -125,6 +125,7 @@ class AtlasBetheHeitlerApprox {
   constexpr static double m_singleGaussianLimit = 0.002;
   double m_lowLimit = 0.10;
   double m_highLimit = 0.20;
+  bool m_clampToRange = false;
 
  public:
   /// Construct the Bethe-Heitler approximation description with two
@@ -138,16 +139,19 @@ class AtlasBetheHeitlerApprox {
   /// @param highTransform whether the high data need to be transformed
   /// @param lowLimit the upper limit for the low data
   /// @param highLimit the upper limit for the high data
+  /// @param clampToRange whether to clamp the input x/x0 to the allowed range
   constexpr AtlasBetheHeitlerApprox(const Data &lowData, const Data &highData,
                                     bool lowTransform, bool highTransform,
                                     double lowLimit = 0.1,
-                                    double highLimit = 0.2)
+                                    double highLimit = 0.2,
+                                    bool clampToRange = false)
       : m_lowData(lowData),
         m_highData(highData),
         m_lowTransform(lowTransform),
         m_highTransform(highTransform),
         m_lowLimit(lowLimit),
-        m_highLimit(highLimit) {}
+        m_highLimit(highLimit),
+        m_clampToRange(clampToRange) {}
 
   /// Returns the number of components the returned mixture will have
   constexpr auto numComponents() const { return NComponents; }
@@ -155,7 +159,13 @@ class AtlasBetheHeitlerApprox {
   /// Checks if an input is valid for the parameterization
   ///
   /// @param x pathlength in terms of the radiation length
-  constexpr bool validXOverX0(ActsScalar x) const { return x < m_highLimit; }
+  constexpr bool validXOverX0(ActsScalar x) const {
+    if (m_clampToRange) {
+      return true;
+    } else {
+      return x < m_highLimit;
+    }
+  }
 
   /// Generates the mixture from the polynomials and reweights them, so
   /// that the sum of all weights is 1
@@ -164,6 +174,11 @@ class AtlasBetheHeitlerApprox {
   auto mixture(ActsScalar x) const {
     using Array =
         boost::container::static_vector<detail::GaussianComponent, NComponents>;
+
+    if (m_clampToRange) {
+      x = std::clamp(x, 0.0, m_highLimit);
+    }
+
     // Build a polynom
     auto poly = [](ActsScalar xx,
                    const std::array<ActsScalar, PolyDegree + 1> &coeffs) {
@@ -238,9 +253,11 @@ class AtlasBetheHeitlerApprox {
   /// the parameterization for high x/x0
   /// @param lowLimit the upper limit for the low x/x0-data
   /// @param highLimit the upper limit for the high x/x0-data
+  /// @param clampToRange forwarded to constructor
   static auto loadFromFiles(const std::string &low_parameters_path,
                             const std::string &high_parameters_path,
-                            double lowLimit = 0.1, double highLimit = 0.2) {
+                            double lowLimit = 0.1, double highLimit = 0.2,
+                            bool clampToRange = false) {
     auto read_file = [](const std::string &filepath) {
       std::ifstream file(filepath);
 
@@ -284,7 +301,8 @@ class AtlasBetheHeitlerApprox {
     const auto [highData, highTransform] = read_file(high_parameters_path);
 
     return AtlasBetheHeitlerApprox(lowData, highData, lowTransform,
-                                   highTransform, lowLimit, highLimit);
+                                   highTransform, lowLimit, highLimit,
+                                   clampToRange);
   }
 };
 
@@ -292,6 +310,7 @@ class AtlasBetheHeitlerApprox {
 /// configuration, that are stored as static data in the source code.
 /// This may not be an optimal configuration, but should allow to run
 /// the GSF without the need to load files
-AtlasBetheHeitlerApprox<6, 5> makeDefaultBetheHeitlerApprox();
+AtlasBetheHeitlerApprox<6, 5> makeDefaultBetheHeitlerApprox(
+    bool clampToRange = false);
 
 }  // namespace Acts
diff --git a/Core/src/Geometry/CylinderVolumeBounds.cpp b/Core/src/Geometry/CylinderVolumeBounds.cpp
index c7bb8192f4a..d9af620ea45 100644
--- a/Core/src/Geometry/CylinderVolumeBounds.cpp
+++ b/Core/src/Geometry/CylinderVolumeBounds.cpp
@@ -309,4 +309,7 @@ void CylinderVolumeBounds::set(
   }
 }
 
+CylinderVolumeBounds::CylinderVolumeBounds(const CylinderVolumeBounds& cylbo) =
+    default;
+
 }  // namespace Acts
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/Core/src/TrackFitting/BetheHeitlerApprox.cpp b/Core/src/TrackFitting/BetheHeitlerApprox.cpp
index b19b9921902..272d6e08c8d 100644
--- a/Core/src/TrackFitting/BetheHeitlerApprox.cpp
+++ b/Core/src/TrackFitting/BetheHeitlerApprox.cpp
@@ -8,7 +8,8 @@
 
 #include "Acts/TrackFitting/BetheHeitlerApprox.hpp"
 
-Acts::AtlasBetheHeitlerApprox<6, 5> Acts::makeDefaultBetheHeitlerApprox() {
+Acts::AtlasBetheHeitlerApprox<6, 5> Acts::makeDefaultBetheHeitlerApprox(
+    bool clampToRange) {
   // Tracking/TrkFitter/TrkGaussianSumFilterUtils/Data/BetheHeitler_cdf_nC6_O5.par
   // clang-format off
   constexpr static AtlasBetheHeitlerApprox<6, 5>::Data cdf_cmps6_order5_data = {{
@@ -51,6 +52,7 @@ Acts::AtlasBetheHeitlerApprox<6, 5> Acts::makeDefaultBetheHeitlerApprox() {
   }};
   // clang-format on
 
-  return AtlasBetheHeitlerApprox<6, 5>(
-      cdf_cmps6_order5_data, cdf_cmps6_order5_data, true, true, 0.2, 0.2);
+  return AtlasBetheHeitlerApprox<6, 5>(cdf_cmps6_order5_data,
+                                       cdf_cmps6_order5_data, true, true, 0.2,
+                                       0.2, clampToRange);
 }
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 100f55cf93c..f6bb69d44c0 100644
--- a/Examples/Framework/src/Framework/Sequencer.cpp
+++ b/Examples/Framework/src/Framework/Sequencer.cpp
@@ -126,6 +126,13 @@ Sequencer::Sequencer(const Sequencer::Config& cfg)
         "ACTS_SEQUENCER_DISABLE_FPEMON");
     m_cfg.trackFpes = false;
   }
+
+  if (m_cfg.trackFpes && !m_cfg.fpeMasks.empty() &&
+      !Acts::FpeMonitor::canSymbolize()) {
+    ACTS_ERROR("FPE monitoring is enabled but symbolization is not available");
+    throw std::runtime_error(
+        "FPE monitoring is enabled but symbolization is not available");
+  }
 }
 
 void Sequencer::addContextDecorator(
@@ -604,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/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py
index 3bb898ca998..371345a17b5 100644
--- a/Examples/Python/python/acts/examples/reconstruction.py
+++ b/Examples/Python/python/acts/examples/reconstruction.py
@@ -1375,8 +1375,13 @@ def addTruthTrackingGsf(
 ) -> None:
     customLogLevel = acts.examples.defaultLogging(s, logLevel)
 
+    # NOTE we specify clampToRange as True to silence warnings in the test about
+    # queries to the loss distribution outside the specified range, since no dedicated
+    # approximation for the ODD is done yet.
+    bha = acts.examples.AtlasBetheHeitlerApprox.makeDefault(clampToRange=True)
+
     gsfOptions = {
-        "betheHeitlerApprox": acts.examples.AtlasBetheHeitlerApprox.makeDefault(),
+        "betheHeitlerApprox": bha,
         "maxComponents": 12,
         "componentMergeMethod": acts.examples.ComponentMergeMethod.maxWeight,
         "mixtureReductionAlgorithm": acts.examples.MixtureReductionAlgorithm.KLDistance,
diff --git a/Examples/Python/src/TrackFitting.cpp b/Examples/Python/src/TrackFitting.cpp
index 20c6010c8ee..c8dc72f9eef 100644
--- a/Examples/Python/src/TrackFitting.cpp
+++ b/Examples/Python/src/TrackFitting.cpp
@@ -39,6 +39,7 @@ namespace py = pybind11;
 
 using namespace ActsExamples;
 using namespace Acts;
+using namespace py::literals;
 
 namespace Acts::Python {
 
@@ -106,12 +107,17 @@ void addTrackFitting(Context& ctx) {
         .value("KLDistance", MixtureReductionAlgorithm::KLDistance);
 
     py::class_<ActsExamples::BetheHeitlerApprox>(mex, "AtlasBetheHeitlerApprox")
-        .def_static("loadFromFiles",
-                    &ActsExamples::BetheHeitlerApprox::loadFromFiles,
-                    py::arg("lowParametersPath"), py::arg("highParametersPath"),
-                    py::arg("lowLimit") = 0.1, py::arg("highLimit") = 0.2)
-        .def_static("makeDefault",
-                    []() { return Acts::makeDefaultBetheHeitlerApprox(); });
+        .def_static(
+            "loadFromFiles", &ActsExamples::BetheHeitlerApprox::loadFromFiles,
+            "lowParametersPath"_a, "highParametersPath"_a, "lowLimit"_a = 0.1,
+            "highLimit"_a = 0.2, "clampToRange"_a = false)
+        .def_static(
+            "makeDefault",
+            [](bool clampToRange) {
+              return Acts::makeDefaultBetheHeitlerApprox(clampToRange);
+            },
+            "clampToRange"_a = false);
+
     mex.def(
         "makeGsfFitterFunction",
         [](std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt
index 6ddb1a405cd..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
@@ -74,10 +74,10 @@ test_exatrkx[cpu-torch]__performance_track_finding.root: 36b3045589c4c17c038dbc8
 test_exatrkx[gpu-onnx]__performance_track_finding.root: 9090de10ffb1489d3f1993e2a3081a3038227e3e5c453e98a9a4f33ea3d6d817
 test_exatrkx[gpu-torch]__performance_track_finding.root: 36b3045589c4c17c038dbc87943366f4af4440f7eea6887afb763871ac149b05
 test_ML_Ambiguity_Solver__performance_ambiML.root: 284ff5c3a08c0b810938e4ac2f8ba8fe2babb17d4c202b624ed69fff731a9006
-test_refitting[odd]__trackstates_gsf_refit.root: a61fe2d80d5d10d3b3505000e8abb589d88302bf5f54b625948793418f2a7fe8
-test_refitting[odd]__tracksummary_gsf_refit.root: 082789fc1a85e578b3cf9a750723d2dcee01b5019e871c6a63e0b271f4e907b1
-test_refitting[generic]__trackstates_gsf_refit.root: 9fa7af9eff12081504c0d648f6debae64f6996e0ca610cf58187d23aa5a13251
-test_refitting[generic]__tracksummary_gsf_refit.root: 35b5ac6f208cae093fff94038d217a2e9915a5ce075da2a95718dda696f2d4a2
+test_refitting[odd]__trackstates_gsf_refit.root: e297749dc1e7eda3b8dea13defa0499986c584740d93e723a901b498b8e90c71
+test_refitting[odd]__tracksummary_gsf_refit.root: d5085882e45a0b699194dff9f40a36e9291227bf65f9aaaf9087f9242ef5ae22
+test_refitting[generic]__trackstates_gsf_refit.root: 4424fdf2f27575db825c1a59f8e53a1595946211cbd5b2c8d3a2f71cdcc77ae9
+test_refitting[generic]__tracksummary_gsf_refit.root: 562deecee4cfb97ceee72eff53d63da079e3249fb62d6bcd556e6f27d495dfd9
 test_truth_tracking_kalman[generic-False-0.0]__trackstates_kf.root: 9f77962b92037cb760b1629a602b1dae61f45e659c45d9a87baa784f6190960e
 test_truth_tracking_kalman[generic-False-0.0]__tracksummary_kf.root: 562deecee4cfb97ceee72eff53d63da079e3249fb62d6bcd556e6f27d495dfd9
 test_truth_tracking_kalman[generic-False-1000.0]__trackstates_kf.root: 56a1bd989b9c1316b9098c65fa75df9e6683e62e35ae68d8f72d27220be0fd7d
diff --git a/Examples/Python/tests/test_examples.py b/Examples/Python/tests/test_examples.py
index 0e34a8a3c1c..d5bb593ca9f 100644
--- a/Examples/Python/tests/test_examples.py
+++ b/Examples/Python/tests/test_examples.py
@@ -689,7 +689,7 @@ def test_refitting(tmp_path, detector_config, assert_root_hash):
     field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T))
 
     seq = Sequencer(
-        events=3,
+        events=10,
         numThreads=1,
     )
 
diff --git a/Examples/Python/tests/test_reader.py b/Examples/Python/tests/test_reader.py
index 3d458cb8009..b83b1317cc8 100644
--- a/Examples/Python/tests/test_reader.py
+++ b/Examples/Python/tests/test_reader.py
@@ -290,8 +290,11 @@ def test_edm4hep_simhit_particle_reader(tmp_path):
     tmp_file = str(tmp_path / "output_edm4hep.root")
     odd_xml_file = str(getOpenDataDetectorDirectory() / "xml" / "OpenDataDetector.xml")
 
-    with multiprocessing.get_context("spawn").Pool() as pool:
-        pool.apply(generate_input_test_edm4hep_simhit_reader, (odd_xml_file, tmp_file))
+    p = multiprocessing.Process(
+        target=generate_input_test_edm4hep_simhit_reader, args=(odd_xml_file, tmp_file)
+    )
+    p.start()
+    p.join()
 
     assert os.path.exists(tmp_file)
 
diff --git a/Examples/Scripts/Python/truth_tracking_gsf_refitting.py b/Examples/Scripts/Python/truth_tracking_gsf_refitting.py
index 797bfc59569..e8b41c9aae8 100755
--- a/Examples/Scripts/Python/truth_tracking_gsf_refitting.py
+++ b/Examples/Scripts/Python/truth_tracking_gsf_refitting.py
@@ -25,8 +25,13 @@ def runRefittingGsf(
         s=s,
     )
 
+    # NOTE we specify clampToRange as True to silence warnings in the test about
+    # queries to the loss distribution outside the specified range, since no dedicated
+    # approximation for the ODD is done yet.
+    bha = acts.examples.AtlasBetheHeitlerApprox.makeDefault(clampToRange=True)
+
     gsfOptions = {
-        "betheHeitlerApprox": acts.examples.AtlasBetheHeitlerApprox.makeDefault(),
+        "betheHeitlerApprox": bha,
         "maxComponents": 12,
         "componentMergeMethod": acts.examples.ComponentMergeMethod.maxWeight,
         "mixtureReductionAlgorithm": acts.examples.MixtureReductionAlgorithm.KLDistance,
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/CMakeLists.txt b/Plugins/FpeMonitoring/CMakeLists.txt
index 05a6126e2f3..ccd7c7f4ceb 100644
--- a/Plugins/FpeMonitoring/CMakeLists.txt
+++ b/Plugins/FpeMonitoring/CMakeLists.txt
@@ -24,6 +24,11 @@ else()
 
         set(_backtrace_setup_complete FALSE)
 
+        find_path(
+            boost_stacktrace_include
+            NAMES "boost/stacktrace.hpp"
+            REQUIRED
+        )
         if(Backtrace_FOUND)
             # check if we need to link against bracktrace or not
             set(backtrace_include "")
@@ -44,6 +49,7 @@ else()
                 "${CMAKE_BINARY_DIR}"
                 "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeTmp/backtrace.cpp"
                 LINK_LIBRARIES ${dl_LIBRARY}
+                CMAKE_FLAGS "-DINCLUDE_DIRECTORIES=${boost_stacktrace_include}"
                 COMPILE_DEFINITIONS -DBOOST_STACKTRACE_USE_BACKTRACE
                 OUTPUT_VARIABLE __OUTPUT
             )
@@ -54,9 +60,9 @@ else()
                 message(CHECK_FAIL "no")
 
                 file(GLOB hints "/usr/lib/gcc/*/*/include")
-                find_file(backtrace_header "backtrace.h" HINTS ${hints})
+                find_file(backtrace_header NAMES "backtrace.h" HINTS ${hints})
 
-                if(${backtrace_header} STREQUAL "backtrcae_header-NOTFOUND")
+                if(${backtrace_header} STREQUAL "backtrace_header-NOTFOUND")
                     message(STATUS "Could not find backtrace header file")
                 else()
                     set(backtrace_include
@@ -82,6 +88,8 @@ else()
                         "${CMAKE_BINARY_DIR}"
                         "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeTmp/backtrace.cpp"
                         LINK_LIBRARIES ${dl_LIBRARY}
+                        CMAKE_FLAGS
+                            "-DINCLUDE_DIRECTORIES=${boost_stacktrace_include}"
                         COMPILE_DEFINITIONS
                             -DBOOST_STACKTRACE_USE_BACKTRACE
                             ${backtrace_include}
@@ -111,6 +119,7 @@ else()
                 "${CMAKE_BINARY_DIR}"
                 "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeTmp/backtrace.cpp"
                 LINK_LIBRARIES ${dl_LIBRARY}
+                CMAKE_FLAGS "-DINCLUDE_DIRECTORIES=${boost_stacktrace_include}"
                 COMPILE_DEFINITIONS
                     -DBOOST_STACKTRACE_USE_BACKTRACE
                     ${backtrace_include}
@@ -137,6 +146,8 @@ else()
                     "${CMAKE_BINARY_DIR}"
                     "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/CMakeTmp/backtrace.cpp"
                     LINK_LIBRARIES backtrace ${dl_LIBRARY}
+                    CMAKE_FLAGS
+                        "-DINCLUDE_DIRECTORIES=${boost_stacktrace_include}"
                     COMPILE_DEFINITIONS
                         -DBOOST_STACKTRACE_USE_BACKTRACE
                         ${backtrace_include}
@@ -158,23 +169,6 @@ else()
             endif()
         endif()
 
-        if(NOT _backtrace_setup_complete)
-            message(CHECK_START "Is addr2line available")
-            if(addr2line_EXECUTABLE)
-                list(APPEND _fpe_options -DBOOST_STACKTRACE_USE_ADDR2LINE)
-                list(
-                    APPEND
-                    _fpe_options
-                    -DBOOST_STACKTRACE_ADDR2LINE_LOCATION=${addr2line_EXECUTABLE}
-                )
-                message(CHECK_PASS "yes")
-
-                set(_backtrace_setup_complete TRUE)
-            else()
-                message(CHECK_FAIL "no")
-            endif()
-        endif()
-
         if(NOT _backtrace_setup_complete)
             message(STATUS "Unable to set up stacktrace setup: use noop")
             list(APPEND _fpe_options -BOOST_STACKTRACE_USE_NOOP)
diff --git a/Plugins/FpeMonitoring/include/Acts/Plugins/FpeMonitoring/FpeMonitor.hpp b/Plugins/FpeMonitoring/include/Acts/Plugins/FpeMonitoring/FpeMonitor.hpp
index 99e39da96d0..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;
@@ -131,6 +131,8 @@ class FpeMonitor {
                                         std::size_t depth);
   static std::string getSourceLocation(const boost::stacktrace::frame &frame);
 
+  static bool canSymbolize();
+
  private:
   void enable();
   void disable();
diff --git a/Plugins/FpeMonitoring/src/FpeMonitor.cpp b/Plugins/FpeMonitoring/src/FpeMonitor.cpp
index 8b4feb5f1d2..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)});
   }
 }
 
@@ -332,4 +332,12 @@ std::string FpeMonitor::getSourceLocation(
   return frame.source_file() + ":" + std::to_string(frame.source_line());
 }
 
+bool FpeMonitor::canSymbolize() {
+#if defined(BOOST_STACKTRACE_USE_NOOP)
+  return false;
+#else
+  return true;
+#endif
+}
+
 }  // namespace Acts
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()