From 69cee36bc73f3f467c5c92991d9de02e21543f22 Mon Sep 17 00:00:00 2001
From: Paul Gessinger <paul.gessinger@cern.ch>
Date: Thu, 8 Mar 2018 11:05:57 +0100
Subject: [PATCH] make CylinderLayer and DiscLayer serializable

still has some caveats, custom approach descriptors are currently not stored
---
 Core/include/ACTS/Layers/CylinderLayer.hpp    |  62 +--
 Core/include/ACTS/Layers/DiscLayer.hpp        |   7 +
 Core/include/ACTS/Surfaces/CylinderBounds.hpp |   2 +-
 Core/include/ACTS/Surfaces/DiamondBounds.hpp  |   4 +-
 Core/include/ACTS/Surfaces/DiscSurface.hpp    |   4 +
 Core/include/ACTS/Surfaces/EllipseBounds.hpp  |   4 +-
 Core/include/ACTS/Surfaces/LineBounds.hpp     |   6 +
 Core/include/ACTS/Surfaces/LineSurface.hpp    |   6 +
 Core/include/ACTS/Surfaces/PerigeeSurface.hpp |   6 +
 Core/include/ACTS/Surfaces/PlanarBounds.hpp   |   3 +-
 Core/include/ACTS/Surfaces/PlaneSurface.hpp   |   6 +-
 .../include/ACTS/Surfaces/RectangleBounds.hpp |   3 +-
 Core/include/ACTS/Surfaces/StrawSurface.hpp   |   4 +
 Core/include/ACTS/Surfaces/Surface.hpp        |   4 +
 Core/include/ACTS/Surfaces/SurfaceArray.hpp   | 389 +++++++++++++++++-
 Core/include/ACTS/Surfaces/SurfaceBounds.hpp  |   9 +-
 .../include/ACTS/Surfaces/TrapezoidBounds.hpp |   4 +-
 Core/include/ACTS/Surfaces/TriangleBounds.hpp |   8 +-
 .../ACTS/Utilities/InstanceFactory.hpp        |  34 +-
 Core/include/ACTS/Utilities/VariantData.hpp   |  28 +-
 .../include/ACTS/Utilities/VariantDataFwd.hpp |  16 +-
 Core/src/Layers/CylinderLayer.cpp             | 112 +++++
 Core/src/Layers/DiscLayer.cpp                 | 109 ++++-
 Core/src/Surfaces/CylinderBounds.cpp          |  26 +-
 Core/src/Surfaces/DiamondBounds.cpp           |   4 +-
 Core/src/Surfaces/EllipseBounds.cpp           |   4 +-
 Core/src/Surfaces/LineBounds.cpp              |  32 ++
 Core/src/Surfaces/LineSurface.cpp             |  47 +++
 Core/src/Surfaces/PerigeeSurface.cpp          |  43 +-
 Core/src/Surfaces/PlaneSurface.cpp            |   4 +-
 Core/src/Surfaces/RectangleBounds.cpp         |  16 +-
 Core/src/Surfaces/TrapezoidBounds.cpp         |  15 +-
 Core/src/Surfaces/TriangleBounds.cpp          |  30 +-
 Core/src/Tools/SurfaceArrayCreator.cpp        |  35 +-
 Tests/Layers/CylinderLayerTests.cpp           |  10 +-
 Tests/Layers/DiscLayerTests.cpp               |  37 ++
 Tests/Surfaces/CylinderBoundsTests.cpp        |   8 +-
 Tests/Surfaces/DiamondBoundsTests.cpp         |   7 +-
 Tests/Surfaces/EllipseBoundsTests.cpp         |   7 +-
 Tests/Surfaces/LineSurfaceStub.hpp            |   7 +
 Tests/Surfaces/LineSurfaceTests.cpp           |  25 ++
 Tests/Surfaces/PerigeeSurfaceTests.cpp        |  19 +
 Tests/Surfaces/PlaneSurfaceTests.cpp          |  25 +-
 Tests/Surfaces/RectangleBoundsTests.cpp       |   5 +-
 Tests/Surfaces/SurfaceArrayTests.cpp          | 158 ++++++-
 Tests/Surfaces/SurfaceStub.hpp                |  10 +
 Tests/Surfaces/TrapezoidBoundsTests.cpp       |  12 +-
 Tests/Surfaces/TriangleBoundsTests.cpp        |  45 +-
 Tests/Tools/LayerCreatorTests.cpp             | 104 ++++-
 49 files changed, 1321 insertions(+), 244 deletions(-)

diff --git a/Core/include/ACTS/Layers/CylinderLayer.hpp b/Core/include/ACTS/Layers/CylinderLayer.hpp
index b29cb202f..08ab43ed0 100644
--- a/Core/include/ACTS/Layers/CylinderLayer.hpp
+++ b/Core/include/ACTS/Layers/CylinderLayer.hpp
@@ -68,35 +68,7 @@ public:
   }
 
   static MutableLayerPtr
-  create(const variant_data& data_)
-  {
-    throw_assert(data_.which() == 4, "Variant data must be map");
-    const variant_map& data = boost::get<variant_map>(data_);
-    std::string        type = data.get<std::string>("type");
-    throw_assert(type == "CylinderLayer", "Type must be CylinderLayer");
-
-    variant_map payload = data.get<variant_map>("payload");
-
-    auto trf = std::make_shared<const Transform3D>(
-        from_variant<Transform3D>(payload.get<variant_map>("transform")));
-
-    LayerType laytyp = passive;
-
-    double thickness = payload.get<double>("thickness");
-
-    InstanceFactory    factory;
-    const variant_map& var_bounds = payload.get<variant_map>("cylinder_bounds");
-
-    auto cbounds = std::dynamic_pointer_cast<const CylinderBounds>(
-        factory.surfaceBounds(var_bounds.get<std::string>("type"), var_bounds));
-
-    return MutableLayerPtr(new CylinderLayer(trf,
-                                             cbounds,
-                                             nullptr,  // SurfaceArray
-                                             thickness,
-                                             nullptr,  // std::move(ad),
-                                             laytyp));
-  }
+  create(const variant_data& data_);
 
   /// Copy constructor - deleted
   CylinderLayer(const CylinderLayer& cla) = delete;
@@ -122,37 +94,7 @@ public:
   surfaceRepresentation() override;
 
   variant_data
-  toVariantData() const
-  {
-    using namespace std::string_literals;
-    variant_map payload;
-
-    payload["transform"] = to_variant(*m_transform);
-
-    // we need to recover the bounds
-    const AbstractVolume* absVol = representingVolume();
-    auto                  cvBounds
-        = dynamic_cast<const CylinderVolumeBounds*>(&absVol->volumeBounds());
-
-    double cylR = cvBounds->innerRadius() + 0.5 * thickness();
-    double hlZ  = cvBounds->halflengthZ();
-
-    CylinderBounds cylBounds(cylR, hlZ);
-
-    const variant_data bounds  = cylBounds.toVariantData();
-    payload["cylinder_bounds"] = bounds;
-
-    payload["thickness"] = thickness();
-
-    // payload["surface_array"] = m_surfaceArray->toVariantData();
-    payload["surface_array"] = 0;
-
-    variant_map data;
-    data["type"]    = "CylinderLayer"s;
-    data["payload"] = payload;
-
-    return data;
-  }
+  toVariantData() const override;
 
 private:
   /// build approach surfaces */
diff --git a/Core/include/ACTS/Layers/DiscLayer.hpp b/Core/include/ACTS/Layers/DiscLayer.hpp
index 99cbad25e..8094c6d79 100644
--- a/Core/include/ACTS/Layers/DiscLayer.hpp
+++ b/Core/include/ACTS/Layers/DiscLayer.hpp
@@ -19,6 +19,7 @@ class MsgStream;
 #include "ACTS/Layers/Layer.hpp"
 #include "ACTS/Surfaces/DiscSurface.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 
 namespace Acts {
 
@@ -63,6 +64,9 @@ public:
                                          laytyp));
   }
 
+  static MutableLayerPtr
+  create(const variant_data& data);
+
   /// Default Constructor
   DiscLayer() = delete;
 
@@ -86,6 +90,9 @@ public:
   DiscSurface&
   surfaceRepresentation() override;
 
+  variant_data
+  toVariantData() const override;
+
 private:
   /// build approach surfaces
   void
diff --git a/Core/include/ACTS/Surfaces/CylinderBounds.hpp b/Core/include/ACTS/Surfaces/CylinderBounds.hpp
index 0073d0d23..e7911d57f 100644
--- a/Core/include/ACTS/Surfaces/CylinderBounds.hpp
+++ b/Core/include/ACTS/Surfaces/CylinderBounds.hpp
@@ -136,7 +136,7 @@ public:
   /// This method returns the halflengthZ
   double
   halflengthZ() const;
-  
+
   variant_data
   toVariantData() const;
 
diff --git a/Core/include/ACTS/Surfaces/DiamondBounds.hpp b/Core/include/ACTS/Surfaces/DiamondBounds.hpp
index 2ccc6f4f9..bf3bd7b9b 100644
--- a/Core/include/ACTS/Surfaces/DiamondBounds.hpp
+++ b/Core/include/ACTS/Surfaces/DiamondBounds.hpp
@@ -54,7 +54,7 @@ public:
                 double haley1,
                 double haley2);
 
-  DiamondBounds(const variant_data &data);
+  DiamondBounds(const variant_data& data);
 
   virtual ~DiamondBounds();
 
@@ -121,7 +121,7 @@ public:
   /// This method returns the halflength in Y of trapezoid at positive Y
   double
   halflengthY2() const;
-  
+
   variant_data
   toVariantData() const override;
 
diff --git a/Core/include/ACTS/Surfaces/DiscSurface.hpp b/Core/include/ACTS/Surfaces/DiscSurface.hpp
index f2b3471e9..a0f6c907b 100644
--- a/Core/include/ACTS/Surfaces/DiscSurface.hpp
+++ b/Core/include/ACTS/Surfaces/DiscSurface.hpp
@@ -18,6 +18,7 @@
 #include "ACTS/Surfaces/Surface.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/Identifier.hpp"
+#include "ACTS/Utilities/VariantDataFwd.hpp"
 
 namespace Acts {
 
@@ -311,6 +312,9 @@ public:
   virtual std::string
   name() const override;
 
+  virtual variant_data
+  toVariantData() const override;
+
 protected:
   std::shared_ptr<const DiscBounds> m_bounds;  ///< bounds (shared)
 };
diff --git a/Core/include/ACTS/Surfaces/EllipseBounds.hpp b/Core/include/ACTS/Surfaces/EllipseBounds.hpp
index cfabfece6..e55b70850 100644
--- a/Core/include/ACTS/Surfaces/EllipseBounds.hpp
+++ b/Core/include/ACTS/Surfaces/EllipseBounds.hpp
@@ -63,7 +63,7 @@ public:
                 double averagePhi = 0.,
                 double halfPhi    = M_PI);
 
-  EllipseBounds(const variant_data &data);
+  EllipseBounds(const variant_data& data);
 
   virtual ~EllipseBounds();
 
@@ -129,7 +129,7 @@ public:
   /// This method returns the halfPhiSector which is covered by the disc
   double
   halfPhiSector() const;
-  
+
   variant_data
   toVariantData() const override;
 
diff --git a/Core/include/ACTS/Surfaces/LineBounds.hpp b/Core/include/ACTS/Surfaces/LineBounds.hpp
index b6fbe6b50..f3d169a23 100644
--- a/Core/include/ACTS/Surfaces/LineBounds.hpp
+++ b/Core/include/ACTS/Surfaces/LineBounds.hpp
@@ -15,6 +15,7 @@
 
 #include "ACTS/Surfaces/SurfaceBounds.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
+#include "ACTS/Utilities/VariantDataFwd.hpp"
 
 namespace Acts {
 
@@ -36,6 +37,8 @@ public:
   /// @param halez is the half length in z, defualt = 0.
   LineBounds(double radius = 0., double halez = 0.);
 
+  LineBounds(const variant_data& data);
+
   virtual ~LineBounds();
 
   virtual LineBounds*
@@ -81,6 +84,9 @@ public:
   virtual std::ostream&
   dump(std::ostream& sl) const final override;
 
+  virtual variant_data
+  toVariantData() const;
+
 private:
   double m_radius, m_halfZ;
 };
diff --git a/Core/include/ACTS/Surfaces/LineSurface.hpp b/Core/include/ACTS/Surfaces/LineSurface.hpp
index 3dcebd306..63eba981b 100644
--- a/Core/include/ACTS/Surfaces/LineSurface.hpp
+++ b/Core/include/ACTS/Surfaces/LineSurface.hpp
@@ -17,6 +17,7 @@
 #include "ACTS/Surfaces/LineBounds.hpp"
 #include "ACTS/Surfaces/Surface.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
+#include "ACTS/Utilities/VariantDataFwd.hpp"
 
 namespace Acts {
 
@@ -74,6 +75,8 @@ public:
   /// @param transf The additional transform applied after copying
   LineSurface(const LineSurface& other, const Transform3D& transf);
 
+  LineSurface(const variant_data& data);
+
   virtual ~LineSurface();
 
   /// Assignment operator
@@ -272,6 +275,9 @@ public:
   virtual std::string
   name() const override;
 
+  virtual variant_data
+  toVariantData() const override;
+
 protected:
   std::shared_ptr<const LineBounds> m_bounds;  ///< bounds (shared)
 
diff --git a/Core/include/ACTS/Surfaces/PerigeeSurface.hpp b/Core/include/ACTS/Surfaces/PerigeeSurface.hpp
index e156718e3..dff4e9046 100644
--- a/Core/include/ACTS/Surfaces/PerigeeSurface.hpp
+++ b/Core/include/ACTS/Surfaces/PerigeeSurface.hpp
@@ -17,6 +17,7 @@
 #include "ACTS/Surfaces/LineSurface.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/GeometryStatics.hpp"
+#include "ACTS/Utilities/VariantDataFwd.hpp"
 
 namespace Acts {
 
@@ -53,6 +54,8 @@ public:
   /// @param transf is the transformed applied after copying
   PerigeeSurface(const PerigeeSurface& other, const Transform3D& transf);
 
+  PerigeeSurface(const variant_data& data);
+
   virtual ~PerigeeSurface();
 
   /// Virtual constructor
@@ -80,6 +83,9 @@ public:
   /// @param sl is the ostream to be dumped into
   virtual std::ostream&
   dump(std::ostream& sl) const final override;
+
+  virtual variant_data
+  toVariantData() const override;
 };
 
 }  // end of namespace
diff --git a/Core/include/ACTS/Surfaces/PlanarBounds.hpp b/Core/include/ACTS/Surfaces/PlanarBounds.hpp
index 9e61ccd83..988f810db 100644
--- a/Core/include/ACTS/Surfaces/PlanarBounds.hpp
+++ b/Core/include/ACTS/Surfaces/PlanarBounds.hpp
@@ -40,8 +40,7 @@ public:
   virtual const RectangleBounds&
   boundingBox() const = 0;
 
-  virtual 
-  variant_data
+  virtual variant_data
   toVariantData() const = 0;
 };
 
diff --git a/Core/include/ACTS/Surfaces/PlaneSurface.hpp b/Core/include/ACTS/Surfaces/PlaneSurface.hpp
index 826fc951d..d3a684356 100644
--- a/Core/include/ACTS/Surfaces/PlaneSurface.hpp
+++ b/Core/include/ACTS/Surfaces/PlaneSurface.hpp
@@ -74,7 +74,7 @@ public:
   PlaneSurface(std::shared_ptr<const Transform3D>  htrans,
                std::shared_ptr<const PlanarBounds> pbounds = nullptr);
 
-  PlaneSurface(const variant_data &data);
+  PlaneSurface(const variant_data& data);
 
   virtual ~PlaneSurface();
 
@@ -207,8 +207,8 @@ public:
   virtual std::string
   name() const override;
 
-  variant_data
-  toVariantData() const;
+  virtual variant_data
+  toVariantData() const override;
 
 protected:
   /// the bounds of this surface
diff --git a/Core/include/ACTS/Surfaces/RectangleBounds.hpp b/Core/include/ACTS/Surfaces/RectangleBounds.hpp
index 43af5380a..e0b5d0bb5 100644
--- a/Core/include/ACTS/Surfaces/RectangleBounds.hpp
+++ b/Core/include/ACTS/Surfaces/RectangleBounds.hpp
@@ -19,7 +19,6 @@
 
 namespace Acts {
 
-
 /// @class RectangleBounds
 ///
 /// Bounds for a rectangular, planar surface.
@@ -96,7 +95,7 @@ public:
   /// Return method for the half length in Y
   double
   halflengthY() const;
-  
+
   variant_data
   toVariantData() const override;
 
diff --git a/Core/include/ACTS/Surfaces/StrawSurface.hpp b/Core/include/ACTS/Surfaces/StrawSurface.hpp
index ce191a511..5eca2263e 100644
--- a/Core/include/ACTS/Surfaces/StrawSurface.hpp
+++ b/Core/include/ACTS/Surfaces/StrawSurface.hpp
@@ -17,6 +17,7 @@
 #include "ACTS/Surfaces/LineSurface.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/Identifier.hpp"
+#include "ACTS/Utilities/VariantDataFwd.hpp"
 
 namespace Acts {
 
@@ -95,6 +96,9 @@ public:
   /// Return properly formatted class name for screen output */
   virtual std::string
   name() const final override;
+
+  virtual variant_data
+  toVariantData() const override;
 };
 
 inline Surface::SurfaceType
diff --git a/Core/include/ACTS/Surfaces/Surface.hpp b/Core/include/ACTS/Surfaces/Surface.hpp
index 1dd3e8323..96632a748 100644
--- a/Core/include/ACTS/Surfaces/Surface.hpp
+++ b/Core/include/ACTS/Surfaces/Surface.hpp
@@ -23,6 +23,7 @@
 #include "ACTS/Utilities/GeometryStatics.hpp"
 #include "ACTS/Utilities/Identifier.hpp"
 #include "ACTS/Utilities/Intersection.hpp"
+#include "ACTS/Utilities/VariantDataFwd.hpp"
 
 namespace Acts {
 
@@ -391,6 +392,9 @@ public:
   virtual std::string
   name() const = 0;
 
+  virtual variant_data
+  toVariantData() const = 0;
+
 protected:
   /// Transform3D definition that positions (translation, rotation) the surface
   /// in global space
diff --git a/Core/include/ACTS/Surfaces/SurfaceArray.hpp b/Core/include/ACTS/Surfaces/SurfaceArray.hpp
index 7b1b78a54..089d09c1b 100644
--- a/Core/include/ACTS/Surfaces/SurfaceArray.hpp
+++ b/Core/include/ACTS/Surfaces/SurfaceArray.hpp
@@ -16,6 +16,8 @@
 #include "ACTS/Utilities/BinningType.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/IAxis.hpp"
+#include "ACTS/Utilities/InstanceFactory.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 #include "ACTS/Utilities/detail/Axis.hpp"
 #include "ACTS/Utilities/detail/Grid.hpp"
 
@@ -122,6 +124,10 @@ public:
     ///       or overflow bin or out of range in any axis.
     virtual bool
     isValidBin(size_t bin) const = 0;
+
+    virtual
+    variant_data
+    toVariantData() const = 0;
   };
 
   /// @brief Lookup helper which encapsulates a @c Grid
@@ -320,6 +326,54 @@ public:
       return true;
     }
 
+    virtual
+    variant_data
+    toVariantData() const override
+    {
+      using namespace std::string_literals;
+      variant_map payload;
+
+      payload["dimensions"] = int(DIM);
+      variant_vector axes;
+
+      for (const auto& axis : m_grid.getAxes()) {
+        variant_map ax_pl;
+
+        ax_pl["axistype"] = axis->isEquidistant() ? "equidistant"s : "variable"s;
+
+        if (axis->isEquidistant()) {
+          ax_pl["min"]   = axis->getMin();
+          ax_pl["max"]   = axis->getMax();
+          ax_pl["nbins"] = int(axis->getNBins());
+        } else {
+          variant_vector bin_edges;
+          for (const auto& bin_edge : axis->getBinEdges()) {
+            bin_edges.push_back(bin_edge);
+          }
+          ax_pl["bin_edges"] = bin_edges;
+        }
+
+        if (axis->getBoundaryType() == detail::AxisBoundaryType::Open) {
+          ax_pl["axisboundarytype"] = "open"s;
+        } else if (axis->getBoundaryType() == detail::AxisBoundaryType::Bound) {
+          ax_pl["axisboundarytype"] = "bound"s;
+        } else if (axis->getBoundaryType() == detail::AxisBoundaryType::Closed) {
+          ax_pl["axisboundarytype"] = "closed"s;
+        }
+
+        variant_map ax_data;
+        ax_data["type"]    = "Axis"s;
+        ax_data["payload"] = ax_pl;
+        axes.push_back(ax_data);
+      }
+      payload["axes"] = axes;
+
+      variant_map data;
+      data["type"]    = "SurfaceGridLookup"s;
+      data["payload"] = payload;
+      return data;
+    }
+
   private:
     void
     populateNeighborCache()
@@ -479,6 +533,17 @@ public:
     /// @return always true
     virtual bool isValidBin(size_t) const override { return true; }
 
+    virtual
+    variant_data
+    toVariantData() const override
+    {
+      using namespace std::string_literals;
+      variant_map data;
+      data["type"]    = "SingleElementLookup"s;
+      data["payload"] = variant_map();
+      return data;
+    }
+
   private:
     SurfaceVector m_element;
   };
@@ -490,8 +555,11 @@ public:
   /// @param surfaces The input vector of surfaces. This is only for
   /// bookkeeping, so we can ask
   SurfaceArray(std::unique_ptr<ISurfaceGridLookup> gridLookup,
-               SurfaceVector                       surfaces)
-    : p_gridLookup(std::move(gridLookup)), m_surfaces(surfaces)
+               SurfaceVector                       surfaces,
+               std::shared_ptr<const Transform3D>  transform = nullptr)
+    : p_gridLookup(std::move(gridLookup))
+    , m_surfaces(surfaces)
+    , m_transform(transform)
   {
   }
 
@@ -500,9 +568,12 @@ public:
   /// @param surfaces The input vector of surfaces. This is only for
   /// bookkeeping, so we can ask
   template <class SGL>
-  SurfaceArray(std::unique_ptr<SGL> gridLookup, SurfaceVector surfaces)
+  SurfaceArray(std::unique_ptr<SGL>               gridLookup,
+               SurfaceVector                      surfaces,
+               std::shared_ptr<const Transform3D> transform = nullptr)
     : p_gridLookup(static_cast<ISurfaceGridLookup*>(gridLookup.release()))
     , m_surfaces(surfaces)
+    , m_transform(transform)
   {
   }
 
@@ -515,6 +586,220 @@ public:
   {
   }
 
+  SurfaceArray(const variant_data&                      data_,
+               std::function<Vector2D(const Vector3D&)> g2l,
+               std::function<Vector3D(const Vector2D&)> l2g,
+               std::shared_ptr<const Transform3D>       transform = nullptr)
+    : m_gridLookup(nullptr), m_transform(transform)
+  {
+    const variant_map& data = boost::get<variant_map>(data_);
+    throw_assert(data.get<std::string>("type") == "SurfaceArray",
+                 "Type must be SurfaceArray");
+
+    const variant_map& payload = data.get<variant_map>("payload");
+    const variant_map& var_sgl = payload.get<variant_map>("surfacegridlookup");
+    throw_assert(var_sgl.get<std::string>("type") == "SurfaceGridLookup",
+                 "Currently only SurfaceGridLookup can be deserialized");
+    const variant_vector& var_surfaces
+        = payload.get<variant_vector>("surfaces");
+
+    InstanceFactory factory;
+
+    std::vector<const Surface*> surfaces;
+    for (const auto& var_srf_ : var_surfaces) {
+      const variant_map& var_srf = boost::get<variant_map>(var_srf_);
+      const Surface*     surface
+          = factory.surface(var_srf.get<std::string>("type"), var_srf);
+      surfaces.push_back(surface);
+    }
+
+    m_surfaces = surfaces;
+
+    // reproduce axes
+    const variant_vector& var_axes
+        = var_sgl.get<variant_map>("payload").get<variant_vector>("axes");
+
+    throw_assert(var_axes.size() == 2,
+                 "This constructor cannot deserialize DIM!=2 data");
+
+    // two dimensional
+    const variant_map& var_axis_a
+        = var_axes.get<variant_map>(0).get<variant_map>("payload");
+    const variant_map& var_axis_b
+        = var_axes.get<variant_map>(1).get<variant_map>("payload");
+
+    std::string axistype_a = var_axis_a.get<std::string>("axistype");
+    std::string axisbdt_a  = var_axis_a.get<std::string>("axisboundarytype");
+    std::string axistype_b = var_axis_b.get<std::string>("axistype");
+    std::string axisbdt_b  = var_axis_b.get<std::string>("axisboundarytype");
+
+    auto makePAxis = [](const std::string& axistype,
+                        const variant_map& var_axis) -> ProtoAxis {
+      ProtoAxis pAxis;
+      if (axistype == "equidistant") {
+        pAxis.bType = equidistant;
+        pAxis.min   = var_axis.get<double>("min");
+        pAxis.max   = var_axis.get<double>("max");
+        pAxis.nBins = var_axis.get<int>("nbins");
+      } else {
+        std::vector<double>   bin_edges;
+        const variant_vector& var_bin_edges
+            = var_axis.get<variant_vector>("bin_edges");
+        for (size_t i = 0; i < var_bin_edges.size(); i++) {
+          bin_edges.push_back(var_bin_edges.get<double>(i));
+        }
+
+        pAxis.bType    = arbitrary;
+        pAxis.binEdges = bin_edges;
+      }
+
+      return pAxis;
+    };
+
+    ProtoAxis pAxisA = makePAxis(axistype_a, var_axis_a);
+    ProtoAxis pAxisB = makePAxis(axistype_b, var_axis_b);
+
+    if (axisbdt_a == "closed" && axisbdt_b == "closed") {
+      m_gridLookup = SurfaceArray::makeSurfaceGridLookup2D<
+          detail::AxisBoundaryType::Closed,
+          detail::AxisBoundaryType::Closed>(g2l, l2g, pAxisA, pAxisB);
+    } else if (axisbdt_a == "closed" && axisbdt_b == "bound") {
+      m_gridLookup = SurfaceArray::makeSurfaceGridLookup2D<
+          detail::AxisBoundaryType::Closed,
+          detail::AxisBoundaryType::Bound>(g2l, l2g, pAxisA, pAxisB);
+    } else if (axisbdt_a == "closed" && axisbdt_b == "open") {
+      m_gridLookup = SurfaceArray::makeSurfaceGridLookup2D<
+          detail::AxisBoundaryType::Closed,
+          detail::AxisBoundaryType::Open>(g2l, l2g, pAxisA, pAxisB);
+    } else if (axisbdt_a == "open" && axisbdt_b == "closed") {
+      m_gridLookup = SurfaceArray::makeSurfaceGridLookup2D<
+          detail::AxisBoundaryType::Open,
+          detail::AxisBoundaryType::Closed>(g2l, l2g, pAxisA, pAxisB);
+    } else if (axisbdt_a == "open" && axisbdt_b == "bound") {
+      m_gridLookup = SurfaceArray::makeSurfaceGridLookup2D<
+          detail::AxisBoundaryType::Open,
+          detail::AxisBoundaryType::Bound>(g2l, l2g, pAxisA, pAxisB);
+    } else if (axisbdt_a == "open" && axisbdt_b == "open") {
+      m_gridLookup = SurfaceArray::makeSurfaceGridLookup2D<
+          detail::AxisBoundaryType::Open,
+          detail::AxisBoundaryType::Open>(g2l, l2g, pAxisA, pAxisB);
+    } else if (axisbdt_a == "bound" && axisbdt_b == "closed") {
+      m_gridLookup = SurfaceArray::makeSurfaceGridLookup2D<
+          detail::AxisBoundaryType::Bound,
+          detail::AxisBoundaryType::Closed>(g2l, l2g, pAxisA, pAxisB);
+    } else if (axisbdt_a == "bound" && axisbdt_b == "bound") {
+      m_gridLookup = SurfaceArray::makeSurfaceGridLookup2D<
+          detail::AxisBoundaryType::Bound,
+          detail::AxisBoundaryType::Bound>(g2l, l2g, pAxisA, pAxisB);
+    } else if (axisbdt_a == "bound" && axisbdt_b == "open") {
+      m_gridLookup = SurfaceArray::makeSurfaceGridLookup2D<
+          detail::AxisBoundaryType::Bound,
+          detail::AxisBoundaryType::Open>(g2l, l2g, pAxisA, pAxisB);
+    }
+
+    m_gridLookup->fill(surfaces);
+  }
+
+  SurfaceArray(const variant_data&                                   data_,
+               std::function<std::array<double, 1>(const Vector3D&)> g2l,
+               std::function<Vector3D(const std::array<double, 1>&)> l2g)
+    : m_gridLookup(nullptr)
+  {
+    const variant_map& data = boost::get<variant_map>(data_);
+    throw_assert(data.get<std::string>("type") == "SurfaceArray",
+                 "Type must be SurfaceArray");
+
+    const variant_map& payload = data.get<variant_map>("payload");
+    const variant_map& var_sgl = payload.get<variant_map>("surfacegridlookup");
+    throw_assert(var_sgl.get<std::string>("type") == "SurfaceGridLookup",
+                 "Currently only SurfaceGridLookup can be deserialized");
+    const variant_vector& var_surfaces
+        = payload.get<variant_vector>("surfaces");
+
+    InstanceFactory factory;
+
+    std::vector<const Surface*> surfaces;
+    for (const auto& var_srf_ : var_surfaces) {
+      const variant_map& var_srf = boost::get<variant_map>(var_srf_);
+      const Surface*     surface
+          = factory.surface(var_srf.get<std::string>("type"), var_srf);
+      surfaces.push_back(surface);
+    }
+
+    m_surfaces = surfaces;
+
+    // reproduce axes
+    const variant_vector& var_axes
+        = var_sgl.get<variant_map>("payload").get<variant_vector>("axes");
+
+    throw_assert(var_axes.size() == 1,
+                 "This constructor cannot deserialize DIM!=2 data");
+
+    const variant_map& var_axis
+        = var_axes.get<variant_map>(0).get<variant_map>("payload");
+
+    std::string axistype = var_axis.get<std::string>("axistype");
+    std::string axisbdt  = var_axis.get<std::string>("axisboundarytype");
+
+    if (axistype == "equidistant" && axisbdt == "bound") {
+      detail::Axis<detail::AxisType::Equidistant,
+                   detail::AxisBoundaryType::Bound>
+      axis(var_axis.get<double>("min"),
+           var_axis.get<double>("max"),
+           var_axis.get<int>("nbins"));
+
+      m_gridLookup = std::make_unique<SurfaceGridLookup<decltype(axis)>>(
+          g2l, l2g, std::make_tuple(axis));
+    } else if (axistype == "equidistant" && axisbdt == "closed") {
+      detail::Axis<detail::AxisType::Equidistant,
+                   detail::AxisBoundaryType::Closed>
+      axis(var_axis.get<double>("min"),
+           var_axis.get<double>("max"),
+           var_axis.get<int>("nbins"));
+
+      m_gridLookup = std::make_unique<SurfaceGridLookup<decltype(axis)>>(
+          g2l, l2g, std::make_tuple(axis));
+    } else if (axistype == "equidistant" && axisbdt == "open") {
+      detail::Axis<detail::AxisType::Equidistant,
+                   detail::AxisBoundaryType::Open>
+      axis(var_axis.get<double>("min"),
+           var_axis.get<double>("max"),
+           var_axis.get<int>("nbins"));
+
+      m_gridLookup = std::make_unique<SurfaceGridLookup<decltype(axis)>>(
+          g2l, l2g, std::make_tuple(axis));
+    } else if (axistype == "variable") {
+
+      std::vector<double>   bin_edges;
+      const variant_vector& var_bin_edges
+          = var_axis.get<variant_vector>("bin_edges");
+      for (size_t i = 0; i < var_bin_edges.size(); i++) {
+        bin_edges.push_back(var_bin_edges.get<double>(i));
+      }
+
+      if (axisbdt == "bound") {
+        detail::Axis<detail::AxisType::Variable,
+                     detail::AxisBoundaryType::Bound>
+            axis(bin_edges);
+        m_gridLookup = std::make_unique<SurfaceGridLookup<decltype(axis)>>(
+            g2l, l2g, std::make_tuple(axis));
+      } else if (axisbdt == "closed") {
+        detail::Axis<detail::AxisType::Variable,
+                     detail::AxisBoundaryType::Closed>
+            axis(bin_edges);
+        m_gridLookup = std::make_unique<SurfaceGridLookup<decltype(axis)>>(
+            g2l, l2g, std::make_tuple(axis));
+      } else if (axisbdt == "open") {
+        detail::Axis<detail::AxisType::Variable, detail::AxisBoundaryType::Open>
+            axis(bin_edges);
+        m_gridLookup = std::make_unique<SurfaceGridLookup<decltype(axis)>>(
+            g2l, l2g, std::make_tuple(axis));
+      }
+    }
+
+    m_gridLookup->fill(surfaces);
+  }
+
   /// @brief Get all surfaces in bin given by position.
   /// @param pos the lookup position
   /// @return reference to @c SurfaceVector contained in bin at that position
@@ -615,6 +900,12 @@ public:
     return p_gridLookup->isValidBin(bin);
   }
 
+  const Transform3D&
+  transform() const
+  {
+    return *m_transform;
+  }
+
   /// @brief String representation of this @c SurfaceArray
   /// @param sl Output stream to write to
   /// @return the output stream given as @p sl
@@ -650,9 +941,101 @@ public:
     return sl;
   }
 
+  variant_data
+  toVariantData() const
+  {
+    using namespace std::string_literals;
+    variant_map payload;
+
+    payload["surfacegridlookup"] = m_gridLookup->toVariantData();
+
+    variant_vector surfaces;
+
+    for (const auto& srf : m_surfaces) {
+      surfaces.push_back(srf->toVariantData());
+    }
+
+    payload["surfaces"] = surfaces;
+
+    variant_map data;
+    data["type"]    = "SurfaceArray"s;
+    data["payload"] = payload;
+    return data;
+  }
+
 private:
   std::unique_ptr<ISurfaceGridLookup> p_gridLookup;
   SurfaceVector                       m_surfaces;
+  // this is only used to keep info on transform applied
+  // by l2g and g2l
+  std::shared_ptr<const Transform3D> m_transform;
+
+  struct ProtoAxis
+  {
+    BinningType         bType;
+    BinningValue        bValue;
+    size_t              nBins;
+    double              min;
+    double              max;
+    std::vector<double> binEdges;
+  };
+
+  template <detail::AxisBoundaryType bdtA,
+            detail::AxisBoundaryType bdtB,
+            typename F1,
+            typename F2>
+  static std::unique_ptr<SurfaceArray::ISurfaceGridLookup>
+  makeSurfaceGridLookup2D(F1                      globalToLocal,
+                          F2                      localToGlobal,
+                          SurfaceArray::ProtoAxis pAxisA,
+                          SurfaceArray::ProtoAxis pAxisB)
+  {
+
+    using ISGL = SurfaceArray::ISurfaceGridLookup;
+    std::unique_ptr<ISGL> ptr;
+
+    // this becomes completely unreadable otherwise
+    // clang-format off
+    if (pAxisA.bType == equidistant && pAxisB.bType == equidistant) {
+
+      detail::Axis<detail::AxisType::Equidistant, bdtA> axisA(pAxisA.min, pAxisA.max, pAxisA.nBins);
+      detail::Axis<detail::AxisType::Equidistant, bdtB> axisB(pAxisB.min, pAxisB.max, pAxisB.nBins);
+
+      using SGL = SurfaceArray::SurfaceGridLookup<decltype(axisA), decltype(axisB)>;
+      ptr = std::unique_ptr<ISGL>(static_cast<ISGL*>(
+            new SGL(globalToLocal, localToGlobal, std::make_tuple(axisA, axisB))));
+
+    } else if (pAxisA.bType == equidistant && pAxisB.bType == arbitrary) {
+
+      detail::Axis<detail::AxisType::Equidistant, bdtA> axisA(pAxisA.min, pAxisA.max, pAxisA.nBins);
+      detail::Axis<detail::AxisType::Variable, bdtB> axisB(pAxisB.binEdges);
+
+      using SGL = SurfaceArray::SurfaceGridLookup<decltype(axisA), decltype(axisB)>;
+      ptr = std::unique_ptr<ISGL>(static_cast<ISGL*>(
+            new SGL(globalToLocal, localToGlobal, std::make_tuple(axisA, axisB))));
+      
+    } else if (pAxisA.bType == arbitrary && pAxisB.bType == equidistant) {
+
+      detail::Axis<detail::AxisType::Variable, bdtA> axisA(pAxisA.binEdges);
+      detail::Axis<detail::AxisType::Equidistant, bdtB> axisB(pAxisB.min, pAxisB.max, pAxisB.nBins);
+
+      using SGL = SurfaceArray::SurfaceGridLookup<decltype(axisA), decltype(axisB)>;
+      ptr = std::unique_ptr<ISGL>(static_cast<ISGL*>(
+            new SGL(globalToLocal, localToGlobal, std::make_tuple(axisA, axisB))));
+
+    } else /*if (pAxisA.bType == arbitrary && pAxisB.bType == arbitrary)*/ {
+
+      detail::Axis<detail::AxisType::Variable, bdtA> axisA(pAxisA.binEdges);
+      detail::Axis<detail::AxisType::Variable, bdtB> axisB(pAxisB.binEdges);
+
+      using SGL = SurfaceArray::SurfaceGridLookup<decltype(axisA), decltype(axisB)>;
+      ptr = std::unique_ptr<ISGL>(static_cast<ISGL*>(
+            new SGL(globalToLocal, localToGlobal, std::make_tuple(axisA, axisB))));
+    }
+    // clang-format on
+
+    return ptr;
+  }
 };
 
 }  // namespace Acts
diff --git a/Core/include/ACTS/Surfaces/SurfaceBounds.hpp b/Core/include/ACTS/Surfaces/SurfaceBounds.hpp
index 5aa219c7f..ceada2186 100644
--- a/Core/include/ACTS/Surfaces/SurfaceBounds.hpp
+++ b/Core/include/ACTS/Surfaces/SurfaceBounds.hpp
@@ -19,7 +19,6 @@
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/VariantDataFwd.hpp"
 
-
 namespace Acts {
 
 /// @class SurfaceBounds
@@ -98,10 +97,10 @@ public:
   /// @param sl is the outstream in which the string dump is done
   virtual std::ostream&
   dump(std::ostream& os) const = 0;
-  
-  //virtual 
-  //variant_data
-  //toVariantData() const = 0;
+
+  // virtual
+  // variant_data
+  // toVariantData() const = 0;
 };
 
 inline bool
diff --git a/Core/include/ACTS/Surfaces/TrapezoidBounds.hpp b/Core/include/ACTS/Surfaces/TrapezoidBounds.hpp
index 646171151..cbd74f2ea 100644
--- a/Core/include/ACTS/Surfaces/TrapezoidBounds.hpp
+++ b/Core/include/ACTS/Surfaces/TrapezoidBounds.hpp
@@ -52,7 +52,7 @@ public:
   /// @param haley half length Y - defined at x=0
   TrapezoidBounds(double minhalex, double maxhalex, double haley);
 
-  TrapezoidBounds(const variant_data &data);
+  TrapezoidBounds(const variant_data& data);
 
   virtual ~TrapezoidBounds();
 
@@ -148,7 +148,7 @@ public:
   /// (second coordinate of local surface frame)
   double
   halflengthY() const;
-  
+
   variant_data
   toVariantData() const override;
 
diff --git a/Core/include/ACTS/Surfaces/TriangleBounds.hpp b/Core/include/ACTS/Surfaces/TriangleBounds.hpp
index 5fc123ae8..53d609adb 100644
--- a/Core/include/ACTS/Surfaces/TriangleBounds.hpp
+++ b/Core/include/ACTS/Surfaces/TriangleBounds.hpp
@@ -13,8 +13,8 @@
 #ifndef ACTS_SURFACESTRIANGLEBOUNDS_H
 #define ACTS_SURFACESTRIANGLEBOUNDS_H
 
-#include <utility>
 #include <array>
+#include <utility>
 
 #include "ACTS/Surfaces/PlanarBounds.hpp"
 #include "ACTS/Surfaces/RectangleBounds.hpp"
@@ -51,7 +51,7 @@ public:
   /// @param vertices is the vector of vertices
   TriangleBounds(const std::array<Vector2D, 3>& vertices);
 
-  TriangleBounds(const variant_data &data);
+  TriangleBounds(const variant_data& data);
 
   virtual ~TriangleBounds();
 
@@ -94,12 +94,12 @@ public:
   /// @param sl is the ostream to be dumped into
   virtual std::ostream&
   dump(std::ostream& sl) const final override;
-  
+
   variant_data
   toVariantData() const override;
 
 private:
-  std::array<Vector2D, 3>       m_vertices;
+  std::array<Vector2D, 3> m_vertices;
   RectangleBounds m_boundingBox;  ///< internal bounding box cache
 };
 
diff --git a/Core/include/ACTS/Utilities/InstanceFactory.hpp b/Core/include/ACTS/Utilities/InstanceFactory.hpp
index 876b186bf..71792b86c 100644
--- a/Core/include/ACTS/Utilities/InstanceFactory.hpp
+++ b/Core/include/ACTS/Utilities/InstanceFactory.hpp
@@ -11,21 +11,24 @@
 
 #include <boost/functional/value_factory.hpp>
 #include <map>
-#include "ACTS/Utilities/VariantData.hpp"
 #include "ACTS/Utilities/ThrowAssert.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 
-#include "ACTS/Surfaces/RectangleBounds.hpp"
-#include "ACTS/Surfaces/TriangleBounds.hpp"
+#include "ACTS/Surfaces/CylinderBounds.hpp"
 #include "ACTS/Surfaces/DiamondBounds.hpp"
 #include "ACTS/Surfaces/EllipseBounds.hpp"
-#include "ACTS/Surfaces/CylinderBounds.hpp"
+#include "ACTS/Surfaces/PlaneSurface.hpp"
+#include "ACTS/Surfaces/RectangleBounds.hpp"
 #include "ACTS/Surfaces/TrapezoidBounds.hpp"
+#include "ACTS/Surfaces/TriangleBounds.hpp"
 
 namespace Acts {
 
-using SurfaceBoundsPtr     = std::shared_ptr<const SurfaceBounds>;
-using PlanarBoundsPtr     = std::shared_ptr<const PlanarBounds>;
-using SurfaceBoundsFactory = std::function<SurfaceBoundsPtr(const variant_data&)>;
+using SurfaceBoundsPtr = std::shared_ptr<const SurfaceBounds>;
+using PlanarBoundsPtr  = std::shared_ptr<const PlanarBounds>;
+using SurfaceBoundsFactory
+    = std::function<SurfaceBoundsPtr(const variant_data&)>;
+using SurfaceFactory = std::function<const Surface*(const variant_data&)>;
 
 class InstanceFactory
 {
@@ -52,13 +55,16 @@ public:
       return std::make_shared<const CylinderBounds>(data);
     };
 
+    m_surfaces["PlaneSurface"]
+        = [](auto const& data) { return new PlaneSurface(data); };
   }
 
   PlanarBoundsPtr
   planarBounds(const std::string& cname, const variant_data& data)
   {
     SurfaceBoundsPtr srfBnd_ptr = surfaceBounds(cname, data);
-    PlanarBoundsPtr plnBnd_ptr = std::dynamic_pointer_cast<const PlanarBounds>(srfBnd_ptr);
+    PlanarBoundsPtr  plnBnd_ptr
+        = std::dynamic_pointer_cast<const PlanarBounds>(srfBnd_ptr);
     throw_assert(plnBnd_ptr, "Conversion to PlanarBounds failed");
     return plnBnd_ptr;
   }
@@ -66,12 +72,22 @@ public:
   SurfaceBoundsPtr
   surfaceBounds(const std::string& cname, const variant_data& data)
   {
-    throw_assert(m_surfaceBounds.count(cname), "No factory found for class "+cname);
+    throw_assert(m_surfaceBounds.count(cname),
+                 "No factory found for class " + cname);
     return m_surfaceBounds.at(cname)(data);
   }
 
+  const Surface*
+  surface(const std::string& cname, const variant_data& data)
+  {
+    throw_assert(m_surfaces.count(cname),
+                 "No factory found for class " + cname);
+    return m_surfaces.at(cname)(data);
+  }
+
 private:
   std::map<std::string, SurfaceBoundsFactory> m_surfaceBounds;
+  std::map<std::string, SurfaceFactory>       m_surfaces;
 };
 
 }  // namespace Acts
diff --git a/Core/include/ACTS/Utilities/VariantData.hpp b/Core/include/ACTS/Utilities/VariantData.hpp
index 44c171f68..38ddfe91e 100644
--- a/Core/include/ACTS/Utilities/VariantData.hpp
+++ b/Core/include/ACTS/Utilities/VariantData.hpp
@@ -68,6 +68,8 @@ public:
   T&
   get(const std::string& key)
   {
+    if (!m_map.count(key))
+      throw std::out_of_range("variant_map key " + key + " not found");
     return boost::get<T>(m_map.at(key));
   }
 
@@ -75,6 +77,8 @@ public:
   const T&
   get(const std::string& key) const
   {
+    if (!m_map.count(key))
+      throw std::out_of_range("variant_map key " + key + " not found");
     return boost::get<T>(m_map.at(key));
   }
 
@@ -142,7 +146,7 @@ public:
   {
     m_vector.push_back(data);
   }
-  
+
   template <class T>
   T&
   get(const size_t& idx)
@@ -359,11 +363,11 @@ to_variant(const ActsMatrixD<4, 4>& matrix)
   return data;
 }
 
-template<typename T>
+template <typename T>
 inline T
 from_variant(const variant_data& data_);
 
-template<>
+template <>
 inline Transform3D
 from_variant<Transform3D>(const variant_data& data_)
 {
@@ -372,14 +376,14 @@ from_variant<Transform3D>(const variant_data& data_)
   throw_assert(data.get<std::string>("type") == "Transform3D",
                "Must be type Transform3D");
 
-  const variant_map &payload = data.get<variant_map>("payload");
+  const variant_map& payload = data.get<variant_map>("payload");
 
-  const variant_vector &matrix_data = payload.get<variant_vector>("data");
-  Transform3D trf;
-  for(size_t i=0;i<4;i++) {
-    for(size_t j=0;j<4;j++) {
+  const variant_vector& matrix_data = payload.get<variant_vector>("data");
+  Transform3D           trf;
+  for (size_t i = 0; i < 4; i++) {
+    for (size_t j = 0; j < 4; j++) {
 
-      size_t k = i*4+j;
+      size_t k     = i * 4 + j;
       double value = matrix_data.get<double>(k);
       trf(i, j) = value;
     }
@@ -388,7 +392,7 @@ from_variant<Transform3D>(const variant_data& data_)
   return trf;
 }
 
-template<>
+template <>
 inline Vector2D
 from_variant<Vector2D>(const variant_data& data_)
 {
@@ -397,10 +401,10 @@ from_variant<Vector2D>(const variant_data& data_)
   throw_assert(data.get<std::string>("type") == "Vector2D",
                "Must be type Vector2D");
 
-  const variant_vector &vector_data = data.get<variant_vector>("payload");
+  const variant_vector& vector_data = data.get<variant_vector>("payload");
 
   Vector2D vec;
-  for(size_t i=0;i<2;i++) {
+  for (size_t i = 0; i < 2; i++) {
     vec[i] = vector_data.get<double>(i);
   }
 
diff --git a/Core/include/ACTS/Utilities/VariantDataFwd.hpp b/Core/include/ACTS/Utilities/VariantDataFwd.hpp
index ba9ee6c44..90815390c 100644
--- a/Core/include/ACTS/Utilities/VariantDataFwd.hpp
+++ b/Core/include/ACTS/Utilities/VariantDataFwd.hpp
@@ -9,13 +9,13 @@
 #ifndef ACTS_UTILITIES_VARIANTDATA_FWD_H
 #define ACTS_UTILITIES_VARIANTDATA_FWD_H 1
 
-#include <boost/variant/variant_fwd.hpp>
 #include <boost/variant/recursive_wrapper.hpp>
+#include <boost/variant/variant_fwd.hpp>
 #include <iostream>
 #include <map>
 #include <sstream>
-#include <vector>
 #include <string>
+#include <vector>
 
 namespace Acts {
 
@@ -23,12 +23,10 @@ class variant_map;
 class variant_vector;
 
 using variant_data = boost::variant<int,
-                                     double,
-                                     std::string,
-                                     bool,
-                                     boost::recursive_wrapper<variant_map>,
-                                     boost::recursive_wrapper<variant_vector>>;
-
-
+                                    double,
+                                    std::string,
+                                    bool,
+                                    boost::recursive_wrapper<variant_map>,
+                                    boost::recursive_wrapper<variant_vector>>;
 }
 #endif
diff --git a/Core/src/Layers/CylinderLayer.cpp b/Core/src/Layers/CylinderLayer.cpp
index 5726b6e17..4eff2e998 100644
--- a/Core/src/Layers/CylinderLayer.cpp
+++ b/Core/src/Layers/CylinderLayer.cpp
@@ -45,6 +45,76 @@ Acts::CylinderLayer::CylinderLayer(
   if (m_approachDescriptor) approachDescriptor()->registerLayer(*this);
 }
 
+std::shared_ptr<Acts::Layer>
+Acts::CylinderLayer::create(const variant_data& data_)
+{
+  throw_assert(data_.which() == 4, "Variant data must be map");
+  const variant_map& data = boost::get<variant_map>(data_);
+  std::string        type = data.get<std::string>("type");
+  throw_assert(type == "CylinderLayer", "Type must be CylinderLayer");
+
+  variant_map payload = data.get<variant_map>("payload");
+
+  auto trf = std::make_shared<const Transform3D>(
+      from_variant<Transform3D>(payload.get<variant_map>("transform")));
+
+  LayerType   laytyp;
+  std::string laytyp_str = payload.get<std::string>("layer_type");
+  if (laytyp_str == "active")
+    laytyp = active;
+  else if (laytyp_str == "passive")
+    laytyp = passive;
+  else /*laytyp_str == "navigation"*/
+    laytyp = navigation;
+
+  double thickness = payload.get<double>("thickness");
+
+  InstanceFactory    factory;
+  const variant_map& var_bounds = payload.get<variant_map>("cylinder_bounds");
+
+  auto cbounds = std::dynamic_pointer_cast<const CylinderBounds>(
+      factory.surfaceBounds(var_bounds.get<std::string>("type"), var_bounds));
+
+  std::unique_ptr<SurfaceArray> sArray = nullptr;
+
+  // only attempt to reover surface array if present
+  if (payload.count("surfacearray")) {
+
+    // get surface array transform
+    const Transform3D& sa_trf = from_variant<Transform3D>(
+        payload.get<variant_map>("surfacearray_transform"));
+    const Transform3D& sa_itrf = sa_trf.inverse();
+
+    // we need to reproduce the coordinate conversions
+    double R      = cbounds->r();
+    double avgPhi = cbounds->averagePhi();
+    auto g2l      = [sa_trf, avgPhi](const Vector3D& pos) -> Vector2D {
+      // @TODO: Check if - is right here, might be the other way round
+      Vector3D loc = sa_trf * pos;
+      return Vector2D(loc.phi() - avgPhi, loc.z());
+    };
+    auto l2g = [sa_itrf, R, avgPhi](const Vector2D& loc) -> Vector3D {
+      return sa_itrf * Vector3D(R * std::cos(loc[0] + avgPhi),
+                                R * std::sin(loc[0] + avgPhi),
+                                loc[1]);
+    };
+
+    sArray = std::make_unique<SurfaceArray>(
+        payload.at("surfacearray"),
+        g2l,
+        l2g,
+        std::make_shared<const Transform3D>(sa_trf));
+  }
+
+  // @TODO: Implement ApproachDescriptor serialization
+  return MutableLayerPtr(new CylinderLayer(trf,
+                                           cbounds,
+                                           std::move(sArray),
+                                           thickness,
+                                           nullptr,  // std::move(ad),
+                                           laytyp));
+}
+
 const Acts::CylinderSurface&
 Acts::CylinderLayer::surfaceRepresentation() const
 {
@@ -100,3 +170,45 @@ Acts::CylinderLayer::buildApproachDescriptor()
     }
   }
 }
+
+Acts::variant_data
+Acts::CylinderLayer::toVariantData() const
+{
+  using namespace std::string_literals;
+  variant_map payload;
+
+  payload["transform"] = to_variant(*m_transform);
+
+  // we need to recover the bounds
+  const AbstractVolume* absVol = representingVolume();
+  auto                  cvBounds
+      = dynamic_cast<const CylinderVolumeBounds*>(&absVol->volumeBounds());
+
+  double cylR = cvBounds->innerRadius() + 0.5 * thickness();
+  double hlZ  = cvBounds->halflengthZ();
+
+  CylinderBounds cylBounds(cylR, hlZ);
+
+  const variant_data bounds  = cylBounds.toVariantData();
+  payload["cylinder_bounds"] = bounds;
+
+  payload["thickness"] = thickness();
+
+  if (layerType() == active)
+    payload["layer_type"] = "active"s;
+  else if (layerType() == passive)
+    payload["layer_type"] = "passive"s;
+  else /*layerType() == navigation*/
+    payload["layer_type"] = "navigation"s;
+
+  // this lacks localToGlobal and globalToLocal
+  if (m_surfaceArray) {
+    payload["surfacearray"]           = m_surfaceArray->toVariantData();
+    payload["surfacearray_transform"] = to_variant(m_surfaceArray->transform());
+  }
+
+  variant_map data;
+  data["type"]    = "CylinderLayer"s;
+  data["payload"] = payload;
+  return data;
+}
diff --git a/Core/src/Layers/DiscLayer.cpp b/Core/src/Layers/DiscLayer.cpp
index 2ef760a1b..a1f505223 100644
--- a/Core/src/Layers/DiscLayer.cpp
+++ b/Core/src/Layers/DiscLayer.cpp
@@ -12,11 +12,13 @@
 
 #include "ACTS/Layers/DiscLayer.hpp"
 #include "ACTS/Layers/GenericApproachDescriptor.hpp"
+#include "ACTS/Layers/Layer.hpp"
 #include "ACTS/Material/SurfaceMaterial.hpp"
 #include "ACTS/Surfaces/DiscBounds.hpp"
 #include "ACTS/Surfaces/RadialBounds.hpp"
 #include "ACTS/Utilities/BinUtility.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 #include "ACTS/Volumes/AbstractVolume.hpp"
 #include "ACTS/Volumes/BoundarySurfaceFace.hpp"
 #include "ACTS/Volumes/CylinderVolumeBounds.hpp"
@@ -48,6 +50,72 @@ Acts::DiscLayer::DiscLayer(std::shared_ptr<const Transform3D>  transform,
   if (m_approachDescriptor) approachDescriptor()->registerLayer(*this);
 }
 
+std::shared_ptr<Acts::Layer>
+Acts::DiscLayer::create(const variant_data& data_)
+{
+  throw_assert(data_.which() == 4, "Variant data must be map");
+  const variant_map& data = boost::get<variant_map>(data_);
+  std::string        type = data.get<std::string>("type");
+  throw_assert(type == "DiscLayer", "Type must be DiscLayer");
+
+  variant_map payload = data.get<variant_map>("payload");
+
+  auto trf = std::make_shared<const Transform3D>(
+      from_variant<Transform3D>(payload.get<variant_map>("transform")));
+
+  LayerType   laytyp;
+  std::string laytyp_str = payload.get<std::string>("layer_type");
+  if (laytyp_str == "active")
+    laytyp = active;
+  else if (laytyp_str == "passive")
+    laytyp = passive;
+  else /*laytyp_str == "navigation"*/
+    laytyp = navigation;
+
+  double thickness = payload.get<double>("thickness");
+  double minR      = payload.get<double>("minR");
+  double maxR      = payload.get<double>("maxR");
+  double R         = (minR + maxR) / 2.;
+
+  auto rbounds = std::make_shared<const RadialBounds>(minR, maxR);
+
+  std::unique_ptr<SurfaceArray> sArray = nullptr;
+
+  // only attempt to reover surface array if present
+  if (payload.count("surfacearray")) {
+
+    // get surface array transform
+    const Transform3D& sa_trf = from_variant<Transform3D>(
+        payload.get<variant_map>("surfacearray_transform"));
+    const Transform3D& sa_itrf = sa_trf.inverse();
+
+    // we need to reproduce the coordinate conversions
+    auto g2l = [sa_trf](const Vector3D& pos) -> Vector2D {
+      // @TODO: Check if - is right here, might be the other way round
+      Vector3D loc = sa_trf * pos;
+      return Vector2D(loc.perp(), loc.phi());
+    };
+    auto l2g = [sa_itrf, R](const Vector2D& loc) -> Vector3D {
+      return sa_itrf
+          * Vector3D(R * std::cos(loc[0]), R * std::sin(loc[0]), loc[1]);
+    };
+
+    sArray = std::make_unique<SurfaceArray>(
+        payload.at("surfacearray"),
+        g2l,
+        l2g,
+        std::make_shared<const Transform3D>(sa_trf));
+  }
+
+  // @TODO: Implement ApproachDescriptor serialization
+  return MutableLayerPtr(new DiscLayer(trf,
+                                       rbounds,
+                                       std::move(sArray),
+                                       thickness,
+                                       nullptr,  // std::move(ad),
+                                       laytyp));
+}
+
 const Acts::DiscSurface&
 Acts::DiscLayer::surfaceRepresentation() const
 {
@@ -92,8 +160,8 @@ Acts::DiscLayer::buildApproachDescriptor()
     std::vector<const Surface*> aSurfaces;
     aSurfaces.push_back(new DiscSurface(asnTransform, m_bounds));
     aSurfaces.push_back(new DiscSurface(aspTransform, m_bounds));
-    // create an ApproachDescriptor with standard surfaces surfaces - these will
-    // be deleted by the approach descriptor
+    // create an ApproachDescriptor with standard surfaces surfaces - these
+    // will be deleted by the approach descriptor
     m_approachDescriptor
         = std::make_unique<const GenericApproachDescriptor<Surface>>(aSurfaces);
   }
@@ -105,3 +173,40 @@ Acts::DiscLayer::buildApproachDescriptor()
     }
   }
 }
+
+Acts::variant_data
+Acts::DiscLayer::toVariantData() const
+{
+  using namespace std::string_literals;
+  variant_map payload;
+
+  payload["transform"] = to_variant(*m_transform);
+
+  // we need to recover the bounds
+  const AbstractVolume* absVol = representingVolume();
+  throw_assert(absVol,
+               "Cannot serialize DiscLayer without representing volume");
+  auto cvBounds
+      = dynamic_cast<const CylinderVolumeBounds*>(&absVol->volumeBounds());
+
+  payload["minR"]      = cvBounds->innerRadius();
+  payload["maxR"]      = cvBounds->outerRadius();
+  payload["thickness"] = thickness();
+
+  if (layerType() == active)
+    payload["layer_type"] = "active"s;
+  else if (layerType() == passive)
+    payload["layer_type"] = "passive"s;
+  else /*layerType() == navigation*/
+    payload["layer_type"] = "navigation"s;
+
+  if (m_surfaceArray) {
+    payload["surfacearray"]           = m_surfaceArray->toVariantData();
+    payload["surfacearray_transform"] = to_variant(m_surfaceArray->transform());
+  }
+
+  variant_map data;
+  data["type"]    = "DiscLayer"s;
+  data["payload"] = payload;
+  return data;
+}
diff --git a/Core/src/Surfaces/CylinderBounds.cpp b/Core/src/Surfaces/CylinderBounds.cpp
index e4801365c..6ad646b07 100644
--- a/Core/src/Surfaces/CylinderBounds.cpp
+++ b/Core/src/Surfaces/CylinderBounds.cpp
@@ -16,8 +16,8 @@
 #include <iomanip>
 #include <iostream>
 
-#include "ACTS/Utilities/detail/periodic.hpp"
 #include "ACTS/Utilities/VariantData.hpp"
+#include "ACTS/Utilities/detail/periodic.hpp"
 
 Acts::CylinderBounds::CylinderBounds(double radius, double halfZ)
   : CylinderBounds(radius, 0, M_PI, halfZ)
@@ -46,16 +46,16 @@ Acts::CylinderBounds::CylinderBounds(const variant_data& data_)
   : m_radius(0), m_avgPhi(0), m_halfPhi(0), m_halfZ(0)
 {
   throw_assert(data_.which() == 4, "Variant data must be map");
-  const variant_map &data = boost::get<variant_map>(data_);
-  std::string type = data.get<std::string>("type");
+  const variant_map& data = boost::get<variant_map>(data_);
+  std::string        type = data.get<std::string>("type");
   throw_assert(type == "CylinderBounds", "Type must be CylinderBounds");
 
-  const variant_map &payload = data.get<variant_map>("payload");
+  const variant_map& payload = data.get<variant_map>("payload");
 
-  m_radius = payload.get<double>("radius");
-  m_avgPhi = payload.get<double>("avgPhi");
+  m_radius  = payload.get<double>("radius");
+  m_avgPhi  = payload.get<double>("avgPhi");
   m_halfPhi = payload.get<double>("halfPhi");
-  m_halfZ = payload.get<double>("halfZ");
+  m_halfZ   = payload.get<double>("halfZ");
 }
 
 Acts::CylinderBounds::~CylinderBounds()
@@ -153,16 +153,12 @@ Acts::CylinderBounds::toVariantData() const
   using namespace std::string_literals;
 
   variant_map payload;
-  payload["radius"] = m_radius;
-  payload["avgPhi"] = m_avgPhi;
+  payload["radius"]  = m_radius;
+  payload["avgPhi"]  = m_avgPhi;
   payload["halfPhi"] = m_halfPhi;
-  payload["halfZ"] = m_halfZ;
+  payload["halfZ"]   = m_halfZ;
 
-  variant_map data({
-    {"type", "CylinderBounds"s},
-    {"payload", payload}
-  });
+  variant_map data({{"type", "CylinderBounds"s}, {"payload", payload}});
 
-  
   return data;
 }
diff --git a/Core/src/Surfaces/DiamondBounds.cpp b/Core/src/Surfaces/DiamondBounds.cpp
index ff9025290..f0e87c7c7 100644
--- a/Core/src/Surfaces/DiamondBounds.cpp
+++ b/Core/src/Surfaces/DiamondBounds.cpp
@@ -56,7 +56,9 @@ Acts::DiamondBounds::DiamondBounds(const variant_data& data_)
                         std::max(m_minY, m_maxY));
 }
 
-Acts::DiamondBounds::~DiamondBounds() {}
+Acts::DiamondBounds::~DiamondBounds()
+{
+}
 
 Acts::DiamondBounds*
 Acts::DiamondBounds::clone() const
diff --git a/Core/src/Surfaces/EllipseBounds.cpp b/Core/src/Surfaces/EllipseBounds.cpp
index 19a778fa5..5c5c1074d 100644
--- a/Core/src/Surfaces/EllipseBounds.cpp
+++ b/Core/src/Surfaces/EllipseBounds.cpp
@@ -57,7 +57,9 @@ Acts::EllipseBounds::EllipseBounds(const variant_data& data_)
       = RectangleBounds(std::max(m_rMinX, m_rMaxX), std::max(m_rMinY, m_rMaxY));
 }
 
-Acts::EllipseBounds::~EllipseBounds() {}
+Acts::EllipseBounds::~EllipseBounds()
+{
+}
 
 Acts::EllipseBounds*
 Acts::EllipseBounds::clone() const
diff --git a/Core/src/Surfaces/LineBounds.cpp b/Core/src/Surfaces/LineBounds.cpp
index d34fca1bd..019d793f2 100644
--- a/Core/src/Surfaces/LineBounds.cpp
+++ b/Core/src/Surfaces/LineBounds.cpp
@@ -11,6 +11,7 @@
 ///////////////////////////////////////////////////////////////////
 
 #include "ACTS/Surfaces/LineBounds.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 
 #include <iomanip>
 #include <iostream>
@@ -20,6 +21,22 @@ Acts::LineBounds::LineBounds(double radius, double halez)
 {
 }
 
+Acts::LineBounds::LineBounds(const variant_data& data_)
+{
+
+  throw_assert(data_.which() == 4, "Variant data must be map");
+  variant_map data = boost::get<variant_map>(data_);
+  throw_assert(data.count("type"), "Variant data must have type.");
+  // std::string type = boost::get<std::string>(data["type"]);
+  std::string type = data.get<std::string>("type");
+  throw_assert(type == "LineBounds", "Variant data type must be LineBounds");
+
+  variant_map payload = data.get<variant_map>("payload");
+
+  m_radius = payload.get<double>("radius");
+  m_halfZ  = payload.get<double>("halfZ");
+}
+
 Acts::LineBounds::~LineBounds()
 {
 }
@@ -70,3 +87,18 @@ Acts::LineBounds::dump(std::ostream& sl) const
   sl << std::setprecision(-1);
   return sl;
 }
+
+Acts::variant_data
+Acts::LineBounds::toVariantData() const
+{
+  using namespace std::string_literals;
+  variant_map payload;
+
+  payload["radius"] = m_radius;
+  payload["halfZ"]  = m_halfZ;
+
+  variant_map data;
+  data["type"]    = "LineBounds"s;
+  data["payload"] = payload;
+  return data;
+}
diff --git a/Core/src/Surfaces/LineSurface.cpp b/Core/src/Surfaces/LineSurface.cpp
index 9a1cfe6bc..1d1fcc951 100644
--- a/Core/src/Surfaces/LineSurface.cpp
+++ b/Core/src/Surfaces/LineSurface.cpp
@@ -17,6 +17,7 @@
 #include <iostream>
 
 #include "ACTS/Utilities/ThrowAssert.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 
 Acts::LineSurface::LineSurface(std::shared_ptr<const Transform3D> htrans,
                                double                             radius,
@@ -52,6 +53,32 @@ Acts::LineSurface::LineSurface(const LineSurface& other,
 {
 }
 
+Acts::LineSurface::LineSurface(const variant_data& data_) : GeometryObject()
+{
+  throw_assert(data_.which() == 4, "Variant data must be map");
+  variant_map data = boost::get<variant_map>(data_);
+  throw_assert(data.count("type"), "Variant data must have type.");
+  // std::string type = boost::get<std::string>(data["type"]);
+  std::string type = data.get<std::string>("type");
+  throw_assert(type == "LineSurface", "Variant data type must be LineSurface");
+
+  variant_map payload    = data.get<variant_map>("payload");
+  variant_map bounds     = payload.get<variant_map>("bounds");
+  std::string boundsType = bounds.get<std::string>("type");
+
+  throw_assert(boundsType == "LineBounds",
+               "Can only construct LineSurface from LineBounds");
+
+  m_bounds = std::make_shared<const LineBounds>(bounds);
+
+  if (payload.count("transform")) {
+    // we have a transform
+    auto trf = std::make_shared<const Transform3D>(
+        from_variant<Transform3D>(payload.get<variant_map>("transform")));
+    m_transform = trf;
+  }
+}
+
 Acts::LineSurface::~LineSurface()
 {
 }
@@ -192,3 +219,23 @@ Acts::LineSurface::bounds() const
   if (m_bounds) return (*m_bounds.get());
   return s_noBounds;
 }
+
+Acts::variant_data
+Acts::LineSurface::toVariantData() const
+{
+  using namespace std::string_literals;
+
+  variant_map payload;
+
+  variant_data bounds = m_bounds->toVariantData();
+  payload["bounds"]   = bounds;
+
+  if (m_transform) {
+    payload["transform"] = to_variant(*m_transform);
+  }
+
+  variant_map data;
+  data["type"]    = "LineSurface"s;
+  data["payload"] = payload;
+  return data;
+}
diff --git a/Core/src/Surfaces/PerigeeSurface.cpp b/Core/src/Surfaces/PerigeeSurface.cpp
index 5191fcb3f..61cae0e67 100644
--- a/Core/src/Surfaces/PerigeeSurface.cpp
+++ b/Core/src/Surfaces/PerigeeSurface.cpp
@@ -11,11 +11,13 @@
 ///////////////////////////////////////////////////////////////////
 
 #include "ACTS/Surfaces/PerigeeSurface.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 
 #include <iomanip>
 #include <iostream>
 
-Acts::PerigeeSurface::PerigeeSurface(const Vector3D& gp) : LineSurface(nullptr)
+Acts::PerigeeSurface::PerigeeSurface(const Vector3D& gp)
+  : LineSurface(nullptr, nullptr)
 {
   Surface::m_transform = std::make_shared<const Transform3D>(
       Translation3D(gp.x(), gp.y(), gp.z()));
@@ -38,6 +40,28 @@ Acts::PerigeeSurface::PerigeeSurface(const PerigeeSurface& other,
 {
 }
 
+Acts::PerigeeSurface::PerigeeSurface(const variant_data& data_)
+  : GeometryObject(), LineSurface(nullptr, nullptr)
+{
+
+  throw_assert(data_.which() == 4, "Variant data must be map");
+  variant_map data = boost::get<variant_map>(data_);
+  throw_assert(data.count("type"), "Variant data must have type.");
+  // std::string type = boost::get<std::string>(data["type"]);
+  std::string type = data.get<std::string>("type");
+  throw_assert(type == "PerigeeSurface",
+               "Variant data type must be PerigeeSurface");
+
+  variant_map payload = data.get<variant_map>("payload");
+
+  if (payload.count("transform")) {
+    // we have a transform
+    auto trf = std::make_shared<const Transform3D>(
+        from_variant<Transform3D>(payload.get<variant_map>("transform")));
+    m_transform = trf;
+  }
+}
+
 Acts::PerigeeSurface::~PerigeeSurface()
 {
 }
@@ -81,3 +105,20 @@ Acts::PerigeeSurface::dump(std::ostream& sl) const
   sl << std::setprecision(-1);
   return sl;
 }
+
+Acts::variant_data
+Acts::PerigeeSurface::toVariantData() const
+{
+  using namespace std::string_literals;
+
+  variant_map payload;
+
+  if (m_transform) {
+    payload["transform"] = to_variant(*m_transform);
+  }
+
+  variant_map data;
+  data["type"]    = "PerigeeSurface"s;
+  data["payload"] = payload;
+  return data;
+}
diff --git a/Core/src/Surfaces/PlaneSurface.cpp b/Core/src/Surfaces/PlaneSurface.cpp
index 1484895e3..de763173f 100644
--- a/Core/src/Surfaces/PlaneSurface.cpp
+++ b/Core/src/Surfaces/PlaneSurface.cpp
@@ -102,7 +102,9 @@ Acts::PlaneSurface::PlaneSurface(const variant_data& data_)
   }
 }
 
-Acts::PlaneSurface::~PlaneSurface() {}
+Acts::PlaneSurface::~PlaneSurface()
+{
+}
 
 Acts::PlaneSurface&
 Acts::PlaneSurface::operator=(const PlaneSurface& other)
diff --git a/Core/src/Surfaces/RectangleBounds.cpp b/Core/src/Surfaces/RectangleBounds.cpp
index b2df2556a..66e8b3bfd 100644
--- a/Core/src/Surfaces/RectangleBounds.cpp
+++ b/Core/src/Surfaces/RectangleBounds.cpp
@@ -11,8 +11,8 @@
 ///////////////////////////////////////////////////////////////////
 
 #include "ACTS/Surfaces/RectangleBounds.hpp"
-#include "ACTS/Utilities/VariantData.hpp"
 #include "ACTS/Utilities/ThrowAssert.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 
 #include <cmath>
 #include <iomanip>
@@ -26,11 +26,11 @@ Acts::RectangleBounds::RectangleBounds(double halex, double haley)
 Acts::RectangleBounds::RectangleBounds(const variant_data& data_)
 {
   throw_assert(data_.which() == 4, "Variant data must be map");
-  const variant_map &data = boost::get<variant_map>(data_);
-  std::string type = data.get<std::string>("type");
+  const variant_map& data = boost::get<variant_map>(data_);
+  std::string        type = data.get<std::string>("type");
   throw_assert(type == "RectangleBounds", "Type must be RectangleBounds");
 
-  const variant_map &payload = data.get<variant_map>("payload");
+  const variant_map& payload = data.get<variant_map>("payload");
 
   m_halfX = payload.get<double>("halflengthX");
   m_halfY = payload.get<double>("halflengthY");
@@ -86,8 +86,6 @@ Acts::RectangleBounds::boundingBox() const
   return (*this);
 }
 
-
-
 // ostream operator overload
 std::ostream&
 Acts::RectangleBounds::dump(std::ostream& sl) const
@@ -109,11 +107,7 @@ Acts::RectangleBounds::toVariantData() const
   payload["halflengthX"] = m_halfX;
   payload["halflengthY"] = m_halfY;
 
-  variant_map data({
-    {"type", "RectangleBounds"s},
-    {"payload", payload}
-  });
+  variant_map data({{"type", "RectangleBounds"s}, {"payload", payload}});
 
-  
   return data;
 }
diff --git a/Core/src/Surfaces/TrapezoidBounds.cpp b/Core/src/Surfaces/TrapezoidBounds.cpp
index 831374138..1e30ef0a8 100644
--- a/Core/src/Surfaces/TrapezoidBounds.cpp
+++ b/Core/src/Surfaces/TrapezoidBounds.cpp
@@ -31,15 +31,15 @@ Acts::TrapezoidBounds::TrapezoidBounds(const variant_data& data_)
   : m_boundingBox(0, 0)
 {
   throw_assert(data_.which() == 4, "Variant data must be map");
-  const variant_map &data = boost::get<variant_map>(data_);
-  std::string type = data.get<std::string>("type");
+  const variant_map& data = boost::get<variant_map>(data_);
+  std::string        type = data.get<std::string>("type");
   throw_assert(type == "TrapezoidBounds", "Type must be TrapezoidBounds");
 
-  const variant_map &payload = data.get<variant_map>("payload");
+  const variant_map& payload = data.get<variant_map>("payload");
 
   m_minHalfX = payload.get<double>("minHalfX");
   m_maxHalfX = payload.get<double>("maxHalfX");
-  m_halfY = payload.get<double>("halfY");
+  m_halfY    = payload.get<double>("halfY");
 
   m_boundingBox = RectangleBounds(m_maxHalfX, m_halfY);
 }
@@ -112,16 +112,17 @@ Acts::TrapezoidBounds::dump(std::ostream& sl) const
 }
 
 Acts::variant_data
-Acts::TrapezoidBounds::toVariantData() const {
+Acts::TrapezoidBounds::toVariantData() const
+{
   using namespace std::string_literals;
 
   variant_map payload;
   payload["minHalfX"] = m_minHalfX;
   payload["maxHalfX"] = m_maxHalfX;
   payload["halfY"]    = m_halfY;
-  
+
   variant_map data;
-  data["type"] = "TrapezoidBounds"s;
+  data["type"]    = "TrapezoidBounds"s;
   data["payload"] = payload;
 
   return data;
diff --git a/Core/src/Surfaces/TriangleBounds.cpp b/Core/src/Surfaces/TriangleBounds.cpp
index 48e9fe4d9..e8c66df31 100644
--- a/Core/src/Surfaces/TriangleBounds.cpp
+++ b/Core/src/Surfaces/TriangleBounds.cpp
@@ -32,26 +32,26 @@ Acts::TriangleBounds::TriangleBounds(const variant_data& data_)
   : m_vertices(), m_boundingBox(0, 0)
 {
   throw_assert(data_.which() == 4, "Variant data must be map");
-  const variant_map &data = boost::get<variant_map>(data_);
-  std::string type = data.get<std::string>("type");
+  const variant_map& data = boost::get<variant_map>(data_);
+  std::string        type = data.get<std::string>("type");
   throw_assert(type == "TriangleBounds", "Type must be TriangleBounds");
 
-  const variant_map &payload = data.get<variant_map>("payload");
-  const variant_vector &vertices = payload.get<variant_vector>("vertices");
-  throw_assert(vertices.size() == 3, "Vertices for triangle must be exactly 3.");
+  const variant_map&    payload  = data.get<variant_map>("payload");
+  const variant_vector& vertices = payload.get<variant_vector>("vertices");
+  throw_assert(vertices.size() == 3,
+               "Vertices for triangle must be exactly 3.");
 
   double mx = 0, my = 0;
-  for(size_t i=0;i<3;i++) {
-    Vector2D vtx = from_variant<Vector2D>(vertices.at(i));
-    mx = std::max(mx, std::abs(vtx.x()));
-    my = std::max(my, std::abs(vtx.y()));
+  for (size_t i = 0; i < 3; i++) {
+    Vector2D vtx     = from_variant<Vector2D>(vertices.at(i));
+    mx               = std::max(mx, std::abs(vtx.x()));
+    my               = std::max(my, std::abs(vtx.y()));
     m_vertices.at(i) = vtx;
   }
 
   m_boundingBox = RectangleBounds(mx, my);
 }
 
-
 Acts::TriangleBounds::~TriangleBounds()
 {
 }
@@ -121,20 +121,20 @@ Acts::TriangleBounds::dump(std::ostream& sl) const
 }
 
 Acts::variant_data
-Acts::TriangleBounds::toVariantData() const {
+Acts::TriangleBounds::toVariantData() const
+{
   using namespace std::string_literals;
 
   variant_map payload;
-  
+
   variant_vector vertices;
-  for(const auto &vtx : m_vertices) {
+  for (const auto& vtx : m_vertices) {
     vertices.push_back(to_variant(vtx));
   }
 
   variant_map data;
-  data["type"] = "TriangleBounds"s;
+  data["type"]    = "TriangleBounds"s;
   data["payload"] = variant_map({{"vertices", vertices}});
 
   return data;
-
 }
diff --git a/Core/src/Tools/SurfaceArrayCreator.cpp b/Core/src/Tools/SurfaceArrayCreator.cpp
index 273fd0286..1dba01c30 100644
--- a/Core/src/Tools/SurfaceArrayCreator.cpp
+++ b/Core/src/Tools/SurfaceArrayCreator.cpp
@@ -36,8 +36,7 @@ Acts::SurfaceArrayCreator::surfaceArrayOnCylinder(
   ACTS_VERBOSE("Creating a SurfaceArray on a cylinder");
   ACTS_VERBOSE(" -- with " << surfaces.size() << " surfaces.")
   ACTS_VERBOSE(" -- with phi x z  = " << binsPhi << " x " << binsZ << " = "
-                                      << binsPhi * binsZ
-                                      << " bins.");
+                                      << binsPhi * binsZ << " bins.");
 
   Transform3D transform
       = _transform != nullptr ? *_transform : Transform3D::Identity();
@@ -68,7 +67,8 @@ Acts::SurfaceArrayCreator::surfaceArrayOnCylinder(
   sl->fill(surfaces);
   completeBinning(*sl, surfaces);
 
-  return std::make_unique<SurfaceArray>(std::move(sl), surfaces);
+  return std::make_unique<SurfaceArray>(
+      std::move(sl), surfaces, std::make_shared<const Transform3D>(transform));
 }
 
 std::unique_ptr<Acts::SurfaceArray>
@@ -100,8 +100,8 @@ Acts::SurfaceArrayCreator::surfaceArrayOnCylinder(
   else
     pAxisZ = createVariableAxis(surfaces, binZ, protoLayer, transform);
 
-  Transform3D itransform = transform.inverse();
-  auto globalToLocal     = [transform](const Vector3D& pos) {
+  Transform3D itransform    = transform.inverse();
+  auto        globalToLocal = [transform](const Vector3D& pos) {
     Vector3D loc = transform * pos;
     return Vector2D(loc.phi(), loc.z());
   };
@@ -126,10 +126,10 @@ Acts::SurfaceArrayCreator::surfaceArrayOnCylinder(
   ACTS_VERBOSE("Creating a SurfaceArray on a cylinder");
   ACTS_VERBOSE(" -- with " << surfaces.size() << " surfaces.")
   ACTS_VERBOSE(" -- with phi x z  = " << bins0 << " x " << bins1 << " = "
-                                      << bins0 * bins1
-                                      << " bins.");
+                                      << bins0 * bins1 << " bins.");
 
-  return std::make_unique<SurfaceArray>(std::move(sl), surfaces);
+  return std::make_unique<SurfaceArray>(
+      std::move(sl), surfaces, std::make_shared<const Transform3D>(transform));
 }
 
 std::unique_ptr<Acts::SurfaceArray>
@@ -179,12 +179,12 @@ Acts::SurfaceArrayCreator::surfaceArrayOnDisc(
 
   ACTS_VERBOSE(" -- with " << surfaces.size() << " surfaces.")
   ACTS_VERBOSE(" -- with r x phi  = " << bins0 << " x " << bins1 << " = "
-                                      << bins0 * bins1
-                                      << " bins.");
+                                      << bins0 * bins1 << " bins.");
   sl->fill(surfaces);
   completeBinning(*sl, surfaces);
 
-  return std::make_unique<SurfaceArray>(std::move(sl), surfaces);
+  return std::make_unique<SurfaceArray>(
+      std::move(sl), surfaces, std::make_shared<const Transform3D>(transform));
 }
 
 std::unique_ptr<Acts::SurfaceArray>
@@ -263,13 +263,13 @@ Acts::SurfaceArrayCreator::surfaceArrayOnDisc(
 
   ACTS_VERBOSE(" -- with " << surfaces.size() << " surfaces.")
   ACTS_VERBOSE(" -- with r x phi  = " << bins0 << " x " << bins1 << " = "
-                                      << bins0 * bins1
-                                      << " bins.");
+                                      << bins0 * bins1 << " bins.");
 
   sl->fill(surfaces);
   completeBinning(*sl, surfaces);
 
-  return std::make_unique<SurfaceArray>(std::move(sl), surfaces);
+  return std::make_unique<SurfaceArray>(
+      std::move(sl), surfaces, std::make_shared<const Transform3D>(transform));
 }
 
 /// SurfaceArrayCreator interface method - create an array on a plane
@@ -288,7 +288,7 @@ Acts::SurfaceArrayCreator::surfaceArrayOnPlane(
 
 std::vector<const Acts::Surface*>
 Acts::SurfaceArrayCreator::findKeySurfaces(
-    const std::vector<const Surface*>& surfaces,
+    const std::vector<const Surface*>&                  surfaces,
     std::function<bool(const Surface*, const Surface*)> equal) const
 {
   std::vector<const Surface*> keys;
@@ -353,8 +353,9 @@ Acts::SurfaceArrayCreator::createVariableAxis(
                                < b->binningPosition(binPhi).phi());
                      });
 
-    double maxPhi = 0.5 * (keys.at(0)->binningPosition(binPhi).phi()
-                           + keys.at(1)->binningPosition(binPhi).phi());
+    double maxPhi = 0.5
+        * (keys.at(0)->binningPosition(binPhi).phi()
+           + keys.at(1)->binningPosition(binPhi).phi());
 
     // create rotation, so that maxPhi is +pi
     double angle = -(M_PI + maxPhi);
diff --git a/Tests/Layers/CylinderLayerTests.cpp b/Tests/Layers/CylinderLayerTests.cpp
index 8e2272441..bd2c1db2f 100644
--- a/Tests/Layers/CylinderLayerTests.cpp
+++ b/Tests/Layers/CylinderLayerTests.cpp
@@ -24,6 +24,7 @@
 #include "ACTS/EventData/SingleTrackParameters.hpp"
 #include "ACTS/Layers/GenericApproachDescriptor.hpp"
 #include "ACTS/Tools/SurfaceArrayCreator.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 #include "ACTS/Volumes/CuboidVolumeBounds.hpp"
 #include "LayerStub.hpp"
 
@@ -51,10 +52,10 @@ namespace Test {
       // next level: need an array of Surfaces;
       const std::vector<const Surface*> aSurfaces{new SurfaceStub(),
                                                   new SurfaceStub()};
-      const double                      thickness(1.0);
-      SurfaceArrayCreator               sac;
-      size_t                            binsX(2), binsY(4);
-      auto                              pSurfaceArray
+      const double        thickness(1.0);
+      SurfaceArrayCreator sac;
+      size_t              binsX(2), binsY(4);
+      auto                pSurfaceArray
           = sac.surfaceArrayOnPlane(aSurfaces, 10, 20, binsX, binsY);
       auto pCylinderLayerFromSurfaces = CylinderLayer::create(
           pTransform, pCylinder, std::move(pSurfaceArray));
@@ -138,7 +139,6 @@ namespace Test {
       BOOST_TEST(cvBoundsExp->outerRadius() == cvBoundsAct->outerRadius());
       BOOST_TEST(cvBoundsExp->halfPhiSector() == cvBoundsAct->halfPhiSector());
       BOOST_TEST(cvBoundsExp->halflengthZ() == cvBoundsAct->halflengthZ());
-
     }
 
     BOOST_AUTO_TEST_SUITE_END()
diff --git a/Tests/Layers/DiscLayerTests.cpp b/Tests/Layers/DiscLayerTests.cpp
index c455b2ffc..51ada50c2 100644
--- a/Tests/Layers/DiscLayerTests.cpp
+++ b/Tests/Layers/DiscLayerTests.cpp
@@ -24,7 +24,9 @@
 #include "ACTS/EventData/SingleTrackParameters.hpp"
 #include "ACTS/Layers/GenericApproachDescriptor.hpp"
 #include "ACTS/Tools/SurfaceArrayCreator.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 #include "ACTS/Volumes/CuboidVolumeBounds.hpp"
+#include "ACTS/Volumes/CylinderVolumeBounds.hpp"
 #include "LayerStub.hpp"
 
 using boost::test_tools::output_test_stream;
@@ -98,6 +100,41 @@ namespace Test {
                  == std::string("Acts::DiscSurface"));
     }
 
+    BOOST_AUTO_TEST_CASE(DiscLayer_toVariantData)
+    {
+      Translation3D translation{0., 1., 2.};
+      auto pTransform = std::make_shared<const Transform3D>(translation);
+      const double minRad(10.), maxRad(5.);  // 20 x 10 disc
+      auto         pDisc = std::make_shared<const RadialBounds>(minRad, maxRad);
+      auto         pDiscLayer = std::dynamic_pointer_cast<DiscLayer>(
+          DiscLayer::create(pTransform, pDisc, nullptr, 6));
+
+      variant_data var_data = pDiscLayer->toVariantData();
+      std::cout << var_data << std::endl;
+
+      variant_map var_map = boost::get<variant_map>(var_data);
+      variant_map pl      = var_map.get<variant_map>("payload");
+      BOOST_TEST(pl.get<double>("thickness") == 6);
+      Transform3D act = from_variant<Transform3D>(pl.at("transform"));
+      BOOST_TEST(pTransform->isApprox(act));
+
+      auto pDiscLayer2
+          = std::dynamic_pointer_cast<DiscLayer>(DiscLayer::create(var_data));
+
+      BOOST_TEST(pDiscLayer->thickness() == pDiscLayer2->thickness());
+      BOOST_TEST(pDiscLayer->transform().isApprox(pDiscLayer2->transform()));
+
+      auto cvBoundsExp = dynamic_cast<const CylinderVolumeBounds*>(
+          &(pDiscLayer->representingVolume()->volumeBounds()));
+      auto cvBoundsAct = dynamic_cast<const CylinderVolumeBounds*>(
+          &(pDiscLayer2->representingVolume()->volumeBounds()));
+
+      BOOST_TEST(cvBoundsExp->innerRadius() == cvBoundsAct->innerRadius());
+      BOOST_TEST(cvBoundsExp->outerRadius() == cvBoundsAct->outerRadius());
+      BOOST_TEST(cvBoundsExp->halfPhiSector() == cvBoundsAct->halfPhiSector());
+      BOOST_TEST(cvBoundsExp->halflengthZ() == cvBoundsAct->halflengthZ());
+    }
+
     BOOST_AUTO_TEST_SUITE_END()
   }  // end of namespace Layers
 }  // end of namespace Test
diff --git a/Tests/Surfaces/CylinderBoundsTests.cpp b/Tests/Surfaces/CylinderBoundsTests.cpp
index 874e26da5..dbace02f2 100644
--- a/Tests/Surfaces/CylinderBoundsTests.cpp
+++ b/Tests/Surfaces/CylinderBoundsTests.cpp
@@ -134,9 +134,9 @@ namespace Test {
     BOOST_TEST(assignedCylinderBounds == cylinderBoundsObject);
   }
 
-
-  BOOST_AUTO_TEST_CASE(RectangleBounds_toVariantData) {
-    double radius = 10, avgPhi = 0, halfPhi = M_PI, halfZ = 15;
+  BOOST_AUTO_TEST_CASE(RectangleBounds_toVariantData)
+  {
+    double         radius = 10, avgPhi = 0, halfPhi = M_PI, halfZ = 15;
     CylinderBounds cylBounds(radius, avgPhi, halfPhi, halfZ);
 
     variant_data var_data = cylBounds.toVariantData();
@@ -158,8 +158,6 @@ namespace Test {
     BOOST_TEST(cylBounds.halflengthZ() == cylBounds2.halflengthZ());
   }
 
-
-
   BOOST_AUTO_TEST_SUITE_END()
 
 }  // end of namespace Test
diff --git a/Tests/Surfaces/DiamondBoundsTests.cpp b/Tests/Surfaces/DiamondBoundsTests.cpp
index 45a0157be..39769a395 100644
--- a/Tests/Surfaces/DiamondBoundsTests.cpp
+++ b/Tests/Surfaces/DiamondBoundsTests.cpp
@@ -146,12 +146,12 @@ namespace Test {
     BOOST_TEST(assignedDiamondBoundsObject == diamondBoundsObject);
   }
 
-
-  BOOST_AUTO_TEST_CASE(DiamondBounds_toVariantData) {
+  BOOST_AUTO_TEST_CASE(DiamondBounds_toVariantData)
+  {
     double minHalfX(10.), midHalfX(50.), maxHalfX(30.), halfY1(10.),
         halfY2(20.);
     DiamondBounds diam(minHalfX, midHalfX, maxHalfX, halfY1, halfY2);
-    variant_data var_data = diam.toVariantData();
+    variant_data  var_data = diam.toVariantData();
 
     std::cout << var_data << std::endl;
 
@@ -172,7 +172,6 @@ namespace Test {
     BOOST_TEST(diam.halflengthY2() == diam2.halflengthY2());
   }
 
-
   BOOST_AUTO_TEST_SUITE_END()
 
 }  // end of namespace Test
diff --git a/Tests/Surfaces/EllipseBoundsTests.cpp b/Tests/Surfaces/EllipseBoundsTests.cpp
index fb9aa0eea..e7b82272b 100644
--- a/Tests/Surfaces/EllipseBoundsTests.cpp
+++ b/Tests/Surfaces/EllipseBoundsTests.cpp
@@ -139,10 +139,12 @@ namespace Test {
     BOOST_TEST(assignedEllipseBoundsObject == ellipseBoundsObject);
   }
 
-  BOOST_AUTO_TEST_CASE(EllipseBounds_toVariantData) {
+  BOOST_AUTO_TEST_CASE(EllipseBounds_toVariantData)
+  {
     double minRad1(10.), minRad2(15.), maxRad1(15.), maxRad2(20.),
         averagePhi(0.), phiSector(M_PI / 2.);
-    EllipseBounds ell(minRad1, minRad2, maxRad1, maxRad2, averagePhi, phiSector);
+    EllipseBounds ell(
+        minRad1, minRad2, maxRad1, maxRad2, averagePhi, phiSector);
     variant_data var_data = ell.toVariantData();
 
     std::cout << var_data << std::endl;
@@ -166,7 +168,6 @@ namespace Test {
     BOOST_TEST(ell.halfPhiSector() == ell2.halfPhiSector());
   }
 
-
   BOOST_AUTO_TEST_SUITE_END()
 
 }  // end of namespace Test
diff --git a/Tests/Surfaces/LineSurfaceStub.hpp b/Tests/Surfaces/LineSurfaceStub.hpp
index e84b92a7c..ed2bf5383 100644
--- a/Tests/Surfaces/LineSurfaceStub.hpp
+++ b/Tests/Surfaces/LineSurfaceStub.hpp
@@ -11,6 +11,7 @@
 //
 #include "ACTS/Surfaces/LineSurface.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
+#include "ACTS/Utilities/VariantDataFwd.hpp"
 //
 //
 #include <limits>
@@ -40,6 +41,12 @@ public:
     : GeometryObject(), LineSurface(lbounds, detelement, identifier)
   { /* nop */
   }
+
+  LineSurfaceStub(const variant_data& data)
+    : GeometryObject(), LineSurface(data)
+  { /* nop */
+  }
+
   //
   LineSurfaceStub(const LineSurfaceStub& ls) : GeometryObject(), LineSurface(ls)
   { /* nop */
diff --git a/Tests/Surfaces/LineSurfaceTests.cpp b/Tests/Surfaces/LineSurfaceTests.cpp
index 1bf639ba1..c4e675d8a 100644
--- a/Tests/Surfaces/LineSurfaceTests.cpp
+++ b/Tests/Surfaces/LineSurfaceTests.cpp
@@ -21,6 +21,7 @@
 #include "ACTS/Material/HomogeneousSurfaceMaterial.hpp"
 #include "ACTS/Surfaces/LineSurface.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 //
 #include "DetectorElementStub.hpp"
 #include "LineSurfaceStub.hpp"
@@ -174,6 +175,30 @@ namespace Test {
                "LineSurfaces are equal value after assignment");  // operator ==
                                                                   // from base
   }
+
+  BOOST_AUTO_TEST_CASE(LineSurface_toVariantData)
+  {
+    double        radius = 2.0, hlZ = 20;
+    Translation3D translation{0., 1., 2.};
+    Transform3D   transform(translation);
+    auto          pTransform = std::make_shared<const Transform3D>(translation);
+    LineSurfaceStub line(pTransform, radius, hlZ);
+    variant_data    var_line = line.toVariantData();
+    std::cout << var_line << std::endl;
+
+    const variant_map& pl
+        = boost::get<variant_map>(var_line).get<variant_map>("payload");
+    const variant_map& bounds_pl
+        = pl.get<variant_map>("bounds").get<variant_map>("payload");
+    BOOST_TEST(bounds_pl.get<double>("radius") == radius);
+    BOOST_TEST(bounds_pl.get<double>("halfZ") == hlZ);
+
+    LineSurfaceStub line2(var_line);
+    auto            lbounds = dynamic_cast<const LineBounds*>(&line2.bounds());
+    BOOST_TEST(lbounds->r() == radius);
+    BOOST_TEST(lbounds->halflengthZ() == hlZ);
+  }
+
   BOOST_AUTO_TEST_SUITE_END()
 
 }  // end of namespace Test
diff --git a/Tests/Surfaces/PerigeeSurfaceTests.cpp b/Tests/Surfaces/PerigeeSurfaceTests.cpp
index 6897a8dfa..9bbfd53ed 100644
--- a/Tests/Surfaces/PerigeeSurfaceTests.cpp
+++ b/Tests/Surfaces/PerigeeSurfaceTests.cpp
@@ -25,6 +25,7 @@
 #include "ACTS/Surfaces/PerigeeSurface.hpp"
 #include "ACTS/Surfaces/RectangleBounds.hpp"  //to get s_noBounds
 #include "ACTS/Utilities/Definitions.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 
 using boost::test_tools::output_test_stream;
 namespace utf    = boost::unit_test;
@@ -100,6 +101,24 @@ namespace Test {
     /// Test equality of assigned to original
     BOOST_TEST(assignedPerigeeSurface == perigeeSurfaceObject);
   }
+
+  /// Unit test for testing PerigeeSurface properties
+  BOOST_AUTO_TEST_CASE(PerigeeSurface_toVariantData)
+  {
+    Vector3D       unitXYZ{1., 1., 1.};
+    PerigeeSurface perigee(unitXYZ);
+
+    variant_data var_data = perigee.toVariantData();
+    std::cout << var_data << std::endl;
+
+    // const variant_map &var_pl =
+    // boost::get<variant_map>(var_data).get<variant_map>("payload");
+    BOOST_TEST(boost::get<variant_map>(var_data).get<std::string>("type")
+               == "PerigeeSurface");
+
+    PerigeeSurface perigee2(var_data);
+    BOOST_TEST(perigee.transform().isApprox(perigee2.transform()));
+  }
   BOOST_AUTO_TEST_SUITE_END()
 
 }  // end of namespace Test
diff --git a/Tests/Surfaces/PlaneSurfaceTests.cpp b/Tests/Surfaces/PlaneSurfaceTests.cpp
index b31ce94f7..90cdcb04a 100644
--- a/Tests/Surfaces/PlaneSurfaceTests.cpp
+++ b/Tests/Surfaces/PlaneSurfaceTests.cpp
@@ -22,8 +22,8 @@
 #include "ACTS/Surfaces/RectangleBounds.hpp"
 #include "ACTS/Surfaces/TriangleBounds.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
-#include "DetectorElementStub.hpp"
 #include "ACTS/Utilities/VariantData.hpp"
+#include "DetectorElementStub.hpp"
 
 namespace tt = boost::test_tools;
 using boost::test_tools::output_test_stream;
@@ -201,20 +201,24 @@ namespace Test {
     // build
     auto rectBounds = std::make_shared<const RectangleBounds>(5, 10);
     auto idTrf = std::make_shared<const Transform3D>(Transform3D::Identity());
-    auto rot = std::make_shared<const Transform3D>(AngleAxis3D(M_PI/4., Vector3D::UnitZ()));
+    auto rot   = std::make_shared<const Transform3D>(
+        AngleAxis3D(M_PI / 4., Vector3D::UnitZ()));
     PlaneSurface rectSrf(rot, rectBounds);
     variant_data rectVariant = rectSrf.toVariantData();
     std::cout << rectVariant << std::endl;
 
     // rebuild from variant
     PlaneSurface rectSrfRec(rectVariant);
-    auto rectBoundsRec = dynamic_cast<const RectangleBounds*>(&rectSrfRec.bounds());
-    BOOST_CHECK_CLOSE(rectBounds->halflengthX(), rectBoundsRec->halflengthX(), 1e-4);
-    BOOST_CHECK_CLOSE(rectBounds->halflengthY(), rectBoundsRec->halflengthY(), 1e-4);
+    auto         rectBoundsRec
+        = dynamic_cast<const RectangleBounds*>(&rectSrfRec.bounds());
+    BOOST_CHECK_CLOSE(
+        rectBounds->halflengthX(), rectBoundsRec->halflengthX(), 1e-4);
+    BOOST_CHECK_CLOSE(
+        rectBounds->halflengthY(), rectBoundsRec->halflengthY(), 1e-4);
     BOOST_TEST(rot->isApprox(rectSrfRec.transform(), 1e-4));
 
-
-    std::array<Vector2D, 3> vertices = {{Vector2D(1,1), Vector2D(1, -1), Vector2D(-1, 1)}};
+    std::array<Vector2D, 3> vertices
+        = {{Vector2D(1, 1), Vector2D(1, -1), Vector2D(-1, 1)}};
     auto triangleBounds = std::make_shared<const TriangleBounds>(vertices);
     PlaneSurface triangleSrf(rot, triangleBounds);
     variant_data triangleVariant = triangleSrf.toVariantData();
@@ -222,17 +226,16 @@ namespace Test {
 
     // rebuild
     PlaneSurface triangleSrfRec(triangleVariant);
-    auto triangleBoundsRec = dynamic_cast<const TriangleBounds*>(&triangleSrfRec.bounds());
-    for(size_t i=0;i<3;i++) {
+    auto         triangleBoundsRec
+        = dynamic_cast<const TriangleBounds*>(&triangleSrfRec.bounds());
+    for (size_t i = 0; i < 3; i++) {
       Vector2D exp = triangleBounds->vertices().at(i);
       Vector2D act = triangleBoundsRec->vertices().at(i);
       BOOST_CHECK_CLOSE(exp.x(), act.x(), 1e-4);
       BOOST_CHECK_CLOSE(exp.y(), act.y(), 1e-4);
     }
     BOOST_TEST(rot->isApprox(triangleSrfRec.transform(), 1e-4));
-
   }
-  
 
   BOOST_AUTO_TEST_SUITE_END()
 
diff --git a/Tests/Surfaces/RectangleBoundsTests.cpp b/Tests/Surfaces/RectangleBoundsTests.cpp
index 0a609587b..e3298b61d 100644
--- a/Tests/Surfaces/RectangleBoundsTests.cpp
+++ b/Tests/Surfaces/RectangleBoundsTests.cpp
@@ -123,9 +123,10 @@ namespace Test {
     delete rectB;
   }
 
-  BOOST_AUTO_TEST_CASE(RectangleBounds_toVariantData) {
+  BOOST_AUTO_TEST_CASE(RectangleBounds_toVariantData)
+  {
     RectangleBounds rect(10, 15);
-    variant_data var_data = rect.toVariantData();
+    variant_data    var_data = rect.toVariantData();
 
     std::cout << var_data << std::endl;
     variant_map var_map = boost::get<variant_map>(var_data);
diff --git a/Tests/Surfaces/SurfaceArrayTests.cpp b/Tests/Surfaces/SurfaceArrayTests.cpp
index d6e702e4f..076105831 100644
--- a/Tests/Surfaces/SurfaceArrayTests.cpp
+++ b/Tests/Surfaces/SurfaceArrayTests.cpp
@@ -21,6 +21,7 @@
 #include "ACTS/Tools/SurfaceArrayCreator.hpp"
 #include "ACTS/Utilities/BinningType.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 #include "ACTS/Utilities/detail/Grid.hpp"
 
 namespace bdata = boost::unit_test::data;
@@ -106,6 +107,35 @@ namespace Test {
       return res;
     }
 
+    SrfVec
+    straightLineSurfaces(size_t      n        = 10.,
+                         double      step     = 3,
+                         Vector3D    origin   = {0, 0, 1.5},
+                         Transform3D pretrans = Transform3D::Identity(),
+                         Vector3D    dir      = {0, 0, 1})
+    {
+      SrfVec res;
+      for (size_t i = 0; i < n; ++i) {
+        Transform3D trans;
+        trans.setIdentity();
+        trans.translate(origin + dir * step * i);
+        // trans.rotate(AngleAxis3D(M_PI/9., Vector3D(0, 0, 1)));
+        trans.rotate(AngleAxis3D(M_PI / 2., Vector3D(1, 0, 0)));
+        trans = trans * pretrans;
+
+        auto bounds = std::make_shared<const RectangleBounds>(2, 1.5);
+
+        auto transptr = std::make_shared<const Transform3D>(trans);
+        auto srf      = std::make_unique<const PlaneSurface>(transptr, bounds);
+
+        res.push_back(srf.get());  // use raw pointer
+        m_surfaces.push_back(
+            std::move(srf));  // keep unique, will get destroyed at the end
+      }
+
+      return res;
+    }
+
     SrfVec
     makeBarrel(int nPhi, int nZ, double w, double h)
     {
@@ -240,7 +270,133 @@ namespace Test {
     BOOST_TEST(sa.surfaces().at(0) == srf.get());
   }
 
+  BOOST_FIXTURE_TEST_CASE(SurfaceArray_toVariantData, SurfaceArrayFixture)
+  {
+    SrfVec brl = makeBarrel(30, 7, 2, 1);
+
+    detail::Axis<detail::AxisType::Equidistant,
+                 detail::AxisBoundaryType::Closed>
+                        phiAxis(-M_PI, M_PI, 30u);
+    std::vector<double> zAxis_bin_edges_exp = {-14, -10, 3, 5, 8, 14};
+    detail::Axis<detail::AxisType::Variable, detail::AxisBoundaryType::Bound>
+        zAxis(zAxis_bin_edges_exp);
+
+    double angleShift = 2 * M_PI / 30. / 2.;
+    auto transform    = [angleShift](const Vector3D& pos) {
+      return Vector2D(pos.phi() + angleShift, pos.z());
+    };
+    double R        = 10;
+    auto itransform = [angleShift, R](const Vector2D& loc) {
+      return Vector3D(R * std::cos(loc[0] - angleShift),
+                      R * std::sin(loc[0] - angleShift),
+                      loc[1]);
+    };
+    SurfaceArray::SurfaceGridLookup<decltype(phiAxis), decltype(zAxis)> sl(
+        transform,
+        itransform,
+        std::make_tuple(std::move(phiAxis), std::move(zAxis)));
+    sl.fill(brl);
+    SurfaceArray sa(sl, brl);
+    sa.dump(std::cout);
+
+    variant_data data = sa.toVariantData();
+    // std::cout << data << std::endl;
+
+    const variant_map& var_map = boost::get<variant_map>(data);
+    BOOST_TEST(var_map.get<std::string>("type") == "SurfaceArray");
+    const variant_map& sa_var_pl = var_map.get<variant_map>("payload");
+    BOOST_TEST(sa_var_pl.count("surfacegridlookup") == 1);
+    const variant_map& sgl_var_pl
+        = sa_var_pl.get<variant_map>("surfacegridlookup")
+              .get<variant_map>("payload");
+    BOOST_TEST(sgl_var_pl.get<int>("dimensions") == 2);
+    const variant_vector& axes = sgl_var_pl.get<variant_vector>("axes");
+    BOOST_TEST(axes.size() == 2);
+
+    const variant_map& phiAxis_pl
+        = axes.get<variant_map>(0).get<variant_map>("payload");
+    BOOST_TEST(phiAxis_pl.get<std::string>("axisboundarytype") == "closed");
+    BOOST_TEST(phiAxis_pl.get<std::string>("axistype") == "equidistant");
+    BOOST_TEST(phiAxis_pl.get<double>("min") == -M_PI);
+    BOOST_TEST(phiAxis_pl.get<double>("max") == M_PI);
+    BOOST_TEST(phiAxis_pl.get<int>("nbins") == 30);
+
+    const variant_map& zAxis_pl
+        = axes.get<variant_map>(1).get<variant_map>("payload");
+    BOOST_TEST(zAxis_pl.get<std::string>("axisboundarytype") == "bound");
+    BOOST_TEST(zAxis_pl.get<std::string>("axistype") == "variable");
+    const variant_vector& zAxis_bin_edges
+        = zAxis_pl.get<variant_vector>("bin_edges");
+    BOOST_TEST(zAxis_bin_edges.size() == 6);
+    for (size_t i = 0; i < zAxis_bin_edges.size(); i++) {
+      BOOST_TEST(zAxis_bin_edges.get<double>(i) == zAxis_bin_edges_exp.at(i));
+    }
+
+    SurfaceArray sa2(data, transform, itransform);
+    sa2.dump(std::cout);
+
+    std::ostringstream dumpExp_os;
+    sa.dump(dumpExp_os);
+    std::string                           dumpExp = dumpExp_os.str();
+    boost::test_tools::output_test_stream dumpAct;
+    sa2.dump(dumpAct);
+    BOOST_TEST(dumpAct.is_equal(dumpExp));
+  }
+
+  BOOST_FIXTURE_TEST_CASE(SurfaceArray_toVariantData_1D, SurfaceArrayFixture)
+  {
+    detail::Axis<detail::AxisType::Equidistant, detail::AxisBoundaryType::Bound>
+         zAxis(0, 30, 10);
+    auto transform = [](const Vector3D& pos) {
+      return std::array<double, 1>({{pos.z()}});
+    };
+    auto itransform = [](const std::array<double, 1>& loc) {
+      return Vector3D(0, 0, loc[0]);
+    };
+    SurfaceArray::SurfaceGridLookup<decltype(zAxis)> sl(
+        transform, itransform, std::make_tuple(zAxis));
+
+    // same thing in 1D
+    SrfVec line = straightLineSurfaces();
+    sl.fill(line);
+    SurfaceArray sa(sl, line);
+
+    sa.dump(std::cout);
+
+    variant_data data = sa.toVariantData();
+    // std::cout << data << std::endl;
+
+    const variant_map& var_map = boost::get<variant_map>(data);
+    BOOST_TEST(var_map.get<std::string>("type") == "SurfaceArray");
+    const variant_map& sa_var_pl = var_map.get<variant_map>("payload");
+    BOOST_TEST(sa_var_pl.count("surfacegridlookup") == 1);
+    const variant_map& sgl_var_pl
+        = sa_var_pl.get<variant_map>("surfacegridlookup")
+              .get<variant_map>("payload");
+    BOOST_TEST(sgl_var_pl.get<int>("dimensions") == 1);
+    const variant_vector& axes = sgl_var_pl.get<variant_vector>("axes");
+    BOOST_TEST(axes.size() == 1);
+
+    const variant_map& zAxis_pl
+        = axes.get<variant_map>(0).get<variant_map>("payload");
+    BOOST_TEST(zAxis_pl.get<std::string>("axisboundarytype") == "bound");
+    BOOST_TEST(zAxis_pl.get<std::string>("axistype") == "equidistant");
+    BOOST_TEST(zAxis_pl.get<double>("min") == 0);
+    BOOST_TEST(zAxis_pl.get<double>("max") == 30);
+    BOOST_TEST(zAxis_pl.get<int>("nbins") == 10);
+
+    SurfaceArray sa2(data, transform, itransform);
+    sa2.dump(std::cout);
+
+    std::ostringstream dumpExp_os;
+    sa.dump(dumpExp_os);
+    std::string                           dumpExp = dumpExp_os.str();
+    boost::test_tools::output_test_stream dumpAct;
+    sa2.dump(dumpAct);
+    BOOST_TEST(dumpAct.is_equal(dumpExp));
+  }
+
   BOOST_AUTO_TEST_SUITE_END()
-}  // end of namespace Test
+}  // namespace Test
 
 }  // end of namespace Acts
diff --git a/Tests/Surfaces/SurfaceStub.hpp b/Tests/Surfaces/SurfaceStub.hpp
index a5af4c685..42ca8bcc9 100644
--- a/Tests/Surfaces/SurfaceStub.hpp
+++ b/Tests/Surfaces/SurfaceStub.hpp
@@ -12,6 +12,7 @@
 #include "ACTS/Surfaces/PlanarBounds.hpp"
 #include "ACTS/Surfaces/Surface.hpp"
 #include "ACTS/Utilities/Definitions.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 
 namespace Acts {
 /// Surface derived class stub
@@ -120,6 +121,15 @@ public:
     return true;
   }
 
+  virtual variant_data
+  toVariantData() const override
+  {
+    using namespace std::string_literals;
+    variant_map data;
+    data["type"] = "SurfaceStub"s;
+    return data;
+  }
+
 private:
   /// the bounds of this surface
   std::shared_ptr<const PlanarBounds> m_bounds;
diff --git a/Tests/Surfaces/TrapezoidBoundsTests.cpp b/Tests/Surfaces/TrapezoidBoundsTests.cpp
index a3f537a06..a27b08a07 100644
--- a/Tests/Surfaces/TrapezoidBoundsTests.cpp
+++ b/Tests/Surfaces/TrapezoidBoundsTests.cpp
@@ -120,13 +120,13 @@ namespace Test {
     BOOST_TEST(assignedTrapezoidBoundsObject == trapezoidBoundsObject);
   }
 
-
-  BOOST_AUTO_TEST_CASE(TrapezoidBounds_toVariantData) {
-    double minHlX = 10;
-    double maxHlX = 15;
-    double hlY = 5;
+  BOOST_AUTO_TEST_CASE(TrapezoidBounds_toVariantData)
+  {
+    double          minHlX = 10;
+    double          maxHlX = 15;
+    double          hlY    = 5;
     TrapezoidBounds trap(minHlX, maxHlX, hlY);
-    variant_data var_data = trap.toVariantData();
+    variant_data    var_data = trap.toVariantData();
 
     std::cout << var_data << std::endl;
 
diff --git a/Tests/Surfaces/TriangleBoundsTests.cpp b/Tests/Surfaces/TriangleBoundsTests.cpp
index e27b9a45b..6d56908c7 100644
--- a/Tests/Surfaces/TriangleBoundsTests.cpp
+++ b/Tests/Surfaces/TriangleBoundsTests.cpp
@@ -23,8 +23,8 @@
 #include "ACTS/Utilities/Definitions.hpp"
 #include "ACTS/Utilities/VariantData.hpp"
 
-namespace utf    = boost::unit_test;
-//const double NaN = std::numeric_limits<double>::quiet_NaN();
+namespace utf = boost::unit_test;
+// const double NaN = std::numeric_limits<double>::quiet_NaN();
 
 namespace Acts {
 
@@ -33,7 +33,9 @@ namespace Test {
   /// Unit test for creating compliant/non-compliant TriangleBounds object
   BOOST_AUTO_TEST_CASE(TriangleBoundsConstruction)
   {
-    std::array<Vector2D, 3> vertices({{Vector2D(1., 1.), Vector2D(4., 1.), Vector2D(4., 5.)}});  // 3-4-5 triangle
+    std::array<Vector2D, 3> vertices({{Vector2D(1., 1.),
+                                       Vector2D(4., 1.),
+                                       Vector2D(4., 5.)}});  // 3-4-5 triangle
     // test default construction
     // TriangleBounds defaultConstructedTriangleBounds;  //deleted
     //
@@ -49,8 +51,9 @@ namespace Test {
   /// Unit tests for TriangleBounds properties
   BOOST_AUTO_TEST_CASE(TriangleBoundsProperties)
   {
-    std::array<Vector2D, 3> vertices({{
-        Vector2D(1., 1.), Vector2D(4., 1.), Vector2D(4., 5.)}});  // 3-4-5 triangle
+    std::array<Vector2D, 3> vertices({{Vector2D(1., 1.),
+                                       Vector2D(4., 1.),
+                                       Vector2D(4., 5.)}});  // 3-4-5 triangle
     /// Test clone
     TriangleBounds triangleBoundsObject(vertices);
     auto           pClonedTriangleBounds = triangleBoundsObject.clone();
@@ -73,7 +76,7 @@ namespace Test {
     BOOST_TEST_MESSAGE(
         "Following two tests fail because the triangle has six vertices");
     BOOST_TEST(triangleBoundsObject.vertices().size() == (size_t)3);
-    for(size_t i=0;i<3;i++) {
+    for (size_t i = 0; i < 3; i++) {
       Vector2D act = triangleBoundsObject.vertices().at(i);
       Vector2D exp = expectedVertices.at(i);
       BOOST_CHECK_CLOSE(act[0], exp[0], 1e-6);
@@ -101,24 +104,29 @@ namespace Test {
   /// Unit test for testing TriangleBounds assignment
   BOOST_AUTO_TEST_CASE(TriangleBoundsAssignment)
   {
-    std::array<Vector2D, 3> vertices({{
-        Vector2D(1., 1.), Vector2D(4., 1.), Vector2D(4., 5)}});  // 3-4-5 triangle
-    std::array<Vector2D, 3> invalid({{Vector2D(-1, -1), Vector2D(-1, -1), Vector2D(-1, -1)}});
-    TriangleBounds        triangleBoundsObject(vertices);
+    std::array<Vector2D, 3> vertices({{Vector2D(1., 1.),
+                                       Vector2D(4., 1.),
+                                       Vector2D(4., 5)}});  // 3-4-5 triangle
+    std::array<Vector2D, 3> invalid(
+        {{Vector2D(-1, -1), Vector2D(-1, -1), Vector2D(-1, -1)}});
+    TriangleBounds triangleBoundsObject(vertices);
     // operator == not implemented in this class
     //
     /// Test assignment
     TriangleBounds assignedTriangleBoundsObject(
         invalid);  // invalid object, in some sense
     assignedTriangleBoundsObject = triangleBoundsObject;
-    BOOST_TEST(assignedTriangleBoundsObject.vertices() == triangleBoundsObject.vertices());
+    BOOST_TEST(assignedTriangleBoundsObject.vertices()
+               == triangleBoundsObject.vertices());
   }
 
-  BOOST_AUTO_TEST_CASE(TriangleBounds_toVariantData) {
-    std::array<Vector2D, 3> vertices({{
-        Vector2D(1., 1.), Vector2D(4., 1.), Vector2D(4., 5.)}});  // 3-4-5 triangle
+  BOOST_AUTO_TEST_CASE(TriangleBounds_toVariantData)
+  {
+    std::array<Vector2D, 3> vertices({{Vector2D(1., 1.),
+                                       Vector2D(4., 1.),
+                                       Vector2D(4., 5.)}});  // 3-4-5 triangle
     TriangleBounds triangle(vertices);
-    variant_data var_data = triangle.toVariantData();
+    variant_data   var_data = triangle.toVariantData();
 
     std::cout << var_data << std::endl;
 
@@ -128,9 +136,9 @@ namespace Test {
 
     variant_vector var_vertices = pl.get<variant_vector>("vertices");
     BOOST_TEST(var_vertices.size() == 3);
-    
-    for(size_t i=0;i<3;i++) {
-      Vector2D exp = vertices.at(i);
+
+    for (size_t i = 0; i < 3; i++) {
+      Vector2D    exp = vertices.at(i);
       variant_map var = var_vertices.get<variant_map>(i);
       BOOST_TEST(var.get<std::string>("type") == "Vector2D");
       variant_vector coords = var.get<variant_vector>("payload");
@@ -143,7 +151,6 @@ namespace Test {
     BOOST_TEST(triangle2.vertices().size() == 3);
   }
 
-
   BOOST_AUTO_TEST_SUITE_END()
 
 }  // end of namespace Test
diff --git a/Tests/Tools/LayerCreatorTests.cpp b/Tests/Tools/LayerCreatorTests.cpp
index df79b069c..fa57b033a 100644
--- a/Tests/Tools/LayerCreatorTests.cpp
+++ b/Tests/Tools/LayerCreatorTests.cpp
@@ -6,7 +6,7 @@
 // License, v. 2.0. If a copy of the MPL was not distributed with this
 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
-#define BOOST_TEST_MODULE SurfaceArrayCreator
+#define BOOST_TEST_MODULE LayerCreator
 #include <boost/test/included/unit_test.hpp>
 
 #include <boost/format.hpp>
@@ -22,6 +22,7 @@
 #include "ACTS/Tools/LayerCreator.hpp"
 #include "ACTS/Tools/SurfaceArrayCreator.hpp"
 #include "ACTS/Utilities/BinningType.hpp"
+#include "ACTS/Utilities/VariantData.hpp"
 
 #include "ACTS/Utilities/Definitions.hpp"
 
@@ -490,6 +491,107 @@ namespace Test {
     checkBinning(*layer->surfaceArray());
   }
 
+  BOOST_FIXTURE_TEST_CASE(LayerCreator_Cylinder_toVariantData,
+                          LayerCreatorFixture)
+  {
+
+    auto barrel = makeBarrelStagger(30, 7, 0, M_PI / 9.);
+    auto brl    = barrel.first;
+    draw_surfaces(brl, "LayerCreator_barrelStagger.obj");
+
+    ProtoLayer                     pl(brl);
+    std::shared_ptr<CylinderLayer> layer
+        = std::dynamic_pointer_cast<CylinderLayer>(
+            p_LC->cylinderLayer(brl, equidistant, equidistant, pl));
+
+    std::cout << (*layer->surfaceArray()) << std::endl;
+
+    const variant_data var_layer = layer->toVariantData();
+    // std::cout << var_layer << std::endl;
+
+    auto layer2 = std::dynamic_pointer_cast<CylinderLayer>(
+        CylinderLayer::create(var_layer));
+    std::cout << (*layer2->surfaceArray()) << std::endl;
+
+    auto sa  = layer->surfaceArray();
+    auto sa2 = layer2->surfaceArray();
+
+    BOOST_TEST(sa);
+    BOOST_TEST(sa2);
+
+    BOOST_TEST(sa->transform().isApprox(sa2->transform()));
+
+    // let's make sure the binning is really ok
+    // we check that lookup at center of input surface centers
+    // gives the same number of bin content surfaces
+    // which also have compatible transforms.
+    // This is as close to "ok" as we can get I think.
+    for (const auto& pr : barrel.second) {
+      auto A = pr.first;
+
+      Vector3D ctr = A->binningPosition(binR);
+
+      std::vector<const Surface*> bc1 = sa->at(ctr);
+      std::vector<const Surface*> bc2 = sa2->at(ctr);
+
+      BOOST_TEST(bc1.size() == bc2.size());
+
+      for (size_t i = 0; i < bc1.size(); i++) {
+        auto srf1 = bc1.at(i);
+        auto srf2 = bc2.at(i);
+
+        BOOST_TEST(srf1->transform().isApprox(srf2->transform()));
+      }
+    }
+  }
+
+  BOOST_FIXTURE_TEST_CASE(LayerCreator_Disc_toVariantData, LayerCreatorFixture)
+  {
+    std::vector<const Surface*> surfaces;
+    auto                        ringa = fullPhiTestSurfacesEC(30, 0, 0, 10);
+    surfaces.insert(surfaces.end(), ringa.begin(), ringa.end());
+    auto ringb = fullPhiTestSurfacesEC(30, 0, 0, 15);
+    surfaces.insert(surfaces.end(), ringb.begin(), ringb.end());
+    auto ringc = fullPhiTestSurfacesEC(30, 0, 0, 20);
+    surfaces.insert(surfaces.end(), ringc.begin(), ringc.end());
+
+    ProtoLayer                 pl(surfaces);
+    std::shared_ptr<DiscLayer> layer = std::dynamic_pointer_cast<DiscLayer>(
+        p_LC->discLayer(surfaces, equidistant, equidistant, pl));
+
+    const variant_data var_layer = layer->toVariantData();
+
+    std::cout << (*layer->surfaceArray()) << std::endl;
+    auto layer2
+        = std::dynamic_pointer_cast<DiscLayer>(DiscLayer::create(var_layer));
+    std::cout << (*layer2->surfaceArray()) << std::endl;
+
+    auto sa  = layer->surfaceArray();
+    auto sa2 = layer2->surfaceArray();
+
+    BOOST_TEST(sa);
+    BOOST_TEST(sa2);
+
+    BOOST_TEST(sa->transform().isApprox(sa2->transform()));
+
+    for (const auto& srfRef : surfaces) {
+
+      Vector3D ctr = srfRef->binningPosition(binR);
+
+      std::vector<const Surface*> bc1 = sa->at(ctr);
+      std::vector<const Surface*> bc2 = sa2->at(ctr);
+
+      BOOST_TEST(bc1.size() == bc2.size());
+
+      for (size_t i = 0; i < bc1.size(); i++) {
+        auto srf1 = bc1.at(i);
+        auto srf2 = bc2.at(i);
+
+        BOOST_TEST(srf1->transform().isApprox(srf2->transform()));
+      }
+    }
+  }
+
   BOOST_AUTO_TEST_SUITE_END()
 }  // end of namespace Test
 
-- 
GitLab