diff --git a/CI/physmon/phys_perf_mon.sh b/CI/physmon/phys_perf_mon.sh
index fa3027bb8f6..841fad5c376 100755
--- a/CI/physmon/phys_perf_mon.sh
+++ b/CI/physmon/phys_perf_mon.sh
@@ -20,8 +20,8 @@ shopt -s extglob
 
 
 mode=${1:-all}
-if ! [[ $mode = @(all|kf|gsf|gx2f|fullchains|simulation) ]]; then
-    echo "Usage: $0 <all|kf|gsf|gx2f|fullchains|simulation> (outdir)"
+if ! [[ $mode = @(all|kf|gsf|gx2f|refit_kf|refit_gsf|fullchains|simulation) ]]; then
+    echo "Usage: $0 <all|kf|gsf|gx2f|refit_kf|refit_gsf|fullchains|simulation> (outdir)"
     exit 1
 fi
 
@@ -152,6 +152,12 @@ fi
 if [[ "$mode" == "all" || "$mode" == "gx2f" ]]; then
     run_physmon_gen "Truth Tracking GX2F" "trackfitting_gx2f"
 fi
+if [[ "$mode" == "all" || "$mode" == "refit_kf" ]]; then
+    run_physmon_gen "Truth Tracking KF refit" "trackrefitting_kf"
+fi
+if [[ "$mode" == "all" || "$mode" == "refit_gsf" ]]; then
+    run_physmon_gen "Truth Tracking GSF refit" "trackrefitting_gsf"
+fi
 if [[ "$mode" == "all" || "$mode" == "fullchains" ]]; then
     run_physmon_gen "CKF single muon" "trackfinding_1muon"
     run_physmon_gen "CKF muon 50" "trackfinding_4muon_50vertices"
@@ -411,6 +417,26 @@ if [[ "$mode" == "all" || "$mode" == "gx2f" ]]; then
         --config CI/physmon/config/trackfitting_gx2f.yml
 fi
 
+if [[ "$mode" == "all" || "$mode" == "kf_refit" ]]; then
+    run_histcmp \
+        $outdir/data/trackrefitting_kf/performance_trackrefitting.root \
+        $refdir/trackrefitting_kf/performance_trackrefitting.root \
+        "Truth tracking (KF refit)" \
+        trackrefitting_kf/performance_trackrefitting.html \
+        trackrefitting_kf/performance_trackrefitting_plots \
+        --config CI/physmon/config/trackfitting_kf.yml
+fi
+
+if [[ "$mode" == "all" || "$mode" == "gsf_refit" ]]; then
+    run_histcmp \
+        $outdir/data/trackrefitting_gsf/performance_trackrefitting.root \
+        $refdir/trackrefitting_gsf/performance_trackrefitting.root \
+        "Truth tracking (GSF refit)" \
+        trackrefitting_gsf/performance_trackrefitting.html \
+        trackrefitting_gsf/performance_trackrefitting_plots \
+        --config CI/physmon/config/trackfitting_gsf.yml
+fi
+
 if [[ "$mode" == "all" || "$mode" == "fullchains" ]]; then
     trackfinding "trackfinding | single muon | truth smeared seeding" trackfinding_1muon/truth_smeared
     trackfinding "trackfinding | single muon | truth estimated seeding" trackfinding_1muon/truth_estimated
diff --git a/CI/physmon/reference/trackrefitting_gsf/performance_trackrefitting.root b/CI/physmon/reference/trackrefitting_gsf/performance_trackrefitting.root
new file mode 100644
index 00000000000..d888ad406d8
Binary files /dev/null and b/CI/physmon/reference/trackrefitting_gsf/performance_trackrefitting.root differ
diff --git a/CI/physmon/reference/trackrefitting_kf/performance_trackrefitting.root b/CI/physmon/reference/trackrefitting_kf/performance_trackrefitting.root
new file mode 100644
index 00000000000..aade25ba5a2
Binary files /dev/null and b/CI/physmon/reference/trackrefitting_kf/performance_trackrefitting.root differ
diff --git a/CI/physmon/workflows/physmon_trackrefitting_gsf.py b/CI/physmon/workflows/physmon_trackrefitting_gsf.py
new file mode 100755
index 00000000000..5477121942a
--- /dev/null
+++ b/CI/physmon/workflows/physmon_trackrefitting_gsf.py
@@ -0,0 +1,34 @@
+#!/usr/bin/env python3
+
+import tempfile
+from pathlib import Path
+import shutil
+
+import acts
+from truth_tracking_gsf_refitting import runRefittingGsf
+
+from physmon_common import makeSetup
+
+setup = makeSetup()
+
+with tempfile.TemporaryDirectory() as temp:
+    s = acts.examples.Sequencer(
+        events=10000,
+        numThreads=-1,
+        logLevel=acts.logging.INFO,
+    )
+
+    tp = Path(temp)
+    runRefittingGsf(
+        trackingGeometry=setup.trackingGeometry,
+        field=setup.field,
+        digiConfigFile=setup.digiConfig,
+        outputDir=tp,
+        s=s,
+    )
+
+    s.run()
+
+    perf_file = tp / "performance_gsf_refit.root"
+    assert perf_file.exists(), "Performance file not found"
+    shutil.copy(perf_file, setup.outdir / "performance_trackrefitting.root")
diff --git a/CI/physmon/workflows/physmon_trackrefitting_kf.py b/CI/physmon/workflows/physmon_trackrefitting_kf.py
new file mode 100755
index 00000000000..ea82d26c701
--- /dev/null
+++ b/CI/physmon/workflows/physmon_trackrefitting_kf.py
@@ -0,0 +1,34 @@
+#!/usr/bin/env python3
+
+import tempfile
+from pathlib import Path
+import shutil
+
+import acts
+from truth_tracking_kalman_refitting import runRefittingKf
+
+from physmon_common import makeSetup
+
+setup = makeSetup()
+
+with tempfile.TemporaryDirectory() as temp:
+    s = acts.examples.Sequencer(
+        events=10000,
+        numThreads=-1,
+        logLevel=acts.logging.INFO,
+    )
+
+    tp = Path(temp)
+    runRefittingKf(
+        trackingGeometry=setup.trackingGeometry,
+        field=setup.field,
+        digiConfigFile=setup.digiConfig,
+        outputDir=tp,
+        s=s,
+    )
+
+    s.run()
+
+    perf_file = tp / "performance_kf_refit.root"
+    assert perf_file.exists(), "Performance file not found"
+    shutil.copy(perf_file, setup.outdir / "performance_trackrefitting.root")
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0bbfae3efa8..e8e9ab3f583 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -31,7 +31,6 @@ option(ACTS_FORCE_ASSERTIONS "Force assertions regardless of build type" OFF)
 option(ACTS_USE_SYSTEM_LIBS "Use system libraries by default" OFF)
 # plugins related options
 option(ACTS_USE_SYSTEM_ACTSVG "Use the ActSVG system library" ${ACTS_USE_SYSTEM_LIBS})
-option(ACTS_USE_SYSTEM_GEOMODEL "Use a system-provided GeoModel installation" ${ACTS_USE_SYSTEM_LIBS})
 option(ACTS_USE_SYSTEM_COVFIE "Use a system-provided covfie installation" ${ACTS_USE_SYSTEM_LIBS})
 option(ACTS_USE_SYSTEM_DETRAY "Use a system-provided detray installation" ${ACTS_USE_SYSTEM_LIBS})
 option(ACTS_USE_SYSTEM_TRACCC "Use a system-provided traccc installation" ${ACTS_USE_SYSTEM_LIBS})
@@ -229,7 +228,6 @@ set(_acts_actsvg_version 0.4.50)
 set(_acts_boost_version 1.71.0)
 set(_acts_dd4hep_version 1.21)
 set(_acts_edm4hep_version 0.7)
-set(_acts_geomodel_version 6.3.0)
 set(_acts_eigen3_version 3.4.0)
 set(_acts_podio_version 1.0.1) # will try this first
 set(_acts_podio_fallback_version 0.16) # if not found, will try this one
@@ -374,8 +372,24 @@ if(ACTS_BUILD_PLUGIN_JSON)
     endif()
 endif()
 if(ACTS_BUILD_PLUGIN_GEOMODEL)
-    find_package(GeoModelCore ${_acts_geomodel_version} REQUIRED CONFIG)
-    find_package(GeoModelIO ${_acts_geomodel_version} REQUIRED CONFIG)
+    find_package(GeoModelCore CONFIG)
+    if(NOT GeoModelCore_FOUND)
+        message(
+            FATAL_ERROR
+            "GeoModel not found. Please install GeoModel or set ACTS_BUILD_PLUGIN_GEOMODEL to OFF."
+        )
+    endif()
+
+    set(_gm_ver_min 6.3.0)
+
+    if(GeoModelCore_VERSION VERSION_LESS _gm_ver_min)
+        message(
+            FATAL_ERROR
+            "GeoModel version ${GeoModelCore_VERSION} is insufficient. Please install GeoModel version ${_gm_ver_min} or newer."
+        )
+    endif()
+    # find other GeoModel components of EXACT same version
+    find_package(GeoModelIO ${GeoModelCore_VERSION} REQUIRED EXACT CONFIG)
 endif()
 if(ACTS_BUILD_PLUGIN_TGEO)
     find_package(
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 c92e001fa9f..002af04cd59 100644
--- a/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
@@ -9,6 +9,7 @@
 #pragma once
 
 #include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Geometry/BoundarySurfaceFace.hpp"
 #include "Acts/Geometry/Volume.hpp"
 #include "Acts/Geometry/VolumeBounds.hpp"
 #include "Acts/Utilities/BinningType.hpp"
@@ -80,6 +81,18 @@ class CylinderVolumeBounds : public VolumeBounds {
     eSize
   };
 
+  /// Enum describing the possible faces of a cylinder volume
+  /// @note These values are synchronized with the BoundarySurfaceFace enum.
+  ///       Once Gen1 is removed, this can be changed.
+  enum class Face : unsigned int {
+    PositiveDisc = BoundarySurfaceFace::positiveFaceXY,
+    NegativeDisc = BoundarySurfaceFace::negativeFaceXY,
+    OuterCylinder = BoundarySurfaceFace::tubeOuterCover,
+    InnerCylinder = BoundarySurfaceFace::tubeInnerCover,
+    NegativePhiPlane = BoundarySurfaceFace::tubeSectorNegativePhi,
+    PositivePhiPlane = BoundarySurfaceFace::tubeSectorPositivePhi
+  };
+
   CylinderVolumeBounds() = delete;
 
   /// Constructor
@@ -115,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/Geometry/PortalShell.hpp b/Core/include/Acts/Geometry/PortalShell.hpp
index 39e60b8774b..7405a64ad0a 100644
--- a/Core/include/Acts/Geometry/PortalShell.hpp
+++ b/Core/include/Acts/Geometry/PortalShell.hpp
@@ -8,7 +8,7 @@
 
 #pragma once
 
-#include "Acts/Geometry/BoundarySurfaceFace.hpp"
+#include "Acts/Geometry/CylinderVolumeBounds.hpp"
 #include "Acts/Utilities/BinningType.hpp"
 #include "Acts/Utilities/Logger.hpp"
 
@@ -68,19 +68,9 @@ class PortalShellBase {
 /// volumes
 class CylinderPortalShell : public PortalShellBase {
  public:
-  /// Enum describing the possible faces of a cylinder portal shell
-  /// @note These values are synchronized with the BoundarySurfaceFace enum.
-  ///       Once Gen1 is removed, this can be changed.
-  enum class Face : unsigned int {
-    PositiveDisc = BoundarySurfaceFace::positiveFaceXY,
-    NegativeDisc = BoundarySurfaceFace::negativeFaceXY,
-    OuterCylinder = BoundarySurfaceFace::tubeOuterCover,
-    InnerCylinder = BoundarySurfaceFace::tubeInnerCover,
-    NegativePhiPlane = BoundarySurfaceFace::tubeSectorNegativePhi,
-    PositivePhiPlane = BoundarySurfaceFace::tubeSectorPositivePhi
-  };
-
-  using enum Face;
+  using Face = CylinderVolumeBounds::Face;
+
+  using enum CylinderVolumeBounds::Face;
 
   /// Retrieve the portal associated to the given face. Can be nullptr if unset.
   /// @param face The face to retrieve the portal for
diff --git a/Core/include/Acts/Propagator/DirectNavigator.hpp b/Core/include/Acts/Propagator/DirectNavigator.hpp
index 8df2afdfe8e..eb9db91e71f 100644
--- a/Core/include/Acts/Propagator/DirectNavigator.hpp
+++ b/Core/include/Acts/Propagator/DirectNavigator.hpp
@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include "Acts/Definitions/Direction.hpp"
 #include "Acts/Geometry/BoundarySurfaceT.hpp"
 #include "Acts/Geometry/Layer.hpp"
 #include "Acts/Geometry/TrackingGeometry.hpp"
@@ -17,6 +18,7 @@
 #include "Acts/Surfaces/BoundaryTolerance.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Intersection.hpp"
+#include "Acts/Utilities/Logger.hpp"
 
 #include <algorithm>
 #include <iterator>
@@ -63,7 +65,7 @@ class DirectNavigator {
     Options options;
 
     /// Index of the next surface to try
-    std::size_t surfaceIndex = 0;
+    int surfaceIndex = 0;
 
     /// Navigation state - external interface: the current surface
     const Surface* currentSurface = nullptr;
@@ -72,6 +74,35 @@ class DirectNavigator {
     bool targetReached = false;
     /// Navigation state - external interface: a break has been detected
     bool navigationBreak = false;
+
+    const Surface* navSurface() const {
+      return options.surfaces.at(surfaceIndex);
+    }
+
+    void nextSurface(Direction direction) {
+      if (direction == Direction::Forward) {
+        ++surfaceIndex;
+      } else {
+        --surfaceIndex;
+      }
+    }
+
+    bool endOfSurfaces() const {
+      return surfaceIndex < 0 ||
+             surfaceIndex >= static_cast<int>(options.surfaces.size());
+    }
+
+    int remainingSurfaces(Direction direction) const {
+      if (direction == Direction::Forward) {
+        return options.surfaces.size() - surfaceIndex;
+      }
+      return surfaceIndex + 1;
+    }
+
+    void resetSurfaceIndex(Direction direction) {
+      surfaceIndex =
+          direction == Direction::Forward ? 0 : options.surfaces.size() - 1;
+    }
   };
 
   explicit DirectNavigator(std::unique_ptr<const Logger> _logger =
@@ -135,17 +166,41 @@ class DirectNavigator {
   void initialize(propagator_state_t& state,
                   const stepper_t& /*stepper*/) const {
     ACTS_VERBOSE("Initialize. Surface sequence for navigation:");
-    for (auto surface : state.navigation.options.surfaces) {
+    for (const Surface* surface : state.navigation.options.surfaces) {
       ACTS_VERBOSE(surface->geometryId()
                    << " - " << surface->center(state.geoContext).transpose());
     }
 
     // We set the current surface to the start surface
     state.navigation.currentSurface = state.navigation.options.startSurface;
-    if (state.navigation.currentSurface) {
+    if (state.navigation.currentSurface != nullptr) {
       ACTS_VERBOSE("Current surface set to start surface "
                    << state.navigation.currentSurface->geometryId());
+    } else {
+      ACTS_VERBOSE("Current surface set to nullptr");
+    }
+
+    // Reset the surface index
+    state.navigation.resetSurfaceIndex(state.options.direction);
+    for (const Surface* surface : state.navigation.options.surfaces) {
+      // make sure we skip over the start surface
+      state.navigation.nextSurface(state.options.direction);
+      if (surface == state.navigation.currentSurface) {
+        break;
+      }
+    }
+    ACTS_VERBOSE("Start surface index set to "
+                 << state.navigation.surfaceIndex);
+    if (state.navigation.endOfSurfaces()) {
+      ACTS_DEBUG(
+          "Did not find the start surface in the sequence. Assuming it is not "
+          "part of the sequence. Trusting the correctness of the input "
+          "sequence. Resetting the surface index.");
+      state.navigation.resetSurfaceIndex(state.options.direction);
     }
+
+    state.navigation.navigationBreak = false;
+    state.navigation.targetReached = false;
   }
 
   /// @brief Navigator pre step call
@@ -157,21 +212,25 @@ class DirectNavigator {
   /// @param [in] stepper Stepper in use
   template <typename propagator_state_t, typename stepper_t>
   void preStep(propagator_state_t& state, const stepper_t& stepper) const {
+    if (state.navigation.navigationBreak) {
+      return;
+    }
+
     ACTS_VERBOSE("pre step");
 
     // Navigator target always resets the current surface
     state.navigation.currentSurface = nullptr;
 
     // Output the position in the sequence
-    ACTS_VERBOSE((state.navigation.options.surfaces.size() -
-                  state.navigation.surfaceIndex)
+    ACTS_VERBOSE(state.navigation.remainingSurfaces(state.options.direction)
                  << " out of " << state.navigation.options.surfaces.size()
                  << " surfaces remain to try.");
 
-    if (state.navigation.surfaceIndex >=
-        state.navigation.options.surfaces.size()) {
+    if (state.navigation.endOfSurfaces()) {
       // Set the navigation break
+      ACTS_VERBOSE("End of surfaces reached, navigation break.");
       state.navigation.navigationBreak = true;
+      stepper.releaseStepSize(state.stepping, ConstrainedStep::actor);
       // If no externally provided target is given, the target is reached
       if (state.navigation.options.targetSurface == nullptr) {
         state.navigation.targetReached = true;
@@ -183,8 +242,7 @@ class DirectNavigator {
 
     // Establish & update the surface status
     // TODO we do not know the intersection index - passing the closer one
-    const auto& surface =
-        *state.navigation.options.surfaces.at(state.navigation.surfaceIndex);
+    const auto& surface = *state.navigation.navSurface();
     const double farLimit = std::numeric_limits<double>::max();
     const auto index =
         chooseIntersection(
@@ -202,7 +260,7 @@ class DirectNavigator {
           "Surface not reachable anymore, switching to next one in "
           "sequence");
       // Move the sequence to the next surface
-      ++state.navigation.surfaceIndex;
+      state.navigation.nextSurface(state.options.direction);
     } else {
       ACTS_VERBOSE("Navigation stepSize set to "
                    << stepper.outputStepSize(state.stepping));
@@ -218,26 +276,27 @@ class DirectNavigator {
   /// @param [in] stepper Stepper in use
   template <typename propagator_state_t, typename stepper_t>
   void postStep(propagator_state_t& state, const stepper_t& stepper) const {
+    if (state.navigation.navigationBreak) {
+      return;
+    }
+
     ACTS_VERBOSE("post step");
 
     // Navigator post step always resets the current surface
     state.navigation.currentSurface = nullptr;
 
     // Output the position in the sequence
-    ACTS_VERBOSE((state.navigation.options.surfaces.size() -
-                  state.navigation.surfaceIndex)
+    ACTS_VERBOSE(state.navigation.remainingSurfaces(state.options.direction)
                  << " out of " << state.navigation.options.surfaces.size()
                  << " surfaces remain to try.");
 
-    if (state.navigation.surfaceIndex >=
-        state.navigation.options.surfaces.size()) {
+    if (state.navigation.endOfSurfaces()) {
       return;
     }
 
     // Establish the surface status
     // TODO we do not know the intersection index - passing the closer one
-    const auto& surface =
-        *state.navigation.options.surfaces.at(state.navigation.surfaceIndex);
+    const auto& surface = *state.navigation.navSurface();
     const double farLimit = std::numeric_limits<double>::max();
     const auto index =
         chooseIntersection(
@@ -252,14 +311,12 @@ class DirectNavigator {
         *m_logger);
     if (surfaceStatus == Intersection3D::Status::onSurface) {
       // Set the current surface
-      state.navigation.currentSurface =
-          state.navigation.options.surfaces.at(state.navigation.surfaceIndex);
+      state.navigation.currentSurface = state.navigation.navSurface();
       ACTS_VERBOSE("Current surface set to  "
                    << state.navigation.currentSurface->geometryId());
       // Move the sequence to the next surface
-      ++state.navigation.surfaceIndex;
-      if (state.navigation.surfaceIndex <
-          state.navigation.options.surfaces.size()) {
+      state.navigation.nextSurface(state.options.direction);
+      if (!state.navigation.endOfSurfaces()) {
         ACTS_VERBOSE("Next surface candidate is  "
                      << state.navigation.options.surfaces
                             .at(state.navigation.surfaceIndex)
diff --git a/Core/include/Acts/Propagator/MultiStepperAborters.hpp b/Core/include/Acts/Propagator/MultiStepperAborters.hpp
index f675200975e..0e28375df94 100644
--- a/Core/include/Acts/Propagator/MultiStepperAborters.hpp
+++ b/Core/include/Acts/Propagator/MultiStepperAborters.hpp
@@ -71,8 +71,8 @@ struct MultiStepperSurfaceReached : public SurfaceReached {
       }
 
       ACTS_VERBOSE(
-          "MultiStepperSurfaceReached aborter | "
-          "Target intersection not found. Maybe next time?");
+          "MultiStepperSurfaceReached aborter | Average distance to target: "
+          << sIntersection.pathLength());
     }
 
     bool reached = true;
@@ -84,6 +84,7 @@ struct MultiStepperSurfaceReached : public SurfaceReached {
 
       if (!SurfaceReached::checkAbort(singleState, singleStepper, navigator,
                                       logger)) {
+        cmp.status() = Acts::Intersection3D::Status::reachable;
         reached = false;
       } else {
         cmp.status() = Acts::Intersection3D::Status::onSurface;
diff --git a/Core/include/Acts/Propagator/Navigator.hpp b/Core/include/Acts/Propagator/Navigator.hpp
index b2bea757749..c2a9427257c 100644
--- a/Core/include/Acts/Propagator/Navigator.hpp
+++ b/Core/include/Acts/Propagator/Navigator.hpp
@@ -205,8 +205,6 @@ class Navigator {
       : m_cfg{std::move(cfg)}, m_logger{std::move(_logger)} {}
 
   State makeState(const Options& options) const {
-    assert(options.startSurface != nullptr && "Start surface must be set");
-
     State state;
     state.options = options;
     state.startSurface = options.startSurface;
diff --git a/Core/include/Acts/Propagator/TryAllNavigator.hpp b/Core/include/Acts/Propagator/TryAllNavigator.hpp
index 9a3a2d78475..83d7cbc64a4 100644
--- a/Core/include/Acts/Propagator/TryAllNavigator.hpp
+++ b/Core/include/Acts/Propagator/TryAllNavigator.hpp
@@ -96,8 +96,6 @@ class TryAllNavigatorBase {
       : m_cfg(std::move(cfg)), m_logger{std::move(_logger)} {}
 
   State makeState(const Options& options) const {
-    assert(options.startSurface != nullptr && "Start surface must be set");
-
     State state;
     state.options = options;
     state.startSurface = options.startSurface;
@@ -587,8 +585,6 @@ class TryAllOverstepNavigator : public TryAllNavigatorBase {
       : TryAllNavigatorBase(std::move(cfg), std::move(logger)) {}
 
   State makeState(const Options& options) const {
-    assert(options.startSurface != nullptr && "Start surface must be set");
-
     State state;
     state.options = options;
     state.startSurface = options.startSurface;
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/SeedFinder.hpp b/Core/include/Acts/Seeding/SeedFinder.hpp
index 120610b36a1..cfe9908574c 100644
--- a/Core/include/Acts/Seeding/SeedFinder.hpp
+++ b/Core/include/Acts/Seeding/SeedFinder.hpp
@@ -17,12 +17,14 @@
 #include "Acts/Seeding/SeedFinderConfig.hpp"
 #include "Acts/Seeding/SeedFinderUtils.hpp"
 #include "Acts/Seeding/SpacePointGrid.hpp"
+#include "Acts/Seeding/detail/UtilityFunctions.hpp"
 
 #include <array>
 #include <limits>
 #include <list>
 #include <map>
 #include <memory>
+#include <ranges>
 #include <set>
 #include <string>
 #include <utility>
@@ -30,6 +32,17 @@
 
 namespace Acts {
 
+template <typename Coll>
+concept GridBinCollection =
+    std::ranges::random_access_range<Coll> &&
+    std::same_as<typename Coll::value_type, std::size_t>;
+
+template <typename Coll, typename external_t, std::size_t N = 3ul>
+concept CollectionStoresSeedsTo = requires(Coll coll, external_t sp) {
+  Acts::detail::pushBackOrInsertAtEnd(coll,
+                                      Acts::Seed<external_t, N>(sp, sp, sp));
+};
+
 enum class SpacePointCandidateType : short { eBottom, eTop };
 
 enum class DetectorMeasurementInfo : short { eDefault, eDetailed };
@@ -44,32 +57,32 @@ class SeedFinder {
  public:
   struct SeedingState {
     // bottom space point
-    std::vector<const external_spacepoint_t*> compatBottomSP;
-    std::vector<const external_spacepoint_t*> compatTopSP;
+    std::vector<const external_spacepoint_t*> compatBottomSP{};
+    std::vector<const external_spacepoint_t*> compatTopSP{};
     // contains parameters required to calculate circle with linear equation
     // ...for bottom-middle
-    std::vector<LinCircle> linCircleBottom;
+    std::vector<LinCircle> linCircleBottom{};
     // ...for middle-top
-    std::vector<LinCircle> linCircleTop;
+    std::vector<LinCircle> linCircleTop{};
 
     // create vectors here to avoid reallocation in each loop
-    std::vector<const external_spacepoint_t*> topSpVec;
-    std::vector<float> curvatures;
-    std::vector<float> impactParameters;
+    std::vector<const external_spacepoint_t*> topSpVec{};
+    std::vector<float> curvatures{};
+    std::vector<float> impactParameters{};
 
     // managing seed candidates for SpM
-    CandidatesForMiddleSp<const external_spacepoint_t> candidates_collector;
+    CandidatesForMiddleSp<const external_spacepoint_t> candidates_collector{};
 
     // managing doublet candidates
     boost::container::small_vector<Acts::Neighbour<grid_t>,
                                    Acts::detail::ipow(3, grid_t::DIM)>
-        bottomNeighbours;
+        bottomNeighbours{};
     boost::container::small_vector<Acts::Neighbour<grid_t>,
                                    Acts::detail::ipow(3, grid_t::DIM)>
-        topNeighbours;
+        topNeighbours{};
 
     // Mutable variables for Space points used in the seeding
-    Acts::SpacePointMutableData spacePointMutableData;
+    Acts::SpacePointMutableData spacePointMutableData{};
   };
 
   /// The only constructor. Requires a config object.
@@ -97,7 +110,9 @@ class SeedFinder {
   /// @param rMiddleSPRange range object containing the minimum and maximum r for middle SP for a certain z bin.
   /// @note Ranges must return pointers.
   /// @note Ranges must be separate objects for each parallel call.
-  template <typename container_t, typename sp_range_t>
+  template <typename container_t, Acts::GridBinCollection sp_range_t>
+    requires Acts::CollectionStoresSeedsTo<container_t, external_spacepoint_t,
+                                           3ul>
   void createSeedsForGroup(const Acts::SeedFinderOptions& options,
                            SeedingState& state, const grid_t& grid,
                            container_t& outputCollection,
diff --git a/Core/include/Acts/Seeding/SeedFinder.ipp b/Core/include/Acts/Seeding/SeedFinder.ipp
index f05e14fa0fb..45b95ec05d5 100644
--- a/Core/include/Acts/Seeding/SeedFinder.ipp
+++ b/Core/include/Acts/Seeding/SeedFinder.ipp
@@ -36,7 +36,9 @@ SeedFinder<external_spacepoint_t, grid_t, platform_t>::SeedFinder(
 }
 
 template <typename external_spacepoint_t, typename grid_t, typename platform_t>
-template <typename container_t, typename sp_range_t>
+template <typename container_t, Acts::GridBinCollection sp_range_t>
+  requires Acts::CollectionStoresSeedsTo<container_t, external_spacepoint_t,
+                                         3ul>
 void SeedFinder<external_spacepoint_t, grid_t, platform_t>::createSeedsForGroup(
     const Acts::SeedFinderOptions& options, SeedingState& state,
     const grid_t& grid, container_t& outputCollection,
@@ -366,8 +368,15 @@ SeedFinder<external_spacepoint_t, grid_t, platform_t>::getCompatibleDoublets(
       const float uT = xNewFrame * iDeltaR2;
       const float vT = yNewFrame * iDeltaR2;
 
-      // interactionPointCut == true we apply this cut first cuts before
-      // coordinate transformation to avoid unnecessary calculations
+      // We check the interaction point by evaluating the minimal distance
+      // between the origin and the straight line connecting the two points in
+      // the doublets. Using a geometric similarity, the Im is given by
+      // yNewFrame * rM / deltaR <= m_config.impactMax
+      // However, we make here an approximation of the impact parameter
+      // which is valid under the assumption yNewFrame / xNewFrame is small
+      // The correct computation would be:
+      // yNewFrame * yNewFrame * rM * rM <= m_config.impactMax *
+      // m_config.impactMax * deltaR2
       if (std::abs(rM * yNewFrame) <= impactMax * xNewFrame) {
         // check if duplet cotTheta is within the region of interest
         // cotTheta is defined as (deltaZ / deltaR) but instead we multiply
diff --git a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp
index 7b020010771..4410f1c125b 100644
--- a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp
+++ b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.hpp
@@ -17,9 +17,18 @@
 
 namespace Acts {
 
+/// Concept to check the provided external space point type
+/// can be used to fill the space point grid
+template <typename external_spacepoint_t>
+concept CylindricalGridElement = requires(external_spacepoint_t sp) {
+  { sp.phi() } -> std::same_as<float>;
+  { sp.z() } -> std::same_as<float>;
+  { sp.radius() } -> std::same_as<float>;
+};
+
 /// Cylindrical Space Point bin is a 2D grid with (phi, z) bins
 /// It stores a vector of internal space points to external space points
-template <typename external_spacepoint_t>
+template <Acts::CylindricalGridElement external_spacepoint_t>
 using CylindricalSpacePointGrid = Acts::Grid<
     std::vector<const external_spacepoint_t*>,
     Acts::Axis<Acts::AxisType::Equidistant, Acts::AxisBoundaryType::Closed>,
@@ -90,6 +99,13 @@ struct CylindricalSpacePointGridConfig {
     config.zMin /= 1_mm;
     config.deltaRMax /= 1_mm;
 
+    for (float& val : config.zBinEdges) {
+      val /= 1_mm;
+    }
+    for (float& val : config.rBinEdges) {
+      val /= 1_mm;
+    }
+
     if (config.phiMin < -std::numbers::pi_v<float> ||
         config.phiMax > std::numbers::pi_v<float>) {
       throw std::runtime_error(
@@ -119,7 +135,7 @@ struct CylindricalSpacePointGridConfig {
 
 struct CylindricalSpacePointGridOptions {
   // magnetic field
-  float bFieldInZ = 0.;
+  float bFieldInZ = 0. * Acts::UnitConstants::T;
   bool isInInternalUnits = false;
   CylindricalSpacePointGridOptions toInternalUnits() const {
     if (isInInternalUnits) {
@@ -152,6 +168,16 @@ class CylindricalSpacePointGridCreator {
       Acts::CylindricalSpacePointGrid<external_spacepoint_t>& grid,
       external_spacepoint_iterator_t spBegin,
       external_spacepoint_iterator_t spEnd);
+
+  template <typename external_spacepoint_t, typename external_collection_t>
+    requires std::ranges::range<external_collection_t> &&
+             std::same_as<typename external_collection_t::value_type,
+                          external_spacepoint_t>
+  static void fillGrid(
+      const Acts::SeedFinderConfig<external_spacepoint_t>& config,
+      const Acts::SeedFinderOptions& options,
+      Acts::CylindricalSpacePointGrid<external_spacepoint_t>& grid,
+      const external_collection_t& collection);
 };
 
 }  // namespace Acts
diff --git a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp
index 062aa7595ef..523d470eae5 100644
--- a/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp
+++ b/Core/include/Acts/Seeding/detail/CylindricalSpacePointGrid.ipp
@@ -214,3 +214,17 @@ void Acts::CylindricalSpacePointGridCreator::fillGrid(
     std::ranges::sort(rbin, {}, [](const auto& rb) { return rb->radius(); });
   }
 }
+
+template <typename external_spacepoint_t, typename external_collection_t>
+  requires std::ranges::range<external_collection_t> &&
+           std::same_as<typename external_collection_t::value_type,
+                        external_spacepoint_t>
+void Acts::CylindricalSpacePointGridCreator::fillGrid(
+    const Acts::SeedFinderConfig<external_spacepoint_t>& config,
+    const Acts::SeedFinderOptions& options,
+    Acts::CylindricalSpacePointGrid<external_spacepoint_t>& grid,
+    const external_collection_t& collection) {
+  Acts::CylindricalSpacePointGridCreator::fillGrid<external_spacepoint_t>(
+      config, options, grid, std::ranges::begin(collection),
+      std::ranges::end(collection));
+}
diff --git a/Core/include/Acts/Seeding/detail/UtilityFunctions.hpp b/Core/include/Acts/Seeding/detail/UtilityFunctions.hpp
index 369906d7f46..9cc61c4fe3f 100644
--- a/Core/include/Acts/Seeding/detail/UtilityFunctions.hpp
+++ b/Core/include/Acts/Seeding/detail/UtilityFunctions.hpp
@@ -30,16 +30,23 @@ concept isCollectionThatSupportsInsert =
     };
 
 // Define some functions
-template <typename value_t>
-void pushBackOrInsertAtEnd(
-    Acts::detail::isCollectionThatSupportsPushBack auto& storage,
-    value_t&& value) {
+template <Acts::detail::isCollectionThatSupportsPushBack storage_t,
+          typename value_t>
+  requires requires(storage_t coll, value_t value) {
+    coll.push_back(value);
+    coll.push_back(std::move(value));
+  }
+void pushBackOrInsertAtEnd(storage_t& storage, value_t&& value) {
   storage.push_back(std::forward<value_t>(value));
 }
 
 template <std::ranges::range storage_t, typename value_t>
   requires(!Acts::detail::isCollectionThatSupportsPushBack<storage_t> &&
-           Acts::detail::isCollectionThatSupportsInsert<storage_t>)
+           Acts::detail::isCollectionThatSupportsInsert<storage_t> &&
+           requires(storage_t coll, value_t value) {
+             coll.insert(std::ranges::end(coll), value);
+             coll.insert(std::ranges::end(coll), std::move(value));
+           })
 void pushBackOrInsertAtEnd(storage_t& storage, value_t&& value) {
   storage.insert(std::ranges::end(storage), std::forward<value_t>(value));
 }
diff --git a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp
index afb74140b8d..05e236e1e29 100644
--- a/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp
+++ b/Core/include/Acts/TrackFinding/CombinatorialKalmanFilter.hpp
@@ -34,6 +34,7 @@
 #include <functional>
 #include <limits>
 #include <memory>
+#include <optional>
 #include <type_traits>
 
 namespace Acts {
@@ -579,8 +580,14 @@ class CombinatorialKalmanFilter {
           ACTS_ERROR("Error in filter: " << res.error());
           result.lastError = res.error();
         }
+
+        if (result.finished) {
+          return;
+        }
       }
 
+      assert(!result.activeBranches.empty() && "No active branches");
+
       const bool isEndOfWorldReached =
           endOfWorldReached.checkAbort(state, stepper, navigator, logger());
       const bool isVolumeConstraintReached = volumeConstraintAborter.checkAbort(
@@ -589,9 +596,8 @@ class CombinatorialKalmanFilter {
           state, stepper, navigator, logger());
       const bool isTargetReached =
           targetReached.checkAbort(state, stepper, navigator, logger());
-      const bool allBranchesStopped = result.activeBranches.empty();
-      if (isEndOfWorldReached || isPathLimitReached || isTargetReached ||
-          allBranchesStopped || isVolumeConstraintReached) {
+      if (isEndOfWorldReached || isVolumeConstraintReached ||
+          isPathLimitReached || isTargetReached) {
         if (isEndOfWorldReached) {
           ACTS_VERBOSE("End of world reached");
         } else if (isVolumeConstraintReached) {
@@ -620,26 +626,12 @@ class CombinatorialKalmanFilter {
           stepper.releaseStepSize(state.stepping, ConstrainedStep::actor);
         }
 
-        if (!allBranchesStopped) {
-          // Record the active branch and remove it from the list
-          storeLastActiveBranch(result);
-          result.activeBranches.pop_back();
-        } else {
-          // This can happen if we stopped all branches in the filter step
-          ACTS_VERBOSE("All branches stopped");
-        }
+        // Record the active branch and remove it from the list
+        storeLastActiveBranch(result);
+        result.activeBranches.pop_back();
 
-        // If no more active branches, done with filtering; Otherwise, reset
-        // propagation state to track state at next active branch
-        if (!result.activeBranches.empty()) {
-          ACTS_VERBOSE("Propagation jumps to branch with tip = "
-                       << result.activeBranches.back().tipIndex());
-          reset(state, stepper, navigator, result);
-        } else {
-          ACTS_VERBOSE("Stop Kalman filtering with "
-                       << result.collectedTracks.size() << " found tracks");
-          result.finished = true;
-        }
+        // Reset propagation state to track state at next active branch
+        reset(state, stepper, navigator, result);
       }
     }
 
@@ -665,9 +657,20 @@ class CombinatorialKalmanFilter {
               typename navigator_t>
     void reset(propagator_state_t& state, const stepper_t& stepper,
                const navigator_t& navigator, result_type& result) const {
+      if (result.activeBranches.empty()) {
+        ACTS_VERBOSE("Stop CKF with " << result.collectedTracks.size()
+                                      << " found tracks");
+        result.finished = true;
+
+        return;
+      }
+
       auto currentBranch = result.activeBranches.back();
       auto currentState = currentBranch.outermostTrackState();
 
+      ACTS_VERBOSE("Propagation jumps to branch with tip = "
+                   << currentBranch.tipIndex());
+
       // Reset the stepping state
       stepper.resetState(state.stepping, currentState.filtered(),
                          currentState.filteredCovariance(),
@@ -713,40 +716,67 @@ class CombinatorialKalmanFilter {
       using PM = TrackStatePropMask;
 
       bool isSensitive = surface->associatedDetectorElement() != nullptr;
-      bool isMaterial = surface->surfaceMaterial() != nullptr;
-
-      std::size_t nBranchesOnSurface = 0;
+      bool hasMaterial = surface->surfaceMaterial() != nullptr;
+      bool isMaterialOnly = hasMaterial && !isSensitive;
+      bool expectMeasurements = isSensitive;
 
-      if (auto [slBegin, slEnd] = m_sourceLinkAccessor(*surface);
-          slBegin != slEnd) {
-        // Screen output message
+      if (isSensitive) {
         ACTS_VERBOSE("Measurement surface " << surface->geometryId()
                                             << " detected.");
+      } else if (isMaterialOnly) {
+        ACTS_VERBOSE("Material surface " << surface->geometryId()
+                                         << " detected.");
+      } else {
+        ACTS_VERBOSE("Passive surface " << surface->geometryId()
+                                        << " detected.");
+        return Result<void>::success();
+      }
 
-        // Transport the covariance to the surface
+      using SourceLinkRange = decltype(m_sourceLinkAccessor(*surface));
+      std::optional<SourceLinkRange> slRange = std::nullopt;
+      bool hasMeasurements = false;
+      if (isSensitive) {
+        slRange = m_sourceLinkAccessor(*surface);
+        hasMeasurements = slRange->first != slRange->second;
+      }
+      bool isHole = isSensitive && !hasMeasurements;
+
+      if (isHole) {
+        ACTS_VERBOSE("Detected hole before measurement selection on surface "
+                     << surface->geometryId());
+      }
+
+      // Transport the covariance to the surface
+      if (isHole || isMaterialOnly) {
+        stepper.transportCovarianceToCurvilinear(state.stepping);
+      } else {
         stepper.transportCovarianceToBound(state.stepping, *surface);
+      }
 
-        // Update state and stepper with pre material effects
-        materialInteractor(surface, state, stepper, navigator,
-                           MaterialUpdateStage::PreUpdate);
+      // Update state and stepper with pre material effects
+      materialInteractor(surface, state, stepper, navigator,
+                         MaterialUpdateStage::PreUpdate);
 
-        // Bind the transported state to the current surface
-        auto boundStateRes =
-            stepper.boundState(state.stepping, *surface, false);
-        if (!boundStateRes.ok()) {
-          return boundStateRes.error();
-        }
-        auto& boundState = *boundStateRes;
-        auto& [boundParams, jacobian, pathLength] = boundState;
-        boundParams.covariance() = state.stepping.cov;
+      // Bind the transported state to the current surface
+      auto boundStateRes = stepper.boundState(state.stepping, *surface, false);
+      if (!boundStateRes.ok()) {
+        return boundStateRes.error();
+      }
+      auto& boundState = *boundStateRes;
+      auto& [boundParams, jacobian, pathLength] = boundState;
+      boundParams.covariance() = state.stepping.cov;
+
+      auto currentBranch = result.activeBranches.back();
+      TrackIndexType prevTip = currentBranch.tipIndex();
 
-        auto currentBranch = result.activeBranches.back();
-        TrackIndexType prevTip = currentBranch.tipIndex();
+      // Create trackstates for all source links (will be filtered later)
+      using TrackStatesResult =
+          Acts::Result<CkfTypes::BranchVector<TrackIndexType>>;
+      TrackStatesResult tsRes = TrackStatesResult::success({});
+      if (hasMeasurements) {
+        auto [slBegin, slEnd] = *slRange;
 
-        // Create trackstates for all source links (will be filtered later)
-        using TrackStatesResult =
-            Acts::Result<CkfTypes::BranchVector<TrackIndexType>>;
-        TrackStatesResult tsRes = trackStateCandidateCreator(
+        tsRes = trackStateCandidateCreator(
             state.geoContext, *calibrationContextPtr, *surface, boundState,
             slBegin, slEnd, prevTip, *result.trackStates,
             result.trackStateCandidates, *result.trackStates, logger());
@@ -755,175 +785,98 @@ class CombinatorialKalmanFilter {
               "Processing of selected track states failed: " << tsRes.error());
           return tsRes.error();
         }
-        const CkfTypes::BranchVector<TrackIndexType>& newTrackStateList =
-            *tsRes;
-
-        if (newTrackStateList.empty()) {
-          ACTS_VERBOSE("Detected hole after measurement selection on surface "
-                       << surface->geometryId());
-
-          // Setting the number of branches on the surface to 1 as the hole
-          // still counts as a branch
-          nBranchesOnSurface = 1;
-
-          auto stateMask = PM::Predicted | PM::Jacobian;
-
-          // Add a hole track state to the multitrajectory
-          TrackIndexType currentTip = addNonSourcelinkState(
-              stateMask, boundState, result, true, prevTip);
-          auto nonSourcelinkState =
-              result.trackStates->getTrackState(currentTip);
-          currentBranch.tipIndex() = currentTip;
-          currentBranch.nHoles()++;
+      }
+      const CkfTypes::BranchVector<TrackIndexType>& newTrackStateList = *tsRes;
+
+      if (!newTrackStateList.empty()) {
+        Result<unsigned int> procRes =
+            processNewTrackStates(state.geoContext, newTrackStateList, result);
+        if (!procRes.ok()) {
+          ACTS_ERROR("Processing of selected track states failed: "
+                     << procRes.error());
+          return procRes.error();
+        }
+        unsigned int nBranchesOnSurface = *procRes;
 
-          BranchStopperResult branchStopperResult =
-              m_extensions.branchStopper(currentBranch, nonSourcelinkState);
+        if (nBranchesOnSurface == 0) {
+          ACTS_VERBOSE("All branches on surface " << surface->geometryId()
+                                                  << " have been stopped");
 
-          // Check the branch
-          if (branchStopperResult == BranchStopperResult::Continue) {
-            // Remembered the active branch and its state
-          } else {
-            // No branch on this surface
-            nBranchesOnSurface = 0;
-            if (branchStopperResult == BranchStopperResult::StopAndKeep) {
-              storeLastActiveBranch(result);
-            }
-            // Remove the branch from list
-            result.activeBranches.pop_back();
-          }
-        } else {
-          Result<unsigned int> procRes = processNewTrackStates(
-              state.geoContext, newTrackStateList, result);
-          if (!procRes.ok()) {
-            ACTS_ERROR("Processing of selected track states failed: "
-                       << procRes.error());
-            return procRes.error();
-          }
-          nBranchesOnSurface = *procRes;
+          reset(state, stepper, navigator, result);
 
-          if (nBranchesOnSurface == 0) {
-            // All branches on the surface have been stopped. Reset happens at
-            // the end of the function
-            ACTS_VERBOSE("All branches on surface " << surface->geometryId()
-                                                    << " have been stopped");
-          } else {
-            // `currentBranch` is invalidated after `processNewTrackStates`
-            currentBranch = result.activeBranches.back();
-            prevTip = currentBranch.tipIndex();
-
-            auto currentState = currentBranch.outermostTrackState();
-
-            if (currentState.typeFlags().test(TrackStateFlag::OutlierFlag)) {
-              // We don't need to update the stepper given an outlier state
-              ACTS_VERBOSE("Outlier state detected on surface "
-                           << surface->geometryId());
-            } else {
-              // If there are measurement track states on this surface
-              ACTS_VERBOSE("Filtering step successful with "
-                           << nBranchesOnSurface << " branches");
-              // Update stepping state using filtered parameters of last track
-              // state on this surface
-              stepper.update(state.stepping,
-                             MultiTrajectoryHelpers::freeFiltered(
-                                 state.options.geoContext, currentState),
-                             currentState.filtered(),
-                             currentState.filteredCovariance(), *surface);
-              ACTS_VERBOSE(
-                  "Stepping state is updated with filtered parameter:");
-              ACTS_VERBOSE("-> " << currentState.filtered().transpose()
-                                 << " of track state with tip = "
-                                 << currentState.index());
-            }
-          }
+          return Result<void>::success();
         }
 
-        // Update state and stepper with post material effects
-        materialInteractor(surface, state, stepper, navigator,
-                           MaterialUpdateStage::PostUpdate);
-      } else if (isSensitive || isMaterial) {
-        ACTS_VERBOSE("Handle " << (isSensitive ? "sensitive" : "passive")
-                               << " surface: " << surface->geometryId()
-                               << " without measurements");
-
-        // No splitting on the surface without source links. Set it to one
-        // first, but could be changed later
-        nBranchesOnSurface = 1;
-
-        auto currentBranch = result.activeBranches.back();
-        TrackIndexType prevTip = currentBranch.tipIndex();
+        // `currentBranch` is invalidated after `processNewTrackStates`
+        currentBranch = result.activeBranches.back();
+        prevTip = currentBranch.tipIndex();
+      } else {
+        if (expectMeasurements) {
+          ACTS_VERBOSE("Detected hole after measurement selection on surface "
+                       << surface->geometryId());
+        }
 
-        // No source links on surface, add either hole or passive material
-        // TrackState. No storage allocation for uncalibrated/calibrated
-        // measurement and filtered parameter
         auto stateMask = PM::Predicted | PM::Jacobian;
 
-        // Transport the covariance to a curvilinear surface
-        stepper.transportCovarianceToCurvilinear(state.stepping);
-
-        // Update state and stepper with pre material effects
-        materialInteractor(surface, state, stepper, navigator,
-                           MaterialUpdateStage::PreUpdate);
-
-        // Transport & bind the state to the current surface
-        auto boundStateRes =
-            stepper.boundState(state.stepping, *surface, false);
-        if (!boundStateRes.ok()) {
-          return boundStateRes.error();
-        }
-        auto& boundState = *boundStateRes;
-        auto& [boundParams, jacobian, pathLength] = boundState;
-        boundParams.covariance() = state.stepping.cov;
-
         // Add a hole or material track state to the multitrajectory
         TrackIndexType currentTip = addNonSourcelinkState(
-            stateMask, boundState, result, isSensitive, prevTip);
-        auto nonSourcelinkState = result.trackStates->getTrackState(currentTip);
+            stateMask, boundState, result, expectMeasurements, prevTip);
         currentBranch.tipIndex() = currentTip;
-
-        if (isSensitive) {
+        auto currentState = currentBranch.outermostTrackState();
+        if (expectMeasurements) {
           currentBranch.nHoles()++;
         }
 
         BranchStopperResult branchStopperResult =
-            m_extensions.branchStopper(currentBranch, nonSourcelinkState);
+            m_extensions.branchStopper(currentBranch, currentState);
 
         // Check the branch
         if (branchStopperResult == BranchStopperResult::Continue) {
           // Remembered the active branch and its state
         } else {
           // No branch on this surface
-          nBranchesOnSurface = 0;
           if (branchStopperResult == BranchStopperResult::StopAndKeep) {
             storeLastActiveBranch(result);
           }
           // Remove the branch from list
           result.activeBranches.pop_back();
-        }
 
-        // Update state and stepper with post material effects
-        materialInteractor(surface, state, stepper, navigator,
-                           MaterialUpdateStage::PostUpdate);
-      } else {
-        // Neither measurement nor material on surface, this branch is still
-        // valid. Count the branch on current surface
-        nBranchesOnSurface = 1;
-      }
+          // Branch on the surface has been stopped - reset
+          ACTS_VERBOSE("Branch on surface " << surface->geometryId()
+                                            << " has been stopped");
 
-      // Reset current branch if there is no branch on current surface
-      if (nBranchesOnSurface == 0) {
-        ACTS_DEBUG("Branch on surface " << surface->geometryId()
-                                        << " is stopped");
-        if (!result.activeBranches.empty()) {
-          ACTS_VERBOSE("Propagation jumps to branch with tip = "
-                       << result.activeBranches.back().tipIndex());
           reset(state, stepper, navigator, result);
-        } else {
-          ACTS_VERBOSE("Stop Kalman filtering with "
-                       << result.collectedTracks.size() << " found tracks");
-          result.finished = true;
+
+          return Result<void>::success();
         }
       }
 
+      auto currentState = currentBranch.outermostTrackState();
+
+      if (currentState.typeFlags().test(TrackStateFlag::OutlierFlag)) {
+        // We don't need to update the stepper given an outlier state
+        ACTS_VERBOSE("Outlier state detected on surface "
+                     << surface->geometryId());
+      } else if (currentState.typeFlags().test(
+                     TrackStateFlag::MeasurementFlag)) {
+        // If there are measurement track states on this surface
+        // Update stepping state using filtered parameters of last track
+        // state on this surface
+        stepper.update(state.stepping,
+                       MultiTrajectoryHelpers::freeFiltered(
+                           state.options.geoContext, currentState),
+                       currentState.filtered(),
+                       currentState.filteredCovariance(), *surface);
+        ACTS_VERBOSE("Stepping state is updated with filtered parameter:");
+        ACTS_VERBOSE("-> " << currentState.filtered().transpose()
+                           << " of track state with tip = "
+                           << currentState.index());
+      }
+
+      // Update state and stepper with post material effects
+      materialInteractor(surface, state, stepper, navigator,
+                         MaterialUpdateStage::PostUpdate);
+
       return Result<void>::success();
     }
 
@@ -947,12 +900,13 @@ class CombinatorialKalmanFilter {
 
       auto rootBranch = result.activeBranches.back();
 
-      // Build the new branches by forking the root branch. Reverse the order to
-      // process the best candidate first
+      // Build the new branches by forking the root branch. Reverse the order
+      // to process the best candidate first
       CkfTypes::BranchVector<TrackProxy> newBranches;
       for (auto it = newTrackStateList.rbegin(); it != newTrackStateList.rend();
            ++it) {
-        // Keep the root branch as the first branch, make a copy for the others
+        // Keep the root branch as the first branch, make a copy for the
+        // others
         auto newBranch = (it == newTrackStateList.rbegin())
                              ? rootBranch
                              : rootBranch.shallowCopy();
diff --git a/Core/include/Acts/TrackFinding/MeasurementSelector.ipp b/Core/include/Acts/TrackFinding/MeasurementSelector.ipp
index 3e32c884585..30d09c002be 100644
--- a/Core/include/Acts/TrackFinding/MeasurementSelector.ipp
+++ b/Core/include/Acts/TrackFinding/MeasurementSelector.ipp
@@ -21,9 +21,9 @@ MeasurementSelector::select(
 
   ACTS_VERBOSE("Invoked MeasurementSelector");
 
-  // Return error if no measurement
+  // Return if no measurement
   if (candidates.empty()) {
-    return CombinatorialKalmanFilterError::MeasurementSelectionFailed;
+    return Result::success(std::pair(candidates.begin(), candidates.end()));
   }
 
   // Get geoID of this surface
@@ -92,20 +92,19 @@ MeasurementSelector::select(
           "No measurement candidate. Return an outlier measurement chi2="
           << minChi2);
       isOutlier = true;
-      return Result::success(std::make_pair(candidates.begin() + minIndex,
-                                            candidates.begin() + minIndex + 1));
+      return Result::success(std::pair(candidates.begin() + minIndex,
+                                       candidates.begin() + minIndex + 1));
     } else {
       ACTS_VERBOSE("No measurement candidate. Return empty chi2=" << minChi2);
-      return Result::success(
-          std::make_pair(candidates.begin(), candidates.begin()));
+      return Result::success(std::pair(candidates.begin(), candidates.begin()));
     }
   }
 
   if (passedCandidates <= 1ul) {
     // return single item range, no sorting necessary
     ACTS_VERBOSE("Returning only 1 element chi2=" << minChi2);
-    return Result::success(std::make_pair(candidates.begin() + minIndex,
-                                          candidates.begin() + minIndex + 1));
+    return Result::success(std::pair(candidates.begin() + minIndex,
+                                     candidates.begin() + minIndex + 1));
   }
 
   std::sort(
@@ -115,7 +114,7 @@ MeasurementSelector::select(
   ACTS_VERBOSE("Number of selected measurements: "
                << passedCandidates << ", max: " << cuts.numMeasurements);
 
-  return Result::success(std::make_pair(
+  return Result::success(std::pair(
       candidates.begin(),
       candidates.begin() + std::min(cuts.numMeasurements, passedCandidates)));
 }
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/include/Acts/TrackFitting/GaussianSumFitter.hpp b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp
index efe21b58a60..f6c990a696d 100644
--- a/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp
+++ b/Core/include/Acts/TrackFitting/GaussianSumFitter.hpp
@@ -123,15 +123,11 @@ struct GaussianSumFitter {
       using Actors = ActorList<GsfActor>;
       using PropagatorOptions = typename propagator_t::template Options<Actors>;
 
-      std::vector<const Surface*> backwardSequence(
-          std::next(sSequence.rbegin()), sSequence.rend());
-      backwardSequence.push_back(opts.referenceSurface);
-
       PropagatorOptions propOptions(opts.geoContext, opts.magFieldContext);
 
       propOptions.setPlainOptions(opts.propagatorPlainOptions);
 
-      propOptions.navigation.surfaces = backwardSequence;
+      propOptions.navigation.surfaces = sSequence;
       propOptions.actorList.template get<GsfActor>()
           .m_cfg.bethe_heitler_approx = &m_betheHeitlerApproximation;
 
diff --git a/Core/src/Detector/detail/BlueprintDrawer.cpp b/Core/src/Detector/detail/BlueprintDrawer.cpp
index bc17c4c5922..b518266dc4d 100644
--- a/Core/src/Detector/detail/BlueprintDrawer.cpp
+++ b/Core/src/Detector/detail/BlueprintDrawer.cpp
@@ -57,28 +57,35 @@ std::string labelStr(
 void Acts::Experimental::detail::BlueprintDrawer::dotStream(
     std::ostream& ss, const Acts::Experimental::Blueprint::Node& node,
     const Options& options) {
+  // Replace the "/" in node names
+  std::string nodeName = node.name;
+  std::replace(nodeName.begin(), nodeName.end(), '/', '_');
+
   // Root / leaf or branch
   if (node.isRoot()) {
     ss << "digraph " << options.graphName << " {" << '\n';
-    ss << node.name << " " << labelStr(options.root, node.name, node.auxiliary)
+    ss << nodeName << " " << labelStr(options.root, nodeName, node.auxiliary)
        << '\n';
-    ss << node.name << " " << shapeStr(options.root) << '\n';
+    ss << nodeName << " " << shapeStr(options.root) << '\n';
 
   } else if (node.isLeaf()) {
-    ss << node.name << " " << labelStr(options.leaf, node.name, node.auxiliary)
+    ss << nodeName << " " << labelStr(options.leaf, nodeName, node.auxiliary)
        << '\n';
-    ss << node.name << " "
+    ss << nodeName << " "
        << ((node.internalsBuilder != nullptr) ? shapeStr(options.leaf)
                                               : shapeStr(options.gap))
        << '\n';
   } else {
-    ss << node.name << " "
-       << labelStr(options.branch, node.name, node.auxiliary) << '\n';
-    ss << node.name << " " << shapeStr(options.branch) << '\n';
+    ss << nodeName << " " << labelStr(options.branch, nodeName, node.auxiliary)
+       << '\n';
+    ss << nodeName << " " << shapeStr(options.branch) << '\n';
   }
   // Recursive for children
   for (const auto& c : node.children) {
-    ss << node.name << " -> " << c->name << ";" << '\n';
+    // Replace the "/" in node names
+    std::string childName = c->name;
+    std::replace(childName.begin(), childName.end(), '/', '_');
+    ss << nodeName << " -> " << childName << ";" << '\n';
     dotStream(ss, *c, options);
   }
 
@@ -86,29 +93,29 @@ void Acts::Experimental::detail::BlueprintDrawer::dotStream(
   Options::Node shape = node.isLeaf() ? options.shape : options.virtualShape;
   std::stringstream bts;
   bts << node.boundsType;
-  ss << node.name + "_shape " << shapeStr(shape) << '\n';
-  ss << node.name + "_shape "
+  ss << nodeName + "_shape " << shapeStr(shape) << '\n';
+  ss << nodeName + "_shape "
      << labelStr(shape, bts.str(),
                  {"t = " + toString(node.transform.translation(), 1),
                   "b = " + toString(node.boundaryValues, 1)})
      << '\n';
-  ss << node.name << " -> " << node.name + "_shape [ arrowhead = \"none\" ];"
+  ss << nodeName << " -> " << nodeName + "_shape [ arrowhead = \"none\" ];"
      << '\n';
 
   // Sub node detection
   if (node.internalsBuilder != nullptr) {
-    ss << node.name + "_int " << shapeStr(options.internals) << '\n';
-    ss << node.name << " -> " << node.name + "_int;" << '\n';
+    ss << nodeName + "_int " << shapeStr(options.internals) << '\n';
+    ss << nodeName << " -> " << nodeName + "_int;" << '\n';
   }
 
   if (node.geoIdGenerator != nullptr) {
-    ss << node.name + "_geoID " << shapeStr(options.geoID) << '\n';
-    ss << node.name << " -> " << node.name + "_geoID;" << '\n';
+    ss << nodeName + "_geoID " << shapeStr(options.geoID) << '\n';
+    ss << nodeName << " -> " << nodeName + "_geoID;" << '\n';
   }
 
   if (node.rootVolumeFinderBuilder != nullptr) {
-    ss << node.name + "_roots " << shapeStr(options.roots) << '\n';
-    ss << node.name << " -> " << node.name + "_roots;" << '\n';
+    ss << nodeName + "_roots " << shapeStr(options.roots) << '\n';
+    ss << nodeName << " -> " << nodeName + "_roots;" << '\n';
   }
 
   if (node.isRoot()) {
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/Geometry/PortalShell.cpp b/Core/src/Geometry/PortalShell.cpp
index 59c83c1529f..b9baf403812 100644
--- a/Core/src/Geometry/PortalShell.cpp
+++ b/Core/src/Geometry/PortalShell.cpp
@@ -367,7 +367,7 @@ std::string CylinderStackPortalShell::label() const {
 
 std::ostream& operator<<(std::ostream& os, CylinderPortalShell::Face face) {
   switch (face) {
-    using enum CylinderPortalShell::Face;
+    using enum CylinderVolumeBounds::Face;
     case PositiveDisc:
       return os << "PositiveDisc";
     case NegativeDisc:
diff --git a/Core/src/Geometry/TrackingVolume.cpp b/Core/src/Geometry/TrackingVolume.cpp
index ad2c939fcbd..e7d069b80c5 100644
--- a/Core/src/Geometry/TrackingVolume.cpp
+++ b/Core/src/Geometry/TrackingVolume.cpp
@@ -347,15 +347,6 @@ void TrackingVolume::closeGeometry(
     std::unordered_map<GeometryIdentifier, const TrackingVolume*>& volumeMap,
     std::size_t& vol, const GeometryIdentifierHook& hook,
     const Logger& logger) {
-  if (!boundarySurfaces().empty() && !portals().empty()) {
-    ACTS_ERROR(
-        "TrackingVolume::closeGeometry: Volume "
-        << volumeName()
-        << " has both boundary surfaces and portals. This is not supported.");
-    throw std::invalid_argument(
-        "Volume has both boundary surfaces and portals");
-  }
-
   if (m_confinedVolumes && !volumes().empty()) {
     ACTS_ERROR(
         "TrackingVolume::closeGeometry: Volume "
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/Alignment/include/ActsExamples/Alignment/AlignmentAlgorithm.hpp b/Examples/Algorithms/Alignment/include/ActsExamples/Alignment/AlignmentAlgorithm.hpp
index fe3442b1a29..db5833e5549 100644
--- a/Examples/Algorithms/Alignment/include/ActsExamples/Alignment/AlignmentAlgorithm.hpp
+++ b/Examples/Algorithms/Alignment/include/ActsExamples/Alignment/AlignmentAlgorithm.hpp
@@ -11,7 +11,6 @@
 #include "Acts/Geometry/GeometryHierarchyMap.hpp"
 #include "Acts/TrackFitting/KalmanFitter.hpp"
 #include "ActsAlignment/Kernel/Alignment.hpp"
-#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 #include "ActsExamples/EventData/ProtoTrack.hpp"
 #include "ActsExamples/EventData/Track.hpp"
@@ -88,8 +87,6 @@ class AlignmentAlgorithm final : public IAlgorithm {
   struct Config {
     /// Input measurements collection.
     std::string inputMeasurements;
-    /// Input source links collection.
-    std::string inputSourceLinks;
     /// Input proto tracks collection, i.e. groups of hit indices.
     std::string inputProtoTracks;
     /// Input initial track parameter estimates for for each proto track.
@@ -133,8 +130,6 @@ class AlignmentAlgorithm final : public IAlgorithm {
 
   ReadDataHandle<MeasurementContainer> m_inputMeasurements{this,
                                                            "InputMeasurements"};
-  ReadDataHandle<IndexSourceLinkContainer> m_inputSourceLinks{
-      this, "InputSourceLinks"};
   ReadDataHandle<TrackParametersContainer> m_inputInitialTrackParameters{
       this, "InputInitialTrackParameters"};
   ReadDataHandle<ProtoTrackContainer> m_inputProtoTracks{this,
diff --git a/Examples/Algorithms/Alignment/src/AlignmentAlgorithm.cpp b/Examples/Algorithms/Alignment/src/AlignmentAlgorithm.cpp
index aa685217f05..2f1f450e23c 100644
--- a/Examples/Algorithms/Alignment/src/AlignmentAlgorithm.cpp
+++ b/Examples/Algorithms/Alignment/src/AlignmentAlgorithm.cpp
@@ -11,6 +11,7 @@
 #include "Acts/Surfaces/PerigeeSurface.hpp"
 #include "Acts/TrackFitting/GainMatrixSmoother.hpp"
 #include "Acts/TrackFitting/GainMatrixUpdater.hpp"
+#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/MeasurementCalibration.hpp"
 #include "ActsExamples/EventData/ProtoTrack.hpp"
 #include "ActsExamples/EventData/Trajectories.hpp"
@@ -23,9 +24,6 @@ ActsExamples::AlignmentAlgorithm::AlignmentAlgorithm(Config cfg,
   if (m_cfg.inputMeasurements.empty()) {
     throw std::invalid_argument("Missing input measurement collection");
   }
-  if (m_cfg.inputSourceLinks.empty()) {
-    throw std::invalid_argument("Missing input source links collection");
-  }
   if (m_cfg.inputProtoTracks.empty()) {
     throw std::invalid_argument("Missing input proto tracks collection");
   }
@@ -39,7 +37,6 @@ ActsExamples::AlignmentAlgorithm::AlignmentAlgorithm(Config cfg,
   }
 
   m_inputMeasurements.initialize(m_cfg.inputMeasurements);
-  m_inputSourceLinks.initialize(m_cfg.inputSourceLinks);
   m_inputProtoTracks.initialize(m_cfg.inputProtoTracks);
   m_inputInitialTrackParameters.initialize(m_cfg.inputInitialTrackParameters);
   m_outputAlignmentParameters.initialize(m_cfg.outputAlignmentParameters);
@@ -49,7 +46,6 @@ ActsExamples::ProcessCode ActsExamples::AlignmentAlgorithm::execute(
     const ActsExamples::AlgorithmContext& ctx) const {
   // Read input data
   const auto& measurements = m_inputMeasurements(ctx);
-  const auto& sourceLinks = m_inputSourceLinks(ctx);
   const auto& protoTracks = m_inputProtoTracks(ctx);
   const auto& initialParameters = m_inputInitialTrackParameters(ctx);
 
@@ -79,14 +75,11 @@ ActsExamples::ProcessCode ActsExamples::AlignmentAlgorithm::execute(
     trackSourceLinks.reserve(protoTrack.size());
 
     // Fill the source links via their indices from the container
-    for (auto hitIndex : protoTrack) {
-      auto sourceLink = sourceLinks.nth(hitIndex);
-      if (sourceLink == sourceLinks.end()) {
-        ACTS_FATAL("Proto track " << itrack << " contains invalid hit index"
-                                  << hitIndex);
-        return ProcessCode::ABORT;
-      }
-      trackSourceLinks.push_back(*sourceLink);
+    for (auto measIndex : protoTrack) {
+      const ConstVariableBoundMeasurementProxy measurement =
+          measurements.getMeasurement(measIndex);
+      IndexSourceLink sourceLink(measurement.geometryId(), measIndex);
+      trackSourceLinks.push_back(sourceLink);
     }
     sourceLinkTrackContainer.push_back(trackSourceLinks);
   }
diff --git a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp
index 560fcc9b3e5..fd8937bd300 100644
--- a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp
+++ b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp
@@ -18,7 +18,6 @@
 #include "ActsExamples/Digitization/SmearingConfig.hpp"
 #include "ActsExamples/EventData/Cluster.hpp"
 #include "ActsExamples/EventData/Index.hpp"
-#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 #include "ActsExamples/EventData/SimHit.hpp"
 #include "ActsExamples/Framework/DataHandle.hpp"
@@ -102,8 +101,6 @@ class DigitizationAlgorithm final : public IAlgorithm {
   ReadDataHandle<SimHitContainer> m_simContainerReadHandle{this,
                                                            "SimHitContainer"};
 
-  WriteDataHandle<IndexSourceLinkContainer> m_sourceLinkWriteHandle{
-      this, "SourceLinks"};
   WriteDataHandle<MeasurementContainer> m_measurementWriteHandle{
       this, "Measurements"};
   WriteDataHandle<CellsMap> m_cellsWriteHandle{this, "Cells"};
diff --git a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationConfig.hpp b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationConfig.hpp
index 5fa97d8142c..6a933f50075 100644
--- a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationConfig.hpp
+++ b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationConfig.hpp
@@ -135,8 +135,6 @@ class DigitizationConfig {
 
   /// Input collection of simulated hits.
   std::string inputSimHits = "simhits";
-  /// Output source links collection.
-  std::string outputSourceLinks = "sourcelinks";
   /// Output measurements collection.
   std::string outputMeasurements = "measurements";
   /// Output cells map (geoID -> collection of cells).
diff --git a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/MeasurementCreation.hpp b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/MeasurementCreation.hpp
index b2c76a9820e..667e0c9f4c3 100644
--- a/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/MeasurementCreation.hpp
+++ b/Examples/Algorithms/Digitization/include/ActsExamples/Digitization/MeasurementCreation.hpp
@@ -10,8 +10,8 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "ActsExamples/EventData/Cluster.hpp"
-#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 
 #include <array>
@@ -20,7 +20,6 @@
 #include <vector>
 
 namespace ActsExamples {
-class IndexSourceLink;
 
 /// Struct to identify digitized parameters
 ///
@@ -35,15 +34,16 @@ struct DigitizedParameters {
 
 /// Helper method for created a measurement from digitized parameters
 ///
+/// @param container The measurement container to insert into
+/// @param geometryId The geometry ID of the measurement surface
 /// @param dParams The digitized parameters of variable size
-/// @param isl The indexed source link for the measurement
 ///
 /// To be used also by the e I/O system
 ///
 /// @return the measurement proxy
 ActsExamples::VariableBoundMeasurementProxy createMeasurement(
-    MeasurementContainer& container, const DigitizedParameters& dParams,
-    const IndexSourceLink& isl) noexcept(false);
+    MeasurementContainer& container, Acts::GeometryIdentifier geometryId,
+    const DigitizedParameters& dParams) noexcept(false);
 
 /// Construct the constituents of a measurement.
 ///
diff --git a/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp b/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp
index d6babfe24d0..72f3905ab57 100644
--- a/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp
+++ b/Examples/Algorithms/Digitization/src/DigitizationAlgorithm.cpp
@@ -18,7 +18,6 @@
 #include "ActsExamples/Digitization/ModuleClusters.hpp"
 #include "ActsExamples/EventData/GeometryContainers.hpp"
 #include "ActsExamples/EventData/Index.hpp"
-#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/SimHit.hpp"
 #include "ActsExamples/Framework/AlgorithmContext.hpp"
 #include "ActsExamples/Utilities/GroupBy.hpp"
@@ -55,9 +54,6 @@ ActsExamples::DigitizationAlgorithm::DigitizationAlgorithm(
   }
 
   if (m_cfg.doClusterization) {
-    if (m_cfg.outputSourceLinks.empty()) {
-      throw std::invalid_argument("Missing source links output collection");
-    }
     if (m_cfg.outputMeasurements.empty()) {
       throw std::invalid_argument("Missing measurements output collection");
     }
@@ -73,7 +69,6 @@ ActsExamples::DigitizationAlgorithm::DigitizationAlgorithm(
           "Missing hit-to-simulated-hits map output collection");
     }
 
-    m_sourceLinkWriteHandle.initialize(m_cfg.outputSourceLinks);
     m_measurementWriteHandle.initialize(m_cfg.outputMeasurements);
     m_clusterWriteHandle.initialize(m_cfg.outputClusters);
     m_measurementParticlesMapWriteHandle.initialize(
@@ -152,12 +147,10 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
 
   // Prepare output containers
   // need list here for stable addresses
-  IndexSourceLinkContainer sourceLinks;
   MeasurementContainer measurements;
   ClusterContainer clusters;
   IndexMultimap<ActsFatras::Barcode> measurementParticlesMap;
   IndexMultimap<Index> measurementSimHitsMap;
-  sourceLinks.reserve(simHits.size());
   measurements.reserve(simHits.size());
   measurementParticlesMap.reserve(simHits.size());
   measurementSimHitsMap.reserve(simHits.size());
@@ -292,15 +285,8 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
               // The measurement container is unordered and the index under
               // which the measurement will be stored is known before adding it.
               Index measurementIdx = measurements.size();
-              IndexSourceLink sourceLink{moduleGeoId, measurementIdx};
-
-              // Add to output containers:
-              // index map and source link container are geometry-ordered.
-              // since the input is also geometry-ordered, new items can
-              // be added at the end.
-              sourceLinks.insert(sourceLinks.end(), sourceLink);
 
-              createMeasurement(measurements, dParameters, sourceLink);
+              createMeasurement(measurements, moduleGeoId, dParameters);
               clusters.emplace_back(std::move(dParameters.cluster));
               // this digitization does hit merging so there can be more than
               // one mapping entry for each digitized hit.
@@ -328,7 +314,6 @@ ActsExamples::ProcessCode ActsExamples::DigitizationAlgorithm::execute(
   }
 
   if (m_cfg.doClusterization) {
-    m_sourceLinkWriteHandle(ctx, std::move(sourceLinks));
     m_measurementWriteHandle(ctx, std::move(measurements));
     m_clusterWriteHandle(ctx, std::move(clusters));
     m_measurementParticlesMapWriteHandle(ctx,
diff --git a/Examples/Algorithms/Digitization/src/MeasurementCreation.cpp b/Examples/Algorithms/Digitization/src/MeasurementCreation.cpp
index c91e95a2299..066901cc0cf 100644
--- a/Examples/Algorithms/Digitization/src/MeasurementCreation.cpp
+++ b/Examples/Algorithms/Digitization/src/MeasurementCreation.cpp
@@ -10,6 +10,7 @@
 
 #include "Acts/EventData/MeasurementHelpers.hpp"
 #include "Acts/EventData/SourceLink.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 
@@ -18,10 +19,8 @@
 #include <utility>
 
 ActsExamples::VariableBoundMeasurementProxy ActsExamples::createMeasurement(
-    MeasurementContainer& container, const DigitizedParameters& dParams,
-    const IndexSourceLink& isl) {
-  Acts::SourceLink sl{isl};
-
+    MeasurementContainer& container, Acts::GeometryIdentifier geometryId,
+    const DigitizedParameters& dParams) {
   if (dParams.indices.size() > 4u) {
     std::string errorMsg = "Invalid/mismatching measurement dimension: " +
                            std::to_string(dParams.indices.size());
@@ -32,6 +31,6 @@ ActsExamples::VariableBoundMeasurementProxy ActsExamples::createMeasurement(
       dParams.indices.size(), [&](auto dim) -> VariableBoundMeasurementProxy {
         auto [indices, par, cov] = measurementConstituents<dim>(dParams);
         return VariableBoundMeasurementProxy{
-            container.emplaceMeasurement<dim>(sl, indices, par, cov)};
+            container.emplaceMeasurement<dim>(geometryId, indices, par, cov)};
       });
 }
diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/GbtsSeedingAlgorithm.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/GbtsSeedingAlgorithm.hpp
index b931c72f979..b422eab9e85 100644
--- a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/GbtsSeedingAlgorithm.hpp
+++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/GbtsSeedingAlgorithm.hpp
@@ -44,8 +44,6 @@ class GbtsSeedingAlgorithm final : public IAlgorithm {
 
     std::vector<Acts::GeometryIdentifier> geometrySelection;
 
-    std::string inputSourceLinks;
-
     std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry;
 
     std::map<std::pair<int, int>, std::pair<int, int>> ActsGbtsMap;
@@ -88,9 +86,6 @@ class GbtsSeedingAlgorithm final : public IAlgorithm {
 
   WriteDataHandle<SimSeedContainer> m_outputSeeds{this, "OutputSeeds"};
 
-  ReadDataHandle<IndexSourceLinkContainer> m_inputSourceLinks{
-      this, "InputSourceLinks"};
-
   ReadDataHandle<ClusterContainer> m_inputClusters{this, "InputClusters"};
 };
 
diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/HoughTransformSeeder.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/HoughTransformSeeder.hpp
index 6ea2400a712..b74047ca6a2 100644
--- a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/HoughTransformSeeder.hpp
+++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/HoughTransformSeeder.hpp
@@ -76,7 +76,6 @@
 #include "Acts/Utilities/Logger.hpp"
 #include "Acts/Utilities/Result.hpp"
 #include "ActsExamples/EventData/Index.hpp"
-#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 #include "ActsExamples/EventData/ProtoTrack.hpp"
 #include "ActsExamples/EventData/SimSpacePoint.hpp"
@@ -168,8 +167,6 @@ class HoughTransformSeeder final : public IAlgorithm {
     std::string outputSeeds;
     /// Output hough track collection.
     std::string outputProtoTracks;
-    /// Input source links collection.
-    std::string inputSourceLinks;
     /// Tracking geometry required to access global-to-local transforms.
     std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry;
     /// For which part of the detector geometry should space points be created.
@@ -281,9 +278,6 @@ class HoughTransformSeeder final : public IAlgorithm {
   ReadDataHandle<MeasurementContainer> m_inputMeasurements{this,
                                                            "InputMeasurements"};
 
-  ReadDataHandle<IndexSourceLinkContainer> m_inputSourceLinks{
-      this, "InputSourceLinks"};
-
   ////////////////////////////////////////////////////////////////////////
   /// Convenience
 
diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SpacePointMaker.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SpacePointMaker.hpp
index 0d5604e5ec4..db4678a6200 100644
--- a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SpacePointMaker.hpp
+++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/SpacePointMaker.hpp
@@ -11,7 +11,6 @@
 #include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "Acts/SpacePointFormation/SpacePointBuilder.hpp"
 #include "Acts/Utilities/Logger.hpp"
-#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 #include "ActsExamples/EventData/SimSpacePoint.hpp"
 #include "ActsExamples/Framework/DataHandle.hpp"
@@ -46,8 +45,6 @@ struct AlgorithmContext;
 class SpacePointMaker final : public IAlgorithm {
  public:
   struct Config {
-    /// Input source links collection.
-    std::string inputSourceLinks;
     /// Input measurements collection.
     std::string inputMeasurements;
     /// Output space points collection.
@@ -86,9 +83,6 @@ class SpacePointMaker final : public IAlgorithm {
 
   Acts::SpacePointBuilder<SimSpacePoint> m_spacePointBuilder;
 
-  ReadDataHandle<IndexSourceLinkContainer> m_inputSourceLinks{
-      this, "InputSourceLinks"};
-
   ReadDataHandle<MeasurementContainer> m_inputMeasurements{this,
                                                            "InputMeasurements"};
 
diff --git a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp
index f63581470e6..deb15750a5f 100644
--- a/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp
+++ b/Examples/Algorithms/TrackFinding/include/ActsExamples/TrackFinding/TrackFindingAlgorithm.hpp
@@ -80,8 +80,6 @@ class TrackFindingAlgorithm final : public IAlgorithm {
   struct Config {
     /// Input measurements collection.
     std::string inputMeasurements;
-    /// Input source links collection.
-    std::string inputSourceLinks;
     /// Input initial track parameter estimates for for each proto track.
     std::string inputInitialTrackParameters;
     /// Input seeds. These are optional and allow for seed deduplication.
@@ -156,9 +154,8 @@ class TrackFindingAlgorithm final : public IAlgorithm {
   const Config& config() const { return m_cfg; }
 
  private:
-  template <typename source_link_accessor_container_t>
-  void computeSharedHits(const source_link_accessor_container_t& sourceLinks,
-                         TrackContainer& tracks) const;
+  void computeSharedHits(TrackContainer& tracks,
+                         const MeasurementContainer& measurements) const;
 
   ActsExamples::ProcessCode finalize() override;
 
@@ -168,8 +165,6 @@ class TrackFindingAlgorithm final : public IAlgorithm {
 
   ReadDataHandle<MeasurementContainer> m_inputMeasurements{this,
                                                            "InputMeasurements"};
-  ReadDataHandle<IndexSourceLinkContainer> m_inputSourceLinks{
-      this, "InputSourceLinks"};
   ReadDataHandle<TrackParametersContainer> m_inputInitialTrackParameters{
       this, "InputInitialTrackParameters"};
   ReadDataHandle<SimSeedContainer> m_inputSeeds{this, "InputSeeds"};
@@ -192,56 +187,4 @@ class TrackFindingAlgorithm final : public IAlgorithm {
       }};
 };
 
-// TODO this is somewhat duplicated in AmbiguityResolutionAlgorithm.cpp
-// TODO we should make a common implementation in the core at some point
-template <typename source_link_accessor_container_t>
-void TrackFindingAlgorithm::computeSharedHits(
-    const source_link_accessor_container_t& sourceLinks,
-    TrackContainer& tracks) const {
-  // Compute shared hits from all the reconstructed tracks
-  // Compute nSharedhits and Update ckf results
-  // hit index -> list of multi traj indexes [traj, meas]
-
-  std::vector<std::size_t> firstTrackOnTheHit(
-      sourceLinks.size(), std::numeric_limits<std::size_t>::max());
-  std::vector<std::size_t> firstStateOnTheHit(
-      sourceLinks.size(), std::numeric_limits<std::size_t>::max());
-
-  for (auto track : tracks) {
-    for (auto state : track.trackStatesReversed()) {
-      if (!state.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) {
-        continue;
-      }
-
-      std::size_t hitIndex = state.getUncalibratedSourceLink()
-                                 .template get<IndexSourceLink>()
-                                 .index();
-
-      // Check if hit not already used
-      if (firstTrackOnTheHit.at(hitIndex) ==
-          std::numeric_limits<std::size_t>::max()) {
-        firstTrackOnTheHit.at(hitIndex) = track.index();
-        firstStateOnTheHit.at(hitIndex) = state.index();
-        continue;
-      }
-
-      // if already used, control if first track state has been marked
-      // as shared
-      int indexFirstTrack = firstTrackOnTheHit.at(hitIndex);
-      int indexFirstState = firstStateOnTheHit.at(hitIndex);
-
-      auto firstState = tracks.getTrack(indexFirstTrack)
-                            .container()
-                            .trackStateContainer()
-                            .getTrackState(indexFirstState);
-      if (!firstState.typeFlags().test(Acts::TrackStateFlag::SharedHitFlag)) {
-        firstState.typeFlags().set(Acts::TrackStateFlag::SharedHitFlag);
-      }
-
-      // Decorate this track state
-      state.typeFlags().set(Acts::TrackStateFlag::SharedHitFlag);
-    }
-  }
-}
-
 }  // namespace ActsExamples
diff --git a/Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp
index 835b027862f..385a3db7ec1 100644
--- a/Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp
+++ b/Examples/Algorithms/TrackFinding/src/GbtsSeedingAlgorithm.cpp
@@ -61,8 +61,6 @@ ActsExamples::GbtsSeedingAlgorithm::GbtsSeedingAlgorithm(
 
   m_outputSeeds.initialize(m_cfg.outputSeeds);
 
-  m_inputSourceLinks.initialize(m_cfg.inputSourceLinks);
-
   m_inputClusters.initialize(m_cfg.inputClusters);
 
   // map
diff --git a/Examples/Algorithms/TrackFinding/src/HoughTransformSeeder.cpp b/Examples/Algorithms/TrackFinding/src/HoughTransformSeeder.cpp
index 02b112fe6b4..c8126773194 100644
--- a/Examples/Algorithms/TrackFinding/src/HoughTransformSeeder.cpp
+++ b/Examples/Algorithms/TrackFinding/src/HoughTransformSeeder.cpp
@@ -76,13 +76,7 @@ ActsExamples::HoughTransformSeeder::HoughTransformSeeder(
         "HoughTransformSeeder: Missing hough track seeds output collection");
   }
 
-  if (m_cfg.inputSourceLinks.empty()) {
-    throw std::invalid_argument(
-        "HoughTransformSeeder: Missing source link input collection");
-  }
-
   m_outputProtoTracks.initialize(m_cfg.outputProtoTracks);
-  m_inputSourceLinks.initialize(m_cfg.inputSourceLinks);
   m_inputMeasurements.initialize(m_cfg.inputMeasurements);
 
   if (!m_cfg.trackingGeometry) {
@@ -503,14 +497,14 @@ void ActsExamples::HoughTransformSeeder::addSpacePoints(
 void ActsExamples::HoughTransformSeeder::addMeasurements(
     const AlgorithmContext& ctx) const {
   const auto& measurements = m_inputMeasurements(ctx);
-  const auto& sourceLinks = m_inputSourceLinks(ctx);
 
   ACTS_DEBUG("Inserting " << measurements.size() << " space points from "
                           << m_cfg.inputMeasurements);
 
   for (Acts::GeometryIdentifier geoId : m_cfg.geometrySelection) {
     // select volume/layer depending on what is set in the geometry id
-    auto range = selectLowestNonZeroGeometryObject(sourceLinks, geoId);
+    auto range =
+        selectLowestNonZeroGeometryObject(measurements.orderedIndices(), geoId);
     // groupByModule only works with geometry containers, not with an
     // arbitrary range. do the equivalent grouping manually
     auto groupedByModule = makeGroupBy(range, detail::GeometryIdGetter());
diff --git a/Examples/Algorithms/TrackFinding/src/SeedingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/SeedingAlgorithm.cpp
index d8bd6265f27..ed5b11e4df2 100644
--- a/Examples/Algorithms/TrackFinding/src/SeedingAlgorithm.cpp
+++ b/Examples/Algorithms/TrackFinding/src/SeedingAlgorithm.cpp
@@ -246,9 +246,8 @@ ActsExamples::ProcessCode ActsExamples::SeedingAlgorithm::execute(
       Acts::CylindricalSpacePointGridCreator::createGrid<value_type>(
           m_cfg.gridConfig, m_cfg.gridOptions);
 
-  Acts::CylindricalSpacePointGridCreator::fillGrid(
-      m_cfg.seedFinderConfig, m_cfg.seedFinderOptions, grid,
-      spContainer.begin(), spContainer.end());
+  Acts::CylindricalSpacePointGridCreator::fillGrid<value_type>(
+      m_cfg.seedFinderConfig, m_cfg.seedFinderOptions, grid, spContainer);
 
   // Compute radius Range
   // we rely on the fact the grid is storing the proxies
diff --git a/Examples/Algorithms/TrackFinding/src/SpacePointMaker.cpp b/Examples/Algorithms/TrackFinding/src/SpacePointMaker.cpp
index e6fed997a02..673ffa8b103 100644
--- a/Examples/Algorithms/TrackFinding/src/SpacePointMaker.cpp
+++ b/Examples/Algorithms/TrackFinding/src/SpacePointMaker.cpp
@@ -32,9 +32,6 @@
 ActsExamples::SpacePointMaker::SpacePointMaker(Config cfg,
                                                Acts::Logging::Level lvl)
     : IAlgorithm("SpacePointMaker", lvl), m_cfg(std::move(cfg)) {
-  if (m_cfg.inputSourceLinks.empty()) {
-    throw std::invalid_argument("Missing source link input collection");
-  }
   if (m_cfg.inputMeasurements.empty()) {
     throw std::invalid_argument("Missing measurement input collection");
   }
@@ -48,7 +45,6 @@ ActsExamples::SpacePointMaker::SpacePointMaker(Config cfg,
     throw std::invalid_argument("Missing space point maker geometry selection");
   }
 
-  m_inputSourceLinks.initialize(m_cfg.inputSourceLinks);
   m_inputMeasurements.initialize(m_cfg.inputMeasurements);
   m_outputSpacePoints.initialize(m_cfg.outputSpacePoints);
 
@@ -119,7 +115,6 @@ ActsExamples::SpacePointMaker::SpacePointMaker(Config cfg,
 
 ActsExamples::ProcessCode ActsExamples::SpacePointMaker::execute(
     const AlgorithmContext& ctx) const {
-  const auto& sourceLinks = m_inputSourceLinks(ctx);
   const auto& measurements = m_inputMeasurements(ctx);
 
   // TODO Support strip measurements
@@ -136,7 +131,8 @@ ActsExamples::ProcessCode ActsExamples::SpacePointMaker::execute(
   SimSpacePointContainer spacePoints;
   for (Acts::GeometryIdentifier geoId : m_cfg.geometrySelection) {
     // select volume/layer depending on what is set in the geometry id
-    auto range = selectLowestNonZeroGeometryObject(sourceLinks, geoId);
+    auto range =
+        selectLowestNonZeroGeometryObject(measurements.orderedIndices(), geoId);
     // groupByModule only works with geometry containers, not with an
     // arbitrary range. do the equivalent grouping manually
     auto groupedByModule = makeGroupBy(range, detail::GeometryIdGetter());
diff --git a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp
index ed9918c4f84..2a36701dc3b 100644
--- a/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp
+++ b/Examples/Algorithms/TrackFinding/src/TrackFindingAlgorithm.cpp
@@ -32,6 +32,7 @@
 #include "Acts/Utilities/Logger.hpp"
 #include "Acts/Utilities/TrackHelpers.hpp"
 #include "ActsExamples/EventData/IndexSourceLink.hpp"
+#include "ActsExamples/EventData/Measurement.hpp"
 #include "ActsExamples/EventData/MeasurementCalibration.hpp"
 #include "ActsExamples/EventData/SimSeed.hpp"
 #include "ActsExamples/EventData/Track.hpp"
@@ -262,9 +263,6 @@ TrackFindingAlgorithm::TrackFindingAlgorithm(Config config,
   if (m_cfg.inputMeasurements.empty()) {
     throw std::invalid_argument("Missing measurements input collection");
   }
-  if (m_cfg.inputSourceLinks.empty()) {
-    throw std::invalid_argument("Missing source links input collection");
-  }
   if (m_cfg.inputInitialTrackParameters.empty()) {
     throw std::invalid_argument(
         "Missing initial track parameters input collection");
@@ -293,7 +291,6 @@ TrackFindingAlgorithm::TrackFindingAlgorithm(Config config,
   }
 
   m_inputMeasurements.initialize(m_cfg.inputMeasurements);
-  m_inputSourceLinks.initialize(m_cfg.inputSourceLinks);
   m_inputInitialTrackParameters.initialize(m_cfg.inputInitialTrackParameters);
   m_inputSeeds.maybeInitialize(m_cfg.inputSeeds);
   m_outputTracks.initialize(m_cfg.outputTracks);
@@ -302,7 +299,6 @@ TrackFindingAlgorithm::TrackFindingAlgorithm(Config config,
 ProcessCode TrackFindingAlgorithm::execute(const AlgorithmContext& ctx) const {
   // Read input data
   const auto& measurements = m_inputMeasurements(ctx);
-  const auto& sourceLinks = m_inputSourceLinks(ctx);
   const auto& initialParameters = m_inputInitialTrackParameters(ctx);
   const SimSeedContainer* seeds = nullptr;
 
@@ -339,7 +335,7 @@ ProcessCode TrackFindingAlgorithm::execute(const AlgorithmContext& ctx) const {
   extensions.branchStopper.connect<&BranchStopper::operator()>(&branchStopper);
 
   IndexSourceLinkAccessor slAccessor;
-  slAccessor.container = &sourceLinks;
+  slAccessor.container = &measurements.orderedIndices();
   Acts::SourceLinkAccessorDelegate<IndexSourceLinkAccessor::Iterator>
       slAccessorDelegate;
   slAccessorDelegate.connect<&IndexSourceLinkAccessor::range>(&slAccessor);
@@ -648,7 +644,7 @@ ProcessCode TrackFindingAlgorithm::execute(const AlgorithmContext& ctx) const {
 
   // Compute shared hits from all the reconstructed tracks
   if (m_cfg.computeSharedHits) {
-    computeSharedHits(sourceLinks, tracks);
+    computeSharedHits(tracks, measurements);
   }
 
   ACTS_DEBUG("Finalized track finding with " << tracks.size()
@@ -698,4 +694,54 @@ ProcessCode TrackFindingAlgorithm::finalize() {
   return ProcessCode::SUCCESS;
 }
 
+// TODO this is somewhat duplicated in AmbiguityResolutionAlgorithm.cpp
+// TODO we should make a common implementation in the core at some point
+void TrackFindingAlgorithm::computeSharedHits(
+    TrackContainer& tracks, const MeasurementContainer& measurements) const {
+  // Compute shared hits from all the reconstructed tracks
+  // Compute nSharedhits and Update ckf results
+  // hit index -> list of multi traj indexes [traj, meas]
+
+  std::vector<std::size_t> firstTrackOnTheHit(
+      measurements.size(), std::numeric_limits<std::size_t>::max());
+  std::vector<std::size_t> firstStateOnTheHit(
+      measurements.size(), std::numeric_limits<std::size_t>::max());
+
+  for (auto track : tracks) {
+    for (auto state : track.trackStatesReversed()) {
+      if (!state.typeFlags().test(Acts::TrackStateFlag::MeasurementFlag)) {
+        continue;
+      }
+
+      std::size_t hitIndex = state.getUncalibratedSourceLink()
+                                 .template get<IndexSourceLink>()
+                                 .index();
+
+      // Check if hit not already used
+      if (firstTrackOnTheHit.at(hitIndex) ==
+          std::numeric_limits<std::size_t>::max()) {
+        firstTrackOnTheHit.at(hitIndex) = track.index();
+        firstStateOnTheHit.at(hitIndex) = state.index();
+        continue;
+      }
+
+      // if already used, control if first track state has been marked
+      // as shared
+      int indexFirstTrack = firstTrackOnTheHit.at(hitIndex);
+      int indexFirstState = firstStateOnTheHit.at(hitIndex);
+
+      auto firstState = tracks.getTrack(indexFirstTrack)
+                            .container()
+                            .trackStateContainer()
+                            .getTrackState(indexFirstState);
+      if (!firstState.typeFlags().test(Acts::TrackStateFlag::SharedHitFlag)) {
+        firstState.typeFlags().set(Acts::TrackStateFlag::SharedHitFlag);
+      }
+
+      // Decorate this track state
+      state.typeFlags().set(Acts::TrackStateFlag::SharedHitFlag);
+    }
+  }
+}
+
 }  // namespace ActsExamples
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/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/TrackFindingFromPrototrackAlgorithm.hpp b/Examples/Algorithms/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/TrackFindingFromPrototrackAlgorithm.hpp
index 91b3ce43cdb..abfc29b462c 100644
--- a/Examples/Algorithms/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/TrackFindingFromPrototrackAlgorithm.hpp
+++ b/Examples/Algorithms/TrackFindingExaTrkX/include/ActsExamples/TrackFindingExaTrkX/TrackFindingFromPrototrackAlgorithm.hpp
@@ -16,7 +16,6 @@
 #include "Acts/TrackFitting/GainMatrixSmoother.hpp"
 #include "Acts/TrackFitting/GainMatrixUpdater.hpp"
 #include "Acts/Utilities/Zip.hpp"
-#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 #include "ActsExamples/EventData/ProtoTrack.hpp"
 #include "ActsExamples/EventData/Trajectories.hpp"
@@ -44,9 +43,6 @@ class TrackFindingFromPrototrackAlgorithm final : public IAlgorithm {
     /// Input measurements
     std::string inputMeasurements;
 
-    /// Input source links
-    std::string inputSourceLinks;
-
     /// Input track parameters
     std::string inputInitialTrackParameters;
 
@@ -98,8 +94,6 @@ class TrackFindingFromPrototrackAlgorithm final : public IAlgorithm {
                                                          "InputProtoTracks"};
   ReadDataHandle<MeasurementContainer> m_inputMeasurements{this,
                                                            "InputMeasurements"};
-  ReadDataHandle<IndexSourceLinkContainer> m_inputSourceLinks{
-      this, "InputSourceLinks"};
   ReadDataHandle<TrackParametersContainer> m_inputInitialTrackParameters{
       this, "InputInitialTrackParameters"};
 
diff --git a/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingFromPrototrackAlgorithm.cpp b/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingFromPrototrackAlgorithm.cpp
index 4eef98a6f83..bf28b1743ba 100644
--- a/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingFromPrototrackAlgorithm.cpp
+++ b/Examples/Algorithms/TrackFindingExaTrkX/src/TrackFindingFromPrototrackAlgorithm.cpp
@@ -48,6 +48,7 @@ struct ProtoTrackSourceLinkAccessor
     return {Iterator{begin}, Iterator{end}};
   }
 };
+
 }  // namespace
 
 namespace ActsExamples {
@@ -58,14 +59,12 @@ TrackFindingFromPrototrackAlgorithm::TrackFindingFromPrototrackAlgorithm(
   m_inputInitialTrackParameters.initialize(m_cfg.inputInitialTrackParameters);
   m_inputMeasurements.initialize(m_cfg.inputMeasurements);
   m_inputProtoTracks.initialize(m_cfg.inputProtoTracks);
-  m_inputSourceLinks.initialize(m_cfg.inputSourceLinks);
   m_outputTracks.initialize(m_cfg.outputTracks);
 }
 
 ActsExamples::ProcessCode TrackFindingFromPrototrackAlgorithm::execute(
     const ActsExamples::AlgorithmContext& ctx) const {
   const auto& measurements = m_inputMeasurements(ctx);
-  const auto& sourceLinks = m_inputSourceLinks(ctx);
   const auto& protoTracks = m_inputProtoTracks(ctx);
   const auto& initialParameters = m_inputInitialTrackParameters(ctx);
 
@@ -98,7 +97,7 @@ ActsExamples::ProcessCode TrackFindingFromPrototrackAlgorithm::execute(
   // The source link accessor
   ProtoTrackSourceLinkAccessor sourceLinkAccessor;
   sourceLinkAccessor.loggerPtr = logger().clone("SourceLinkAccessor");
-  sourceLinkAccessor.container = &sourceLinks;
+  sourceLinkAccessor.container = &measurements.orderedIndices();
 
   Acts::SourceLinkAccessorDelegate<IndexSourceLinkAccessor::Iterator>
       slAccessorDelegate;
@@ -133,7 +132,8 @@ ActsExamples::ProcessCode TrackFindingFromPrototrackAlgorithm::execute(
 
     // Fill the source links via their indices from the container
     for (const auto hitIndex : protoTracks.at(i)) {
-      if (auto it = sourceLinks.nth(hitIndex); it != sourceLinks.end()) {
+      if (auto it = measurements.orderedIndices().nth(hitIndex);
+          it != measurements.orderedIndices().end()) {
         sourceLinkAccessor.protoTrackSourceLinks.insert(*it);
       } else {
         ACTS_FATAL("Proto track " << i << " contains invalid hit index"
@@ -177,7 +177,7 @@ ActsExamples::ProcessCode TrackFindingFromPrototrackAlgorithm::execute(
   // once this is done.
   // Compute shared hits from all the reconstructed tracks if
   // (m_cfg.computeSharedHits) {
-  //   computeSharedHits(sourceLinks, tracks);
+  //   computeSharedHits(measurements, tracks);
   // }
 
   ACTS_INFO("Event " << ctx.eventNumber << ": " << nFailed << " / " << nSeed
diff --git a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFittingAlgorithm.hpp b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFittingAlgorithm.hpp
index 3c3a62c4c4c..934203aeeaa 100644
--- a/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFittingAlgorithm.hpp
+++ b/Examples/Algorithms/TrackFitting/include/ActsExamples/TrackFitting/TrackFittingAlgorithm.hpp
@@ -72,8 +72,6 @@ class TrackFittingAlgorithm final : public IAlgorithm {
 
   ReadDataHandle<MeasurementContainer> m_inputMeasurements{this,
                                                            "InputMeasurements"};
-  ReadDataHandle<IndexSourceLinkContainer> m_inputSourceLinks{
-      this, "InputSourceLinks"};
   ReadDataHandle<ProtoTrackContainer> m_inputProtoTracks{this,
                                                          "InputProtoTracks"};
   ReadDataHandle<TrackParametersContainer> m_inputInitialTrackParameters{
diff --git a/Examples/Algorithms/TrackFitting/src/TrackFittingAlgorithm.cpp b/Examples/Algorithms/TrackFitting/src/TrackFittingAlgorithm.cpp
index 49ed2c285cd..50360e3d3b4 100644
--- a/Examples/Algorithms/TrackFitting/src/TrackFittingAlgorithm.cpp
+++ b/Examples/Algorithms/TrackFitting/src/TrackFittingAlgorithm.cpp
@@ -18,6 +18,7 @@
 #include "Acts/Surfaces/PerigeeSurface.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/Result.hpp"
+#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 #include "ActsExamples/EventData/MeasurementCalibration.hpp"
 #include "ActsExamples/EventData/ProtoTrack.hpp"
@@ -126,10 +127,11 @@ ActsExamples::ProcessCode ActsExamples::TrackFittingAlgorithm::execute(
     trackSourceLinks.reserve(protoTrack.size());
 
     // Fill the source links via their indices from the container
-    for (auto hitIndex : protoTrack) {
+    for (auto measIndex : protoTrack) {
       ConstVariableBoundMeasurementProxy measurement =
-          measurements.getMeasurement(hitIndex);
-      trackSourceLinks.push_back(measurement.sourceLink());
+          measurements.getMeasurement(measIndex);
+      IndexSourceLink sourceLink(measurement.geometryId(), measIndex);
+      trackSourceLinks.push_back(Acts::SourceLink(sourceLink));
     }
 
     ACTS_DEBUG("Invoke direct fitter for track " << itrack);
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/Algorithms/Utilities/include/ActsExamples/Utilities/MeasurementMapSelector.hpp b/Examples/Algorithms/Utilities/include/ActsExamples/Utilities/MeasurementMapSelector.hpp
index 3cf5eace6d9..ef495b3a9a9 100644
--- a/Examples/Algorithms/Utilities/include/ActsExamples/Utilities/MeasurementMapSelector.hpp
+++ b/Examples/Algorithms/Utilities/include/ActsExamples/Utilities/MeasurementMapSelector.hpp
@@ -17,6 +17,7 @@
 #include "ActsFatras/EventData/Barcode.hpp"
 
 #include <string>
+#include <utility>
 #include <vector>
 
 namespace ActsExamples {
@@ -30,12 +31,12 @@ class MeasurementMapSelector final : public IAlgorithm {
 
  public:
   struct Config {
+    /// Input measurements
+    std::string inputMeasurements;
+
     /// Input spacepoints collection.
     std::string inputMeasurementParticleMap;
 
-    /// Input source links
-    std::string inputSourceLinks;
-
     /// Output protoTracks collection.
     std::string outputMeasurementParticleMap;
 
@@ -48,27 +49,26 @@ class MeasurementMapSelector final : public IAlgorithm {
   /// @param cfg is the config struct to configure the algorithm
   /// @param level is the logging level
   MeasurementMapSelector(Config cfg, Acts::Logging::Level lvl)
-      : IAlgorithm("MeasurementMapSelector", lvl), m_cfg(cfg) {
-    m_inputSourceLinks.initialize(m_cfg.inputSourceLinks);
+      : IAlgorithm("MeasurementMapSelector", lvl), m_cfg(std::move(cfg)) {
+    m_inputMeasurements.initialize(m_cfg.inputMeasurements);
     m_inputMap.initialize(m_cfg.inputMeasurementParticleMap);
     m_outputMap.initialize(m_cfg.outputMeasurementParticleMap);
   }
 
-  virtual ~MeasurementMapSelector() = default;
-
   /// Filter the measurements
   ///
   /// @param ctx is the algorithm context that holds event-wise information
   /// @return a process code to steer the algorithm flow
   ActsExamples::ProcessCode execute(
       const ActsExamples::AlgorithmContext& ctx) const final {
-    const auto& inputSourceLinks = m_inputSourceLinks(ctx);
+    const auto& inputMeasurements = m_inputMeasurements(ctx);
     const auto& inputMap = m_inputMap(ctx);
 
     Map outputMap;
 
     for (const auto geoId : m_cfg.geometrySelection) {
-      auto range = selectLowestNonZeroGeometryObject(inputSourceLinks, geoId);
+      auto range = selectLowestNonZeroGeometryObject(
+          inputMeasurements.orderedIndices(), geoId);
       for (const auto& sl : range) {
         const auto [begin, end] = inputMap.equal_range(sl.index());
         outputMap.insert(begin, end);
@@ -86,8 +86,8 @@ class MeasurementMapSelector final : public IAlgorithm {
   // configuration
   Config m_cfg;
 
-  ReadDataHandle<IndexSourceLinkContainer> m_inputSourceLinks{
-      this, "InputSourceLinks"};
+  ReadDataHandle<MeasurementContainer> m_inputMeasurements{this,
+                                                           "InputMeasurements"};
   ReadDataHandle<Map> m_inputMap{this, "InputMeasurementParticleMap"};
   WriteDataHandle<Map> m_outputMap{this, "OutputMeasurementParticleMap"};
 };
diff --git a/Examples/Algorithms/Utilities/src/PrototracksToTracks.cpp b/Examples/Algorithms/Utilities/src/PrototracksToTracks.cpp
index 01e2fac6a38..12a2b70e6cf 100644
--- a/Examples/Algorithms/Utilities/src/PrototracksToTracks.cpp
+++ b/Examples/Algorithms/Utilities/src/PrototracksToTracks.cpp
@@ -8,6 +8,7 @@
 
 #include "ActsExamples/Utilities/ProtoTracksToTracks.hpp"
 
+#include "Acts/EventData/SourceLink.hpp"
 #include "Acts/Surfaces/PerigeeSurface.hpp"
 #include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/ProtoTrack.hpp"
@@ -27,16 +28,12 @@ PrototracksToTracks::PrototracksToTracks(Config cfg, Acts::Logging::Level lvl)
 }
 
 ProcessCode PrototracksToTracks::execute(const AlgorithmContext& ctx) const {
+  const auto& measurements = m_inputMeasurements(ctx);
+
   auto trackContainer = std::make_shared<Acts::VectorTrackContainer>();
   auto mtj = std::make_shared<Acts::VectorMultiTrajectory>();
   TrackContainer tracks(trackContainer, mtj);
 
-  boost::container::flat_map<Index, Acts::SourceLink> slMap;
-  for (const auto& m : m_inputMeasurements(ctx)) {
-    const auto idx = m.sourceLink().template get<IndexSourceLink>().index();
-    slMap.insert(std::pair<Index, Acts::SourceLink>{idx, m.sourceLink()});
-  }
-
   const auto& prototracks = m_inputProtoTracks(ctx);
   ACTS_DEBUG("Received " << prototracks.size() << " prototracks");
 
@@ -54,12 +51,15 @@ ProcessCode PrototracksToTracks::execute(const AlgorithmContext& ctx) const {
     maxSize = std::max(maxSize, protoTrack.size());
 
     auto track = tracks.makeTrack();
-    for (auto idx : protoTrack) {
+    for (auto measIndex : protoTrack) {
+      ConstVariableBoundMeasurementProxy measurement =
+          measurements.getMeasurement(measIndex);
+      IndexSourceLink sourceLink(measurement.geometryId(), measIndex);
+
       auto trackStateProxy =
           track.appendTrackState(Acts::TrackStatePropMask::None);
       trackStateProxy.typeFlags().set(Acts::TrackStateFlag::MeasurementFlag);
-      trackStateProxy.setUncalibratedSourceLink(
-          Acts::SourceLink{slMap.at(idx)});
+      trackStateProxy.setUncalibratedSourceLink(Acts::SourceLink(sourceLink));
     }
 
     track.nMeasurements() = static_cast<std::uint32_t>(protoTrack.size());
diff --git a/Examples/Algorithms/Utilities/src/TrajectoriesToPrototracks.cpp b/Examples/Algorithms/Utilities/src/TrajectoriesToPrototracks.cpp
index cd961ab6141..c92d985d492 100644
--- a/Examples/Algorithms/Utilities/src/TrajectoriesToPrototracks.cpp
+++ b/Examples/Algorithms/Utilities/src/TrajectoriesToPrototracks.cpp
@@ -9,6 +9,7 @@
 #include "ActsExamples/Utilities/TrajectoriesToPrototracks.hpp"
 
 #include "Acts/EventData/MultiTrajectory.hpp"
+#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/ProtoTrack.hpp"
 #include "ActsExamples/EventData/Trajectories.hpp"
 
@@ -16,7 +17,6 @@
 #include <vector>
 
 namespace ActsExamples {
-class IndexSourceLink;
 struct AlgorithmContext;
 
 TrajectoriesToPrototracks::TrajectoriesToPrototracks(Config cfg,
diff --git a/Examples/Framework/include/ActsExamples/EventData/IndexSourceLink.hpp b/Examples/Framework/include/ActsExamples/EventData/IndexSourceLink.hpp
index cc2e6249a2e..ea60a030b96 100644
--- a/Examples/Framework/include/ActsExamples/EventData/IndexSourceLink.hpp
+++ b/Examples/Framework/include/ActsExamples/EventData/IndexSourceLink.hpp
@@ -84,11 +84,6 @@ struct IndexSourceLinkSurfaceAccessor {
 
 }  // namespace Experimental
 
-/// Container of index source links.
-///
-/// Since the source links provide a `.geometryId()` accessor, they can be
-/// stored in an ordered geometry container.
-using IndexSourceLinkContainer = GeometryIdMultiset<IndexSourceLink>;
 /// Accessor for the above source link container
 ///
 /// It wraps up a few lookup methods to be used in the Combinatorial Kalman
@@ -105,4 +100,5 @@ struct IndexSourceLinkAccessor : GeometryIdMultisetAccessor<IndexSourceLink> {
     return {Iterator{begin}, Iterator{end}};
   }
 };
+
 }  // namespace ActsExamples
diff --git a/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp b/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp
index 6721740fbfa..a0bf4559509 100644
--- a/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp
+++ b/Examples/Framework/include/ActsExamples/EventData/Measurement.hpp
@@ -11,13 +11,15 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
 #include "Acts/EventData/MeasurementHelpers.hpp"
-#include "Acts/EventData/SourceLink.hpp"
 #include "Acts/EventData/SubspaceHelpers.hpp"
 #include "Acts/EventData/Types.hpp"
 #include "Acts/EventData/detail/CalculateResiduals.hpp"
 #include "Acts/EventData/detail/ParameterTraits.hpp"
 #include "Acts/EventData/detail/PrintParameters.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "Acts/Utilities/Iterator.hpp"
+#include "ActsExamples/EventData/GeometryContainers.hpp"
+#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/MeasurementConcept.hpp"
 
 #include <array>
@@ -69,6 +71,7 @@ class MeasurementContainer {
   using ConstFixedProxy = FixedMeasurementProxy<Acts::eBoundSize, Size, true>;
   using VariableProxy = VariableMeasurementProxy<Acts::eBoundSize, false>;
   using ConstVariableProxy = VariableMeasurementProxy<Acts::eBoundSize, true>;
+  using OrderedIndices = GeometryIdMultiset<IndexSourceLink>;
 
   MeasurementContainer();
 
@@ -82,8 +85,9 @@ class MeasurementContainer {
 
   /// @brief Add a measurement of a given size
   /// @param size The size of the measurement
+  /// @param geometryId The geometry identifier of the measurement surface
   /// @return The index of the added measurement
-  Index addMeasurement(std::uint8_t size);
+  Index addMeasurement(std::uint8_t size, Acts::GeometryIdentifier geometryId);
 
   /// @brief Get a variable-size measurement proxy
   /// @param index The index of the measurement
@@ -122,21 +126,29 @@ class MeasurementContainer {
 
   /// @brief Make a measurement of a given size
   /// @param size The size of the measurement
+  /// @param geometryId The geometry identifier of the measurement surface
   /// @return The variable-size measurement proxy
-  VariableProxy makeMeasurement(std::uint8_t size);
+  VariableProxy makeMeasurement(std::uint8_t size,
+                                Acts::GeometryIdentifier geometryId);
   /// @brief Make a fixed-size measurement
   /// @tparam Size The size of the measurement
+  /// @param geometryId The geometry identifier of the measurement surface
   /// @return The fixed-size measurement proxy
   template <std::size_t Size>
-  FixedProxy<Size> makeMeasurement() {
-    return getMeasurement<Size>(addMeasurement(Size));
+  FixedProxy<Size> makeMeasurement(Acts::GeometryIdentifier geometryId) {
+    return getMeasurement<Size>(addMeasurement(Size, geometryId));
   }
 
   template <typename... Args>
-  VariableProxy emplaceMeasurement(std::uint8_t size, Args&&... args);
+  VariableProxy emplaceMeasurement(std::uint8_t size,
+                                   Acts::GeometryIdentifier geometryId,
+                                   Args&&... args);
 
   template <std::size_t Size, typename... Args>
-  FixedProxy<Size> emplaceMeasurement(Args&&... args);
+  FixedProxy<Size> emplaceMeasurement(Acts::GeometryIdentifier geometryId,
+                                      Args&&... args);
+
+  const OrderedIndices& orderedIndices() const;
 
   using iterator =
       Acts::ContainerIndexIterator<MeasurementContainer, VariableProxy, false>;
@@ -161,10 +173,12 @@ class MeasurementContainer {
 
   std::vector<MeasurementEntry> m_entries;
 
-  std::vector<std::optional<Acts::SourceLink>> m_sourceLinks;
+  std::vector<Acts::GeometryIdentifier> m_geometryIds;
   std::vector<std::uint8_t> m_subspaceIndices;
   std::vector<double> m_parameters;
   std::vector<double> m_covariances;
+
+  OrderedIndices m_orderedIndices;
 };
 
 /// @brief Base class for measurement proxies
@@ -223,18 +237,10 @@ class MeasurementProxyBase {
     return self().subspaceHelper().indexOf(i);
   }
 
-  /// @brief Set the source link of the measurement
-  /// @param sourceLink The source link
-  void setSourceLink(const Acts::SourceLink& sourceLink)
-    requires(!ReadOnly)
-  {
-    container().m_sourceLinks.at(m_index) = sourceLink;
-  }
-
-  /// @brief Get the source link of the measurement
-  /// @return The source link
-  const Acts::SourceLink& sourceLink() const {
-    return container().m_sourceLinks.at(m_index).value();
+  /// @brief Get the geometry ID of the measurement
+  /// @return The geometry ID
+  Acts::GeometryIdentifier geometryId() const {
+    return container().m_geometryIds.at(m_index);
   }
 
   /// @brief Set the subspace indices of the measurement
@@ -262,23 +268,22 @@ class MeasurementProxyBase {
     return self().subspaceHelper().expandMatrix(self().covariance());
   }
 
-  /// @brief Construct the measurement from a sourcelink, subspace vector,
+  /// @brief Construct the measurement from a subspace vector,
   /// parameters, and covariance.
   ///
   template <typename Subspace, typename ParameterDerived,
             typename CovarianceDerived>
-  void fill(const Acts::SourceLink& source_link, Subspace&& subspace,
+  void fill(Subspace&& subspace,
             const Eigen::DenseBase<ParameterDerived>& parameters,
             const Eigen::DenseBase<CovarianceDerived>& covariance)
     requires(!ReadOnly)
   {
-    setSourceLink(source_link);
     self().setSubspaceIndices(std::forward<Subspace>(subspace));
     self().parameters() = parameters;
     self().covariance() = covariance;
   }
 
-  /// @brief Construct the measurement from a sourcelink, subspace vector,
+  /// @brief Construct the measurement from a subspace vector,
   /// parameters, and covariance.
   ///
   template <MeasurementConcept OtherDerived>
@@ -286,8 +291,7 @@ class MeasurementProxyBase {
     requires(!ReadOnly)
   {
     assert(size() == other.size() && "Size mismatch");
-    fill(other.sourceLink(), other.subspaceIndexVector(), other.parameters(),
-         other.covariance());
+    fill(other.subspaceIndexVector(), other.parameters(), other.covariance());
   }
 
   /// @brief Copy the data from another measurement
@@ -493,8 +497,8 @@ class VariableMeasurementProxy
 
 template <typename... Args>
 MeasurementContainer::VariableProxy MeasurementContainer::emplaceMeasurement(
-    std::uint8_t size, Args&&... args) {
-  VariableProxy meas = makeMeasurement(size);
+    std::uint8_t size, Acts::GeometryIdentifier geometryId, Args&&... args) {
+  VariableProxy meas = makeMeasurement(size, geometryId);
 
   meas.fill(std::forward<Args>(args)...);
 
@@ -503,8 +507,8 @@ MeasurementContainer::VariableProxy MeasurementContainer::emplaceMeasurement(
 
 template <std::size_t Size, typename... Args>
 MeasurementContainer::FixedProxy<Size> MeasurementContainer::emplaceMeasurement(
-    Args&&... args) {
-  FixedProxy<Size> meas = makeMeasurement<Size>();
+    Acts::GeometryIdentifier geometryId, Args&&... args) {
+  FixedProxy<Size> meas = makeMeasurement<Size>(geometryId);
 
   meas.fill(std::forward<Args>(args)...);
 
@@ -514,4 +518,5 @@ MeasurementContainer::FixedProxy<Size> MeasurementContainer::emplaceMeasurement(
 static_assert(
     std::random_access_iterator<MeasurementContainer::iterator> &&
     std::random_access_iterator<MeasurementContainer::const_iterator>);
+
 }  // namespace ActsExamples
diff --git a/Examples/Framework/include/ActsExamples/EventData/MeasurementConcept.hpp b/Examples/Framework/include/ActsExamples/EventData/MeasurementConcept.hpp
index f95ab117728..14627a78c53 100644
--- a/Examples/Framework/include/ActsExamples/EventData/MeasurementConcept.hpp
+++ b/Examples/Framework/include/ActsExamples/EventData/MeasurementConcept.hpp
@@ -9,6 +9,7 @@
 #pragma once
 
 #include "Acts/EventData/SourceLink.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
 
 #include <concepts>
 
@@ -17,7 +18,7 @@ namespace ActsExamples {
 template <typename T>
 concept MeasurementConcept = requires(const T& m) {
   { m.size() } -> std::integral;
-  { m.sourceLink() } -> Acts::Concepts::decayed_same_as<Acts::SourceLink>;
+  { m.geometryId() } -> std::same_as<Acts::GeometryIdentifier>;
   { m.subspaceIndexVector() };
   { m.parameters() };
   { m.covariance() };
diff --git a/Examples/Framework/src/EventData/Measurement.cpp b/Examples/Framework/src/EventData/Measurement.cpp
index f37661948db..ab2dad07489 100644
--- a/Examples/Framework/src/EventData/Measurement.cpp
+++ b/Examples/Framework/src/EventData/Measurement.cpp
@@ -8,6 +8,9 @@
 
 #include "ActsExamples/EventData/Measurement.hpp"
 
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "ActsExamples/EventData/IndexSourceLink.hpp"
+
 namespace ActsExamples {
 
 MeasurementContainer::MeasurementContainer() = default;
@@ -17,20 +20,26 @@ std::size_t MeasurementContainer::size() const {
 }
 
 void MeasurementContainer::reserve(std::size_t size) {
-  m_sourceLinks.reserve(size);
+  m_geometryIds.reserve(size);
   m_subspaceIndices.reserve(size * 2);
   m_parameters.reserve(size * 2);
   m_covariances.reserve(size * 2 * 2);
 }
 
-std::size_t MeasurementContainer::addMeasurement(std::uint8_t size) {
+std::size_t MeasurementContainer::addMeasurement(
+    std::uint8_t size, Acts::GeometryIdentifier geometryId) {
   m_entries.push_back({m_subspaceIndices.size(), m_parameters.size(),
                        m_covariances.size(), size});
-  m_sourceLinks.emplace_back();
+  m_geometryIds.emplace_back(geometryId);
   m_subspaceIndices.resize(m_subspaceIndices.size() + size);
   m_parameters.resize(m_parameters.size() + size);
   m_covariances.resize(m_covariances.size() + size * size);
-  return m_entries.size() - 1;
+
+  std::size_t index = m_entries.size() - 1;
+  IndexSourceLink sourceLink(geometryId, index);
+  m_orderedIndices.emplace_hint(m_orderedIndices.end(), sourceLink);
+
+  return index;
 }
 
 MeasurementContainer::VariableProxy MeasurementContainer::at(
@@ -54,8 +63,13 @@ MeasurementContainer::ConstVariableProxy MeasurementContainer::getMeasurement(
 }
 
 MeasurementContainer::VariableProxy MeasurementContainer::makeMeasurement(
-    std::uint8_t size) {
-  return getMeasurement(addMeasurement(size));
+    std::uint8_t size, Acts::GeometryIdentifier geometryId) {
+  return getMeasurement(addMeasurement(size, geometryId));
+}
+
+const MeasurementContainer::OrderedIndices&
+MeasurementContainer::orderedIndices() const {
+  return m_orderedIndices;
 }
 
 MeasurementContainer::iterator MeasurementContainer::begin() {
diff --git a/Examples/Framework/src/Framework/Sequencer.cpp b/Examples/Framework/src/Framework/Sequencer.cpp
index 4d734da7e85..5c1a76b5bd4 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/Io/Csv/include/ActsExamples/Io/Csv/CsvMeasurementReader.hpp b/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvMeasurementReader.hpp
index c26387e1260..d3c5359d16e 100644
--- a/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvMeasurementReader.hpp
+++ b/Examples/Io/Csv/include/ActsExamples/Io/Csv/CsvMeasurementReader.hpp
@@ -13,7 +13,6 @@
 #include "ActsExamples/EventData/Cluster.hpp"
 #include "ActsExamples/EventData/GeometryContainers.hpp"
 #include "ActsExamples/EventData/Index.hpp"
-#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 #include "ActsExamples/EventData/SimHit.hpp"
 #include "ActsExamples/EventData/SimParticle.hpp"
@@ -32,7 +31,6 @@ class Surface;
 }
 
 namespace ActsExamples {
-class IndexSourceLink;
 struct AlgorithmContext;
 
 /// Read in a measurement cluster collection in comma-separated-value format.
@@ -60,8 +58,6 @@ class CsvMeasurementReader final : public IReader {
     std::string outputMeasurements;
     /// Output measurement to sim hit collection.
     std::string outputMeasurementSimHitsMap;
-    /// Output source links collection.
-    std::string outputSourceLinks;
     /// Output cluster collection (optional).
     std::string outputClusters;
 
@@ -102,9 +98,6 @@ class CsvMeasurementReader final : public IReader {
   WriteDataHandle<IndexMultimap<Index>> m_outputMeasurementSimHitsMap{
       this, "OutputMeasurementSimHitsMap"};
 
-  WriteDataHandle<GeometryIdMultiset<IndexSourceLink>> m_outputSourceLinks{
-      this, "OutputSourceLinks"};
-
   WriteDataHandle<ClusterContainer> m_outputClusters{this, "OutputClusters"};
 
   WriteDataHandle<IndexMultimap<ActsFatras::Barcode>>
diff --git a/Examples/Io/Csv/src/CsvMeasurementReader.cpp b/Examples/Io/Csv/src/CsvMeasurementReader.cpp
index cbc030c6939..f426784ea30 100644
--- a/Examples/Io/Csv/src/CsvMeasurementReader.cpp
+++ b/Examples/Io/Csv/src/CsvMeasurementReader.cpp
@@ -45,7 +45,6 @@ ActsExamples::CsvMeasurementReader::CsvMeasurementReader(
 
   m_outputMeasurements.initialize(m_cfg.outputMeasurements);
   m_outputMeasurementSimHitsMap.initialize(m_cfg.outputMeasurementSimHitsMap);
-  m_outputSourceLinks.initialize(m_cfg.outputSourceLinks);
   m_outputClusters.maybeInitialize(m_cfg.outputClusters);
   m_outputMeasurementParticlesMap.maybeInitialize(
       m_cfg.outputMeasurementParticlesMap);
@@ -196,15 +195,11 @@ ActsExamples::ProcessCode ActsExamples::CsvMeasurementReader::read(
   MeasurementContainer tmpMeasurements;
   GeometryIdMultimap<ConstVariableBoundMeasurementProxy> orderedMeasurements;
   IndexMultimap<Index> measurementSimHitsMap;
-  IndexSourceLinkContainer sourceLinks;
-  // need list here for stable addresses
-  std::list<IndexSourceLink> sourceLinkStorage;
 
   tmpMeasurements.reserve(measurementData.size());
   orderedMeasurements.reserve(measurementData.size());
   // Safe long as we have single particle to sim hit association
   measurementSimHitsMap.reserve(measurementData.size());
-  sourceLinks.reserve(measurementData.size());
 
   auto measurementSimHitLinkData =
       readEverything<ActsExamples::MeasurementSimHitLink>(
@@ -252,10 +247,7 @@ ActsExamples::ProcessCode ActsExamples::CsvMeasurementReader::read(
 
     // The measurement container is unordered and the index under which
     // the measurement will be stored is known before adding it.
-    const Index index = orderedMeasurements.size();
-    IndexSourceLink& sourceLink = sourceLinkStorage.emplace_back(geoId, index);
-    auto measurement =
-        createMeasurement(tmpMeasurements, dParameters, sourceLink);
+    auto measurement = createMeasurement(tmpMeasurements, geoId, dParameters);
 
     // Due to the previous sorting of the raw hit data by geometry id, new
     // measurements should always end up at the end of the container. previous
@@ -267,13 +259,11 @@ ActsExamples::ProcessCode ActsExamples::CsvMeasurementReader::read(
       ACTS_FATAL("Something went horribly wrong with the hit sorting");
       return ProcessCode::ABORT;
     }
-
-    sourceLinks.insert(sourceLinks.end(), std::cref(sourceLink));
   }
 
   MeasurementContainer measurements;
   for (auto& [_, meas] : orderedMeasurements) {
-    measurements.emplaceMeasurement(meas.size(), meas);
+    measurements.emplaceMeasurement(meas.size(), meas.geometryId(), meas);
   }
 
   // Generate measurement-particles-map
@@ -294,7 +284,6 @@ ActsExamples::ProcessCode ActsExamples::CsvMeasurementReader::read(
   // Write the data to the EventStore
   m_outputMeasurements(ctx, std::move(measurements));
   m_outputMeasurementSimHitsMap(ctx, std::move(measurementSimHitsMap));
-  m_outputSourceLinks(ctx, std::move(sourceLinks));
 
   /////////////////////////
   // Cluster information //
diff --git a/Examples/Io/Csv/src/CsvMeasurementWriter.cpp b/Examples/Io/Csv/src/CsvMeasurementWriter.cpp
index f421c8c2501..851966e8fba 100644
--- a/Examples/Io/Csv/src/CsvMeasurementWriter.cpp
+++ b/Examples/Io/Csv/src/CsvMeasurementWriter.cpp
@@ -13,7 +13,6 @@
 #include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "ActsExamples/EventData/Cluster.hpp"
 #include "ActsExamples/EventData/Index.hpp"
-#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 #include "ActsExamples/Framework/AlgorithmContext.hpp"
 #include "ActsExamples/Io/Csv/CsvInputOutput.hpp"
@@ -101,8 +100,7 @@ ActsExamples::ProcessCode ActsExamples::CsvMeasurementWriter::writeT(
       writerMeasurementSimHitMap.append({measIdx, simHitIdx});
     }
 
-    Acts::GeometryIdentifier geoId =
-        measurement.sourceLink().template get<IndexSourceLink>().geometryId();
+    Acts::GeometryIdentifier geoId = measurement.geometryId();
     // MEASUREMENT information ------------------------------------
 
     // Encoded geometry identifier. same for all hits on the module
diff --git a/Examples/Io/Csv/src/CsvTrackWriter.cpp b/Examples/Io/Csv/src/CsvTrackWriter.cpp
index 26ac3ad2e6b..9f6d74ec730 100644
--- a/Examples/Io/Csv/src/CsvTrackWriter.cpp
+++ b/Examples/Io/Csv/src/CsvTrackWriter.cpp
@@ -33,10 +33,6 @@
 #include <unordered_set>
 #include <utility>
 
-namespace ActsExamples {
-class IndexSourceLink;
-}  // namespace ActsExamples
-
 using namespace ActsExamples;
 
 CsvTrackWriter::CsvTrackWriter(const CsvTrackWriter::Config& config,
diff --git a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementReader.hpp b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementReader.hpp
index a81341f0edc..34a95d3709a 100644
--- a/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementReader.hpp
+++ b/Examples/Io/EDM4hep/include/ActsExamples/Io/EDM4hep/EDM4hepMeasurementReader.hpp
@@ -41,8 +41,6 @@ class EDM4hepMeasurementReader final : public IReader {
     std::string outputMeasurements;
     /// Output measurement to sim hit collection.
     std::string outputMeasurementSimHitsMap;
-    /// Output source links collection.
-    std::string outputSourceLinks;
     /// Output cluster collection (optional).
     std::string outputClusters;
   };
@@ -79,9 +77,6 @@ class EDM4hepMeasurementReader final : public IReader {
   WriteDataHandle<IndexMultimap<Index>> m_outputMeasurementSimHitsMap{
       this, "OutputMeasurementSimHitsMap"};
 
-  WriteDataHandle<GeometryIdMultiset<IndexSourceLink>> m_outputSourceLinks{
-      this, "OutputSourceLinks"};
-
   WriteDataHandle<ClusterContainer> m_outputClusters{this, "OutputClusters"};
 };
 
diff --git a/Examples/Io/EDM4hep/src/EDM4hepMeasurementReader.cpp b/Examples/Io/EDM4hep/src/EDM4hepMeasurementReader.cpp
index 5a9d6840999..237ff19e3b6 100644
--- a/Examples/Io/EDM4hep/src/EDM4hepMeasurementReader.cpp
+++ b/Examples/Io/EDM4hep/src/EDM4hepMeasurementReader.cpp
@@ -37,7 +37,6 @@ EDM4hepMeasurementReader::EDM4hepMeasurementReader(
 
   m_outputMeasurements.initialize(m_cfg.outputMeasurements);
   m_outputMeasurementSimHitsMap.initialize(m_cfg.outputMeasurementSimHitsMap);
-  m_outputSourceLinks.initialize(m_cfg.outputSourceLinks);
   m_outputClusters.maybeInitialize(m_cfg.outputClusters);
 }
 
@@ -55,7 +54,6 @@ ProcessCode EDM4hepMeasurementReader::read(const AlgorithmContext& ctx) {
   ClusterContainer clusters;
   // TODO what about those?
   IndexMultimap<Index> measurementSimHitsMap;
-  IndexSourceLinkContainer sourceLinks;
 
   podio::Frame frame = reader().readEntry("events", ctx.eventNumber);
 
@@ -76,7 +74,6 @@ ProcessCode EDM4hepMeasurementReader::read(const AlgorithmContext& ctx) {
   // Write the data to the EventStore
   m_outputMeasurements(ctx, std::move(measurements));
   m_outputMeasurementSimHitsMap(ctx, std::move(measurementSimHitsMap));
-  m_outputSourceLinks(ctx, std::move(sourceLinks));
   if (!m_cfg.outputClusters.empty()) {
     m_outputClusters(ctx, std::move(clusters));
   }
diff --git a/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp b/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp
index ef6b794cc4c..09dcbc6331b 100644
--- a/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp
+++ b/Examples/Io/EDM4hep/src/EDM4hepUtil.cpp
@@ -152,9 +152,6 @@ VariableBoundMeasurementProxy EDM4hepUtil::readMeasurement(
   // no need for digitization as we only want to identify the sensor
   Acts::GeometryIdentifier geometryId = geometryMapper(from.getCellID());
 
-  IndexSourceLink sourceLink{
-      geometryId, static_cast<Index>(podioObjectIDToInteger(from.id()))};
-
   auto pos = from.getPosition();
   auto cov = from.getCovMatrix();
 
@@ -173,7 +170,7 @@ VariableBoundMeasurementProxy EDM4hepUtil::readMeasurement(
   dParameters.values.push_back(pos.z);
   dParameters.variances.push_back(cov[5]);
 
-  auto to = createMeasurement(container, dParameters, sourceLink);
+  auto to = createMeasurement(container, geometryId, dParameters);
 
   if (fromClusters != nullptr) {
     for (const auto objectId : from.getRawHits()) {
@@ -202,8 +199,7 @@ void EDM4hepUtil::writeMeasurement(
     edm4hep::MutableTrackerHitPlane to, const Cluster* fromCluster,
     edm4hep::TrackerHitCollection& toClusters,
     const MapGeometryIdTo& geometryMapper) {
-  Acts::GeometryIdentifier geoId =
-      from.sourceLink().template get<IndexSourceLink>().geometryId();
+  Acts::GeometryIdentifier geoId = from.geometryId();
 
   if (geometryMapper) {
     // no need for digitization as we only want to identify the sensor
diff --git a/Examples/Io/Root/src/RootAthenaDumpReader.cpp b/Examples/Io/Root/src/RootAthenaDumpReader.cpp
index 3ea44e78016..eea524fe544 100644
--- a/Examples/Io/Root/src/RootAthenaDumpReader.cpp
+++ b/Examples/Io/Root/src/RootAthenaDumpReader.cpp
@@ -299,7 +299,6 @@ RootAthenaDumpReader::readMeasurements(
   IndexMultimap<ActsFatras::Barcode> measPartMap;
 
   // We cannot use im for the index since we might skip measurements
-  std::size_t idx = 0;
   std::unordered_map<int, std::size_t> imIdxMap;
 
   for (int im = 0; im < nCL; im++) {
@@ -375,7 +374,7 @@ RootAthenaDumpReader::readMeasurements(
                                 << CLloc_direction3[im]);
     const auto& locCov = CLlocal_cov->at(im);
 
-    std::optional<IndexSourceLink> sl;
+    Acts::GeometryIdentifier geoId;
     std::vector<double> localParams;
     if (m_cfg.geometryIdMap && m_cfg.trackingGeometry) {
       const auto& geoIdMap = m_cfg.geometryIdMap->left;
@@ -384,8 +383,7 @@ RootAthenaDumpReader::readMeasurements(
         continue;
       }
 
-      auto geoId = m_cfg.geometryIdMap->left.at(CLmoduleID[im]);
-      sl = IndexSourceLink(geoId, idx);
+      geoId = m_cfg.geometryIdMap->left.at(CLmoduleID[im]);
 
       auto surface = m_cfg.trackingGeometry->findSurface(geoId);
       if (surface == nullptr) {
@@ -427,7 +425,7 @@ RootAthenaDumpReader::readMeasurements(
       // bounds?
       localParams = std::vector<double>(loc->begin(), loc->end());
     } else {
-      sl = IndexSourceLink(Acts::GeometryIdentifier(CLmoduleID[im]), idx);
+      geoId = Acts::GeometryIdentifier(CLmoduleID[im]);
       localParams = {CLloc_direction1[im], CLloc_direction2[im]};
     }
 
@@ -444,7 +442,8 @@ RootAthenaDumpReader::readMeasurements(
       digiPars.values = {localParams[0]};
     }
 
-    createMeasurement(measurements, digiPars, *sl);
+    std::size_t measIndex = measurements.size();
+    createMeasurement(measurements, geoId, digiPars);
 
     // Create measurement particles map and particles container
     for (const auto& [subevt, barcode] :
@@ -459,12 +458,10 @@ RootAthenaDumpReader::readMeasurements(
         particles.emplace(dummyBarcode, Acts::PdgParticle::eInvalid);
       }
       measPartMap.insert(
-          std::pair<Index, ActsFatras::Barcode>{idx, dummyBarcode});
+          std::pair<Index, ActsFatras::Barcode>{measIndex, dummyBarcode});
     }
 
-    // Finally increment the measurement index
-    imIdxMap.emplace(im, idx);
-    ++idx;
+    imIdxMap.emplace(im, measIndex);
   }
 
   if (measurements.size() < static_cast<std::size_t>(nCL)) {
diff --git a/Examples/Io/Root/src/RootMeasurementWriter.cpp b/Examples/Io/Root/src/RootMeasurementWriter.cpp
index a716cb721a3..b7ab332fdea 100644
--- a/Examples/Io/Root/src/RootMeasurementWriter.cpp
+++ b/Examples/Io/Root/src/RootMeasurementWriter.cpp
@@ -12,7 +12,6 @@
 #include "Acts/Utilities/Enumerate.hpp"
 #include "ActsExamples/EventData/AverageSimHits.hpp"
 #include "ActsExamples/EventData/Index.hpp"
-#include "ActsExamples/EventData/IndexSourceLink.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 #include "ActsExamples/Framework/AlgorithmContext.hpp"
 #include "ActsExamples/Utilities/Range.hpp"
@@ -270,8 +269,7 @@ ProcessCode RootMeasurementWriter::writeT(
     const ConstVariableBoundMeasurementProxy meas =
         measurements.getMeasurement(hitIdx);
 
-    Acts::GeometryIdentifier geoId =
-        meas.sourceLink().template get<IndexSourceLink>().geometryId();
+    Acts::GeometryIdentifier geoId = meas.geometryId();
     // find the corresponding surface
     auto surfaceItr = m_cfg.surfaceByIdentifier.find(geoId);
     if (surfaceItr == m_cfg.surfaceByIdentifier.end()) {
diff --git a/Examples/Io/Root/src/RootTrackStatesWriter.cpp b/Examples/Io/Root/src/RootTrackStatesWriter.cpp
index ad8467f239c..160eada776b 100644
--- a/Examples/Io/Root/src/RootTrackStatesWriter.cpp
+++ b/Examples/Io/Root/src/RootTrackStatesWriter.cpp
@@ -44,8 +44,6 @@
 
 namespace ActsExamples {
 
-class IndexSourceLink;
-
 using Acts::VectorHelpers::eta;
 using Acts::VectorHelpers::perp;
 using Acts::VectorHelpers::phi;
diff --git a/Examples/Python/python/acts/examples/reconstruction.py b/Examples/Python/python/acts/examples/reconstruction.py
index 7ac0989744f..371345a17b5 100644
--- a/Examples/Python/python/acts/examples/reconstruction.py
+++ b/Examples/Python/python/acts/examples/reconstruction.py
@@ -403,7 +403,6 @@ def addSeeding(
             logger.info("Using Hough Transform seeding")
             houghTransformConfig.inputSpacePoints = [spacePoints]
             houghTransformConfig.inputMeasurements = "measurements"
-            houghTransformConfig.inputSourceLinks = "sourcelinks"
             houghTransformConfig.outputProtoTracks = "prototracks"
             houghTransformConfig.outputSeeds = "seeds"
             houghTransformConfig.trackingGeometry = trackingGeometry
@@ -635,7 +634,6 @@ def addSpacePointsMaking(
     logLevel = acts.examples.defaultLogging(sequence, logLevel)()
     spAlg = acts.examples.SpacePointMaker(
         level=logLevel,
-        inputSourceLinks="sourcelinks",
         inputMeasurements="measurements",
         outputSpacePoints="spacepoints",
         trackingGeometry=trackingGeometry,
@@ -1168,7 +1166,6 @@ def addGbtsSeeding(
         geometrySelection=acts.examples.readJsonGeometryList(
             str(geoSelectionConfigFile)
         ),
-        inputSourceLinks="sourcelinks",
         trackingGeometry=trackingGeometry,
         fill_module_csv=False,
         inputClusters="clusters",
@@ -1378,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,
@@ -1528,7 +1530,6 @@ def addCKFTracks(
             ]
         ),
         inputMeasurements="measurements",
-        inputSourceLinks="sourcelinks",
         inputInitialTrackParameters="estimatedparameters",
         inputSeeds=(
             "estimatedseeds"
@@ -1813,7 +1814,6 @@ def addExaTrkX(
     s.addAlgorithm(
         acts.examples.SpacePointMaker(
             level=customLogLevel(),
-            inputSourceLinks="sourcelinks",
             inputMeasurements="measurements",
             outputSpacePoints="spacepoints",
             trackingGeometry=trackingGeometry,
diff --git a/Examples/Python/python/acts/examples/simulation.py b/Examples/Python/python/acts/examples/simulation.py
index 125ce422c16..b3b1046dc70 100644
--- a/Examples/Python/python/acts/examples/simulation.py
+++ b/Examples/Python/python/acts/examples/simulation.py
@@ -818,7 +818,6 @@ def addDigitization(
         surfaceByIdentifier=trackingGeometry.geoIdSurfaceMap(),
         randomNumbers=rnd,
         inputSimHits="simhits",
-        outputSourceLinks="sourcelinks",
         outputMeasurements="measurements",
         outputMeasurementParticlesMap="measurement_particles_map",
         outputMeasurementSimHitsMap="measurement_simhits_map",
diff --git a/Examples/Python/src/Base.cpp b/Examples/Python/src/Base.cpp
index 4d5eec557bd..cc8a5f67dc1 100644
--- a/Examples/Python/src/Base.cpp
+++ b/Examples/Python/src/Base.cpp
@@ -364,12 +364,16 @@ void addBinning(Context& ctx) {
                           .value("binY", Acts::BinningValue::binY)
                           .value("binZ", Acts::BinningValue::binZ)
                           .value("binR", Acts::BinningValue::binR)
-                          .value("binPhi", Acts::BinningValue::binPhi);
+                          .value("binPhi", Acts::BinningValue::binPhi)
+                          .value("binRPhi", Acts::BinningValue::binRPhi)
+                          .value("binH", Acts::BinningValue::binH)
+                          .value("binEta", Acts::BinningValue::binEta)
+                          .value("binMag", Acts::BinningValue::binMag);
 
   auto boundaryType = py::enum_<Acts::AxisBoundaryType>(m, "AxisBoundaryType")
-                          .value("bound", Acts::AxisBoundaryType::Bound)
-                          .value("closed", Acts::AxisBoundaryType::Closed)
-                          .value("open", Acts::AxisBoundaryType::Open);
+                          .value("Bound", Acts::AxisBoundaryType::Bound)
+                          .value("Closed", Acts::AxisBoundaryType::Closed)
+                          .value("Open", Acts::AxisBoundaryType::Open);
 
   auto axisType = py::enum_<Acts::AxisType>(m, "AxisType")
                       .value("equidistant", Acts::AxisType::Equidistant)
diff --git a/Examples/Python/src/Digitization.cpp b/Examples/Python/src/Digitization.cpp
index a68d4871724..c392d6dcc8c 100644
--- a/Examples/Python/src/Digitization.cpp
+++ b/Examples/Python/src/Digitization.cpp
@@ -62,7 +62,6 @@ void addDigitization(Context& ctx) {
 
     ACTS_PYTHON_STRUCT_BEGIN(c, Config);
     ACTS_PYTHON_MEMBER(inputSimHits);
-    ACTS_PYTHON_MEMBER(outputSourceLinks);
     ACTS_PYTHON_MEMBER(outputMeasurements);
     ACTS_PYTHON_MEMBER(outputClusters);
     ACTS_PYTHON_MEMBER(outputMeasurementParticlesMap);
diff --git a/Examples/Python/src/EDM4hepComponent.cpp b/Examples/Python/src/EDM4hepComponent.cpp
index 5017596b7db..75746587be8 100644
--- a/Examples/Python/src/EDM4hepComponent.cpp
+++ b/Examples/Python/src/EDM4hepComponent.cpp
@@ -42,7 +42,7 @@ PYBIND11_MODULE(ActsPythonBindingsEDM4hep, m) {
   ACTS_PYTHON_DECLARE_READER(ActsExamples::EDM4hepMeasurementReader, m,
                              "EDM4hepMeasurementReader", inputPath,
                              outputMeasurements, outputMeasurementSimHitsMap,
-                             outputSourceLinks, outputClusters);
+                             outputClusters);
 
   ACTS_PYTHON_DECLARE_WRITER(ActsExamples::EDM4hepMeasurementWriter, m,
                              "EDM4hepMeasurementWriter", inputMeasurements,
diff --git a/Examples/Python/src/ExaTrkXTrackFinding.cpp b/Examples/Python/src/ExaTrkXTrackFinding.cpp
index db9190e7f66..3778e8576f5 100644
--- a/Examples/Python/src/ExaTrkXTrackFinding.cpp
+++ b/Examples/Python/src/ExaTrkXTrackFinding.cpp
@@ -251,9 +251,8 @@ void addExaTrkXTrackFinding(Context &ctx) {
   ACTS_PYTHON_DECLARE_ALGORITHM(
       ActsExamples::TrackFindingFromPrototrackAlgorithm, mex,
       "TrackFindingFromPrototrackAlgorithm", inputProtoTracks,
-      inputMeasurements, inputSourceLinks, inputInitialTrackParameters,
-      outputTracks, measurementSelectorCfg, trackingGeometry, magneticField,
-      findTracks, tag);
+      inputMeasurements, inputInitialTrackParameters, outputTracks,
+      measurementSelectorCfg, trackingGeometry, magneticField, findTracks, tag);
 }
 
 }  // namespace Acts::Python
diff --git a/Examples/Python/src/GeoModel.cpp b/Examples/Python/src/GeoModel.cpp
index eb9b4406808..58609a573d5 100644
--- a/Examples/Python/src/GeoModel.cpp
+++ b/Examples/Python/src/GeoModel.cpp
@@ -192,8 +192,9 @@ void addGeoModel(Context& ctx) {
         .def_readwrite(
             "topBoundsOverride",
             &Acts::GeoModelBlueprintCreater::Options::topBoundsOverride)
-        .def_readwrite("table",
-                       &Acts::GeoModelBlueprintCreater::Options::table);
+        .def_readwrite("table", &Acts::GeoModelBlueprintCreater::Options::table)
+        .def_readwrite("dotGraph",
+                       &Acts::GeoModelBlueprintCreater::Options::dotGraph);
   }
 
   gm.def(
diff --git a/Examples/Python/src/Geometry.cpp b/Examples/Python/src/Geometry.cpp
index 8ec8ed97ef1..d676dbe6243 100644
--- a/Examples/Python/src/Geometry.cpp
+++ b/Examples/Python/src/Geometry.cpp
@@ -25,6 +25,7 @@
 #include "Acts/Detector/interface/IInternalStructureBuilder.hpp"
 #include "Acts/Detector/interface/IRootVolumeFinderBuilder.hpp"
 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
+#include "Acts/Geometry/CylinderVolumeStack.hpp"
 #include "Acts/Geometry/Extent.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Geometry/GeometryHierarchyMap.hpp"
@@ -36,6 +37,7 @@
 #include "Acts/Plugins/Python/Utilities.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Surfaces/SurfaceArray.hpp"
+#include "Acts/Utilities/BinningType.hpp"
 #include "Acts/Utilities/Helpers.hpp"
 #include "Acts/Utilities/RangeXD.hpp"
 #include "Acts/Visualization/ViewConfig.hpp"
@@ -46,6 +48,7 @@
 #include <unordered_map>
 #include <vector>
 
+#include <boost/algorithm/string/join.hpp>
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
 
@@ -106,7 +109,12 @@ void addGeometry(Context& ctx) {
         .def("approach", &Acts::GeometryIdentifier::approach)
         .def("sensitive", &Acts::GeometryIdentifier::sensitive)
         .def("extra", &Acts::GeometryIdentifier::extra)
-        .def("value", &Acts::GeometryIdentifier::value);
+        .def("value", &Acts::GeometryIdentifier::value)
+        .def("__str__", [](const Acts::GeometryIdentifier& self) {
+          std::stringstream ss;
+          ss << self;
+          return ss.str();
+        });
   }
 
   {
@@ -178,15 +186,31 @@ void addGeometry(Context& ctx) {
 
   {
     py::class_<Acts::VolumeBounds, std::shared_ptr<Acts::VolumeBounds>>(
-        m, "VolumeBounds");
-
-    py::class_<Acts::CylinderVolumeBounds,
-               std::shared_ptr<Acts::CylinderVolumeBounds>, Acts::VolumeBounds>(
-        m, "CylinderVolumeBounds")
-        .def(py::init<ActsScalar, ActsScalar, ActsScalar, ActsScalar,
-                      ActsScalar, ActsScalar, ActsScalar>(),
-             "rmin"_a, "rmax"_a, "halfz"_a, "halfphi"_a = M_PI, "avgphi"_a = 0.,
-             "bevelMinZ"_a = 0., "bevelMaxZ"_a = 0.);
+        m, "VolumeBounds")
+        .def("type", &Acts::VolumeBounds::type)
+        .def("__str__", [](const Acts::VolumeBounds& self) {
+          std::stringstream ss;
+          ss << self;
+          return ss.str();
+        });
+
+    auto cvb =
+        py::class_<Acts::CylinderVolumeBounds,
+                   std::shared_ptr<Acts::CylinderVolumeBounds>,
+                   Acts::VolumeBounds>(m, "CylinderVolumeBounds")
+            .def(py::init<ActsScalar, ActsScalar, ActsScalar, ActsScalar,
+                          ActsScalar, ActsScalar, ActsScalar>(),
+                 "rmin"_a, "rmax"_a, "halfz"_a, "halfphi"_a = M_PI,
+                 "avgphi"_a = 0., "bevelMinZ"_a = 0., "bevelMaxZ"_a = 0.);
+
+    py::enum_<CylinderVolumeBounds::Face>(cvb, "Face")
+        .value("PositiveDisc", CylinderVolumeBounds::Face::PositiveDisc)
+        .value("NegativeDisc", CylinderVolumeBounds::Face::NegativeDisc)
+        .value("OuterCylinder", CylinderVolumeBounds::Face::OuterCylinder)
+        .value("InnerCylinder", CylinderVolumeBounds::Face::InnerCylinder)
+        .value("NegativePhiPlane", CylinderVolumeBounds::Face::NegativePhiPlane)
+        .value("PositivePhiPlane",
+               CylinderVolumeBounds::Face::PositivePhiPlane);
   }
 
   {
@@ -209,22 +233,72 @@ void addGeometry(Context& ctx) {
         }));
   }
 
+  py::class_<ExtentEnvelope>(m, "ExtentEnvelope")
+      .def(py::init<>())
+      .def(py::init<const Envelope&>())
+      .def(py::init([](Envelope x, Envelope y, Envelope z, Envelope r,
+                       Envelope phi, Envelope rPhi, Envelope h, Envelope eta,
+                       Envelope mag) {
+             return ExtentEnvelope({.x = x,
+                                    .y = y,
+                                    .z = z,
+                                    .r = r,
+                                    .phi = phi,
+                                    .rPhi = rPhi,
+                                    .h = h,
+                                    .eta = eta,
+                                    .mag = mag});
+           }),
+           py::arg("x") = zeroEnvelope, py::arg("y") = zeroEnvelope,
+           py::arg("z") = zeroEnvelope, py::arg("r") = zeroEnvelope,
+           py::arg("phi") = zeroEnvelope, py::arg("rPhi") = zeroEnvelope,
+           py::arg("h") = zeroEnvelope, py::arg("eta") = zeroEnvelope,
+           py::arg("mag") = zeroEnvelope)
+      .def_static("Zero", &ExtentEnvelope::Zero)
+      .def("__getitem__", [](ExtentEnvelope& self,
+                             BinningValue bValue) { return self[bValue]; })
+      .def("__setitem__", [](ExtentEnvelope& self, BinningValue bValue,
+                             const Envelope& value) { self[bValue] = value; })
+      .def("__str__", [](const ExtentEnvelope& self) {
+        std::array<std::string, numBinningValues()> values;
+
+        std::stringstream ss;
+        for (BinningValue val : allBinningValues()) {
+          ss << val << "=(" << self[val][0] << ", " << self[val][1] << ")";
+          values.at(toUnderlying(val)) = ss.str();
+          ss.str("");
+        }
+
+        ss.str("");
+        ss << "ExtentEnvelope(";
+        ss << boost::algorithm::join(values, ", ");
+        ss << ")";
+        return ss.str();
+      });
+
+  py::class_<Extent>(m, "Extent")
+      .def(py::init<const ExtentEnvelope&>(),
+           py::arg("envelope") = ExtentEnvelope::Zero())
+      .def("range",
+           [](const Acts::Extent& self,
+              Acts::BinningValue bval) -> std::array<ActsScalar, 2> {
+             return {self.min(bval), self.max(bval)};
+           })
+      .def("__str__", &Extent::toString);
+
   {
-    py::class_<Acts::Extent>(m, "Extent")
-        .def(py::init(
-            [](const std::vector<std::tuple<Acts::BinningValue,
-                                            std::array<Acts::ActsScalar, 2u>>>&
-                   franges) {
-              Acts::Extent extent;
-              for (const auto& [bval, frange] : franges) {
-                extent.set(bval, frange[0], frange[1]);
-              }
-              return extent;
-            }))
-        .def("range", [](const Acts::Extent& self, Acts::BinningValue bval) {
-          return std::array<Acts::ActsScalar, 2u>{self.min(bval),
-                                                  self.max(bval)};
-        });
+    auto cylStack = py::class_<CylinderVolumeStack>(m, "CylinderVolumeStack");
+
+    py::enum_<CylinderVolumeStack::AttachmentStrategy>(cylStack,
+                                                       "AttachmentStrategy")
+        .value("Gap", CylinderVolumeStack::AttachmentStrategy::Gap)
+        .value("First", CylinderVolumeStack::AttachmentStrategy::First)
+        .value("Second", CylinderVolumeStack::AttachmentStrategy::Second)
+        .value("Midpoint", CylinderVolumeStack::AttachmentStrategy::Midpoint);
+
+    py::enum_<CylinderVolumeStack::ResizeStrategy>(cylStack, "ResizeStrategy")
+        .value("Gap", CylinderVolumeStack::ResizeStrategy::Gap)
+        .value("Expand", CylinderVolumeStack::ResizeStrategy::Expand);
   }
 }
 
@@ -307,10 +381,15 @@ void addExperimentalGeometry(Context& ctx) {
     // Be able to construct a proto binning
     py::class_<ProtoBinning>(m, "ProtoBinning")
         .def(py::init<Acts::BinningValue, Acts::AxisBoundaryType,
-                      const std::vector<Acts::ActsScalar>&, std::size_t>())
+                      const std::vector<Acts::ActsScalar>&, std::size_t>(),
+             "bValue"_a, "bType"_a, "e"_a, "exp"_a = 0u)
         .def(py::init<Acts::BinningValue, Acts::AxisBoundaryType,
                       Acts::ActsScalar, Acts::ActsScalar, std::size_t,
-                      std::size_t>());
+                      std::size_t>(),
+             "bValue"_a, "bType"_a, "minE"_a, "maxE"_a, "nbins"_a, "exp"_a = 0u)
+        .def(py::init<Acts::BinningValue, Acts::AxisBoundaryType, std::size_t,
+                      std::size_t>(),
+             "bValue"_a, "bType"_a, "nbins"_a, "exp"_a = 0u);
   }
 
   {
diff --git a/Examples/Python/src/Input.cpp b/Examples/Python/src/Input.cpp
index e876efaa69c..b3ad9e1f657 100644
--- a/Examples/Python/src/Input.cpp
+++ b/Examples/Python/src/Input.cpp
@@ -63,8 +63,8 @@ void addInput(Context& ctx) {
 
   ACTS_PYTHON_DECLARE_READER(
       ActsExamples::CsvMeasurementReader, mex, "CsvMeasurementReader", inputDir,
-      outputMeasurements, outputMeasurementSimHitsMap, outputSourceLinks,
-      outputClusters, outputMeasurementParticlesMap, inputSimHits);
+      outputMeasurements, outputMeasurementSimHitsMap, outputClusters,
+      outputMeasurementParticlesMap, inputSimHits);
 
   ACTS_PYTHON_DECLARE_READER(ActsExamples::CsvSimHitReader, mex,
                              "CsvSimHitReader", inputDir, inputStem,
diff --git a/Examples/Python/src/TrackFinding.cpp b/Examples/Python/src/TrackFinding.cpp
index 4179d58d94a..ad3ad364ce7 100644
--- a/Examples/Python/src/TrackFinding.cpp
+++ b/Examples/Python/src/TrackFinding.cpp
@@ -58,10 +58,9 @@ namespace Acts::Python {
 void addTrackFinding(Context& ctx) {
   auto [m, mex] = ctx.get("main", "examples");
 
-  ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::SpacePointMaker, mex,
-                                "SpacePointMaker", inputSourceLinks,
-                                inputMeasurements, outputSpacePoints,
-                                trackingGeometry, geometrySelection);
+  ACTS_PYTHON_DECLARE_ALGORITHM(
+      ActsExamples::SpacePointMaker, mex, "SpacePointMaker", inputMeasurements,
+      outputSpacePoints, trackingGeometry, geometrySelection);
 
   {
     using Config = Acts::SeedFilterConfig;
@@ -274,14 +273,14 @@ void addTrackFinding(Context& ctx) {
   ACTS_PYTHON_DECLARE_ALGORITHM(
       ActsExamples::GbtsSeedingAlgorithm, mex, "GbtsSeedingAlgorithm",
       inputSpacePoints, outputSeeds, seedFilterConfig, seedFinderConfig,
-      seedFinderOptions, layerMappingFile, geometrySelection, inputSourceLinks,
-      trackingGeometry, ActsGbtsMap, fill_module_csv, inputClusters);
+      seedFinderOptions, layerMappingFile, geometrySelection, trackingGeometry,
+      ActsGbtsMap, fill_module_csv, inputClusters);
 
   ACTS_PYTHON_DECLARE_ALGORITHM(
       ActsExamples::HoughTransformSeeder, mex, "HoughTransformSeeder",
-      inputSpacePoints, outputProtoTracks, inputSourceLinks, trackingGeometry,
-      geometrySelection, inputMeasurements, subRegions, nLayers, xMin, xMax,
-      yMin, yMax, houghHistSize_x, houghHistSize_y, hitExtend_x, threshold,
+      inputSpacePoints, outputProtoTracks, trackingGeometry, geometrySelection,
+      inputMeasurements, subRegions, nLayers, xMin, xMax, yMin, yMax,
+      houghHistSize_x, houghHistSize_y, hitExtend_x, threshold,
       localMaxWindowSize, kA);
 
   ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::MuonHoughSeeder, mex,
@@ -323,7 +322,6 @@ void addTrackFinding(Context& ctx) {
     auto c = py::class_<Config>(alg, "Config").def(py::init<>());
     ACTS_PYTHON_STRUCT_BEGIN(c, Config);
     ACTS_PYTHON_MEMBER(inputMeasurements);
-    ACTS_PYTHON_MEMBER(inputSourceLinks);
     ACTS_PYTHON_MEMBER(inputInitialTrackParameters);
     ACTS_PYTHON_MEMBER(inputSeeds);
     ACTS_PYTHON_MEMBER(outputTracks);
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/src/Utilities.cpp b/Examples/Python/src/Utilities.cpp
index 7a68a6a3bdb..8993da1c46a 100644
--- a/Examples/Python/src/Utilities.cpp
+++ b/Examples/Python/src/Utilities.cpp
@@ -51,7 +51,7 @@ void addUtilities(Context& ctx) {
 
   ACTS_PYTHON_DECLARE_ALGORITHM(
       ActsExamples::MeasurementMapSelector, mex, "MeasurementMapSelector",
-      inputMeasurementParticleMap, inputSourceLinks,
+      inputMeasurements, inputMeasurementParticleMap,
       outputMeasurementParticleMap, geometrySelection);
 
   ACTS_PYTHON_DECLARE_ALGORITHM(ActsExamples::PrototracksToTracks, mex,
diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt
index 5eb50dee1b7..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,8 +74,8 @@ 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: 75cac3d6c0d9ce02169637cff739ab1e72b00c8fd3db9e99945d40b34c35b70b
-test_refitting[odd]__tracksummary_gsf_refit.root: 90b9c9b3e24441066376f8cf529df10119b31f83563f001ac438fa70874dcb16
+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
diff --git a/Examples/Python/tests/test_examples.py b/Examples/Python/tests/test_examples.py
index 46044a9513c..d5bb593ca9f 100644
--- a/Examples/Python/tests/test_examples.py
+++ b/Examples/Python/tests/test_examples.py
@@ -699,6 +699,7 @@ def test_refitting(tmp_path, detector_config, assert_root_hash):
         runRefittingGsf(
             trackingGeometry=trackingGeometry,
             field=field,
+            digiConfigFile=detector_config.digiConfigFile,
             outputDir=tmp_path,
             s=seq,
         ).run()
diff --git a/Examples/Python/tests/test_reader.py b/Examples/Python/tests/test_reader.py
index f4e803776d5..b83b1317cc8 100644
--- a/Examples/Python/tests/test_reader.py
+++ b/Examples/Python/tests/test_reader.py
@@ -206,7 +206,6 @@ def test_csv_meas_reader(tmp_path, fatras, trk_geo, conf_const):
             level=acts.logging.WARNING,
             outputMeasurements="measurements",
             outputMeasurementSimHitsMap="simhitsmap",
-            outputSourceLinks="sourcelinks",
             outputMeasurementParticlesMap="meas_ptcl_map",
             inputSimHits=simAlg.config.outputSimHits,
             inputDir=str(out),
@@ -215,7 +214,7 @@ def test_csv_meas_reader(tmp_path, fatras, trk_geo, conf_const):
 
     algs = [
         AssertCollectionExistsAlg(k, f"check_alg_{k}", acts.logging.WARNING)
-        for k in ("measurements", "simhitsmap", "sourcelinks", "meas_ptcl_map")
+        for k in ("measurements", "simhitsmap", "meas_ptcl_map")
     ]
     for alg in algs:
         s.addAlgorithm(alg)
@@ -291,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)
 
@@ -363,14 +365,13 @@ def test_edm4hep_measurement_reader(tmp_path, fatras, conf_const):
             level=acts.logging.WARNING,
             outputMeasurements="measurements",
             outputMeasurementSimHitsMap="simhitsmap",
-            outputSourceLinks="sourcelinks",
             inputPath=str(out),
         )
     )
 
     algs = [
         AssertCollectionExistsAlg(k, f"check_alg_{k}", acts.logging.WARNING)
-        for k in ("measurements", "simhitsmap", "sourcelinks")
+        for k in ("measurements", "simhitsmap")
     ]
     for alg in algs:
         s.addAlgorithm(alg)
diff --git a/Examples/Scripts/Python/truth_tracking_gsf_refitting.py b/Examples/Scripts/Python/truth_tracking_gsf_refitting.py
index 138540db2af..e8b41c9aae8 100755
--- a/Examples/Scripts/Python/truth_tracking_gsf_refitting.py
+++ b/Examples/Scripts/Python/truth_tracking_gsf_refitting.py
@@ -11,25 +11,27 @@
 
 
 def runRefittingGsf(
-    trackingGeometry,
-    field,
-    outputDir,
+    trackingGeometry: acts.TrackingGeometry,
+    field: acts.MagneticFieldProvider,
+    digiConfigFile: Path,
+    outputDir: Path,
     s: acts.examples.Sequencer = None,
 ):
-    srcdir = Path(__file__).resolve().parent.parent.parent.parent
-
     s = runTruthTrackingKalman(
         trackingGeometry,
         field,
-        digiConfigFile=srcdir
-        / "Examples/Algorithms/Digitization/share/default-smearing-config-generic.json",
-        # "thirdparty/OpenDataDetector/config/odd-digi-smearing-config.json",
+        digiConfigFile=digiConfigFile,
         outputDir=outputDir,
         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,
@@ -41,7 +43,7 @@ def runRefittingGsf(
         acts.examples.RefittingAlgorithm(
             acts.logging.INFO,
             inputTracks="kf_tracks",
-            outputTracks="gsf_tracks",
+            outputTracks="gsf_refit_tracks",
             fit=acts.examples.makeGsfFitterFunction(
                 trackingGeometry, field, **gsfOptions
             ),
@@ -51,7 +53,7 @@ def runRefittingGsf(
     s.addAlgorithm(
         acts.examples.TrackTruthMatcher(
             level=acts.logging.INFO,
-            inputTracks="gsf_tracks",
+            inputTracks="gsf_refit_tracks",
             inputParticles="truth_seeds_selected",
             inputMeasurementParticlesMap="measurement_particles_map",
             outputTrackParticleMatching="refit_track_particle_matching",
@@ -62,7 +64,7 @@ def runRefittingGsf(
     s.addWriter(
         acts.examples.RootTrackStatesWriter(
             level=acts.logging.INFO,
-            inputTracks="gsf_tracks",
+            inputTracks="gsf_refit_tracks",
             inputParticles="truth_seeds_selected",
             inputTrackParticleMatching="refit_track_particle_matching",
             inputSimHits="simhits",
@@ -87,7 +89,7 @@ def runRefittingGsf(
             inputTracks="tracks",
             inputParticles="truth_seeds_selected",
             inputTrackParticleMatching="track_particle_matching",
-            filePath=str(outputDir / "performance_refitter.root"),
+            filePath=str(outputDir / "performance_gsf_refit.root"),
         )
     )
 
@@ -95,10 +97,28 @@ def runRefittingGsf(
 
 
 if __name__ == "__main__":
-    outputDir = Path.cwd()
+    srcdir = Path(__file__).resolve().parent.parent.parent.parent
+
+    # ODD
+    from acts.examples.odd import getOpenDataDetector
+
+    detector, trackingGeometry, decorators = getOpenDataDetector()
+    digiConfigFile = (
+        srcdir / "thirdparty/OpenDataDetector/config/odd-digi-smearing-config.json"
+    )
+
+    ## GenericDetector
+    # detector, trackingGeometry, _ = acts.examples.GenericDetector.create()
+    # digiConfigFile = (
+    #     srcdir
+    #     / "Examples/Algorithms/Digitization/share/default-smearing-config-generic.json"
+    # )
 
-    # detector, trackingGeometry, decorators = getOpenDataDetector()
-    detector, trackingGeometry, decorators = acts.examples.GenericDetector.create()
     field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T))
 
-    runRefittingGsf(trackingGeometry, field, outputDir).run()
+    runRefittingGsf(
+        trackingGeometry=trackingGeometry,
+        field=field,
+        digiConfigFile=digiConfigFile,
+        outputDir=Path.cwd(),
+    ).run()
diff --git a/Examples/Scripts/Python/truth_tracking_kalman_refitting.py b/Examples/Scripts/Python/truth_tracking_kalman_refitting.py
new file mode 100755
index 00000000000..7e795b9e310
--- /dev/null
+++ b/Examples/Scripts/Python/truth_tracking_kalman_refitting.py
@@ -0,0 +1,122 @@
+#!/usr/bin/env python3
+
+from pathlib import Path
+
+import acts
+import acts.examples
+
+from truth_tracking_kalman import runTruthTrackingKalman
+
+u = acts.UnitConstants
+
+
+def runRefittingKf(
+    trackingGeometry: acts.TrackingGeometry,
+    field: acts.MagneticFieldProvider,
+    digiConfigFile: Path,
+    outputDir: Path,
+    multipleScattering: bool = True,
+    energyLoss: bool = True,
+    reverseFilteringMomThreshold=0 * u.GeV,
+    s: acts.examples.Sequencer = None,
+):
+    s = runTruthTrackingKalman(
+        trackingGeometry,
+        field,
+        digiConfigFile=digiConfigFile,
+        outputDir=outputDir,
+        s=s,
+    )
+
+    kalmanOptions = {
+        "multipleScattering": multipleScattering,
+        "energyLoss": energyLoss,
+        "reverseFilteringMomThreshold": reverseFilteringMomThreshold,
+        "freeToBoundCorrection": acts.examples.FreeToBoundCorrection(False),
+        "level": acts.logging.INFO,
+    }
+
+    s.addAlgorithm(
+        acts.examples.RefittingAlgorithm(
+            level=acts.logging.INFO,
+            inputTracks="kf_tracks",
+            outputTracks="kf_refit_tracks",
+            fit=acts.examples.makeKalmanFitterFunction(
+                trackingGeometry, field, **kalmanOptions
+            ),
+        )
+    )
+
+    s.addAlgorithm(
+        acts.examples.TrackTruthMatcher(
+            level=acts.logging.INFO,
+            inputTracks="kf_refit_tracks",
+            inputParticles="truth_seeds_selected",
+            inputMeasurementParticlesMap="measurement_particles_map",
+            outputTrackParticleMatching="refit_track_particle_matching",
+            outputParticleTrackMatching="refit_particle_track_matching",
+        )
+    )
+
+    s.addWriter(
+        acts.examples.RootTrackStatesWriter(
+            level=acts.logging.INFO,
+            inputTracks="kf_refit_tracks",
+            inputParticles="truth_seeds_selected",
+            inputTrackParticleMatching="refit_track_particle_matching",
+            inputSimHits="simhits",
+            inputMeasurementSimHitsMap="measurement_simhits_map",
+            filePath=str(outputDir / "trackstates_kf_refit.root"),
+        )
+    )
+
+    s.addWriter(
+        acts.examples.RootTrackSummaryWriter(
+            level=acts.logging.INFO,
+            inputTracks="tracks",
+            inputParticles="truth_seeds_selected",
+            inputTrackParticleMatching="refit_track_particle_matching",
+            filePath=str(outputDir / "tracksummary_kf_refit.root"),
+        )
+    )
+
+    s.addWriter(
+        acts.examples.TrackFitterPerformanceWriter(
+            level=acts.logging.INFO,
+            inputTracks="tracks",
+            inputParticles="truth_seeds_selected",
+            inputTrackParticleMatching="track_particle_matching",
+            filePath=str(outputDir / "performance_kf_refit.root"),
+        )
+    )
+
+    return s
+
+
+if __name__ == "__main__":
+    srcdir = Path(__file__).resolve().parent.parent.parent.parent
+    outputDir = Path.cwd()
+
+    # ODD
+    from acts.examples.odd import getOpenDataDetector
+
+    detector, trackingGeometry, decorators = getOpenDataDetector()
+    digiConfigFile = (
+        srcdir / "thirdparty/OpenDataDetector/config/odd-digi-smearing-config.json"
+    )
+
+    ## GenericDetector
+    # detector, trackingGeometry, _ = acts.examples.GenericDetector.create()
+    # digiConfigFile = (
+    #     srcdir
+    #     / "Examples/Algorithms/Digitization/share/default-smearing-config-generic.json"
+    # )
+
+    field = acts.ConstantBField(acts.Vector3(0, 0, 2 * u.T))
+
+    runRefittingKf(
+        trackingGeometry=trackingGeometry,
+        field=field,
+        digiConfigFile=digiConfigFile,
+        outputDir=Path.cwd(),
+    ).run()
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/GeoModel/include/Acts/Plugins/GeoModel/GeoModelBlueprintCreater.hpp b/Plugins/GeoModel/include/Acts/Plugins/GeoModel/GeoModelBlueprintCreater.hpp
index df895eb0741..c6358bdbec2 100644
--- a/Plugins/GeoModel/include/Acts/Plugins/GeoModel/GeoModelBlueprintCreater.hpp
+++ b/Plugins/GeoModel/include/Acts/Plugins/GeoModel/GeoModelBlueprintCreater.hpp
@@ -58,6 +58,8 @@ class GeoModelBlueprintCreater {
     std::string topEntry;
     /// Optionally override the top node bounds
     std::string topBoundsOverride = "";
+    /// Export dot graph
+    std::string dotGraph = "";
   };
 
   /// The Blueprint return object
diff --git a/Plugins/GeoModel/src/GeoModelBlueprintCreater.cpp b/Plugins/GeoModel/src/GeoModelBlueprintCreater.cpp
index 93712d5d911..0b3d2f66e42 100644
--- a/Plugins/GeoModel/src/GeoModelBlueprintCreater.cpp
+++ b/Plugins/GeoModel/src/GeoModelBlueprintCreater.cpp
@@ -10,6 +10,7 @@
 
 #include "Acts/Detector/GeometryIdGenerator.hpp"
 #include "Acts/Detector/LayerStructureBuilder.hpp"
+#include "Acts/Detector/detail/BlueprintDrawer.hpp"
 #include "Acts/Detector/detail/BlueprintHelper.hpp"
 #include "Acts/Detector/interface/IGeometryIdGenerator.hpp"
 #include "Acts/Plugins/GeoModel/GeoModelTree.hpp"
@@ -20,6 +21,8 @@
 #include "Acts/Utilities/Helpers.hpp"
 #include "Acts/Utilities/RangeXD.hpp"
 
+#include <fstream>
+
 #include <boost/algorithm/string.hpp>
 
 using namespace Acts::detail;
@@ -128,6 +131,14 @@ Acts::GeoModelBlueprintCreater::create(const GeometryContext& gctx,
   blueprint.topNode =
       createNode(cache, gctx, topEntry->second, blueprintTableMap, Extent());
 
+  // Export to dot graph if configured
+  if (!options.dotGraph.empty()) {
+    std::ofstream dotFile(options.dotGraph);
+    Experimental::detail::BlueprintDrawer::dotStream(dotFile,
+                                                     *blueprint.topNode);
+    dotFile.close();
+  }
+
   // Return the ready-to-use blueprint
   return blueprint;
 }
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/CommonHelpers/Acts/Tests/CommonHelpers/WhiteBoardUtilities.hpp b/Tests/CommonHelpers/Acts/Tests/CommonHelpers/WhiteBoardUtilities.hpp
index f40f658e63e..85a760bd6af 100644
--- a/Tests/CommonHelpers/Acts/Tests/CommonHelpers/WhiteBoardUtilities.hpp
+++ b/Tests/CommonHelpers/Acts/Tests/CommonHelpers/WhiteBoardUtilities.hpp
@@ -15,6 +15,7 @@
 #include <tuple>
 
 namespace Acts::Test {
+
 struct DummySequenceElement : public ActsExamples::SequenceElement {
   ActsExamples::ProcessCode initialize() override { return {}; };
   ActsExamples::ProcessCode finalize() override { return {}; };
@@ -102,4 +103,5 @@ struct GenericReadWriteTool {
     return get(get, std::tuple<>{}, std::integral_constant<std::size_t, 0>{});
   }
 };
+
 }  // namespace Acts::Test
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/Core/Geometry/PortalShellTests.cpp b/Tests/UnitTests/Core/Geometry/PortalShellTests.cpp
index 4c44617524a..15c5a9c7953 100644
--- a/Tests/UnitTests/Core/Geometry/PortalShellTests.cpp
+++ b/Tests/UnitTests/Core/Geometry/PortalShellTests.cpp
@@ -68,7 +68,7 @@ BOOST_AUTO_TEST_CASE(ConstructionFromVolume) {
   SingleCylinderPortalShell shell1{cyl1};
   BOOST_CHECK_EQUAL(shell1.size(), 4);
 
-  using enum CylinderPortalShell::Face;
+  using enum CylinderVolumeBounds::Face;
 
   const auto* pDisc = shell1.portal(PositiveDisc);
   BOOST_REQUIRE_NE(pDisc, nullptr);
@@ -299,7 +299,7 @@ BOOST_AUTO_TEST_CASE(ConstructionFromVolume) {
 //              inner cylinder
 
 BOOST_AUTO_TEST_CASE(PortalAssignment) {
-  using enum CylinderPortalShell::Face;
+  using enum CylinderVolumeBounds::Face;
   TrackingVolume vol(
       Transform3::Identity(),
       std::make_shared<CylinderVolumeBounds>(30_mm, 100_mm, 100_mm));
@@ -350,7 +350,7 @@ BOOST_AUTO_TEST_CASE(PortalAssignment) {
 
 BOOST_AUTO_TEST_SUITE(CylinderStack)
 BOOST_AUTO_TEST_CASE(ZDirection) {
-  using enum CylinderPortalShell::Face;
+  using enum CylinderVolumeBounds::Face;
   BOOST_TEST_CONTEXT("rMin>0") {
     TrackingVolume vol1(
         Transform3{Translation3{Vector3::UnitZ() * -100_mm}},
@@ -447,7 +447,7 @@ BOOST_AUTO_TEST_CASE(ZDirection) {
 }
 
 BOOST_AUTO_TEST_CASE(RDirection) {
-  using enum CylinderPortalShell::Face;
+  using enum CylinderVolumeBounds::Face;
   BOOST_TEST_CONTEXT("rMin>0") {
     TrackingVolume vol1(
         Transform3::Identity(),
@@ -610,7 +610,7 @@ BOOST_AUTO_TEST_CASE(NestedStacks) {
       gctx, {&stack, &shell3}, BinningValue::binZ, *logger};
   BOOST_CHECK(stack2.isValid());
 
-  using enum CylinderPortalShell::Face;
+  using enum CylinderVolumeBounds::Face;
 
   auto lookup = [](auto& shell, CylinderPortalShell::Face face,
                    Vector3 position,
@@ -726,7 +726,7 @@ BOOST_AUTO_TEST_CASE(ConnectOuter) {
 
   SingleCylinderPortalShell shell{cyl1};
 
-  using enum CylinderPortalShell::Face;
+  using enum CylinderVolumeBounds::Face;
   BOOST_CHECK_EQUAL(
       shell.portal(OuterCylinder)->getLink(Direction::AlongNormal), nullptr);
   BOOST_CHECK_EQUAL(
@@ -749,7 +749,7 @@ BOOST_AUTO_TEST_CASE(ConnectOuter) {
 }
 
 BOOST_AUTO_TEST_CASE(RegisterInto) {
-  using enum CylinderPortalShell::Face;
+  using enum CylinderVolumeBounds::Face;
   TrackingVolume vol1(
       Transform3::Identity(),
       std::make_shared<CylinderVolumeBounds>(0_mm, 100_mm, 100_mm));
diff --git a/Tests/UnitTests/Examples/EventData/MeasurementTests.cpp b/Tests/UnitTests/Examples/EventData/MeasurementTests.cpp
index fc3d33d801f..d59ff4a135e 100644
--- a/Tests/UnitTests/Examples/EventData/MeasurementTests.cpp
+++ b/Tests/UnitTests/Examples/EventData/MeasurementTests.cpp
@@ -11,9 +11,8 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/TrackParametrization.hpp"
-#include "Acts/EventData/SourceLink.hpp"
 #include "Acts/EventData/detail/GenerateParameters.hpp"
-#include "Acts/EventData/detail/TestSourceLink.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
 #include "ActsExamples/EventData/Measurement.hpp"
 
@@ -30,15 +29,13 @@
 using namespace Acts;
 using namespace Acts::detail::Test;
 using namespace ActsExamples;
-using SourceLink = Acts::detail::Test::TestSourceLink;
 namespace bd = boost::unit_test::data;
 
 namespace {
 constexpr BoundIndices boundIndices[] = {
     eBoundLoc0, eBoundLoc1, eBoundTime, eBoundPhi, eBoundTheta, eBoundQOverP,
 };
-const TestSourceLink sourceOrig;
-const Acts::SourceLink source{sourceOrig};
+constexpr Acts::GeometryIdentifier geoId = 1;
 // fix seed for reproducible tests
 std::default_random_engine rng(123);
 }  // namespace
@@ -54,8 +51,7 @@ BOOST_DATA_TEST_CASE(VariableBoundOne, bd::make(boundIndices), index) {
 
   auto [params, cov] = generateParametersCovariance<ActsScalar, 1u>(rng);
 
-  FixedBoundMeasurementProxy<1> meas = container.makeMeasurement<1>();
-  meas.setSourceLink(source);
+  FixedBoundMeasurementProxy<1> meas = container.makeMeasurement<1>(geoId);
   meas.setSubspaceIndices(std::array{index});
   meas.parameters() = params;
   meas.covariance() = cov;
@@ -66,8 +62,7 @@ BOOST_DATA_TEST_CASE(VariableBoundOne, bd::make(boundIndices), index) {
   }
   BOOST_CHECK_EQUAL(meas.parameters(), params);
   BOOST_CHECK_EQUAL(meas.covariance(), cov);
-  BOOST_CHECK_EQUAL(meas.sourceLink().template get<TestSourceLink>(),
-                    sourceOrig);
+  BOOST_CHECK_EQUAL(meas.geometryId(), geoId);
 }
 
 BOOST_DATA_TEST_CASE(VariableBoundOneEmplace, bd::make(boundIndices), index) {
@@ -76,7 +71,7 @@ BOOST_DATA_TEST_CASE(VariableBoundOneEmplace, bd::make(boundIndices), index) {
   auto [params, cov] = generateParametersCovariance<ActsScalar, 1u>(rng);
 
   FixedBoundMeasurementProxy<1> meas =
-      container.emplaceMeasurement<1>(source, std::array{index}, params, cov);
+      container.emplaceMeasurement<1>(geoId, std::array{index}, params, cov);
 
   BOOST_CHECK_EQUAL(meas.size(), 1);
   for (auto i : boundIndices) {
@@ -84,8 +79,7 @@ BOOST_DATA_TEST_CASE(VariableBoundOneEmplace, bd::make(boundIndices), index) {
   }
   BOOST_CHECK_EQUAL(meas.parameters(), params);
   BOOST_CHECK_EQUAL(meas.covariance(), cov);
-  BOOST_CHECK_EQUAL(meas.sourceLink().template get<TestSourceLink>(),
-                    sourceOrig);
+  BOOST_CHECK_EQUAL(meas.geometryId(), geoId);
 }
 
 BOOST_AUTO_TEST_CASE(VariableBoundAll) {
@@ -94,8 +88,7 @@ BOOST_AUTO_TEST_CASE(VariableBoundAll) {
   auto [params, cov] = generateBoundParametersCovariance(rng);
 
   FixedBoundMeasurementProxy<eBoundSize> meas =
-      container.makeMeasurement<eBoundSize>();
-  meas.setSourceLink(source);
+      container.makeMeasurement<eBoundSize>(geoId);
   meas.setSubspaceIndices(std::array{eBoundLoc0, eBoundLoc1, eBoundTime,
                                      eBoundPhi, eBoundTheta, eBoundQOverP});
   meas.parameters() = params;
@@ -107,7 +100,7 @@ BOOST_AUTO_TEST_CASE(VariableBoundAll) {
   }
   BOOST_CHECK_EQUAL(meas.parameters(), params);
   BOOST_CHECK_EQUAL(meas.covariance(), cov);
-  BOOST_CHECK_EQUAL(meas.sourceLink().get<TestSourceLink>(), sourceOrig);
+  BOOST_CHECK_EQUAL(meas.geometryId(), geoId);
 }
 
 BOOST_AUTO_TEST_CASE(VariableBoundAllEmplace) {
@@ -117,7 +110,7 @@ BOOST_AUTO_TEST_CASE(VariableBoundAllEmplace) {
 
   FixedBoundMeasurementProxy<eBoundSize> meas =
       container.emplaceMeasurement<eBoundSize>(
-          source,
+          geoId,
           std::array{eBoundLoc0, eBoundLoc1, eBoundTime, eBoundPhi, eBoundTheta,
                      eBoundQOverP},
           params, cov);
@@ -128,7 +121,7 @@ BOOST_AUTO_TEST_CASE(VariableBoundAllEmplace) {
   }
   BOOST_CHECK_EQUAL(meas.parameters(), params);
   BOOST_CHECK_EQUAL(meas.covariance(), cov);
-  BOOST_CHECK_EQUAL(meas.sourceLink().get<TestSourceLink>(), sourceOrig);
+  BOOST_CHECK_EQUAL(meas.geometryId(), geoId);
 }
 
 BOOST_AUTO_TEST_CASE(VariableBoundReassign) {
@@ -137,8 +130,7 @@ BOOST_AUTO_TEST_CASE(VariableBoundReassign) {
   // generate w/ two parameter
   auto [params1, cov1] = generateParametersCovariance<ActsScalar, 2u>(rng);
 
-  VariableBoundMeasurementProxy meas = container.makeMeasurement(2);
-  meas.setSourceLink(source);
+  VariableBoundMeasurementProxy meas = container.makeMeasurement(2, geoId);
   meas.setSubspaceIndices(std::array{eBoundPhi, eBoundTheta});
   meas.parameters() = params1;
   meas.covariance() = cov1;
@@ -154,8 +146,7 @@ BOOST_AUTO_TEST_CASE(VariableBoundReassign) {
   // reassign w/ all parameters
   auto [paramsN, covN] = generateBoundParametersCovariance(rng);
 
-  meas = container.makeMeasurement(eBoundSize);
-  meas.setSourceLink(source);
+  meas = container.makeMeasurement(eBoundSize, geoId);
   meas.setSubspaceIndices(std::array{eBoundLoc0, eBoundLoc1, eBoundTime,
                                      eBoundPhi, eBoundTheta, eBoundQOverP});
   meas.parameters() = paramsN;
diff --git a/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp b/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp
index 7786a5ec8a1..c9c21b0983a 100644
--- a/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp
+++ b/Tests/UnitTests/Examples/Io/Csv/MeasurementReaderWriterTests.cpp
@@ -29,7 +29,6 @@ using namespace ActsExamples;
 using namespace Acts::Test;
 
 BOOST_AUTO_TEST_CASE(CsvMeasurementRoundTrip) {
-  IndexSourceLinkContainer sourceLinksOriginal;
   MeasurementContainer measOriginal;
   ClusterContainer clusterOriginal;
   IndexMultimap<Index> mapOriginal;
@@ -45,14 +44,11 @@ BOOST_AUTO_TEST_CASE(CsvMeasurementRoundTrip) {
   std::uniform_real_distribution<double> distf(0.0, 1.0);
 
   for (auto i = 0ul; i < nMeasurements; ++i) {
-    IndexSourceLink sl(someGeoId, static_cast<Index>(i));
-    sourceLinksOriginal.insert(sl);
-
     Acts::Vector2 p = Acts::Vector2::Random();
     Acts::SquareMatrix2 c = Acts::SquareMatrix2::Random();
 
-    FixedBoundMeasurementProxy<2> m = measOriginal.makeMeasurement<2>();
-    m.setSourceLink(Acts::SourceLink(sl));
+    FixedBoundMeasurementProxy<2> m =
+        measOriginal.makeMeasurement<2>(someGeoId);
     m.setSubspaceIndices(std::array{Acts::eBoundLoc0, Acts::eBoundLoc1});
     m.parameters() = p;
     m.covariance() = c;
@@ -118,14 +114,12 @@ BOOST_AUTO_TEST_CASE(CsvMeasurementRoundTrip) {
   readerConfig.outputMeasurementSimHitsMap =
       writerConfig.inputMeasurementSimHitsMap;
   readerConfig.outputClusters = writerConfig.inputClusters;
-  readerConfig.outputSourceLinks = "sourcelinks";
 
   CsvMeasurementReader reader(readerConfig, Acts::Logging::WARNING);
 
-  auto readTool =
-      writeTool.add(readerConfig.outputSourceLinks, sourceLinksOriginal);
+  auto readTool = writeTool.add(readerConfig.outputMeasurements, measOriginal);
 
-  const auto [measRead, clusterRead, mapRead, sourceLinksRead] =
+  const auto [measRead, clusterRead, mapRead, measRead2] =
       readTool.read(reader);
 
   ///////////
@@ -164,10 +158,10 @@ BOOST_AUTO_TEST_CASE(CsvMeasurementRoundTrip) {
     BOOST_REQUIRE(a == b);
   }
 
-  static_assert(std::is_same_v<std::decay_t<decltype(sourceLinksRead)>,
-                               decltype(sourceLinksOriginal)>);
-  BOOST_REQUIRE(sourceLinksRead.size() == sourceLinksOriginal.size());
-  for (const auto &[a, b] : Acts::zip(sourceLinksRead, sourceLinksOriginal)) {
+  static_assert(
+      std::is_same_v<std::decay_t<decltype(measRead)>, decltype(measOriginal)>);
+  BOOST_REQUIRE(measRead.size() == measOriginal.size());
+  for (const auto &[a, b] : Acts::zip(measRead, measOriginal)) {
     BOOST_REQUIRE(a.geometryId() == b.geometryId());
     BOOST_REQUIRE(a.index() == b.index());
   }
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()
diff --git a/docs/getting_started.md b/docs/getting_started.md
index fbc247a3e05..f17ec5abcf1 100644
--- a/docs/getting_started.md
+++ b/docs/getting_started.md
@@ -9,9 +9,9 @@ following commands will clone the repository, configure, and build the core
 library:
 
 ```console
-$ git clone https://github.com/acts-project/acts <source>
-$ cmake -B <build> -S <source>
-$ cmake --build <build>
+git clone https://github.com/acts-project/acts <source>
+cmake -B <build> -S <source>
+cmake --build <build>
 ```
 
 For a full list of dependencies, including specific versions, see the
@@ -23,27 +23,27 @@ section.
 
 The following dependencies are required to build the ACTS core library:
 
--   A C++17 compatible compiler (recent versions of either gcc and clang should work)
--   [CMake](https://cmake.org) >= 3.14
--   [Boost](https://www.boost.org) >= 1.71 with `filesystem`, `program_options`, and `unit_test_framework`
--   [Eigen](https://eigen.tuxfamily.org) >= 3.3.7
+- A C++17 compatible compiler (recent versions of either gcc and clang should work)
+- [CMake](https://cmake.org) >= 3.14
+- [Boost](https://www.boost.org) >= 1.71 with `filesystem`, `program_options`, and `unit_test_framework`
+- [Eigen](https://eigen.tuxfamily.org) >= 3.3.7
 
 The following dependencies are optional and are needed to build additional
 components:
 
--   [CUDA](https://developer.nvidia.com/cuda-zone) for the CUDA plugin and the Exa.TrkX plugin and its examples
--   [DD4hep](http://dd4hep.cern.ch) >= 1.11 for the DD4hep plugin and some examples
--   [Doxygen](http://doxygen.org) >= 1.8.15 for the documentation
--   [Geant4](https://geant4.org/) for some examples
--   [HepMC](https://gitlab.cern.ch/hepmc/HepMC3) >= 3.2.1 for some examples
--   [Intel Threading Building Blocks](https://github.com/oneapi-src/oneTBB) >= 2020.1 for the examples
--   [ONNX Runtime](https://onnxruntime.ai/) >= 1.12.0 for the ONNX plugin, the Exa.TrkX plugin and some examples
--   [Pythia8](https://pythia.org) for some examples
--   [ROOT](https://root.cern.ch) >= 6.20 for the TGeo plugin and the examples
--   [Sphinx](https://www.sphinx-doc.org) >= 2.0 with [Breathe](https://breathe.readthedocs.io/en/latest/), [Exhale](https://exhale.readthedocs.io/en/latest/), and [recommonmark](https://recommonmark.readthedocs.io/en/latest/index.html) extensions for the documentation
--   [cugraph](https://github.com/rapidsai/cugraph) for the Exa.TrkX plugin
--   [libtorch](https://pytorch.org/cppdocs/installing.html) for the Exa.TrkX plugin
--   [Pybind11](https://github.com/pybind/pybind11) for the Python bindings of the examples
+- [CUDA](https://developer.nvidia.com/cuda-zone) for the CUDA plugin and the Exa.TrkX plugin and its examples
+- [DD4hep](http://dd4hep.cern.ch) >= 1.11 for the DD4hep plugin and some examples
+- [Doxygen](http://doxygen.org) >= 1.8.15 for the documentation
+- [Geant4](https://geant4.org/) for some examples
+- [HepMC](https://gitlab.cern.ch/hepmc/HepMC3) >= 3.2.1 for some examples
+- [Intel Threading Building Blocks](https://github.com/oneapi-src/oneTBB) >= 2020.1 for the examples
+- [ONNX Runtime](https://onnxruntime.ai/) >= 1.12.0 for the ONNX plugin, the Exa.TrkX plugin and some examples
+- [Pythia8](https://pythia.org) for some examples
+- [ROOT](https://root.cern.ch) >= 6.20 for the TGeo plugin and the examples
+- [Sphinx](https://www.sphinx-doc.org) >= 2.0 with [Breathe](https://breathe.readthedocs.io/en/latest/), [Exhale](https://exhale.readthedocs.io/en/latest/), and [recommonmark](https://recommonmark.readthedocs.io/en/latest/index.html) extensions for the documentation
+- [cugraph](https://github.com/rapidsai/cugraph) for the Exa.TrkX plugin
+- [libtorch](https://pytorch.org/cppdocs/installing.html) for the Exa.TrkX plugin
+- [Pybind11](https://github.com/pybind/pybind11) for the Python bindings of the examples
 
 There are some additional dependencies that are automatically provided as part of
 the build system.
@@ -69,7 +69,7 @@ runs the configuration and searches for the dependencies. The `<build>`
 directory is automatically created.
 
 ```console
-$ cmake -B <build> -S <source>
+cmake -B <build> -S <source>
 ```
 
 The build can be configured via various options that are listed in detail in the
@@ -77,19 +77,19 @@ The build can be configured via various options that are listed in detail in the
 The previous command could be e.g. modified to
 
 ```console
-$ cmake -B <build> -S <source> -DACTS_BUILD_UNITTESTS=on -DACTS_BUILD_FATRAS=on
+cmake -B <build> -S <source> -DACTS_BUILD_UNITTESTS=on -DACTS_BUILD_FATRAS=on
 ```
 
 After the configuration succeeded, the software is build. This is also done with cmake via the following command
 
 ```console
-$ cmake --build <build>
+cmake --build <build>
 ```
 
 This automatically calls the configure build tool, e.g. Make or Ninja. To build only a specific target, the target names has to be separated from the CMake options by `--`, i.e.
 
 ```console
-$ cmake --build <build> -- ActsFatras # to build the Fatras library
+cmake --build <build> -- ActsFatras # to build the Fatras library
 ```
 
 The build commands are the same regardless of where you are building the
@@ -103,8 +103,8 @@ e.g. CERNs lxplus login machines, the dependencies can be easily satisfied
 via a LCG releases available through CVMFS. A setup script is provided to activate a compatible releases that can be used as follows:
 
 ```console
-$ cd <source>
-$ source CI/setup_cvmfs_lcg.sh
+cd <source>
+source CI/setup_cvmfs_lcg.sh
 ```
 
 After sourcing the setup script, you can build ACTS as described above. The
@@ -112,10 +112,10 @@ following commands will build ACTS in the `<source>/build` directory with the
 Fatras component.
 
 ```console
-$ cd <source>
-$ source CI/setup_cvmfs_lcg.sh
-$ cmake -B build -S . -DACTS_BUILD_FATRAS=on
-$ cmake --build build
+cd <source>
+source CI/setup_cvmfs_lcg.sh
+cmake -B build -S . -DACTS_BUILD_FATRAS=on
+cmake --build build
 ```
 
 ### In a container
@@ -144,13 +144,13 @@ available tags, e.g. for the `ubuntu2004` image, you can use the following
 command:
 
 ```console
-$ docker search --list-tags ghcr.io/acts-project/ubuntu2404
+docker search --list-tags ghcr.io/acts-project/ubuntu2404
 ```
 
 The following command then downloads a stable tag of the `ubuntu2404` image:
 
 ```console
-$ docker pull ghcr.io/acts-project/ubuntu2404:51
+docker pull ghcr.io/acts-project/ubuntu2404:51
 ```
 
 This should print the image id as part of the output. You can also find out the
@@ -163,7 +163,7 @@ following command will make the source directory available as `/acts` in the
 container and start an interactive `bash` shell
 
 ```console
-$ docker run --volume=<source>:/acts:ro --interactive --tty <image> /bin/bash
+docker run --volume=<source>:/acts:ro --interactive --tty <image> /bin/bash
 ```
 
 where `<image>` is the image id that was previously mentioned. If you are using the Ubuntu-based image you are already good to go. For the images based on LCG releases, you can now activate the LCG release in the container shell by sourcing a setup script:
@@ -191,6 +191,7 @@ install ACTS' dependencies; see the [building with Spack](misc/spack) page for
 more information.
 
 (build_docs)=
+
 ## Building the documentation
 
 The documentation uses [Doxygen][doxygen] to extract the source code
@@ -201,8 +202,8 @@ need to have [Doxygen][doxygen] version `1.9.5` or newer installed.
 package manager `pip`:
 
 ```console
-$ cd <source>
-$ pip install -r docs/requirements.txt
+cd <source>
+pip install -r docs/requirements.txt
 ```
 
 :::{tip}
@@ -211,8 +212,8 @@ environment](https://realpython.com/python-virtual-environments-a-primer/) for
 this purpose! For example, run
 
 ```console
-$ python -m venv docvenv
-$ source docvenv/bin/activate
+python -m venv docvenv
+source docvenv/bin/activate
 ```
 
 to create a local virtual environment, and then run the `pip` command above.
@@ -221,13 +222,13 @@ to create a local virtual environment, and then run the `pip` command above.
 To activate the documentation build targets, the `ACTS_BUILD_DOCS` option has to be set
 
 ```console
-$ cmake -B <build> -S <source> -DACTS_BUILD_DOCS=on
+cmake -B <build> -S <source> -DACTS_BUILD_DOCS=on
 ```
 
 Then the documentation can be build with this target
 
 ```console
-$ cmake --build <build> --target docs
+cmake --build <build> --target docs
 ```
 
 The default option includes the Doxygen, Sphinx, and the Breathe extension,
@@ -239,8 +240,6 @@ of errors you will need to manually pull in symbols to be documented.
 [doxygen]: https://doxygen.nl/
 [sphinx]: https://www.sphinx-doc.org
 [breathe]: https://breathe.readthedocs.io
-[exhale]: https://exhale.readthedocs.io
-[rtd_acts]: https://acts.readthedocs.io
 
 ## Build options
 
@@ -248,7 +247,7 @@ CMake options can be set by adding `-D<OPTION>=<VALUE>` to the configuration
 command. The following command would e.g. enable the unit tests
 
 ```console
-$ cmake -B <build> -S <source> -DACTS_BUILD_UNITTESTS=ON
+cmake -B <build> -S <source> -DACTS_BUILD_UNITTESTS=ON
 ```
 
 Multiple options can be given. `cmake` caches the options so that only changed
@@ -264,7 +263,6 @@ components.
 | ACTS_FORCE_ASSERTIONS               | Force assertions regardless of build<br>type<br> type: `bool`, default: `OFF`                                                                                                                                                      |
 | ACTS_USE_SYSTEM_LIBS                | Use system libraries by default<br> type: `bool`, default: `OFF`                                                                                                                                                                   |
 | ACTS_USE_SYSTEM_ACTSVG              | Use the ActSVG system library<br> type: `bool`, default: `ACTS_USE_SYSTEM_LIBS -> OFF`                                                                                                                                             |
-| ACTS_USE_SYSTEM_GEOMODEL            | Use a system-provided GeoModel<br>installation<br> type: `bool`, default: `ACTS_USE_SYSTEM_LIBS -> OFF`                                                                                                                            |
 | ACTS_USE_SYSTEM_COVFIE              | Use a system-provided covfie<br>installation<br> type: `bool`, default: `ACTS_USE_SYSTEM_LIBS -> OFF`                                                                                                                              |
 | ACTS_USE_SYSTEM_DETRAY              | Use a system-provided detray<br>installation<br> type: `bool`, default: `ACTS_USE_SYSTEM_LIBS -> OFF`                                                                                                                              |
 | ACTS_USE_SYSTEM_TRACCC              | Use a system-provided traccc<br>installation<br> type: `bool`, default: `ACTS_USE_SYSTEM_LIBS -> OFF`                                                                                                                              |
@@ -350,7 +348,7 @@ documentation](https://cmake.org/documentation/).
 The build is also affected by some environment variables. They can be set by prepending them to the configuration call:
 
 ```console
-$ DD4hep_DIR=<path/to/dd4hep> cmake -B <build> -S <source>
+DD4hep_DIR=<path/to/dd4hep> cmake -B <build> -S <source>
 ```
 
 The following environment variables might be useful.
@@ -367,8 +365,8 @@ ACTS comes packaged with a detector modeled using DD4hep that can be used to tes
 It is available via the git submodule feature by performing the following steps ([`git lfs`](https://git-lfs.com/) need to be installed on your machine):
 
 ```console
-$ git submodule init
-$ git submodule update
+git submodule init
+git submodule update
 ```
 
 To use it, you will then need to build ACTS with the `ACTS_BUILD_ODD` option and then point either `LD_LIBRARY_PATH` on Linux or
@@ -381,7 +379,6 @@ oddMaterialDeco = acts.IMaterialDecorator.fromFile("PATH_TO_Acts/thirdparty/Open
 detector, trackingGeometry, decorators = getOpenDataDetector(oddMaterialDeco)
 ```
 
-
 ## Using ACTS
 
 When using ACTS in your own CMake-based project, you need to include the