diff --git a/CI/physmon/reference/performance_amvf_gridseeder_seeded_hist.root b/CI/physmon/reference/performance_amvf_gridseeder_seeded_hist.root
index df4b1db76b0..882a5fc6d30 100644
Binary files a/CI/physmon/reference/performance_amvf_gridseeder_seeded_hist.root and b/CI/physmon/reference/performance_amvf_gridseeder_seeded_hist.root differ
diff --git a/Core/include/Acts/Detector/DetectorVolume.hpp b/Core/include/Acts/Detector/DetectorVolume.hpp
index d16a06ca9cc..c0e6a91fb8c 100644
--- a/Core/include/Acts/Detector/DetectorVolume.hpp
+++ b/Core/include/Acts/Detector/DetectorVolume.hpp
@@ -484,6 +484,9 @@ class DetectorVolume : public std::enable_shared_from_this<DetectorVolume> {
 /// @brief  A detector volume factory which first constructs the detector volume
 /// and then constructs the portals. This ensures that the std::shared_ptr
 /// holding the detector volume is not weak when assigning to the portals.
+///
+/// @note Optional containment check is invoked by setting the number
+/// of segments nSeg to be greater than 0
 class DetectorVolumeFactory {
  public:
   /// Create a detector volume - from factory
@@ -494,11 +497,19 @@ class DetectorVolumeFactory {
       const std::vector<std::shared_ptr<Surface>>& surfaces,
       const std::vector<std::shared_ptr<DetectorVolume>>& volumes,
       DetectorVolumeUpdater detectorVolumeUpdater,
-      SurfaceCandidatesUpdater surfaceCandidateUpdater) {
+      SurfaceCandidatesUpdater surfaceCandidateUpdater, int nSeg = -1) {
     auto dVolume = DetectorVolume::makeShared(
         gctx, name, transform, std::move(bounds), surfaces, volumes,
         std::move(detectorVolumeUpdater), std::move(surfaceCandidateUpdater));
     dVolume->construct(gctx, portalGenerator);
+
+    /// Volume extent is constructed from the portals
+    /// So the surface/subvolume containment
+    /// check has to happen here
+    if (nSeg > 0 && !dVolume->checkContainment(gctx, nSeg)) {
+      throw std::invalid_argument(
+          "DetectorVolume: surfaces or subvolumes are not contained by volume");
+    }
     return dVolume;
   }
 
diff --git a/Core/include/Acts/Geometry/ConeVolumeBounds.hpp b/Core/include/Acts/Geometry/ConeVolumeBounds.hpp
index 93b7283e9ae..c875f27b6b8 100644
--- a/Core/include/Acts/Geometry/ConeVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/ConeVolumeBounds.hpp
@@ -56,9 +56,10 @@ class ConeVolumeBounds : public VolumeBounds {
   /// @param halflengthZ The minimum z value of the inner and outer cones
   /// @param averagePhi The phi orientation of the sector
   /// @param halfPhiSector The opening angle phi sector
-  ConeVolumeBounds(double innerAlpha, double innerOffsetZ, double outerAlpha,
-                   double outerOffsetZ, double halflengthZ, double averagePhi,
-                   double halfPhiSector) noexcept(false);
+  ConeVolumeBounds(ActsScalar innerAlpha, ActsScalar innerOffsetZ,
+                   ActsScalar outerAlpha, ActsScalar outerOffsetZ,
+                   ActsScalar halflengthZ, ActsScalar averagePhi,
+                   ActsScalar halfPhiSector) noexcept(false);
 
   /// Constructor - for general cylidner-cone setups
   ///
@@ -71,14 +72,14 @@ class ConeVolumeBounds : public VolumeBounds {
   ///
   /// @note depending on cylinderR > coneR it is constructing a cone with
   /// cylindrical cutout or a cylinder with conical cutout
-  ConeVolumeBounds(double cylinderR, double alpha, double offsetZ,
-                   double halflengthZ, double averagePhi,
-                   double halfPhiSector) noexcept(false);
+  ConeVolumeBounds(ActsScalar cylinderR, ActsScalar alpha, ActsScalar offsetZ,
+                   ActsScalar halflengthZ, ActsScalar averagePhi,
+                   ActsScalar halfPhiSector) noexcept(false);
 
   /// Constructor - from a fixed size array
   ///
   /// @param values The bound values
-  ConeVolumeBounds(const std::array<double, eSize>& values) noexcept(false)
+  ConeVolumeBounds(const std::array<ActsScalar, eSize>& values) noexcept(false)
       : m_values(values) {
     checkConsistency();
     buildSurfaceBounds();
@@ -93,14 +94,14 @@ class ConeVolumeBounds : public VolumeBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// This method checks if position in the 3D volume
   /// frame is inside the cylinder
   ///
   /// @param pos is the position in volume frame to be checked
   /// @param tol is the absolute tolerance to be applied
-  bool inside(const Vector3& pos, double tol = 0.) const final;
+  bool inside(const Vector3& pos, ActsScalar tol = 0.) const final;
 
   /// Oriented surfaces, i.e. the decomposed boundary surfaces and the
   /// according navigation direction into the volume given the normal
@@ -126,30 +127,30 @@ class ConeVolumeBounds : public VolumeBounds {
 
   /// Access to the bound values
   /// @param bValue the class nested enum for the array access
-  double get(BoundValues bValue) const { return m_values[bValue]; }
+  ActsScalar get(BoundValues bValue) const { return m_values[bValue]; }
 
   // Return the derived innerRmin
-  double innerRmin() const;
+  ActsScalar innerRmin() const;
 
   // Return the derived innerRmin
-  double innerRmax() const;
+  ActsScalar innerRmax() const;
 
   // Return the derived inner tan(alpha)
-  double innerTanAlpha() const;
+  ActsScalar innerTanAlpha() const;
 
   // Return the derived outerRmin
-  double outerRmin() const;
+  ActsScalar outerRmin() const;
 
   // Return the derived outerRmax
-  double outerRmax() const;
+  ActsScalar outerRmax() const;
 
   // Return the derived outer tan(alpha)
-  double outerTanAlpha() const;
+  ActsScalar outerTanAlpha() const;
 
   /// Output Method for std::ostream
   ///
-  /// @param sl is ostream operator to be dumped into
-  std::ostream& toStream(std::ostream& sl) const final;
+  /// @param os is ostream operator to be dumped into
+  std::ostream& toStream(std::ostream& os) const final;
 
  private:
   /// Check the input values for consistency,
@@ -159,14 +160,8 @@ class ConeVolumeBounds : public VolumeBounds {
   /// Create the surface bounds
   void buildSurfaceBounds();
 
-  /// Templated dump methods
-  /// @tparam stream_t The type of the stream for dumping
-  /// @param dt The stream object
-  template <class stream_t>
-  stream_t& dumpT(stream_t& dt) const;
-
   /// The bound values
-  std::array<double, eSize> m_values;
+  std::array<ActsScalar, eSize> m_values;
   std::shared_ptr<CylinderBounds> m_innerCylinderBounds{nullptr};
   std::shared_ptr<ConeBounds> m_innerConeBounds{nullptr};
   std::shared_ptr<ConeBounds> m_outerConeBounds{nullptr};
@@ -176,54 +171,12 @@ class ConeVolumeBounds : public VolumeBounds {
   std::shared_ptr<PlanarBounds> m_sectorBounds{nullptr};
 
   /// Derived values
-  double m_innerRmin = 0.;
-  double m_innerRmax = 0.;
-  double m_innerTanAlpha = 0.;
-  double m_outerRmin = 0.;
-  double m_outerRmax = 0.;
-  double m_outerTanAlpha = 0.;
+  ActsScalar m_innerRmin = 0.;
+  ActsScalar m_innerRmax = 0.;
+  ActsScalar m_innerTanAlpha = 0.;
+  ActsScalar m_outerRmin = 0.;
+  ActsScalar m_outerRmax = 0.;
+  ActsScalar m_outerTanAlpha = 0.;
 };
 
-inline double ConeVolumeBounds::innerRmin() const {
-  return m_innerRmin;
-}
-
-inline double ConeVolumeBounds::innerRmax() const {
-  return m_innerRmax;
-}
-
-inline double ConeVolumeBounds::innerTanAlpha() const {
-  return m_innerTanAlpha;
-}
-
-inline double ConeVolumeBounds::outerRmin() const {
-  return m_outerRmin;
-}
-
-inline double ConeVolumeBounds::outerRmax() const {
-  return m_outerRmax;
-}
-
-inline double ConeVolumeBounds::outerTanAlpha() const {
-  return m_outerTanAlpha;
-}
-
-inline std::vector<double> ConeVolumeBounds::values() const {
-  std::vector<double> valvector;
-  valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
-  return valvector;
-}
-
-template <class stream_t>
-stream_t& ConeVolumeBounds::dumpT(stream_t& dt) const {
-  dt << std::setiosflags(std::ios::fixed);
-  dt << std::setprecision(5);
-  dt << "Acts::ConeVolumeBounds : (innerAlpha, innerOffsetZ, outerAlpha,";
-  dt << "  outerOffsetZ, halflenghZ, averagePhi, halfPhiSector) = ";
-  dt << get(eInnerAlpha) << ", " << get(eInnerOffsetZ) << ", ";
-  dt << get(eOuterAlpha) << ", " << get(eOuterOffsetZ) << ", ";
-  dt << get(eHalfLengthZ) << ", " << get(eAveragePhi) << std::endl;
-  return dt;
-}
-
 }  // namespace Acts
diff --git a/Core/include/Acts/Geometry/CuboidVolumeBounds.hpp b/Core/include/Acts/Geometry/CuboidVolumeBounds.hpp
index 52f8697688b..5605dd7532d 100644
--- a/Core/include/Acts/Geometry/CuboidVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/CuboidVolumeBounds.hpp
@@ -63,12 +63,14 @@ class CuboidVolumeBounds : public VolumeBounds {
   /// @param halex is the half length of the cube in x
   /// @param haley is the half length of the cube in y
   /// @param halez is the half length of the cube in z
-  CuboidVolumeBounds(double halex, double haley, double halez) noexcept(false);
+  CuboidVolumeBounds(ActsScalar halex, ActsScalar haley,
+                     ActsScalar halez) noexcept(false);
 
   /// Constructor - from a fixed size array
   ///
   /// @param values iw the bound values
-  CuboidVolumeBounds(const std::array<double, eSize>& values) noexcept(false)
+  CuboidVolumeBounds(const std::array<ActsScalar, eSize>& values) noexcept(
+      false)
       : m_values(values) {
     checkConsistency();
     buildSurfaceBounds();
@@ -91,14 +93,14 @@ class CuboidVolumeBounds : public VolumeBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// This method checks if position in the 3D volume
   /// frame is inside the cylinder
   ///
   /// @param pos is the position in volume frame to be checked
   /// @param tol is the absolute tolerance to be applied
-  bool inside(const Vector3& pos, double tol = 0.) const override;
+  bool inside(const Vector3& pos, ActsScalar tol = 0.) const override;
 
   /// Oriented surfaces, i.e. the decomposed boundary surfaces and the
   /// according navigation direction into the volume given the normal
@@ -122,31 +124,33 @@ class CuboidVolumeBounds : public VolumeBounds {
                                   const Vector3& envelope = {0, 0, 0},
                                   const Volume* entity = nullptr) const final;
 
-  /// Binning borders in double
+  /// Get the canonical binning values, i.e. the binning values
+  /// for that fully describe the shape's extent
+  ///
+  /// @return vector of canonical binning values
+  std::vector<Acts::BinningValue> canonicalBinning() const override {
+    return {Acts::binX, Acts::binY, Acts::binZ};
+  };
+
+  /// Binning borders in ActsScalar
   ///
   /// @param bValue is the binning schema used
   ///
   /// @return float offset to be used for the binning
-  double binningBorder(BinningValue bValue) const final;
+  ActsScalar binningBorder(BinningValue bValue) const final;
 
   /// Access to the bound values
   /// @param bValue the class nested enum for the array access
-  double get(BoundValues bValue) const { return m_values[bValue]; }
+  ActsScalar get(BoundValues bValue) const { return m_values[bValue]; }
 
   /// Output Method for std::ostream
   ///
-  /// @param sl is ostream operator to be dumped into
-  std::ostream& toStream(std::ostream& sl) const override;
+  /// @param os is ostream operator to be dumped into
+  std::ostream& toStream(std::ostream& os) const override;
 
  private:
-  /// Templated dumpT method
-  /// @tparam stream_t The type for the dump stream
-  /// @param dt The dump stream object
-  template <class stream_t>
-  stream_t& dumpT(stream_t& dt) const;
-
   /// The bound values ordered in a fixed size array
-  std::array<double, eSize> m_values;
+  std::array<ActsScalar, eSize> m_values;
 
   std::shared_ptr<const RectangleBounds> m_xyBounds{nullptr};
   std::shared_ptr<const RectangleBounds> m_yzBounds{nullptr};
@@ -159,34 +163,4 @@ class CuboidVolumeBounds : public VolumeBounds {
   /// will throw a logic_exception if consistency is not given
   void checkConsistency() noexcept(false);
 };
-
-inline bool CuboidVolumeBounds::inside(const Vector3& pos, double tol) const {
-  return (std::abs(pos.x()) <= get(eHalfLengthX) + tol &&
-          std::abs(pos.y()) <= get(eHalfLengthY) + tol &&
-          std::abs(pos.z()) <= get(eHalfLengthZ) + tol);
-}
-
-inline std::vector<double> CuboidVolumeBounds::values() const {
-  std::vector<double> valvector;
-  valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
-  return valvector;
-}
-
-inline void CuboidVolumeBounds::checkConsistency() noexcept(false) {
-  if (get(eHalfLengthX) <= 0 || get(eHalfLengthY) <= 0 ||
-      get(eHalfLengthZ) <= 0.) {
-    throw std::invalid_argument(
-        "CuboidVolumeBounds: invalid input, zero or negative.");
-  }
-}
-
-template <class stream_t>
-stream_t& CuboidVolumeBounds::dumpT(stream_t& dt) const {
-  dt << std::setiosflags(std::ios::fixed);
-  dt << std::setprecision(5);
-  dt << "Acts::CuboidVolumeBounds: (halfLengthX, halfLengthY, halfLengthZ) = ";
-  dt << "(" << get(eHalfLengthX) << ", " << get(eHalfLengthY) << ", "
-     << get(eHalfLengthZ) << ")";
-  return dt;
-}
 }  // namespace Acts
diff --git a/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp b/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp
index 52031d2514f..3033b83595e 100644
--- a/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp
@@ -114,6 +114,14 @@ class CutoutCylinderVolumeBounds : public VolumeBounds {
                                   const Vector3& envelope = {0, 0, 0},
                                   const Volume* entity = nullptr) const final;
 
+  /// Get the canonical binning values, i.e. the binning values
+  /// for that fully describe the shape's extent
+  ///
+  /// @return vector of canonical binning values
+  std::vector<Acts::BinningValue> canonicalBinning() const override {
+    return {Acts::binR, Acts::binPhi, Acts::binZ};
+  };
+
   /// Write information about this instance to an outstream
   ///
   /// @param sl The outstream
diff --git a/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp b/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
index 4f8eb9dddca..6676659ebc0 100644
--- a/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
@@ -94,10 +94,10 @@ class CylinderVolumeBounds : public VolumeBounds {
   /// @param avgphi The average phi value
   /// @param bevelMinZ The bevel angle, in radians, for the negative side
   /// @param bevelMaxZ The bevel angle, in radians, for the positive side
-  CylinderVolumeBounds(double rmin, double rmax, double halfz,
-                       double halfphi = M_PI, double avgphi = 0.,
-                       double bevelMinZ = 0.,
-                       double bevelMaxZ = 0.) noexcept(false)
+  CylinderVolumeBounds(ActsScalar rmin, ActsScalar rmax, ActsScalar halfz,
+                       ActsScalar halfphi = M_PI, ActsScalar avgphi = 0.,
+                       ActsScalar bevelMinZ = 0.,
+                       ActsScalar bevelMaxZ = 0.) noexcept(false)
       : m_values({rmin, rmax, halfz, halfphi, avgphi, bevelMinZ, bevelMaxZ}) {
     checkConsistency();
     buildSurfaceBounds();
@@ -106,7 +106,8 @@ class CylinderVolumeBounds : public VolumeBounds {
   /// Constructor - from a fixed size array
   ///
   /// @param values The bound values
-  CylinderVolumeBounds(const std::array<double, eSize>& values) noexcept(false)
+  CylinderVolumeBounds(const std::array<ActsScalar, eSize>& values) noexcept(
+      false)
       : m_values(values) {
     checkConsistency();
     buildSurfaceBounds();
@@ -117,14 +118,14 @@ class CylinderVolumeBounds : public VolumeBounds {
   /// @param cBounds the cylinder bounds
   /// @param thickness of the extrusion
   CylinderVolumeBounds(const CylinderBounds& cBounds,
-                       double thickness) noexcept(false);
+                       ActsScalar thickness) noexcept(false);
 
   /// Constructor - extruded from radial bounds and thickness
   ///
   /// @param rBounds the Radial bounds
   /// @param thickness
   CylinderVolumeBounds(const RadialBounds& rBounds,
-                       double thickness) noexcept(false);
+                       ActsScalar thickness) noexcept(false);
 
   /// Copy Constructor
   ///
@@ -141,14 +142,14 @@ class CylinderVolumeBounds : public VolumeBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// This method checks if position in the 3D volume
   /// frame is inside the cylinder
   ///
   /// @param pos is a global position to be checked
   /// @param tol is the tolerance for the check
-  bool inside(const Vector3& pos, double tol = 0.) const override;
+  bool inside(const Vector3& pos, ActsScalar tol = 0.) const override;
 
   /// Oriented surfaces, i.e. the decomposed boundary surfaces and the
   /// according navigation direction into the volume given the normal
@@ -172,26 +173,35 @@ class CylinderVolumeBounds : public VolumeBounds {
                                   const Vector3& envelope = {0, 0, 0},
                                   const Volume* entity = nullptr) const final;
 
+  /// Get the canonical binning values, i.e. the binning values
+  /// for that fully describe the shape's extent
+  ///
+  /// @return vector of canonical binning values
+  std::vector<Acts::BinningValue> canonicalBinning() const override {
+    return {Acts::binR, Acts::binPhi, Acts::binZ};
+  };
+
   /// Binning offset - overloaded for some R-binning types
   ///
   /// @param bValue is the type used for the binning
   Vector3 binningOffset(BinningValue bValue) const override;
 
-  /// Binning borders in double
+  /// Binning borders in ActsScalar
   ///
   /// @param bValue is the type used for the binning
-  double binningBorder(BinningValue bValue) const override;
+  ActsScalar binningBorder(BinningValue bValue) const override;
 
   /// Output Method for std::ostream
-  std::ostream& toStream(std::ostream& sl) const override;
+  /// @param os is the output stream
+  std::ostream& toStream(std::ostream& os) const override;
 
   /// Access to the bound values
   /// @param bValue the class nested enum for the array access
-  double get(BoundValues bValue) const { return m_values[bValue]; }
+  ActsScalar get(BoundValues bValue) const { return m_values[bValue]; }
 
  private:
-  /// The internal version of the bounds can be float/double
-  std::array<double, eSize> m_values{};
+  /// The internal version of the bounds can be float/ActsScalar
+  std::array<ActsScalar, eSize> m_values{};
   /// Bounds of the inner CylinderBounds
   std::shared_ptr<const CylinderBounds> m_innerCylinderBounds{nullptr};
   /// Bounds of the inner CylinderBounds
@@ -207,85 +217,6 @@ class CylinderVolumeBounds : public VolumeBounds {
 
   /// Helper method to create the surface bounds
   void buildSurfaceBounds();
-
-  /// Templated dumpT method
-  /// @tparam stream_t The type for the dump stream
-  /// @param dt The dump stream object
-  template <class stream_t>
-  stream_t& dumpT(stream_t& dt) const;
 };
 
-inline bool CylinderVolumeBounds::inside(const Vector3& pos, double tol) const {
-  using VectorHelpers::perp;
-  using VectorHelpers::phi;
-  double ros = perp(pos);
-  bool insidePhi = cos(phi(pos)) >= cos(get(eHalfPhiSector)) - tol;
-  bool insideR = insidePhi
-                     ? ((ros >= get(eMinR) - tol) && (ros <= get(eMaxR) + tol))
-                     : false;
-  bool insideZ =
-      insideR ? (std::abs(pos.z()) <= get(eHalfLengthZ) + tol) : false;
-  return (insideZ && insideR && insidePhi);
-}
-
-inline Vector3 CylinderVolumeBounds::binningOffset(BinningValue bValue)
-    const {  // the medium radius is taken for r-type binning
-  if (bValue == Acts::binR || bValue == Acts::binRPhi) {
-    return Vector3(0.5 * (get(eMinR) + get(eMaxR)), 0., 0.);
-  }
-  return VolumeBounds::binningOffset(bValue);
-}
-
-inline double CylinderVolumeBounds::binningBorder(BinningValue bValue) const {
-  if (bValue == Acts::binR) {
-    return 0.5 * (get(eMaxR) - get(eMinR));
-  }
-  if (bValue == Acts::binZ) {
-    return get(eHalfLengthZ);
-  }
-  return VolumeBounds::binningBorder(bValue);
-}
-
-template <class stream_t>
-stream_t& CylinderVolumeBounds::dumpT(stream_t& dt) const {
-  dt << std::setiosflags(std::ios::fixed);
-  dt << std::setprecision(5);
-  dt << "Acts::CylinderVolumeBounds: (rMin, rMax, halfZ, halfPhi, "
-        "averagePhi, minBevelZ, maxBevelZ) = ";
-  dt << get(eMinR) << ", " << get(eMaxR) << ", " << get(eHalfLengthZ) << ", "
-     << get(eHalfPhiSector) << ", " << get(eAveragePhi) << ", "
-     << get(eBevelMinZ) << ", " << get(eBevelMaxZ);
-  return dt;
-}
-
-inline std::vector<double> CylinderVolumeBounds::values() const {
-  std::vector<double> valvector;
-  valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
-  return valvector;
-}
-
-inline void CylinderVolumeBounds::checkConsistency() noexcept(false) {
-  if (get(eMinR) < 0. || get(eMaxR) <= 0. || get(eMinR) >= get(eMaxR)) {
-    throw std::invalid_argument("CylinderVolumeBounds: invalid radial input.");
-  }
-  if (get(eHalfLengthZ) <= 0) {
-    throw std::invalid_argument(
-        "CylinderVolumeBounds: invalid longitudinal input.");
-  }
-  if (get(eHalfPhiSector) < 0. || get(eHalfPhiSector) > M_PI) {
-    throw std::invalid_argument(
-        "CylinderVolumeBounds: invalid phi sector setup.");
-  }
-  if (get(eAveragePhi) != detail::radian_sym(get(eAveragePhi))) {
-    throw std::invalid_argument(
-        "CylinderVolumeBounds: invalid phi positioning.");
-  }
-  if (get(eBevelMinZ) != detail::radian_sym(get(eBevelMinZ))) {
-    throw std::invalid_argument("CylinderBounds: invalid bevel at min Z.");
-  }
-  if (get(eBevelMaxZ) != detail::radian_sym(get(eBevelMaxZ))) {
-    throw std::invalid_argument("CylinderBounds: invalid bevel at max Z.");
-  }
-}
-
 }  // namespace Acts
diff --git a/Core/include/Acts/Geometry/CylinderVolumeBuilder.hpp b/Core/include/Acts/Geometry/CylinderVolumeBuilder.hpp
index 8b242e139f8..c2b0e52c2ea 100644
--- a/Core/include/Acts/Geometry/CylinderVolumeBuilder.hpp
+++ b/Core/include/Acts/Geometry/CylinderVolumeBuilder.hpp
@@ -527,7 +527,8 @@ class CylinderVolumeBuilder : public ITrackingVolumeBuilder {
   ///         processing
   MutableTrackingVolumePtr trackingVolume(
       const GeometryContext& gctx, TrackingVolumePtr existingVolume = nullptr,
-      VolumeBoundsPtr externalBounds = nullptr) const override;
+      std::shared_ptr<const VolumeBounds> externalBounds =
+          nullptr) const override;
 
   /// Set configuration method
   ///
diff --git a/Core/include/Acts/Geometry/CylinderVolumeHelper.hpp b/Core/include/Acts/Geometry/CylinderVolumeHelper.hpp
index 5762552d7ec..84ce33f41cb 100644
--- a/Core/include/Acts/Geometry/CylinderVolumeHelper.hpp
+++ b/Core/include/Acts/Geometry/CylinderVolumeHelper.hpp
@@ -79,7 +79,8 @@ class CylinderVolumeHelper : public ITrackingVolumeHelper {
   MutableTrackingVolumePtr createTrackingVolume(
       const GeometryContext& gctx, const LayerVector& layers,
       std::shared_ptr<const IVolumeMaterial> volumeMaterial,
-      VolumeBoundsPtr volumeBounds, MutableTrackingVolumeVector mtvVector = {},
+      std::shared_ptr<const VolumeBounds> volumeBounds,
+      MutableTrackingVolumeVector mtvVector = {},
       const Transform3& transform = Transform3::Identity(),
       const std::string& volumeName = "UndefinedVolume",
       BinningType bType = arbitrary) const override;
diff --git a/Core/include/Acts/Geometry/Extent.hpp b/Core/include/Acts/Geometry/Extent.hpp
index 9d650264ea3..5d22e39081a 100644
--- a/Core/include/Acts/Geometry/Extent.hpp
+++ b/Core/include/Acts/Geometry/Extent.hpp
@@ -13,7 +13,6 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Utilities/BinningType.hpp"
 #include "Acts/Utilities/Enumerate.hpp"
-#include "Acts/Utilities/Range1D.hpp"
 #include "Acts/Utilities/RangeXD.hpp"
 
 #include <array>
@@ -122,17 +121,17 @@ class Extent {
   /// @param bValue is the binning value to be returned
   ///
   /// @return a one dimensional arrange
-  Range1D<ActsScalar>& range(BinningValue bValue);
+  auto range(BinningValue bValue) { return m_range[bValue]; }
 
   /// Return the individual 1-dimensional range
   ///
   /// @param bValue is the binning value to be returned
   ///
   /// @return a one dimensional arrange
-  const Range1D<ActsScalar>& range(BinningValue bValue) const;
+  Range1D<ActsScalar> range(BinningValue bValue) const;
 
   /// Return the N-dimension range
-  const RangeXD<binValues, ActsScalar> range() const;
+  const RangeXD<binValues, ActsScalar>& range() const;
 
   /// Return an D-dimensional sub range according to the
   /// the given @param binValues
@@ -227,16 +226,11 @@ class Extent {
   std::array<std::vector<ActsScalar>, binValues> m_valueHistograms;
 };
 
-inline Range1D<ActsScalar>& Acts::Extent::range(BinningValue bValue) {
+inline Range1D<ActsScalar> Acts::Extent::range(BinningValue bValue) const {
   return m_range[bValue];
 }
 
-inline const Range1D<ActsScalar>& Acts::Extent::range(
-    BinningValue bValue) const {
-  return m_range[bValue];
-}
-
-inline const RangeXD<binValues, ActsScalar> Extent::range() const {
+inline const RangeXD<binValues, ActsScalar>& Extent::range() const {
   return m_range;
 }
 
diff --git a/Core/include/Acts/Geometry/GenericCuboidVolumeBounds.hpp b/Core/include/Acts/Geometry/GenericCuboidVolumeBounds.hpp
index d205f51d92a..138ab0d7692 100644
--- a/Core/include/Acts/Geometry/GenericCuboidVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/GenericCuboidVolumeBounds.hpp
@@ -89,6 +89,14 @@ class GenericCuboidVolumeBounds : public VolumeBounds {
                                   const Vector3& envelope = {0, 0, 0},
                                   const Volume* entity = nullptr) const final;
 
+  /// Get the canonical binning values, i.e. the binning values
+  /// for that fully describe the shape's extent
+  ///
+  /// @return vector of canonical binning values
+  std::vector<Acts::BinningValue> canonicalBinning() const override {
+    return {Acts::binX, Acts::binY, Acts::binZ};
+  };
+
   /// @param sl is the output stream to be written into
   std::ostream& toStream(std::ostream& sl) const override;
 
diff --git a/Core/include/Acts/Geometry/ITrackingVolumeBuilder.hpp b/Core/include/Acts/Geometry/ITrackingVolumeBuilder.hpp
index d858ab76817..eb477c203f8 100644
--- a/Core/include/Acts/Geometry/ITrackingVolumeBuilder.hpp
+++ b/Core/include/Acts/Geometry/ITrackingVolumeBuilder.hpp
@@ -22,7 +22,6 @@ class Layer;
 class Volume;
 using TrackingVolumePtr = std::shared_ptr<const TrackingVolume>;
 using MutableTrackingVolumePtr = std::shared_ptr<TrackingVolume>;
-using VolumeBoundsPtr = std::shared_ptr<const VolumeBounds>;
 using LayerPtr = std::shared_ptr<const Layer>;
 using LayerVector = std::vector<LayerPtr>;
 ///  @class ITrackingVolumeBuilder
@@ -51,7 +50,7 @@ class ITrackingVolumeBuilder {
   /// @return shared pointer to a newly created TrackingVolume
   virtual MutableTrackingVolumePtr trackingVolume(
       const GeometryContext& gctx, TrackingVolumePtr oppositeVolume = nullptr,
-      VolumeBoundsPtr outsideBounds = nullptr) const = 0;
+      std::shared_ptr<const VolumeBounds> outsideBounds = nullptr) const = 0;
 };
 
 }  // namespace Acts
diff --git a/Core/include/Acts/Geometry/ITrackingVolumeHelper.hpp b/Core/include/Acts/Geometry/ITrackingVolumeHelper.hpp
index 664915b519f..df3e75f2549 100644
--- a/Core/include/Acts/Geometry/ITrackingVolumeHelper.hpp
+++ b/Core/include/Acts/Geometry/ITrackingVolumeHelper.hpp
@@ -26,7 +26,6 @@ class IVolumeMaterial;
 using LayerPtr = std::shared_ptr<const Layer>;
 using TrackingVolumePtr = std::shared_ptr<const TrackingVolume>;
 using MutableTrackingVolumePtr = std::shared_ptr<TrackingVolume>;
-using VolumeBoundsPtr = std::shared_ptr<const VolumeBounds>;
 
 using LayerVector = std::vector<LayerPtr>;
 using TrackingVolumeVector = std::vector<TrackingVolumePtr>;
@@ -64,7 +63,8 @@ class ITrackingVolumeHelper {
   virtual MutableTrackingVolumePtr createTrackingVolume(
       const GeometryContext& gctx, const LayerVector& layers,
       std::shared_ptr<const IVolumeMaterial> volumeMaterial,
-      VolumeBoundsPtr volumeBounds, MutableTrackingVolumeVector mtvVector = {},
+      std::shared_ptr<const VolumeBounds> volumeBounds,
+      MutableTrackingVolumeVector mtvVector = {},
       const Transform3& transform = Transform3::Identity(),
       const std::string& volumeName = "UndefinedVolume",
       BinningType btype = arbitrary) const = 0;
diff --git a/Core/include/Acts/Geometry/TrackingGeometryBuilder.hpp b/Core/include/Acts/Geometry/TrackingGeometryBuilder.hpp
index e7f20ac649a..a3d1bf4a011 100644
--- a/Core/include/Acts/Geometry/TrackingGeometryBuilder.hpp
+++ b/Core/include/Acts/Geometry/TrackingGeometryBuilder.hpp
@@ -41,7 +41,7 @@ class TrackingGeometryBuilder : public ITrackingGeometryBuilder {
     /// The list of tracking volume builders
     std::vector<std::function<std::shared_ptr<TrackingVolume>(
         const GeometryContext& gctx, const TrackingVolumePtr&,
-        const VolumeBoundsPtr&)>>
+        const std::shared_ptr<const VolumeBounds>&)>>
         trackingVolumeBuilders;
 
     /// The tracking volume helper for detector construction
diff --git a/Core/include/Acts/Geometry/TrackingVolume.hpp b/Core/include/Acts/Geometry/TrackingVolume.hpp
index 00a4e020dd4..2b12114c42f 100644
--- a/Core/include/Acts/Geometry/TrackingVolume.hpp
+++ b/Core/include/Acts/Geometry/TrackingVolume.hpp
@@ -125,7 +125,8 @@ class TrackingVolume : public Volume {
   ///
   /// @return shared pointer to a new TrackingVolume
   static MutableTrackingVolumePtr create(
-      const Transform3& transform, VolumeBoundsPtr volumeBounds,
+      const Transform3& transform,
+      std::shared_ptr<const VolumeBounds> volumeBounds,
       const std::shared_ptr<const TrackingVolumeArray>& containedVolumes =
           nullptr,
       const std::string& volumeName = "undefined") {
@@ -146,7 +147,8 @@ class TrackingVolume : public Volume {
   ///
   /// @return shared pointer to a new TrackingVolume
   static MutableTrackingVolumePtr create(
-      const Transform3& transform, VolumeBoundsPtr volumeBounds,
+      const Transform3& transform,
+      std::shared_ptr<const VolumeBounds> volumeBounds,
       std::shared_ptr<const IVolumeMaterial> volumeMaterial,
       std::unique_ptr<const LayerArray> containedLayers = nullptr,
       std::shared_ptr<const TrackingVolumeArray> containedVolumes = nullptr,
@@ -410,14 +412,6 @@ class TrackingVolume : public Volume {
   ///  - positiveFaceXY
   GlueVolumesDescriptor& glueVolumesDescriptor();
 
-  /// Register the color code
-  ///
-  /// @param icolor is a color number
-  void registerColorCode(unsigned int icolor);
-
-  /// Get the color code
-  unsigned int colorCode() const;
-
   /// Return the MotherVolume - if it exists
   const TrackingVolume* motherVolume() const;
 
@@ -435,7 +429,8 @@ class TrackingVolume : public Volume {
   /// @param volbounds is the description of the volume boundaries
   /// @param containedVolumeArray are the static volumes that fill this volume
   /// @param volumeName is a string identifier
-  TrackingVolume(const Transform3& transform, VolumeBoundsPtr volbounds,
+  TrackingVolume(const Transform3& transform,
+                 std::shared_ptr<const VolumeBounds> volbounds,
                  const std::shared_ptr<const TrackingVolumeArray>&
                      containedVolumeArray = nullptr,
                  const std::string& volumeName = "undefined");
@@ -452,7 +447,8 @@ class TrackingVolume : public Volume {
   /// @param denseVolumeVector  The contained dense volumes
   /// @param volumeName is a string identifier
   TrackingVolume(
-      const Transform3& transform, VolumeBoundsPtr volumeBounds,
+      const Transform3& transform,
+      std::shared_ptr<const VolumeBounds> volumeBounds,
       std::shared_ptr<const IVolumeMaterial> volumeMaterial,
       std::unique_ptr<const LayerArray> staticLayerArray = nullptr,
       std::shared_ptr<const TrackingVolumeArray> containedVolumeArray = nullptr,
@@ -515,9 +511,6 @@ class TrackingVolume : public Volume {
 
   /// Volume name for debug reasons & screen output
   std::string m_name;
-
-  /// color code for displaying
-  unsigned int m_colorCode{20};
 };
 
 inline const std::string& TrackingVolume::volumeName() const {
@@ -551,14 +544,6 @@ TrackingVolume::confinedVolumes() const {
   return m_confinedVolumes;
 }
 
-inline void TrackingVolume::registerColorCode(unsigned int icolor) {
-  m_colorCode = icolor;
-}
-
-inline unsigned int TrackingVolume::colorCode() const {
-  return m_colorCode;
-}
-
 inline const TrackingVolume* TrackingVolume::motherVolume() const {
   return m_motherVolume;
 }
diff --git a/Core/include/Acts/Geometry/TrapezoidVolumeBounds.hpp b/Core/include/Acts/Geometry/TrapezoidVolumeBounds.hpp
index d4999183ca1..bc00685e228 100644
--- a/Core/include/Acts/Geometry/TrapezoidVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/TrapezoidVolumeBounds.hpp
@@ -68,8 +68,8 @@ class TrapezoidVolumeBounds : public VolumeBounds {
   /// @param maxhalex is the half length in x at maximal y
   /// @param haley is the half length in y
   /// @param halez is the half length in z
-  TrapezoidVolumeBounds(double minhalex, double maxhalex, double haley,
-                        double halez) noexcept(false);
+  TrapezoidVolumeBounds(ActsScalar minhalex, ActsScalar maxhalex,
+                        ActsScalar haley, ActsScalar halez) noexcept(false);
 
   /// Constructor - the trapezoid boundaries (arbitrary trapezoid)
   ///
@@ -78,13 +78,14 @@ class TrapezoidVolumeBounds : public VolumeBounds {
   /// @param halez is the half length in z
   /// @param alpha is the opening angle at -x,-y
   /// @param beta is the opening angle at +x,-y
-  TrapezoidVolumeBounds(double minhalex, double haley, double halez,
-                        double alpha, double beta) noexcept(false);
+  TrapezoidVolumeBounds(ActsScalar minhalex, ActsScalar haley, ActsScalar halez,
+                        ActsScalar alpha, ActsScalar beta) noexcept(false);
 
   /// Constructor - from a fixed size array
   ///
   /// @param values The bound values
-  TrapezoidVolumeBounds(const std::array<double, eSize>& values) noexcept(false)
+  TrapezoidVolumeBounds(const std::array<ActsScalar, eSize>& values) noexcept(
+      false)
       : m_values(values) {
     checkConsistency();
     buildSurfaceBounds();
@@ -102,7 +103,7 @@ class TrapezoidVolumeBounds : public VolumeBounds {
   /// Return the bound values as dynamically sized vector
   ///
   /// @return this returns a copy of the internal values
-  std::vector<double> values() const final;
+  std::vector<ActsScalar> values() const final;
 
   /// This method checks if position in the 3D volume frame
   /// is inside the cylinder
@@ -111,7 +112,7 @@ class TrapezoidVolumeBounds : public VolumeBounds {
   /// @param tol is the tolerance applied
   ///
   /// @return boolean indicator if position is inside
-  bool inside(const Vector3& pos, double tol = 0.) const override;
+  bool inside(const Vector3& pos, ActsScalar tol = 0.) const override;
 
   /// Oriented surfaces, i.e. the decomposed boundary surfaces and the
   /// according navigation direction into the volume given the normal
@@ -136,15 +137,16 @@ class TrapezoidVolumeBounds : public VolumeBounds {
                                   const Volume* entity = nullptr) const final;
 
   /// Output Method for std::ostream
-  std::ostream& toStream(std::ostream& sl) const override;
+  /// @param os is the output stream
+  std::ostream& toStream(std::ostream& os) const override;
 
   /// Access to the bound values
   /// @param bValue the class nested enum for the array access
-  double get(BoundValues bValue) const { return m_values[bValue]; }
+  ActsScalar get(BoundValues bValue) const { return m_values[bValue]; }
 
  private:
-  /// The internal version of the bounds can be float/double
-  std::array<double, eSize> m_values{};
+  /// The internal version of the bounds can be float/ActsScalar
+  std::array<ActsScalar, eSize> m_values{};
   /// The face PlaneSurface parallel to local xy plane
   std::shared_ptr<const TrapezoidBounds> m_faceXYTrapezoidBounds{nullptr};
   /// Thhe face PlaneSurface attached to alpha (negative local x)
@@ -162,43 +164,6 @@ class TrapezoidVolumeBounds : public VolumeBounds {
 
   /// Helper method to create the surface bounds
   void buildSurfaceBounds();
-
-  /// Templated dump methods
-  /// @tparam stream_t The type of the stream for dumping
-  /// @param dt The stream object
-  template <class stream_t>
-  stream_t& dumpT(stream_t& dt) const;
 };
 
-template <class stream_t>
-stream_t& TrapezoidVolumeBounds::dumpT(stream_t& dt) const {
-  dt << std::setiosflags(std::ios::fixed);
-  dt << std::setprecision(5);
-  dt << "Acts::TrapezoidVolumeBounds: (minhalfX, halfY, halfZ, alpha, beta) "
-        "= ";
-  dt << "(" << get(eHalfLengthXnegY) << ", " << get(eHalfLengthXposY) << ", "
-     << get(eHalfLengthY) << ", " << get(eHalfLengthZ);
-  dt << ", " << get(eAlpha) << ", " << get(eBeta) << ")";
-  return dt;
-}
-
-inline std::vector<double> TrapezoidVolumeBounds::values() const {
-  std::vector<double> valvector;
-  valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
-  return valvector;
-}
-
-inline void TrapezoidVolumeBounds::checkConsistency() noexcept(false) {
-  if (get(eHalfLengthXnegY) < 0. || get(eHalfLengthXposY) < 0.) {
-    throw std::invalid_argument(
-        "TrapezoidVolumeBounds: invalid trapezoid parameters in x.");
-  }
-  if (get(eHalfLengthY) <= 0.) {
-    throw std::invalid_argument("TrapezoidVolumeBounds: invalid y extrusion.");
-  }
-  if (get(eHalfLengthZ) <= 0.) {
-    throw std::invalid_argument("TrapezoidVolumeBounds: invalid z extrusion.");
-  }
-}
-
 }  // namespace Acts
diff --git a/Core/include/Acts/Geometry/Volume.hpp b/Core/include/Acts/Geometry/Volume.hpp
index c8159118e11..037bfdfaeb7 100644
--- a/Core/include/Acts/Geometry/Volume.hpp
+++ b/Core/include/Acts/Geometry/Volume.hpp
@@ -20,7 +20,6 @@
 namespace Acts {
 
 class VolumeBounds;
-using VolumeBoundsPtr = std::shared_ptr<const VolumeBounds>;
 
 /// @class Volume
 ///
@@ -40,7 +39,8 @@ class Volume : public virtual GeometryObject {
   /// @param volbounds is the volume boundary definitions
   /// @note This will automatically build an oriented bounding box with an
   /// envelope value of (0.05, 0.05, 0.05)mm
-  Volume(const Transform3& transform, VolumeBoundsPtr volbounds);
+  Volume(const Transform3& transform,
+         std::shared_ptr<const VolumeBounds> volbounds);
 
   /// Copy Constructor - with optional shift
   ///
@@ -70,6 +70,9 @@ class Volume : public virtual GeometryObject {
   /// returns the volumeBounds()
   const VolumeBounds& volumeBounds() const;
 
+  /// Set volume bounds and update volume bounding boxes implicitly
+  void assignVolumeBounds(std::shared_ptr<const VolumeBounds> volbounds);
+
   /// Construct bounding box for this shape
   /// @param envelope Optional envelope to add / subtract from min/max
   /// @return Constructed bounding box pointing to this volume
@@ -85,7 +88,7 @@ class Volume : public virtual GeometryObject {
   /// @param tol is the tolerance parameter
   ///
   /// @return boolean indicator if the position is inside
-  bool inside(const Vector3& gpos, double tol = 0.) const;
+  bool inside(const Vector3& gpos, ActsScalar tol = 0.) const;
 
   /// The binning position method
   /// - as default the center is given, but may be overloaded
@@ -101,26 +104,10 @@ class Volume : public virtual GeometryObject {
   Transform3 m_transform;
   Transform3 m_itransform;
   Vector3 m_center;
-  VolumeBoundsPtr m_volumeBounds;
+  std::shared_ptr<const VolumeBounds> m_volumeBounds;
   BoundingBox m_orientedBoundingBox;
 };
 
-inline const Transform3& Volume::transform() const {
-  return m_transform;
-}
-
-inline const Transform3& Volume::itransform() const {
-  return m_itransform;
-}
-
-inline const Vector3& Volume::center() const {
-  return m_center;
-}
-
-inline const VolumeBounds& Volume::volumeBounds() const {
-  return (*(m_volumeBounds.get()));
-}
-
 /**Overload of << operator for std::ostream for debug output*/
 std::ostream& operator<<(std::ostream& sl, const Volume& vol);
 
diff --git a/Core/include/Acts/Geometry/VolumeBounds.hpp b/Core/include/Acts/Geometry/VolumeBounds.hpp
index 860977172f7..23bdcd4f819 100644
--- a/Core/include/Acts/Geometry/VolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/VolumeBounds.hpp
@@ -26,8 +26,6 @@ class Surface;
 class VolumeBounds;
 class Direction;
 
-using VolumeBoundsPtr = std::shared_ptr<const VolumeBounds>;
-
 using OrientedSurface = std::pair<std::shared_ptr<RegularSurface>, Direction>;
 using OrientedSurfaces = std::vector<OrientedSurface>;
 
@@ -115,6 +113,18 @@ class VolumeBounds {
       const Transform3* trf = nullptr, const Vector3& envelope = {0, 0, 0},
       const Volume* entity = nullptr) const = 0;
 
+  /// Get the canonical binning values, i.e. the binning values
+  /// for that fully describe the shape's extent
+  ///
+  /// @return vector of canonical binning values
+  ///
+  /// @note This is the default implementation that
+  /// returns the bounding box binning. Individual shapes
+  /// should override this method
+  virtual std::vector<Acts::BinningValue> canonicalBinning() const {
+    return {Acts::binX, Acts::binY, Acts::binZ};
+  };
+
   /// Binning offset - overloaded for some R-binning types
   ///
   /// @param bValue is the binning schema used
diff --git a/Core/include/Acts/Utilities/DBScan.hpp b/Core/include/Acts/Utilities/DBScan.hpp
index f978dde06a2..79f862ce573 100644
--- a/Core/include/Acts/Utilities/DBScan.hpp
+++ b/Core/include/Acts/Utilities/DBScan.hpp
@@ -159,8 +159,8 @@ class DBScan {
       // hypercube with a length of 2 epsilon).
       typename Tree::range_t range;
       for (std::size_t dim = 0; dim < kDims; dim++) {
-        range[dim] = std::make_pair(currentPoint[dim] - m_eps,
-                                    currentPoint[dim] + m_eps);
+        range[dim] = Range1D<scalar_t>(currentPoint[dim] - m_eps,
+                                       currentPoint[dim] + m_eps);
       }
       // We use the KDTree to find the neighbours.
       // An extra cut needs to be applied to only keep the neighbours that
diff --git a/Core/include/Acts/Utilities/KDTree.hpp b/Core/include/Acts/Utilities/KDTree.hpp
index d78f8cf021c..ddf7e1877af 100644
--- a/Core/include/Acts/Utilities/KDTree.hpp
+++ b/Core/include/Acts/Utilities/KDTree.hpp
@@ -49,7 +49,7 @@ class KDTree {
   using value_t = Type;
 
   /// @brief The type describing a multi-dimensional orthogonal range.
-  using range_t = RangeXD<Dims, Scalar, Vector>;
+  using range_t = RangeXD<Dims, Scalar>;
 
   /// @brief The type of coordinates for points.
   using coordinate_t = Vector<Scalar, Dims>;
@@ -303,7 +303,7 @@ class KDTree {
     range_t r;
 
     for (std::size_t j = 0; j < Dims; ++j) {
-      r[j] = {min_v[j], nextRepresentable(max_v[j])};
+      r[j] = Range1D<Scalar>{min_v[j], nextRepresentable(max_v[j])};
     }
 
     return r;
diff --git a/Core/include/Acts/Utilities/Range1D.hpp b/Core/include/Acts/Utilities/Range1D.hpp
deleted file mode 100644
index 3a9e24f7d64..00000000000
--- a/Core/include/Acts/Utilities/Range1D.hpp
+++ /dev/null
@@ -1,323 +0,0 @@
-// This file is part of the Acts project.
-//
-// Copyright (C) 2021 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 http://mozilla.org/MPL/2.0/.
-
-#pragma once
-
-#include <algorithm>
-#include <limits>
-
-namespace Acts {
-/// @brief A one-dimensional range between two points.
-///
-/// This type describes a one-demensional range of values, designed to be used
-/// in the construction of more complex multi-dimensional types. The class
-/// provides functionality for growing and shrinking the range, as well as some
-/// other utilities. These ranges are half-open, including the minimum but
-/// excluding the maximum.
-///
-/// @tparam Type The scalar type of the values contained in this range.
-template <typename Type>
-class Range1D {
- public:
-  /// @brief Construct a new degenerate range object
-  ///
-  /// This constructor coonstructs a degenerate range object with a maximum
-  /// lower than the minimum. In other words, this range is empty.
-  Range1D()
-      : m_min(std::numeric_limits<Type>::lowest()),
-        m_max(std::numeric_limits<Type>::max()) {}
-
-  /// @brief Construct a new range object from a lower and upper bound
-  ///
-  /// Construct a new range object given the values for the minimum and
-  /// maximum. Note that it is perfectly possible to construct a degenerate
-  /// range in this way.
-  ///
-  /// @param min The minimum value in the range (inclusive)
-  /// @param max The maximum value in the range (exclusive)
-  Range1D(Type min, Type max) : m_min(min), m_max(max) {}
-
-  /// @brief Construct a new range object from a pair of bounds
-  ///
-  /// Construct a new range object from a pair of values, the first of which
-  /// is taken to be the minimum, and the second of which is taken to be the
-  /// maximum.
-  ///
-  /// @param p The pair of values to use as the minimum and maximum
-  Range1D(const std::pair<Type, Type>& p) : m_min(p.first), m_max(p.second) {}
-
-  /// @brief Construct a new range object from an existing range
-  ///
-  /// This simply copies the values from the existing range to the new one.
-  /// It's the copy constructor.
-  ///
-  /// @param o The range to copy
-  Range1D(const Range1D<Type>& o) : m_min(o.min()), m_max(o.max()) {}
-
-  /// @brief Set the minimum value
-  ///
-  /// Override the minimum value of the range, regardless of what was already
-  /// set.
-  ///
-  /// @note If you want to shrink or expand the range, use the shrink and
-  /// expand methods.
-  ///
-  /// @param v The value to use as the new minimum
-  void setMin(const Type& v) { m_min = v; }
-
-  /// @brief Set the maximum value
-  ///
-  /// Override the maximum value of the range, regardless of what was already
-  /// set.
-  ///
-  /// @note If you want to shrink or expand the range, use the shrink and
-  /// expand methods.
-  ///
-  /// @param v The value to use as the new maximum
-  void setMax(const Type& v) { m_max = v; }
-
-  /// @brief Set the minimum and maximum value
-  ///
-  /// Override both the minimum and maximum value of the range, regardless of
-  /// what they were set to.
-  ///
-  /// @note If you want to shrink or expand the range, use the shrink and
-  /// expand methods.
-  ///
-  /// @note After this operation, the range should be exactly equal to [min,
-  /// max]
-  ///
-  /// @param min The new minimum value of the range
-  /// @param max The new maximum value of the range
-  void set(const Type& min, const Type& max) {
-    m_min = min;
-    m_max = max;
-  }
-
-  /// @brief Shrink a range by increasing the minimum value
-  ///
-  /// Shrink the range by increasing the minimum value. If the given value is
-  /// smaller than the current minimum (in other words, if the proposed new
-  /// range would be larger than the current range), this is a no-op.
-  ///
-  /// @param v The proposed new minimum for the range
-  void shrinkMin(const Type& v) { m_min = std::max(m_min, v); }
-
-  /// @brief Shrink a range by decreasing the maximum value
-  ///
-  /// Shrink the range by decreasing the maximum value. If the given value is
-  /// larger than the current maximum (in other words, if the proposed new
-  /// range would be larger than the current range), this is a no-op.
-  ///
-  /// @param v The proposed new maximum for the range
-  void shrinkMax(const Type& v) { m_max = std::min(m_max, v); }
-
-  /// @brief Shrink a range on both ends
-  ///
-  /// Shrink a range by increasing the minimum value as well as decreasing the
-  /// maximum value. If either of the values are already smaller or larger
-  /// (respectively) than the proposed values, then that particular boundary of
-  /// the interval is not shrunk.
-  ///
-  /// @note After this operation, the range is always equal to or smaller than
-  /// [min, max].
-  ///
-  /// @param min The proposed new minimum for the range
-  /// @param max The proposed new maximum for the range
-  void shrink(const Type& min, const Type& max) {
-    shrinkMin(min);
-    shrinkMax(max);
-  }
-
-  /// @brief Expand a range by decreasing the minimum value
-  ///
-  /// Expand the range by decreasing the minimum value. If the given value is
-  /// larger than the current minimum (in other words, if the proposed new
-  /// range would be smaller than the current range), this is a no-op.
-  ///
-  /// @param v The proposed new minimum for the range
-  void expandMin(const Type& v) { m_min = std::min(m_min, v); }
-
-  /// @brief Expand a range by increasing the maximum value
-  ///
-  /// Expand the range by increasing the maximum value. If the given value is
-  /// smaller than the current maximum (in other words, if the proposed new
-  /// range would be smaller than the current range), this is a no-op.
-  ///
-  /// @param v The proposed new maximum for the range
-  void expandMax(const Type& v) { m_max = std::max(m_max, v); }
-
-  /// @brief Expand a range on both ends
-  ///
-  /// Expand a range by decreasing the minimum value as well as increasing the
-  /// maximum value. If either of the values are already larger or smaller
-  /// (respectively) than the proposed values, then that particular boundary of
-  /// the interval is not expanded.
-  ///
-  /// @note After this operation, the range is always equal to or larger than
-  /// [min, max].
-  ///
-  /// @param min The proposed new minimum for the range
-  /// @param max The proposed new maximum for the range
-  void expand(const Type& min, const Type& max) {
-    expandMin(min);
-    expandMax(max);
-  }
-
-  /// @brief Return the minimum value of the range (inclusive)
-  Type min(void) const { return m_min; }
-
-  /// @brief Return the maximum value of the range (inclusive)
-  Type max(void) const { return m_max; }
-
-  /// @brief Compute the size of the range
-  ///
-  /// The size of a range is defined as the difference between the minimum and
-  /// the maximum. For degenerate ranges, this is zero.
-  ///
-  /// @warning Due to the nature of numbers, the result of this function can be
-  /// somewhat ambiguous. For natural numbers, you could argue that the range
-  /// [n, n] has size 0 or size 1. In this case we say it has size 0. The
-  /// uncountable nature of the reals means this doesn't matter for them, but
-  /// this can be awkward when working with integers.
-  ///
-  /// @return The size of the range
-  Type size(void) const {
-    return std::max(static_cast<Type>(0), m_max - m_min);
-  }
-
-  /// @brief Determine if this range is degenerate or not
-  ///
-  /// A degenerate range has a minimum higher than the maximum, and thus
-  /// cannot contain any values.
-  ///
-  /// @return true The range is degenerate and has size zero
-  /// @return false The range is not degenerate
-  bool degenerate(void) const { return m_min >= m_max; }
-
-  /// @brief Determine if the range contains a given value
-  ///
-  /// A value is inside a range if and only if it is greater than the minimum
-  /// and smaller than the maximum.
-  ///
-  /// @param v The value to check
-  ///
-  /// @return true The value is inside the range
-  /// @return false The value is not inside the range
-  bool contains(const Type& v) const { return m_min <= v && v < m_max; }
-
-  /// @brief Determine whether the range intersects another range
-  ///
-  /// The intersection of a range is the space where both ranges overlap. If
-  /// the ranges overlap at all, they are said to intersect. This operation
-  /// is commutative.
-  ///
-  /// @param o The other range to check
-  ///
-  /// @return true The ranges intersect
-  /// @return false The ranges do not intersect
-  bool operator&&(const Range1D<Type>& o) const {
-    return m_min < o.max() && o.min() < m_max;
-  }
-
-  /// @brief Determine whether the range is equal to another range
-  ///
-  /// Two ranges are equal if and only if their minima and maxima are the
-  /// same.
-  ///
-  /// @warning This method relies on the existence of a well-defined notion
-  /// of equality for the underlying types. Using this method on floating
-  /// ranges may have unintended effecrs.
-  ///
-  /// @param o The other range to check
-  ///
-  /// @return true The ranges are equal
-  /// @return false The ranges are not equal
-  bool operator==(const Range1D<Type>& o) const {
-    return min() == o.min() && max() == o.max();
-  }
-
-  /// @brief Determine whether the left-hand range is a subset of the
-  /// right-hand range
-  ///
-  /// A range is a subset of another range if and only if all values
-  /// contained in the first range are also contained in the second range.
-  ///
-  /// @param o The other range to check
-  ///
-  /// @return true The left-hand range is a subset of the right-hand range
-  /// @return false The left-hand range is not a subset of the right-hand
-  /// range
-  bool operator<=(const Range1D<Type>& o) const {
-    return min() >= o.min() && max() <= o.max();
-  }
-
-  /// @brief Determine whether the left-hand range is a superset of the
-  /// right-hand range
-  ///
-  /// A range is a superset of another range if and only if all values
-  /// contained in the second range are also contained in the first range.
-  ///
-  /// @param o The other range to check
-  ///
-  /// @return true The left-hand range is a superset of thr right-hand range
-  /// @return false The left-hand range is not a superset of the right-hand
-  /// range
-  bool operator>=(const Range1D<Type>& o) const {
-    return min() <= o.min() && max() >= o.max();
-  }
-
-  /// @brief Assignment operator
-  ///
-  /// Copy the right-hand range into the left-hand range, which means setting
-  /// the minimum and maximum to equal the minimum and maximum of the
-  /// right-hand side.
-  ///
-  /// @param o The range of values to copy
-  ///
-  /// @return This range
-  Range1D<Type>& operator=(const Range1D<Type>& o) {
-    m_min = o.min();
-    m_max = o.max();
-
-    return *this;
-  }
-
-  /// @brief Compute the intersection of two ranges
-  ///
-  /// The intersection of two ranges is the range containing all values
-  /// contained in both ranges. If the two ranges do not intersect, the
-  /// intersection is a degenerate range. This operation is commutative.
-  ///
-  /// @param o The range to compute the intersection with
-  ///
-  /// @return The intersection range between the two ranges
-  Range1D<Type> operator&(const Range1D<Type>& o) const {
-    return Range1D<Type>(std::max(m_min, o.min()), std::min(m_max, o.max()));
-  }
-
-  /// @brief Set the range to the intersection of itself and another range
-  ///
-  /// This is an assignment version of operator&, which updates the range on
-  /// which it is called to ensure that the new range is the intersection of
-  /// the old range and the new range.
-  ///
-  /// @param o The range to compute the intersection with
-  ///
-  /// @return This object
-  Range1D<Type>& operator&=(const Range1D<Type>& o) {
-    m_min = std::max(m_min, o.min());
-    m_max = std::min(m_max, o.max());
-
-    return *this;
-  }
-
- private:
-  Type m_min, m_max;
-};
-}  // namespace Acts
diff --git a/Core/include/Acts/Utilities/RangeXD.hpp b/Core/include/Acts/Utilities/RangeXD.hpp
index 2f6acaaa56f..d614004234e 100644
--- a/Core/include/Acts/Utilities/RangeXD.hpp
+++ b/Core/include/Acts/Utilities/RangeXD.hpp
@@ -8,10 +8,9 @@
 
 #pragma once
 
-#include "Acts/Utilities/Range1D.hpp"
-
 #include <algorithm>
 #include <array>
+#include <cassert>
 #include <limits>
 #include <sstream>
 #include <string>
@@ -29,20 +28,59 @@ namespace Acts {
 template <std::size_t Dims, typename Type,
           template <typename, std::size_t> typename Vector = std::array>
 class RangeXD {
+ private:
+  // @TODO: Replace with std::span or boost::span once available
+  template <typename, std::size_t>
+  struct SingleElementContainer {
+    Type* element;
+
+    Type& operator[](std::size_t i) {
+      (void)i;
+      assert(i == 0);
+
+      return *element;
+    }
+  };
+
  public:
-  /// @brief The type used to describe coordinates in our range
-  using coordinate_t = Vector<Type, Dims>;
+  RangeXD() {
+    for (std::size_t i = 0; i < Dims; ++i) {
+      min(i) = std::numeric_limits<Type>::lowest();
+      max(i) = std::numeric_limits<Type>::max();
+    }
+  }
+
+  /// @brief Construct a range from a pair of minimum and maximum values
+  /// @param minima The minimum values of the range
+  /// @param maxima The maximum values of the range
+  RangeXD(Vector<Type, Dims> minima, Vector<Type, Dims> maxima)
+      : m_minima(minima), m_maxima(maxima) {}
+
+  /// @brief Construct a range from a pair of single minimum and maximum values
+  /// @note Only available for one-dimensional ranges
+  /// @param minimum The minimum value of the range
+  /// @param maximum The maximum value of the range
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  RangeXD(Type minimum, Type maximum)
+      : m_minima({minimum}), m_maxima({maximum}) {}
+
+  /// @brief Construct a range from a pair of minimum and maximum values
+  /// @note Only available for one-dimensional ranges
+  /// @param p The pair of minimum and maximum values
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  RangeXD(const std::pair<Type, Type>& p)
+      : m_minima({p.first}), m_maxima({p.second}) {}
 
   /// @brief Determine whether this range is degenerate
   ///
-  /// A degenerate multi-dimensional range has no volume and cannot contain any
-  /// values. This is the case if any of its dimensions are degenerate.
+  /// A degenerate multi-dimensional range has no volume and cannot contain
+  /// any values. This is the case if any of its dimensions are degenerate.
   ///
   /// @return true The range is degenerate
   /// @return false The range is not degenerate
-  bool degenerate(void) const {
+  bool degenerate() const {
     for (std::size_t i = 0; i < Dims; ++i) {
-      if (m_dims[i].degenerate()) {
+      if (min(i) >= max(i)) {
         return true;
       }
     }
@@ -58,9 +96,10 @@ class RangeXD {
   ///
   /// @return true The coordinate is inside the range
   /// @return false The coordinate is outside the range
-  bool contains(const coordinate_t& v) const {
+  template <template <typename, std::size_t> typename coordinate_t = std::array>
+  bool contains(const coordinate_t<Type, Dims>& v) const {
     for (std::size_t i = 0; i < Dims; ++i) {
-      if (!m_dims[i].contains(v[i])) {
+      if (!(min(i) <= v[i] && v[i] < max(i))) {
         return false;
       }
     }
@@ -72,14 +111,35 @@ class RangeXD {
   ///
   /// @param i The index of the dimension to access
   /// @return A reference to the dimension contained in this range
-  Range1D<Type>& operator[](const std::size_t& i) { return m_dims[i]; }
+  RangeXD<1, Type, SingleElementContainer> operator[](const std::size_t& i) {
+    return RangeXD<1, Type, SingleElementContainer>{{&min(i)}, {&max(i)}};
+  }
 
   /// @brief Access one of the dimensional ranges of the volume
   ///
   /// @param i The index of the dimension to access
   /// @return A reference to the dimension contained in this range
-  const Range1D<Type>& operator[](const std::size_t& i) const {
-    return m_dims[i];
+  RangeXD<1, Type> operator[](const std::size_t& i) const {
+    return RangeXD<1, Type>{{min(i)}, {max(i)}};
+  }
+
+  /// @brief Assignment operator
+  ///
+  /// Copy the right-hand range into the left-hand range, which means setting
+  /// the minimum and maximum to equal the minimum and maximum of the
+  /// right-hand side.
+  ///
+  /// @param o The range of values to copy
+  ///
+  /// @return This range
+  template <template <typename, std::size_t> typename V>
+  RangeXD& operator=(const RangeXD<Dims, Type, V>& o) {
+    for (std::size_t i = 0; i < Dims; ++i) {
+      min(i) = o.min(i);
+      max(i) = o.max(i);
+    }
+
+    return *this;
   }
 
   /// @brief Determine whether two ranges are equal
@@ -93,7 +153,7 @@ class RangeXD {
   /// @return false The ranges are not equal
   bool operator==(const RangeXD<Dims, Type, Vector>& o) const {
     for (std::size_t i = 0; i < Dims; ++i) {
-      if (!(this->m_dims[i] == o[i])) {
+      if (!(min(i) == o.min(i) && max(i) == o.max(i))) {
         return false;
       }
     }
@@ -115,7 +175,7 @@ class RangeXD {
   /// @return false The first range is not a subset of the second range
   bool operator<=(const RangeXD<Dims, Type, Vector>& o) const {
     for (std::size_t i = 0; i < Dims; ++i) {
-      if (!(this->m_dims[i] <= o[i])) {
+      if (!(min(i) >= o.min(i) && max(i) <= o.max(i))) {
         return false;
       }
     }
@@ -138,7 +198,7 @@ class RangeXD {
   /// range
   bool operator>=(const RangeXD<Dims, Type, Vector>& o) const {
     for (std::size_t i = 0; i < Dims; ++i) {
-      if (!(this->m_dims[i] >= o[i])) {
+      if (!(min(i) <= o.min(i) && max(i) >= o.max(i))) {
         return false;
       }
     }
@@ -148,8 +208,8 @@ class RangeXD {
 
   /// @brief Compute the intersection of this range with another range
   ///
-  /// The intersection of one orthogonal range with another orthogonal range is
-  /// in itself an orthogonal range. This operation is commutative. This
+  /// The intersection of one orthogonal range with another orthogonal range
+  /// is in itself an orthogonal range. This operation is commutative. This
   /// intersection between two n-dimensional ranges is defined simply as the
   /// intersection in each dimension of the two ranges.
   ///
@@ -161,7 +221,8 @@ class RangeXD {
     RangeXD<Dims, Type> res;
 
     for (std::size_t i = 0; i < Dims; ++i) {
-      res[i] = m_dims[i] & o[i];
+      res.min(i) = std::max(min(i), o.min(i));
+      res.max(i) = std::min(max(i), o.max(i));
     }
 
     return res;
@@ -179,7 +240,8 @@ class RangeXD {
   RangeXD<Dims, Type, Vector>& operator&=(
       const RangeXD<Dims, Type, Vector>& o) {
     for (std::size_t i = 0; i < Dims; ++i) {
-      m_dims[i] &= o[i];
+      min(i) = std::max(min(i), o.min(i));
+      max(i) = std::min(max(i), o.max(i));
     }
 
     return *this;
@@ -187,8 +249,8 @@ class RangeXD {
 
   /// @brief Determine whether this range intersects another
   ///
-  /// Two n-dimensional ranges intersect if and only if they intersect in every
-  /// one of their n dimensions. Otherwise, they are disjoint.
+  /// Two n-dimensional ranges intersect if and only if they intersect in
+  /// every one of their n dimensions. Otherwise, they are disjoint.
   ///
   /// @param r The other range to check
   ///
@@ -196,7 +258,7 @@ class RangeXD {
   /// @return false The ranges do not intersect
   bool operator&&(const RangeXD<Dims, Type, Vector>& r) const {
     for (std::size_t i = 0; i < Dims; ++i) {
-      if (!(m_dims[i] && r[i])) {
+      if (!(min(i) < r.max(i) && r.min(i) < max(i))) {
         return false;
       }
     }
@@ -206,15 +268,15 @@ class RangeXD {
 
   /// @brief Represent the range as a string
   ///
-  /// This method produces a helpful string that can be used to debug the range
-  /// if needed. Not really designed to be used in production code.
+  /// This method produces a helpful string that can be used to debug the
+  /// range if needed. Not really designed to be used in production code.
   ///
   /// @return A string representing the range
   std::string toString(void) const {
     std::stringstream s;
 
     for (std::size_t i = 0; i < Dims; ++i) {
-      s << m_dims[i].min() << " <= v[" << i << "] <= " << m_dims[i].max();
+      s << min(i) << " <= v[" << i << "] <= " << max(i);
       if (i != Dims - 1) {
         s << ", ";
       }
@@ -223,7 +285,339 @@ class RangeXD {
     return s.str();
   }
 
+  /// @brief Shrink a range by increasing the minimum value
+  ///
+  /// Shrink the range by increasing the minimum value. If the given value is
+  /// smaller than the current minimum (in other words, if the proposed new
+  /// range would be larger than the current range), this is a no-op.
+  ///
+  /// @param i The index of the dimension to shrink
+  /// @param v The proposed new minimum for the range
+  void shrinkMin(std::size_t i, const Type& v) { min(i) = std::max(min(i), v); }
+
+  /// @brief Shrink a range by decreasing the maximum value
+  ///
+  /// Shrink the range by decreasing the maximum value. If the given value is
+  /// larger than the current maximum (in other words, if the proposed new
+  /// range would be larger than the current range), this is a no-op.
+  ///
+  /// @param i The index of the dimension to shrink
+  /// @param v The proposed new maximum for the range
+  void shrinkMax(std::size_t i, const Type& v) { max(i) = std::min(max(i), v); }
+
+  /// @brief Shrink a range on both ends
+  ///
+  /// Shrink a range by increasing the minimum value as well as decreasing the
+  /// maximum value. If either of the values are already smaller or larger
+  /// (respectively) than the proposed values, then that particular boundary
+  /// of the interval is not shrunk.
+  ///
+  /// @note After this operation, the range is always equal to or smaller than
+  /// [min, max].
+  ///
+  /// @param i The index of the dimension to shrink
+  /// @param min The proposed new minimum for the range
+  /// @param max The proposed new maximum for the range
+  void shrink(std::size_t i, const Type& min, const Type& max) {
+    shrinkMin(i, min);
+    shrinkMax(i, max);
+  }
+
+  /// @brief Expand a range by decreasing the minimum value
+  ///
+  /// Expand the range by decreasing the minimum value. If the given value is
+  /// larger than the current minimum (in other words, if the proposed new
+  /// range would be smaller than the current range), this is a no-op.
+  ///
+  /// @param i The index of the dimension to expand
+  /// @param v The proposed new minimum for the range
+  void expandMin(std::size_t i, const Type& v) { min(i) = std::min(min(i), v); }
+
+  /// @brief Expand a range by increasing the maximum value
+  ///
+  /// Expand the range by increasing the maximum value. If the given value is
+  /// smaller than the current maximum (in other words, if the proposed new
+  /// range would be smaller than the current range), this is a no-op.
+  ///
+  /// @param i The index of the dimension to expand
+  /// @param v The proposed new maximum for the range
+  void expandMax(std::size_t i, const Type& v) { max(i) = std::max(max(i), v); }
+
+  /// @brief Expand a range on both ends
+  ///
+  /// Expand a range by decreasing the minimum value as well as increasing the
+  /// maximum value. If either of the values are already larger or smaller
+  /// (respectively) than the proposed values, then that particular boundary
+  /// of the interval is not expanded.
+  ///
+  /// @note After this operation, the range is always equal to or larger than
+  /// [min, max].
+  ///
+  /// @param i The index of the dimension to expand
+  /// @param min The proposed new minimum for the range
+  /// @param max The proposed new maximum for the range
+  void expand(std::size_t i, const Type& min, const Type& max) {
+    expandMin(i, min);
+    expandMax(i, max);
+  }
+
+  /// @brief Set the minimum value
+  ///
+  /// Override the minimum value of the range, regardless of what was already
+  /// set.
+  ///
+  /// @note If you want to shrink or expand the range, use the shrink and
+  /// expand methods.
+  ///
+  /// @param i The index of the dimension to set
+  /// @param v The value to use as the new minimum
+  void setMin(std::size_t i, const Type& v) { min(i) = v; }
+
+  /// @brief Set the maximum value
+  ///
+  /// Override the maximum value of the range, regardless of what was already
+  /// set.
+  ///
+  /// @note If you want to shrink or expand the range, use the shrink and
+  /// expand methods.
+  ///
+  /// @param i The index of the dimension to set
+  /// @param v The value to use as the new maximum
+  void setMax(std::size_t i, const Type& v) { max(i) = v; }
+
+  /// @brief Set the minimum and maximum value
+  ///
+  /// Override both the minimum and maximum value of the range, regardless of
+  /// what they were set to.
+  ///
+  /// @note If you want to shrink or expand the range, use the shrink and
+  /// expand methods.
+  ///
+  /// @note After this operation, the range should be exactly equal to [min,
+  /// max]
+  ///
+  /// @param i The index of the dimension to set
+  /// @param min The new minimum value of the range
+  /// @param max The new maximum value of the range
+  void set(std::size_t i, const Type& min, const Type& max) {
+    setMin(i, min);
+    setMax(i, max);
+  }
+
+  /// @brief Return the minimum value of the range @p i (inclusive)
+  /// @param i The index of the dimension to access
+  Type& min(std::size_t i) { return m_minima[i]; }
+
+  /// @brief Return the maximum value of the range @p i (inclusive)
+  /// @param i The index of the dimension to access
+  Type& max(std::size_t i) { return m_maxima[i]; }
+
+  /// @brief Return the minimum value of the range @p i (inclusive)
+  /// @param i The index of the dimension to access
+  Type min(std::size_t i) const { return m_minima[i]; }
+
+  /// @brief Return the maximum value of the range @p i (inclusive)
+  /// @param i The index of the dimension to access
+  Type max(std::size_t i) const { return m_maxima[i]; }
+
+  /// Methods for manipulating a range of dimension 1
+  /// @{
+
+  /// @brief Expand a range by decreasing the minimum value
+  ///
+  /// Expand the range by decreasing the minimum value. If the given value is
+  /// larger than the current minimum (in other words, if the proposed new
+  /// range would be smaller than the current range), this is a no-op.
+  ///
+  /// @param v The proposed new minimum for the range
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  void expandMin(const Type& v) {
+    min() = std::min(min(), v);
+  }
+
+  /// @brief Expand a range by increasing the maximum value
+  ///
+  /// Expand the range by increasing the maximum value. If the given value is
+  /// smaller than the current maximum (in other words, if the proposed new
+  /// range would be smaller than the current range), this is a no-op.
+  ///
+  /// @param v The proposed new maximum for the range
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  void expandMax(const Type& v) {
+    max() = std::max(max(), v);
+  }
+
+  /// @brief Expand a range on both ends
+  ///
+  /// Expand a range by decreasing the minimum value as well as increasing the
+  /// maximum value. If either of the values are already larger or smaller
+  /// (respectively) than the proposed values, then that particular boundary
+  /// of the interval is not expanded.
+  ///
+  /// @note After this operation, the range is always equal to or larger than
+  /// [min, max].
+  ///
+  /// @param min The proposed new minimum for the range
+  /// @param max The proposed new maximum for the range
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  void expand(const Type& min, const Type& max) {
+    expandMin(min);
+    expandMax(max);
+  }
+
+  /// @brief Shrink a range by increasing the minimum value
+  ///
+  /// Shrink the range by increasing the minimum value. If the given value is
+  /// smaller than the current minimum (in other words, if the proposed new
+  /// range would be larger than the current range), this is a no-op.
+  ///
+  /// @param v The proposed new minimum for the range
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  void shrinkMin(const Type& v) {
+    min() = std::max(min(), v);
+  }
+
+  /// @brief Shrink a range by decreasing the maximum value
+  ///
+  /// Shrink the range by decreasing the maximum value. If the given value is
+  /// larger than the current maximum (in other words, if the proposed new
+  /// range would be larger than the current range), this is a no-op.
+  ///
+  /// @param v The proposed new maximum for the range
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  void shrinkMax(const Type& v) {
+    max() = std::min(max(), v);
+  }
+
+  /// @brief Shrink a range on both ends
+  ///
+  /// Shrink a range by increasing the minimum value as well as decreasing the
+  /// maximum value. If either of the values are already smaller or larger
+  /// (respectively) than the proposed values, then that particular boundary
+  /// of the interval is not shrunk.
+  ///
+  /// @note After this operation, the range is always equal to or smaller than
+  /// [min, max].
+  ///
+  /// @param min The proposed new minimum for the range
+  /// @param max The proposed new maximum for the range
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  void shrink(const Type& min, const Type& max) {
+    shrinkMin(min);
+    shrinkMax(max);
+  }
+
+  /// @brief Set the minimum value
+  ///
+  /// Override the minimum value of the range, regardless of what was already
+  /// set.
+  ///
+  /// @note If you want to shrink or expand the range, use the shrink and
+  /// expand methods.
+  ///
+  /// @param v The value to use as the new minimum
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  void setMin(const Type& v) {
+    min() = v;
+  }
+
+  /// @brief Set the maximum value
+  ///
+  /// Override the maximum value of the range, regardless of what was already
+  /// set.
+  ///
+  /// @note If you want to shrink or expand the range, use the shrink and
+  /// expand methods.
+  ///
+  /// @param v The value to use as the new maximum
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  void setMax(const Type& v) {
+    max() = v;
+  }
+
+  /// @brief Set the minimum and maximum value
+  ///
+  /// Override both the minimum and maximum value of the range, regardless of
+  /// what they were set to.
+  ///
+  /// @note If you want to shrink or expand the range, use the shrink and
+  /// expand methods.
+  ///
+  /// @note After this operation, the range should be exactly equal to [min,
+  /// max]
+  ///
+  /// @param min The new minimum value of the range
+  /// @param max The new maximum value of the range
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  void set(const Type& min, const Type& max) {
+    setMin(min);
+    setMax(max);
+  }
+
+  /// @brief Return the minimum value of the range (inclusive)
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  Type min() const {
+    return min(0);
+  }
+
+  /// @brief Return the minimum value of the range (inclusive)
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  Type& min() {
+    return min(0);
+  }
+
+  /// @brief Return the maximum value of the range (inclusive)
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  Type max() const {
+    return max(0);
+  }
+
+  /// @brief Return the maximum value of the range (inclusive)
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  Type& max() {
+    return max(0);
+  }
+
+  /// @brief Compute the size of the range
+  ///
+  /// The size of a range is defined as the difference between the minimum and
+  /// the maximum. For degenerate ranges, this is zero.
+  ///
+  /// @warning Due to the nature of numbers, the result of this function can be
+  /// somewhat ambiguous. For natural numbers, you could argue that the range
+  /// [n, n] has size 0 or size 1. In this case we say it has size 0. The
+  /// uncountable nature of the reals means this doesn't matter for them, but
+  /// this can be awkward when working with integers.
+  ///
+  /// @return The size of the range
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  Type size() const {
+    return std::max(static_cast<Type>(0), max() - min());
+  }
+
+  /// @brief Determine if the range contains a given value
+  ///
+  /// A value is inside a range if and only if it is greater than the minimum
+  /// and smaller than the maximum.
+  ///
+  /// @param v The value to check
+  ///
+  /// @return true The value is inside the range
+  /// @return false The value is not inside the range
+  template <std::size_t I = Dims, typename = std::enable_if_t<I == 1>>
+  bool contains(const Type& v) const {
+    return min() <= v && v < max();
+  }
+
+  /// @}
+
  private:
-  std::array<Range1D<Type>, Dims> m_dims;
+  Vector<Type, Dims> m_minima{};
+  Vector<Type, Dims> m_maxima{};
 };
+
+template <typename Type,
+          template <typename, std::size_t> typename Vector = std::array>
+using Range1D = RangeXD<1, Type, Vector>;
+
 }  // namespace Acts
diff --git a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp
index 23f6759085c..0bf05086b67 100644
--- a/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp
+++ b/Core/include/Acts/Vertexing/AdaptiveMultiVertexFinder.hpp
@@ -75,6 +75,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder {
     // track as compatible to vertex. If useTime is set to true, the time
     // coordinate also contributes to the significance and tracksMaxSignificance
     // needs to be increased.
+    // 5 corresponds to a p-value of ~0.92 using `chi2(x=5,ndf=2)`
     double tracksMaxSignificance = 5.;
 
     // Max chi2 value for which tracks are considered compatible with
@@ -98,6 +99,7 @@ class AdaptiveMultiVertexFinder final : public IVertexFinder {
 
     // Maximum significance on the distance between two vertices
     // to allow merging of two vertices.
+    // 3 corresponds to a p-value of ~0.92 using `chi2(x=3,ndf=1)`
     double maxMergeVertexSignificance = 3.;
 
     // Minimum weight a track has to have to be considered a compatible
diff --git a/Core/src/Detector/DetectorVolume.cpp b/Core/src/Detector/DetectorVolume.cpp
index c233f36ae3b..c0edaf2d11d 100644
--- a/Core/src/Detector/DetectorVolume.cpp
+++ b/Core/src/Detector/DetectorVolume.cpp
@@ -55,7 +55,6 @@ Acts::Experimental::DetectorVolume::DetectorVolume(
   }
 
   [[maybe_unused]] const auto& gctx_ref = gctx;
-  assert(checkContainment(gctx) && "Objects are not contained by volume.");
 }
 
 Acts::Experimental::DetectorVolume::DetectorVolume(
@@ -264,20 +263,27 @@ Acts::Extent Acts::Experimental::DetectorVolume::extent(
 
 bool Acts::Experimental::DetectorVolume::checkContainment(
     const GeometryContext& gctx, std::size_t nseg) const {
+  // We don't have a logging instance here
+  // so can't throw a warning for shapes that are
+  // using the bounding box
+  auto binningValues = volumeBounds().canonicalBinning();
+
   // Create the volume extent
   auto volumeExtent = extent(gctx, nseg);
   // Check surfaces
-  for (const auto* s : surfaces()) {
-    auto sExtent = s->polyhedronRepresentation(gctx, nseg).extent();
-    if (!volumeExtent.contains(sExtent)) {
-      return false;
+  for (auto b : binningValues) {
+    for (const auto* s : surfaces()) {
+      auto sExtent = s->polyhedronRepresentation(gctx, nseg).extent();
+      if (!volumeExtent.contains(sExtent, b)) {
+        return false;
+      }
     }
-  }
-  // Check volumes
-  for (const auto* v : volumes()) {
-    auto vExtent = v->extent(gctx, nseg);
-    if (!volumeExtent.contains(vExtent)) {
-      return false;
+    // Check volumes
+    for (const auto* v : volumes()) {
+      auto vExtent = v->extent(gctx, nseg);
+      if (!volumeExtent.contains(vExtent, b)) {
+        return false;
+      }
     }
   }
   // All contained
diff --git a/Core/src/Geometry/ConeVolumeBounds.cpp b/Core/src/Geometry/ConeVolumeBounds.cpp
index c339cd317d6..879a590df97 100644
--- a/Core/src/Geometry/ConeVolumeBounds.cpp
+++ b/Core/src/Geometry/ConeVolumeBounds.cpp
@@ -28,10 +28,11 @@
 #include <type_traits>
 #include <utility>
 
-Acts::ConeVolumeBounds::ConeVolumeBounds(double innerAlpha, double innerOffsetZ,
-                                         double outerAlpha, double outerOffsetZ,
-                                         double halflengthZ, double averagePhi,
-                                         double halfPhiSector) noexcept(false)
+namespace Acts {
+ConeVolumeBounds::ConeVolumeBounds(
+    ActsScalar innerAlpha, ActsScalar innerOffsetZ, ActsScalar outerAlpha,
+    ActsScalar outerOffsetZ, ActsScalar halflengthZ, ActsScalar averagePhi,
+    ActsScalar halfPhiSector) noexcept(false)
     : VolumeBounds(), m_values() {
   m_values[eInnerAlpha] = innerAlpha;
   m_values[eInnerOffsetZ] = innerOffsetZ;
@@ -44,10 +45,10 @@ Acts::ConeVolumeBounds::ConeVolumeBounds(double innerAlpha, double innerOffsetZ,
   checkConsistency();
 }
 
-Acts::ConeVolumeBounds::ConeVolumeBounds(double cylinderR, double alpha,
-                                         double offsetZ, double halflengthZ,
-                                         double averagePhi,
-                                         double halfPhiSector) noexcept(false)
+ConeVolumeBounds::ConeVolumeBounds(ActsScalar cylinderR, ActsScalar alpha,
+                                   ActsScalar offsetZ, ActsScalar halflengthZ,
+                                   ActsScalar averagePhi,
+                                   ActsScalar halfPhiSector) noexcept(false)
     : VolumeBounds(), m_values() {
   m_values[eInnerAlpha] = 0.;
   m_values[eInnerOffsetZ] = 0.;
@@ -58,11 +59,11 @@ Acts::ConeVolumeBounds::ConeVolumeBounds(double cylinderR, double alpha,
   m_values[eHalfPhiSector] = halfPhiSector;
 
   // Cone parameters
-  double tanAlpha = std::tan(alpha);
-  double zmin = offsetZ - halflengthZ;
-  double zmax = offsetZ + halflengthZ;
-  double rmin = std::abs(zmin) * tanAlpha;
-  double rmax = std::abs(zmax) * tanAlpha;
+  ActsScalar tanAlpha = std::tan(alpha);
+  ActsScalar zmin = offsetZ - halflengthZ;
+  ActsScalar zmax = offsetZ + halflengthZ;
+  ActsScalar rmin = std::abs(zmin) * tanAlpha;
+  ActsScalar rmax = std::abs(zmax) * tanAlpha;
 
   if (rmin >= cylinderR) {
     // Cylindrical cut-out of a cone
@@ -170,7 +171,7 @@ Acts::OrientedSurfaces Acts::ConeVolumeBounds::orientedSurfaces(
   return oSurfaces;
 }
 
-void Acts::ConeVolumeBounds::checkConsistency() noexcept(false) {
+void ConeVolumeBounds::checkConsistency() noexcept(false) {
   if (innerRmin() > outerRmin() || innerRmax() > outerRmax()) {
     throw std::invalid_argument("ConeVolumeBounds: invalid radial input.");
   }
@@ -190,35 +191,36 @@ void Acts::ConeVolumeBounds::checkConsistency() noexcept(false) {
   }
 }
 
-bool Acts::ConeVolumeBounds::inside(const Vector3& pos, double tol) const {
-  double z = pos.z();
-  double zmin = z + tol;
-  double zmax = z - tol;
+bool ConeVolumeBounds::inside(const Vector3& pos, ActsScalar tol) const {
+  ActsScalar z = pos.z();
+  ActsScalar zmin = z + tol;
+  ActsScalar zmax = z - tol;
   // Quick check outside z
   if (zmin < -get(eHalfLengthZ) || zmax > get(eHalfLengthZ)) {
     return false;
   }
-  double r = VectorHelpers::perp(pos);
+  ActsScalar r = VectorHelpers::perp(pos);
   if (std::abs(get(eHalfPhiSector) - M_PI) > s_onSurfaceTolerance) {
     // need to check the phi sector - approximate phi tolerance
-    double phitol = tol / r;
-    double phi = VectorHelpers::phi(pos);
-    double phimin = phi - phitol;
-    double phimax = phi + phitol;
+    ActsScalar phitol = tol / r;
+    ActsScalar phi = VectorHelpers::phi(pos);
+    ActsScalar phimin = phi - phitol;
+    ActsScalar phimax = phi + phitol;
     if (phimin < get(eAveragePhi) - get(eHalfPhiSector) ||
         phimax > get(eAveragePhi) + get(eHalfPhiSector)) {
       return false;
     }
   }
   // We are within phi sector check box r quickly
-  double rmin = r + tol;
-  double rmax = r - tol;
+  ActsScalar rmin = r + tol;
+  ActsScalar rmax = r - tol;
   if (rmin > innerRmax() && rmax < outerRmin()) {
     return true;
   }
   // Finally we need to check the cone
   if (m_innerConeBounds != nullptr) {
-    double innerConeR = m_innerConeBounds->r(std::abs(z + get(eInnerOffsetZ)));
+    ActsScalar innerConeR =
+        m_innerConeBounds->r(std::abs(z + get(eInnerOffsetZ)));
     if (innerConeR > rmin) {
       return false;
     }
@@ -227,7 +229,8 @@ bool Acts::ConeVolumeBounds::inside(const Vector3& pos, double tol) const {
   }
   // And the outer cone
   if (m_outerConeBounds != nullptr) {
-    double outerConeR = m_outerConeBounds->r(std::abs(z + get(eOuterOffsetZ)));
+    ActsScalar outerConeR =
+        m_outerConeBounds->r(std::abs(z + get(eOuterOffsetZ)));
     if (outerConeR < rmax) {
       return false;
     }
@@ -237,12 +240,12 @@ bool Acts::ConeVolumeBounds::inside(const Vector3& pos, double tol) const {
   return true;
 }
 
-void Acts::ConeVolumeBounds::buildSurfaceBounds() {
+void ConeVolumeBounds::buildSurfaceBounds() {
   // Build inner cone or inner cylinder
   if (get(eInnerAlpha) > s_epsilon) {
     m_innerTanAlpha = std::tan(get(eInnerAlpha));
-    double innerZmin = get(eInnerOffsetZ) - get(eHalfLengthZ);
-    double innerZmax = get(eInnerOffsetZ) + get(eHalfLengthZ);
+    ActsScalar innerZmin = get(eInnerOffsetZ) - get(eHalfLengthZ);
+    ActsScalar innerZmax = get(eInnerOffsetZ) + get(eHalfLengthZ);
     m_innerRmin = std::abs(innerZmin) * m_innerTanAlpha;
     m_innerRmax = std::abs(innerZmax) * m_innerTanAlpha;
     m_innerConeBounds =
@@ -255,8 +258,8 @@ void Acts::ConeVolumeBounds::buildSurfaceBounds() {
 
   if (get(eOuterAlpha) > s_epsilon) {
     m_outerTanAlpha = std::tan(get(eOuterAlpha));
-    double outerZmin = get(eOuterOffsetZ) - get(eHalfLengthZ);
-    double outerZmax = get(eOuterOffsetZ) + get(eHalfLengthZ);
+    ActsScalar outerZmin = get(eOuterOffsetZ) - get(eHalfLengthZ);
+    ActsScalar outerZmax = get(eOuterOffsetZ) + get(eHalfLengthZ);
     m_outerRmin = std::abs(outerZmin) * m_outerTanAlpha;
     m_outerRmax = std::abs(outerZmax) * m_outerTanAlpha;
     m_outerConeBounds =
@@ -288,16 +291,54 @@ void Acts::ConeVolumeBounds::buildSurfaceBounds() {
   }
 }
 
-// ostream operator overload
-std::ostream& Acts::ConeVolumeBounds::toStream(std::ostream& sl) const {
-  return dumpT(sl);
+std::ostream& ConeVolumeBounds::toStream(std::ostream& os) const {
+  os << std::setiosflags(std::ios::fixed);
+  os << std::setprecision(5);
+  os << "Acts::ConeVolumeBounds : (innerAlpha, innerOffsetZ, outerAlpha,";
+  os << "  outerOffsetZ, halflenghZ, averagePhi, halfPhiSector) = ";
+  os << get(eInnerAlpha) << ", " << get(eInnerOffsetZ) << ", ";
+  os << get(eOuterAlpha) << ", " << get(eOuterOffsetZ) << ", ";
+  os << get(eHalfLengthZ) << ", " << get(eAveragePhi) << std::endl;
+  return os;
 }
 
-Acts::Volume::BoundingBox Acts::ConeVolumeBounds::boundingBox(
-    const Acts::Transform3* trf, const Vector3& envelope,
-    const Volume* entity) const {
+Volume::BoundingBox ConeVolumeBounds::boundingBox(const Transform3* trf,
+                                                  const Vector3& envelope,
+                                                  const Volume* entity) const {
   Vector3 vmin(-outerRmax(), -outerRmax(), -0.5 * get(eHalfLengthZ));
   Vector3 vmax(outerRmax(), outerRmax(), 0.5 * get(eHalfLengthZ));
   Volume::BoundingBox box(entity, vmin - envelope, vmax + envelope);
   return trf == nullptr ? box : box.transformed(*trf);
 }
+
+ActsScalar ConeVolumeBounds::innerRmin() const {
+  return m_innerRmin;
+}
+
+ActsScalar ConeVolumeBounds::innerRmax() const {
+  return m_innerRmax;
+}
+
+ActsScalar ConeVolumeBounds::innerTanAlpha() const {
+  return m_innerTanAlpha;
+}
+
+ActsScalar ConeVolumeBounds::outerRmin() const {
+  return m_outerRmin;
+}
+
+ActsScalar ConeVolumeBounds::outerRmax() const {
+  return m_outerRmax;
+}
+
+ActsScalar ConeVolumeBounds::outerTanAlpha() const {
+  return m_outerTanAlpha;
+}
+
+std::vector<ActsScalar> ConeVolumeBounds::values() const {
+  std::vector<ActsScalar> valvector;
+  valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
+  return valvector;
+}
+
+}  // namespace Acts
diff --git a/Core/src/Geometry/CuboidVolumeBounds.cpp b/Core/src/Geometry/CuboidVolumeBounds.cpp
index 7323d30b0b9..c3b570b1481 100644
--- a/Core/src/Geometry/CuboidVolumeBounds.cpp
+++ b/Core/src/Geometry/CuboidVolumeBounds.cpp
@@ -17,21 +17,23 @@
 #include <type_traits>
 #include <utility>
 
-Acts::CuboidVolumeBounds::CuboidVolumeBounds(double halex, double haley,
-                                             double halez)
+namespace Acts {
+
+CuboidVolumeBounds::CuboidVolumeBounds(ActsScalar halex, ActsScalar haley,
+                                       ActsScalar halez)
     : VolumeBounds(), m_values({halex, haley, halez}) {
   checkConsistency();
   buildSurfaceBounds();
 }
 
-Acts::CuboidVolumeBounds::CuboidVolumeBounds(const CuboidVolumeBounds& bobo)
+CuboidVolumeBounds::CuboidVolumeBounds(const CuboidVolumeBounds& bobo)
     : VolumeBounds(),
       m_values(bobo.m_values),
       m_xyBounds(bobo.m_xyBounds),
       m_yzBounds(bobo.m_yzBounds),
       m_zxBounds(bobo.m_zxBounds) {}
 
-Acts::CuboidVolumeBounds& Acts::CuboidVolumeBounds::operator=(
+CuboidVolumeBounds& CuboidVolumeBounds::operator=(
     const CuboidVolumeBounds& bobo) {
   if (this != &bobo) {
     m_values = bobo.m_values;
@@ -81,12 +83,17 @@ Acts::OrientedSurfaces Acts::CuboidVolumeBounds::orientedSurfaces(
   return oSurfaces;
 }
 
-std::ostream& Acts::CuboidVolumeBounds::toStream(std::ostream& sl) const {
-  return dumpT(sl);
+std::ostream& CuboidVolumeBounds::toStream(std::ostream& os) const {
+  os << std::setiosflags(std::ios::fixed);
+  os << std::setprecision(5);
+  os << "Acts::CuboidVolumeBounds: (halfLengthX, halfLengthY, halfLengthZ) = ";
+  os << "(" << get(eHalfLengthX) << ", " << get(eHalfLengthY) << ", "
+     << get(eHalfLengthZ) << ")";
+  return os;
 }
 
-Acts::Volume::BoundingBox Acts::CuboidVolumeBounds::boundingBox(
-    const Acts::Transform3* trf, const Vector3& envelope,
+Volume::BoundingBox CuboidVolumeBounds::boundingBox(
+    const Transform3* trf, const Vector3& envelope,
     const Volume* entity) const {
   Vector3 vmin(-get(eHalfLengthX), -get(eHalfLengthY), -get(eHalfLengthZ));
   Vector3 vmax(get(eHalfLengthX), get(eHalfLengthY), get(eHalfLengthZ));
@@ -95,7 +102,7 @@ Acts::Volume::BoundingBox Acts::CuboidVolumeBounds::boundingBox(
   return trf == nullptr ? box : box.transformed(*trf);
 }
 
-void Acts::CuboidVolumeBounds::buildSurfaceBounds() {
+void CuboidVolumeBounds::buildSurfaceBounds() {
   m_xyBounds = std::make_shared<const RectangleBounds>(get(eHalfLengthX),
                                                        get(eHalfLengthY));
   m_yzBounds = std::make_shared<const RectangleBounds>(get(eHalfLengthY),
@@ -104,7 +111,7 @@ void Acts::CuboidVolumeBounds::buildSurfaceBounds() {
                                                        get(eHalfLengthX));
 }
 
-double Acts::CuboidVolumeBounds::binningBorder(BinningValue bValue) const {
+ActsScalar CuboidVolumeBounds::binningBorder(BinningValue bValue) const {
   if (bValue <= binZ) {
     return m_values[bValue];
   }
@@ -114,3 +121,24 @@ double Acts::CuboidVolumeBounds::binningBorder(BinningValue bValue) const {
   }
   return 0.0;
 }
+
+bool CuboidVolumeBounds::inside(const Vector3& pos, ActsScalar tol) const {
+  return (std::abs(pos.x()) <= get(eHalfLengthX) + tol &&
+          std::abs(pos.y()) <= get(eHalfLengthY) + tol &&
+          std::abs(pos.z()) <= get(eHalfLengthZ) + tol);
+}
+
+std::vector<ActsScalar> CuboidVolumeBounds::values() const {
+  std::vector<ActsScalar> valvector;
+  valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
+  return valvector;
+}
+
+void CuboidVolumeBounds::checkConsistency() noexcept(false) {
+  if (get(eHalfLengthX) <= 0 || get(eHalfLengthY) <= 0 ||
+      get(eHalfLengthZ) <= 0.) {
+    throw std::invalid_argument(
+        "CuboidVolumeBounds: invalid input, zero or negative.");
+  }
+}
+}  // namespace Acts
diff --git a/Core/src/Geometry/CuboidVolumeBuilder.cpp b/Core/src/Geometry/CuboidVolumeBuilder.cpp
index 7c53f557ded..10df4c7a6a6 100644
--- a/Core/src/Geometry/CuboidVolumeBuilder.cpp
+++ b/Core/src/Geometry/CuboidVolumeBuilder.cpp
@@ -200,7 +200,7 @@ std::shared_ptr<Acts::TrackingVolume> Acts::CuboidVolumeBuilder::buildVolume(
 
 Acts::MutableTrackingVolumePtr Acts::CuboidVolumeBuilder::trackingVolume(
     const GeometryContext& gctx, Acts::TrackingVolumePtr /*gctx*/,
-    Acts::VolumeBoundsPtr /*bounds*/) const {
+    std::shared_ptr<const VolumeBounds> /*bounds*/) const {
   // Build volumes
   std::vector<std::shared_ptr<TrackingVolume>> volumes;
   volumes.reserve(m_cfg.volumeCfg.size());
diff --git a/Core/src/Geometry/CylinderVolumeBounds.cpp b/Core/src/Geometry/CylinderVolumeBounds.cpp
index ff60b6fac77..261b2446b56 100644
--- a/Core/src/Geometry/CylinderVolumeBounds.cpp
+++ b/Core/src/Geometry/CylinderVolumeBounds.cpp
@@ -23,10 +23,12 @@
 #include <type_traits>
 #include <utility>
 
-Acts::CylinderVolumeBounds::CylinderVolumeBounds(
-    const CylinderBounds& cBounds, double thickness) noexcept(false)
+namespace Acts {
+
+CylinderVolumeBounds::CylinderVolumeBounds(const CylinderBounds& cBounds,
+                                           ActsScalar thickness) noexcept(false)
     : VolumeBounds() {
-  double cR = cBounds.get(CylinderBounds::eR);
+  ActsScalar cR = cBounds.get(CylinderBounds::eR);
   if (thickness <= 0. || (cR - 0.5 * thickness) < 0.) {
     throw(std::invalid_argument(
         "CylinderVolumeBounds: invalid extrusion thickness."));
@@ -41,8 +43,8 @@ Acts::CylinderVolumeBounds::CylinderVolumeBounds(
   buildSurfaceBounds();
 }
 
-Acts::CylinderVolumeBounds::CylinderVolumeBounds(
-    const RadialBounds& rBounds, double thickness) noexcept(false)
+CylinderVolumeBounds::CylinderVolumeBounds(const RadialBounds& rBounds,
+                                           ActsScalar thickness) noexcept(false)
     : VolumeBounds() {
   if (thickness <= 0.) {
     throw(std::invalid_argument(
@@ -53,12 +55,12 @@ Acts::CylinderVolumeBounds::CylinderVolumeBounds(
   m_values[eHalfLengthZ] = 0.5 * thickness;
   m_values[eHalfPhiSector] = rBounds.get(RadialBounds::eHalfPhiSector);
   m_values[eAveragePhi] = rBounds.get(RadialBounds::eAveragePhi);
-  m_values[eBevelMinZ] = (double)0.;
-  m_values[eBevelMaxZ] = (double)0.;
+  m_values[eBevelMinZ] = (ActsScalar)0.;
+  m_values[eBevelMaxZ] = (ActsScalar)0.;
   buildSurfaceBounds();
 }
 
-Acts::OrientedSurfaces Acts::CylinderVolumeBounds::orientedSurfaces(
+std::vector<OrientedSurface> CylinderVolumeBounds::orientedSurfaces(
     const Transform3& transform) const {
   OrientedSurfaces oSurfaces;
   oSurfaces.reserve(6);
@@ -66,11 +68,11 @@ Acts::OrientedSurfaces Acts::CylinderVolumeBounds::orientedSurfaces(
   Translation3 vMinZ(0., 0., -get(eHalfLengthZ));
   Translation3 vMaxZ(0., 0., get(eHalfLengthZ));
   // Set up transform for beveled edges if they are defined
-  double bevelMinZ = get(eBevelMinZ);
-  double bevelMaxZ = get(eBevelMaxZ);
+  ActsScalar bevelMinZ = get(eBevelMinZ);
+  ActsScalar bevelMaxZ = get(eBevelMaxZ);
   Transform3 transMinZ, transMaxZ;
   if (bevelMinZ != 0.) {
-    double sy = 1 - 1 / std::cos(bevelMinZ);
+    ActsScalar sy = 1 - 1 / std::cos(bevelMinZ);
     transMinZ = transform * vMinZ *
                 Eigen::AngleAxisd(-bevelMinZ, Eigen::Vector3d(1., 0., 0.)) *
                 Eigen::Scaling(1., 1. + sy, 1.);
@@ -78,7 +80,7 @@ Acts::OrientedSurfaces Acts::CylinderVolumeBounds::orientedSurfaces(
     transMinZ = transform * vMinZ;
   }
   if (bevelMaxZ != 0.) {
-    double sy = 1 - 1 / std::cos(bevelMaxZ);
+    ActsScalar sy = 1 - 1 / std::cos(bevelMaxZ);
     transMaxZ = transform * vMaxZ *
                 Eigen::AngleAxisd(bevelMaxZ, Eigen::Vector3d(1., 0., 0.)) *
                 Eigen::Scaling(1., 1. + sy, 1.);
@@ -136,7 +138,7 @@ Acts::OrientedSurfaces Acts::CylinderVolumeBounds::orientedSurfaces(
   return oSurfaces;
 }
 
-void Acts::CylinderVolumeBounds::buildSurfaceBounds() {
+void CylinderVolumeBounds::buildSurfaceBounds() {
   if (get(eMinR) > s_epsilon) {
     m_innerCylinderBounds = std::make_shared<const CylinderBounds>(
         get(eMinR), get(eHalfLengthZ), get(eHalfPhiSector), get(eAveragePhi),
@@ -154,14 +156,21 @@ void Acts::CylinderVolumeBounds::buildSurfaceBounds() {
   }
 }
 
-std::ostream& Acts::CylinderVolumeBounds::toStream(std::ostream& sl) const {
-  return dumpT<std::ostream>(sl);
+std::ostream& CylinderVolumeBounds::toStream(std::ostream& os) const {
+  os << std::setiosflags(std::ios::fixed);
+  os << std::setprecision(5);
+  os << "CylinderVolumeBounds: (rMin, rMax, halfZ, halfPhi, "
+        "averagePhi, minBevelZ, maxBevelZ) = ";
+  os << get(eMinR) << ", " << get(eMaxR) << ", " << get(eHalfLengthZ) << ", "
+     << get(eHalfPhiSector) << ", " << get(eAveragePhi) << ", "
+     << get(eBevelMinZ) << ", " << get(eBevelMaxZ);
+  return os;
 }
 
-Acts::Volume::BoundingBox Acts::CylinderVolumeBounds::boundingBox(
+Volume::BoundingBox CylinderVolumeBounds::boundingBox(
     const Transform3* trf, const Vector3& envelope,
     const Volume* entity) const {
-  double xmax = 0, xmin = 0, ymax = 0, ymin = 0;
+  ActsScalar xmax = 0, xmin = 0, ymax = 0, ymin = 0;
   xmax = get(eMaxR);
 
   if (get(eHalfPhiSector) > M_PI / 2.) {
@@ -184,3 +193,66 @@ Acts::Volume::BoundingBox Acts::CylinderVolumeBounds::boundingBox(
   Volume::BoundingBox box{entity, vmin - envelope, vmax + envelope};
   return trf == nullptr ? box : box.transformed(*trf);
 }
+
+bool CylinderVolumeBounds::inside(const Vector3& pos, ActsScalar tol) const {
+  using VectorHelpers::perp;
+  using VectorHelpers::phi;
+  ActsScalar ros = perp(pos);
+  bool insidePhi = cos(phi(pos)) >= cos(get(eHalfPhiSector)) - tol;
+  bool insideR = insidePhi
+                     ? ((ros >= get(eMinR) - tol) && (ros <= get(eMaxR) + tol))
+                     : false;
+  bool insideZ =
+      insideR ? (std::abs(pos.z()) <= get(eHalfLengthZ) + tol) : false;
+  return (insideZ && insideR && insidePhi);
+}
+
+Vector3 CylinderVolumeBounds::binningOffset(BinningValue bValue)
+    const {  // the medium radius is taken for r-type binning
+  if (bValue == Acts::binR || bValue == Acts::binRPhi) {
+    return Vector3(0.5 * (get(eMinR) + get(eMaxR)), 0., 0.);
+  }
+  return VolumeBounds::binningOffset(bValue);
+}
+
+ActsScalar CylinderVolumeBounds::binningBorder(BinningValue bValue) const {
+  if (bValue == Acts::binR) {
+    return 0.5 * (get(eMaxR) - get(eMinR));
+  }
+  if (bValue == Acts::binZ) {
+    return get(eHalfLengthZ);
+  }
+  return VolumeBounds::binningBorder(bValue);
+}
+
+std::vector<ActsScalar> CylinderVolumeBounds::values() const {
+  std::vector<ActsScalar> valvector;
+  valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
+  return valvector;
+}
+
+void CylinderVolumeBounds::checkConsistency() noexcept(false) {
+  if (get(eMinR) < 0. || get(eMaxR) <= 0. || get(eMinR) >= get(eMaxR)) {
+    throw std::invalid_argument("CylinderVolumeBounds: invalid radial input.");
+  }
+  if (get(eHalfLengthZ) <= 0) {
+    throw std::invalid_argument(
+        "CylinderVolumeBounds: invalid longitudinal input.");
+  }
+  if (get(eHalfPhiSector) < 0. || get(eHalfPhiSector) > M_PI) {
+    throw std::invalid_argument(
+        "CylinderVolumeBounds: invalid phi sector setup.");
+  }
+  if (get(eAveragePhi) != detail::radian_sym(get(eAveragePhi))) {
+    throw std::invalid_argument(
+        "CylinderVolumeBounds: invalid phi positioning.");
+  }
+  if (get(eBevelMinZ) != detail::radian_sym(get(eBevelMinZ))) {
+    throw std::invalid_argument("CylinderBounds: invalid bevel at min Z.");
+  }
+  if (get(eBevelMaxZ) != detail::radian_sym(get(eBevelMaxZ))) {
+    throw std::invalid_argument("CylinderBounds: invalid bevel at max Z.");
+  }
+}
+
+}  // namespace Acts
diff --git a/Core/src/Geometry/CylinderVolumeBuilder.cpp b/Core/src/Geometry/CylinderVolumeBuilder.cpp
index a1b532f8b53..66f87529797 100644
--- a/Core/src/Geometry/CylinderVolumeBuilder.cpp
+++ b/Core/src/Geometry/CylinderVolumeBuilder.cpp
@@ -58,7 +58,7 @@ void Acts::CylinderVolumeBuilder::setLogger(
 std::shared_ptr<Acts::TrackingVolume>
 Acts::CylinderVolumeBuilder::trackingVolume(
     const GeometryContext& gctx, TrackingVolumePtr existingVolume,
-    VolumeBoundsPtr externalBounds) const {
+    std::shared_ptr<const VolumeBounds> externalBounds) const {
   ACTS_DEBUG("Configured to build volume : " << m_cfg.volumeName);
   if (existingVolume) {
     ACTS_DEBUG("- will wrap/enclose : " << existingVolume->volumeName());
diff --git a/Core/src/Geometry/CylinderVolumeHelper.cpp b/Core/src/Geometry/CylinderVolumeHelper.cpp
index a67d0991384..5b969ec0d92 100644
--- a/Core/src/Geometry/CylinderVolumeHelper.cpp
+++ b/Core/src/Geometry/CylinderVolumeHelper.cpp
@@ -202,8 +202,8 @@ Acts::CylinderVolumeHelper::createTrackingVolume(
   const Transform3 transform = Transform3(Translation3(0., 0., zPosition));
   // call to the creation method with Bounds & Translation3
   return createTrackingVolume(gctx, layers, volumeMaterial,
-                              VolumeBoundsPtr(cBounds), mtvVector, transform,
-                              volumeName, bType);
+                              std::shared_ptr<const VolumeBounds>(cBounds),
+                              mtvVector, transform, volumeName, bType);
 }
 
 std::shared_ptr<Acts::TrackingVolume>
diff --git a/Core/src/Geometry/Extent.cpp b/Core/src/Geometry/Extent.cpp
index aaf9759228d..a38b7769e62 100644
--- a/Core/src/Geometry/Extent.cpp
+++ b/Core/src/Geometry/Extent.cpp
@@ -96,7 +96,7 @@ void Acts::Extent::set(BinningValue bValue, ActsScalar min, ActsScalar max) {
   if (bValue == binR && minval < 0.) {
     minval = 0.;
   }
-  m_range[bValue] = Range1D{minval, max};
+  m_range[bValue] = Range1D<ActsScalar>{minval, max};
   m_constrains.set(bValue);
 }
 
diff --git a/Core/src/Geometry/KDTreeTrackingGeometryBuilder.cpp b/Core/src/Geometry/KDTreeTrackingGeometryBuilder.cpp
index 008bc327dbf..ed2969c1e56 100644
--- a/Core/src/Geometry/KDTreeTrackingGeometryBuilder.cpp
+++ b/Core/src/Geometry/KDTreeTrackingGeometryBuilder.cpp
@@ -26,7 +26,6 @@
 #include "Acts/Surfaces/SurfaceArray.hpp"
 #include "Acts/Surfaces/SurfaceBounds.hpp"
 #include "Acts/Utilities/BinningType.hpp"
-#include "Acts/Utilities/Range1D.hpp"
 #include "Acts/Utilities/RangeXD.hpp"
 
 #include <cstddef>
diff --git a/Core/src/Geometry/Polyhedron.cpp b/Core/src/Geometry/Polyhedron.cpp
index b6970ecc7fc..c2af9153349 100644
--- a/Core/src/Geometry/Polyhedron.cpp
+++ b/Core/src/Geometry/Polyhedron.cpp
@@ -10,7 +10,6 @@
 
 #include "Acts/Surfaces/detail/VerticesHelper.hpp"
 #include "Acts/Utilities/BinningType.hpp"
-#include "Acts/Utilities/Range1D.hpp"
 
 #include <algorithm>
 #include <cmath>
diff --git a/Core/src/Geometry/TrackingVolume.cpp b/Core/src/Geometry/TrackingVolume.cpp
index 2531b688689..89fbb7eda26 100644
--- a/Core/src/Geometry/TrackingVolume.cpp
+++ b/Core/src/Geometry/TrackingVolume.cpp
@@ -37,7 +37,7 @@ class ISurfaceMaterial;
 }  // namespace Acts
 
 Acts::TrackingVolume::TrackingVolume(
-    const Transform3& transform, VolumeBoundsPtr volbounds,
+    const Transform3& transform, std::shared_ptr<const VolumeBounds> volbounds,
     const std::shared_ptr<const TrackingVolumeArray>& containedVolumeArray,
     const std::string& volumeName)
     : Volume(transform, std::move(volbounds)),
@@ -52,7 +52,8 @@ Acts::TrackingVolume::TrackingVolume(
 
 // constructor for arguments
 Acts::TrackingVolume::TrackingVolume(
-    const Transform3& transform, VolumeBoundsPtr volumeBounds,
+    const Transform3& transform,
+    std::shared_ptr<const VolumeBounds> volumeBounds,
     std::shared_ptr<const IVolumeMaterial> volumeMaterial,
     std::unique_ptr<const LayerArray> staticLayerArray,
     std::shared_ptr<const TrackingVolumeArray> containedVolumeArray,
diff --git a/Core/src/Geometry/TrapezoidVolumeBounds.cpp b/Core/src/Geometry/TrapezoidVolumeBounds.cpp
index d70a5a126b0..49e4eb5c163 100644
--- a/Core/src/Geometry/TrapezoidVolumeBounds.cpp
+++ b/Core/src/Geometry/TrapezoidVolumeBounds.cpp
@@ -20,9 +20,11 @@
 #include <cstddef>
 #include <utility>
 
-Acts::TrapezoidVolumeBounds::TrapezoidVolumeBounds(double minhalex,
-                                                   double maxhalex,
-                                                   double haley, double halez)
+namespace Acts {
+
+TrapezoidVolumeBounds::TrapezoidVolumeBounds(ActsScalar minhalex,
+                                             ActsScalar maxhalex,
+                                             ActsScalar haley, ActsScalar halez)
     : VolumeBounds() {
   m_values[eHalfLengthXnegY] = minhalex;
   m_values[eHalfLengthXposY] = maxhalex;
@@ -34,9 +36,9 @@ Acts::TrapezoidVolumeBounds::TrapezoidVolumeBounds(double minhalex,
   buildSurfaceBounds();
 }
 
-Acts::TrapezoidVolumeBounds::TrapezoidVolumeBounds(double minhalex,
-                                                   double haley, double halez,
-                                                   double alpha, double beta)
+TrapezoidVolumeBounds::TrapezoidVolumeBounds(ActsScalar minhalex,
+                                             ActsScalar haley, ActsScalar halez,
+                                             ActsScalar alpha, ActsScalar beta)
     : VolumeBounds() {
   m_values[eHalfLengthXnegY] = minhalex;
   m_values[eHalfLengthY] = haley;
@@ -44,14 +46,15 @@ Acts::TrapezoidVolumeBounds::TrapezoidVolumeBounds(double minhalex,
   m_values[eAlpha] = alpha;
   m_values[eBeta] = beta;
   // now calculate the remaining max half X
-  double gamma = (alpha > beta) ? (alpha - 0.5 * M_PI) : (beta - 0.5 * M_PI);
+  ActsScalar gamma =
+      (alpha > beta) ? (alpha - 0.5 * M_PI) : (beta - 0.5 * M_PI);
   m_values[eHalfLengthXposY] = minhalex + (2. * haley) * tan(gamma);
 
   checkConsistency();
   buildSurfaceBounds();
 }
 
-Acts::OrientedSurfaces Acts::TrapezoidVolumeBounds::orientedSurfaces(
+std::vector<OrientedSurface> TrapezoidVolumeBounds::orientedSurfaces(
     const Transform3& transform) const {
   OrientedSurfaces oSurfaces;
   oSurfaces.reserve(6);
@@ -73,9 +76,9 @@ Acts::OrientedSurfaces Acts::TrapezoidVolumeBounds::orientedSurfaces(
   sf = Surface::makeShared<PlaneSurface>(pzTransform, m_faceXYTrapezoidBounds);
   oSurfaces.push_back(OrientedSurface(std::move(sf), Direction::Negative));
 
-  double poshOffset = get(eHalfLengthY) / std::tan(get(eAlpha));
-  double neghOffset = get(eHalfLengthY) / std::tan(get(eBeta));
-  double topShift = poshOffset + neghOffset;
+  ActsScalar poshOffset = get(eHalfLengthY) / std::tan(get(eAlpha));
+  ActsScalar neghOffset = get(eHalfLengthY) / std::tan(get(eBeta));
+  ActsScalar topShift = poshOffset + neghOffset;
 
   // Face surfaces yz
   // (3) - At point B, attached to beta opening angle
@@ -113,7 +116,7 @@ Acts::OrientedSurfaces Acts::TrapezoidVolumeBounds::orientedSurfaces(
   return oSurfaces;
 }
 
-void Acts::TrapezoidVolumeBounds::buildSurfaceBounds() {
+void TrapezoidVolumeBounds::buildSurfaceBounds() {
   m_faceXYTrapezoidBounds = std::make_shared<const TrapezoidBounds>(
       get(eHalfLengthXnegY), get(eHalfLengthXposY), get(eHalfLengthY));
 
@@ -130,7 +133,7 @@ void Acts::TrapezoidVolumeBounds::buildSurfaceBounds() {
       get(eHalfLengthZ), get(eHalfLengthXposY));
 }
 
-bool Acts::TrapezoidVolumeBounds::inside(const Vector3& pos, double tol) const {
+bool TrapezoidVolumeBounds::inside(const Vector3& pos, ActsScalar tol) const {
   if (std::abs(pos.z()) > get(eHalfLengthZ) + tol) {
     return false;
   }
@@ -143,17 +146,24 @@ bool Acts::TrapezoidVolumeBounds::inside(const Vector3& pos, double tol) const {
   return inside;
 }
 
-std::ostream& Acts::TrapezoidVolumeBounds::toStream(std::ostream& sl) const {
-  return dumpT<std::ostream>(sl);
+std::ostream& TrapezoidVolumeBounds::toStream(std::ostream& os) const {
+  os << std::setiosflags(std::ios::fixed);
+  os << std::setprecision(5);
+  os << "TrapezoidVolumeBounds: (minhalfX, halfY, halfZ, alpha, beta) "
+        "= ";
+  os << "(" << get(eHalfLengthXnegY) << ", " << get(eHalfLengthXposY) << ", "
+     << get(eHalfLengthY) << ", " << get(eHalfLengthZ);
+  os << ", " << get(eAlpha) << ", " << get(eBeta) << ")";
+  return os;
 }
 
-Acts::Volume::BoundingBox Acts::TrapezoidVolumeBounds::boundingBox(
-    const Acts::Transform3* trf, const Vector3& envelope,
+Volume::BoundingBox TrapezoidVolumeBounds::boundingBox(
+    const Transform3* trf, const Vector3& envelope,
     const Volume* entity) const {
-  double minx = get(eHalfLengthXnegY);
-  double maxx = get(eHalfLengthXposY);
-  double haley = get(eHalfLengthY);
-  double halez = get(eHalfLengthZ);
+  ActsScalar minx = get(eHalfLengthXnegY);
+  ActsScalar maxx = get(eHalfLengthXposY);
+  ActsScalar haley = get(eHalfLengthY);
+  ActsScalar halez = get(eHalfLengthZ);
 
   std::array<Vector3, 8> vertices = {{{-minx, -haley, -halez},
                                       {+minx, -haley, -halez},
@@ -180,3 +190,24 @@ Acts::Volume::BoundingBox Acts::TrapezoidVolumeBounds::boundingBox(
 
   return {entity, vmin - envelope, vmax + envelope};
 }
+
+std::vector<ActsScalar> TrapezoidVolumeBounds::values() const {
+  std::vector<ActsScalar> valvector;
+  valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
+  return valvector;
+}
+
+void TrapezoidVolumeBounds::checkConsistency() noexcept(false) {
+  if (get(eHalfLengthXnegY) < 0. || get(eHalfLengthXposY) < 0.) {
+    throw std::invalid_argument(
+        "TrapezoidVolumeBounds: invalid trapezoid parameters in x.");
+  }
+  if (get(eHalfLengthY) <= 0.) {
+    throw std::invalid_argument("TrapezoidVolumeBounds: invalid y extrusion.");
+  }
+  if (get(eHalfLengthZ) <= 0.) {
+    throw std::invalid_argument("TrapezoidVolumeBounds: invalid z extrusion.");
+  }
+}
+
+}  // namespace Acts
diff --git a/Core/src/Geometry/Volume.cpp b/Core/src/Geometry/Volume.cpp
index 98599825dd9..d87da05a010 100644
--- a/Core/src/Geometry/Volume.cpp
+++ b/Core/src/Geometry/Volume.cpp
@@ -16,8 +16,10 @@
 
 using namespace Acts::UnitLiterals;
 
-Acts::Volume::Volume(const Transform3& transform,
-                     std::shared_ptr<const VolumeBounds> volbounds)
+namespace Acts {
+
+Volume::Volume(const Transform3& transform,
+               std::shared_ptr<const VolumeBounds> volbounds)
     : GeometryObject(),
       m_transform(transform),
       m_itransform(m_transform.inverse()),
@@ -26,7 +28,7 @@ Acts::Volume::Volume(const Transform3& transform,
       m_orientedBoundingBox(m_volumeBounds->boundingBox(
           nullptr, {0.05_mm, 0.05_mm, 0.05_mm}, this)) {}
 
-Acts::Volume::Volume(const Volume& vol, const Transform3& shift)
+Volume::Volume(const Volume& vol, const Transform3& shift)
     : GeometryObject(),
       m_transform(shift * vol.m_transform),
       m_itransform(m_transform.inverse()),
@@ -35,11 +37,11 @@ Acts::Volume::Volume(const Volume& vol, const Transform3& shift)
       m_orientedBoundingBox(m_volumeBounds->boundingBox(
           nullptr, {0.05_mm, 0.05_mm, 0.05_mm}, this)) {}
 
-Acts::Vector3 Acts::Volume::binningPosition(const GeometryContext& /*gctx*/,
-                                            Acts::BinningValue bValue) const {
+Vector3 Volume::binningPosition(const GeometryContext& /*gctx*/,
+                                BinningValue bValue) const {
   // for most of the binning types it is actually the center,
   // just for R-binning types the
-  if (bValue == Acts::binR || bValue == Acts::binRPhi) {
+  if (bValue == binR || bValue == binRPhi) {
     // the binning Position for R-type may have an offset
     return (center() + m_volumeBounds->binningOffset(bValue));
   }
@@ -48,7 +50,7 @@ Acts::Vector3 Acts::Volume::binningPosition(const GeometryContext& /*gctx*/,
 }
 
 // assignment operator
-Acts::Volume& Acts::Volume::operator=(const Acts::Volume& vol) {
+Volume& Volume::operator=(const Volume& vol) {
   if (this != &vol) {
     m_transform = vol.m_transform;
     m_center = vol.m_center;
@@ -57,21 +59,44 @@ Acts::Volume& Acts::Volume::operator=(const Acts::Volume& vol) {
   return *this;
 }
 
-bool Acts::Volume::inside(const Acts::Vector3& gpos, double tol) const {
-  Acts::Vector3 posInVolFrame((transform().inverse()) * gpos);
+bool Volume::inside(const Vector3& gpos, ActsScalar tol) const {
+  Vector3 posInVolFrame((transform().inverse()) * gpos);
   return (volumeBounds()).inside(posInVolFrame, tol);
 }
 
-std::ostream& Acts::operator<<(std::ostream& sl, const Acts::Volume& vol) {
+std::ostream& operator<<(std::ostream& sl, const Volume& vol) {
   sl << "Volume with " << vol.volumeBounds() << std::endl;
   return sl;
 }
 
-Acts::Volume::BoundingBox Acts::Volume::boundingBox(
-    const Vector3& envelope) const {
+Volume::BoundingBox Volume::boundingBox(const Vector3& envelope) const {
   return m_volumeBounds->boundingBox(&m_transform, envelope, this);
 }
 
-const Acts::Volume::BoundingBox& Acts::Volume::orientedBoundingBox() const {
+const Volume::BoundingBox& Volume::orientedBoundingBox() const {
   return m_orientedBoundingBox;
 }
+
+void Volume::assignVolumeBounds(std::shared_ptr<const VolumeBounds> volbounds) {
+  m_volumeBounds = std::move(volbounds);
+  m_orientedBoundingBox =
+      m_volumeBounds->boundingBox(nullptr, {0.05_mm, 0.05_mm, 0.05_mm}, this);
+}
+
+const Transform3& Volume::transform() const {
+  return m_transform;
+}
+
+const Transform3& Volume::itransform() const {
+  return m_itransform;
+}
+
+const Vector3& Volume::center() const {
+  return m_center;
+}
+
+const VolumeBounds& Volume::volumeBounds() const {
+  return *m_volumeBounds;
+}
+
+}  // namespace Acts
diff --git a/Examples/Algorithms/TrackFinding/src/SeedingAlgorithm.cpp b/Examples/Algorithms/TrackFinding/src/SeedingAlgorithm.cpp
index 3bc1e6e329a..55100f4e6a0 100644
--- a/Examples/Algorithms/TrackFinding/src/SeedingAlgorithm.cpp
+++ b/Examples/Algorithms/TrackFinding/src/SeedingAlgorithm.cpp
@@ -19,7 +19,6 @@
 #include "Acts/Utilities/Grid.hpp"
 #include "Acts/Utilities/GridBinFinder.hpp"
 #include "Acts/Utilities/Helpers.hpp"
-#include "Acts/Utilities/Range1D.hpp"
 #include "ActsExamples/EventData/SimSeed.hpp"
 
 #include <cmath>
diff --git a/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp b/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp
index 5e236f9fa40..d5e97411c34 100644
--- a/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp
+++ b/Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp
@@ -133,17 +133,22 @@ auto ActsExamples::AdaptiveMultiVertexFinderAlgorithm::makeVertexFinder() const
   finderConfig.tracksMaxZinterval = 1. * Acts::UnitConstants::mm;
   finderConfig.maxIterations = 200;
   finderConfig.useTime = m_cfg.useTime;
+  // 5 corresponds to a p-value of ~0.92 using `chi2(x=5,ndf=2)`
+  finderConfig.tracksMaxSignificance = 5;
+  // This should be used consistently with and without time
+  finderConfig.doFullSplitting = false;
+  // 3 corresponds to a p-value of ~0.92 using `chi2(x=3,ndf=1)`
+  finderConfig.maxMergeVertexSignificance = 3;
   if (m_cfg.useTime) {
     // When using time, we have an extra contribution to the chi2 by the time
     // coordinate. We thus need to increase tracksMaxSignificance (i.e., the
     // maximum chi2 that a track can have to be associated with a vertex).
-    finderConfig.tracksMaxSignificance = 7.5;
-    finderConfig.doFullSplitting = true;
-    // Reset the maximum significance that two vertices can have before they
-    // are considered as merged. The default value 3 is tuned for comparing
-    // the vertices' z-coordinates. Since we consider 4 dimensions here, we
-    // need to multiply the value by 4 and thus we set it to 3 * 4 = 12.
-    finderConfig.maxMergeVertexSignificance = 12.;
+    // Using the same p-value for 3 dof instead of 2.
+    // 6.7 corresponds to a p-value of ~0.92 using `chi2(x=6.7,ndf=3)`
+    finderConfig.tracksMaxSignificance = 6.7;
+    // Using the same p-value for 2 dof instead of 1.
+    // 5 corresponds to a p-value of ~0.92 using `chi2(x=5,ndf=2)`
+    finderConfig.maxMergeVertexSignificance = 5;
   }
   finderConfig.extractParameters
       .template connect<&Acts::InputTrack::extractParameters>();
diff --git a/Examples/Detectors/TelescopeDetector/src/BuildTelescopeDetector.cpp b/Examples/Detectors/TelescopeDetector/src/BuildTelescopeDetector.cpp
index 11b2244ba10..0f8d2732ef1 100644
--- a/Examples/Detectors/TelescopeDetector/src/BuildTelescopeDetector.cpp
+++ b/Examples/Detectors/TelescopeDetector/src/BuildTelescopeDetector.cpp
@@ -128,7 +128,7 @@ ActsExamples::Telescope::buildDetector(
   // The volume bounds is set to be a bit larger than either cubic with planes
   // or cylinder with discs
   auto length = positions.back() - positions.front();
-  Acts::VolumeBoundsPtr boundsVol = nullptr;
+  std::shared_ptr<const Acts::VolumeBounds> boundsVol = nullptr;
   if (surfaceType == TelescopeSurfaceType::Plane) {
     boundsVol = std::make_shared<const Acts::CuboidVolumeBounds>(
         bounds[0] + 5._mm, bounds[1] + 5._mm, length + 10._mm);
diff --git a/Examples/Scripts/Python/full_chain_odd.py b/Examples/Scripts/Python/full_chain_odd.py
index 9a8040b37b7..185a6b0f808 100755
--- a/Examples/Scripts/Python/full_chain_odd.py
+++ b/Examples/Scripts/Python/full_chain_odd.py
@@ -1,5 +1,11 @@
 #!/usr/bin/env python3
-import os, argparse, pathlib, acts, acts.examples
+
+import os
+import argparse
+import pathlib
+
+import acts
+import acts.examples
 from acts.examples.simulation import (
     addParticleGun,
     MomentumConfig,
@@ -29,12 +35,20 @@
 )
 from common import getOpenDataDetectorDirectory
 from acts.examples.odd import getOpenDataDetector
-import acts.examples.edm4hep
 
-parser = argparse.ArgumentParser(description="Full chain with the OpenDataDetector")
 
-parser.add_argument("--events", "-n", help="Number of events", type=int, default=100)
+u = acts.UnitConstants
+
 
+parser = argparse.ArgumentParser(description="Full chain with the OpenDataDetector")
+parser.add_argument(
+    "--output",
+    "-o",
+    help="Output directory",
+    type=pathlib.Path,
+    default=pathlib.Path.cwd() / "odd_output",
+)
+parser.add_argument("--events", "-n", help="Number of events", type=int, default=100)
 parser.add_argument("--skip", "-s", help="Number of events", type=int, default=0)
 parser.add_argument("--edm4hep", help="Use edm4hep inputs", type=pathlib.Path)
 parser.add_argument(
@@ -58,13 +72,12 @@
 
 args = vars(parser.parse_args())
 
+outputDir = args["output"]
 ttbar = args["ttbar"]
 g4_simulation = args["geant4"]
 ambiguity_MLSolver = args["MLSolver"]
 seedFilter_ML = args["MLSeedFilter"]
-u = acts.UnitConstants
 geoDir = getOpenDataDetectorDirectory()
-outputDir = pathlib.Path.cwd() / "odd_output"
 # acts.examples.dump_args_calls(locals())  # show python binding calls
 
 oddMaterialMap = geoDir / "data/odd-material-maps.root"
@@ -86,6 +99,8 @@
 )
 
 if args["edm4hep"]:
+    import acts.examples.edm4hep
+
     edm4hepReader = acts.examples.edm4hep.EDM4hepReader(
         inputPath=str(args["edm4hep"]),
         inputSimHits=[
diff --git a/Plugins/DD4hep/src/ConvertDD4hepDetector.cpp b/Plugins/DD4hep/src/ConvertDD4hepDetector.cpp
index 5fcd56d1355..c4432b88161 100644
--- a/Plugins/DD4hep/src/ConvertDD4hepDetector.cpp
+++ b/Plugins/DD4hep/src/ConvertDD4hepDetector.cpp
@@ -128,14 +128,14 @@ std::unique_ptr<const TrackingGeometry> convertDD4hepDetector(
 
   std::vector<std::function<std::shared_ptr<TrackingVolume>(
       const GeometryContext&, const TrackingVolumePtr&,
-      const VolumeBoundsPtr&)>>
+      const std::shared_ptr<const VolumeBounds>&)>>
       volumeFactories;
 
   for (const auto& vb : volumeBuilders) {
     volumeFactories.push_back(
         [vb](const GeometryContext& vgctx,
              const std::shared_ptr<const TrackingVolume>& inner,
-             const VolumeBoundsPtr&) {
+             const std::shared_ptr<const VolumeBounds>&) {
           return vb->trackingVolume(vgctx, inner);
         });
   }
diff --git a/Plugins/DD4hep/src/DD4hepLayerBuilder.cpp b/Plugins/DD4hep/src/DD4hepLayerBuilder.cpp
index d7a616ed40a..607c106ca61 100644
--- a/Plugins/DD4hep/src/DD4hepLayerBuilder.cpp
+++ b/Plugins/DD4hep/src/DD4hepLayerBuilder.cpp
@@ -27,7 +27,6 @@
 #include "Acts/Surfaces/SurfaceArray.hpp"
 #include "Acts/Utilities/Helpers.hpp"
 #include "Acts/Utilities/Logger.hpp"
-#include "Acts/Utilities/Range1D.hpp"
 
 #include <algorithm>
 #include <array>
diff --git a/Plugins/ExaTrkX/src/buildEdges.cpp b/Plugins/ExaTrkX/src/buildEdges.cpp
index 1b133291d7f..79f5d110f46 100644
--- a/Plugins/ExaTrkX/src/buildEdges.cpp
+++ b/Plugins/ExaTrkX/src/buildEdges.cpp
@@ -191,7 +191,7 @@ struct Span {
   const_iterator cbegin() const { return ptr; }
   const_iterator cend() const { return ptr + S; }
 
-  auto operator[](std::size_t i) const { return ptr[i]; }
+  auto &operator[](std::size_t i) const { return ptr[i]; }
 };
 
 template <std::size_t Dim>
@@ -235,9 +235,9 @@ struct BuildEdgesKDTree {
     for (int iself = 0; iself < embedFeatures.size(0); ++iself) {
       const Span<float, Dim> self{dataPtr + iself * Dim};
 
-      Acts::RangeXD<Dim, float, Span> range;
+      Acts::RangeXD<Dim, float> range;
       for (auto j = 0ul; j < Dim; ++j) {
-        range[j] = Acts::Range1D(self[j] - rVal, self[j] + rVal);
+        range[j] = Acts::Range1D<float>(self[j] - rVal, self[j] + rVal);
       }
 
       tree.rangeSearchMapDiscard(
diff --git a/Plugins/Json/include/Acts/Plugins/Json/AlgebraJsonConverter.hpp b/Plugins/Json/include/Acts/Plugins/Json/AlgebraJsonConverter.hpp
index b54cef52eca..b93a8703795 100644
--- a/Plugins/Json/include/Acts/Plugins/Json/AlgebraJsonConverter.hpp
+++ b/Plugins/Json/include/Acts/Plugins/Json/AlgebraJsonConverter.hpp
@@ -28,7 +28,7 @@ struct Options {
   /// Write the identity transform explicitly
   bool writeIdentity = false;
   /// Apply a transpose to flip column/row order
-  bool transpose = true;
+  bool transpose = false;
 };
 
 /// @brief The Transform converter to json
diff --git a/Plugins/Json/include/Acts/Plugins/Json/UtilitiesJsonConverter.hpp b/Plugins/Json/include/Acts/Plugins/Json/UtilitiesJsonConverter.hpp
index 65cc312d6a2..8f6e546b936 100644
--- a/Plugins/Json/include/Acts/Plugins/Json/UtilitiesJsonConverter.hpp
+++ b/Plugins/Json/include/Acts/Plugins/Json/UtilitiesJsonConverter.hpp
@@ -11,7 +11,7 @@
 #include "Acts/Plugins/Json/ActsJson.hpp"
 #include "Acts/Utilities/BinUtility.hpp"
 #include "Acts/Utilities/BinningData.hpp"
-#include "Acts/Utilities/Range1D.hpp"
+#include "Acts/Utilities/RangeXD.hpp"
 
 #include <nlohmann/json.hpp>
 
@@ -20,8 +20,6 @@
 
 namespace Acts {
 class BinningData;
-template <typename Type>
-class Range1D;
 
 void to_json(nlohmann::json& j, const BinningData& bd);
 
diff --git a/Plugins/Json/src/ExtentJsonConverter.cpp b/Plugins/Json/src/ExtentJsonConverter.cpp
index a0ab9322cf8..8f644ce765e 100644
--- a/Plugins/Json/src/ExtentJsonConverter.cpp
+++ b/Plugins/Json/src/ExtentJsonConverter.cpp
@@ -12,7 +12,6 @@
 #include "Acts/Plugins/Json/UtilitiesJsonConverter.hpp"
 #include "Acts/Utilities/BinningType.hpp"
 #include "Acts/Utilities/Enumerate.hpp"
-#include "Acts/Utilities/Range1D.hpp"
 #include "Acts/Utilities/RangeXD.hpp"
 
 #include <array>
diff --git a/Tests/UnitTests/Alignment/Kernel/AlignmentTests.cpp b/Tests/UnitTests/Alignment/Kernel/AlignmentTests.cpp
index 3acf3133c3c..20e2f8c676a 100644
--- a/Tests/UnitTests/Alignment/Kernel/AlignmentTests.cpp
+++ b/Tests/UnitTests/Alignment/Kernel/AlignmentTests.cpp
@@ -143,9 +143,10 @@ struct TelescopeDetector {
     // The volume transform
     Translation3 transVol(0, 0, 0);
     Transform3 trafoVol(rotation * transVol);
-    VolumeBoundsPtr boundsVol = std::make_shared<const CuboidVolumeBounds>(
-        rBounds->halfLengthX() + 10._mm, rBounds->halfLengthY() + 10._mm,
-        length + 10._mm);
+    std::shared_ptr<const VolumeBounds> boundsVol =
+        std::make_shared<const CuboidVolumeBounds>(
+            rBounds->halfLengthX() + 10._mm, rBounds->halfLengthY() + 10._mm,
+            length + 10._mm);
 
     LayerArrayCreator::Config lacConfig;
     LayerArrayCreator layArrCreator(
diff --git a/Tests/UnitTests/Core/Detector/DetectorVolumeTests.cpp b/Tests/UnitTests/Core/Detector/DetectorVolumeTests.cpp
index 91b861a1f43..c93ebd345e8 100644
--- a/Tests/UnitTests/Core/Detector/DetectorVolumeTests.cpp
+++ b/Tests/UnitTests/Core/Detector/DetectorVolumeTests.cpp
@@ -25,6 +25,8 @@
 #include "Acts/Navigation/SurfaceCandidatesUpdaters.hpp"
 #include "Acts/Surfaces/CylinderBounds.hpp"
 #include "Acts/Surfaces/CylinderSurface.hpp"
+#include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Surfaces/RectangleBounds.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
 #include "Acts/Utilities/BinningType.hpp"
@@ -53,6 +55,75 @@ Acts::GeometryContext tContext;
 
 BOOST_AUTO_TEST_SUITE(Detector)
 
+BOOST_AUTO_TEST_CASE(SurfaceVolumeContainment) {
+  // Create a surface that is placed way off
+  auto surfaceOutOfBounds = Acts::Surface::makeShared<Acts::CylinderSurface>(
+      Acts::Transform3::Identity() * Acts::Translation3(-1000, 0., 0.),
+      std::make_shared<Acts::CylinderBounds>(1, 1));
+
+  auto vBounds = Acts::CuboidVolumeBounds(10.0, 10.0, 10.0);
+  auto portalGenerator =
+      Acts::Experimental::defaultPortalAndSubPortalGenerator();
+  BOOST_CHECK_THROW(
+      Acts::Experimental::DetectorVolumeFactory::construct(
+          portalGenerator, Acts::GeometryContext(),
+          "CubeWithOutofBoundsSurface", Acts::Transform3::Identity(),
+          std::make_shared<Acts::CuboidVolumeBounds>(vBounds),
+          {surfaceOutOfBounds}, {}, Acts::Experimental::tryAllSubVolumes(),
+          Acts::Experimental::tryAllPortalsAndSurfaces(), 1000),
+      std::invalid_argument);
+
+  // Create a surface that is too big
+  auto surfaceTooBig = Acts::Surface::makeShared<Acts::PlaneSurface>(
+      Acts::Transform3::Identity() * Acts::Translation3(0, 0., 0.),
+      std::make_shared<Acts::RectangleBounds>(1, 100));
+
+  BOOST_CHECK_THROW(
+      Acts::Experimental::DetectorVolumeFactory::construct(
+          portalGenerator, Acts::GeometryContext(), "CubeWithSurfaceTooBig",
+          Acts::Transform3::Identity(),
+          std::make_shared<Acts::CuboidVolumeBounds>(vBounds), {surfaceTooBig},
+          {}, Acts::Experimental::tryAllSubVolumes(),
+          Acts::Experimental::tryAllPortalsAndSurfaces(), 1000),
+      std::invalid_argument);
+
+  // Envelope a bigger volume into a smaller one
+  auto bigVolume = Acts::Experimental::DetectorVolumeFactory::construct(
+      portalGenerator, Acts::GeometryContext(), "BigCube",
+      Acts::Transform3::Identity(),
+      std::make_shared<Acts::CuboidVolumeBounds>(vBounds), {}, {},
+      Acts::Experimental::tryAllSubVolumes(),
+      Acts::Experimental::tryAllPortalsAndSurfaces(), 1000);
+
+  auto smallBounds = Acts::CuboidVolumeBounds(1.0, 1.0, 1.0);
+  BOOST_CHECK_THROW(
+      Acts::Experimental::DetectorVolumeFactory::construct(
+          portalGenerator, Acts::GeometryContext(),
+          "SmallCubeWithBigCubeInside", Acts::Transform3::Identity(),
+          std::make_shared<Acts::CuboidVolumeBounds>(smallBounds), {},
+          {bigVolume}, Acts::Experimental::tryAllSubVolumes(),
+          Acts::Experimental::tryAllPortalsAndSurfaces(), 1000),
+      std::invalid_argument);
+
+  // Envelope a misaligned subvolume
+  auto smallVolumeMisaligned =
+      Acts::Experimental::DetectorVolumeFactory::construct(
+          portalGenerator, Acts::GeometryContext(), "SmallCubeMisaligned",
+          Acts::Transform3::Identity() * Acts::Translation3(9.5, 0., 0.),
+          std::make_shared<Acts::CuboidVolumeBounds>(smallBounds), {}, {},
+          Acts::Experimental::tryAllSubVolumes(),
+          Acts::Experimental::tryAllPortalsAndSurfaces(), 1000);
+
+  BOOST_CHECK_THROW(
+      Acts::Experimental::DetectorVolumeFactory::construct(
+          portalGenerator, Acts::GeometryContext(), "CubeWithMisalignedVolume",
+          Acts::Transform3::Identity(),
+          std::make_shared<Acts::CuboidVolumeBounds>(vBounds), {},
+          {smallVolumeMisaligned}, Acts::Experimental::tryAllSubVolumes(),
+          Acts::Experimental::tryAllPortalsAndSurfaces(), 1000),
+      std::invalid_argument);
+}
+
 BOOST_AUTO_TEST_CASE(CylindricalDetectorVolumePortals) {
   Acts::ActsScalar rInner = 10.;
   Acts::ActsScalar rOuter = 100.;
diff --git a/Tests/UnitTests/Core/Geometry/LayerCreatorTests.cpp b/Tests/UnitTests/Core/Geometry/LayerCreatorTests.cpp
index 14a76a8195e..f150d5734bd 100644
--- a/Tests/UnitTests/Core/Geometry/LayerCreatorTests.cpp
+++ b/Tests/UnitTests/Core/Geometry/LayerCreatorTests.cpp
@@ -29,7 +29,6 @@
 #include "Acts/Utilities/BinningType.hpp"
 #include "Acts/Utilities/IAxis.hpp"
 #include "Acts/Utilities/Logger.hpp"
-#include "Acts/Utilities/Range1D.hpp"
 
 #include <algorithm>
 #include <array>
diff --git a/Tests/UnitTests/Core/Seeding/SeedFinderTest.cpp b/Tests/UnitTests/Core/Seeding/SeedFinderTest.cpp
index 456033ad723..ff616fa0e7c 100644
--- a/Tests/UnitTests/Core/Seeding/SeedFinderTest.cpp
+++ b/Tests/UnitTests/Core/Seeding/SeedFinderTest.cpp
@@ -18,7 +18,6 @@
 #include "Acts/Seeding/SpacePointGrid.hpp"
 #include "Acts/Utilities/Grid.hpp"
 #include "Acts/Utilities/GridBinFinder.hpp"
-#include "Acts/Utilities/Range1D.hpp"
 
 #include <algorithm>
 #include <chrono>
diff --git a/Tests/UnitTests/Core/Utilities/KDTreeTests.cpp b/Tests/UnitTests/Core/Utilities/KDTreeTests.cpp
index fc4453c366a..20979d4dab7 100644
--- a/Tests/UnitTests/Core/Utilities/KDTreeTests.cpp
+++ b/Tests/UnitTests/Core/Utilities/KDTreeTests.cpp
@@ -9,7 +9,6 @@
 #include <boost/test/unit_test.hpp>
 
 #include "Acts/Utilities/KDTree.hpp"
-#include "Acts/Utilities/Range1D.hpp"
 #include "Acts/Utilities/RangeXD.hpp"
 
 #include <algorithm>
diff --git a/Tests/UnitTests/Core/Utilities/Range1DTests.cpp b/Tests/UnitTests/Core/Utilities/Range1DTests.cpp
index 6ec47b31974..17c14bb754c 100644
--- a/Tests/UnitTests/Core/Utilities/Range1DTests.cpp
+++ b/Tests/UnitTests/Core/Utilities/Range1DTests.cpp
@@ -8,7 +8,7 @@
 
 #include <boost/test/unit_test.hpp>
 
-#include "Acts/Utilities/Range1D.hpp"
+#include "Acts/Utilities/RangeXD.hpp"
 
 #include <limits>
 #include <utility>
diff --git a/Tests/UnitTests/Core/Utilities/RangeXDTests.cpp b/Tests/UnitTests/Core/Utilities/RangeXDTests.cpp
index 035237e5cf4..b864b46b229 100644
--- a/Tests/UnitTests/Core/Utilities/RangeXDTests.cpp
+++ b/Tests/UnitTests/Core/Utilities/RangeXDTests.cpp
@@ -8,7 +8,6 @@
 
 #include <boost/test/unit_test.hpp>
 
-#include "Acts/Utilities/Range1D.hpp"
 #include "Acts/Utilities/RangeXD.hpp"
 
 #include <limits>
diff --git a/Tests/UnitTests/Plugins/DD4hep/DD4hepCylindricalDetectorTests.cpp b/Tests/UnitTests/Plugins/DD4hep/DD4hepCylindricalDetectorTests.cpp
index a565bad3dca..894bfb89cce 100644
--- a/Tests/UnitTests/Plugins/DD4hep/DD4hepCylindricalDetectorTests.cpp
+++ b/Tests/UnitTests/Plugins/DD4hep/DD4hepCylindricalDetectorTests.cpp
@@ -160,10 +160,10 @@ Acts::Test::CylindricalTrackingGeometry::DetectorStore generateXML() {
   std::vector<std::array<Acts::ActsScalar, 2u>> innerOuter = {
       {25., 35.}, {65., 75.}, {110., 120.}};
   auto b0Surfaces = cGeometry.surfacesCylinder(dStore, 8.4, 36., 0.15, 0.14,
-                                               32., 3., 2., {16, 14});
+                                               31., 3., 2., {16, 14});
 
   auto b1Surfaces = cGeometry.surfacesCylinder(dStore, 8.4, 36., 0.15, 0.14,
-                                               72., 3., 2., {32, 14});
+                                               71., 3., 2., {32, 14});
 
   auto b2Surfaces = cGeometry.surfacesCylinder(dStore, 8.4, 36., 0.15, 0.14,
                                                116., 3., 2., {52, 14});
@@ -192,7 +192,7 @@ Acts::Test::CylindricalTrackingGeometry::DetectorStore generateXML() {
   cxml << beampipe_head_xml << '\n';
   cxml << indent_12_xml << "<acts_passive_surface>" << '\n';
   cxml << indent_12_xml
-       << "<tubs rmin=\"25*mm\" rmax=\"25.8*mm\" dz=\"1800*mm\" "
+       << "<tubs rmin=\"19*mm\" rmax=\"20*mm\" dz=\"800*mm\" "
           "material=\"Air\"/>"
        << '\n';
   cxml << indent_12_xml << "</acts_passive_surface>" << '\n';
diff --git a/Tests/UnitTests/Plugins/Json/DetectorJsonConverterTests.cpp b/Tests/UnitTests/Plugins/Json/DetectorJsonConverterTests.cpp
index b4508076a2c..39ede9a27f8 100644
--- a/Tests/UnitTests/Plugins/Json/DetectorJsonConverterTests.cpp
+++ b/Tests/UnitTests/Plugins/Json/DetectorJsonConverterTests.cpp
@@ -178,7 +178,7 @@ BOOST_AUTO_TEST_CASE(BeamPipeEndcapBarrelDetector) {
                                              Acts::Logging::VERBOSE));
 
     Acts::Experimental::VolumeStructureBuilder::Config shapeConfig;
-    shapeConfig.boundValues = {20, 100, 10., M_PI, 0.};
+    shapeConfig.boundValues = {18, 100, 10., M_PI, 0.};
     shapeConfig.transform = Acts::Transform3(Acts::Transform3::Identity())
                                 .pretranslate(Acts::Vector3(0., 0., ep));
     shapeConfig.boundsType = Acts::VolumeBounds::BoundsType::eCylinder;
@@ -202,7 +202,7 @@ BOOST_AUTO_TEST_CASE(BeamPipeEndcapBarrelDetector) {
 
   // Central barrel
   Acts::Experimental::VolumeStructureBuilder::Config innerShapeConfig;
-  innerShapeConfig.boundValues = {20., 60., 700., M_PI, 0.};
+  innerShapeConfig.boundValues = {18., 60., 700., M_PI, 0.};
   innerShapeConfig.boundsType = Acts::VolumeBounds::BoundsType::eCylinder;
 
   auto innerShapeBuilder =
@@ -298,7 +298,7 @@ BOOST_AUTO_TEST_CASE(BeamPipeEndcapBarrelDetector) {
 
   // Beam Pipe
   Acts::Experimental::VolumeStructureBuilder::Config bpShapeConfig;
-  bpShapeConfig.boundValues = {0., 20., 720., M_PI, 0.};
+  bpShapeConfig.boundValues = {0., 18., 720., M_PI, 0.};
   bpShapeConfig.boundsType = Acts::VolumeBounds::BoundsType::eCylinder;
 
   auto bpShapeBuilder =
diff --git a/Tests/UnitTests/Plugins/Json/DetectorVolumeJsonConverterTests.cpp b/Tests/UnitTests/Plugins/Json/DetectorVolumeJsonConverterTests.cpp
index 2cf710fc7e8..27ed95d7bf2 100644
--- a/Tests/UnitTests/Plugins/Json/DetectorVolumeJsonConverterTests.cpp
+++ b/Tests/UnitTests/Plugins/Json/DetectorVolumeJsonConverterTests.cpp
@@ -14,6 +14,7 @@
 #include "Acts/Detector/LayerStructureBuilder.hpp"
 #include "Acts/Detector/PortalGenerators.hpp"
 #include "Acts/Detector/VolumeStructureBuilder.hpp"
+#include "Acts/Geometry/ConeVolumeBounds.hpp"
 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Navigation/DetectorVolumeFinders.hpp"
diff --git a/Tests/UnitTests/Plugins/Json/UtilitiesJsonConverterTests.cpp b/Tests/UnitTests/Plugins/Json/UtilitiesJsonConverterTests.cpp
index 3a911304eb2..a2b0e56ea68 100644
--- a/Tests/UnitTests/Plugins/Json/UtilitiesJsonConverterTests.cpp
+++ b/Tests/UnitTests/Plugins/Json/UtilitiesJsonConverterTests.cpp
@@ -13,7 +13,6 @@
 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
 #include "Acts/Utilities/BinUtility.hpp"
 #include "Acts/Utilities/BinningType.hpp"
-#include "Acts/Utilities/Range1D.hpp"
 
 #include <cmath>
 #include <fstream>