diff --git a/Core/include/Acts/Geometry/CuboidVolumeBounds.hpp b/Core/include/Acts/Geometry/CuboidVolumeBounds.hpp
index 7c9622f3f19c2b83cecb60ddfed8007fff68249d..6e8f60f13d2633bcc444019b2e6104e7e6395741 100644
--- a/Core/include/Acts/Geometry/CuboidVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/CuboidVolumeBounds.hpp
@@ -1,6 +1,6 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 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
@@ -8,12 +8,16 @@
 
 #pragma once
 
-#include <cmath>
 #include "Acts/Geometry/Volume.hpp"
 #include "Acts/Geometry/VolumeBounds.hpp"
 #include "Acts/Utilities/BoundingBox.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 
+#include <array>
+#include <cmath>
+#include <exception>
+#include <vector>
+
 namespace Acts {
 
 class RectangleBounds;
@@ -44,35 +48,52 @@ class Surface;
 
 class CuboidVolumeBounds : public VolumeBounds {
  public:
-  /// @enum BoundValues for readability
-  enum BoundValues { bv_halfX = 0, bv_halfY = 1, bv_halfZ = 2, bv_length = 3 };
+  /// @enum BoundValues for streaming and access
+  enum BoundValues : int {
+    eHalfLengthX = 0,
+    eHalfLengthY = 1,
+    eHalfLengthZ = 2,
+    eSize = 3
+  };
 
-  /// Default Constructor
-  CuboidVolumeBounds();
+  CuboidVolumeBounds() = delete;
 
   /// Constructor - the box boundaries
   ///
   /// @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);
+  CuboidVolumeBounds(double halex, double haley, double halez) noexcept(false);
+
+  /// Constructor - from a fixed size array
+  ///
+  /// @param values iw the bound values
+  CuboidVolumeBounds(const std::array<double, eSize>& values) noexcept(false)
+      : m_values(values) {
+    checkConsistency();
+  }
 
   /// Copy Constructor
   ///
   /// @param bobo is the source volume bounds to be copied
   CuboidVolumeBounds(const CuboidVolumeBounds& bobo);
 
-  /// Destructor
-  ~CuboidVolumeBounds() override;
-
   /// Assignment operator
   ///
   /// @param bobo is the source volume bounds to be assigned
   CuboidVolumeBounds& operator=(const CuboidVolumeBounds& bobo);
 
-  /// Virtual constructor
+  ~CuboidVolumeBounds() override = default;
+
   CuboidVolumeBounds* clone() const override;
 
+  VolumeBounds::BoundsType type() const final { return VolumeBounds::eCuboid; }
+
+  /// Return the bound values as dynamically sized vector
+  ///
+  /// @return this returns a copy of the internal values
+  std::vector<double> values() const final;
+
   /// This method checks if position in the 3D volume
   /// frame is inside the cylinder
   ///
@@ -95,42 +116,30 @@ class CuboidVolumeBounds : public VolumeBounds {
                                   const Vector3D& envelope = {0, 0, 0},
                                   const Volume* entity = nullptr) const final;
 
-  /// This method returns the halflength in local x
-  double halflengthX() const;
-
-  /// This method returns the halflength in local y
-  double halflengthY() const;
-
-  /// This method returns the halflength in local z
-  double halflengthZ() const;
-
   /// Output Method for std::ostream
   ///
   /// @param sl is ostream operator to be dumped into
   std::ostream& toStream(std::ostream& sl) 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]; }
+
  private:
   /// Templated dumpT method
   template <class T>
   T& dumpT(T& dt) const;
 
-  /// This method returns the associated RecantleBounds of the face PlaneSurface
-  /// parallel to local xy plane
-  std::shared_ptr<const RectangleBounds> faceXYRectangleBounds() const;
-
-  /// This method returns the associated RecantleBounds of the face PlaneSurface
-  /// parallel to local yz plane
-  std::shared_ptr<const RectangleBounds> faceYZRectangleBounds() const;
+  /// The bound values ordered in a fixed size array
+  std::array<double, eSize> m_values;
 
-  /// This method returns the associated RecantleBounds of the face PlaneSurface
-  // parallel to local zx plane
-  std::shared_ptr<const RectangleBounds> faceZXRectangleBounds() const;
+  std::shared_ptr<const RectangleBounds> m_xyBounds = nullptr;
+  std::shared_ptr<const RectangleBounds> m_yzBounds = nullptr;
+  std::shared_ptr<const RectangleBounds> m_zxBounds = nullptr;
 
-  /// The bound values
-  std::vector<double> m_values;
-  std::shared_ptr<const RectangleBounds> m_xyBounds;
-  std::shared_ptr<const RectangleBounds> m_yzBounds;
-  std::shared_ptr<const RectangleBounds> m_zxBounds;
+  /// Check the input values for consistency,
+  /// will throw a logic_exception if consistency is not given
+  void checkConsistency() noexcept(false);
 };
 
 inline CuboidVolumeBounds* CuboidVolumeBounds::clone() const {
@@ -138,30 +147,32 @@ inline CuboidVolumeBounds* CuboidVolumeBounds::clone() const {
 }
 
 inline bool CuboidVolumeBounds::inside(const Vector3D& pos, double tol) const {
-  return (std::abs(pos.x()) <= m_values.at(bv_halfX) + tol &&
-          std::abs(pos.y()) <= m_values.at(bv_halfY) + tol &&
-          std::abs(pos.z()) <= m_values.at(bv_halfZ) + tol);
-}
-
-inline double CuboidVolumeBounds::halflengthX() const {
-  return m_values.at(bv_halfX);
+  return (std::abs(pos.x()) <= get(eHalfLengthX) + tol &&
+          std::abs(pos.y()) <= get(eHalfLengthY) + tol &&
+          std::abs(pos.z()) <= get(eHalfLengthZ) + tol);
 }
 
-inline double CuboidVolumeBounds::halflengthY() const {
-  return m_values.at(bv_halfY);
+inline std::vector<double> CuboidVolumeBounds::values() const {
+  std::vector<double> valvector;
+  valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
+  return valvector;
 }
 
-inline double CuboidVolumeBounds::halflengthZ() const {
-  return m_values.at(bv_halfZ);
+inline void CuboidVolumeBounds::checkConsistency() noexcept(false) {
+  if (get(eHalfLengthX) <= 0 or get(eHalfLengthY) <= 0 or
+      get(eHalfLengthZ) <= 0.) {
+    throw std::invalid_argument(
+        "CuboidVolumeBounds: invalid input, zero or negative.");
+  }
 }
 
 template <class T>
 T& CuboidVolumeBounds::dumpT(T& dt) const {
   dt << std::setiosflags(std::ios::fixed);
   dt << std::setprecision(5);
-  dt << "Acts::CuboidVolumeBounds: (halfX, halfY, halfZ) = ";
-  dt << "(" << m_values.at(bv_halfX) << ", " << m_values.at(bv_halfY) << ", "
-     << m_values.at(bv_halfZ) << ")";
+  dt << "Acts::CuboidVolumeBounds: (halfLengthX, halfLengthY, halfLengthZ) = ";
+  dt << "(" << get(eHalfLengthX) << ", " << get(eHalfLengthY) << ", "
+     << get(eHalfLengthZ) << ")";
   return dt;
 }
 }  // namespace Acts
\ No newline at end of file
diff --git a/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp b/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp
index 90c9aff5611f13812463de1e612b539b4296b1c6..eaaf81dd2fd47d10cb718f8ceee53c2941e5da26 100644
--- a/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/CutoutCylinderVolumeBounds.hpp
@@ -16,6 +16,10 @@
 #include "Acts/Utilities/BoundingBox.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 
+#include <array>
+#include <exception>
+#include <vector>
+
 namespace Acts {
 
 class IVisualization;
@@ -27,30 +31,59 @@ class IVisualization;
 /// |    |---------|    | rmed
 /// |    |         |    |
 /// ------         ------ rmin
-///       -- dz2 --
-/// -------- dz1 -------
+///       -- hlZc --
+/// --------- hlZ -------
 ///
 ///
 class CutoutCylinderVolumeBounds : public VolumeBounds {
  public:
+  /// @enum BoundValues for streaming and access
+  enum BoundValues : int {
+    eMinR = 0,
+    eMedR = 1,
+    eMaxR = 2,
+    eHalfLengthZ = 3,
+    eHalfLengthZcutout = 4,
+    eSize = 5
+  };
+
+  CutoutCylinderVolumeBounds() = delete;
+
   /// Constructor from defining parameters
   ///
   /// @param rmin Minimum radius at the "choke points"
   /// @param rmed The medium radius (outer radius of the cutout)
   /// @param rmax The outer radius of the overall shape
-  /// @param dz1 The longer halflength of the shape
-  /// @param dz2 The shorter halflength of the shape
-  CutoutCylinderVolumeBounds(double rmin, double rmed, double rmax, double dz1,
-                             double dz2)
-      : m_rmin(rmin), m_rmed(rmed), m_rmax(rmax), m_dz1(dz1), m_dz2(dz2) {}
+  /// @param hlZ The longer halflength of the shape
+  /// @param hlZc The cutout halflength of the shape
+  CutoutCylinderVolumeBounds(double rmin, double rmed, double rmax, double hlZ,
+                             double hlZc) noexcept(false)
+      : m_values({rmin, rmed, rmax, hlZ, hlZc}) {
+    checkConsistency();
+  }
+
+  /// Constructor - from a fixed size array
+  ///
+  /// @param values The bound values
+  CutoutCylinderVolumeBounds(const std::array<double, eSize>& values) noexcept(
+      false)
+      : m_values(values) {
+    checkConsistency();
+  }
 
-  /// Virtual default constructor
   ~CutoutCylinderVolumeBounds() override = default;
 
-  /// Clone method.
-  /// @return Pointer to a copy of the shape
   VolumeBounds* clone() const override;
 
+  VolumeBounds::BoundsType type() const final {
+    return VolumeBounds::eCutoutCylinder;
+  }
+
+  /// Return the bound values as dynamically sized vector
+  ///
+  /// @return this returns a copy of the internal values
+  std::vector<double> values() const final;
+
   /// Inside method to test whether a point is inside the shape
   ///
   /// @param gpos The point to test
@@ -83,32 +116,36 @@ class CutoutCylinderVolumeBounds : public VolumeBounds {
   /// @return The outstream
   std::ostream& toStream(std::ostream& sl) const override;
 
-  /// Return the minimum radius
-  /// @return The minimum radius
-  double rMin() const { return m_rmin; }
-
-  /// Return the medium radius
-  /// @return The medium radius
-  double rMed() const { return m_rmed; }
-
-  /// Return the maximum radius
-  /// @return The maximum radius
-  double rMax() const { return m_rmax; }
-
-  /// Return the longer halflength in z.
-  /// @return The halflength
-  double dZ1() const { return m_dz1; }
-
-  /// Return the shorter halflength in z.
-  /// @return The halflength
-  double dZ2() const { return m_dz2; }
+  /// Access to the bound values
+  /// @param bValue the class nested enum for the array access
+  double get(BoundValues bValue) const { return m_values[bValue]; }
 
  private:
-  double m_rmin;
-  double m_rmed;
-  double m_rmax;
-  double m_dz1;
-  double m_dz2;
+  std::array<double, eSize> m_values;
+
+  /// Check the input values for consistency,
+  /// will throw a logic_exception if consistency is not given
+  void checkConsistency() noexcept(false);
 };
 
+inline std::vector<double> CutoutCylinderVolumeBounds::values() const {
+  std::vector<double> valvector;
+  valvector.insert(valvector.begin(), m_values.begin(), m_values.end());
+  return valvector;
+}
+
+inline void CutoutCylinderVolumeBounds::checkConsistency() noexcept(false) {
+  if (get(eMinR) < 0. or get(eMedR) <= 0. or get(eMaxR) <= 0. or
+      get(eMinR) >= get(eMedR) or get(eMinR) >= get(eMaxR) or
+      get(eMedR) >= get(eMaxR)) {
+    throw std::invalid_argument(
+        "CutoutCylinderVolumeBounds: invalid radial input.");
+  }
+  if (get(eHalfLengthZ) <= 0 or get(eHalfLengthZcutout) <= 0. or
+      get(eHalfLengthZcutout) > get(eHalfLengthZ)) {
+    throw std::invalid_argument(
+        "CutoutCylinderVolumeBounds: invalid longitudinal input.");
+  }
+}
+
 }  // namespace Acts
diff --git a/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp b/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
index 1ed77ae6263438077f4c977db25f25e81f8851a9..a29d239a22a2adf87b3d44942482d57e7137d940 100644
--- a/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/CylinderVolumeBounds.hpp
@@ -7,7 +7,7 @@
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 #pragma once
-#include <cmath>
+
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Geometry/Volume.hpp"
 #include "Acts/Geometry/VolumeBounds.hpp"
@@ -16,12 +16,17 @@
 #include "Acts/Utilities/BoundingBox.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 #include "Acts/Utilities/Helpers.hpp"
+#include "Acts/Utilities/detail/periodic.hpp"
+
+#include <array>
+#include <cmath>
+#include <exception>
+#include <vector>
 
 namespace Acts {
 
 class Surface;
 class CylinderBounds;
-class DiscBounds;
 class RadialBounds;
 class PlanarBounds;
 class IVisualization;
@@ -70,66 +75,76 @@ class IVisualization;
 
 class CylinderVolumeBounds : public VolumeBounds {
  public:
-  /// @enum BoundValues for readability
+  /// @enum BoundValues for streaming and access
   enum BoundValues {
-    bv_innerRadius = 0,
-    bv_outerRadius = 1,
-    bv_halfPhiSector = 2,
-    bv_halfZ = 3,
-    bv_length = 4
+    eMinR = 0,
+    eMaxR = 1,
+    eHalfLengthZ = 2,
+    eHalfPhiSector = 3,
+    eAveragePhi = 4,
+    eSize = 5
   };
 
-  /// Default Constructor
-  CylinderVolumeBounds();
+  CylinderVolumeBounds() = delete;
 
-  /// Constructor - full cylinder
+  /// Constructor
   ///
-  /// @param radius is the outer radius of the cylinder
-  /// @param halez is the half length in z
-  CylinderVolumeBounds(double radius, double halez);
+  /// @param rmin The inner radius of the cylinder
+  /// @param rmax The outer radius of the cylinder
+  /// @param halfz The half length in z
+  /// @param halfphi The half lopening angle
+  /// @param avgphi The average phi value
+  CylinderVolumeBounds(double rmin, double rmax, double halfz,
+                       double halfphi = M_PI,
+                       double avgphi = 0.) noexcept(false)
+      : m_values({rmin, rmax, halfz, halfphi, avgphi}) {
+    checkConsistency();
+    buildSurfaceBounds();
+  }
 
-  /// Constructor - extruded cylinder
+  /// Constructor - from a fixed size array
   ///
-  /// @param rinner is the inner radius of the cylinder
-  /// @param router is the outer radius of the cylinder
-  /// @param halez is the half length in z
-  CylinderVolumeBounds(double rinner, double router, double halez);
+  /// @param values The bound values
+  CylinderVolumeBounds(const std::array<double, eSize>& values) noexcept(false)
+      : m_values(values) {
+    checkConsistency();
+    buildSurfaceBounds();
+  }
 
-  /// Constructor - extruded cylinder
-  ///
-  /// @param rinner is the inner radius of the cylinder
-  /// @param router is the outer radius of the cylinder
-  /// @param haphi is the half opening angle
-  /// @param halez is the half length in z
-  CylinderVolumeBounds(double rinner, double router, double haphi,
-                       double halez);
-
-  /// Constructor - from cylinder bounds and thickness
+  /// Constructor - extruded from cylinder bounds and thickness
   ///
   /// @param cbounds the cylinder bounds
-  /// @param thickness
-  CylinderVolumeBounds(const CylinderBounds& cBounds, double thickness);
+  /// @param thickness of the extrusion
+  CylinderVolumeBounds(const CylinderBounds& cBounds,
+                       double thickness) noexcept(false);
 
-  /// Constructor - from radial bounds and thickness
+  /// Constructor - extruded from radial bounds and thickness
   ///
   /// @param rbounds the Radial bounds
   /// @param thickness
-  CylinderVolumeBounds(const RadialBounds& rBounds, double thickness);
+  CylinderVolumeBounds(const RadialBounds& rBounds,
+                       double thickness) noexcept(false);
 
   /// Copy Constructor
   ///
   /// @param cylbo is the source cylinder volume bounds for the copy
-  CylinderVolumeBounds(const CylinderVolumeBounds& cylbo);
+  CylinderVolumeBounds(const CylinderVolumeBounds& cylbo) = default;
 
-  /// Destructor
-  ~CylinderVolumeBounds() override;
+  ~CylinderVolumeBounds() override = default;
 
-  /// Assignment operator
-  CylinderVolumeBounds& operator=(const CylinderVolumeBounds& cylbo);
+  CylinderVolumeBounds& operator=(const CylinderVolumeBounds& cylbo) = default;
 
-  /// Virtual constructor
   CylinderVolumeBounds* clone() const override;
 
+  VolumeBounds::BoundsType type() const final {
+    return VolumeBounds::eCylinder;
+  }
+
+  /// Return the bound values as dynamically sized vector
+  ///
+  /// @return this returns a copy of the internal values
+  std::vector<double> values() const final;
+
   /// This method checks if position in the 3D volume
   /// frame is inside the cylinder
   ///
@@ -162,54 +177,35 @@ class CylinderVolumeBounds : public VolumeBounds {
   /// @param bValue is the type used for the binning
   double binningBorder(BinningValue bValue) const override;
 
-  /// This method returns the inner radius
-  double innerRadius() const;
-
-  /// This method returns the outer radius
-  double outerRadius() const;
-
-  /// This method returns the medium radius
-  double mediumRadius() const;
-
-  /// This method returns the delta radius
-  double deltaRadius() const;
-
-  /// This method returns the halfPhiSector angle
-  double halfPhiSector() const;
-
-  /// This method returns the halflengthZ
-  double halflengthZ() const;
-
   /// Output Method for std::ostream
   std::ostream& toStream(std::ostream& sl) 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]; }
+
  private:
+  /// The internal version of the bounds can be float/double
+  std::array<double, eSize> m_values;
+  /// Bounds of the inner CylinderSurfaces
+  std::shared_ptr<const CylinderBounds> m_innerCylinderBounds{nullptr};
+  /// Bounds of the inner CylinderSurfaces
+  std::shared_ptr<const CylinderBounds> m_outerCylinderBounds{nullptr};
+  /// Bounds of the bottom/top DiscSurface
+  std::shared_ptr<const RadialBounds> m_discBounds{nullptr};
+  /// Bounds of the sector planes
+  std::shared_ptr<const PlanarBounds> m_sectorPlaneBounds{nullptr};
+
+  /// Check the input values for consistency,
+  /// will throw a logic_exception if consistency is not given
+  void checkConsistency() noexcept(false);
+
+  /// Helper method to create the surface bounds
+  void buildSurfaceBounds();
+
   /// templated dumpT method
   template <class T>
   T& dumpT(T& tstream) const;
-
-  /// This method returns the associated CylinderBounds
-  /// of the inner CylinderSurfaces.
-  std::shared_ptr<const CylinderBounds> innerCylinderBounds() const;
-
-  /// This method returns the associated CylinderBounds
-  /// of the inner CylinderSurfaces.
-  std::shared_ptr<const CylinderBounds> outerCylinderBounds() const;
-
-  /// This method returns the associated RadialBounds
-  /// for the bottom/top DiscSurface
-  std::shared_ptr<const DiscBounds> discBounds() const;
-
-  /// This method returns the associated PlaneBounds
-  /// limiting a sectoral CylinderVolume
-  std::shared_ptr<const PlanarBounds> sectorPlaneBounds() const;
-
-  /// The internal version of the bounds can be float/double
-  std::vector<double> m_values;
-
-  /// numerical stability
-  /// @todo unify the numerical stability checks
-  static const double s_numericalStable;
 };
 
 inline CylinderVolumeBounds* CylinderVolumeBounds::clone() const {
@@ -221,66 +217,66 @@ inline bool CylinderVolumeBounds::inside(const Vector3D& pos,
   using VectorHelpers::perp;
   using VectorHelpers::phi;
   double ros = perp(pos);
-  bool insidePhi = cos(phi(pos)) >= cos(m_values[bv_halfPhiSector]) - tol;
-  bool insideR = insidePhi ? ((ros >= m_values[bv_innerRadius] - tol) &&
-                              (ros <= m_values[bv_outerRadius] + tol))
-                           : false;
+  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()) <= m_values[bv_halfZ] + tol) : false;
+      insideR ? (std::abs(pos.z()) <= get(eHalfLengthZ) + tol) : false;
   return (insideZ && insideR && insidePhi);
 }
 
 inline Vector3D CylinderVolumeBounds::binningOffset(BinningValue bValue)
     const {  // the medium radius is taken for r-type binning
   if (bValue == Acts::binR || bValue == Acts::binRPhi) {
-    return Vector3D(mediumRadius(), 0., 0.);
+    return Vector3D(0.5 * (get(eMinR) + get(eMaxR)), 0., 0.);
   }
   return VolumeBounds::binningOffset(bValue);
 }
 
-inline double CylinderVolumeBounds::binningBorder(BinningValue bValue)
-    const {  // the medium radius is taken for r-type binning
+inline double CylinderVolumeBounds::binningBorder(BinningValue bValue) const {
   if (bValue == Acts::binR) {
-    return 0.5 * deltaRadius();
+    return 0.5 * (get(eMaxR) - get(eMinR));
   }
   if (bValue == Acts::binZ) {
-    return halflengthZ();
+    return get(eHalfLengthZ);
   }
   return VolumeBounds::binningBorder(bValue);
 }
 
-inline double CylinderVolumeBounds::innerRadius() const {
-  return m_values.at(bv_innerRadius);
-}
-
-inline double CylinderVolumeBounds::outerRadius() const {
-  return m_values.at(bv_outerRadius);
-}
-
-inline double CylinderVolumeBounds::mediumRadius() const {
-  return 0.5 * (m_values.at(bv_innerRadius) + m_values.at(bv_outerRadius));
-}
-
-inline double CylinderVolumeBounds::deltaRadius() const {
-  return (m_values.at(bv_outerRadius) - m_values.at(bv_innerRadius));
-}
-
-inline double CylinderVolumeBounds::halfPhiSector() const {
-  return m_values.at(bv_halfPhiSector);
-}
-
-inline double CylinderVolumeBounds::halflengthZ() const {
-  return m_values.at(bv_halfZ);
-}
-
 template <class T>
 T& CylinderVolumeBounds::dumpT(T& tstream) const {
   tstream << std::setiosflags(std::ios::fixed);
   tstream << std::setprecision(5);
-  tstream << "Acts::CylinderVolumeBounds: (rMin, rMax, halfPhi, halfZ) = ";
-  tstream << m_values.at(bv_innerRadius) << ", " << m_values.at(bv_outerRadius)
-          << ", " << m_values.at(bv_halfPhiSector) << ", "
-          << m_values.at(bv_halfZ);
+  tstream << "Acts::CylinderVolumeBounds: (rMin, rMax, halfZ, halfPhi, "
+             "averagePhi) = ";
+  tstream << get(eMinR) << ", " << get(eMaxR) << ", " << get(eHalfLengthZ)
+          << ", " << get(eHalfPhiSector) << get(eAveragePhi);
   return tstream;
 }
+
+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. or get(eMaxR) <= 0. or 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. or 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.");
+  }
+}
+
 }  // namespace Acts
\ No newline at end of file
diff --git a/Core/include/Acts/Geometry/DoubleTrapezoidVolumeBounds.hpp b/Core/include/Acts/Geometry/DoubleTrapezoidVolumeBounds.hpp
deleted file mode 100644
index 3355393eb070bb8f6330d84b8c75ae1810935fb5..0000000000000000000000000000000000000000
--- a/Core/include/Acts/Geometry/DoubleTrapezoidVolumeBounds.hpp
+++ /dev/null
@@ -1,242 +0,0 @@
-// This file is part of the Acts project.
-//
-// Copyright (C) 2016-2018 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/.
-
-///////////////////////////////////////////////////////////////////
-// DoubleTrapezoidVolumeBounds.h, Acts project
-///////////////////////////////////////////////////////////////////
-
-#pragma once
-#include "Acts/Geometry/Volume.hpp"
-#include "Acts/Geometry/VolumeBounds.hpp"
-#include "Acts/Utilities/BoundingBox.hpp"
-#include "Acts/Utilities/Definitions.hpp"
-
-namespace Acts {
-
-class Surface;
-class RectangleBounds;
-class TrapezoidBounds;
-class DiamondBounds;
-
-/// @class DoubleTrapezoidVolumeBounds
-///
-/// Bounds for a double trapezoidal shaped Volume, the decomposeToSurfaces
-/// method
-/// creates a
-/// vector of 8 surfaces:
-///
-///  BoundarySurfaceFace [index]:
-///
-///  - negativeFaceXY     [0] : Diamond Acts::PlaneSurface,
-///                             parallel to \f$ xy \f$ plane at negative \f$z\f$
-///  - positiveFaceXY     [1] : Diamond Acts::PlaneSurface,
-///                             parallel to \f$ xy \f$ plane at positive \f$z\f$
-///  - trapezoidFaceAlpha1 [2] : Rectangular  Acts::PlaneSurface,
-///                              attached to [0] and [1] at negative \f$ x \f$
-/// (associated to alpha1)
-///  - trapezoidFaceBeta1  [3] : Rectangular  Acts::PlaneSurface,
-///                             attached to [0] and [1] at positive \f$ x \f$
-/// (associated to beta1)
-///  - trapezoidFaceAlpha2 [5] : Rectangular  Acts::PlaneSurface,
-///                              attached to [0] and [1] at negative \f$ x \f$
-/// (associated to alpha2)
-///  - trapezoidFaceBeta2  [6] : Rectangular  Acts::PlaneSurface,
-///                              attached to [0] and [1] at positive \f$ x \f$
-/// (associated to beta2)
-///  - negativeFaceZX     [4] : Rectangular  Acts::PlaneSurface,
-///                             parallel to \f$ zx \f$ plane at negative \f$y\f$
-///  - positiveFaceZX     [5] : Rectangular  Acts::PlaneSurface,
-///                             parallel to \f$ zx \f$ plane at positive \f$y\f$
-///
-///  @image html DoubleTrapezoidVolumeBounds_decomp.gif
-
-class DoubleTrapezoidVolumeBounds : public VolumeBounds {
- public:
-  /// @enum BoundValues for readability
-  enum BoundValues {
-    bv_minHalfX = 0,
-    bv_medHalfX = 1,
-    bv_maxHalfX = 2,
-    bv_halfY1 = 3,
-    bv_halfY2 = 4,
-    bv_halfZ = 5,
-    bv_alpha1 = 6,
-    bv_alpha2 = 7,
-    bv_length = 8
-  };
-
-  /// Default Constructor
-  DoubleTrapezoidVolumeBounds();
-
-  /// Constructor - the double trapezoid boundaries (
-  /// symmetric trapezoid/diamond
-  ///
-  /// @param minhalex half length in x at minimum y
-  /// @param medhalex half length in x aty = 0
-  /// @param maxhalex half length in x at maximum x
-  /// @param haley1 first half length in y (to negative)
-  /// @param haley2 second half length in y (to positive)
-  /// @param halez half length in z
-  DoubleTrapezoidVolumeBounds(double minhalex, double medhalex, double maxhalex,
-                              double haley1, double haley2, double halez);
-
-  /// Copy Constructor
-  ///
-  /// @param trabo is the source bounds
-  DoubleTrapezoidVolumeBounds(const DoubleTrapezoidVolumeBounds& trabo);
-
-  /// Destructor
-  ~DoubleTrapezoidVolumeBounds() override;
-
-  /// Assignment operator
-  ///
-  /// @param trabo is the source bounds
-  DoubleTrapezoidVolumeBounds& operator=(
-      const DoubleTrapezoidVolumeBounds& trabo);
-
-  /// Virtual constructor
-  DoubleTrapezoidVolumeBounds* clone() const override;
-
-  /// This method checks if position in the 3D volume frame
-  /// is inside the cylinder
-  ///
-  /// @param pos is the global position to be checked for inside
-  /// @param tol is the tolerance parametere
-  bool inside(const Vector3D& pos, double tol = 0.) const override;
-
-  /// decompose into boundary surfaces
-  ///
-  /// @param transformPtr is the transform applied by the volume
-  ///
-  /// @return a vector of surfaces to be used as boundary surfaces
-  std::vector<std::shared_ptr<const Surface>> decomposeToSurfaces(
-      const Transform3D* transformPtr) const override;
-
-  /// Construct bounding box for this shape
-  /// @param trf Optional transform
-  /// @param envelope Optional envelope to add / subtract from min/max
-  /// @param entity Entity to associate this bounding box with
-  /// @return Constructed bounding box
-  Volume::BoundingBox boundingBox(const Transform3D* trf = nullptr,
-                                  const Vector3D& envelope = {0, 0, 0},
-                                  const Volume* entity = nullptr) const final;
-
-  /// This method returns the X halflength at minimal Y
-  double minHalflengthX() const;
-
-  /// This method returns the (maximal) halflength in local x
-  double medHalflengthX() const;
-
-  /// This method returns the X halflength at maximal Y (local coordinates)
-  double maxHalflengthX() const;
-
-  /// This method returns the halflength1 in local y
-  double halflengthY1() const;
-
-  /// This method returns the halflength2 in local y
-  double halflengthY2() const;
-
-  /// This method returns the halflength in local z
-  double halflengthZ() const;
-
-  /// This method returns the opening angle in point A (negative local x)
-  double alpha1() const;
-
-  /// This method returns the opening angle in point A' (negative local x)
-  double alpha2() const;
-
-  /// Output Method for std::ostream
-  std::ostream& toStream(std::ostream& sl) const override;
-
- private:
-  /// dump method
-  ///
-  /// @tparam dT is the output stream to be dumped into
-  template <class T>
-  T& dumpT(T& dT) const;
-
-  /// This method returns the associated DoubleTrapezoidBounds of the face
-  /// PlaneSurface parallel to local xy plane
-  DiamondBounds* faceXYDiamondBounds() const;
-
-  /// This method returns the associated RecantleBounds of the face PlaneSurface
-  /// attached to alpha (negative local x)
-  RectangleBounds* faceAlpha1RectangleBounds() const;
-
-  /// This method returns the associated RecantleBounds of the face PlaneSurface
-  /// attached to alpha (negative local x)
-  RectangleBounds* faceAlpha2RectangleBounds() const;
-
-  /// This method returns the associated RecantleBounds of the face PlaneSurface
-  /// attached to beta (positive local x)
-  RectangleBounds* faceBeta1RectangleBounds() const;
-
-  /// This method returns the associated RecantleBounds of the face PlaneSurface
-  /// attached to beta (positive local x)
-  RectangleBounds* faceBeta2RectangleBounds() const;
-
-  /// This method returns the associated RecantleBounds of the face PlaneSurface
-  /// parallel to local zx plane, negative local y
-  RectangleBounds* faceZXRectangleBoundsBottom() const;
-
-  /// This method returns the associated RecantleBounds of the face PlaneSurface
-  /// parallel to local zx plane, positive local y
-  RectangleBounds* faceZXRectangleBoundsTop() const;
-
-  std::vector<double> m_values;  ///< the internal store
-};
-
-inline DoubleTrapezoidVolumeBounds* DoubleTrapezoidVolumeBounds::clone() const {
-  return new DoubleTrapezoidVolumeBounds(*this);
-}
-
-inline double DoubleTrapezoidVolumeBounds::minHalflengthX() const {
-  return m_values.at(bv_minHalfX);
-}
-
-inline double DoubleTrapezoidVolumeBounds::medHalflengthX() const {
-  return m_values.at(bv_medHalfX);
-}
-
-inline double DoubleTrapezoidVolumeBounds::maxHalflengthX() const {
-  return m_values.at(bv_maxHalfX);
-}
-
-inline double DoubleTrapezoidVolumeBounds::halflengthY1() const {
-  return m_values.at(bv_halfY1);
-}
-
-inline double DoubleTrapezoidVolumeBounds::halflengthY2() const {
-  return m_values.at(bv_halfY2);
-}
-
-inline double DoubleTrapezoidVolumeBounds::halflengthZ() const {
-  return m_values.at(bv_halfZ);
-}
-
-inline double DoubleTrapezoidVolumeBounds::alpha1() const {
-  return m_values.at(bv_alpha1);
-}
-
-inline double DoubleTrapezoidVolumeBounds::alpha2() const {
-  return m_values.at(bv_alpha2);
-}
-
-template <class T>
-T& DoubleTrapezoidVolumeBounds::dumpT(T& dT) const {
-  dT << std::setiosflags(std::ios::fixed);
-  dT << std::setprecision(5);
-  dT << "Acts::DoubleTrapezoidVolumeBounds: (minhalfX, medhalfX, maxhalfX, "
-        "halfY1, halfY2, halfZ) = ";
-  dT << "(" << m_values.at(bv_minHalfX) << ", " << m_values.at(bv_medHalfX)
-     << ", " << m_values.at(bv_maxHalfX);
-  dT << ", " << m_values.at(bv_halfY1) << ", " << m_values.at(bv_halfY2) << ", "
-     << m_values.at(bv_halfZ) << ")";
-  return dT;
-}
-}  // namespace Acts
\ No newline at end of file
diff --git a/Core/include/Acts/Geometry/GenericCuboidVolumeBounds.hpp b/Core/include/Acts/Geometry/GenericCuboidVolumeBounds.hpp
index 9d0df1e93b056149f8eace94795db1ad2ef237de..7171d1be869c34f0f8b7a103ea7357c2ed317e08 100644
--- a/Core/include/Acts/Geometry/GenericCuboidVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/GenericCuboidVolumeBounds.hpp
@@ -8,35 +8,55 @@
 
 #pragma once
 
-#include <array>
-#include <ostream>
-
 #include "Acts/Geometry/Volume.hpp"
 #include "Acts/Geometry/VolumeBounds.hpp"
 #include "Acts/Utilities/BoundingBox.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 
+#include <array>
+#include <ostream>
+#include <vector>
+
 namespace Acts {
 
 class IVisualization;
 
 class GenericCuboidVolumeBounds : public VolumeBounds {
  public:
+  static constexpr size_t eSize = 24;
+
   GenericCuboidVolumeBounds() = delete;
 
   /// Constructor from a set of vertices
+  ///
+  /// @param vertices The set of input vertices
+  ///
   /// The ordering is considered to be:
   /// - the first 4 vertices are the "top" face
   /// - the second 4 vertices are the "bottom" face
   /// - both faces are given in counter clock wise order
-  GenericCuboidVolumeBounds(const std::array<Acts::Vector3D, 8>& vertices);
+  GenericCuboidVolumeBounds(
+      const std::array<Acts::Vector3D, 8>& vertices) noexcept(false);
+
+  /// Constructor from a fixed size array
+  ///
+  /// @param values The input values
+  GenericCuboidVolumeBounds(const std::array<double, eSize>& values) noexcept(
+      false);
 
   ~GenericCuboidVolumeBounds() override = default;
 
-  ///  clone() method to make deep copy in Volume copy constructor and for
-  /// assigment operator  of the Surface class.
   VolumeBounds* clone() const override;
 
+  VolumeBounds::BoundsType type() const final {
+    return VolumeBounds::eGenericCuboid;
+  }
+
+  /// Return the bound values as dynamically sized vector
+  ///
+  /// @return this returns a copy of the internal values
+  std::vector<double> values() const final;
+
   /// Checking if position given in volume frame is inside
   ///
   /// @param gpos is the global position to be checked
@@ -56,31 +76,32 @@ class GenericCuboidVolumeBounds : public VolumeBounds {
   std::vector<std::shared_ptr<const Surface>> decomposeToSurfaces(
       const Transform3D* transform) const override;
 
-  /**
-   * Construct bounding box for this shape
-   * @param trf Optional transform
-   * @param envelope Optional envelope to add / subtract from min/max
-   * @param entity Entity to associate this bounding box with
-   * @return Constructed bounding box
-   */
+  /// Construct bounding box for this shape
+  /// @param trf Optional transform
+  /// @param envelope Optional envelope to add / subtract from min/max
+  /// @param entity Entity to associate this bounding box with
+  /// @return Constructed bounding box
   Volume::BoundingBox boundingBox(const Transform3D* trf = nullptr,
                                   const Vector3D& envelope = {0, 0, 0},
                                   const Volume* entity = nullptr) const final;
 
-  ///
   /// @param sl is the output stream to be written into
   std::ostream& toStream(std::ostream& sl) const override;
 
-  /**
-   * Draw this shape using a visualization helper
-   * @param helper The visualizatin helper
-   * @param transform Optional transformation matrix
-   */
+  /// Draw this shape using a visualization helper
+  /// @param helper The visualizatin helper
+  /// @param transform Optional transformation matrix
+  ///
   void draw(IVisualization& helper,
             const Transform3D& transform = Transform3D::Identity()) const;
 
  private:
   std::array<Vector3D, 8> m_vertices;
   std::array<Vector3D, 6> m_normals;
+
+  /// Private helper method to contruct the Volume bounds
+  /// to be called by the constructors, from the ordered input vertices
+  void construct() noexcept(false);
+  ;
 };
 }  // namespace Acts
diff --git a/Core/include/Acts/Geometry/TrapezoidVolumeBounds.hpp b/Core/include/Acts/Geometry/TrapezoidVolumeBounds.hpp
index eb2320d96f12a9c0bc5446986c7f52cdd2ab2aee..568db1dfa921b910a863662a9543b99a19b7c2b3 100644
--- a/Core/include/Acts/Geometry/TrapezoidVolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/TrapezoidVolumeBounds.hpp
@@ -11,11 +11,15 @@
 ///////////////////////////////////////////////////////////////////
 
 #pragma once
+
 #include "Acts/Geometry/Volume.hpp"
 #include "Acts/Geometry/VolumeBounds.hpp"
 #include "Acts/Utilities/BoundingBox.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 
+#include <array>
+#include <vector>
+
 namespace Acts {
 
 class Surface;
@@ -48,20 +52,18 @@ class TrapezoidBounds;
 
 class TrapezoidVolumeBounds : public VolumeBounds {
  public:
-  /// @enum BoundValues for readability
+  /// @enum BoundValues for acces / streaming
   enum BoundValues {
-    bv_minHalfX = 0,  //!< minimal halflength in x
-    bv_maxHalfX = 1,  //!< maximal halflength in x
-    bv_halfY = 2,     //!< halflength in y
-    bv_halfZ = 3,     //!< halflength in z
-    bv_alpha = 4,     //!< opening angle alpha (in point A)
-    bv_beta = 5,      //!< opening angle beta  (in point B)
-    bv_length = 6     // length of the bounds vector
-
+    eHalfLengthXnegY = 0,  //!< halflength in x at negative y
+    eHalfLengthXposY = 1,  //!< halflength in x at positive y
+    eHalfLengthY = 2,      //!< halflength in y
+    eHalfLengthZ = 3,      //!< halflength in z
+    eAlpha = 4,            //!< opening angle alpha (in point A)
+    eBeta = 5,             //!< opening angle beta  (in point B)
+    eSize = 6              //!< length of the bounds vector
   };
 
-  /// Default Constructor
-  TrapezoidVolumeBounds();
+  TrapezoidVolumeBounds() = delete;
 
   /// Constructor - the trapezoid boundaries (symmetric trapezoid)
   ///
@@ -70,7 +72,7 @@ class TrapezoidVolumeBounds : public VolumeBounds {
   /// @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);
+                        double halez) noexcept(false);
 
   /// Constructor - the trapezoid boundaries (arbitrary trapezoid)
   ///
@@ -80,22 +82,35 @@ class TrapezoidVolumeBounds : public VolumeBounds {
   /// @param alpha is the openeing angle at -x,-y
   /// @param beta is the openeing angle at +x,-y
   TrapezoidVolumeBounds(double minhalex, double haley, double halez,
-                        double alpha, double beta);
+                        double alpha, double beta) noexcept(false);
+
+  /// Constructor - from a fixed size array
+  ///
+  /// @param values The bound values
+  TrapezoidVolumeBounds(const std::array<double, eSize>& values) noexcept(false)
+      : m_values(values) {
+    checkConsistency();
+    buildSurfaceBounds();
+  }
 
-  /// Copy Constructor
-  /// @param trabo The object to be copied
-  TrapezoidVolumeBounds(const TrapezoidVolumeBounds& trabo);
+  TrapezoidVolumeBounds(const TrapezoidVolumeBounds& trabo) = default;
 
-  /// Destructor
-  ~TrapezoidVolumeBounds() override;
+  TrapezoidVolumeBounds& operator=(const TrapezoidVolumeBounds& trabo) =
+      default;
 
-  /// Assignment operator
-  /// @param trabo The object to be assigned
-  TrapezoidVolumeBounds& operator=(const TrapezoidVolumeBounds& trabo);
+  ~TrapezoidVolumeBounds() override = default;
 
-  /// Virtual constructor
   TrapezoidVolumeBounds* clone() const override;
 
+  VolumeBounds::BoundsType type() const final {
+    return VolumeBounds::eTrapezoid;
+  }
+
+  /// Return the bound values as dynamically sized vector
+  ///
+  /// @return this returns a copy of the internal values
+  std::vector<double> values() const final;
+
   /// This method checks if position in the 3D volume frame
   /// is inside the cylinder
   ///
@@ -123,88 +138,70 @@ class TrapezoidVolumeBounds : public VolumeBounds {
                                   const Vector3D& envelope = {0, 0, 0},
                                   const Volume* entity = nullptr) const final;
 
-  /// This method returns the minimal halflength in local x
-  double minHalflengthX() const;
-
-  /// This method returns the maximal halflength in local x
-  double maxHalflengthX() const;
-
-  /// This method returns the halflength in local y
-  double halflengthY() const;
-
-  /// This method returns the halflength in local z
-  double halflengthZ() const;
-
-  /// This method returns the opening angle in point A (negative local x)
-  double alpha() const;
-
-  /// This method returns the opening angle in point B (negative local x)
-  double beta() const;
-
   /// Output Method for std::ostream
   std::ostream& toStream(std::ostream& sl) 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]; }
+
  private:
+  /// The internal version of the bounds can be float/double
+  std::array<double, 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)
+  std::shared_ptr<const RectangleBounds> m_faceAlphaRectangleBounds = nullptr;
+  /// The face PlaneSurface attached to beta (positive local x)
+  std::shared_ptr<const RectangleBounds> m_faceBetaRectangleBounds = nullptr;
+  /// The face PlaneSurface parallel to local zx plane, negative local y
+  std::shared_ptr<const RectangleBounds> m_faceZXRectangleBounds = nullptr;
+
+  /// Check the input values for consistency,
+  /// will throw a logic_exception if consistency is not given
+  void checkConsistency() noexcept(false);
+
+  /// Helper method to create the surface bounds
+  void buildSurfaceBounds();
+
   /// Templated dump methos
   template <class T>
   T& dumpT(T& dt) const;
-
-  /// This method returns the associated TrapezoidBounds
-  /// of the face PlaneSurface parallel to local xy plane
-  TrapezoidBounds* faceXYTrapezoidBounds() const;
-
-  /// This method returns the associated RecangleBounds
-  /// of the face PlaneSurface attached to alpha (negative local x)
-  RectangleBounds* faceAlphaRectangleBounds() const;
-
-  /// This method returns the associated RectangleBounds
-  /// of the face PlaneSurface attached to beta (positive local x)
-  RectangleBounds* faceBetaRectangleBounds() const;
-
-  /// This method returns the associated RectangleBounds
-  /// of the face PlaneSurface parallel to local zx plane, negative local y
-  RectangleBounds* faceZXRectangleBoundsBottom() const;
-
-  /// This method returns the associated RectangleBounds
-  /// of the face PlaneSurface parallel to local zx plane, positive local y
-  RectangleBounds* faceZXRectangleBoundsTop() const;
-
-  /// the bounds values
-  std::vector<double> m_values;
 };
 
 inline TrapezoidVolumeBounds* TrapezoidVolumeBounds::clone() const {
   return new TrapezoidVolumeBounds(*this);
 }
 
-inline double TrapezoidVolumeBounds::minHalflengthX() const {
-  return m_values.at(bv_minHalfX);
-}
-inline double TrapezoidVolumeBounds::maxHalflengthX() const {
-  return m_values.at(bv_maxHalfX);
-}
-inline double TrapezoidVolumeBounds::halflengthY() const {
-  return m_values.at(bv_halfY);
-}
-inline double TrapezoidVolumeBounds::halflengthZ() const {
-  return m_values.at(bv_halfZ);
-}
-inline double TrapezoidVolumeBounds::alpha() const {
-  return m_values.at(bv_alpha);
-}
-inline double TrapezoidVolumeBounds::beta() const {
-  return m_values.at(bv_beta);
-}
-
 template <class T>
 T& TrapezoidVolumeBounds::dumpT(T& dt) const {
   dt << std::setiosflags(std::ios::fixed);
   dt << std::setprecision(5);
   dt << "Acts::TrapezoidVolumeBounds: (minhalfX, halfY, halfZ, alpha, beta) "
         "= ";
-  dt << "(" << m_values.at(bv_minHalfX) << ", " << m_values.at(bv_halfY) << ", "
-     << m_values.at(bv_halfZ);
-  dt << ", " << m_values.at(bv_alpha) << ", " << m_values.at(bv_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. or 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
\ No newline at end of file
diff --git a/Core/include/Acts/Geometry/VolumeBounds.hpp b/Core/include/Acts/Geometry/VolumeBounds.hpp
index c39a141b9d44371931bd090cc5025a7a0d0c3c2a..968b61e8c45893330ca5bbccabaf99c767b9ee8d 100644
--- a/Core/include/Acts/Geometry/VolumeBounds.hpp
+++ b/Core/include/Acts/Geometry/VolumeBounds.hpp
@@ -15,20 +15,6 @@
 #include "Acts/Utilities/BoundingBox.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 
-#ifndef VOLUMEBOUNDS_boundValues_FILL
-#define VOLUMEBOUNDS_boundValues_FILL(val) m_values[bv_##val] = val
-#endif
-
-#ifndef VOLUMEBOUNDS_boundValues_ACCESS
-#define VOLUMEBOUNDS_boundValues_ACCESS(val) \
-  double val() const { return m_values[bv_##val]; }
-#endif
-
-#ifndef VOLUMEBOUNDS_DERIVED_ACCESS
-#define VOLUMEBOUNDS_DERIVED_ACCESS(derived) \
-  double derived() const { return m_##derived; }
-#endif
-
 namespace Acts {
 
 class Surface;
@@ -39,6 +25,7 @@ using VolumeBoundsPtr = std::shared_ptr<const VolumeBounds>;
 
 using SurfacePtr = std::shared_ptr<const Surface>;
 using SurfacePtrVector = std::vector<SurfacePtr>;
+
 /// @class VolumeBounds
 ///
 /// Pure Absract Base Class for Volume bounds.
@@ -55,16 +42,38 @@ using SurfacePtrVector = std::vector<SurfacePtr>;
 /// Surfaces into BoundarySurfaces.
 class VolumeBounds {
  public:
-  /// Default Constructor*/
+  // @enum BoundsType
+  /// This is nested to the VolumeBounds, as also SurfaceBounds will have
+  /// Bounds Type.
+  enum BoundsType : int {
+    eCone = 0,
+    eCuboid = 1,
+    eCutoutCylinder = 2,
+    eCylinder = 3,
+    eGenericCuboid = 4,
+    eTrapezoid = 5,
+    eOther = 6
+  };
+
   VolumeBounds() = default;
 
-  /// Destructor
   virtual ~VolumeBounds() = default;
 
   ///  clone() method to make deep copy in Volume copy constructor and for
   /// assigment operator  of the Surface class.
   virtual VolumeBounds* clone() const = 0;
 
+  /// Return the bounds type - for persistency optimization
+  ///
+  /// @return is a BoundsType enum
+  virtual BoundsType type() const = 0;
+
+  /// Access method for bound values, this is a dynamically sized
+  /// vector containing the parameters needed to describe these bounds
+  ///
+  /// @return of the stored values for this SurfaceBounds object
+  virtual std::vector<double> values() const = 0;
+
   /// Checking if position given in volume frame is inside
   ///
   /// @param gpos is the global position to be checked
@@ -125,4 +134,11 @@ inline double VolumeBounds::binningBorder(BinningValue /*bValue*/) const {
 /// Overload of << operator for std::ostream for debug output
 std::ostream& operator<<(std::ostream& sl, const VolumeBounds& vb);
 
+inline bool operator==(const VolumeBounds& lhs, const VolumeBounds& rhs) {
+  if (&lhs == &rhs) {
+    return true;
+  }
+  return (lhs.type() == rhs.type()) && (lhs.values() == rhs.values());
+}
+
 }  // namespace Acts
diff --git a/Core/include/Acts/Surfaces/CylinderBounds.hpp b/Core/include/Acts/Surfaces/CylinderBounds.hpp
index ffc4a545b0e5907fe7ba1a367f4d2980dd0ad424..d5c9a10b24fe57ad60131296cf3bf2692a9fa031 100644
--- a/Core/include/Acts/Surfaces/CylinderBounds.hpp
+++ b/Core/include/Acts/Surfaces/CylinderBounds.hpp
@@ -14,6 +14,7 @@
 
 #include <array>
 #include <cmath>
+#include <iostream>
 #include <vector>
 
 namespace Acts {
diff --git a/Core/include/Acts/Surfaces/SurfaceBounds.hpp b/Core/include/Acts/Surfaces/SurfaceBounds.hpp
index 75a72e7846ce732a5e2555c5847375fd37b65196..6e16f30e28046373aa24443775e64fcc4033d731 100644
--- a/Core/include/Acts/Surfaces/SurfaceBounds.hpp
+++ b/Core/include/Acts/Surfaces/SurfaceBounds.hpp
@@ -26,8 +26,8 @@ namespace Acts {
 class SurfaceBounds {
  public:
   /// @enum BoundsType
-  ///
-  /// This enumerator simplifies the persistency and bounds identification
+  /// This is nested to the SurfaceBounds, as also VolumeBounds will have
+  /// Bounds Type.
   enum BoundsType : int {
     eCone = 0,
     eCylinder = 1,
diff --git a/Core/src/Geometry/CMakeLists.txt b/Core/src/Geometry/CMakeLists.txt
index f4669b6d20d9680efcd5c20c5d3fe1dafd3bb685..3e053e26f5adfa6d77623a3a20eb69e4bc6fec66 100644
--- a/Core/src/Geometry/CMakeLists.txt
+++ b/Core/src/Geometry/CMakeLists.txt
@@ -13,7 +13,6 @@ target_sources_local(
     CylinderVolumeHelper.cpp
     Extent.cpp
     DiscLayer.cpp
-    DoubleTrapezoidVolumeBounds.cpp
     GenericApproachDescriptor.cpp
     GenericCuboidVolumeBounds.cpp
     GeometryID.cpp
diff --git a/Core/src/Geometry/CuboidVolumeBounds.cpp b/Core/src/Geometry/CuboidVolumeBounds.cpp
index 219c233efa070a52e59c5dc263679b614bffdbdf..b36a9987c4b5e7c7e0d11baab7d87470cabe1221 100644
--- a/Core/src/Geometry/CuboidVolumeBounds.cpp
+++ b/Core/src/Geometry/CuboidVolumeBounds.cpp
@@ -1,36 +1,28 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 CERN for the benefit of the Acts project
+// Copyright (C) 2016-2020 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/.
 
 #include "Acts/Geometry/CuboidVolumeBounds.hpp"
-#include <cmath>
-#include <iostream>
 #include "Acts/Surfaces/PlaneSurface.hpp"
 #include "Acts/Surfaces/RectangleBounds.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/BoundingBox.hpp"
 
-Acts::CuboidVolumeBounds::CuboidVolumeBounds()
-    : VolumeBounds(), m_values(bv_length, 0.) {}
+#include <cmath>
+#include <iostream>
 
 Acts::CuboidVolumeBounds::CuboidVolumeBounds(double halex, double haley,
                                              double halez)
     : VolumeBounds(),
-      m_values(bv_length, 0.),
-      m_xyBounds(nullptr),
-      m_yzBounds(nullptr),
-      m_zxBounds(nullptr) {
-  m_values.at(bv_halfX) = halex;
-  m_values.at(bv_halfY) = haley;
-  m_values.at(bv_halfZ) = halez;
-
-  m_xyBounds = faceXYRectangleBounds();
-  m_yzBounds = faceYZRectangleBounds();
-  m_zxBounds = faceZXRectangleBounds();
+      m_values({halex, haley, halez}),
+      m_xyBounds(std::make_shared<const RectangleBounds>(halex, haley)),
+      m_yzBounds(std::make_shared<const RectangleBounds>(haley, halez)),
+      m_zxBounds(std::make_shared<const RectangleBounds>(halez, halex)) {
+  checkConsistency();
 }
 
 Acts::CuboidVolumeBounds::CuboidVolumeBounds(const CuboidVolumeBounds& bobo)
@@ -40,8 +32,6 @@ Acts::CuboidVolumeBounds::CuboidVolumeBounds(const CuboidVolumeBounds& bobo)
       m_yzBounds(bobo.m_yzBounds),
       m_zxBounds(bobo.m_zxBounds) {}
 
-Acts::CuboidVolumeBounds::~CuboidVolumeBounds() = default;
-
 Acts::CuboidVolumeBounds& Acts::CuboidVolumeBounds::operator=(
     const CuboidVolumeBounds& bobo) {
   if (this != &bobo) {
@@ -66,12 +56,12 @@ Acts::SurfacePtrVector Acts::CuboidVolumeBounds::decomposeToSurfaces(
   //   (1) - at negative local z
   tTransform =
       new Transform3D(transform * AngleAxis3D(M_PI, Vector3D(0., 1., 0.)) *
-                      Translation3D(Vector3D(0., 0., halflengthZ())));
+                      Translation3D(Vector3D(0., 0., get(eHalfLengthZ))));
   rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
       std::shared_ptr<const Transform3D>(tTransform), m_xyBounds));
   //   (2) - at positive local z
-  tTransform = new Transform3D(transform *
-                               Translation3D(Vector3D(0., 0., halflengthZ())));
+  tTransform = new Transform3D(
+      transform * Translation3D(Vector3D(0., 0., get(eHalfLengthZ))));
   rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
       std::shared_ptr<const Transform3D>(tTransform), m_xyBounds));
   // face surfaces yz -------------------------------------
@@ -79,57 +69,38 @@ Acts::SurfacePtrVector Acts::CuboidVolumeBounds::decomposeToSurfaces(
   //   (3) - at negative local x
   tTransform =
       new Transform3D(transform * AngleAxis3D(M_PI, Vector3D(0., 0., 1.)) *
-                      Translation3D(Vector3D(halflengthX(), 0., 0)) *
+                      Translation3D(Vector3D(get(eHalfLengthX), 0., 0)) *
                       AngleAxis3D(0.5 * M_PI, Vector3D(0., 1., 0)) *
                       AngleAxis3D(0.5 * M_PI, Vector3D(0., 0., 1.)));
   rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
       std::shared_ptr<const Transform3D>(tTransform), m_yzBounds));
   //   (4) - at positive local x
-  tTransform = new Transform3D(transform *
-                               Translation3D(Vector3D(halflengthX(), 0., 0.)) *
-                               AngleAxis3D(0.5 * M_PI, Vector3D(0., 1., 0.)) *
-                               AngleAxis3D(0.5 * M_PI, Vector3D(0., 0., 1.)));
+  tTransform = new Transform3D(
+      transform * Translation3D(Vector3D(get(eHalfLengthX), 0., 0.)) *
+      AngleAxis3D(0.5 * M_PI, Vector3D(0., 1., 0.)) *
+      AngleAxis3D(0.5 * M_PI, Vector3D(0., 0., 1.)));
   rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
       std::shared_ptr<const Transform3D>(tTransform), m_yzBounds));
   // face surfaces zx -------------------------------------
   //   (5) - at negative local y
   tTransform =
       new Transform3D(transform * AngleAxis3D(M_PI, Vector3D(1., 0., 0.)) *
-                      Translation3D(Vector3D(0., halflengthY(), 0.)) *
+                      Translation3D(Vector3D(0., get(eHalfLengthY), 0.)) *
                       AngleAxis3D(-0.5 * M_PI, Vector3D(0., 1., 0.)) *
                       AngleAxis3D(-0.5 * M_PI, Vector3D(1., 0., 0.)));
   rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
       std::shared_ptr<const Transform3D>(tTransform), m_zxBounds));
   //   (6) - at positive local y
-  tTransform = new Transform3D(transform *
-                               Translation3D(Vector3D(0., halflengthY(), 0.)) *
-                               AngleAxis3D(-0.5 * M_PI, Vector3D(0., 1., 0.)) *
-                               AngleAxis3D(-0.5 * M_PI, Vector3D(1., 0., 0.)));
+  tTransform = new Transform3D(
+      transform * Translation3D(Vector3D(0., get(eHalfLengthY), 0.)) *
+      AngleAxis3D(-0.5 * M_PI, Vector3D(0., 1., 0.)) *
+      AngleAxis3D(-0.5 * M_PI, Vector3D(1., 0., 0.)));
   rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
       std::shared_ptr<const Transform3D>(tTransform), m_zxBounds));
   // return the surfaces
   return rSurfaces;
 }
 
-std::shared_ptr<const Acts::RectangleBounds>
-Acts::CuboidVolumeBounds::faceXYRectangleBounds() const {
-  return std::make_shared<const RectangleBounds>(m_values.at(bv_halfX),
-                                                 m_values.at(bv_halfY));
-}
-
-std::shared_ptr<const Acts::RectangleBounds>
-Acts::CuboidVolumeBounds::faceYZRectangleBounds() const {
-  return std::make_shared<const RectangleBounds>(m_values.at(bv_halfY),
-                                                 m_values.at(bv_halfZ));
-}
-
-std::shared_ptr<const Acts::RectangleBounds>
-Acts::CuboidVolumeBounds::faceZXRectangleBounds() const {
-  return std::make_shared<const RectangleBounds>(m_values.at(bv_halfZ),
-                                                 m_values.at(bv_halfX));
-}
-
-// ostream operator overload
 std::ostream& Acts::CuboidVolumeBounds::toStream(std::ostream& sl) const {
   return dumpT(sl);
 }
@@ -137,8 +108,8 @@ std::ostream& Acts::CuboidVolumeBounds::toStream(std::ostream& sl) const {
 Acts::Volume::BoundingBox Acts::CuboidVolumeBounds::boundingBox(
     const Acts::Transform3D* trf, const Vector3D& envelope,
     const Volume* entity) const {
-  Vector3D vmin(-halflengthX(), -halflengthY(), -halflengthZ());
-  Vector3D vmax(halflengthX(), halflengthY(), halflengthZ());
+  Vector3D vmin(-get(eHalfLengthX), -get(eHalfLengthY), -get(eHalfLengthZ));
+  Vector3D vmax(get(eHalfLengthX), get(eHalfLengthY), get(eHalfLengthZ));
 
   Volume::BoundingBox box(entity, vmin - envelope, vmax + envelope);
   return trf == nullptr ? box : box.transformed(*trf);
diff --git a/Core/src/Geometry/CutoutCylinderVolumeBounds.cpp b/Core/src/Geometry/CutoutCylinderVolumeBounds.cpp
index 81d89b835e9ff279bf186d23f43949cba670cb28..6152315bb83c134543d1d29a235801dcd50f009b 100644
--- a/Core/src/Geometry/CutoutCylinderVolumeBounds.cpp
+++ b/Core/src/Geometry/CutoutCylinderVolumeBounds.cpp
@@ -30,8 +30,8 @@ bool Acts::CutoutCylinderVolumeBounds::inside(const Acts::Vector3D& gpos,
   using VectorHelpers::phi;
   double ros = perp(gpos);
 
-  bool insideR = (ros >= m_rmin - tol) && (ros <= m_rmax + tol);
-  bool insideZ = std::abs(gpos.z()) <= m_dz1 + tol;
+  bool insideR = (ros >= get(eMinR) - tol) && (ros <= get(eMaxR) + tol);
+  bool insideZ = std::abs(gpos.z()) <= get(eHalfLengthZ) + tol;
 
   if (!insideR || !insideZ) {
     return false;
@@ -39,8 +39,8 @@ bool Acts::CutoutCylinderVolumeBounds::inside(const Acts::Vector3D& gpos,
 
   // we're inside the outer volume, but we might be in inside the
   // cutout section in the middle
-  bool insideRInner = ros <= m_rmed - tol;
-  bool insideZInner = std::abs(gpos.z()) < m_dz2 - tol;
+  bool insideRInner = ros <= get(eMedR) - tol;
+  bool insideZInner = std::abs(gpos.z()) < get(eHalfLengthZcutout) - tol;
 
   return !insideRInner || !insideZInner;  // we are not, inside bounds
 }
@@ -58,66 +58,68 @@ Acts::CutoutCylinderVolumeBounds::decomposeToSurfaces(
     trf = std::make_shared<const Transform3D>(Transform3D::Identity());
   }
 
-  if (m_rmin == 0.) {
+  if (get(eMinR) == 0.) {
     surfaces.resize(6);  // exactly six surfaces (no choke inner cover)
   } else {
     surfaces.resize(8);  // exactly eight surfaces
   }
 
   // outer cylinder envelope
-  auto outerBounds = std::make_shared<CylinderBounds>(m_rmax, m_dz1);
+  auto outerBounds =
+      std::make_shared<CylinderBounds>(get(eMaxR), get(eHalfLengthZ));
   auto outer = Surface::makeShared<CylinderSurface>(trf, outerBounds);
   surfaces.at(tubeOuterCover) = outer;
 
   // inner (small) cylinder envelope
-  auto ctr_innerBounds = std::make_shared<CylinderBounds>(m_rmed, m_dz2);
+  auto ctr_innerBounds =
+      std::make_shared<CylinderBounds>(get(eMedR), get(eHalfLengthZcutout));
   auto ctr_inner = Surface::makeShared<CylinderSurface>(trf, ctr_innerBounds);
   surfaces.at(tubeInnerCover) = ctr_inner;
 
   // z position of the pos and neg choke points
-  double hlChoke = (m_dz1 - m_dz2) * 0.5;
-  double zChoke = m_dz2 + hlChoke;
+  double hlChoke = (get(eHalfLengthZ) - get(eHalfLengthZcutout)) * 0.5;
+  double zChoke = get(eHalfLengthZcutout) + hlChoke;
 
-  if (m_rmin > 0.) {
+  if (get(eMinR) > 0.) {
     auto posChokeTrf = std::make_shared<const Transform3D>(
         *trf * Translation3D(Vector3D(0, 0, zChoke)));
     auto posInner =
-        Surface::makeShared<CylinderSurface>(posChokeTrf, m_rmin, hlChoke);
+        Surface::makeShared<CylinderSurface>(posChokeTrf, get(eMinR), hlChoke);
     surfaces.at(index7) = posInner;
 
     auto negChokeTrf = std::make_shared<const Transform3D>(
         *trf * Translation3D(Vector3D(0, 0, -zChoke)));
     auto negInner =
-        Surface::makeShared<CylinderSurface>(negChokeTrf, m_rmin, hlChoke);
+        Surface::makeShared<CylinderSurface>(negChokeTrf, get(eMinR), hlChoke);
     surfaces.at(index6) = negInner;
   }
 
   // outer disks
   auto posOutDiscTrf = std::make_shared<const Transform3D>(
-      *trf * Translation3D(Vector3D(0, 0, m_dz1)));
+      *trf * Translation3D(Vector3D(0, 0, get(eHalfLengthZ))));
   auto posOutDisc =
-      Surface::makeShared<DiscSurface>(posOutDiscTrf, m_rmin, m_rmax);
+      Surface::makeShared<DiscSurface>(posOutDiscTrf, get(eMinR), get(eMaxR));
   surfaces.at(positiveFaceXY) = posOutDisc;
 
   auto negOutDiscTrf = std::make_shared<const Transform3D>(
-      *trf * Translation3D(Vector3D(0, 0, -m_dz1)) *
+      *trf * Translation3D(Vector3D(0, 0, -get(eHalfLengthZ))) *
       AngleAxis3D(M_PI, Vector3D::UnitX()));
   auto negOutDisc =
-      Surface::makeShared<DiscSurface>(negOutDiscTrf, m_rmin, m_rmax);
+      Surface::makeShared<DiscSurface>(negOutDiscTrf, get(eMinR), get(eMaxR));
   surfaces.at(negativeFaceXY) = negOutDisc;
 
   // inner disks
   auto posInDiscTrf = std::make_shared<const Transform3D>(
-      *trf * Translation3D(Vector3D(0, 0, m_dz2)));
+      *trf * Translation3D(Vector3D(0, 0, get(eHalfLengthZcutout))));
   auto posInDisc =
-      Surface::makeShared<DiscSurface>(posInDiscTrf, m_rmin, m_rmed);
+      Surface::makeShared<DiscSurface>(posInDiscTrf, get(eMinR), get(eMedR));
   surfaces.at(index5) = posInDisc;
 
   auto negInDiscTrf = std::make_shared<const Transform3D>(
-      *trf * Translation3D(Vector3D(0, 0, -m_dz2)) *
+      *trf * Translation3D(Vector3D(0, 0, -get(eHalfLengthZcutout))) *
       AngleAxis3D(M_PI, Vector3D::UnitX()));
   auto negInDisc =
-      Surface::makeShared<DiscSurface>(negInDiscTrf, m_rmin, m_rmed);
+      Surface::makeShared<DiscSurface>(negInDiscTrf, get(eMinR), get(eMedR));
   surfaces.at(index4) = negInDisc;
 
   return surfaces;
@@ -131,8 +133,8 @@ Acts::Volume::BoundingBox Acts::CutoutCylinderVolumeBounds::boundingBox(
   // no phi sector is possible, so this is just the outer size of
   // the cylinder
 
-  vmax = {m_rmax, m_rmax, m_dz1};
-  vmin = {-m_rmax, -m_rmax, -m_dz1};
+  vmax = {get(eMaxR), get(eMaxR), get(eHalfLengthZ)};
+  vmin = {-get(eMaxR), -get(eMaxR), -get(eHalfLengthZ)};
 
   Acts::Volume::BoundingBox box(entity, vmin - envelope, vmax + envelope);
   // transform at the very end, if required
@@ -142,8 +144,8 @@ Acts::Volume::BoundingBox Acts::CutoutCylinderVolumeBounds::boundingBox(
 std::ostream& Acts::CutoutCylinderVolumeBounds::toStream(
     std::ostream& sl) const {
   sl << "Acts::CutoutCylinderVolumeBounds(\n";
-  sl << "rmin = " << m_rmin << " rmed = " << m_rmed << " rmax = " << m_rmax
-     << "\n";
-  sl << "dz1 = " << m_dz1 << " dz2 = " << m_dz2;
+  sl << "rmin = " << get(eMinR) << " rmed = " << get(eMedR)
+     << " rmax = " << get(eMaxR) << "\n";
+  sl << "dz1 = " << get(eHalfLengthZ) << " dz2 = " << get(eHalfLengthZcutout);
   return sl;
 }
diff --git a/Core/src/Geometry/CylinderVolumeBounds.cpp b/Core/src/Geometry/CylinderVolumeBounds.cpp
index 6e0dd1c2287405443c9cae28a029557de968f206..c1548b547e00a97f79fd36ee3c580d51941c58c0 100644
--- a/Core/src/Geometry/CylinderVolumeBounds.cpp
+++ b/Core/src/Geometry/CylinderVolumeBounds.cpp
@@ -11,8 +11,6 @@
 ///////////////////////////////////////////////////////////////////
 
 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
-#include <cmath>
-#include <iostream>
 #include "Acts/Surfaces/CylinderBounds.hpp"
 #include "Acts/Surfaces/CylinderSurface.hpp"
 #include "Acts/Surfaces/DiscSurface.hpp"
@@ -22,68 +20,38 @@
 #include "Acts/Utilities/BoundingBox.hpp"
 #include "Acts/Utilities/IVisualization.hpp"
 
-const double Acts::CylinderVolumeBounds::s_numericalStable = 10e-2;
-
-Acts::CylinderVolumeBounds::CylinderVolumeBounds()
-    : VolumeBounds(), m_values(4, 0.) {}
-
-Acts::CylinderVolumeBounds::CylinderVolumeBounds(double radius, double halez)
-    : VolumeBounds(), m_values(4, 0.) {
-  m_values.at(bv_innerRadius) = 0.;
-  m_values.at(bv_outerRadius) = std::abs(radius);
-  m_values.at(bv_halfPhiSector) = M_PI;
-  m_values.at(bv_halfZ) = std::abs(halez);
-}
-
-Acts::CylinderVolumeBounds::CylinderVolumeBounds(double rinner, double router,
-                                                 double halez)
-    : VolumeBounds(), m_values(4, 0.) {
-  m_values.at(bv_innerRadius) = std::abs(rinner);
-  m_values.at(bv_outerRadius) = std::abs(router);
-  m_values.at(bv_halfPhiSector) = M_PI;
-  m_values.at(bv_halfZ) = std::abs(halez);
-}
-
-Acts::CylinderVolumeBounds::CylinderVolumeBounds(double rinner, double router,
-                                                 double haphi, double halez)
-    : VolumeBounds(), m_values(4, 0.) {
-  m_values.at(bv_innerRadius) = std::abs(rinner);
-  m_values.at(bv_outerRadius) = std::abs(router);
-  m_values.at(bv_halfPhiSector) = std::abs(haphi);
-  m_values.at(bv_halfZ) = std::abs(halez);
-}
+#include <cmath>
+#include <iostream>
 
-Acts::CylinderVolumeBounds::CylinderVolumeBounds(const CylinderBounds& cBounds,
-                                                 double thickness)
-    : VolumeBounds(), m_values(4, 0.) {
+Acts::CylinderVolumeBounds::CylinderVolumeBounds(
+    const CylinderBounds& cBounds, double thickness) noexcept(false)
+    : VolumeBounds() {
   double cR = cBounds.get(CylinderBounds::eR);
-  m_values.at(bv_innerRadius) = cR - 0.5 * thickness;
-  m_values.at(bv_outerRadius) = cR + 0.5 * thickness;
-  m_values.at(bv_halfPhiSector) = cBounds.get(CylinderBounds::eHalfPhiSector);
-  m_values.at(bv_halfZ) = cBounds.get(CylinderBounds::eHalfLengthZ);
-}
-
-Acts::CylinderVolumeBounds::CylinderVolumeBounds(const RadialBounds& rBounds,
-                                                 double thickness)
-    : VolumeBounds(), m_values(4, 0.) {
-  m_values.at(bv_innerRadius) = rBounds.get(RadialBounds::eMinR);
-  m_values.at(bv_outerRadius) = rBounds.get(RadialBounds::eMaxR);
-  m_values.at(bv_halfPhiSector) = rBounds.get(RadialBounds::eHalfPhiSector);
-  m_values.at(bv_halfZ) = 0.5 * thickness;
+  if (thickness <= 0. or (cR - 0.5 * thickness) < 0.) {
+    throw(std::invalid_argument(
+        "CylinderVolumeBounds: invalid extrusion thickness."));
+  }
+  m_values[eMinR] = cR - 0.5 * thickness;
+  m_values[eMaxR] = cR + 0.5 * thickness;
+  m_values[eHalfLengthZ] = cBounds.get(CylinderBounds::eHalfLengthZ);
+  m_values[eHalfPhiSector] = cBounds.get(CylinderBounds::eHalfPhiSector);
+  m_values[eAveragePhi] = cBounds.get(CylinderBounds::eAveragePhi);
+  buildSurfaceBounds();
 }
 
 Acts::CylinderVolumeBounds::CylinderVolumeBounds(
-    const CylinderVolumeBounds& cylbo)
-    : VolumeBounds(), m_values(cylbo.m_values) {}
-
-Acts::CylinderVolumeBounds::~CylinderVolumeBounds() = default;
-
-Acts::CylinderVolumeBounds& Acts::CylinderVolumeBounds::operator=(
-    const CylinderVolumeBounds& cylbo) {
-  if (this != &cylbo) {
-    m_values = cylbo.m_values;
+    const RadialBounds& rBounds, double thickness) noexcept(false)
+    : VolumeBounds() {
+  if (thickness <= 0.) {
+    throw(std::invalid_argument(
+        "CylinderVolumeBounds: invalid extrusion thickness."));
   }
-  return *this;
+  m_values[eMinR] = rBounds.get(RadialBounds::eMinR);
+  m_values[eMaxR] = rBounds.get(RadialBounds::eMaxR);
+  m_values[eHalfLengthZ] = 0.5 * thickness;
+  m_values[eHalfPhiSector] = rBounds.get(RadialBounds::eHalfPhiSector);
+  m_values[eAveragePhi] = rBounds.get(RadialBounds::eAveragePhi);
+  buildSurfaceBounds();
 }
 
 std::vector<std::shared_ptr<const Acts::Surface>>
@@ -100,76 +68,62 @@ Acts::CylinderVolumeBounds::decomposeToSurfaces(
   const Transform3D* tTransform = nullptr;
   Vector3D cylCenter(transform.translation());
 
-  std::shared_ptr<const DiscBounds> dBounds = discBounds();
   // bottom Disc (negative z)
   tTransform =
       new Transform3D(transform * AngleAxis3D(M_PI, Vector3D(1., 0., 0.)) *
-                      Translation3D(Vector3D(0., 0., halflengthZ())));
+                      Translation3D(Vector3D(0., 0., get(eHalfLengthZ))));
   rSurfaces.push_back(Surface::makeShared<DiscSurface>(
-      std::shared_ptr<const Transform3D>(tTransform), dBounds));
+      std::shared_ptr<const Transform3D>(tTransform), m_discBounds));
   // top Disc (positive z)
-  tTransform = new Transform3D(transform *
-                               Translation3D(Vector3D(0., 0., halflengthZ())));
+  tTransform = new Transform3D(
+      transform * Translation3D(Vector3D(0., 0., get(eHalfLengthZ))));
   rSurfaces.push_back(Surface::makeShared<DiscSurface>(
-      std::shared_ptr<const Transform3D>(tTransform), dBounds));
+      std::shared_ptr<const Transform3D>(tTransform), m_discBounds));
 
   // outer Cylinder - shares the transform
   rSurfaces.push_back(
-      Surface::makeShared<CylinderSurface>(trfShared, outerCylinderBounds()));
+      Surface::makeShared<CylinderSurface>(trfShared, m_outerCylinderBounds));
 
   // innermost Cylinder
-  if (innerRadius() > s_numericalStable) {
+  if (m_innerCylinderBounds != nullptr) {
     rSurfaces.push_back(
-        Surface::makeShared<CylinderSurface>(trfShared, innerCylinderBounds()));
+        Surface::makeShared<CylinderSurface>(trfShared, m_innerCylinderBounds));
   }
 
   // the cylinder is sectoral
-  if (std::abs(halfPhiSector() - M_PI) > s_numericalStable) {
-    std::shared_ptr<const PlanarBounds> sp12Bounds = sectorPlaneBounds();
+  if (m_sectorPlaneBounds != nullptr) {
     // sectorPlane 1 (negative phi)
     const Transform3D* sp1Transform = new Transform3D(
-        transform * AngleAxis3D(-halfPhiSector(), Vector3D(0., 0., 1.)) *
-        Translation3D(Vector3D(mediumRadius(), 0., 0.)) *
+        transform * AngleAxis3D(-get(eHalfPhiSector), Vector3D(0., 0., 1.)) *
+        Translation3D(Vector3D(0.5 * (get(eMinR) + get(eMaxR)), 0., 0.)) *
         AngleAxis3D(M_PI / 2, Vector3D(1., 0., 0.)));
     rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-        std::shared_ptr<const Transform3D>(sp1Transform), sp12Bounds));
+        std::shared_ptr<const Transform3D>(sp1Transform), m_sectorPlaneBounds));
     // sectorPlane 2 (positive phi)
     const Transform3D* sp2Transform = new Transform3D(
-        transform * AngleAxis3D(halfPhiSector(), Vector3D(0., 0., 1.)) *
-        Translation3D(Vector3D(mediumRadius(), 0., 0.)) *
+        transform * AngleAxis3D(get(eHalfPhiSector), Vector3D(0., 0., 1.)) *
+        Translation3D(Vector3D(0.5 * (get(eMinR) + get(eMaxR)), 0., 0.)) *
         AngleAxis3D(-M_PI / 2, Vector3D(1., 0., 0.)));
     rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-        std::shared_ptr<const Transform3D>(sp2Transform), sp12Bounds));
+        std::shared_ptr<const Transform3D>(sp2Transform), m_sectorPlaneBounds));
   }
   return rSurfaces;
 }
 
-std::shared_ptr<const Acts::CylinderBounds>
-Acts::CylinderVolumeBounds::innerCylinderBounds() const {
-  return std::make_shared<const CylinderBounds>(m_values.at(bv_innerRadius),
-                                                m_values.at(bv_halfZ),
-                                                m_values.at(bv_halfPhiSector));
-}
-
-std::shared_ptr<const Acts::CylinderBounds>
-Acts::CylinderVolumeBounds::outerCylinderBounds() const {
-  return std::make_shared<const CylinderBounds>(m_values.at(bv_outerRadius),
-                                                m_values.at(bv_halfZ),
-                                                m_values.at(bv_halfPhiSector));
-}
-
-std::shared_ptr<const Acts::DiscBounds> Acts::CylinderVolumeBounds::discBounds()
-    const {
-  return std::shared_ptr<const DiscBounds>(
-      new RadialBounds(m_values.at(bv_innerRadius), m_values.at(bv_outerRadius),
-                       m_values.at(bv_halfPhiSector)));
-}
-
-std::shared_ptr<const Acts::PlanarBounds>
-Acts::CylinderVolumeBounds::sectorPlaneBounds() const {
-  return std::shared_ptr<const PlanarBounds>(new RectangleBounds(
-      0.5 * (m_values.at(bv_outerRadius) - m_values.at(bv_innerRadius)),
-      m_values.at(bv_halfZ)));
+void Acts::CylinderVolumeBounds::buildSurfaceBounds() {
+  if (get(eMinR) > s_epsilon) {
+    m_innerCylinderBounds = std::make_shared<const CylinderBounds>(
+        get(eMinR), get(eHalfLengthZ), get(eHalfPhiSector), get(eAveragePhi));
+  }
+  m_outerCylinderBounds = std::make_shared<const CylinderBounds>(
+      get(eMaxR), get(eHalfLengthZ), get(eHalfPhiSector), get(eAveragePhi));
+  m_discBounds = std::make_shared<const RadialBounds>(
+      get(eMinR), get(eMaxR), get(eHalfPhiSector), get(eAveragePhi));
+
+  if (std::abs(get(eHalfPhiSector) - M_PI) > s_epsilon) {
+    m_sectorPlaneBounds = std::make_shared<const RectangleBounds>(
+        0.5 * (get(eMaxR) - get(eMinR)), get(eHalfLengthZ));
+  }
 }
 
 std::ostream& Acts::CylinderVolumeBounds::toStream(std::ostream& sl) const {
@@ -180,23 +134,23 @@ Acts::Volume::BoundingBox Acts::CylinderVolumeBounds::boundingBox(
     const Transform3D* trf, const Vector3D& envelope,
     const Volume* entity) const {
   double xmax, xmin, ymax, ymin;
-  xmax = outerRadius();
+  xmax = get(eMaxR);
 
-  if (halfPhiSector() > M_PI / 2.) {
+  if (get(eHalfPhiSector) > M_PI / 2.) {
     // more than half
-    ymax = outerRadius();
-    ymin = -outerRadius();
-    xmin = outerRadius() * std::cos(halfPhiSector());
+    ymax = xmax;
+    ymin = -xmax;
+    xmin = xmax * std::cos(get(eHalfPhiSector));
   } else {
     // less than half
-    ymax = outerRadius() * std::sin(halfPhiSector());
+    ymax = get(eMaxR) * std::sin(get(eHalfPhiSector));
     ymin = -ymax;
     // in this case, xmin is given by the inner radius
-    xmin = innerRadius() * std::cos(halfPhiSector());
+    xmin = get(eMinR) * std::cos(get(eHalfPhiSector));
   }
 
-  Vector3D vmin(xmin, ymin, -halflengthZ());
-  Vector3D vmax(xmax, ymax, halflengthZ());
+  Vector3D vmin(xmin, ymin, -get(eHalfLengthZ));
+  Vector3D vmax(xmax, ymax, get(eHalfLengthZ));
 
   // this is probably not perfect, but at least conservative
   Volume::BoundingBox box{entity, vmin - envelope, vmax + envelope};
diff --git a/Core/src/Geometry/CylinderVolumeBuilder.cpp b/Core/src/Geometry/CylinderVolumeBuilder.cpp
index e7c017bf633a944d2a3bf32a12e28ef53c0c821c..5aa6dbe01bdd38cb3b9163405fc3b041a45bb7f9 100644
--- a/Core/src/Geometry/CylinderVolumeBuilder.cpp
+++ b/Core/src/Geometry/CylinderVolumeBuilder.cpp
@@ -91,12 +91,16 @@ Acts::CylinderVolumeBuilder::trackingVolume(
         &existingVolume->volumeBounds());
     // set the inside values
     wConfig.existingVolumeConfig.present = true;
-    wConfig.existingVolumeConfig.rMin = existingBounds->innerRadius();
-    wConfig.existingVolumeConfig.rMax = existingBounds->outerRadius();
+    wConfig.existingVolumeConfig.rMin =
+        existingBounds->get(CylinderVolumeBounds::eMinR);
+    wConfig.existingVolumeConfig.rMax =
+        existingBounds->get(CylinderVolumeBounds::eMaxR);
     wConfig.existingVolumeConfig.zMin =
-        existingVolume->center().z() - existingBounds->halflengthZ();
+        existingVolume->center().z() -
+        existingBounds->get(CylinderVolumeBounds::eHalfLengthZ);
     wConfig.existingVolumeConfig.zMax =
-        existingVolume->center().z() + existingBounds->halflengthZ();
+        existingVolume->center().z() +
+        existingBounds->get(CylinderVolumeBounds::eHalfLengthZ);
   }
   //
   // b) outside config
@@ -109,10 +113,14 @@ Acts::CylinderVolumeBuilder::trackingVolume(
     if (ocvBounds != nullptr) {
       // get values from the out bounds
       wConfig.externalVolumeConfig.present = true;
-      wConfig.externalVolumeConfig.rMin = ocvBounds->innerRadius();
-      wConfig.externalVolumeConfig.rMax = ocvBounds->outerRadius();
-      wConfig.externalVolumeConfig.zMin = -ocvBounds->halflengthZ();
-      wConfig.externalVolumeConfig.zMax = ocvBounds->halflengthZ();
+      wConfig.externalVolumeConfig.rMin =
+          ocvBounds->get(CylinderVolumeBounds::eMinR);
+      wConfig.externalVolumeConfig.rMax =
+          ocvBounds->get(CylinderVolumeBounds::eMaxR);
+      wConfig.externalVolumeConfig.zMin =
+          -ocvBounds->get(CylinderVolumeBounds::eHalfLengthZ);
+      wConfig.externalVolumeConfig.zMax =
+          ocvBounds->get(CylinderVolumeBounds::eHalfLengthZ);
     }
   }
 
@@ -528,10 +536,12 @@ Acts::VolumeConfig Acts::CylinderVolumeBuilder::analyzeContent(
       const CylinderVolumeBounds* cvBounds =
           dynamic_cast<const CylinderVolumeBounds*>(&volume->volumeBounds());
       if (cvBounds != nullptr) {
-        takeSmaller(lConfig.rMin, cvBounds->innerRadius());
-        takeBigger(lConfig.rMax, cvBounds->outerRadius());
-        takeSmaller(lConfig.zMin, -cvBounds->halflengthZ());
-        takeBigger(lConfig.zMax, cvBounds->halflengthZ());
+        takeSmaller(lConfig.rMin, cvBounds->get(CylinderVolumeBounds::eMinR));
+        takeBigger(lConfig.rMax, cvBounds->get(CylinderVolumeBounds::eMaxR));
+        takeSmaller(lConfig.zMin,
+                    -cvBounds->get(CylinderVolumeBounds::eHalfLengthZ));
+        takeBigger(lConfig.zMax,
+                   cvBounds->get(CylinderVolumeBounds::eHalfLengthZ));
       }
     }
   }
diff --git a/Core/src/Geometry/CylinderVolumeHelper.cpp b/Core/src/Geometry/CylinderVolumeHelper.cpp
index c0e14aa79bf1631403ea3ef403bee9ec3da03c0f..65df57e28e77105f57d1441c5f6cdb2a74f49ba7 100644
--- a/Core/src/Geometry/CylinderVolumeHelper.cpp
+++ b/Core/src/Geometry/CylinderVolumeHelper.cpp
@@ -103,15 +103,20 @@ Acts::CylinderVolumeHelper::createTrackingVolume(
     // get the zMin/Max
     double zMin =
         (transform ? transform->translation().z() : 0.) +
-        (cylinderBounds != nullptr ? -cylinderBounds->halflengthZ() : 0.);
-    double zMax =
-        (transform ? transform->translation().z() : 0.) +
-        (cylinderBounds != nullptr ? cylinderBounds->halflengthZ() : 0.);
+        (cylinderBounds != nullptr
+             ? -cylinderBounds->get(CylinderVolumeBounds::eHalfLengthZ)
+             : 0.);
+    double zMax = (transform ? transform->translation().z() : 0.) +
+                  (cylinderBounds != nullptr
+                       ? cylinderBounds->get(CylinderVolumeBounds::eHalfLengthZ)
+                       : 0.);
     // get the rMin/rmAx
-    double rMin =
-        cylinderBounds != nullptr ? cylinderBounds->innerRadius() : rMinRaw;
-    double rMax =
-        cylinderBounds != nullptr ? cylinderBounds->outerRadius() : rMaxRaw;
+    double rMin = cylinderBounds != nullptr
+                      ? cylinderBounds->get(CylinderVolumeBounds::eMinR)
+                      : rMinRaw;
+    double rMax = cylinderBounds != nullptr
+                      ? cylinderBounds->get(CylinderVolumeBounds::eMaxR)
+                      : rMaxRaw;
 
     ACTS_VERBOSE(
         "Filling the layers into an appropriate layer array - with "
@@ -173,8 +178,8 @@ Acts::CylinderVolumeHelper::createTrackingVolume(
   zPosition = std::abs(zPosition) < 0.1 ? 0. : zPosition;
 
   // now create the cylinder volume bounds
-  cBounds = rMin > 0.1 ? new CylinderVolumeBounds(rMin, rMax, halflengthZ)
-                       : new CylinderVolumeBounds(rMax, halflengthZ);
+  cBounds = new CylinderVolumeBounds(rMin, rMax, halflengthZ);
+
   // transform
   std::shared_ptr<const Transform3D> transform =
       (zPosition != 0) ? std::make_shared<const Transform3D>(
@@ -324,8 +329,9 @@ Acts::CylinderVolumeHelper::createContainerTrackingVolume(
     return nullptr;
   }
   // Check whether it is a r-binned case or a z-binned case
-  bool rCase = std::abs(firstVolumeBounds->innerRadius() -
-                        lastVolumeBounds->innerRadius()) > 0.1;
+  bool rCase =
+      std::abs(firstVolumeBounds->get(CylinderVolumeBounds::eMinR) -
+               lastVolumeBounds->get(CylinderVolumeBounds::eMinR)) > 0.1;
 
   // Fill these ones depending on the rCase though assignment
   double zMin = 0.;
@@ -336,20 +342,25 @@ Acts::CylinderVolumeHelper::createContainerTrackingVolume(
   double zSep1 = 0.;
   double zSep2 = 0.;
   if (rCase) {
-    zMin = (*firstVolume)->center().z() - firstVolumeBounds->halflengthZ();
-    zMax = (*firstVolume)->center().z() + firstVolumeBounds->halflengthZ();
+    zMin = (*firstVolume)->center().z() -
+           firstVolumeBounds->get(CylinderVolumeBounds::eHalfLengthZ);
+    zMax = (*firstVolume)->center().z() +
+           firstVolumeBounds->get(CylinderVolumeBounds::eHalfLengthZ);
     zSep1 = zMin;
     zSep2 = zMax;
-    rMin = firstVolumeBounds->innerRadius();
-    rGlueMin = firstVolumeBounds->outerRadius();
-    rMax = lastVolumeBounds->outerRadius();
+    rMin = firstVolumeBounds->get(CylinderVolumeBounds::eMinR);
+    rGlueMin = firstVolumeBounds->get(CylinderVolumeBounds::eMaxR);
+    rMax = lastVolumeBounds->get(CylinderVolumeBounds::eMaxR);
   } else {
-    zMin = (*firstVolume)->center().z() - firstVolumeBounds->halflengthZ();
-    zMax = (*lastVolume)->center().z() + lastVolumeBounds->halflengthZ();
-    zSep1 = (*firstVolume)->center().z() + firstVolumeBounds->halflengthZ();
+    zMin = (*firstVolume)->center().z() -
+           firstVolumeBounds->get(CylinderVolumeBounds::eHalfLengthZ);
+    zMax = (*lastVolume)->center().z() +
+           lastVolumeBounds->get(CylinderVolumeBounds::eHalfLengthZ);
+    zSep1 = (*firstVolume)->center().z() +
+            firstVolumeBounds->get(CylinderVolumeBounds::eHalfLengthZ);
     zSep2 = zSep1;
-    rMin = firstVolumeBounds->innerRadius();
-    rMax = firstVolumeBounds->outerRadius();
+    rMin = firstVolumeBounds->get(CylinderVolumeBounds::eMinR);
+    rMax = firstVolumeBounds->get(CylinderVolumeBounds::eMaxR);
   }
   // Estimate the z - position
   double zPos = 0.5 * (zMin + zMax);
@@ -360,9 +371,7 @@ Acts::CylinderVolumeHelper::createContainerTrackingVolume(
           : nullptr;
   // Create the bounds from the information gathered so far
   CylinderVolumeBounds* topVolumeBounds =
-      std::abs(rMin) > 0.1
-          ? new CylinderVolumeBounds(rMin, rMax, 0.5 * std::abs(zMax - zMin))
-          : new CylinderVolumeBounds(rMax, 0.5 * std::abs(zMax - zMin));
+      new CylinderVolumeBounds(rMin, rMax, 0.5 * std::abs(zMax - zMin));
 
   // some screen output
   ACTS_VERBOSE("Container volume bounds are " << (*topVolumeBounds));
@@ -413,9 +422,11 @@ bool Acts::CylinderVolumeHelper::estimateAndCheckDimension(
                               << " layers to gather overall dimensions");
   if (cylinderVolumeBounds != nullptr)
     ACTS_DEBUG("Cylinder volume bounds are given: (rmin/rmax/dz) = "
-               << "(" << cylinderVolumeBounds->innerRadius() << "/"
-               << cylinderVolumeBounds->outerRadius() << "/"
-               << cylinderVolumeBounds->halflengthZ() << ")");
+               << "(" << cylinderVolumeBounds->get(CylinderVolumeBounds::eMinR)
+               << "/" << cylinderVolumeBounds->get(CylinderVolumeBounds::eMaxR)
+               << "/"
+               << cylinderVolumeBounds->get(CylinderVolumeBounds::eHalfLengthZ)
+               << ")");
 
   // prepare for parsing the layers
   double layerRmin = 10e10;
@@ -509,34 +520,46 @@ bool Acts::CylinderVolumeHelper::estimateAndCheckDimension(
                << layerZmax);
 
   double zFromTransform = transform ? transform->translation().z() : 0.;
-  ACTS_VERBOSE("    -> while created bounds are (rMin/rMax/zMin/zMax) = "
-               << cylinderVolumeBounds->innerRadius() << " / "
-               << cylinderVolumeBounds->outerRadius() << " / "
-               << zFromTransform - cylinderVolumeBounds->halflengthZ() << " / "
-               << zFromTransform + cylinderVolumeBounds->halflengthZ());
+  ACTS_VERBOSE(
+      "    -> while created bounds are (rMin/rMax/zMin/zMax) = "
+      << cylinderVolumeBounds->get(CylinderVolumeBounds::eMinR) << " / "
+      << cylinderVolumeBounds->get(CylinderVolumeBounds::eMaxR) << " / "
+      << zFromTransform -
+             cylinderVolumeBounds->get(CylinderVolumeBounds::eHalfLengthZ)
+      << " / "
+      << zFromTransform +
+             cylinderVolumeBounds->get(CylinderVolumeBounds::eHalfLengthZ));
 
   // both is NOW given --- check it -----------------------------
   if (cylinderVolumeBounds != nullptr) {
     // only check
-    if (zFromTransform - cylinderVolumeBounds->halflengthZ() <= layerZmin &&
-        zFromTransform + cylinderVolumeBounds->halflengthZ() >= layerZmax &&
-        cylinderVolumeBounds->innerRadius() <= layerRmin &&
-        cylinderVolumeBounds->outerRadius() >= layerRmax) {
+    if (zFromTransform -
+                cylinderVolumeBounds->get(CylinderVolumeBounds::eHalfLengthZ) <=
+            layerZmin &&
+        zFromTransform +
+                cylinderVolumeBounds->get(CylinderVolumeBounds::eHalfLengthZ) >=
+            layerZmax &&
+        cylinderVolumeBounds->get(CylinderVolumeBounds::eMinR) <= layerRmin &&
+        cylinderVolumeBounds->get(CylinderVolumeBounds::eMaxR) >= layerRmax) {
       return true;
     } else {
       ACTS_WARNING(
           "Provided layers are not contained by volume ! Bailing out. ");
       ACTS_WARNING("- zFromTransform: " << zFromTransform);
       ACTS_WARNING("- volumeZmin:"
-                   << zFromTransform - cylinderVolumeBounds->halflengthZ()
+                   << zFromTransform - cylinderVolumeBounds->get(
+                                           CylinderVolumeBounds::eHalfLengthZ)
                    << ", layerZmin: " << layerZmin);
       ACTS_WARNING("- volumeZmax: "
-                   << zFromTransform + cylinderVolumeBounds->halflengthZ()
+                   << zFromTransform + cylinderVolumeBounds->get(
+                                           CylinderVolumeBounds::eHalfLengthZ)
                    << ", layerZmax: " << layerZmax);
-      ACTS_WARNING("- volumeRmin: " << cylinderVolumeBounds->innerRadius()
-                                    << ", layerRmin: " << layerRmin);
-      ACTS_WARNING("- volumeRmax: " << cylinderVolumeBounds->outerRadius()
-                                    << ", layerRmax: " << layerRmax);
+      ACTS_WARNING("- volumeRmin: "
+                   << cylinderVolumeBounds->get(CylinderVolumeBounds::eMinR)
+                   << ", layerRmin: " << layerRmin);
+      ACTS_WARNING("- volumeRmax: "
+                   << cylinderVolumeBounds->get(CylinderVolumeBounds::eMaxR)
+                   << ", layerRmax: " << layerRmax);
       return false;
     }
   }
@@ -809,7 +832,7 @@ void Acts::CylinderVolumeHelper::glueTrackingVolumes(
       auto cylVolBounds = dynamic_cast<const Acts::CylinderVolumeBounds*>(
           &tvolOne->volumeBounds());
       double zPos = tvolOne->center().z();
-      double zHL = cylVolBounds->halflengthZ();
+      double zHL = cylVolBounds->get(CylinderVolumeBounds::eHalfLengthZ);
       transform =
           std::make_shared<const Transform3D>(Translation3D(0, 0, zPos + zHL));
       // this puts the surface on the positive z side of the cyl vol bounds
diff --git a/Core/src/Geometry/DoubleTrapezoidVolumeBounds.cpp b/Core/src/Geometry/DoubleTrapezoidVolumeBounds.cpp
deleted file mode 100644
index 9df847635b2aab3cd87e89b4f98bddb0adf0995a..0000000000000000000000000000000000000000
--- a/Core/src/Geometry/DoubleTrapezoidVolumeBounds.cpp
+++ /dev/null
@@ -1,287 +0,0 @@
-// This file is part of the Acts project.
-//
-// Copyright (C) 2016-2018 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/.
-
-///////////////////////////////////////////////////////////////////
-// DoubleTrapezoidVolumeBounds.cpp, Acts project
-///////////////////////////////////////////////////////////////////
-
-#include "Acts/Geometry/DoubleTrapezoidVolumeBounds.hpp"
-
-#include <cmath>
-#include <iomanip>
-#include <iostream>
-#include "Acts/Surfaces/DiamondBounds.hpp"
-#include "Acts/Surfaces/PlaneSurface.hpp"
-#include "Acts/Surfaces/RectangleBounds.hpp"
-
-Acts::DoubleTrapezoidVolumeBounds::DoubleTrapezoidVolumeBounds()
-    : VolumeBounds(), m_values(bv_length, 0.) {}
-
-Acts::DoubleTrapezoidVolumeBounds::DoubleTrapezoidVolumeBounds(
-    double minhalex, double medhalex, double maxhalex, double haley1,
-    double haley2, double halez)
-    : VolumeBounds(), m_values(bv_length, 0.) {
-  m_values.at(bv_minHalfX) = minhalex;
-  m_values.at(bv_medHalfX) = medhalex;
-  m_values.at(bv_maxHalfX) = maxhalex;
-  m_values.at(bv_halfY1) = haley1;
-  m_values.at(bv_halfY2) = haley2;
-  m_values.at(bv_halfZ) = halez;
-  m_values.at(bv_alpha1) =
-      atan2(m_values.at(bv_medHalfX) - m_values.at(bv_minHalfX),
-            2. * m_values.at(bv_halfY1));
-  m_values.at(bv_alpha2) =
-      atan2(m_values.at(bv_medHalfX) - m_values.at(bv_maxHalfX),
-            2. * m_values.at(bv_halfY2));
-}
-
-Acts::DoubleTrapezoidVolumeBounds::DoubleTrapezoidVolumeBounds(
-    const Acts::DoubleTrapezoidVolumeBounds& trabo)
-    : VolumeBounds(), m_values(trabo.m_values) {}
-
-Acts::DoubleTrapezoidVolumeBounds::~DoubleTrapezoidVolumeBounds() = default;
-
-Acts::DoubleTrapezoidVolumeBounds& Acts::DoubleTrapezoidVolumeBounds::operator=(
-    const Acts::DoubleTrapezoidVolumeBounds& trabo) {
-  if (this != &trabo) {
-    m_values = trabo.m_values;
-  }
-  return *this;
-}
-
-std::vector<std::shared_ptr<const Acts::Surface>>
-Acts::DoubleTrapezoidVolumeBounds::decomposeToSurfaces(
-    const Transform3D* transformPtr) const {
-  std::vector<std::shared_ptr<const Surface>> rSurfaces;
-
-  // the transform
-  Transform3D transform =
-      (transformPtr == nullptr) ? Transform3D::Identity() : (*transformPtr);
-
-  // face surfaces xy
-  RotationMatrix3D diamondRotation(transform.rotation());
-  Vector3D diamondX(diamondRotation.col(0));
-  Vector3D diamondY(diamondRotation.col(1));
-  Vector3D diamondZ(diamondRotation.col(2));
-  Vector3D diamondCenter(transform.translation());
-
-  const Transform3D* tTransform = nullptr;
-
-  //   (1) - at negative local z
-  tTransform =
-      new Transform3D(transform * AngleAxis3D(M_PI, Vector3D(0., 1., 0.)) *
-                      Translation3D(Vector3D(0., 0., halflengthZ())));
-  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-      std::shared_ptr<const Transform3D>(tTransform),
-      std::shared_ptr<const PlanarBounds>(faceXYDiamondBounds())));
-  //   (2) - at positive local z
-  tTransform = new Transform3D(transform *
-                               Translation3D(Vector3D(0., 0., halflengthZ())));
-  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-      std::shared_ptr<const Transform3D>(tTransform),
-      std::shared_ptr<const PlanarBounds>(faceXYDiamondBounds())));
-  // face surfaces yz
-  // transmute cyclical
-  //   (3) - at point A, attached to alpha opening angle
-  Vector3D A(diamondCenter - minHalflengthX() * diamondX -
-             2 * halflengthY1() * diamondY);
-  AngleAxis3D alpha1ZRotation(alpha1(), Vector3D(0., 0., 1.));
-  RotationMatrix3D alpha1Rotation(
-      diamondRotation * alpha1ZRotation *
-      AngleAxis3D(-0.5 * M_PI, Vector3D(0., 1., 0.)) *
-      AngleAxis3D(0.5 * M_PI, Vector3D(0., 0., 1.)));
-  RectangleBounds* faceAlpha1Bounds = faceAlpha1RectangleBounds();
-  Vector3D faceAlpha1Position(A + alpha1Rotation.col(0) *
-                                      faceAlpha1Bounds->halfLengthX());
-  tTransform =
-      new Transform3D(Translation3D(faceAlpha1Position) * alpha1Rotation);
-  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-      std::shared_ptr<const Transform3D>(tTransform),
-      std::shared_ptr<const PlanarBounds>(faceAlpha1Bounds)));
-  //   (4) - at point B, attached to beta opening angle
-  Vector3D B(diamondCenter + minHalflengthX() * diamondX -
-             2 * halflengthY1() * diamondY);
-  AngleAxis3D beta1ZRotation(-alpha1(), Vector3D(0., 0., 1.));
-  RotationMatrix3D beta1Rotation(diamondRotation * beta1ZRotation *
-                                 AngleAxis3D(0.5 * M_PI, Vector3D(0., 1., 0.)) *
-                                 AngleAxis3D(0.5 * M_PI, Vector3D(0., 0., 1.)));
-  RectangleBounds* faceBeta1Bounds = faceBeta1RectangleBounds();
-  Vector3D faceBeta1Position(B + beta1Rotation.col(0) *
-                                     faceBeta1Bounds->halfLengthX());
-  tTransform =
-      new Transform3D(Translation3D(faceBeta1Position) * beta1Rotation);
-  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-      std::shared_ptr<const Transform3D>(tTransform),
-      std::shared_ptr<const PlanarBounds>(faceBeta1Bounds)));
-  // face surfaces yz
-  // transmute cyclical
-  //   (5) - at point A', attached to alpha opening angle
-  Vector3D AA(diamondCenter - maxHalflengthX() * diamondX +
-              2 * halflengthY2() * diamondY);
-  AngleAxis3D alpha2ZRotation(-alpha2(), Vector3D(0., 0., 1.));
-  RotationMatrix3D alpha2Rotation(
-      diamondRotation * alpha2ZRotation *
-      AngleAxis3D(-0.5 * M_PI, Vector3D(0., 1., 0.)) *
-      AngleAxis3D(-0.5 * M_PI, Vector3D(0., 0., 1.)));
-  RectangleBounds* faceAlpha2Bounds = faceAlpha2RectangleBounds();
-  Vector3D faceAlpha2Position(AA + alpha2Rotation.col(0) *
-                                       faceAlpha2Bounds->halfLengthX());
-  tTransform =
-      new Transform3D(Translation3D(faceAlpha2Position) * alpha2Rotation);
-  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-      std::shared_ptr<const Transform3D>(tTransform),
-      std::shared_ptr<const PlanarBounds>(faceAlpha2Bounds)));
-  //   (6) - at point B', attached to beta opening angle
-  Vector3D BB(diamondCenter + maxHalflengthX() * diamondX +
-              2 * halflengthY2() * diamondY);
-  AngleAxis3D beta2ZRotation(alpha2(), Vector3D(0., 0., 1.));
-  RotationMatrix3D beta2Rotation(
-      diamondRotation * beta2ZRotation *
-      AngleAxis3D(0.5 * M_PI, Vector3D(0., 1., 0.)) *
-      AngleAxis3D(-0.5 * M_PI, Vector3D(0., 0., 1.)));
-  RectangleBounds* faceBeta2Bounds = faceBeta2RectangleBounds();
-  Vector3D faceBeta2Position(BB + beta2Rotation.col(0) *
-                                      faceBeta2Bounds->halfLengthX());
-  tTransform =
-      new Transform3D(Translation3D(faceBeta2Position) * beta2Rotation);
-  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-      std::shared_ptr<const Transform3D>(tTransform),
-      std::shared_ptr<const PlanarBounds>(faceBeta2Bounds)));
-  // face surfaces zx
-  //   (7) - at negative local y
-  tTransform =
-      new Transform3D(transform * AngleAxis3D(M_PI, Vector3D(1., 0., 0.)) *
-                      Translation3D(Vector3D(0., 2 * halflengthY1(), 0.)) *
-                      AngleAxis3D(-0.5 * M_PI, Vector3D(0., 1., 0.)) *
-                      AngleAxis3D(-0.5 * M_PI, Vector3D(1., 0., 0.)));
-  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-      std::shared_ptr<const Transform3D>(tTransform),
-      std::shared_ptr<const PlanarBounds>(faceZXRectangleBoundsBottom())));
-  //   (8) - at positive local y
-  tTransform = new Transform3D(
-      transform * Translation3D(Vector3D(0., 2 * halflengthY2(), 0.)) *
-      AngleAxis3D(-0.5 * M_PI, Vector3D(0., 1., 0.)) *
-      AngleAxis3D(-0.5 * M_PI, Vector3D(1., 0., 0.)));
-  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-      std::shared_ptr<const Transform3D>(tTransform),
-      std::shared_ptr<const PlanarBounds>(faceZXRectangleBoundsTop())));
-
-  return rSurfaces;
-}
-
-// faces in xy
-Acts::DiamondBounds* Acts::DoubleTrapezoidVolumeBounds::faceXYDiamondBounds()
-    const {
-  return new DiamondBounds(m_values.at(bv_minHalfX), m_values.at(bv_medHalfX),
-                           m_values.at(bv_maxHalfX), 2 * m_values.at(bv_halfY1),
-                           2 * m_values.at(bv_halfY2));
-}
-
-Acts::RectangleBounds*
-Acts::DoubleTrapezoidVolumeBounds::faceAlpha1RectangleBounds() const {
-  return new RectangleBounds(
-      m_values.at(bv_halfY1) / cos(m_values.at(bv_alpha1)),
-      m_values.at(bv_halfZ));
-}
-
-Acts::RectangleBounds*
-Acts::DoubleTrapezoidVolumeBounds::faceAlpha2RectangleBounds() const {
-  return new RectangleBounds(
-      m_values.at(bv_halfY2) / cos(m_values.at(bv_alpha2)),
-      m_values.at(bv_halfZ));
-}
-
-Acts::RectangleBounds*
-Acts::DoubleTrapezoidVolumeBounds::faceBeta1RectangleBounds() const {
-  return new RectangleBounds(
-      m_values.at(bv_halfY1) / cos(m_values.at(bv_alpha1)),
-      m_values.at(bv_halfZ));
-}
-
-Acts::RectangleBounds*
-Acts::DoubleTrapezoidVolumeBounds::faceBeta2RectangleBounds() const {
-  return new RectangleBounds(
-      m_values.at(bv_halfY2) / cos(m_values.at(bv_alpha2)),
-      m_values.at(bv_halfZ));
-}
-
-Acts::RectangleBounds*
-Acts::DoubleTrapezoidVolumeBounds::faceZXRectangleBoundsBottom() const {
-  return new RectangleBounds(m_values.at(bv_halfZ), m_values.at(bv_minHalfX));
-}
-
-Acts::RectangleBounds*
-Acts::DoubleTrapezoidVolumeBounds::faceZXRectangleBoundsTop() const {
-  return new RectangleBounds(m_values.at(bv_halfZ), m_values.at(bv_maxHalfX));
-}
-
-bool Acts::DoubleTrapezoidVolumeBounds::inside(const Vector3D& pos,
-                                               double tol) const {
-  if (std::abs(pos.z()) > m_values.at(bv_halfZ) + tol) {
-    return false;
-  }
-  if (pos.y() < -2 * m_values.at(bv_halfY1) - tol) {
-    return false;
-  }
-  if (pos.y() > 2 * m_values.at(bv_halfY2) - tol) {
-    return false;
-  }
-  DiamondBounds* faceXYBounds = faceXYDiamondBounds();
-  Vector2D locp(pos.x(), pos.y());
-  bool inside(faceXYBounds->inside(locp, BoundaryCheck(true, true, tol, tol)));
-  delete faceXYBounds;
-  return inside;
-}
-
-// ostream operator overload
-std::ostream& Acts::DoubleTrapezoidVolumeBounds::toStream(
-    std::ostream& sl) const {
-  return dumpT<std::ostream>(sl);
-}
-
-Acts::Volume::BoundingBox Acts::DoubleTrapezoidVolumeBounds::boundingBox(
-    const Transform3D* trf, const Vector3D& envelope,
-    const Volume* entity) const {
-  float minx = minHalflengthX();
-  float medx = medHalflengthX();
-  float maxx = maxHalflengthX();
-  float haley1 = 2 * halflengthY1();
-  float haley2 = 2 * halflengthY2();
-  float halez = halflengthZ();
-
-  std::array<Vector3D, 12> vertices = {{
-      {-minx, -haley1, -halez},
-      {+minx, -haley1, -halez},
-      {+medx, 0, -halez},
-      {-medx, 0, -halez},
-      {-maxx, +haley2, -halez},
-      {+maxx, +haley2, -halez},
-      {-minx, -haley1, +halez},
-      {+minx, -haley1, +halez},
-      {+medx, 0, +halez},
-      {-medx, 0, +halez},
-      {-maxx, +haley2, +halez},
-      {+maxx, +haley2, +halez},
-  }};
-
-  Transform3D transform = Transform3D::Identity();
-  if (trf != nullptr) {
-    transform = *trf;
-  }
-
-  Vector3D vmin = transform * vertices[0];
-  Vector3D vmax = transform * vertices[0];
-
-  for (size_t i = 1; i < 12; i++) {
-    const Vector3D vtx = transform * vertices[i];
-    vmin = vmin.cwiseMin(vtx);
-    vmax = vmax.cwiseMax(vtx);
-  }
-
-  return {entity, vmin - envelope, vmax + envelope};
-}
diff --git a/Core/src/Geometry/GenericCuboidVolumeBounds.cpp b/Core/src/Geometry/GenericCuboidVolumeBounds.cpp
index c6b9d496268cb320b0baabfee0b9500fc8a21bff..139d9e3302e5fa1212e991b34f7380c897d6b360 100644
--- a/Core/src/Geometry/GenericCuboidVolumeBounds.cpp
+++ b/Core/src/Geometry/GenericCuboidVolumeBounds.cpp
@@ -21,49 +21,20 @@
 #include "Acts/Utilities/Definitions.hpp"
 
 Acts::GenericCuboidVolumeBounds::GenericCuboidVolumeBounds(
-    const std::array<Acts::Vector3D, 8>& vertices)
+    const std::array<Acts::Vector3D, 8>& vertices) noexcept(false)
     : m_vertices(vertices) {
-  // calculate approximate center of gravity first, so we can make sure
-  // the normals point inwards
-  Vector3D cog(0, 0, 0);
+  construct();
+}
 
-  for (size_t i = 0; i < 8; i++) {
-    cog += m_vertices[i];
+Acts::GenericCuboidVolumeBounds::GenericCuboidVolumeBounds(
+    const std::array<double, GenericCuboidVolumeBounds::eSize>&
+        values) noexcept(false)
+    : m_vertices() {
+  for (size_t iv = 0; iv < 8; ++iv) {
+    m_vertices[iv] =
+        Vector3D(values[iv * 3], values[iv * 3 + 1], values[iv * 3 + 2]);
   }
-
-  cog *= 0.125;  // 1/8.
-
-  size_t idx = 0;
-
-  auto handle_face = [&](const auto& a, const auto& b, const auto& c,
-                         const auto& d) {
-    // we assume a b c d to be counter clockwise
-    const Vector3D ab = b - a, ac = c - a;
-    Vector3D normal = ab.cross(ac).normalized();
-
-    if ((cog - a).dot(normal) < 0) {
-      // normal points outwards, flip normal
-      normal *= -1.;
-    }
-
-    // get rid of -0 values if present
-    normal += Vector3D::Zero();
-
-    // check if d is on the surface
-    throw_assert((std::abs((a - d).dot(normal)) < 1e-6),
-                 "Four points do not lie on the same plane!");
-
-    m_normals[idx] = normal;
-    idx++;
-  };
-
-  // handle faces
-  handle_face(m_vertices[0], m_vertices[1], m_vertices[2], m_vertices[3]);
-  handle_face(m_vertices[4], m_vertices[5], m_vertices[6], m_vertices[7]);
-  handle_face(m_vertices[0], m_vertices[3], m_vertices[7], m_vertices[4]);
-  handle_face(m_vertices[1], m_vertices[2], m_vertices[6], m_vertices[5]);
-  handle_face(m_vertices[2], m_vertices[3], m_vertices[7], m_vertices[6]);
-  handle_face(m_vertices[1], m_vertices[0], m_vertices[4], m_vertices[5]);
+  construct();
 }
 
 Acts::VolumeBounds* Acts::GenericCuboidVolumeBounds::clone() const {
@@ -175,6 +146,63 @@ std::ostream& Acts::GenericCuboidVolumeBounds::toStream(
   return sl;
 }
 
+void Acts::GenericCuboidVolumeBounds::construct() noexcept(false) {
+  // calculate approximate center of gravity first, so we can make sure
+  // the normals point inwards
+  Vector3D cog(0, 0, 0);
+
+  for (size_t i = 0; i < 8; i++) {
+    cog += m_vertices[i];
+  }
+
+  cog *= 0.125;  // 1/8.
+
+  size_t idx = 0;
+
+  auto handle_face = [&](const auto& a, const auto& b, const auto& c,
+                         const auto& d) {
+    // we assume a b c d to be counter clockwise
+    const Vector3D ab = b - a, ac = c - a;
+    Vector3D normal = ab.cross(ac).normalized();
+
+    if ((cog - a).dot(normal) < 0) {
+      // normal points outwards, flip normal
+      normal *= -1.;
+    }
+
+    // get rid of -0 values if present
+    normal += Vector3D::Zero();
+
+    // check if d is on the surface
+    if (std::abs((a - d).dot(normal)) > 1e-6) {
+      throw(std::invalid_argument(
+          "GenericCuboidBounds: Four points do not lie on the same plane!"));
+    }
+
+    m_normals[idx] = normal;
+    idx++;
+  };
+
+  // handle faces
+  handle_face(m_vertices[0], m_vertices[1], m_vertices[2], m_vertices[3]);
+  handle_face(m_vertices[4], m_vertices[5], m_vertices[6], m_vertices[7]);
+  handle_face(m_vertices[0], m_vertices[3], m_vertices[7], m_vertices[4]);
+  handle_face(m_vertices[1], m_vertices[2], m_vertices[6], m_vertices[5]);
+  handle_face(m_vertices[2], m_vertices[3], m_vertices[7], m_vertices[6]);
+  handle_face(m_vertices[1], m_vertices[0], m_vertices[4], m_vertices[5]);
+}
+
+std::vector<double> Acts::GenericCuboidVolumeBounds::values() const {
+  std::vector<double> rvalues;
+  rvalues.reserve(eSize);
+  for (size_t iv = 0; iv < 8; ++iv) {
+    for (size_t ic = 0; ic < 3; ++ic) {
+      rvalues.push_back(m_vertices[iv][ic]);
+    }
+  }
+  return rvalues;
+}
+
 Acts::Volume::BoundingBox Acts::GenericCuboidVolumeBounds::boundingBox(
     const Acts::Transform3D* trf, const Vector3D& envelope,
     const Volume* entity) const {
diff --git a/Core/src/Geometry/TrapezoidVolumeBounds.cpp b/Core/src/Geometry/TrapezoidVolumeBounds.cpp
index 564a85450185153b2d49340f6cc6b6b9263c3ef1..f51049f1be21d5c73c261b7f30f862f648e0c7bc 100644
--- a/Core/src/Geometry/TrapezoidVolumeBounds.cpp
+++ b/Core/src/Geometry/TrapezoidVolumeBounds.cpp
@@ -12,59 +12,48 @@
 
 #include "Acts/Geometry/TrapezoidVolumeBounds.hpp"
 
-#include <cmath>
-#include <iomanip>
-#include <iostream>
 #include "Acts/Geometry/GeometryStatics.hpp"
 #include "Acts/Surfaces/PlaneSurface.hpp"
 #include "Acts/Surfaces/RectangleBounds.hpp"
 #include "Acts/Surfaces/TrapezoidBounds.hpp"
 #include "Acts/Utilities/BoundingBox.hpp"
 
-Acts::TrapezoidVolumeBounds::TrapezoidVolumeBounds()
-    : VolumeBounds(), m_values(bv_length, 0.) {}
+#include <cmath>
+#include <iomanip>
+#include <iostream>
 
 Acts::TrapezoidVolumeBounds::TrapezoidVolumeBounds(double minhalex,
                                                    double maxhalex,
                                                    double haley, double halez)
-    : VolumeBounds(), m_values(bv_length, 0.) {
-  m_values.at(bv_minHalfX) = minhalex;
-  m_values.at(bv_maxHalfX) = maxhalex;
-  m_values.at(bv_halfY) = haley;
-  m_values.at(bv_halfZ) = halez;
-  m_values.at(bv_alpha) =
-      atan((m_values.at(bv_maxHalfX) - m_values.at(bv_minHalfX)) / 2 /
-           m_values.at(bv_halfY)) +
+    : VolumeBounds() {
+  m_values[eHalfLengthXnegY] = minhalex;
+  m_values[eHalfLengthXposY] = maxhalex;
+  m_values[eHalfLengthY] = haley;
+  m_values[eHalfLengthZ] = halez;
+  m_values[eAlpha] =
+      atan((m_values[eHalfLengthXposY] - m_values[eHalfLengthXnegY]) / 2 /
+           m_values[eHalfLengthY]) +
       0.5 * M_PI;
-  m_values.at(bv_beta) = m_values.at(bv_alpha);
+  m_values[eBeta] = m_values[eAlpha];
+  checkConsistency();
+  buildSurfaceBounds();
 }
 
 Acts::TrapezoidVolumeBounds::TrapezoidVolumeBounds(double minhalex,
                                                    double haley, double halez,
                                                    double alpha, double beta)
-    : VolumeBounds(), m_values(bv_length, 0.) {
-  m_values.at(bv_minHalfX) = minhalex;
-  m_values.at(bv_halfY) = haley;
-  m_values.at(bv_halfZ) = halez;
-  m_values.at(bv_alpha) = alpha;
-  m_values.at(bv_beta) = beta;
+    : VolumeBounds() {
+  m_values[eHalfLengthXnegY] = minhalex;
+  m_values[eHalfLengthY] = haley;
+  m_values[eHalfLengthZ] = halez;
+  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);
-  m_values.at(bv_maxHalfX) = minhalex + (2. * haley) * tan(gamma);
-}
-
-Acts::TrapezoidVolumeBounds::TrapezoidVolumeBounds(
-    const TrapezoidVolumeBounds& trabo)
-    : VolumeBounds(), m_values(trabo.m_values) {}
-
-Acts::TrapezoidVolumeBounds::~TrapezoidVolumeBounds() = default;
+  m_values[eHalfLengthXposY] = minhalex + (2. * haley) * tan(gamma);
 
-Acts::TrapezoidVolumeBounds& Acts::TrapezoidVolumeBounds::operator=(
-    const TrapezoidVolumeBounds& trabo) {
-  if (this != &trabo) {
-    m_values = trabo.m_values;
-  }
-  return *this;
+  checkConsistency();
+  buildSurfaceBounds();
 }
 
 std::vector<std::shared_ptr<const Acts::Surface>>
@@ -83,130 +72,110 @@ Acts::TrapezoidVolumeBounds::decomposeToSurfaces(
   Vector3D trapezoidCenter(transform.translation());
 
   //   (1) - at negative local z
-  std::shared_ptr<const PlanarBounds> xytBounds(faceXYTrapezoidBounds());
   tTransform =
       new Transform3D(transform * AngleAxis3D(M_PI, Vector3D(0., 1., 0.)) *
-                      Translation3D(Vector3D(0., 0., halflengthZ())));
+                      Translation3D(Vector3D(0., 0., get(eHalfLengthZ))));
 
   rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-      std::shared_ptr<const Transform3D>(tTransform), xytBounds));
+      std::shared_ptr<const Transform3D>(tTransform), m_faceXYTrapezoidBounds));
   //   (2) - at positive local z
-  tTransform = new Transform3D(transform *
-                               Translation3D(Vector3D(0., 0., halflengthZ())));
+  tTransform = new Transform3D(
+      transform * Translation3D(Vector3D(0., 0., get(eHalfLengthZ))));
   rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-      std::shared_ptr<const Transform3D>(tTransform), xytBounds));
+      std::shared_ptr<const Transform3D>(tTransform), m_faceXYTrapezoidBounds));
 
   // face surfaces yz
   // transmute cyclical
   //   (3) - at point A, attached to alpha opening angle
-  Vector3D A(minHalflengthX(), halflengthY(), trapezoidCenter.z());
+  Vector3D A(get(eHalfLengthXnegY), get(eHalfLengthY), trapezoidCenter.z());
   RotationMatrix3D alphaZRotation =
-      (s_idRotation * AngleAxis3D(alpha() - 0.5 * M_PI, Vector3D(0., 0., 1.)))
+      (s_idRotation *
+       AngleAxis3D(get(eAlpha) - 0.5 * M_PI, Vector3D(0., 0., 1.)))
           .toRotationMatrix();
   RotationMatrix3D faceAlphaRotation;
   faceAlphaRotation.col(0) = alphaZRotation.col(1);
   faceAlphaRotation.col(1) = -alphaZRotation.col(2);
   faceAlphaRotation.col(2) = -alphaZRotation.col(0);
-  std::shared_ptr<const PlanarBounds> faceAlphaBounds(
-      faceAlphaRectangleBounds());
+
   // Vector3D
-  // faceAlphaPosition(A+faceAlphaRotation.colX()*faceAlphaBounds->halflengthX());
-  Vector3D faceAlphaPosition0(-0.5 * (minHalflengthX() + maxHalflengthX()), 0.,
-                              0.);
+  // faceAlphaPosition(A+faceAlphaRotation.colX()*m_faceAlphaRectangleBounds->halflengthX());
+  Vector3D faceAlphaPosition0(
+      -0.5 * (get(eHalfLengthXnegY) + get(eHalfLengthXposY)), 0., 0.);
   Vector3D faceAlphaPosition = transform * faceAlphaPosition0;
   tTransform = new Transform3D(Translation3D(faceAlphaPosition) *
                                (trapezoidRotation * faceAlphaRotation));
   rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-      std::shared_ptr<const Transform3D>(tTransform), faceAlphaBounds));
+      std::shared_ptr<const Transform3D>(tTransform),
+      m_faceAlphaRectangleBounds));
 
   //   (4) - at point B, attached to beta opening angle
-  Vector3D B(minHalflengthX(), -halflengthY(), trapezoidCenter.z());
+  Vector3D B(get(eHalfLengthXnegY), -get(eHalfLengthY), trapezoidCenter.z());
   RotationMatrix3D betaZRotation =
-      (s_idRotation * AngleAxis3D(-(beta() - 0.5 * M_PI), Vector3D(0., 0., 1.)))
+      (s_idRotation *
+       AngleAxis3D(-(get(eBeta) - 0.5 * M_PI), Vector3D(0., 0., 1.)))
           .toRotationMatrix();
   RotationMatrix3D faceBetaRotation;
   faceBetaRotation.col(0) = betaZRotation.col(1);
   faceBetaRotation.col(1) = betaZRotation.col(2);
   faceBetaRotation.col(2) = betaZRotation.col(0);
-  std::shared_ptr<const PlanarBounds> faceBetaBounds(faceBetaRectangleBounds());
   // Vector3D
-  // faceBetaPosition(B+faceBetaRotation.colX()*faceBetaBounds->halflengthX());
-  Vector3D faceBetaPosition0(0.5 * (minHalflengthX() + maxHalflengthX()), 0.,
-                             0.);
+  // faceBetaPosition(B+faceBetaRotation.colX()*m_faceBetaRectangleBounds->halflengthX());
+  Vector3D faceBetaPosition0(
+      0.5 * (get(eHalfLengthXnegY) + get(eHalfLengthXposY)), 0., 0.);
   Vector3D faceBetaPosition = transform * faceBetaPosition0;
   tTransform = new Transform3D(Translation3D(faceBetaPosition) *
                                (trapezoidRotation * faceBetaRotation));
   rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
-      std::shared_ptr<const Transform3D>(tTransform), faceBetaBounds));
+      std::shared_ptr<const Transform3D>(tTransform),
+      m_faceBetaRectangleBounds));
 
   // face surfaces zx
   //   (5) - at negative local x
   tTransform =
       new Transform3D(transform * AngleAxis3D(M_PI, Vector3D(1., 0., 0.)) *
-                      Translation3D(Vector3D(0., halflengthY(), 0.)) *
+                      Translation3D(Vector3D(0., get(eHalfLengthY), 0.)) *
                       AngleAxis3D(-0.5 * M_PI, Vector3D(0., 1., 0.)) *
                       AngleAxis3D(-0.5 * M_PI, Vector3D(1., 0., 0.)));
   rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
       std::shared_ptr<const Transform3D>(tTransform),
-      std::shared_ptr<const PlanarBounds>(faceZXRectangleBoundsBottom())));
+      std::shared_ptr<const PlanarBounds>(m_faceZXRectangleBounds)));
   //   (6) - at positive local x
-  tTransform = new Transform3D(transform *
-                               Translation3D(Vector3D(0., halflengthY(), 0.)) *
-                               AngleAxis3D(-0.5 * M_PI, Vector3D(0., 1., 0.)) *
-                               AngleAxis3D(-0.5 * M_PI, Vector3D(1., 0., 0.)));
+  tTransform = new Transform3D(
+      transform * Translation3D(Vector3D(0., get(eHalfLengthY), 0.)) *
+      AngleAxis3D(-0.5 * M_PI, Vector3D(0., 1., 0.)) *
+      AngleAxis3D(-0.5 * M_PI, Vector3D(1., 0., 0.)));
   rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
       std::shared_ptr<const Transform3D>(tTransform),
-      std::shared_ptr<const PlanarBounds>(faceZXRectangleBoundsTop())));
+      std::shared_ptr<const PlanarBounds>(m_faceZXRectangleBounds)));
 
   return rSurfaces;
 }
 
-Acts::TrapezoidBounds* Acts::TrapezoidVolumeBounds::faceXYTrapezoidBounds()
-    const {
-  return new TrapezoidBounds(m_values.at(bv_minHalfX), m_values.at(bv_maxHalfX),
-                             m_values.at(bv_halfY));
-}
+void Acts::TrapezoidVolumeBounds::buildSurfaceBounds() {
+  m_faceXYTrapezoidBounds = std::make_shared<const TrapezoidBounds>(
+      get(eHalfLengthXnegY), get(eHalfLengthXposY), get(eHalfLengthY));
 
-Acts::RectangleBounds* Acts::TrapezoidVolumeBounds::faceAlphaRectangleBounds()
-    const {
-  return new RectangleBounds(
-      m_values.at(bv_halfY) / cos(m_values.at(bv_alpha) - 0.5 * M_PI),
-      m_values.at(bv_halfZ));
-}
-
-Acts::RectangleBounds* Acts::TrapezoidVolumeBounds::faceBetaRectangleBounds()
-    const {
-  return new RectangleBounds(
-      m_values.at(bv_halfY) / cos(m_values.at(bv_beta) - 0.5 * M_PI),
-      m_values.at(bv_halfZ));
-}
+  m_faceAlphaRectangleBounds = std::make_shared<const RectangleBounds>(
+      get(eHalfLengthY) / cos(get(eAlpha) - 0.5 * M_PI), get(eHalfLengthZ));
 
-Acts::RectangleBounds*
-Acts::TrapezoidVolumeBounds::faceZXRectangleBoundsBottom() const {
-  return new RectangleBounds(m_values.at(bv_halfZ), m_values.at(bv_minHalfX));
-}
+  m_faceBetaRectangleBounds = std::make_shared<const RectangleBounds>(
+      get(eHalfLengthY) / cos(get(eBeta) - 0.5 * M_PI), get(eHalfLengthZ));
 
-Acts::RectangleBounds* Acts::TrapezoidVolumeBounds::faceZXRectangleBoundsTop()
-    const {
-  // double delta = (m_values.at(bv_alpha) < m_values.at(bv_beta)) ?
-  // m_values.at(bv_alpha) - M_PI/2. : m_values.at(bv_beta) - M_PI/2.;
-  // return new RectangleBounds(m_values.at(bv_halfZ),
-  // 0.5*(m_values.at(bv_minHalfX)+m_values.at(bv_minHalfX)+2.*m_values.at(bv_halfY)/cos(delta)));
-  return new RectangleBounds(m_values.at(bv_halfZ), m_values.at(bv_maxHalfX));
+  m_faceZXRectangleBounds = std::make_shared<const RectangleBounds>(
+      get(eHalfLengthZ), get(eHalfLengthXnegY));
 }
 
 bool Acts::TrapezoidVolumeBounds::inside(const Vector3D& pos,
                                          double tol) const {
-  if (std::abs(pos.z()) > m_values.at(bv_halfZ) + tol) {
+  if (std::abs(pos.z()) > get(eHalfLengthZ) + tol) {
     return false;
   }
-  if (std::abs(pos.y()) > m_values.at(bv_halfY) + tol) {
+  if (std::abs(pos.y()) > get(eHalfLengthY) + tol) {
     return false;
   }
-  TrapezoidBounds* faceXYBounds = faceXYTrapezoidBounds();
   Vector2D locp(pos.x(), pos.y());
-  bool inside(faceXYBounds->inside(locp, BoundaryCheck(true, true, tol, tol)));
-  delete faceXYBounds;
+  bool inside(m_faceXYTrapezoidBounds->inside(
+      locp, BoundaryCheck(true, true, tol, tol)));
   return inside;
 }
 
@@ -217,10 +186,10 @@ std::ostream& Acts::TrapezoidVolumeBounds::toStream(std::ostream& sl) const {
 Acts::Volume::BoundingBox Acts::TrapezoidVolumeBounds::boundingBox(
     const Acts::Transform3D* trf, const Vector3D& envelope,
     const Volume* entity) const {
-  double minx = minHalflengthX();
-  double maxx = maxHalflengthX();
-  double haley = halflengthY();
-  double halez = halflengthZ();
+  double minx = get(eHalfLengthXnegY);
+  double maxx = get(eHalfLengthXposY);
+  double haley = get(eHalfLengthY);
+  double halez = get(eHalfLengthZ);
 
   std::array<Vector3D, 8> vertices = {{{-minx, -haley, -halez},
                                        {+minx, -haley, -halez},
diff --git a/Tests/UnitTests/Core/Geometry/CMakeLists.txt b/Tests/UnitTests/Core/Geometry/CMakeLists.txt
index 7535c5ce500bcb67d3660857861ac11ba2b7f08f..0273c76373f71570ac74e961c6c30170bf155b67 100644
--- a/Tests/UnitTests/Core/Geometry/CMakeLists.txt
+++ b/Tests/UnitTests/Core/Geometry/CMakeLists.txt
@@ -1,11 +1,11 @@
 add_unittest(AlignmentContextTests AlignmentContextTests.cpp)
+add_unittest(CuboidVolumeBoundsTests CuboidVolumeBoundsTests.cpp)
 add_unittest(CuboidVolumeBuilderTests CuboidVolumeBuilderTests.cpp)
 add_unittest(CutoutCylinderVolumeBoundsTests CutoutCylinderVolumeBoundsTests.cpp)
 add_unittest(CylinderLayerTests CylinderLayerTests.cpp)
 add_unittest(CylinderVolumeBoundsTests CylinderVolumeBoundsTests.cpp)
 add_unittest(CylinderVolumeBuilderTests CylinderVolumeBuilderTests.cpp)
 add_unittest(DiscLayerTests DiscLayerTests.cpp)
-add_unittest(DoubleTrapezoidVolumeBoundsTests DoubleTrapezoidVolumeBoundsTests.cpp)
 add_unittest(ExtentTests ExtentTests.cpp)
 add_unittest(GenericApproachDescriptorTests GenericApproachDescriptorTests.cpp)
 add_unittest(GenericCuboidVolumeBoundsTests GenericCuboidVolumeBoundsTests.cpp)
diff --git a/Tests/UnitTests/Core/Geometry/CuboidVolumeBoundsTests.cpp b/Tests/UnitTests/Core/Geometry/CuboidVolumeBoundsTests.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..504f4373d0200fc1adb0f912727e05a332e47479
--- /dev/null
+++ b/Tests/UnitTests/Core/Geometry/CuboidVolumeBoundsTests.cpp
@@ -0,0 +1,96 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2019 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/.
+
+#include <boost/test/unit_test.hpp>
+
+#include "Acts/Geometry/CuboidVolumeBounds.hpp"
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
+#include "Acts/Utilities/Definitions.hpp"
+
+namespace Acts {
+namespace Test {
+
+GeometryContext gctx = GeometryContext();
+
+double hx{10.}, hy{20.}, hz{30.};
+
+BOOST_AUTO_TEST_SUITE(Geometry)
+
+BOOST_AUTO_TEST_CASE(CuboidVolumeConstruction) {
+  // Test Construction
+  CuboidVolumeBounds box(hx, hy, hz);
+
+  // Test copy construction
+  CuboidVolumeBounds copied(box);
+  BOOST_CHECK_EQUAL(box, copied);
+
+  // Test assigned
+  CuboidVolumeBounds assigned = box;
+  BOOST_CHECK_EQUAL(box, assigned);
+}
+
+BOOST_AUTO_TEST_CASE(CuboidVolumeRecreation) {
+  CuboidVolumeBounds original(hx, hy, hz);
+  auto valvector = original.values();
+  std::array<double, CuboidVolumeBounds::eSize> values;
+  std::copy_n(valvector.begin(), CuboidVolumeBounds::eSize, values.begin());
+  CuboidVolumeBounds recreated(values);
+  BOOST_CHECK_EQUAL(original, recreated);
+}
+
+BOOST_AUTO_TEST_CASE(CuboidVolumeException) {
+  // Test exception negative x
+  BOOST_CHECK_THROW(CuboidVolumeBounds(-hx, hy, hz), std::logic_error);
+  // Test exception negative y
+  BOOST_CHECK_THROW(CuboidVolumeBounds(hx, -hy, hz), std::logic_error);
+  // Test exception negative z
+  BOOST_CHECK_THROW(CuboidVolumeBounds(hx, hy, -hz), std::logic_error);
+  // Other iterations 0
+  BOOST_CHECK_THROW(CuboidVolumeBounds(-hx, hy, -hz), std::logic_error);
+  // Other iterations 1
+  BOOST_CHECK_THROW(CuboidVolumeBounds(-hx, -hy, hz), std::logic_error);
+  // Other iterations 2
+  BOOST_CHECK_THROW(CuboidVolumeBounds(hx, -hy, -hz), std::logic_error);
+  // Other iterations : all
+  BOOST_CHECK_THROW(CuboidVolumeBounds(-hx, -hy, -hz), std::logic_error);
+}
+
+BOOST_AUTO_TEST_CASE(CuboidVolumeProperties) {
+  CuboidVolumeBounds box(hx, hy, hz);
+  // Test the type
+  BOOST_TEST(box.type() == VolumeBounds::eCuboid);
+  // Test the halflength x
+  CHECK_CLOSE_ABS(box.get(CuboidVolumeBounds::eHalfLengthX), hx, s_epsilon);
+  // Test the halflength y
+  CHECK_CLOSE_ABS(box.get(CuboidVolumeBounds::eHalfLengthY), hy, s_epsilon);
+  // Test the halflength z
+  CHECK_CLOSE_ABS(box.get(CuboidVolumeBounds::eHalfLengthZ), hz, s_epsilon);
+  // Test the streaming
+  std::vector<double> refvalues = {hx, hy, hz};
+  BOOST_TEST(box.values() == refvalues);
+
+  // Inside position
+  Vector3D inside({5., 10., 8.});
+  // Outside positions  in x, y, z
+  std::vector<Vector3D> outsides = {
+      {20., 1., -2.}, {1., -30., 2.}, {-1., 2., 100.}};
+
+  // Inside position
+  BOOST_TEST(box.inside(inside, s_onSurfaceTolerance));
+
+  // Outside position
+  for (const auto& outside : outsides) {
+    BOOST_TEST(!box.inside(outside, s_onSurfaceTolerance));
+  }
+}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+}  // namespace Test
+}  // namespace Acts
diff --git a/Tests/UnitTests/Core/Geometry/CutoutCylinderVolumeBoundsTests.cpp b/Tests/UnitTests/Core/Geometry/CutoutCylinderVolumeBoundsTests.cpp
index 0fe43355cbf588a10e34f1ea2fde5034f66c7cbc..b9e3e3d737104154db81b960966bc00182c8ab7e 100644
--- a/Tests/UnitTests/Core/Geometry/CutoutCylinderVolumeBoundsTests.cpp
+++ b/Tests/UnitTests/Core/Geometry/CutoutCylinderVolumeBoundsTests.cpp
@@ -24,14 +24,80 @@
 namespace Acts {
 namespace Test {
 
-BOOST_AUTO_TEST_SUITE(Volumes)
+BOOST_AUTO_TEST_SUITE(Geometry)
 
-BOOST_AUTO_TEST_CASE(construction_test) {
+BOOST_AUTO_TEST_CASE(CutoutCylinderVolumeBoundsConstruction) {
   CutoutCylinderVolumeBounds ccvb(5, 10, 15, 30, 25);
   ccvb.toStream(std::cout);
+
+  // Test copy construction
+  CutoutCylinderVolumeBounds copied(ccvb);
+  BOOST_CHECK_EQUAL(ccvb, copied);
+
+  // Test assigned
+  CutoutCylinderVolumeBounds assigned = ccvb;
+  BOOST_CHECK_EQUAL(ccvb, assigned);
+}
+
+BOOST_AUTO_TEST_CASE(CutoutCylinderVolumeBoundsRecreation) {
+  CutoutCylinderVolumeBounds original(5, 10, 15, 30, 25);
+  std::array<double, CutoutCylinderVolumeBounds::eSize> values;
+  std::vector<double> valvector = original.values();
+  std::copy_n(valvector.begin(), CutoutCylinderVolumeBounds::eSize,
+              values.begin());
+  CutoutCylinderVolumeBounds recreated(values);
+  BOOST_CHECK_EQUAL(original, recreated);
+}
+
+BOOST_AUTO_TEST_CASE(CutoutCylinderVolumeBoundsExceptions) {
+  double rmin{5}, rmed{10}, rmax{15}, hz{30}, hzc{25};
+
+  // Test negative rmin
+  BOOST_CHECK_THROW(CutoutCylinderVolumeBounds(-rmin, rmed, rmax, hz, hzc),
+                    std::logic_error);
+
+  // Test negative rmed
+  BOOST_CHECK_THROW(CutoutCylinderVolumeBounds(rmin, -rmed, rmax, hz, hzc),
+                    std::logic_error);
+
+  // Test negative rmax
+  BOOST_CHECK_THROW(CutoutCylinderVolumeBounds(rmin, rmed, -rmax, hz, hzc),
+                    std::logic_error);
+
+  // Test swapped rmin / rmed
+  BOOST_CHECK_THROW(CutoutCylinderVolumeBounds(rmed, rmin, rmax, hz, hzc),
+                    std::logic_error);
+
+  // Test swapped rmin / rmax
+  BOOST_CHECK_THROW(CutoutCylinderVolumeBounds(rmax, rmed, rmin, hz, hzc),
+                    std::logic_error);
+
+  // Test swapped rmed / rmax
+  BOOST_CHECK_THROW(CutoutCylinderVolumeBounds(rmin, rmax, rmed, hz, hzc),
+                    std::logic_error);
+
+  // Test negative hz
+  BOOST_CHECK_THROW(CutoutCylinderVolumeBounds(rmin, rmed, rmax, -hz, hzc),
+                    std::logic_error);
+
+  // Test negative hzc
+  BOOST_CHECK_THROW(CutoutCylinderVolumeBounds(rmin, rmed, rmax, hz, -hzc),
+                    std::logic_error);
+}
+
+BOOST_AUTO_TEST_CASE(CutoutCylinderVolumeBoundsAccess) {
+  double rmin{5}, rmed{10}, rmax{15}, hz{30}, hzc{25};
+  CutoutCylinderVolumeBounds ccvb(rmin, rmed, rmax, hz, hzc);
+
+  BOOST_CHECK_EQUAL(ccvb.get(CutoutCylinderVolumeBounds::eMinR), rmin);
+  BOOST_CHECK_EQUAL(ccvb.get(CutoutCylinderVolumeBounds::eMedR), rmed);
+  BOOST_CHECK_EQUAL(ccvb.get(CutoutCylinderVolumeBounds::eMaxR), rmax);
+  BOOST_CHECK_EQUAL(ccvb.get(CutoutCylinderVolumeBounds::eHalfLengthZ), hz);
+  BOOST_CHECK_EQUAL(ccvb.get(CutoutCylinderVolumeBounds::eHalfLengthZcutout),
+                    hzc);
 }
 
-BOOST_AUTO_TEST_CASE(inside_test) {
+BOOST_AUTO_TEST_CASE(CutoutCylinderVolumeBoundsInside) {
   CutoutCylinderVolumeBounds ccvb(5, 10, 15, 30, 25);
 
   BOOST_CHECK(!ccvb.inside({0, 0, 0}));
@@ -101,7 +167,7 @@ BOOST_AUTO_TEST_CASE(inside_test) {
   BOOST_CHECK(!ccvb.inside({17, 0, -23}));
 }
 
-BOOST_AUTO_TEST_CASE(boundingbox_test) {
+BOOST_AUTO_TEST_CASE(CutoutCylinderVolumeBoundsBoundingBox) {
   GeometryContext tgContext = GeometryContext();
   std::vector<IdentifiedPolyderon> tPolyhedrons;
 
diff --git a/Tests/UnitTests/Core/Geometry/CylinderLayerTests.cpp b/Tests/UnitTests/Core/Geometry/CylinderLayerTests.cpp
index f106448721467396275db5f98f4d9294f0984849..d9a75acfbca1552fba849a8a95f4a67b8539bb81 100644
--- a/Tests/UnitTests/Core/Geometry/CylinderLayerTests.cpp
+++ b/Tests/UnitTests/Core/Geometry/CylinderLayerTests.cpp
@@ -45,7 +45,8 @@ BOOST_AUTO_TEST_CASE(CylinderLayerConstruction) {
   auto pTransform = std::make_shared<const Transform3D>(translation);
   double radius(0.5), halfz(10.);
   auto pCylinder = std::make_shared<const CylinderBounds>(radius, halfz);
-  auto pCylinderLayer = CylinderLayer::create(pTransform, pCylinder);
+  auto pCylinderLayer =
+      CylinderLayer::create(pTransform, pCylinder, nullptr, 1.);
   BOOST_CHECK_EQUAL(pCylinderLayer->layerType(), LayerType::passive);
   // next level: need an array of Surfaces;
   // bounds object, rectangle type
@@ -57,7 +58,7 @@ BOOST_AUTO_TEST_CASE(CylinderLayerConstruction) {
       Surface::makeShared<PlaneSurface>(pNullTransform, rBounds)};
   const double thickness(1.0);
   auto pCylinderLayerFromSurfaces =
-      CylinderLayer::create(pTransform, pCylinder, nullptr);
+      CylinderLayer::create(pTransform, pCylinder, nullptr, thickness);
   BOOST_CHECK_EQUAL(pCylinderLayerFromSurfaces->layerType(),
                     LayerType::passive);
   // construct with thickness:
@@ -86,7 +87,8 @@ BOOST_AUTO_TEST_CASE(CylinderLayerProperties /*, *utf::expected_failures(1)*/) {
   auto pTransform = std::make_shared<const Transform3D>(translation);
   double radius(0.5), halfz(10.);
   auto pCylinder = std::make_shared<const CylinderBounds>(radius, halfz);
-  auto pCylinderLayer = CylinderLayer::create(pTransform, pCylinder);
+  auto pCylinderLayer =
+      CylinderLayer::create(pTransform, pCylinder, nullptr, 1.);
   // auto planeSurface = pCylinderLayer->surfaceRepresentation();
   BOOST_CHECK_EQUAL(pCylinderLayer->surfaceRepresentation().name(),
                     std::string("Acts::CylinderSurface"));
diff --git a/Tests/UnitTests/Core/Geometry/CylinderVolumeBoundsTests.cpp b/Tests/UnitTests/Core/Geometry/CylinderVolumeBoundsTests.cpp
index 90c48854bcec73b6338ddd7211e0c581617f7a51..8925c77c20bff974d73e752e146b2cc605b22c33 100644
--- a/Tests/UnitTests/Core/Geometry/CylinderVolumeBoundsTests.cpp
+++ b/Tests/UnitTests/Core/Geometry/CylinderVolumeBoundsTests.cpp
@@ -11,6 +11,8 @@
 
 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
 #include "Acts/Geometry/Polyhedron.hpp"
+#include "Acts/Surfaces/CylinderBounds.hpp"
+#include "Acts/Surfaces/RadialBounds.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
 #include "Acts/Tests/CommonHelpers/ObjTestWriter.hpp"
@@ -24,10 +26,122 @@ namespace Acts {
 
 namespace Test {
 
-BOOST_AUTO_TEST_SUITE(Volumes)
+BOOST_AUTO_TEST_SUITE(Geometry)
+
+BOOST_AUTO_TEST_CASE(CylinderVolumeBoundsConstruction) {
+  double rmin{10.}, rmax{20.}, halfz{30.}, halfphi{M_PI / 4}, avgphi{0.};
+
+  // Test different construciton modes: solid
+  CylinderVolumeBounds solidCylinder(0., rmax, halfz);
+  BOOST_CHECK_EQUAL(solidCylinder.decomposeToSurfaces().size(), 3);
+
+  // Test different construciton modes: sectoral solid
+  CylinderVolumeBounds solidCylinderSector(0., rmax, halfz, halfphi);
+  BOOST_CHECK_EQUAL(solidCylinderSector.decomposeToSurfaces().size(), 5);
+
+  // Test different construciton modes: tube
+  CylinderVolumeBounds tubeCylinder(rmin, rmax, halfz);
+  BOOST_CHECK_EQUAL(tubeCylinder.decomposeToSurfaces().size(), 4);
+
+  // Test different construciton modes: sectoral tube
+  CylinderVolumeBounds tubeCylinderSector(rmin, rmax, halfz, halfphi);
+  BOOST_CHECK_EQUAL(tubeCylinderSector.decomposeToSurfaces().size(), 6);
+
+  CylinderVolumeBounds original(rmin, rmax, halfz, halfphi, avgphi);
+
+  // Test construction from CylinderBounds and thickness
+  double rmed = 0.5 * (rmin + rmax);
+  double rthickness = (rmax - rmin);
+  CylinderBounds cBounds(rmed, halfz, halfphi, avgphi);
+  CylinderVolumeBounds fromCylinder(cBounds, rthickness);
+  BOOST_CHECK_EQUAL(original, fromCylinder);
+
+  // Test construction from RadialBounds and thickness
+  RadialBounds rBounds(rmin, rmax, halfphi, avgphi);
+  CylinderVolumeBounds fromDisc(rBounds, 2 * halfz);
+  BOOST_CHECK_EQUAL(original, fromDisc);
+
+  // Test the copy construction
+  CylinderVolumeBounds copied(original);
+  BOOST_CHECK_EQUAL(original, copied);
+
+  // Test the assignment
+  CylinderVolumeBounds assigned = original;
+  BOOST_CHECK_EQUAL(original, assigned);
+}
+
+BOOST_AUTO_TEST_CASE(CylinderVolumeBoundsRecreation) {
+  double rmin{10.}, rmax{20.}, halfz{30.}, halfphi{M_PI / 4}, avgphi{0.};
+
+  CylinderVolumeBounds original(rmin, rmax, halfz, halfphi, avgphi);
+  std::array<double, CylinderVolumeBounds::eSize> values;
+  std::vector<double> valvector = original.values();
+  std::copy_n(valvector.begin(), CylinderVolumeBounds::eSize, values.begin());
+  CylinderVolumeBounds recreated(values);
+  BOOST_CHECK_EQUAL(original, recreated);
+}
+
+BOOST_AUTO_TEST_CASE(CylinderVolumeBoundsExceptions) {
+  double rmin{10.}, rmax{20.}, halfz{30.}, halfphi{M_PI / 4}, avgphi{0.};
+
+  // Negative inner radius
+  BOOST_CHECK_THROW(CylinderVolumeBounds(-rmin, rmax, halfz, halfphi, avgphi),
+                    std::logic_error);
+
+  // Negative outer radius
+  BOOST_CHECK_THROW(CylinderVolumeBounds(rmin, -rmax, halfz, halfphi, avgphi),
+                    std::logic_error);
+
+  // Swapped radii
+  BOOST_CHECK_THROW(CylinderVolumeBounds(rmax, rmin, halfz, halfphi, avgphi),
+                    std::logic_error);
+
+  // Zero half length
+  BOOST_CHECK_THROW(CylinderVolumeBounds(rmax, rmin, 0., halfphi, avgphi),
+                    std::logic_error);
+
+  // Negative half length
+  BOOST_CHECK_THROW(CylinderVolumeBounds(rmax, rmin, -halfz, halfphi, avgphi),
+                    std::logic_error);
+
+  // Out of bounds half phi
+  BOOST_CHECK_THROW(CylinderVolumeBounds(rmax, rmin, halfz, -4., avgphi),
+                    std::logic_error);
+
+  // Wrong positioning phi
+  BOOST_CHECK_THROW(CylinderVolumeBounds(rmax, rmin, halfz, halfphi, 4.),
+                    std::logic_error);
+
+  // Test construction from CylinderBounds and thickness
+  double rmed = 0.5 * (rmin + rmax);
+  CylinderBounds cBounds(rmed, halfz, halfphi, avgphi);
+  RadialBounds rBounds(rmin, rmax, halfphi, avgphi);
+
+  // Negative thickness
+  BOOST_CHECK_THROW(CylinderVolumeBounds(cBounds, -1.), std::logic_error);
+
+  // Wrong thickness
+  BOOST_CHECK_THROW(CylinderVolumeBounds(cBounds, 1000.), std::logic_error);
+
+  // Test construction from RadialBounds and thickness
+  BOOST_CHECK_THROW(CylinderVolumeBounds(rBounds, -1), std::logic_error);
+}
+
+BOOST_AUTO_TEST_CASE(CylinderVolumeBoundsAccess) {
+  double rmin{10.}, rmax{20.}, halfz{30.}, halfphi{M_PI / 4}, avgphi{0.};
+  CylinderVolumeBounds cvBounds(rmin, rmax, halfz, halfphi, avgphi);
+
+  // Test the accessors
+  BOOST_CHECK_EQUAL(cvBounds.get(CylinderVolumeBounds::eMinR), rmin);
+  BOOST_CHECK_EQUAL(cvBounds.get(CylinderVolumeBounds::eMaxR), rmax);
+  BOOST_CHECK_EQUAL(cvBounds.get(CylinderVolumeBounds::eHalfLengthZ), halfz);
+  BOOST_CHECK_EQUAL(cvBounds.get(CylinderVolumeBounds::eHalfPhiSector),
+                    halfphi);
+  BOOST_CHECK_EQUAL(cvBounds.get(CylinderVolumeBounds::eAveragePhi), avgphi);
+}
 
 /// Unit test for testing the decomposeToSurfaces() function
-BOOST_DATA_TEST_CASE(CylinderVolumeBounds_decomposeToSurfaces,
+BOOST_DATA_TEST_CASE(CylinderVolumeBoundsDecomposeToSurfaces,
                      bdata::random(-M_PI, M_PI) ^ bdata::random(-M_PI, M_PI) ^
                          bdata::random(-M_PI, M_PI) ^ bdata::random(-10., 10.) ^
                          bdata::random(-10., 10.) ^ bdata::random(-10., 10.) ^
@@ -48,7 +162,10 @@ BOOST_DATA_TEST_CASE(CylinderVolumeBounds_decomposeToSurfaces,
   AngleAxis3D rotZ(gamma, Vector3D(0., 0., 1.));
 
   // create the cylinder bounds
-  CylinderVolumeBounds cylBounds(1., 2., 3.);
+  double rmin = 1.;
+  double rmax = 2.;
+  double halfz = 3.;
+  CylinderVolumeBounds cylBounds(rmin, rmax, halfz);
   // create the transformation matrix
   auto mutableTransformPtr = std::make_shared<Transform3D>(Translation3D(pos));
   (*mutableTransformPtr) *= rotZ;
@@ -63,9 +180,9 @@ BOOST_DATA_TEST_CASE(CylinderVolumeBounds_decomposeToSurfaces,
 
   // check if difference is halfZ - sign and direction independent
   CHECK_CLOSE_REL((pos - boundarySurfaces.at(0)->center(tgContext)).norm(),
-                  cylBounds.halflengthZ(), 1e-12);
+                  cylBounds.get(CylinderVolumeBounds::eHalfLengthZ), 1e-12);
   CHECK_CLOSE_REL((pos - boundarySurfaces.at(1)->center(tgContext)).norm(),
-                  cylBounds.halflengthZ(), 1e-12);
+                  cylBounds.get(CylinderVolumeBounds::eHalfLengthZ), 1e-12);
   // transform to local
   double posDiscPosZ =
       (transformPtr->inverse() * boundarySurfaces.at(1)->center(tgContext)).z();
@@ -77,8 +194,12 @@ BOOST_DATA_TEST_CASE(CylinderVolumeBounds_decomposeToSurfaces,
   BOOST_CHECK_GT(centerPosZ, negDiscPosZ);
   // check positions of disc boundarysurfaces
   // checks for zero value. double precision value is not exact.
-  CHECK_CLOSE_ABS(negDiscPosZ + cylBounds.halflengthZ(), centerPosZ, 1e-12);
-  CHECK_CLOSE_ABS(posDiscPosZ - cylBounds.halflengthZ(), centerPosZ, 1e-12);
+  CHECK_CLOSE_ABS(
+      negDiscPosZ + cylBounds.get(CylinderVolumeBounds::eHalfLengthZ),
+      centerPosZ, 1e-12);
+  CHECK_CLOSE_ABS(
+      posDiscPosZ - cylBounds.get(CylinderVolumeBounds::eHalfLengthZ),
+      centerPosZ, 1e-12);
   // orientation of disc surfaces
   // positive disc durface should point in positive direction in the frame of
   // the volume
@@ -97,7 +218,7 @@ BOOST_DATA_TEST_CASE(CylinderVolumeBounds_decomposeToSurfaces,
   CHECK_CLOSE_REL(boundarySurfaces.at(2)->center(tgContext), pos, 1e-12);
 }
 
-BOOST_AUTO_TEST_CASE(bounding_box_creation) {
+BOOST_AUTO_TEST_CASE(CylinderVolumeBoundsBoundingBox) {
   GeometryContext tgContext = GeometryContext();
   std::vector<IdentifiedPolyderon> tPolyhedrons;
 
@@ -119,7 +240,7 @@ BOOST_AUTO_TEST_CASE(bounding_box_creation) {
 
   float tol = 1e-4;
 
-  CylinderVolumeBounds cvb(5, 10);
+  CylinderVolumeBounds cvb(0., 5, 10);
   auto cvbSurfaces = cvb.decomposeToSurfaces();
   combineAndDecompose(cvbSurfaces, "Solid");
   auto bb = cvb.boundingBox();
@@ -147,7 +268,7 @@ BOOST_AUTO_TEST_CASE(bounding_box_creation) {
   ObjTestWriter::writeObj("CylinderVolumeBounds_Tube_BB", bb);
 
   double angle = M_PI / 8.;
-  cvb = CylinderVolumeBounds(5, 8, angle, 13);
+  cvb = CylinderVolumeBounds(5, 8, 13, angle);
   bb = cvb.boundingBox();
   BOOST_CHECK_EQUAL(bb.entity(), nullptr);
   CHECK_CLOSE_ABS(bb.max(), Vector3D(8, 8 * std::sin(angle), 13), tol);
diff --git a/Tests/UnitTests/Core/Geometry/DoubleTrapezoidVolumeBoundsTests.cpp b/Tests/UnitTests/Core/Geometry/DoubleTrapezoidVolumeBoundsTests.cpp
deleted file mode 100644
index 78e69b8f7b7a72d7e615a8884cb65d52e4acaed4..0000000000000000000000000000000000000000
--- a/Tests/UnitTests/Core/Geometry/DoubleTrapezoidVolumeBoundsTests.cpp
+++ /dev/null
@@ -1,54 +0,0 @@
-// This file is part of the Acts project.
-//
-// Copyright (C) 2019 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/.
-
-#include <boost/test/unit_test.hpp>
-
-#include "Acts/Geometry/DoubleTrapezoidVolumeBounds.hpp"
-#include "Acts/Surfaces/PlanarBounds.hpp"
-#include "Acts/Surfaces/PlaneSurface.hpp"
-#include "Acts/Surfaces/Surface.hpp"
-#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
-#include "Acts/Utilities/BoundingBox.hpp"
-#include "Acts/Utilities/Definitions.hpp"
-
-namespace tt = boost::test_tools;
-
-namespace Acts {
-
-namespace Test {
-BOOST_AUTO_TEST_SUITE(Volumes)
-
-BOOST_AUTO_TEST_CASE(bounding_box_creation) {
-  float tol = 1e-4;
-
-  DoubleTrapezoidVolumeBounds dtvb(5, 10, 8, 4, 7, 3);
-
-  auto bb = dtvb.boundingBox();
-  CHECK_CLOSE_ABS(bb.max(), Vector3D(10, 14, 3), tol);
-  CHECK_CLOSE_ABS(bb.min(), Vector3D(-10, -8, -3), tol);
-
-  Transform3D trf;
-
-  trf = Translation3D(Vector3D(0, 30, 20));
-
-  bb = dtvb.boundingBox(&trf);
-  CHECK_CLOSE_ABS(bb.max(), Vector3D(10, 44, 23), tol);
-  CHECK_CLOSE_ABS(bb.min(), Vector3D(-10, 22, 17), tol);
-
-  trf = AngleAxis3D(M_PI / 2., Vector3D(-2, 4, 5).normalized());
-
-  bb = dtvb.boundingBox(&trf);
-  CHECK_CLOSE_ABS(bb.max(), Vector3D(8.9517, 11.7462, 10.263), tol);
-  CHECK_CLOSE_ABS(bb.min(), Vector3D(-14.7572, -7.9101, -9.85174), tol);
-}
-
-BOOST_AUTO_TEST_SUITE_END()
-
-}  // namespace Test
-
-}  // namespace Acts