diff --git a/Core/include/Acts/MagneticField/MultiRangeBField.hpp b/Core/include/Acts/MagneticField/MultiRangeBField.hpp
new file mode 100644
index 00000000000..60b90097161
--- /dev/null
+++ b/Core/include/Acts/MagneticField/MultiRangeBField.hpp
@@ -0,0 +1,67 @@
+// 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/Definitions/Algebra.hpp"
+#include "Acts/MagneticField/MagneticFieldContext.hpp"
+#include "Acts/MagneticField/MagneticFieldError.hpp"
+#include "Acts/MagneticField/MagneticFieldProvider.hpp"
+#include "Acts/Utilities/RangeXD.hpp"
+
+namespace Acts {
+
+/// @ingroup MagneticField
+///
+/// @brief Magnetic field provider modelling a magnetic field consisting of
+/// several (potentially overlapping) regions of constant values.
+class MultiRangeBField final : public MagneticFieldProvider {
+ private:
+  struct Cache {
+    explicit Cache(const MagneticFieldContext& /*unused*/);
+
+    std::optional<std::size_t> index = {};
+  };
+
+  using BFieldRange = std::pair<RangeXD<3, double>, Vector3>;
+
+  // The different ranges and their corresponding field vectors. Note that
+  // regions positioned _later_ in this vector take priority over earlier
+  // regions.
+  std::vector<BFieldRange> fieldRanges;
+
+ public:
+  /// @brief Construct a magnetic field from a vector of ranges.
+  ///
+  /// @warning These ranges are listed in increasing order of precedence,
+  /// i.e. ranges further along the vector have higher priority.
+  explicit MultiRangeBField(const std::vector<BFieldRange>& ranges);
+
+  explicit MultiRangeBField(std::vector<BFieldRange>&& ranges);
+
+  /// @brief Construct a cache object.
+  MagneticFieldProvider::Cache makeCache(
+      const MagneticFieldContext& mctx) const override;
+
+  /// @brief Request the value of the magnetic field at a given position.
+  ///
+  /// @param [in] position Global 3D position for the lookup.
+  /// @param [in, out] cache Cache object.
+  /// @returns A successful value containing a field vector if the given
+  /// location is contained inside any of the regions, or a failure value
+  /// otherwise.
+  Result<Vector3> getField(const Vector3& position,
+                           MagneticFieldProvider::Cache& cache) const override;
+
+  /// @brief Get the field gradient at a given position.
+  ///
+  /// @warning This is not currently implemented.
+  Result<Vector3> getFieldGradient(
+      const Vector3& position, ActsMatrix<3, 3>& /*unused*/,
+      MagneticFieldProvider::Cache& cache) const override;
+};
+}  // namespace Acts
diff --git a/Core/include/Acts/Utilities/AlgebraHelpers.hpp b/Core/include/Acts/Utilities/AlgebraHelpers.hpp
index dd3937ac341..aa522a7dbc2 100644
--- a/Core/include/Acts/Utilities/AlgebraHelpers.hpp
+++ b/Core/include/Acts/Utilities/AlgebraHelpers.hpp
@@ -178,7 +178,7 @@ inline ActsMatrix<A::RowsAtCompileTime, B::ColsAtCompileTime> blockedMult(
 /// Calculate the inverse of an Eigen matrix after checking if it can be
 /// numerically inverted. This allows to catch potential FPEs before they occur.
 /// For matrices up to 4x4, the inverse is computed directly. For larger
-/// matrices, the FullPivLU is used.
+/// matrices, and dynamic matrices the FullPivLU is used.
 ///
 /// @tparam Derived Eigen derived concrete type
 /// @tparam Result Eigen result type defaulted to input type
diff --git a/Core/src/MagneticField/CMakeLists.txt b/Core/src/MagneticField/CMakeLists.txt
index 4b9ebb0be0f..3b164d39505 100644
--- a/Core/src/MagneticField/CMakeLists.txt
+++ b/Core/src/MagneticField/CMakeLists.txt
@@ -1,4 +1,8 @@
 target_sources(
     ActsCore
-    PRIVATE BFieldMapUtils.cpp SolenoidBField.cpp MagneticFieldError.cpp
+    PRIVATE
+        BFieldMapUtils.cpp
+        SolenoidBField.cpp
+        MagneticFieldError.cpp
+        MultiRangeBField.cpp
 )
diff --git a/Core/src/MagneticField/MultiRangeBField.cpp b/Core/src/MagneticField/MultiRangeBField.cpp
new file mode 100644
index 00000000000..8899b50d802
--- /dev/null
+++ b/Core/src/MagneticField/MultiRangeBField.cpp
@@ -0,0 +1,78 @@
+// 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/MagneticField/MultiRangeBField.hpp"
+
+namespace Acts {
+
+MultiRangeBField::Cache::Cache(const MagneticFieldContext& /*unused*/) {}
+
+MultiRangeBField::MultiRangeBField(const std::vector<BFieldRange>& ranges)
+    : fieldRanges(ranges) {}
+
+MultiRangeBField::MultiRangeBField(
+    std::vector<MultiRangeBField::BFieldRange>&& ranges)
+    : fieldRanges(std::move(ranges)) {}
+
+MagneticFieldProvider::Cache MultiRangeBField::makeCache(
+    const MagneticFieldContext& mctx) const {
+  return MagneticFieldProvider::Cache(std::in_place_type<Cache>, mctx);
+}
+
+Result<Vector3> MultiRangeBField::getField(
+    const Vector3& position, MagneticFieldProvider::Cache& cache) const {
+  // Because we assume that the regions are in increasing order of
+  // precedence, we can iterate over the array, remembering the _last_
+  // region that contained the given point. At the end of the loop, this
+  // region will then be the one we query for its field vector.
+  std::optional<std::size_t> foundRange = {};
+
+  // The cache object for this type contains an optional integer; if the
+  // integer is set, it indicates the index of the region in which the last
+  // access succeeded. This acts as a cache because if the current access is
+  // still within that region, it is guaranteed that none of the regions
+  // with lower priority -- i.e. that are earlier in the region vector --
+  // can be relevant to the current access. Thus, we request the cache index
+  // if it exists and perform a membership check on it; if that succeeds, we
+  // remember the corresponding region as a candidate.
+  if (Cache& lCache = cache.as<Cache>();
+      lCache.index.has_value() &&
+      std::get<0>(fieldRanges[*lCache.index])
+          .contains({position[0], position[1], position[2]})) {
+    foundRange = *lCache.index;
+  }
+
+  // Now, we iterate over the ranges. If we already have a range candidate,
+  // we start iteration just after the existing candidate; otherwise we
+  // start from the beginning.
+  for (std::size_t i = (foundRange.has_value() ? (*foundRange) + 1 : 0);
+       i < fieldRanges.size(); ++i) {
+    if (std::get<0>(fieldRanges[i])
+            .contains({position[0], position[1], position[2]})) {
+      foundRange = i;
+    }
+  }
+
+  // Update the cache using the result of this access.
+  cache.as<Cache>().index = foundRange;
+
+  // If we found a valid range, return the corresponding vector; otherwise
+  // return an out-of-bounds error.
+  if (foundRange.has_value()) {
+    return Result<Vector3>::success(std::get<1>(fieldRanges[*foundRange]));
+  } else {
+    return Result<Vector3>::failure(MagneticFieldError::OutOfBounds);
+  }
+}
+
+Result<Vector3> MultiRangeBField::getFieldGradient(
+    const Vector3& position, ActsMatrix<3, 3>& /*unused*/,
+    MagneticFieldProvider::Cache& cache) const {
+  return getField(position, cache);
+}
+}  // namespace Acts
diff --git a/Examples/Io/Obj/CMakeLists.txt b/Examples/Io/Obj/CMakeLists.txt
index 9ca888f726a..bce7bf19e6b 100644
--- a/Examples/Io/Obj/CMakeLists.txt
+++ b/Examples/Io/Obj/CMakeLists.txt
@@ -2,6 +2,7 @@ add_library(
     ActsExamplesIoObj
     SHARED
     src/ObjTrackingGeometryWriter.cpp
+    src/ObjSimHitWriter.cpp
     src/ObjPropagationStepsWriter.cpp
 )
 
diff --git a/Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjPropagationStepsWriter.hpp b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjPropagationStepsWriter.hpp
similarity index 100%
rename from Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjPropagationStepsWriter.hpp
rename to Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjPropagationStepsWriter.hpp
diff --git a/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjSimHitWriter.hpp b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjSimHitWriter.hpp
new file mode 100644
index 00000000000..8ccd92645b9
--- /dev/null
+++ b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjSimHitWriter.hpp
@@ -0,0 +1,84 @@
+// 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/Definitions/Units.hpp"
+#include "Acts/Utilities/Logger.hpp"
+#include "ActsExamples/EventData/SimHit.hpp"
+#include "ActsExamples/Framework/ProcessCode.hpp"
+#include "ActsExamples/Framework/WriterT.hpp"
+
+#include <cstdint>
+#include <mutex>
+#include <string>
+
+namespace ActsExamples {
+struct AlgorithmContext;
+
+/// Write out a simhit collection before detector digitization as wavefront obj
+/// file(s per event).
+///
+/// This writes two files per event, one for the hits one for the trajectory.
+/// The latter can be smoothed using a spline interpolation.
+///
+///     event000000001-<stem>.obj
+///     event000000001-<stem>_trajectory.obj
+///     event000000002-<stem>.obj
+///     event000000002-<stem>_trajectory.obj
+///
+class ObjSimHitWriter : public WriterT<SimHitContainer> {
+ public:
+  struct Config {
+    /// Input sim hit collection to write.
+    std::string inputSimHits;
+    /// Where to place output files
+    std::string outputDir;
+    /// Output filename stem.
+    std::string outputStem = "simhits";
+    /// Number of decimal digits for floating point precision in output.
+    std::size_t outputPrecision = std::numeric_limits<float>::max_digits10;
+    /// Draw line connections between hits
+    bool drawConnections = true;
+    /// Momentum threshold for hits
+    double momentumThreshold = 0.05 * Acts::UnitConstants::GeV;
+    /// Momentum threshold for trajectories
+    double momentumThresholdTraj = 0.05 * Acts::UnitConstants::GeV;
+    /// Number of points to interpolate between hits
+    std::size_t nInterpolatedPoints = 10;
+  };
+
+  /// Construct the particle writer.
+  ///
+  /// @param config is the configuration object
+  /// @param level is the logging level
+  ObjSimHitWriter(const Config& config, Acts::Logging::Level level);
+
+  /// Ensure underlying file is closed.
+  ~ObjSimHitWriter() override = default;
+
+  /// End-of-run hook
+  ProcessCode finalize() override;
+
+  /// Get readonly access to the config parameters
+  const Config& config() const { return m_cfg; }
+
+ protected:
+  /// Type-specific write implementation.
+  ///
+  /// @param[in] ctx is the algorithm context
+  /// @param[in] hits are the hits to be written
+  ProcessCode writeT(const AlgorithmContext& ctx,
+                     const SimHitContainer& hits) override;
+
+ private:
+  Config m_cfg;
+  std::mutex m_writeMutex;
+};
+
+}  // namespace ActsExamples
diff --git a/Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjTrackingGeometryWriter.hpp b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjTrackingGeometryWriter.hpp
similarity index 100%
rename from Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjTrackingGeometryWriter.hpp
rename to Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjTrackingGeometryWriter.hpp
diff --git a/Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjWriterOptions.hpp b/Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjWriterOptions.hpp
similarity index 100%
rename from Examples/Io/Obj/include/ActsExamples/Plugins/Obj/ObjWriterOptions.hpp
rename to Examples/Io/Obj/include/ActsExamples/Io/Obj/ObjWriterOptions.hpp
diff --git a/Examples/Io/Obj/src/ObjPropagationStepsWriter.cpp b/Examples/Io/Obj/src/ObjPropagationStepsWriter.cpp
index 6d0388984cc..0f196ce162c 100644
--- a/Examples/Io/Obj/src/ObjPropagationStepsWriter.cpp
+++ b/Examples/Io/Obj/src/ObjPropagationStepsWriter.cpp
@@ -6,7 +6,7 @@
 // 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 "ActsExamples/Plugins/Obj/ObjPropagationStepsWriter.hpp"
+#include "ActsExamples/Io/Obj/ObjPropagationStepsWriter.hpp"
 
 namespace ActsExamples {
 
diff --git a/Examples/Io/Obj/src/ObjSimHitWriter.cpp b/Examples/Io/Obj/src/ObjSimHitWriter.cpp
new file mode 100644
index 00000000000..564ff9c5409
--- /dev/null
+++ b/Examples/Io/Obj/src/ObjSimHitWriter.cpp
@@ -0,0 +1,190 @@
+// 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 "ActsExamples/Io/Obj/ObjSimHitWriter.hpp"
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/Common.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "ActsExamples/EventData/SimHit.hpp"
+#include "ActsExamples/Framework/AlgorithmContext.hpp"
+#include "ActsExamples/Utilities/Paths.hpp"
+#include "ActsFatras/EventData/Barcode.hpp"
+#include "ActsFatras/EventData/Hit.hpp"
+
+#include <fstream>
+#include <stdexcept>
+#include <unordered_map>
+#include <vector>
+
+#include <unsupported/Eigen/Splines>
+
+namespace {
+
+/// @brief Helper function to interpolate points
+///
+/// @tparam input_vector_type
+/// @param inputs input vector points
+/// @param nPoints number of interpolation points
+///
+/// @return std::vector<Acts::Vector3> interpolated points
+template <typename input_vector_type>
+std::vector<Acts::Vector3> interpolatedPoints(
+    const std::vector<input_vector_type>& inputs, std::size_t nPoints) {
+  std::vector<Acts::Vector3> output;
+
+  if (nPoints < 2) {
+    // No interpolation done return simply the output vector
+    for (const auto& input : inputs) {
+      output.push_back(input.template head<3>());
+    }
+
+  } else {
+    Eigen::MatrixXd points(3, inputs.size());
+    for (std::size_t i = 0; i < inputs.size(); ++i) {
+      points.col(i) = inputs[i].template head<3>().transpose();
+    }
+    Eigen::Spline<double, 3> spline3D =
+        Eigen::SplineFitting<Eigen::Spline<double, 3>>::Interpolate(points, 2);
+
+    double step = 1. / (nPoints - 1);
+    for (std::size_t i = 0; i < nPoints; ++i) {
+      double t = i * step;
+      output.push_back(spline3D(t));
+    }
+  }
+  return output;
+}
+
+}  // namespace
+
+ActsExamples::ObjSimHitWriter::ObjSimHitWriter(
+    const ActsExamples::ObjSimHitWriter::Config& config,
+    Acts::Logging::Level level)
+    : WriterT(config.inputSimHits, "ObjSimHitWriter", level), m_cfg(config) {
+  // inputSimHits is already checked by base constructor
+  if (m_cfg.outputStem.empty()) {
+    throw std::invalid_argument("Missing output filename stem");
+  }
+}
+
+ActsExamples::ProcessCode ActsExamples::ObjSimHitWriter::writeT(
+    const AlgorithmContext& ctx, const ActsExamples::SimHitContainer& simHits) {
+  // ensure exclusive access to tree/file while writing
+  std::scoped_lock lock(m_writeMutex);
+
+  // open per-event file for all simhit components
+  std::string pathSimHit = perEventFilepath(
+      m_cfg.outputDir, m_cfg.outputStem + ".obj", ctx.eventNumber);
+  std::string pathSimTrajectory = perEventFilepath(
+      m_cfg.outputDir, m_cfg.outputStem + "_trajectory.obj", ctx.eventNumber);
+
+  std::ofstream osHits(pathSimHit, std::ofstream::out | std::ofstream::trunc);
+  std::ofstream osTrajectory(pathSimTrajectory,
+                             std::ofstream::out | std::ofstream::trunc);
+
+  if (!osHits || !osTrajectory) {
+    throw std::ios_base::failure("Could not open '" + pathSimHit + "' or '" +
+                                 pathSimTrajectory + "' to write");
+  }
+
+  // Only hit plotting
+  if (!m_cfg.drawConnections) {
+    // Write data from internal immplementation
+    for (const auto& simHit : simHits) {
+      // local simhit information in global coord.
+      const Acts::Vector4& globalPos4 = simHit.fourPosition();
+
+      osHits << "v " << globalPos4[Acts::ePos0] / Acts::UnitConstants::mm << " "
+             << globalPos4[Acts::ePos1] / Acts::UnitConstants::mm << " "
+             << globalPos4[Acts::ePos2] / Acts::UnitConstants::mm << std::endl;
+
+    }  // end simHit loop
+  } else {
+    // We need to associate first
+    std::unordered_map<std::size_t, std::vector<Acts::Vector4>> particleHits;
+    // Pre-loop over hits ... write those below threshold
+    for (const auto& simHit : simHits) {
+      double momentum = simHit.momentum4Before().head<3>().norm();
+      if (momentum < m_cfg.momentumThreshold) {
+        ACTS_VERBOSE("Skipping: Hit below threshold: " << momentum);
+        continue;
+      } else if (momentum < m_cfg.momentumThresholdTraj) {
+        ACTS_VERBOSE(
+            "Skipping (trajectory): Hit below threshold: " << momentum);
+        osHits << "v "
+               << simHit.fourPosition()[Acts::ePos0] / Acts::UnitConstants::mm
+               << " "
+               << simHit.fourPosition()[Acts::ePos1] / Acts::UnitConstants::mm
+               << " "
+               << simHit.fourPosition()[Acts::ePos2] / Acts::UnitConstants::mm
+               << std::endl;
+        continue;
+      }
+      ACTS_VERBOSE("Accepting: Hit above threshold: " << momentum);
+
+      if (particleHits.find(simHit.particleId().value()) ==
+          particleHits.end()) {
+        particleHits[simHit.particleId().value()] = {};
+      }
+      particleHits[simHit.particleId().value()].push_back(
+          simHit.fourPosition());
+    }
+    // Draw loop
+    std::size_t lOffset = 1;
+    for (auto& [pId, pHits] : particleHits) {
+      // Draw the particle hits
+      std::sort(pHits.begin(), pHits.end(),
+                [](const Acts::Vector4& a, const Acts::Vector4& b) {
+                  return a[Acts::eTime] < b[Acts::eTime];
+                });
+
+      osHits << "o particle_" << pId << std::endl;
+      for (const auto& hit : pHits) {
+        osHits << "v " << hit[Acts::ePos0] / Acts::UnitConstants::mm << " "
+               << hit[Acts::ePos1] / Acts::UnitConstants::mm << " "
+               << hit[Acts::ePos2] / Acts::UnitConstants::mm << std::endl;
+      }
+      osHits << '\n';
+
+      // Interpolate the points
+      std::vector<Acts::Vector3> trajectory;
+      if (pHits.size() < 3) {
+        for (const auto& hit : pHits) {
+          trajectory.push_back(hit.template head<3>());
+        }
+      } else {
+        trajectory =
+            interpolatedPoints(pHits, pHits.size() * m_cfg.nInterpolatedPoints);
+      }
+
+      osTrajectory << "o particle_trajectory_" << pId << std::endl;
+      for (const auto& hit : trajectory) {
+        osTrajectory << "v " << hit[Acts::ePos0] / Acts::UnitConstants::mm
+                     << " " << hit[Acts::ePos1] / Acts::UnitConstants::mm << " "
+                     << hit[Acts::ePos2] / Acts::UnitConstants::mm << std::endl;
+      }
+      // Draw the line
+      for (std::size_t iv = lOffset + 1; iv < lOffset + trajectory.size();
+           ++iv) {
+        osTrajectory << "l " << iv - 1 << " " << iv << '\n';
+      }
+      osTrajectory << '\n';
+      // Increase the offset count
+      lOffset += trajectory.size();
+    }
+  }
+  osHits.close();
+  osTrajectory.close();
+
+  return ActsExamples::ProcessCode::SUCCESS;
+}
+
+ActsExamples::ProcessCode ActsExamples::ObjSimHitWriter::finalize() {
+  return ActsExamples::ProcessCode::SUCCESS;
+}
diff --git a/Examples/Io/Obj/src/ObjTrackingGeometryWriter.cpp b/Examples/Io/Obj/src/ObjTrackingGeometryWriter.cpp
index 4937757c85f..d2dbe52e745 100644
--- a/Examples/Io/Obj/src/ObjTrackingGeometryWriter.cpp
+++ b/Examples/Io/Obj/src/ObjTrackingGeometryWriter.cpp
@@ -6,7 +6,7 @@
 // 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 "ActsExamples/Plugins/Obj/ObjTrackingGeometryWriter.hpp"
+#include "ActsExamples/Io/Obj/ObjTrackingGeometryWriter.hpp"
 
 #include "Acts/Utilities/Logger.hpp"
 #include "ActsExamples/Framework/AlgorithmContext.hpp"
diff --git a/Examples/Python/python/acts/examples/simulation.py b/Examples/Python/python/acts/examples/simulation.py
index b60b05466a0..d076a79b30b 100644
--- a/Examples/Python/python/acts/examples/simulation.py
+++ b/Examples/Python/python/acts/examples/simulation.py
@@ -420,6 +420,7 @@ def addFatras(
     outputSimHits: str = "simhits",
     outputDirCsv: Optional[Union[Path, str]] = None,
     outputDirRoot: Optional[Union[Path, str]] = None,
+    outputDirObj: Optional[Union[Path, str]] = None,
     logLevel: Optional[acts.logging.Level] = None,
 ) -> None:
     """This function steers the detector simulation using Fatras
@@ -444,6 +445,8 @@ def addFatras(
         the output folder for the Csv output, None triggers no output
     outputDirRoot : Path|str, path, None
         the output folder for the Root output, None triggers no output
+    outputDirObj : Path|str, path, None
+        the output folder for the Obj output, None triggers no output
     """
 
     customLogLevel = acts.examples.defaultLogging(s, logLevel)
@@ -507,6 +510,7 @@ def addFatras(
         particlesPostSelected,
         outputDirCsv,
         outputDirRoot,
+        outputDirObj,
         logLevel,
     )
 
@@ -519,6 +523,7 @@ def addSimWriters(
     particlesSimulated: str = "particles_simulated",
     outputDirCsv: Optional[Union[Path, str]] = None,
     outputDirRoot: Optional[Union[Path, str]] = None,
+    outputDirObj: Optional[Union[Path, str]] = None,
     logLevel: Optional[acts.logging.Level] = None,
 ) -> None:
     customLogLevel = acts.examples.defaultLogging(s, logLevel)
@@ -563,6 +568,19 @@ def addSimWriters(
             )
         )
 
+    if outputDirObj is not None:
+        outputDirObj = Path(outputDirObj)
+        if not outputDirObj.exists():
+            outputDirObj.mkdir()
+        s.addWriter(
+            acts.examples.ObjSimHitWriter(
+                level=customLogLevel(),
+                inputSimHits=simHits,
+                outputDir=str(outputDirObj),
+                outputStem="hits",
+            )
+        )
+
 
 def getG4DetectorConstructionFactory(
     detector: Any,
@@ -620,6 +638,7 @@ def addGeant4(
     keepParticlesWithoutHits=True,
     outputDirCsv: Optional[Union[Path, str]] = None,
     outputDirRoot: Optional[Union[Path, str]] = None,
+    outputDirObj: Optional[Union[Path, str]] = None,
     logLevel: Optional[acts.logging.Level] = None,
     killVolume: Optional[acts.Volume] = None,
     killAfterTime: float = float("inf"),
@@ -647,6 +666,8 @@ def addGeant4(
         the output folder for the Csv output, None triggers no output
     outputDirRoot : Path|str, path, None
         the output folder for the Root output, None triggers no output
+    outputDirObj : Path|str, path, None
+        the output folder for the Obj output, None triggers no output
     killVolume: acts.Volume, None
         if given, particles are killed when going outside this volume.
     killAfterTime: float
@@ -737,7 +758,8 @@ def addGeant4(
         particlesPostSelected,
         outputDirCsv,
         outputDirRoot,
-        logLevel,
+        outputDirObj,
+        logLevel=logLevel,
     )
 
     return s
diff --git a/Examples/Python/src/Geometry.cpp b/Examples/Python/src/Geometry.cpp
index d69b35a86ba..4c15e0dcf94 100644
--- a/Examples/Python/src/Geometry.cpp
+++ b/Examples/Python/src/Geometry.cpp
@@ -492,6 +492,21 @@ void addExperimentalGeometry(Context& ctx) {
         .def(py::init<std::shared_ptr<KdtSurfacesDim2Bin100>, const Extent&>());
   }
 
+  {
+    using RangeXDDim3 = Acts::RangeXD<3u, double>;
+
+    py::class_<RangeXDDim3>(m, "RangeXDDim3")
+        .def(py::init([](const std::array<double, 2u>& range0,
+                         const std::array<double, 2u>& range1,
+                         const std::array<double, 2u>& range2) {
+          RangeXDDim3 range;
+          range[0].shrink(range0[0], range0[1]);
+          range[1].shrink(range1[0], range1[1]);
+          range[2].shrink(range2[0], range2[1]);
+          return range;
+        }));
+  }
+
   {
     // The external volume structure builder
     py::class_<Acts::Experimental::IExternalStructureBuilder,
diff --git a/Examples/Python/src/MagneticField.cpp b/Examples/Python/src/MagneticField.cpp
index 6e0c1728941..15463386402 100644
--- a/Examples/Python/src/MagneticField.cpp
+++ b/Examples/Python/src/MagneticField.cpp
@@ -12,6 +12,7 @@
 #include "Acts/MagneticField/BFieldMapUtils.hpp"
 #include "Acts/MagneticField/ConstantBField.hpp"
 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
+#include "Acts/MagneticField/MultiRangeBField.hpp"
 #include "Acts/MagneticField/NullBField.hpp"
 #include "Acts/MagneticField/SolenoidBField.hpp"
 #include "Acts/Plugins/Python/Utilities.hpp"
@@ -86,6 +87,11 @@ void addMagneticField(Context& ctx) {
              std::shared_ptr<Acts::NullBField>>(m, "NullBField")
       .def(py::init<>());
 
+  py::class_<Acts::MultiRangeBField, Acts::MagneticFieldProvider,
+             std::shared_ptr<Acts::MultiRangeBField>>(m, "MultiRangeBField")
+      .def(py::init<
+           std::vector<std::pair<Acts::RangeXD<3, double>, Acts::Vector3>>>());
+
   {
     using Config = Acts::SolenoidBField::Config;
 
diff --git a/Examples/Python/src/Output.cpp b/Examples/Python/src/Output.cpp
index 4281f2db26e..251ca2b6357 100644
--- a/Examples/Python/src/Output.cpp
+++ b/Examples/Python/src/Output.cpp
@@ -25,6 +25,9 @@
 #include "ActsExamples/Io/Csv/CsvTrackParameterWriter.hpp"
 #include "ActsExamples/Io/Csv/CsvTrackWriter.hpp"
 #include "ActsExamples/Io/Csv/CsvTrackingGeometryWriter.hpp"
+#include "ActsExamples/Io/Obj/ObjPropagationStepsWriter.hpp"
+#include "ActsExamples/Io/Obj/ObjSimHitWriter.hpp"
+#include "ActsExamples/Io/Obj/ObjTrackingGeometryWriter.hpp"
 #include "ActsExamples/Io/Root/RootBFieldWriter.hpp"
 #include "ActsExamples/Io/Root/RootMaterialTrackWriter.hpp"
 #include "ActsExamples/Io/Root/RootMaterialWriter.hpp"
@@ -46,8 +49,6 @@
 #include "ActsExamples/Io/Root/TrackFitterPerformanceWriter.hpp"
 #include "ActsExamples/Io/Root/VertexNTupleWriter.hpp"
 #include "ActsExamples/MaterialMapping/IMaterialWriter.hpp"
-#include "ActsExamples/Plugins/Obj/ObjPropagationStepsWriter.hpp"
-#include "ActsExamples/Plugins/Obj/ObjTrackingGeometryWriter.hpp"
 #include "ActsExamples/TrackFinding/ITrackParamsLookupReader.hpp"
 #include "ActsExamples/TrackFinding/ITrackParamsLookupWriter.hpp"
 
@@ -111,6 +112,10 @@ void addOutput(Context& ctx) {
                              "ObjPropagationStepsWriter", collection, outputDir,
                              outputScalor, outputPrecision);
 
+  ACTS_PYTHON_DECLARE_WRITER(ActsExamples::ObjSimHitWriter, mex,
+                             "ObjSimHitWriter", inputSimHits, outputDir,
+                             outputStem, outputPrecision, drawConnections);
+
   {
     auto c = py::class_<ViewConfig>(m, "ViewConfig").def(py::init<>());
 
diff --git a/Examples/Python/tests/test_magnetic_field.py b/Examples/Python/tests/test_magnetic_field.py
index a51ebd55ab6..d1336e3ca8a 100644
--- a/Examples/Python/tests/test_magnetic_field.py
+++ b/Examples/Python/tests/test_magnetic_field.py
@@ -72,3 +72,41 @@ def test_solenoid(conf_const):
     )
 
     assert isinstance(field, acts.examples.InterpolatedMagneticField2)
+
+
+def test_multiregion_bfield():
+    with pytest.raises(TypeError):
+        acts.MultiRangeBField()
+
+    rs = [
+        (acts.RangeXDDim3((0, 3), (0, 3), (0, 3)), acts.Vector3(0.0, 0.0, 2.0)),
+        (acts.RangeXDDim3((1, 2), (1, 2), (1, 10)), acts.Vector3(2.0, 0.0, 0.0)),
+    ]
+    f = acts.MultiRangeBField(rs)
+    assert f
+
+    ctx = acts.MagneticFieldContext()
+    assert ctx
+
+    fc = f.makeCache(ctx)
+    assert fc
+
+    rv = f.getField(acts.Vector3(0.5, 0.5, 0.5), fc)
+    assert rv[0] == pytest.approx(0.0)
+    assert rv[1] == pytest.approx(0.0)
+    assert rv[2] == pytest.approx(2.0)
+
+    rv = f.getField(acts.Vector3(1.5, 1.5, 5.0), fc)
+    assert rv[0] == pytest.approx(2.0)
+    assert rv[1] == pytest.approx(0.0)
+    assert rv[2] == pytest.approx(0.0)
+
+    rv = f.getField(acts.Vector3(2.5, 2.5, 2.5), fc)
+    assert rv[0] == pytest.approx(0.0)
+    assert rv[1] == pytest.approx(0.0)
+    assert rv[2] == pytest.approx(2.0)
+
+    rv = f.getField(acts.Vector3(1.5, 1.5, 1.5), fc)
+    assert rv[0] == pytest.approx(2.0)
+    assert rv[1] == pytest.approx(0.0)
+    assert rv[2] == pytest.approx(0.0)
diff --git a/Examples/Scripts/Python/full_chain_odd.py b/Examples/Scripts/Python/full_chain_odd.py
index 8b86402b820..7ddaa02054f 100755
--- a/Examples/Scripts/Python/full_chain_odd.py
+++ b/Examples/Scripts/Python/full_chain_odd.py
@@ -136,6 +136,12 @@
     default=True,
     action=argparse.BooleanOptionalAction,
 )
+parser.add_argument(
+    "--output-obj",
+    help="Switch obj output on/off",
+    default=True,
+    action=argparse.BooleanOptionalAction,
+)
 
 args = parser.parse_args()
 
@@ -277,6 +283,7 @@
             ),
             outputDirRoot=outputDir if args.output_root else None,
             outputDirCsv=outputDir if args.output_csv else None,
+            outputDirObj=outputDir if args.output_obj else None,
             rnd=rnd,
             killVolume=trackingGeometry.highestTrackingVolume,
             killAfterTime=25 * u.ns,
@@ -305,6 +312,7 @@
             enableInteractions=True,
             outputDirRoot=outputDir if args.output_root else None,
             outputDirCsv=outputDir if args.output_csv else None,
+            outputDirObj=outputDir if args.output_obj else None,
             rnd=rnd,
         )
 
diff --git a/Examples/Scripts/Python/geant4.py b/Examples/Scripts/Python/geant4.py
index 42cb5db95bb..65e154c2704 100755
--- a/Examples/Scripts/Python/geant4.py
+++ b/Examples/Scripts/Python/geant4.py
@@ -36,6 +36,7 @@ def runGeant4(
         field,
         outputDirCsv=outputDir / "csv",
         outputDirRoot=outputDir,
+        outputDirObj=outputDir / "obj",
         rnd=rnd,
         materialMappings=materialMappings,
         volumeMappings=volumeMappings,
diff --git a/Tests/UnitTests/Core/MagneticField/CMakeLists.txt b/Tests/UnitTests/Core/MagneticField/CMakeLists.txt
index 60a829bf548..4600cc8352e 100644
--- a/Tests/UnitTests/Core/MagneticField/CMakeLists.txt
+++ b/Tests/UnitTests/Core/MagneticField/CMakeLists.txt
@@ -1,4 +1,5 @@
 add_unittest(ConstantBField ConstantBFieldTests.cpp)
 add_unittest(InterpolatedBFieldMap InterpolatedBFieldMapTests.cpp)
 add_unittest(SolenoidBField SolenoidBFieldTests.cpp)
+add_unittest(MultiRangeBField MultiRangeBFieldTests.cpp)
 add_unittest(MagneticFieldProvider MagneticFieldProviderTests.cpp)
diff --git a/Tests/UnitTests/Core/MagneticField/MultiRangeBFieldTests.cpp b/Tests/UnitTests/Core/MagneticField/MultiRangeBFieldTests.cpp
new file mode 100644
index 00000000000..3644db058b4
--- /dev/null
+++ b/Tests/UnitTests/Core/MagneticField/MultiRangeBFieldTests.cpp
@@ -0,0 +1,74 @@
+// 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/floating_point_comparison.hpp>
+#include <boost/test/unit_test.hpp>
+
+#include "Acts/MagneticField/MagneticFieldContext.hpp"
+#include "Acts/MagneticField/MultiRangeBField.hpp"
+#include "Acts/Utilities/Result.hpp"
+
+namespace Acts::Test {
+
+MagneticFieldContext mfContext = MagneticFieldContext();
+
+BOOST_AUTO_TEST_CASE(TestMultiRangeBField) {
+  std::vector<std::pair<RangeXD<3, double>, Vector3>> inputs;
+
+  inputs.emplace_back(RangeXD<3, double>{{0., 0., 0.}, {3., 3., 3.}},
+                      Vector3{0., 0., 2.});
+  inputs.emplace_back(RangeXD<3, double>{{1., 1., 1.}, {2., 2., 10.}},
+                      Vector3{2., 0., 0.});
+
+  const MultiRangeBField bfield(std::move(inputs));
+
+  auto bcache = bfield.makeCache(mfContext);
+
+  // Test a point inside the first volume.
+  {
+    Result<Vector3> r = bfield.getField({0.5, 0.5, 0.5}, bcache);
+    BOOST_CHECK(r.ok());
+    BOOST_CHECK_CLOSE((*r)[0], 0., 0.05);
+    BOOST_CHECK_CLOSE((*r)[1], 0., 0.05);
+    BOOST_CHECK_CLOSE((*r)[2], 2., 0.05);
+  }
+
+  // Test a point inside the second volume, not overlapping the first.
+  {
+    Result<Vector3> r = bfield.getField({1.5, 1.5, 5.}, bcache);
+    BOOST_CHECK(r.ok());
+    BOOST_CHECK_CLOSE((*r)[0], 2., 0.05);
+    BOOST_CHECK_CLOSE((*r)[1], 0., 0.05);
+    BOOST_CHECK_CLOSE((*r)[2], 0., 0.05);
+  }
+
+  // Test a point inside the first volume again.
+  {
+    Result<Vector3> r = bfield.getField({2.5, 2.5, 2.5}, bcache);
+    BOOST_CHECK(r.ok());
+    BOOST_CHECK_CLOSE((*r)[0], 0., 0.05);
+    BOOST_CHECK_CLOSE((*r)[1], 0., 0.05);
+    BOOST_CHECK_CLOSE((*r)[2], 2., 0.05);
+  }
+
+  // Test a point inside the second volume in the overlap region.
+  {
+    Result<Vector3> r = bfield.getField({1.5, 1.5, 1.5}, bcache);
+    BOOST_CHECK(r.ok());
+    BOOST_CHECK_CLOSE((*r)[0], 2., 0.05);
+    BOOST_CHECK_CLOSE((*r)[1], 0., 0.05);
+    BOOST_CHECK_CLOSE((*r)[2], 0., 0.05);
+  }
+
+  // Test a point outside all volumes.
+  {
+    Result<Vector3> r = bfield.getField({-1., -1., -1}, bcache);
+    BOOST_CHECK(!r.ok());
+  }
+}
+}  // namespace Acts::Test
diff --git a/Tests/UnitTests/Core/Utilities/AlgebraHelpersTests.cpp b/Tests/UnitTests/Core/Utilities/AlgebraHelpersTests.cpp
index 86165ddb2ba..87e47d3e3a9 100644
--- a/Tests/UnitTests/Core/Utilities/AlgebraHelpersTests.cpp
+++ b/Tests/UnitTests/Core/Utilities/AlgebraHelpersTests.cpp
@@ -51,6 +51,15 @@ BOOST_AUTO_TEST_CASE(safeInverseLargeMatrix) {
   BOOST_CHECK_EQUAL(*identityInv, identity);
 }
 
+BOOST_AUTO_TEST_CASE(safeInverseDynamicMatrix) {
+  Eigen::MatrixXd identity{Eigen::MatrixXd::Identity(2, 2)};
+
+  auto identityInv = Acts::safeInverse(identity);
+
+  BOOST_CHECK(identityInv);
+  BOOST_CHECK_EQUAL(*identityInv, identity);
+}
+
 BOOST_AUTO_TEST_CASE(SafeInverseBadSmallMatrix) {
   Eigen::Matrix<double, 2, 2> m;
   m << 1, 1, 2, 2;
@@ -111,13 +120,6 @@ BOOST_AUTO_TEST_CASE(SafeInverseFPELargeMatrix) {
 //   auto mInv = Acts::safeInverse(m);
 // }
 
-/// This test should not compile
-// BOOST_AUTO_TEST_CASE(SafeInverseDynamicMatrix) {
-//   Eigen::MatrixXd m{Eigen::MatrixXd::Identity(2, 2)};
-//
-//   auto mInv = Acts::safeInverse(m);
-// }
-
 BOOST_AUTO_TEST_SUITE_END()
 
 BOOST_AUTO_TEST_SUITE(SafeExp)
diff --git a/docs/core/magnetic_field.md b/docs/core/magnetic_field.md
index 1924234bde0..488d5890203 100644
--- a/docs/core/magnetic_field.md
+++ b/docs/core/magnetic_field.md
@@ -235,6 +235,28 @@ analytical implementation and is much faster to lookup:
 :::{doxygenfunction} Acts::solenoidFieldMap
 :::
 
+### Multi-range constant field
+
+The multi-range constant field allows modelling cases where a magnetic field
+can be described as multiple (potentially overlapping) regions, each of which
+has its own constant magnetic field. This provides more flexibility than the
+{class}`Acts::ConstantBField` while providing higher performance than
+{class}`Acts::InterpolatedBFieldMap`.
+
+This magnetic field provider is configured using a list of pairs, where each
+pair defines a region in three-dimensional space as well as a field vector.
+Magnetic field lookup then proceeds by finding the _last_ region in the
+user-provided list that contains the requested coordinate and returning the
+corresponding field vector.
+
+The implementation uses a simple caching mechanism to store the last matched
+region, providing improved performance for consecutive lookups within the same
+region. This is thread-safe when each thread uses its own cache instance. The
+field configuration itself is immutable after construction.
+
+:::{doxygenclass} Acts::MultiRangeBField
+:::
+
 ## Full provider interface
 
 :::{doxygenclass} Acts::MagneticFieldProvider
diff --git a/docs/tracking.md b/docs/tracking.md
index 54e582d1e0f..442a3793adc 100644
--- a/docs/tracking.md
+++ b/docs/tracking.md
@@ -97,16 +97,17 @@ position $l$ and the momentum vector $\vec p$ are shown.
 
 Aside from the nominal quantities captured in $\vec x$, the related
 uncertainties and correlations need to be taken into account as well. They
-can be expressed as a $5\times 5$ covariance matrix like
+can be expressed as a $6\times 6$ covariance matrix like
 
 \begin{equation*}
   C =
   \begin{bmatrix}
-   \sigma^2(l_0)& \text{cov}(l_0,l_1) & \text{cov}(l_0, \phi) & \text{cov}(l_0, \theta) & \text{cov}(l_0, q/p) \\
-   . & \sigma^2(l_1) & \text{cov}(l_1, \phi) & \text{cov}(l_1, \theta) & \text{cov}(l_1, q/p) \\
-   . & . &  \sigma^2(\phi) & \text{cov}(\phi,\theta) & \text{cov}(\phi, q/p) \\
-   . & . & . & \sigma^2(\theta) & \text{cov}(\theta, q/p) \\
-   . & . & . & . & \sigma^2(q/p)
+   \sigma^2(l_0)& \text{cov}(l_0,l_1) & \text{cov}(l_0, \phi) & \text{cov}(l_0, \theta) & \text{cov}(l_0, q/p) & \text{cov}(l_0, t) \\
+   . & \sigma^2(l_1) & \text{cov}(l_1, \phi) & \text{cov}(l_1, \theta) & \text{cov}(l_1, q/p) & \text{cov}(l_1, t) \\
+   . & . &  \sigma^2(\phi) & \text{cov}(\phi,\theta) & \text{cov}(\phi, q/p) & \text{cov}(\phi, t) \\
+   . & . & . & \sigma^2(\theta) & \text{cov}(\theta, q/p) & \text{cov}(\theta, t) \\
+   . & . & . & . & \sigma^2(q/p) & \text{cov}(q/p, t) \\
+   . & . & . & . & . & \sigma^2(t)
   \end{bmatrix}
 \end{equation*}