From 5c3071befb8b77fbb408af06b8950d3cb7a82633 Mon Sep 17 00:00:00 2001
From: Joana Niermann <joana.niermann@cern.ch>
Date: Wed, 2 Aug 2023 20:41:49 +0200
Subject: [PATCH] Add the grid reader to the detector reader. It generates all
 possible grid types and then selects the ones to be added to the detector by
 switch case statements on the type ids that are in the payloads. In the grid
 builder, the missing information is added to the surface descriptors, once
 the basic volume geometry has been build.

It contains some other fixes, like correctly adapting the cylinder grid radius to the volume,
updating the index offsets when pushing a grid into a grid collection,
as well as maintaining correct surface sorting during the IO, so that the surface indices
saved in the grid match still match the content of the surface_lookup.O
---
 .../detray/core/detail/multi_store.hpp        |   7 +
 .../detray/core/detail/single_store.hpp       |   7 +-
 core/include/detray/definitions/grid_axis.hpp |   2 +-
 core/include/detray/geometry/barcode.hpp      |   2 +-
 .../geometry/detail/surface_descriptor.hpp    |  11 +
 core/include/detray/geometry/surface.hpp      |   7 +-
 core/include/detray/masks/masks.hpp           |   7 +
 .../surface_finders/brute_force_finder.hpp    |  24 ++
 .../detray/surface_finders/grid/grid.hpp      |  19 +-
 .../surface_finders/grid/grid_collection.hpp  |  43 +++
 core/include/detray/tools/bin_fillers.hpp     |   7 +-
 core/include/detray/tools/bounding_volume.hpp |   1 -
 .../include/detray/tools/detector_builder.hpp |   6 +-
 core/include/detray/tools/grid_builder.hpp    | 145 +++----
 core/include/detray/tools/grid_factory.hpp    |  57 ++-
 core/include/detray/tools/surface_factory.hpp |  47 ++-
 .../tools/surface_factory_interface.hpp       |  26 +-
 core/include/detray/tools/volume_builder.hpp  |   2 +-
 .../detray/utils/consistency_checker.hpp      |   1 +
 core/include/detray/utils/inspectors.hpp      |   1 -
 core/include/detray/utils/type_list.hpp       | 170 +++++++++
 .../common/detail/detector_components_io.hpp  |  12 +-
 .../detray/io/common/detector_reader.hpp      |  13 +-
 .../detray/io/common/geometry_reader.hpp      |  21 +-
 .../detray/io/common/geometry_writer.hpp      |  10 +-
 io/include/detray/io/common/grid_reader.hpp   | 356 ++++++++++++++++++
 io/include/detray/io/common/grid_writer.hpp   |  43 ++-
 .../io/common/homogeneous_material_reader.hpp |   5 +-
 io/include/detray/io/common/payloads.hpp      |  12 +-
 .../detray/io/json/json_geometry_io.hpp       |   6 +
 io/include/detray/io/json/json_grids_io.hpp   |  28 +-
 io/include/detray/io/json/json_reader.hpp     |  16 +-
 io/include/detray/io/json/json_writer.hpp     |  12 +-
 .../plugins/svgtools/conversion/surface.hpp   |   3 +-
 .../plugins/svgtools/utils/mask_utils.hpp     |   3 +-
 .../tests/common/test_toy_detector.hpp        |  51 +--
 .../include/tests/common/tools/hash_tree.hpp  |   1 -
 tests/unit_tests/core/CMakeLists.txt          |   1 +
 tests/unit_tests/core/tools_hash_tree.cpp     |   1 -
 tests/unit_tests/core/utils_type_list.cpp     |  40 ++
 tests/unit_tests/cpu/grid_grid_builder.cpp    |  35 +-
 .../unit_tests/cpu/tools_material_builder.cpp |   3 +-
 .../unit_tests/io/io_json_detector_reader.cpp |  12 +-
 tests/unit_tests/io/io_json_payload.cpp       |  12 +-
 .../detray/detectors/create_toy_geometry.hpp  |  84 +++--
 .../detray/detectors/create_wire_chamber.hpp  |  34 +-
 .../detray/detectors/detector_helper.hpp      |   4 +-
 .../include/detray/detectors/toy_metadata.hpp |   2 +-
 48 files changed, 1131 insertions(+), 281 deletions(-)
 create mode 100644 core/include/detray/utils/type_list.hpp
 create mode 100644 io/include/detray/io/common/grid_reader.hpp
 create mode 100644 tests/unit_tests/core/utils_type_list.cpp

diff --git a/core/include/detray/core/detail/multi_store.hpp b/core/include/detray/core/detail/multi_store.hpp
index 718df5027..83a9fcd62 100644
--- a/core/include/detray/core/detail/multi_store.hpp
+++ b/core/include/detray/core/detail/multi_store.hpp
@@ -207,6 +207,13 @@ class multi_store {
         detail::get<value_types::to_index(id)>(m_tuple_container).reserve(n);
     }
 
+    /// Resize the underlying container to @param n for a collection given by
+    /// @tparam id
+    template <ID id>
+    DETRAY_HOST void resize(std::size_t n, const context_type & /*ctx*/) {
+        detail::get<value_types::to_index(id)>(m_tuple_container).resize(n);
+    }
+
     /// Add a new element to a collection
     ///
     /// @tparam ID is the id of the collection
diff --git a/core/include/detray/core/detail/single_store.hpp b/core/include/detray/core/detail/single_store.hpp
index 9ba70059b..8c7f3893c 100644
--- a/core/include/detray/core/detail/single_store.hpp
+++ b/core/include/detray/core/detail/single_store.hpp
@@ -164,11 +164,16 @@ class single_store {
         m_container.clear();
     }
 
-    /// Reserve memory of size @param n for a collection given by @tparam id
+    /// Reserve memory of size @param n for a given geometry context
     DETRAY_HOST void reserve(std::size_t n, const context_type & /*ctx*/) {
         m_container.reserve(n);
     }
 
+    /// Resize the underlying container to @param n for a given geometry context
+    DETRAY_HOST void resize(std::size_t n, const context_type & /*ctx*/) {
+        m_container.resize(n);
+    }
+
     /// Add a new element to the collection - copy
     ///
     /// @tparam U type that can be converted to T
diff --git a/core/include/detray/definitions/grid_axis.hpp b/core/include/detray/definitions/grid_axis.hpp
index 08e9285cf..0524bc54d 100644
--- a/core/include/detray/definitions/grid_axis.hpp
+++ b/core/include/detray/definitions/grid_axis.hpp
@@ -62,4 +62,4 @@ using bounds_t =
                        n_axis::closed<axis_label>>;
 /// @}
 
-}  // namespace detray::n_axis
\ No newline at end of file
+}  // namespace detray::n_axis
diff --git a/core/include/detray/geometry/barcode.hpp b/core/include/detray/geometry/barcode.hpp
index 7c48d5d4b..e2e48941f 100644
--- a/core/include/detray/geometry/barcode.hpp
+++ b/core/include/detray/geometry/barcode.hpp
@@ -142,7 +142,7 @@ class barcode {
     DETRAY_HOST
     friend std::ostream& operator<<(std::ostream& os, const barcode c) {
         if (c.is_invalid()) {
-            return (os << "undefined");
+            os << "INVALID: ";
         }
 
         static const char* const names[] = {
diff --git a/core/include/detray/geometry/detail/surface_descriptor.hpp b/core/include/detray/geometry/detail/surface_descriptor.hpp
index 3af3cc509..9c4515d13 100644
--- a/core/include/detray/geometry/detail/surface_descriptor.hpp
+++ b/core/include/detray/geometry/detail/surface_descriptor.hpp
@@ -192,6 +192,17 @@ class surface_descriptor {
         return m_barcode.id() == surface_id::e_passive;
     }
 
+    /// @returns a string stream that prints the surface details
+    DETRAY_HOST
+    friend std::ostream &operator<<(std::ostream &os,
+                                    const surface_descriptor &sf) {
+        os << sf.m_barcode;
+        os << " | trf.: " << sf._trf;
+        os << " | mask: " << sf._mask;
+        os << " | mat.: " << sf._material;
+        return os;
+    }
+
     private:
     geometry::barcode m_barcode{};
     mask_link _mask{};
diff --git a/core/include/detray/geometry/surface.hpp b/core/include/detray/geometry/surface.hpp
index 2a2ee1927..a3f66cdd0 100644
--- a/core/include/detray/geometry/surface.hpp
+++ b/core/include/detray/geometry/surface.hpp
@@ -357,12 +357,7 @@ class surface {
     /// @returns a string stream that prints the surface details
     DETRAY_HOST
     friend std::ostream &operator<<(std::ostream &os, const surface &sf) {
-        os << sf.barcode();
-        os << " | trf.: " << sf.m_desc.transform();
-        os << " | mask: " << static_cast<dindex>(sf.m_desc.mask().id()) << ", "
-           << sf.m_desc.mask().index();
-        os << " | mat.: " << static_cast<dindex>(sf.m_desc.material().id())
-           << ", " << sf.m_desc.material().index();
+        os << sf.m_desc;
         return os;
     }
 
diff --git a/core/include/detray/masks/masks.hpp b/core/include/detray/masks/masks.hpp
index fe0ba6ae8..778cc07f0 100644
--- a/core/include/detray/masks/masks.hpp
+++ b/core/include/detray/masks/masks.hpp
@@ -244,6 +244,13 @@ class mask {
         return ss.str();
     }
 
+    /// @returns a string stream that prints the mask details
+    DETRAY_HOST
+    friend std::ostream& operator<<(std::ostream& os, const mask& m) {
+        os << m.to_string();
+        return os;
+    }
+
     private:
     shape _shape;
     mask_values _values;
diff --git a/core/include/detray/surface_finders/brute_force_finder.hpp b/core/include/detray/surface_finders/brute_force_finder.hpp
index 62efa7026..50ad35b0e 100644
--- a/core/include/detray/surface_finders/brute_force_finder.hpp
+++ b/core/include/detray/surface_finders/brute_force_finder.hpp
@@ -11,6 +11,7 @@
 #include "detray/core/detail/container_buffers.hpp"
 #include "detray/core/detail/container_views.hpp"
 #include "detray/definitions/containers.hpp"
+#include "detray/definitions/detail/algorithms.hpp"
 #include "detray/definitions/indexing.hpp"
 #include "detray/definitions/qualifiers.hpp"
 #include "detray/utils/ranges.hpp"
@@ -131,6 +132,12 @@ class brute_force_collection {
         return size() == size_type{0};
     }
 
+    /// @returns access to the volume offsets - const
+    DETRAY_HOST const auto& offsets() const { return m_offsets; }
+
+    /// @returns access to the volume offsets
+    DETRAY_HOST auto& offsets() { return m_offsets; }
+
     /// @return access to the surface container - const.
     DETRAY_HOST_DEVICE
     auto all() const -> const vector_type<surface_t>& { return m_surfaces; }
@@ -161,6 +168,23 @@ class brute_force_collection {
         m_offsets.push_back(static_cast<dindex>(m_surfaces.size()));
     }
 
+    /// Remove surface from collection
+    DETRAY_HOST auto erase(
+        typename vector_type<surface_t>::iterator pos) noexcept(false) {
+        // Remove one element
+        auto next = m_surfaces.erase(pos);
+
+        // Update the upper bound of the range and all following ranges
+        const auto idx{static_cast<std::size_t>(pos - m_surfaces.begin())};
+        auto offset =
+            detail::upper_bound(m_offsets.begin(), m_offsets.end(), idx);
+        for (auto itr = offset; itr != m_offsets.end(); ++itr) {
+            --(*itr);
+        }
+
+        return next;
+    }
+
     /// @return the view on the brute force finders - non-const
     DETRAY_HOST
     constexpr auto get_data() noexcept -> view_type {
diff --git a/core/include/detray/surface_finders/grid/grid.hpp b/core/include/detray/surface_finders/grid/grid.hpp
index 33d07a9bb..17939a75d 100644
--- a/core/include/detray/surface_finders/grid/grid.hpp
+++ b/core/include/detray/surface_finders/grid/grid.hpp
@@ -202,12 +202,18 @@ class grid {
         return bin(n_axis::multi_bin<Dim>{{indices...}});
     }
 
-    /// @param mbin the multi-index of bins over all axes
+    /// @param mbin the multi-index of bins over all axes - const
     DETRAY_HOST_DEVICE
     auto bin(const n_axis::multi_bin<Dim> &mbin) const {
         return bin(m_serializer(m_axes, mbin));
     }
 
+    /// @param mbin the multi-index of bins over all axes
+    DETRAY_HOST_DEVICE
+    auto bin(const n_axis::multi_bin<Dim> &mbin) {
+        return bin(m_serializer(m_axes, mbin));
+    }
+
     /// @param gbin the global bin index - const
     DETRAY_HOST_DEVICE
     auto bin(const dindex gbin) const {
@@ -276,7 +282,7 @@ class grid {
         return all();
     }
 
-    /// Find the value of a single bin
+    /// Find the value of a single bin - const
     ///
     /// @param p is point in the local frame
     ///
@@ -286,6 +292,15 @@ class grid {
         return bin(m_axes.bins(p));
     }
 
+    /// Find the value of a single bin
+    ///
+    /// @param p is point in the local frame
+    ///
+    /// @return the iterable view of the bin content
+    DETRAY_HOST_DEVICE auto search(const typename local_frame_type::point3 &p) {
+        return bin(m_axes.bins(p));
+    }
+
     /// @brief Return a neighborhood of values from the grid
     ///
     /// The lookup is done with a neighborhood around the bin which contains the
diff --git a/core/include/detray/surface_finders/grid/grid_collection.hpp b/core/include/detray/surface_finders/grid/grid_collection.hpp
index 29ea44597..9e86f901f 100644
--- a/core/include/detray/surface_finders/grid/grid_collection.hpp
+++ b/core/include/detray/surface_finders/grid/grid_collection.hpp
@@ -181,21 +181,64 @@ class grid_collection<
     DETRAY_HOST constexpr auto push_back(
         const typename grid_type::template type<true> &gr) noexcept(false)
         -> void {
+        // Current offset into the global bin storage for the new grid
         m_offsets.push_back(static_cast<size_type>(m_bins.size()));
 
+        // Add the bins of the new grid to the collection
         const auto *grid_bins = gr.data().bin_data();
         m_bins.insert(m_bins.end(), grid_bins->begin(), grid_bins->end());
 
+        // Add the axes data of the new grid to the collection
+        // (how to lookup the axis bin edges)
         const auto *axes_data = gr.axes().data().axes_data();
         m_axes_data.insert(m_axes_data.end(), axes_data->begin(),
                            axes_data->end());
 
+        // Current offset into the global bin edges storage
+        dindex bin_edges_offset = static_cast<dindex>(m_bin_edges.size());
+        // The binning types of the axes
+        const auto binnings = get_binning(gr.axes());
+
+        // Update the bin edges index range for the axes in the grid collection
+        for (std::size_t i = m_axes_data.size() - grid_type::Dim;
+             i < m_axes_data.size(); ++i) {
+            auto &bin_entry_range = m_axes_data[i];
+            bin_entry_range[0] += bin_edges_offset;
+            // If the axis has irregular binning, the second entry is the index
+            // of the last bin edge value instead of the number of bins
+            if (binnings[i] == n_axis::binning::e_irregular) {
+                bin_entry_range[1] += bin_edges_offset;
+            }
+        }
+
+        // Add the bin edges of the new grid to the collection
         const auto *bin_edges = gr.axes().data().edges();
+
         m_bin_edges.insert(m_bin_edges.end(), bin_edges->begin(),
                            bin_edges->end());
     }
 
     private:
+    /// @returns an array that contians the binning enum values of the
+    /// corresponding single axes in @param axes
+    template <bool ownership, typename local_frame_t, typename... axis_ts>
+    static auto get_binning(
+        const n_axis::multi_axis<ownership, local_frame_t, axis_ts...> &axes) {
+
+        // Serialize every single axis and construct array from their payloads
+        std::array<n_axis::binning, sizeof...(axis_ts)> binning_types{
+            get_binning(axes.template get_axis<axis_ts>())...};
+
+        return binning_types;
+    }
+
+    /// @returns the binning enum value of a single axis @param axis
+    template <typename bounds_t, typename binning_t>
+    static auto get_binning(
+        const n_axis::single_axis<bounds_t, binning_t> &axis) {
+        return axis.binning();
+    }
+
     /// Offsets for the respective grids into the bin storage
     vector_type<size_type> m_offsets{};
     /// Contains the bin content for all grids
diff --git a/core/include/detray/tools/bin_fillers.hpp b/core/include/detray/tools/bin_fillers.hpp
index f1f8d7c98..7d9c7c742 100644
--- a/core/include/detray/tools/bin_fillers.hpp
+++ b/core/include/detray/tools/bin_fillers.hpp
@@ -69,10 +69,13 @@ struct fill_by_pos {
         // Fill the volumes surfaces into the grid
         for (const auto [idx, sf] : detray::views::enumerate(surfaces)) {
             // no portals in grids allowed
-            assert(not sf.is_portal());
+            if (not sf.is_sensitive() or sf.volume() != vol.index()) {
+                continue;
+            }
 
-            const auto &sf_trf = transforms.at(idx, ctx);
+            const auto &sf_trf = transforms.at(sf.transform(), ctx);
             const auto &t = sf_trf.translation();
+
             // transform to axis coordinate system
             const auto loc_pos = grid.global_to_local(vol.transform(), t, t);
 
diff --git a/core/include/detray/tools/bounding_volume.hpp b/core/include/detray/tools/bounding_volume.hpp
index 5b0661af2..1c973aaa7 100644
--- a/core/include/detray/tools/bounding_volume.hpp
+++ b/core/include/detray/tools/bounding_volume.hpp
@@ -15,7 +15,6 @@
 
 // System include(s)
 #include <cassert>
-#include <iostream>
 #include <limits>
 #include <type_traits>
 #include <vector>
diff --git a/core/include/detray/tools/detector_builder.hpp b/core/include/detray/tools/detector_builder.hpp
index 6528c25c1..f07f1675f 100644
--- a/core/include/detray/tools/detector_builder.hpp
+++ b/core/include/detray/tools/detector_builder.hpp
@@ -59,12 +59,12 @@ class detector_builder {
 
     /// Decorate a volume builder at position @param volume_idx with more
     /// functionality
-    template <template <typename> class builder_t>
+    template <class builder_t>
     DETRAY_HOST auto decorate(dindex volume_idx)
         -> volume_builder_interface<detector_type>* {
 
-        m_volumes[volume_idx] = std::make_unique<builder_t<detector_type>>(
-            std::move(m_volumes[volume_idx]));
+        m_volumes[volume_idx] =
+            std::make_unique<builder_t>(std::move(m_volumes[volume_idx]));
 
         return m_volumes[volume_idx].get();
     }
diff --git a/core/include/detray/tools/grid_builder.hpp b/core/include/detray/tools/grid_builder.hpp
index 4f6b617e4..b97958124 100644
--- a/core/include/detray/tools/grid_builder.hpp
+++ b/core/include/detray/tools/grid_builder.hpp
@@ -39,8 +39,8 @@ class grid_builder final : public volume_decorator<detector_t> {
 
     // TODO: nullptr can lead to exceptions, remove in the full implementation
     DETRAY_HOST
-    grid_builder(std::unique_ptr<volume_builder_interface<detector_t>>
-                     vol_builder = nullptr)
+    grid_builder(
+        std::unique_ptr<volume_builder_interface<detector_t>> vol_builder)
         : volume_decorator<detector_t>(std::move(vol_builder)) {}
 
     /// Should the passive surfaces be added to the grid ?
@@ -49,17 +49,28 @@ class grid_builder final : public volume_decorator<detector_t> {
     }
 
     /// Delegate init call depending on @param span type
-    template <typename axis_span_t>
+    template <typename grid_shape_t>
     DETRAY_HOST void init_grid(
-        const axis_span_t &span,
+        const mask<grid_shape_t> &bounds,
         const std::array<std::size_t, grid_t::Dim> &n_bins,
         const std::array<std::vector<scalar_type>, grid_t::Dim> &ax_bin_edges =
             {{}}) {
-        init_impl(span, n_bins, ax_bin_edges,
+        init_impl(bounds, n_bins, ax_bin_edges,
                   typename grid_t::axes_type::bounds{},
                   typename grid_t::axes_type::binnings{});
     }
 
+    /// Build the empty grid from axis parameters
+    DETRAY_HOST void init_grid(
+        const std::vector<scalar_type> &spans,
+        const std::vector<std::size_t> &n_bins,
+        const std::vector<std::vector<scalar_type>> &ax_bin_edges = {{}}) {
+
+        m_grid = m_factory.template new_grid<typename grid_t::local_frame_type>(
+            spans, n_bins, ax_bin_edges, typename grid_t::axes_type::bounds{},
+            typename grid_t::axes_type::binnings{});
+    }
+
     /// Fill grid from existing volume using a bin filling strategy
     /// This can also be called without a volume builder
     template <typename volume_type, typename... Args>
@@ -83,53 +94,79 @@ class grid_builder final : public volume_decorator<detector_t> {
         bin_filler(m_grid, vol, surfaces, transforms, masks, ctx, args...);
     }
 
-    /// Overwrite, to add the sensitives to the grid, instead of the surface vec
-    DETRAY_HOST
-    void add_sensitives(
-        std::shared_ptr<surface_factory_interface<detector_t>> sf_factory,
-        typename detector_t::geometry_context ctx = {}) override {
-        (*sf_factory)(volume_decorator<detector_t>::operator()(), m_surfaces,
-                      m_transforms, m_masks, ctx);
-    }
-
-    /// Overwrite, to add the passives to the grid, instead of the surface vec
-    DETRAY_HOST
-    void add_passives(
-        std::shared_ptr<surface_factory_interface<detector_t>> ps_factory,
-        typename detector_t::geometry_context ctx = {}) override {
-        if (m_add_passives) {
-            (*ps_factory)(volume_decorator<detector_t>::operator()(),
-                          m_surfaces, m_transforms, m_masks, ctx);
-            ps_factory->clear();
-        } else {
-            volume_decorator<detector_t>::add_passives(std::move(ps_factory),
-                                                       ctx);
-        }
-    }
-
     /// Add the volume and the grid to the detector @param det
     DETRAY_HOST
     auto build(detector_t &det, typename detector_t::geometry_context ctx = {})
         -> typename detector_t::volume_type * override {
+
+        using surface_desc_t = typename detector_t::surface_type;
+
         // Add the surfaces (portals and/or passives) that are owned by the vol
         typename detector_t::volume_type *vol_ptr =
             volume_decorator<detector_t>::build(det, ctx);
-        m_bin_filler(m_grid, detector_volume{det, *vol_ptr}, m_surfaces,
-                     m_transforms, m_masks, ctx);
-
-        // Add the surfaces that were filled into the grid directly to the
-        // detector and update their links
-        const auto trf_offset{det.transform_store().size(ctx)};
-        for (auto &sf_desc : m_grid.all()) {
-            const auto sf = surface{det, sf_desc};
-            sf.template visit_mask<detail::mask_index_update>(sf_desc);
-            sf_desc.set_volume(vol_ptr->index());
-            sf_desc.update_transform(trf_offset);
-        }
 
-        // Add transforms and masks to detector
-        det.append_masks(std::move(m_masks));
-        det.append_transforms(std::move(m_transforms));
+        // Take the surfaces that should be filled into the grid out of the
+        // brute force finder
+        const auto vol_idx{vol_ptr->index()};
+        constexpr auto bf_id{detector_t::sf_finders::id::e_brute_force};
+        auto &bf_search = det.surface_store().template get<bf_id>();
+
+        // Grid has not been filled previously, fill it automatically
+        if (m_grid.size() == 0u) {
+            std::vector<surface_desc_t> surfaces{};
+            for (auto itr = bf_search.all().begin();
+                 itr != bf_search.all().end();) {
+                if (itr->volume() != vol_idx) {
+                    continue;
+                }
+                if (itr->is_sensitive() or
+                    (m_add_passives and itr->is_passive())) {
+                    surfaces.push_back(*itr);
+                    bf_search.erase(itr);
+                } else {
+                    ++itr;
+                }
+            }
+
+            this->fill_grid(
+                detector_volume{det,
+                                volume_decorator<detector_t>::operator()()},
+                surfaces, det.transform_store(), det.mask_store(), ctx);
+        } else {
+            // The grid is prefilled with surface descriptors that contain the
+            // correct surface indices per bin (e.g. from file IO).
+            // Now add the rest of the linking information, which is only
+            // available after the volume builder ran
+            for (auto itr = bf_search.all().begin();
+                 itr != bf_search.all().end();) {
+                if (itr->volume() != vol_idx) {
+                    ++itr;
+                    continue;
+                }
+                if (itr->is_sensitive() or
+                    (m_add_passives and itr->is_passive())) {
+                    const auto vol = det.volume_by_index(itr->volume());
+                    // The current volume is already built, so the surface
+                    // interface is safe to use
+                    const auto sf = surface{det, *itr};
+                    const auto t = sf.center(ctx);
+                    const auto loc_pos =
+                        m_grid.global_to_local(vol.transform(), t, t);
+                    auto bin_content = m_grid.search(loc_pos);
+
+                    for (surface_desc_t &sf_desc : bin_content) {
+                        // Find the correct surface and update all links
+                        if (sf_desc.index() == sf.index()) {
+                            sf_desc = *itr;
+                        }
+                    }
+
+                    bf_search.erase(itr);
+                } else {
+                    ++itr;
+                }
+            }
+        }
 
         // Add the grid to the detector and link it to its volume
         constexpr auto gid{detector_t::sf_finders::template get_id<grid_t>()};
@@ -159,30 +196,10 @@ class grid_builder final : public volume_decorator<detector_t> {
             bounds, n_bins, ax_bin_edges);
     }
 
-    /// Build the empty grid from axis parameters
-    template <typename grid_shape_t, typename... axis_bounds,
-              typename... binning_ts>
-    DETRAY_HOST void init_impl(
-        const std::vector<scalar_type> &spans,
-        const std::vector<std::size_t> *n_bins,
-        const std::vector<std::vector<scalar_type>> &ax_bin_edges,
-        std::tuple<axis_bounds...>, std::tuple<binning_ts...>) {
-
-        m_grid =
-            m_factory
-                .template new_grid<grid_shape_t, axis_bounds..., binning_ts...>(
-                    spans, n_bins, ax_bin_edges);
-    }
-
     grid_factory_t m_factory{};
     typename grid_t::template type<true> m_grid{};
     bin_filler_t m_bin_filler{};
     bool m_add_passives{false};
-
-    // surfaces that are filled into the grid, but not the volume
-    typename detector_t::surface_container_t m_surfaces{};
-    typename detector_t::transform_container m_transforms{};
-    typename detector_t::mask_container m_masks{};
 };
 
 /// Grid builder from single components
diff --git a/core/include/detray/tools/grid_factory.hpp b/core/include/detray/tools/grid_factory.hpp
index 9de78548e..17b94cfb7 100644
--- a/core/include/detray/tools/grid_factory.hpp
+++ b/core/include/detray/tools/grid_factory.hpp
@@ -23,6 +23,7 @@
 #include <cassert>
 #include <iostream>
 #include <string>
+#include <type_traits>
 
 namespace detray {
 
@@ -58,13 +59,6 @@ class grid_factory {
     explicit grid_factory(vecmem::memory_resource &resource)
         : m_resource(&resource) {}
 
-    /// Get new grid collection
-    template <typename grid_t>
-    constexpr auto new_collection() {
-        return grid_collection<typename grid_t::template type<!is_owning>>(
-            m_resource);
-    }
-
     /// Print grid - up to three dimensions
     /// @note will likely become obsolete with the actsvg implementation
     template <typename grid_t>
@@ -305,6 +299,7 @@ class grid_factory {
         // Axes boundaries and local indices
         using boundary = trapezoid2D<>::boundaries;
         using axes = trapezoid2D<>::axes<>;
+
         constexpr auto e_x_axis = static_cast<dindex>(axes::axis_loc0);
         constexpr auto e_y_axis = static_cast<dindex>(axes::axis_loc1);
 
@@ -336,22 +331,54 @@ class grid_factory {
     /// @param ax_bin_edges the explicit bin edges for irregular axes
     ///                     (lower bin edges + the the upper edge of the
     ///                     last bin), otherwise ignored.
-    template <typename grid_shape_t, typename... bound_ts,
-              typename... binning_ts>
+    template <
+        typename grid_shape_t, typename... bound_ts, typename... binning_ts,
+        std::enable_if_t<std::is_enum_v<typename grid_shape_t::boundaries>,
+                         bool> = true>
+    auto new_grid(const std::vector<scalar_type> spans,
+                  const std::vector<std::size_t> n_bins,
+                  const std::vector<std::vector<scalar_type>> &ax_bin_edges,
+                  const std::tuple<bound_ts...> &bounds,
+                  const std::tuple<binning_ts...> &binnings) const {
+
+        return new_grid<
+            typename grid_shape_t::template local_frame_type<algebra_t>>(
+            spans, n_bins, ax_bin_edges, bounds, binnings);
+    }
+
+    /// @brief Create and empty grid with fully initialized axes.
+    ///
+    /// @tparam grid_shape_t the shape of the resulting grid
+    ///         (e.g. cylinder2D).
+    /// @tparam e_bounds the bounds of the regular axes
+    ///         (open vs. closed bounds).
+    /// @tparam binnings the binning types of the axes
+    ///         (regular vs. irregular)
+    ///
+    /// @param spans the span of the axis values for regular axes, otherwise
+    ///              ignored.
+    /// @param n_bins the number of bins for regular axes, otherwise ignored
+    /// @param ax_bin_edges the explicit bin edges for irregular axes
+    ///                     (lower bin edges + the the upper edge of the
+    ///                     last bin), otherwise ignored.
+    template <
+        typename grid_frame_t, typename... bound_ts, typename... binning_ts,
+        std::enable_if_t<
+            std::is_object_v<typename grid_frame_t::rotation_matrix>, bool> =
+            true>
     auto new_grid(const std::vector<scalar_type> spans,
                   const std::vector<std::size_t> n_bins,
                   const std::vector<std::vector<scalar_type>> &ax_bin_edges,
-                  const std::tuple<bound_ts...>,
-                  const std::tuple<binning_ts...>) const {
+                  const std::tuple<bound_ts...> & = {},
+                  const std::tuple<binning_ts...> & = {}) const {
 
         static_assert(sizeof...(bound_ts) == sizeof...(binning_ts),
                       "number of axis bounds and binning types has to match");
 
         // Build the coordinate axes and the grid
-        using axes_t = n_axis::multi_axis<
-            is_owning,
-            typename grid_shape_t::template local_frame_type<algebra_t>,
-            n_axis::single_axis<bound_ts, binning_ts>...>;
+        using axes_t =
+            n_axis::multi_axis<is_owning, grid_frame_t,
+                               n_axis::single_axis<bound_ts, binning_ts>...>;
         using bin_t = typename grid_type<axes_t>::bin_type;
 
         // Prepare data
diff --git a/core/include/detray/tools/surface_factory.hpp b/core/include/detray/tools/surface_factory.hpp
index 078ec16fa..6fdda35b0 100644
--- a/core/include/detray/tools/surface_factory.hpp
+++ b/core/include/detray/tools/surface_factory.hpp
@@ -17,6 +17,7 @@
 #include "detray/utils/ranges.hpp"
 
 // System include(s)
+#include <algorithm>
 #include <cassert>
 #include <limits>
 #include <memory>
@@ -99,7 +100,7 @@ class surface_factory
     DETRAY_HOST
     void push_back(surface_data_t &&sf_data) override {
 
-        auto [trf, vlink, bounds] = sf_data.get_data();
+        auto [vlink, bounds, trf, index] = sf_data.get_data();
 
         assert(bounds.size() == mask_shape_t::boundaries::e_size);
 
@@ -113,6 +114,7 @@ class surface_factory
             m_volume_link.push_back(vlink);
         }
 
+        m_indices.push_back(index);
         m_transforms.push_back(trf);
         m_components.push_back(std::move(bounds));
     }
@@ -127,6 +129,7 @@ class surface_factory
         const auto n_surfaces{
             static_cast<dindex>(size() + surface_data.size())};
 
+        m_indices.reserve(n_surfaces);
         m_transforms.reserve(n_surfaces);
         m_components.reserve(n_surfaces);
 
@@ -145,6 +148,7 @@ class surface_factory
     /// Clear old data
     DETRAY_HOST
     auto clear() -> void override {
+        m_indices.clear();
         m_components.clear();
         m_transforms.clear();
         // cannot clear the single-view
@@ -176,7 +180,7 @@ class surface_factory
         // The material will be added in a later step
         constexpr auto no_material = surface_t::material_id::e_none;
         // In case the surfaces container is prefilled with other surfaces
-        dindex surfaces_offset = static_cast<dindex>(surfaces.size());
+        const dindex surfaces_offset = static_cast<dindex>(surfaces.size());
 
         // Nothing to construct
         if (size() == 0u) {
@@ -185,8 +189,15 @@ class surface_factory
 
         for (const auto [idx, comp] : detray::views::enumerate(m_components)) {
 
+            // Append the surfaces relative to the current number of surfaces
+            // in the stores
+            const dindex sf_idx{is_invalid_value(m_indices[idx])
+                                    ? dindex_invalid
+                                    : m_indices[idx]};
+
             // Add transform
-            transforms.push_back(m_transforms[idx], ctx);
+            const dindex trf_idx =
+                insert_in_container(transforms, m_transforms[idx], sf_idx, ctx);
 
             if constexpr (std::is_same_v<mask_shape_t, unmasked>) {
                 masks.template emplace_back<mask_id>(
@@ -200,20 +211,42 @@ class surface_factory
             // Add surface with all links set (relative to the given containers)
             mask_link_t mask_link{mask_id, masks.template size<mask_id>() - 1u};
             material_link_t material_link{no_material, dindex_invalid};
-            surfaces.emplace_back(transforms.size(ctx) - 1u, mask_link,
-                                  material_link, volume.index(), dindex_invalid,
-                                  sf_id);
+
+            insert_in_container(surfaces,
+                                {trf_idx, mask_link, material_link,
+                                 volume.index(), dindex_invalid, sf_id},
+                                sf_idx);
         }
 
         return {surfaces_offset, static_cast<dindex>(surfaces.size())};
     }
 
     private:
+    template <typename container_t, typename... Args>
+    DETRAY_HOST dindex insert_in_container(
+        container_t &cont, const typename container_t::value_type value,
+        const dindex idx, Args &&... args) const {
+        if (is_invalid_value(idx)) {
+            cont.push_back(value, std::forward<Args>(args)...);
+
+            return static_cast<dindex>(cont.size() - 1u);
+        } else {
+            if (cont.size() < idx + 1) {
+                cont.resize(idx + 1, std::forward<Args>(args)...);
+            }
+
+            cont.at(idx, std::forward<Args>(args)...) = value;
+
+            return idx;
+        }
+    }
+
     /// Check that the containers have the same size
     DETRAY_HOST
     void check() const {
         // This should not happend (need same number and ordering of data)
         assert(m_components.size() == m_transforms.size());
+        assert(m_components.size() == m_indices.size());
         if constexpr (sf_id == surface_id::e_portal) {
             assert(m_components.size() == m_volume_link.size());
         } else {
@@ -221,6 +254,8 @@ class surface_factory
         }
     }
 
+    dindex_range m_range{dindex_invalid, dindex_invalid};
+    std::vector<dindex> m_indices{};
     std::vector<std::vector<scalar_t>> m_components{};
     std::vector<typename detector_t::transform3> m_transforms{};
     volume_link_collection m_volume_link{};
diff --git a/core/include/detray/tools/surface_factory_interface.hpp b/core/include/detray/tools/surface_factory_interface.hpp
index 851c62f5f..a133399d3 100644
--- a/core/include/detray/tools/surface_factory_interface.hpp
+++ b/core/include/detray/tools/surface_factory_interface.hpp
@@ -33,26 +33,32 @@ class surface_data {
     DETRAY_HOST
     surface_data(
         const typename detector_t::transform3 &trf, navigation_link volume_link,
-        const std::vector<typename detector_t::scalar_type> &mask_boundaries)
-        : m_transform{trf},
-          m_volume_link{volume_link},
-          m_boundaries{mask_boundaries} {}
+        const std::vector<typename detector_t::scalar_type> &mask_boundaries,
+        const dindex idx = dindex_invalid)
+        : m_volume_link{volume_link},
+          m_index{idx},
+          m_boundaries{mask_boundaries},
+          m_transform{trf} {}
 
     DETRAY_HOST
     auto get_data()
-        -> std::tuple<typename detector_t::transform3 &, navigation_link &,
-                      std::vector<typename detector_t::scalar_type> &> {
-        return std::tie(m_transform, m_volume_link, m_boundaries);
+        -> std::tuple<navigation_link &,
+                      std::vector<typename detector_t::scalar_type> &,
+                      typename detector_t::transform3 &, dindex &> {
+        return std::tie(m_volume_link, m_boundaries, m_transform, m_index);
     }
 
     private:
-    /// The surface placement
-    typename detector_t::transform3 m_transform;
     /// The index of the volume that this surface links to
     navigation_link m_volume_link;
-    // simple tuple of all mask types in the detector. Only one entry is filled
+    /// The position of the surface in the detector containers, used to match
+    /// the surface to e.g. its material
+    dindex m_index;
+    // Simple tuple of all mask types in the detector. Only one entry is filled
     // with the mask that corresponds to this specific surface.
     std::vector<typename detector_t::scalar_type> m_boundaries;
+    /// The surface placement
+    typename detector_t::transform3 m_transform;
 };
 
 /// @brief How to generate surfaces with their corresponding masks and
diff --git a/core/include/detray/tools/volume_builder.hpp b/core/include/detray/tools/volume_builder.hpp
index e4fa2c7b3..39d6610e4 100644
--- a/core/include/detray/tools/volume_builder.hpp
+++ b/core/include/detray/tools/volume_builder.hpp
@@ -147,7 +147,7 @@ class volume_builder : public volume_builder_interface<detector_t> {
 
         // Update mask and transform index of surfaces and set a
         // unique barcode (index of surface in container)
-        auto sf_offset{static_cast<dindex>(det.portals().size())};
+        auto sf_offset{static_cast<dindex>(det.surface_lookup().size())};
         for (auto& sf_desc : m_surfaces) {
             const auto sf = surface{det, sf_desc};
             sf.template visit_mask<detail::mask_index_update>(sf_desc);
diff --git a/core/include/detray/utils/consistency_checker.hpp b/core/include/detray/utils/consistency_checker.hpp
index 0044f096c..26b91ab92 100644
--- a/core/include/detray/utils/consistency_checker.hpp
+++ b/core/include/detray/utils/consistency_checker.hpp
@@ -209,6 +209,7 @@ inline bool check_consistency(const detector_t &det) {
 
         // Check that nothing is obviously broken
         if (not sf.self_check(err_stream)) {
+            err_stream << "\nat surface no. " << std::to_string(idx);
             throw std::invalid_argument(err_stream.str());
         }
 
diff --git a/core/include/detray/utils/inspectors.hpp b/core/include/detray/utils/inspectors.hpp
index db324dab3..ee33a5654 100644
--- a/core/include/detray/utils/inspectors.hpp
+++ b/core/include/detray/utils/inspectors.hpp
@@ -16,7 +16,6 @@
 
 // System include(s)
 #include <iomanip>
-#include <iostream>
 #include <sstream>
 #include <string>
 #include <tuple>
diff --git a/core/include/detray/utils/type_list.hpp b/core/include/detray/utils/type_list.hpp
new file mode 100644
index 000000000..e5a6cefd2
--- /dev/null
+++ b/core/include/detray/utils/type_list.hpp
@@ -0,0 +1,170 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s)
+#include "detray/utils/tuple.hpp"
+#include "detray/utils/type_traits.hpp"
+
+// System include(s)
+#include <string>
+#include <string_view>
+#include <type_traits>
+
+namespace detray::types {
+
+/// @brief type list implementation
+template <typename... Ts>
+struct type_list {};
+
+/// Number of types in the list
+/// @{
+template <typename = void>
+struct get_size {};
+
+template <typename... Ts>
+struct get_size<type_list<Ts...>>
+    : std::integral_constant<std::size_t, sizeof...(Ts)> {};
+
+template <typename L>
+constexpr inline std::size_t size{get_size<L>()};
+/// @}
+
+/// Access the first type
+/// @{
+template <typename = void>
+struct get_front {};
+
+template <typename T, typename... Ts>
+struct get_front<type_list<T, Ts...>> {
+    using type = T;
+};
+template <typename L>
+using front = typename get_front<L>::type;
+/// @}
+
+/// Access the last type
+/// @{
+template <typename = void>
+struct get_back {};
+
+template <typename T, typename... Ts>
+struct get_back<type_list<T, Ts...>> {
+    using type = std::conditional_t<sizeof...(Ts) == 0, T,
+                                    typename get_back<type_list<Ts...>>::type>;
+};
+// Base case
+template <>
+struct get_back<type_list<>> {
+    using type = void;
+};
+template <typename L>
+using back = typename get_back<L>::type;
+/// @}
+
+/// Access the N-th type
+/// @{
+template <int N, typename = void>
+struct get_at {};
+
+template <int N, typename T, typename... Ts>
+struct get_at<N, type_list<T, Ts...>> {
+    using type = detail::remove_cvref_t<decltype(
+        detray::get<N>(detray::tuple<T, Ts...>{}))>;
+};
+template <typename L, int N>
+using at = typename get_at<N, L>::type;
+/// @}
+
+/// Append a type
+/// @{
+template <typename N, typename = void>
+struct do_push_back {};
+
+template <typename N, typename... Ts>
+struct do_push_back<N, type_list<Ts...>> {
+    using type = type_list<Ts..., N>;
+};
+template <typename L, typename N>
+using push_back = typename do_push_back<N, L>::type;
+/// @}
+
+/// Prepend a type
+/// @{
+template <typename N, typename = void>
+struct do_push_front {};
+
+template <typename N, typename... Ts>
+struct do_push_front<N, type_list<Ts...>> {
+    using type = type_list<N, Ts...>;
+};
+template <typename L, typename N>
+using push_front = typename do_push_front<N, L>::type;
+/// @}
+
+/// Print the type list
+/// @{
+
+/// @see
+/// https://stackoverflow.com/questions/281818/unmangling-the-result-of-stdtype-infoname
+template <typename T>
+std::string demangle_type_name() {
+#if defined(__clang__)
+    constexpr std::string_view prefix{"[T = "};
+    constexpr std::string_view suffix{"; "};
+    constexpr std::string_view function{__PRETTY_FUNCTION__};
+#elif defined(__GNUC__)
+    constexpr std::string_view prefix{"with T = "};
+    constexpr std::string_view suffix{"; "};
+    constexpr std::string_view function{__PRETTY_FUNCTION__};
+#elif defined(_MSC_VER)
+    constexpr std::string_view prefix{"get_type_name<"};
+    constexpr std::string_view suffix{">(void)"};
+    constexpr std::string_view function{__FUNCSIG__};
+#else
+    std::string type_name{typeid(T).name()};
+    type_name += " ";
+    std::string_view function{type_name};
+    constexpr std::string_view prefix{""};
+    constexpr std::string_view suffix{" "};
+#endif
+
+    const std::size_t start{function.find(prefix) + prefix.size()};
+    const std::size_t end{function.find(suffix)};
+    const std::size_t size{end - start};
+
+    return std::string{function.substr(start, size)};
+}
+
+template <typename = void>
+struct print {};
+
+template <typename... Ts>
+struct print<type_list<Ts...>> {
+
+    template <typename P = void, typename... Ps>
+    void print_typeid() {
+
+        std::printf("%s", demangle_type_name<P>().c_str());
+
+        // Keep unrolling the pack
+        if constexpr (sizeof...(Ps) > 0) {
+            std::printf(", ");
+            return print_typeid<Ps...>();
+        }
+    }
+
+    print() {
+        std::printf("type_list<");
+        print_typeid<Ts...>();
+        std::printf(">\n");
+    }
+};
+/// @}
+
+}  // namespace detray::types
diff --git a/io/include/detray/io/common/detail/detector_components_io.hpp b/io/include/detray/io/common/detail/detector_components_io.hpp
index 16a7978b8..4d7caa079 100644
--- a/io/include/detray/io/common/detail/detector_components_io.hpp
+++ b/io/include/detray/io/common/detail/detector_components_io.hpp
@@ -41,12 +41,12 @@ class detector_component_writers final {
     detector_component_writers() = default;
 
     /// Create a new writer of type @tparam writer_t
-    template <template <typename> class writer_t,
+    template <template <typename...> class writer_t, typename... Args,
               std::enable_if_t<std::is_base_of_v<writer_interface<detector_t>,
-                                                 writer_t<detector_t>>,
+                                                 writer_t<detector_t, Args...>>,
                                bool> = true>
     void add() {
-        add(std::make_unique<writer_t<detector_t>>());
+        add(std::make_unique<writer_t<detector_t, Args...>>());
     }
 
     /// Attach an existing writer via @param w_ptr to the writers
@@ -92,12 +92,12 @@ class detector_component_readers final {
     detector_component_readers() = default;
 
     /// Create a new reader of type @tparam reader_t
-    template <template <typename> class reader_t,
+    template <template <typename...> class reader_t, typename... Args,
               std::enable_if_t<std::is_base_of_v<reader_interface<detector_t>,
-                                                 reader_t<detector_t>>,
+                                                 reader_t<detector_t, Args...>>,
                                bool> = true>
     void add(const std::string& file_name) {
-        add(std::make_unique<reader_t<detector_t>>(), file_name);
+        add(std::make_unique<reader_t<detector_t, Args...>>(), file_name);
     }
 
     /// Attach an existing reader via @param w_ptr to the readers
diff --git a/io/include/detray/io/common/detector_reader.hpp b/io/include/detray/io/common/detector_reader.hpp
index b2cb7273a..a7cd4ff46 100644
--- a/io/include/detray/io/common/detector_reader.hpp
+++ b/io/include/detray/io/common/detector_reader.hpp
@@ -78,7 +78,7 @@ namespace detail {
 
 /// From the list of file that are given as part of the config @param cfg,
 /// infer the readers that are needed by peeking into the file header
-template <class detector_t>
+template <class detector_t, std::size_t CAP, std::size_t DIM>
 auto assemble_reader(const io::detector_reader_config& cfg) {
 
     detail::detector_component_readers<detector_t> readers;
@@ -99,6 +99,12 @@ auto assemble_reader(const io::detector_reader_config& cfg) {
             } else if (header.tag == "homogeneous_material") {
                 readers.template add<json_homogeneous_material_reader>(
                     file_name);
+            } else if (header.tag == "surface_grids") {
+                using surface_t = typename detector_t::surface_type;
+                readers.template add<json_grid_reader, surface_t,
+                                     std::integral_constant<std::size_t, CAP>,
+                                     std::integral_constant<std::size_t, DIM>>(
+                    file_name);
             } else {
                 throw std::invalid_argument("Unsupported file tag '" +
                                             header.tag +
@@ -134,7 +140,7 @@ namespace io {
 /// @param cfg the detector reader configuration
 ///
 /// @returns a complete detector object + a map that contains the volume names
-template <class detector_t,
+template <class detector_t, std::size_t CAP = 9u, std::size_t DIM = 2u,
           template <typename> class volume_builder_t = volume_builder>
 auto read_detector(vecmem::memory_resource& resc,
                    const detector_reader_config& cfg) {
@@ -149,7 +155,8 @@ auto read_detector(vecmem::memory_resource& resc,
         det_builder;
 
     // Find all required readers
-    auto [reader, det_name] = detray::detail::assemble_reader<detector_t>(cfg);
+    auto [reader, det_name] =
+        detray::detail::assemble_reader<detector_t, CAP, DIM>(cfg);
 
     // Add a constant b-field (no bfield reader is added in this case)
     if constexpr (std::is_same_v<bfield_bknd_t, bfield::const_bknd_t>) {
diff --git a/io/include/detray/io/common/geometry_reader.hpp b/io/include/detray/io/common/geometry_reader.hpp
index e70daea21..f42b080d5 100644
--- a/io/include/detray/io/common/geometry_reader.hpp
+++ b/io/include/detray/io/common/geometry_reader.hpp
@@ -76,8 +76,6 @@ class geometry_reader : public reader_interface<detector_t> {
         typename detector_t::name_map& name_map,
         const detector_payload& det_data) {
 
-        // @todo Add volume grid
-
         // Deserialize the volumes one-by-one
         for (const auto& vol_data : det_data.volumes) {
             // Get a generic volume builder first and decorate it later
@@ -87,8 +85,8 @@ class geometry_reader : public reader_interface<detector_t> {
             // Set the volume name
             name_map[vbuilder->vol_index() + 1u] = vol_data.name;
 
-            // @todo add the volume placement, once it can be checked for the
-            // test detectors
+            // Volume placement
+            vbuilder->add_volume_placement(deserialize(vol_data.transform));
 
             // Prepare the surface factories (one per shape and surface type)
             std::map<std::pair<surface_id, mask_shape>, sf_factory_ptr_t>
@@ -172,10 +170,14 @@ class geometry_reader : public reader_interface<detector_t> {
                   sf_data.mask.boundaries.end(),
                   std::back_inserter(mask_boundaries));
 
+        const std::size_t sf_idx{sf_data.index_in_coll.has_value()
+                                     ? *(sf_data.index_in_coll)
+                                     : detail::invalid_value<std::size_t>()};
+
         return {deserialize(sf_data.transform),
                 static_cast<nav_link_t>(
                     base_type::deserialize(sf_data.mask.volume_link)),
-                std::move(mask_boundaries)};
+                std::move(mask_boundaries), static_cast<dindex>(sf_idx)};
     }
 
     private:
@@ -214,7 +216,7 @@ class geometry_reader : public reader_interface<detector_t> {
                                             surface_id::e_passive>;
                         return std::make_shared<ps_factory_t>();
                     case surface_id::e_unknown:
-                        throw std::runtime_error(
+                        throw std::invalid_argument(
                             "Unknown surface type in geometry file");
                 };
             }
@@ -227,10 +229,11 @@ class geometry_reader : public reader_interface<detector_t> {
         }
         // Test some edge cases
         if (shape_id == mask_shape::unknown) {
-            throw std::runtime_error("Unknown mask shape in geometry file!");
+            throw std::invalid_argument(
+                "Unknown mask shape id in geometry file!");
         } else {
-            throw std::runtime_error(
-                "Given mask type could not be matched: " +
+            throw std::invalid_argument(
+                "Given shape id could not be matched to a mask type: " +
                 std::to_string(static_cast<std::int64_t>(shape_id)));
 
             // Cannot be reached
diff --git a/io/include/detray/io/common/geometry_writer.hpp b/io/include/detray/io/common/geometry_writer.hpp
index ccd9c1ba7..bfe44fb99 100644
--- a/io/include/detray/io/common/geometry_writer.hpp
+++ b/io/include/detray/io/common/geometry_writer.hpp
@@ -109,9 +109,11 @@ class geometry_writer : public writer_interface<detector_t> {
     }
 
     /// Serialize a detector surface @param sf into its io payload
-    static surface_payload serialize(const surface<detector_t>& sf) {
+    static surface_payload serialize(const surface<detector_t>& sf,
+                                     std::size_t sf_idx) {
         surface_payload sf_data;
 
+        sf_data.index_in_coll = sf_idx;
         sf_data.type = sf.id();
         sf_data.barcode = sf.barcode().value();
         sf_data.transform = serialize(sf.transform({}));
@@ -134,9 +136,13 @@ class geometry_writer : public writer_interface<detector_t> {
             serialize(det.transform_store()[vol_desc.transform()]);
         vol_data.type = vol_desc.id();
 
+        // Count the surfaces belonging to this volume
+        std::size_t sf_idx{0};
+
         for (const auto& sf_desc : det.surface_lookup()) {
             if (sf_desc.volume() == vol_desc.index()) {
-                vol_data.surfaces.push_back(serialize(surface{det, sf_desc}));
+                vol_data.surfaces.push_back(
+                    serialize(surface{det, sf_desc}, sf_idx++));
             }
         }
 
diff --git a/io/include/detray/io/common/grid_reader.hpp b/io/include/detray/io/common/grid_reader.hpp
new file mode 100644
index 000000000..5b819f9be
--- /dev/null
+++ b/io/include/detray/io/common/grid_reader.hpp
@@ -0,0 +1,356 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s)
+#include "detray/definitions/indexing.hpp"
+#include "detray/io/common/detail/type_traits.hpp"
+#include "detray/io/common/io_interface.hpp"
+#include "detray/io/common/payloads.hpp"
+#include "detray/tools/detector_builder.hpp"
+#include "detray/tools/grid_builder.hpp"
+#include "detray/tools/grid_factory.hpp"
+#include "detray/utils/ranges.hpp"
+#include "detray/utils/type_list.hpp"
+
+// System include(s)
+#include <queue>
+#include <stdexcept>
+#include <string>
+#include <type_traits>
+#include <vector>
+
+namespace detray {
+
+/// @brief Abstract base class for surface grid readers
+template <class detector_t,
+          typename value_t = typename detector_t::surface_type,
+          typename CAP = std::integral_constant<std::size_t, 9>,
+          typename DIM = std::integral_constant<std::size_t, 2>>
+class grid_reader : public reader_interface<detector_t> {
+
+    using base_type = reader_interface<detector_t>;
+    /// IO accelerator ids do not need to coincide with the detector ids,
+    /// because they are shared with ACTS
+    using acc_type = io::detail::acc_type;
+    using algebra_t = typename detector_t::transform3;
+    using scalar_t = typename algebra_t::scalar_type;
+
+    // Bin content type to be written into the grid bin payloads
+    // For detector surface descriptors, write only the surface index
+    using content_t = std::conditional_t<
+        std::is_same_v<value_t, typename detector_t::surface_type>, std::size_t,
+        value_t>;
+
+    static constexpr std::size_t Dim{DIM()};
+    static constexpr std::size_t bin_capacity{CAP()};
+
+    protected:
+    /// Tag the reader as "surface_grids"
+    inline static const std::string tag = "surface_grids";
+
+    public:
+    /// Same constructors for this class as for base_type
+    using base_type::base_type;
+
+    protected:
+    /// Deserialize the detector grids @param grids_data from their IO
+    /// payload
+    static void deserialize(
+        detector_builder<typename detector_t::metadata,
+                         typename detector_t::bfield_type::backend_t,
+                         volume_builder> &det_builder,
+        typename detector_t::name_map &,
+        const detector_grids_payload<content_t> &grids_data) {
+
+        // Deserialize the grids volume by volume
+        for (const auto &grid_data : grids_data.grids) {
+
+            std::queue<n_axis::bounds> bounds;
+            std::queue<n_axis::binning> binnings;
+            for (const auto &axis_data : grid_data.axes) {
+                bounds.push(axis_data.bounds);
+                binnings.push(axis_data.binning);
+            }
+            deserialize(bounds, binnings, grid_data, det_builder);
+        }
+    }
+
+    private:
+    /// @brief recursively build the grid: axis bounds (open, closed, circular)
+    ///
+    /// @tparam bounds_ts type list that contains the bounds types that were
+    ///         identified from the IO ids so far (start with empty list)
+    /// @tparam binning_ts type list that contains the binning types that were
+    ///         identified from the IO ids so far (start with empty list)
+    ///
+    /// @param bound_ids runtime queue of bounds type ids (read from file)
+    /// @param binning_ids runtime queue of binning type ids (read from file)
+    template <typename bounds_ts = types::type_list<>,
+              typename binning_ts = types::type_list<>, typename... Ts>
+    static void deserialize(std::queue<n_axis::bounds> &bound_ids,
+                            std::queue<n_axis::binning> &binning_ids,
+                            Ts &&... data) {
+        using namespace n_axis;
+
+        constexpr std::size_t n_bounds_types{types::size<bounds_ts>};
+
+        // Base case: If the bounds types are filled, continue with the binnings
+        if constexpr (n_bounds_types == Dim) {
+            return deserialize<bounds_ts, binning_ts>(
+                binning_ids, std::forward<Ts>(data)...);
+        } else if (!bound_ids.empty()) {
+            // The axis label, e.g. x, y or z by number
+            constexpr auto lb{static_cast<label>(n_bounds_types)};
+
+            const auto first_id{bound_ids.front()};
+            bound_ids.pop();
+
+            // Based on the type id, add the next bounds type to the type list
+            // and continue
+            switch (first_id) {
+                case bounds::e_closed: {
+                    using new_bounds_ts =
+                        types::push_back<bounds_ts, closed<lb>>;
+                    return deserialize<new_bounds_ts, binning_ts>(
+                        bound_ids, binning_ids, std::forward<Ts>(data)...);
+                }
+                case bounds::e_open: {
+                    using new_bounds_ts = types::push_back<bounds_ts, open<lb>>;
+                    return deserialize<new_bounds_ts, binning_ts>(
+                        bound_ids, binning_ids, std::forward<Ts>(data)...);
+                }
+                case bounds::e_circular: {
+                    using new_bounds_ts =
+                        types::push_back<bounds_ts, circular<lb>>;
+                    return deserialize<new_bounds_ts, binning_ts>(
+                        bound_ids, binning_ids, std::forward<Ts>(data)...);
+                }
+                // Test some edge cases
+                default: {
+                    throw std::invalid_argument(
+                        "Given type id could not be matched to a grid boundary "
+                        "type: " +
+                        std::to_string(static_cast<std::int64_t>(first_id)));
+                    break;
+                }
+            };
+        }
+    }
+
+    /// @brief recursively build the grid: axis binning (regular, irregular)
+    ///
+    /// @tparam bounds_ts type list that contains the bounds types
+    /// @tparam binning_ts type list that contains the binning types that were
+    ///         identified from the IO ids so far (start with empty list)
+    ///
+    /// @param binning_ids runtime queue of binning type ids (read from file)
+    template <typename bounds_ts, typename binning_ts, typename... Ts,
+              std::enable_if_t<types::size<bounds_ts> == Dim, bool> = true>
+    static void deserialize(std::queue<n_axis::binning> &binning_ids,
+                            Ts &&... data) {
+
+        using namespace n_axis;
+
+        using regular_binning_t = regular<host_container_types, scalar_t>;
+        using irregular_binning_t = irregular<host_container_types, scalar_t>;
+
+        // Base case: If the binning types are filled, continue with the frame
+        if constexpr (types::size<binning_ts> == Dim) {
+            return deserialize<bounds_ts, binning_ts>(
+                std::forward<Ts>(data)...);
+        } else if (!binning_ids.empty()) {
+
+            const auto first_id{binning_ids.front()};
+            binning_ids.pop();
+
+            switch (first_id) {
+                case binning::e_regular: {
+                    using new_binning_ts =
+                        types::push_back<binning_ts, regular_binning_t>;
+                    return deserialize<bounds_ts, new_binning_ts>(
+                        binning_ids, std::forward<Ts>(data)...);
+                }
+                case binning::e_irregular: {
+                    using new_binning_ts =
+                        types::push_back<binning_ts, irregular_binning_t>;
+                    return deserialize<bounds_ts, new_binning_ts>(
+                        binning_ids, std::forward<Ts>(data)...);
+                }
+                // Test some edge cases
+                default: {
+                    throw std::invalid_argument(
+                        "Given type id could not be matched to a grid binning "
+                        "type: " +
+                        std::to_string(static_cast<std::int64_t>(first_id)));
+                    break;
+                }
+            };
+        }
+    }
+
+    /// @brief recursively build the grid: find the grid geometry (loc. coord.)
+    ///
+    /// @tparam bounds_ts type list that contains the bounds types
+    /// @tparam binning_ts type list that contains the binning types
+    ///
+    /// @param grid_data grid IO payload (read from file)
+    /// @param det_builder gather the grid data and build the final volume
+    template <typename bounds_ts, typename binning_ts,
+              std::enable_if_t<types::size<bounds_ts> == Dim and
+                                   types::size<binning_ts> == Dim,
+                               bool> = true>
+    static void deserialize(
+        const grid_payload<content_t> &grid_data,
+        detector_builder<typename detector_t::metadata,
+                         typename detector_t::bfield_type::backend_t,
+                         volume_builder> &det_builder) {
+
+        // Throw expection if the accelerator link type id is invalid
+        auto print_error = [](io::detail::acc_type acc_link) -> void {
+            if (acc_link == io::detail::acc_type::unknown) {
+                throw std::invalid_argument(
+                    "Unknown accelerator id in geometry file!");
+            } else {
+                throw std::invalid_argument(
+                    "Given accelerator id could not be matched to a grid "
+                    "type: " +
+                    std::to_string(static_cast<std::int64_t>(acc_link)));
+            }
+        };
+
+        // Need to pass actual instances to deduce contained types as
+        // template parameter packs
+        constexpr auto bounds = bounds_ts{};
+        constexpr auto binnings = binning_ts{};
+
+        // Check only 2-dimensional grid types
+        if constexpr (Dim >= 2) {
+            switch (grid_data.acc_link.type) {
+                // rectangle, trapezoid, (triangle) grids
+                case io::detail::acc_type::cartesian2_grid: {
+                    return deserialize<cartesian2<algebra_t>>(
+                        grid_data, det_builder, bounds, binnings);
+                }
+                // ring/disc, annulus grids
+                case io::detail::acc_type::polar2_grid: {
+                    return deserialize<polar2<algebra_t>>(
+                        grid_data, det_builder, bounds, binnings);
+                }
+                // 2D cylinder grid
+                case io::detail::acc_type::cylinder2_grid: {
+                    return deserialize<cylindrical2<algebra_t>>(
+                        grid_data, det_builder, bounds, binnings);
+                }
+                default: {
+                    print_error(grid_data.acc_link.type);
+                    break;
+                }
+            };
+        } else if constexpr (Dim >= 3) {
+            switch (grid_data.acc_link.type) {
+                // cuboid grid
+                case io::detail::acc_type::cuboid3_grid: {
+                    return deserialize<cartesian3<algebra_t>>(
+                        grid_data, det_builder, bounds, binnings);
+                }
+                // 3D cylinder grid
+                case io::detail::acc_type::cylinder3_grid: {
+                    return deserialize<cylindrical3<algebra_t>>(
+                        grid_data, det_builder, bounds, binnings);
+                }
+                default: {
+                    print_error(grid_data.acc_link.type);
+                    break;
+                }
+            };
+        } else {
+            throw std::invalid_argument("No 1D grid type defined in detray");
+        }
+    }
+
+    /// @brief End of recursion: build the grid from the @param grid_data
+    template <typename local_frame_t, typename... bounds_ts,
+              typename... binning_ts,
+              std::enable_if_t<sizeof...(bounds_ts) == Dim and
+                                   sizeof...(binning_ts) == Dim,
+                               bool> = true>
+    static void deserialize(
+        const grid_payload<content_t> &grid_data,
+        detector_builder<typename detector_t::metadata,
+                         typename detector_t::bfield_type::backend_t,
+                         volume_builder> &det_builder,
+        types::type_list<bounds_ts...>, types::type_list<binning_ts...>) {
+        // Assemble the grid type
+        using populator_t = regular_attacher<bin_capacity>;
+        using axes_t =
+            n_axis::multi_axis<false, local_frame_t,
+                               n_axis::single_axis<bounds_ts, binning_ts>...>;
+        // Ok for now: Only have this serializer
+        using grid_t = grid<axes_t, value_t, simple_serializer, populator_t>;
+
+        static_assert(grid_t::Dim == Dim,
+                      "Grid dimension does not meet dimension of grid reader");
+
+        // The compiler will instantiate this function for all possible types of
+        // grids: Only proceed, if the grid type is known by the detector
+        if constexpr (detector_t::sf_finders::template is_defined<grid_t>()) {
+
+            // Decorate the current volume builder with the grid
+            using grid_builder_t =
+                grid_builder<detector_t, grid_t, detail::fill_by_pos>;
+
+            const auto volume_idx{
+                base_type::deserialize(grid_data.volume_link)};
+            auto v_builder =
+                det_builder.template decorate<grid_builder_t>(volume_idx);
+            auto vgr_builder = dynamic_cast<grid_builder_t *>(v_builder);
+
+            // Initialize the grid axes
+            std::vector<std::size_t> n_bins{};
+            std::vector<scalar_t> spans{};
+            std::vector<std::vector<scalar_t>> ax_bin_edges{};
+
+            for (const auto &axis_data : grid_data.axes) {
+                n_bins.push_back(axis_data.bins);
+                std::vector<scalar_t> edges{};
+                std::copy(axis_data.edges.begin(), axis_data.edges.end(),
+                          std::back_inserter(edges));
+                ax_bin_edges.emplace_back(std::move(edges));
+                spans.push_back(static_cast<scalar_t>(axis_data.edges.front()));
+                spans.push_back(static_cast<scalar_t>(axis_data.edges.back()));
+            }
+
+            vgr_builder->init_grid(spans, n_bins, ax_bin_edges);
+
+            value_t empty_sf{};
+            n_axis::multi_bin<Dim> mbin;
+            for (const auto &bin_data : grid_data.bins) {
+
+                assert(Dim == bin_data.loc_index.size() &&
+                       "Numer of local bin indices in input file does not "
+                       "match grid dimension");
+
+                // The local bin indices for the bin to be filled
+                for (const auto &[i, bin_idx] :
+                     detray::views::enumerate(bin_data.loc_index)) {
+                    mbin[i] = bin_idx;
+                }
+
+                // For now assume surfaces ids as the only grid input
+                for (const auto c : bin_data.content) {
+                    empty_sf.set_volume(volume_idx);
+                    empty_sf.set_index(static_cast<dindex>(c));
+                    vgr_builder->get().populate(mbin, empty_sf);
+                }
+            }
+        }
+    }
+};
+
+}  // namespace detray
diff --git a/io/include/detray/io/common/grid_writer.hpp b/io/include/detray/io/common/grid_writer.hpp
index 46d5fc02e..33adf559a 100644
--- a/io/include/detray/io/common/grid_writer.hpp
+++ b/io/include/detray/io/common/grid_writer.hpp
@@ -26,14 +26,21 @@
 namespace detray {
 
 /// @brief Abstract base class for accelerator grid writers
-template <class detector_t>
+template <class detector_t,
+          typename value_t = typename detector_t::surface_type>
 class grid_writer : public writer_interface<detector_t> {
 
     using base_type = writer_interface<detector_t>;
 
+    // Bin content type to be written into the grid bin payloads
+    // For detector surface descriptors, write only the surface index
+    using content_t = std::conditional_t<
+        std::is_same_v<value_t, typename detector_t::surface_type>, std::size_t,
+        value_t>;
+
     protected:
-    /// Tag the writer as "grids"
-    inline static const std::string tag = "grids";
+    /// Tag the writer as "surface_grids"
+    inline static const std::string tag = "surface_grids";
 
     public:
     /// Same constructors for this class as for base_type
@@ -43,6 +50,7 @@ class grid_writer : public writer_interface<detector_t> {
     /// Serialize the header information into its payload
     static grid_header_payload write_header(const detector_t& det,
                                             const std::string_view det_name) {
+
         grid_header_payload header_data;
 
         header_data.common = base_type::serialize(det_name, tag);
@@ -56,10 +64,14 @@ class grid_writer : public writer_interface<detector_t> {
 
     /// Serialize the grid collections of a detector @param det into their io
     /// payload
-    static detector_grids_payload serialize(
+    template <
+        typename V = value_t,
+        std::enable_if_t<std::is_same_v<V, typename detector_t::surface_type>,
+                         bool> = true>
+    static detector_grids_payload<std::size_t> serialize(
         const detector_t& det, const typename detector_t::name_map&) {
 
-        detector_grids_payload grids_data;
+        detector_grids_payload<std::size_t> grids_data;
 
         for (const auto& vol_desc : det.volumes()) {
             // Links to all acceleration data structures in the volume
@@ -84,10 +96,12 @@ class grid_writer : public writer_interface<detector_t> {
     /// Serialize a grid @param gr of type @param type and index @param idx
     /// into its io payload
     template <class grid_t>
-    static grid_payload serialize(std::size_t volume_index,
-                                  io::detail::acc_type type,
-                                  const std::size_t idx, const grid_t& gr) {
-        grid_payload grid_data;
+    static grid_payload<content_t> serialize(std::size_t volume_index,
+                                             io::detail::acc_type type,
+                                             const std::size_t idx,
+                                             const grid_t& gr) {
+
+        grid_payload<content_t> grid_data;
 
         grid_data.volume_link = base_type::serialize(volume_index);
         grid_data.acc_link = base_type::serialize(type, idx);
@@ -147,9 +161,10 @@ class grid_writer : public writer_interface<detector_t> {
 
     /// Serialize a multi-bin @param mbin into its io payload
     template <std::size_t DIM, typename content_range_t>
-    static grid_bin_payload serialize(const n_axis::multi_bin<DIM> mbin,
-                                      const content_range_t& content) {
-        grid_bin_payload bin_data;
+    static grid_bin_payload<content_t> serialize(
+        const n_axis::multi_bin<DIM> mbin, const content_range_t& content) {
+
+        grid_bin_payload<content_t> bin_data;
 
         // Local bin indices are written in the order the grid axis are stored
         for (unsigned int i = 0u; i < DIM; ++i) {
@@ -174,7 +189,9 @@ class grid_writer : public writer_interface<detector_t> {
             [[maybe_unused]] const grid_group_t& coll,
             [[maybe_unused]] const index_t& index,
             [[maybe_unused]] std::size_t volume_index,
-            [[maybe_unused]] detector_grids_payload& grids_data) const {
+            [[maybe_unused]] detector_grids_payload<content_t>& grids_data)
+            const {
+
             using accel_t = typename grid_group_t::value_type;
 
             if constexpr (detail::is_grid_v<accel_t>) {
diff --git a/io/include/detray/io/common/homogeneous_material_reader.hpp b/io/include/detray/io/common/homogeneous_material_reader.hpp
index 3a9c8c34c..5e7d8fd81 100644
--- a/io/include/detray/io/common/homogeneous_material_reader.hpp
+++ b/io/include/detray/io/common/homogeneous_material_reader.hpp
@@ -57,8 +57,9 @@ class homogeneous_material_reader : public reader_interface<detector_t> {
         // Deserialize the material volume by volume
         for (const auto& mv_data : det_mat_data.volumes) {
             // Decorate the current volume builder with material
-            auto vm_builder = det_builder.template decorate<material_builder>(
-                base_type::deserialize(mv_data.volume_link));
+            auto vm_builder =
+                det_builder.template decorate<material_builder<detector_t>>(
+                    base_type::deserialize(mv_data.volume_link));
 
             // Add the material data to the factory
             auto mat_factory = std::make_shared<material_factory<detector_t>>();
diff --git a/io/include/detray/io/common/payloads.hpp b/io/include/detray/io/common/payloads.hpp
index 78ac52529..808bffe9d 100644
--- a/io/include/detray/io/common/payloads.hpp
+++ b/io/include/detray/io/common/payloads.hpp
@@ -85,6 +85,7 @@ struct mask_payload {
 
 /// @brief  A payload for surfaces
 struct surface_payload {
+    std::optional<std::size_t> index_in_coll;
     transform_payload transform{};
     mask_payload mask{};
     std::optional<material_link_payload> material;
@@ -171,12 +172,14 @@ struct axis_payload {
 };
 
 /// @brief A payload for a grid bin
+template <typename content_t = std::size_t>
 struct grid_bin_payload {
     std::vector<unsigned int> loc_index{};
-    std::vector<std::size_t> content{};
+    std::vector<content_t> content{};
 };
 
 /// @brief A payload for a grid definition
+template <typename bin_content_t = std::size_t>
 struct grid_payload {
     using grid_type = io::detail::acc_type;
 
@@ -184,13 +187,14 @@ struct grid_payload {
     acc_links_payload acc_link{};
 
     std::vector<axis_payload> axes{};
-    std::vector<grid_bin_payload> bins{};
+    std::vector<grid_bin_payload<bin_content_t>> bins{};
     std::optional<transform_payload> transform;
 };
 
 /// @brief A payload for the grid collections of a detector
+template <typename bin_content_t = std::size_t>
 struct detector_grids_payload {
-    std::vector<grid_payload> grids = {};
+    std::vector<grid_payload<bin_content_t>> grids = {};
 };
 
 /// @}
@@ -198,7 +202,7 @@ struct detector_grids_payload {
 /// @brief A payload for a detector geometry
 struct detector_payload {
     std::vector<volume_payload> volumes = {};
-    grid_payload volume_grid;
+    grid_payload<std::size_t> volume_grid;
 };
 
 }  // namespace detray
diff --git a/io/include/detray/io/json/json_geometry_io.hpp b/io/include/detray/io/json/json_geometry_io.hpp
index 791dbec63..db11faaf4 100644
--- a/io/include/detray/io/json/json_geometry_io.hpp
+++ b/io/include/detray/io/json/json_geometry_io.hpp
@@ -67,6 +67,9 @@ inline void to_json(nlohmann::ordered_json& j, const surface_payload& s) {
     if (s.material.has_value()) {
         j["material"] = s.material.value();
     }
+    if (s.index_in_coll.has_value()) {
+        j["index_in_coll"] = s.index_in_coll.value();
+    }
 }
 
 inline void from_json(const nlohmann::ordered_json& j, surface_payload& s) {
@@ -78,6 +81,9 @@ inline void from_json(const nlohmann::ordered_json& j, surface_payload& s) {
     if (j.find("material") != j.end()) {
         s.material = j["material"];
     }
+    if (j.find("index_in_coll") != j.end()) {
+        s.index_in_coll = j["index_in_coll"];
+    }
 }
 
 inline void to_json(nlohmann::ordered_json& j, const volume_payload& v) {
diff --git a/io/include/detray/io/json/json_grids_io.hpp b/io/include/detray/io/json/json_grids_io.hpp
index 9ac533de0..a38e844c2 100644
--- a/io/include/detray/io/json/json_grids_io.hpp
+++ b/io/include/detray/io/json/json_grids_io.hpp
@@ -56,17 +56,23 @@ inline void from_json(const nlohmann::ordered_json& j, axis_payload& a) {
     a.edges = j["edges"].get<std::vector<real_io>>();
 }
 
-inline void to_json(nlohmann::ordered_json& j, const grid_bin_payload& g) {
+template <typename content_t>
+inline void to_json(nlohmann::ordered_json& j,
+                    const grid_bin_payload<content_t>& g) {
     j["loc_index"] = g.loc_index;
     j["content"] = g.content;
 }
 
-inline void from_json(const nlohmann::ordered_json& j, grid_bin_payload& g) {
+template <typename content_t>
+inline void from_json(const nlohmann::ordered_json& j,
+                      grid_bin_payload<content_t>& g) {
     g.loc_index = j["loc_index"].get<std::vector<unsigned int>>();
-    g.content = j["content"].get<std::vector<std::size_t>>();
+    g.content = j["content"].get<std::vector<content_t>>();
 }
 
-inline void to_json(nlohmann::ordered_json& j, const grid_payload& g) {
+template <typename content_t>
+inline void to_json(nlohmann::ordered_json& j,
+                    const grid_payload<content_t>& g) {
     j["volume_link"] = g.volume_link;
     j["acc_link"] = g.acc_link;
 
@@ -87,7 +93,9 @@ inline void to_json(nlohmann::ordered_json& j, const grid_payload& g) {
     }
 }
 
-inline void from_json(const nlohmann::ordered_json& j, grid_payload& g) {
+template <typename content_t>
+inline void from_json(const nlohmann::ordered_json& j,
+                      grid_payload<content_t>& g) {
     g.volume_link = j["volume_link"];
     g.acc_link = j["acc_link"];
 
@@ -99,7 +107,7 @@ inline void from_json(const nlohmann::ordered_json& j, grid_payload& g) {
 
     nlohmann::ordered_json jbins = j["bins"];
     for (auto jbin : jbins) {
-        grid_bin_payload b = jbin;
+        grid_bin_payload<content_t> b = jbin;
         g.bins.push_back(std::move(b));
     }
 
@@ -109,8 +117,9 @@ inline void from_json(const nlohmann::ordered_json& j, grid_payload& g) {
     }
 }
 
+template <typename content_t>
 inline void to_json(nlohmann::ordered_json& j,
-                    const detector_grids_payload& d) {
+                    const detector_grids_payload<content_t>& d) {
     if (not d.grids.empty()) {
         nlohmann::ordered_json jgrids;
         for (const auto& gr : d.grids) {
@@ -120,11 +129,12 @@ inline void to_json(nlohmann::ordered_json& j,
     }
 }
 
+template <typename content_t>
 inline void from_json(const nlohmann::ordered_json& j,
-                      detector_grids_payload& d) {
+                      detector_grids_payload<content_t>& d) {
     if (j.find("grids") != j.end()) {
         for (auto jgrid : j["grids"]) {
-            grid_payload grp = jgrid;
+            grid_payload<content_t> grp = jgrid;
             d.grids.push_back(std::move(grp));
         }
     }
diff --git a/io/include/detray/io/json/json_reader.hpp b/io/include/detray/io/json/json_reader.hpp
index 832c4864a..55b548b85 100644
--- a/io/include/detray/io/json/json_reader.hpp
+++ b/io/include/detray/io/json/json_reader.hpp
@@ -11,6 +11,7 @@
 #include "detray/io/common/detail/file_handle.hpp"
 #include "detray/io/common/detail/utils.hpp"
 #include "detray/io/common/geometry_reader.hpp"
+#include "detray/io/common/grid_reader.hpp"
 #include "detray/io/common/homogeneous_material_reader.hpp"
 #include "detray/io/json/json.hpp"
 #include "detray/io/json/json_serializers.hpp"
@@ -56,10 +57,11 @@ inline common_header_payload read_json_header(const std::string& file_name) {
 ///
 /// @note The resulting reader types will fulfill @c reader_interface through
 /// the common readers they are being extended with
-template <class detector_t, template <class> class common_reader_t>
-class json_reader final : public common_reader_t<detector_t> {
+template <class detector_t, template <typename...> class common_reader_t,
+          typename... Args>
+class json_reader final : public common_reader_t<detector_t, Args...> {
 
-    using base_reader = common_reader_t<detector_t>;
+    using base_reader = common_reader_t<detector_t, Args...>;
 
     public:
     /// Set json file extension
@@ -93,4 +95,12 @@ template <typename detector_t>
 using json_homogeneous_material_reader =
     json_reader<detector_t, homogeneous_material_reader>;
 
+/// Reads a homogeneous material descritption from file in json format
+template <typename detector_t,
+          typename value_t = typename detector_t::surface_type,
+          typename CAP = std::integral_constant<std::size_t, 2>,
+          typename DIM = std::integral_constant<std::size_t, 2>>
+using json_grid_reader =
+    json_reader<detector_t, grid_reader, value_t, CAP, DIM>;
+
 }  // namespace detray
diff --git a/io/include/detray/io/json/json_writer.hpp b/io/include/detray/io/json/json_writer.hpp
index e1dd56e86..5345d307f 100644
--- a/io/include/detray/io/json/json_writer.hpp
+++ b/io/include/detray/io/json/json_writer.hpp
@@ -32,10 +32,11 @@ namespace detray {
 ///
 /// @note The resulting writer types will fulfill @c writer_interface through
 /// the common writers they are being extended with
-template <class detector_t, template <class> class common_writer_t>
-class json_writer final : public common_writer_t<detector_t> {
+template <class detector_t, template <typename...> class common_writer_t,
+          typename... Args>
+class json_writer final : public common_writer_t<detector_t, Args...> {
 
-    using base_writer = common_writer_t<detector_t>;
+    using base_writer = common_writer_t<detector_t, Args...>;
 
     public:
     /// File gets created with the json file extension
@@ -88,7 +89,8 @@ using json_homogeneous_material_writer =
     json_writer<detector_t, homogeneous_material_writer>;
 
 /// Write the detector grid collections to file in json format
-template <typename detector_t>
-using json_grid_writer = json_writer<detector_t, grid_writer>;
+template <typename detector_t,
+          typename value_t = typename detector_t::surface_type>
+using json_grid_writer = json_writer<detector_t, grid_writer, value_t>;
 
 }  // namespace detray
diff --git a/plugins/svgtools/include/detray/plugins/svgtools/conversion/surface.hpp b/plugins/svgtools/include/detray/plugins/svgtools/conversion/surface.hpp
index 29686030a..bd7e52e7a 100644
--- a/plugins/svgtools/include/detray/plugins/svgtools/conversion/surface.hpp
+++ b/plugins/svgtools/include/detray/plugins/svgtools/conversion/surface.hpp
@@ -17,7 +17,6 @@
 
 // System include(s)
 #include <algorithm>
-#include <iostream>
 #include <iterator>
 
 namespace detray::svgtools::conversion {
@@ -193,4 +192,4 @@ auto surface(const typename detector_t::geometry_context& context,
         d_surface.transform(context));
 }
 
-}  // namespace detray::svgtools::conversion
\ No newline at end of file
+}  // namespace detray::svgtools::conversion
diff --git a/plugins/svgtools/include/detray/plugins/svgtools/utils/mask_utils.hpp b/plugins/svgtools/include/detray/plugins/svgtools/utils/mask_utils.hpp
index 489ed7602..b6c9b47df 100644
--- a/plugins/svgtools/include/detray/plugins/svgtools/utils/mask_utils.hpp
+++ b/plugins/svgtools/include/detray/plugins/svgtools/utils/mask_utils.hpp
@@ -13,7 +13,6 @@
 
 // System include(s)
 #include <cmath>
-#include <iostream>
 
 namespace detray::svgtools::utils {
 
@@ -41,4 +40,4 @@ inline auto global_vertices(const transform_t& transform, const mask_t& mask) {
     return ret;
 }
 
-}  // namespace detray::svgtools::utils
\ No newline at end of file
+}  // namespace detray::svgtools::utils
diff --git a/tests/common/include/tests/common/test_toy_detector.hpp b/tests/common/include/tests/common/test_toy_detector.hpp
index 75e90905f..790d7b643 100644
--- a/tests/common/include/tests/common/test_toy_detector.hpp
+++ b/tests/common/include/tests/common/test_toy_detector.hpp
@@ -37,7 +37,8 @@ inline void test_finder(const acc_t& finder, const dindex volume_index,
 
     std::vector<dindex> indices;
 
-    // Check if right volume is linked to surface in grid
+    // Check if the correct volume is linked to surface in grid
+    // and record the surface indices
     for (const auto& sf : finder.all()) {
         EXPECT_EQ(sf.volume(), volume_index);
         indices.push_back(sf.index());
@@ -51,13 +52,12 @@ inline void test_finder(const acc_t& finder, const dindex volume_index,
     }
 }
 
+template <typename bfield_t>
 inline bool test_toy_detector(
-    const detector<toy_metadata, covfie::field<bfield::const_bknd_t>>& toy_det,
-    const detector<toy_metadata, covfie::field<bfield::const_bknd_t>>::name_map&
-        names) {
+    const detector<toy_metadata, bfield_t>& toy_det,
+    const typename detector<toy_metadata, bfield_t>::name_map& names) {
 
-    using detector_t =
-        detector<toy_metadata, covfie::field<bfield::const_bknd_t>>;
+    using detector_t = detector<toy_metadata, bfield_t>;
     using geo_obj_ids = typename detector_t::geo_obj_ids;
     using volume_t = typename detector_t::volume_type;
     using nav_link_t = typename detector_t::surface_type::navigation_link;
@@ -93,25 +93,27 @@ inline bool test_toy_detector(
     // Link to outer world (leaving detector)
     constexpr auto leaving_world{detail::invalid_value<nav_link_t>()};
     const bool has_grids =
-        (sf_finders.size<sf_finder_ids::e_cylinder2_grid>() != 0u) ||
-        (sf_finders.size<sf_finder_ids::e_disc_grid>() != 0u);
-    const bool has_material = (materials.size<material_ids::e_slab>() != 0);
+        (sf_finders.template size<sf_finder_ids::e_cylinder2_grid>() != 0u) ||
+        (sf_finders.template size<sf_finder_ids::e_disc_grid>() != 0u);
+    const bool has_material =
+        (materials.template size<material_ids::e_slab>() != 0);
 
     // Check number of geomtery objects
     EXPECT_EQ(volumes.size(), 20u);
     EXPECT_EQ(toy_det.n_surfaces(), 3244u);
     EXPECT_EQ(transforms.size(ctx), 3264u);
-    EXPECT_EQ(masks.size<mask_ids::e_rectangle2>(), 2492u);
-    EXPECT_EQ(masks.size<mask_ids::e_trapezoid2>(), 648u);
-    EXPECT_EQ(masks.size<mask_ids::e_portal_cylinder2>(), 52u);
-    EXPECT_EQ(masks.size<mask_ids::e_portal_ring2>(), 52u);
-    EXPECT_EQ(sf_finders.size<sf_finder_ids::e_brute_force>(), 20u);
+    EXPECT_EQ(masks.template size<mask_ids::e_rectangle2>(), 2492u);
+    EXPECT_EQ(masks.template size<mask_ids::e_trapezoid2>(), 648u);
+    EXPECT_EQ(masks.template size<mask_ids::e_portal_cylinder2>(), 52u);
+    EXPECT_EQ(masks.template size<mask_ids::e_portal_ring2>(), 52u);
+    EXPECT_EQ(sf_finders.template size<sf_finder_ids::e_brute_force>(), 20u);
     if (has_grids) {
-        EXPECT_EQ(sf_finders.size<sf_finder_ids::e_cylinder2_grid>(), 4);
-        EXPECT_EQ(sf_finders.size<sf_finder_ids::e_disc_grid>(), 6);
+        EXPECT_EQ(sf_finders.template size<sf_finder_ids::e_cylinder2_grid>(),
+                  4);
+        EXPECT_EQ(sf_finders.template size<sf_finder_ids::e_disc_grid>(), 6);
     }
     if (has_material) {
-        EXPECT_EQ(materials.size<material_ids::e_slab>(), 3244u);
+        EXPECT_EQ(materials.template size<material_ids::e_slab>(), 3244u);
     }
 
     /** Test the links of a volume.
@@ -164,7 +166,7 @@ inline bool test_toy_detector(
                 if (has_material) {
                     EXPECT_EQ(sf_itr->material(), material_index);
                     EXPECT_EQ(
-                        materials.get<
+                        materials.template get<
                             material_ids::e_slab>()[sf_itr->material().index()],
                         mat);
                 }
@@ -206,7 +208,7 @@ inline bool test_toy_detector(
                 if (has_material) {
                     EXPECT_EQ(sf_itr->material(), material_index);
                     EXPECT_EQ(
-                        materials.get<
+                        materials.template get<
                             material_ids::e_slab>()[sf_itr->material().index()],
                         mat);
                 }
@@ -223,7 +225,7 @@ inline bool test_toy_detector(
     /// @param vol_itr iterator over the volume descriptors
     /// @param sf_finder_store the detectors acceleration data structure store
     /// @param pt_range index range of the portals in the surface lookup
-    /// @param sf_range index range of the portals in the surface lookup
+    /// @param sf_range index range of the surfaces in the surface lookup
     auto test_sf_finders =
         [has_grids](
             decltype(volumes.begin())& vol_itr,
@@ -234,9 +236,8 @@ inline bool test_toy_detector(
             const auto& link = vol_itr->full_link();
 
             // Test the portal search
-            const auto& bf_finder =
-                sf_finder_store
-                    .get<sf_finder_ids::e_brute_force>()[link[0].index()];
+            const auto& bf_finder = sf_finder_store.template get<
+                sf_finder_ids::e_brute_force>()[link[0].index()];
 
             // This means no grids, all surfaces are in the brute force method
             if (not has_grids) {
@@ -249,11 +250,11 @@ inline bool test_toy_detector(
                 // Test the module search if grids were filled
                 if (not link[1].is_invalid()) {
                     if (link[1].id() == sf_finder_ids::e_cylinder2_grid) {
-                        const auto& cyl_grid = sf_finder_store.get<
+                        const auto& cyl_grid = sf_finder_store.template get<
                             sf_finder_ids::e_cylinder2_grid>()[link[1].index()];
                         test_finder(cyl_grid, vol_itr->index(), sf_range);
                     } else {
-                        const auto& disc_grid = sf_finder_store.get<
+                        const auto& disc_grid = sf_finder_store.template get<
                             sf_finder_ids::e_disc_grid>()[link[1].index()];
                         test_finder(disc_grid, vol_itr->index(), sf_range);
                     }
diff --git a/tests/common/include/tests/common/tools/hash_tree.hpp b/tests/common/include/tests/common/tools/hash_tree.hpp
index 439f2ac7b..577ee0bc4 100644
--- a/tests/common/include/tests/common/tools/hash_tree.hpp
+++ b/tests/common/include/tests/common/tools/hash_tree.hpp
@@ -8,7 +8,6 @@
 #pragma once
 
 #include <cmath>
-#include <iostream>
 #include <iterator>
 #include <type_traits>
 
diff --git a/tests/unit_tests/core/CMakeLists.txt b/tests/unit_tests/core/CMakeLists.txt
index 4e5041642..9c91991d4 100644
--- a/tests/unit_tests/core/CMakeLists.txt
+++ b/tests/unit_tests/core/CMakeLists.txt
@@ -18,4 +18,5 @@ detray_add_test( core
    "utils_ranges.cpp"
    "utils_quadratic_equation.cpp"
    "utils_sort.cpp"
+   "utils_type_list.cpp"
    LINK_LIBRARIES GTest::gtest_main detray_tests_common detray::core_array )
diff --git a/tests/unit_tests/core/tools_hash_tree.cpp b/tests/unit_tests/core/tools_hash_tree.cpp
index fc0982c84..097980d21 100644
--- a/tests/unit_tests/core/tools_hash_tree.cpp
+++ b/tests/unit_tests/core/tools_hash_tree.cpp
@@ -8,7 +8,6 @@
 #include <gtest/gtest.h>
 
 #include <cmath>
-#include <iostream>
 
 // detray test
 #include "tests/common/tools/hash_tree.hpp"
diff --git a/tests/unit_tests/core/utils_type_list.cpp b/tests/unit_tests/core/utils_type_list.cpp
new file mode 100644
index 000000000..4f41dbe21
--- /dev/null
+++ b/tests/unit_tests/core/utils_type_list.cpp
@@ -0,0 +1,40 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2023 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s).
+#include "detray/utils/type_list.hpp"
+
+// Google Test include(s).
+#include <gtest/gtest.h>
+
+// Test type list implementation
+TEST(utils, type_list) {
+    using namespace detray;
+
+    using list = types::type_list<float, int, double>;
+
+    static_assert(std::is_same_v<types::front<list>, float>,
+                  "Could not access type list front");
+
+    static_assert(std::is_same_v<types::back<list>, double>,
+                  "Could not access type list back");
+
+    static_assert(std::is_same_v<types::push_back<list, char>,
+                                 types::type_list<float, int, double, char>>,
+                  "Failed to push back new type");
+
+    static_assert(std::is_same_v<types::push_front<list, char>,
+                                 types::type_list<char, float, int, double>>,
+                  "Failed to push front new type");
+
+    static_assert(types::size<list> == 3ul, "Incorrect size");
+
+    static_assert(std::is_same_v<types::at<list, 1>, int>,
+                  "Failed access type");
+
+    types::print<list>();
+}
diff --git a/tests/unit_tests/cpu/grid_grid_builder.cpp b/tests/unit_tests/cpu/grid_grid_builder.cpp
index 1156d00d4..a17470a89 100644
--- a/tests/unit_tests/cpu/grid_grid_builder.cpp
+++ b/tests/unit_tests/cpu/grid_grid_builder.cpp
@@ -118,16 +118,6 @@ GTEST_TEST(detray_tools, grid_factory) {
                                                 closed<label::e_cyl_z>,
                                                 regular<>, irregular<>>(
         cyl2, {bin_edges_z.size() - 1, 10u}, {bin_edges_phi, bin_edges_z});
-
-    // Get grid collection with the correct allocator
-    auto grid_coll = gr_factory.new_collection<decltype(cyl_gr)>();
-    EXPECT_TRUE(grid_coll.size() == 0u);
-    grid_coll.push_back(cyl_gr);
-    grid_coll.push_back(cyl_gr2);
-    EXPECT_TRUE(grid_coll.size() == 2u);
-
-    EXPECT_EQ(grid_coll[0].search(loc_p)[0], 33u);
-    // gr_factory.to_string(grid_coll[0]);
 }
 
 /// Unittest: Test the grid builder
@@ -137,10 +127,10 @@ GTEST_TEST(detray_tools, grid_builder) {
     using cyl_grid_t =
         grid<coordinate_axes<cylinder2D<>::axes<>, false, host_container_types>,
              test_detector_t::surface_type, simple_serializer,
-             regular_attacher<9>>;
+             regular_attacher<1>>;
 
     auto gbuilder = grid_builder<test_detector_t, cyl_grid_t,
-                                 detray::detail::bin_associator>{};
+                                 detray::detail::bin_associator>{nullptr};
 
     // The cylinder portals are at the end of the surface range by construction
     const auto cyl_mask = mask<cylinder2D<>>{0u, 10.f, -500.f, 500.f};
@@ -173,7 +163,7 @@ GTEST_TEST(detray_tools, decorator_grid_builder) {
     using cyl_grid_t =
         grid<coordinate_axes<cylinder2D<>::axes<>, false, host_container_types>,
              test_detector_t::surface_type, simple_serializer,
-             regular_attacher<9>>;
+             regular_attacher<1>>;
 
     using portal_cylinder_factory_t =
         surface_factory<test_detector_t, cylinder2D<>,
@@ -225,18 +215,18 @@ GTEST_TEST(detray_tools, decorator_grid_builder) {
     auto rect_factory = std::make_shared<rectangle_factory>();
 
     typename rectangle_factory::sf_data_collection rect_sf_data;
-    rect_sf_data.emplace_back(transform3(point3{0.f, 0.f, -10.f}), vol_idx,
+    rect_sf_data.emplace_back(transform3(point3{7.07f, 7.07f, -500.f}), vol_idx,
                               std::vector<scalar>{10.f, 8.f});
-    rect_sf_data.emplace_back(transform3(point3{0.f, 0.f, -20.f}), vol_idx,
+    rect_sf_data.emplace_back(transform3(point3{7.07f, 7.07f, -250.f}), vol_idx,
                               std::vector<scalar>{10.f, 8.f});
-    rect_sf_data.emplace_back(transform3(point3{0.f, 0.f, -30.f}), vol_idx,
+    rect_sf_data.emplace_back(transform3(point3{7.07f, 7.07f, 100.f}), vol_idx,
                               std::vector<scalar>{10.f, 8.f});
     rect_factory->push_back(std::move(rect_sf_data));
 
     auto trpz_factory = std::make_shared<trapezoid_factory>();
 
     typename trapezoid_factory::sf_data_collection trpz_sf_data;
-    trpz_sf_data.emplace_back(transform3(point3{0.f, 0.f, 1000.f}), vol_idx,
+    trpz_sf_data.emplace_back(transform3(point3{7.07f, 7.07f, 600.f}), vol_idx,
                               std::vector<scalar>{1.f, 3.f, 2.f, 0.25f});
     trpz_factory->push_back(std::move(trpz_sf_data));
 
@@ -280,8 +270,7 @@ GTEST_TEST(detray_tools, decorator_grid_builder) {
     EXPECT_EQ(vol.template link<geo_obj_id::e_passive>(),
               sf_range[geo_obj_id::e_passive]);
 
-    // Only the portals should be in the detector's surface container now
-    EXPECT_EQ(d.surface_lookup().size(), 3u);
+    EXPECT_EQ(d.surface_lookup().size(), 7u);
     EXPECT_EQ(d.mask_store().template size<mask_id::e_portal_cylinder2>(), 3u);
     EXPECT_EQ(d.mask_store().template size<mask_id::e_portal_ring2>(), 0u);
     EXPECT_EQ(d.mask_store().template size<mask_id::e_cylinder2>(), 3u);
@@ -289,9 +278,13 @@ GTEST_TEST(detray_tools, decorator_grid_builder) {
     EXPECT_EQ(d.mask_store().template size<mask_id::e_trapezoid2>(), 1u);
 
     // check the portals in the detector
-    for (const auto& sf : d.surface_lookup()) {
+    const auto& bf_finder =
+        d.surface_store()
+            .template get<test_detector_t::sf_finders::id::e_brute_force>()[0];
+    for (const auto& sf : bf_finder.all()) {
         EXPECT_TRUE((sf.id() == surface_id::e_portal) or
                     (sf.id() == surface_id::e_passive));
+        EXPECT_EQ(sf.volume(), 0u);
     }
 
     // check the sensitive surfaces in the grid
@@ -299,7 +292,7 @@ GTEST_TEST(detray_tools, decorator_grid_builder) {
         d.surface_store()
             .template get<
                 test_detector_t::sf_finders::id::e_cylinder2_grid>()[0];
-    dindex trf_idx{4u};
+    dindex trf_idx{3u};
     for (const auto& sf : cyl_grid.all()) {
         EXPECT_TRUE(sf.is_sensitive());
         EXPECT_EQ(sf.volume(), 0u);
diff --git a/tests/unit_tests/cpu/tools_material_builder.cpp b/tests/unit_tests/cpu/tools_material_builder.cpp
index 3f6afd257..91e59334d 100644
--- a/tests/unit_tests/cpu/tools_material_builder.cpp
+++ b/tests/unit_tests/cpu/tools_material_builder.cpp
@@ -239,7 +239,8 @@ GTEST_TEST(detray_tools, detector_builder_with_material) {
 
     // Add material
     auto mv_builder =
-        det_builder.template decorate<material_builder>(vbuilder->vol_index());
+        det_builder.template decorate<material_builder<detector_t>>(
+            vbuilder->vol_index());
 
     typename detector_t::point3 t{0.f, 0.f, 20.f};
     mv_builder->add_volume_placement(t);
diff --git a/tests/unit_tests/io/io_json_detector_reader.cpp b/tests/unit_tests/io/io_json_detector_reader.cpp
index 0e816149f..60f6a3ec0 100644
--- a/tests/unit_tests/io/io_json_detector_reader.cpp
+++ b/tests/unit_tests/io/io_json_detector_reader.cpp
@@ -98,12 +98,17 @@ TEST(io, json_toy_detector_reader) {
     io::detector_reader_config reader_cfg{};
     reader_cfg.add_file("toy_detector_geometry.json")
         .add_file("toy_detector_homogeneous_material.json")
+        .add_file("toy_detector_surface_grids.json")
         .add_file(toy_cfg.bfield_file());
 
     const auto [det, names] =
-        io::read_detector<detector_t>(host_mr, reader_cfg);
+        io::read_detector<detector_t, 1u>(host_mr, reader_cfg);
+
+    // Write the result to a different set of files
+    writer_cfg.replace_files(false);
+    io::write_detector(det, names, writer_cfg);
 
-    // EXPECT_TRUE(test_toy_detector(det, names));
+    EXPECT_TRUE(test_toy_detector(det, names));
 }
 
 /// Test the reading and writing of a wire chamber
@@ -124,7 +129,8 @@ TEST(io, json_wire_chamber_reader) {
     // Read the detector back in
     io::detector_reader_config reader_cfg{};
     reader_cfg.add_file("wire_chamber_geometry.json")
-        .add_file("wire_chamber_homogeneous_material.json");
+        .add_file("wire_chamber_homogeneous_material.json")
+        .add_file("wire_chamber_surface_grids.json");
 
     const auto [det, names] =
         io::read_detector<detector_t>(host_mr, reader_cfg);
diff --git a/tests/unit_tests/io/io_json_payload.cpp b/tests/unit_tests/io/io_json_payload.cpp
index f1fc95d89..11bf29770 100644
--- a/tests/unit_tests/io/io_json_payload.cpp
+++ b/tests/unit_tests/io/io_json_payload.cpp
@@ -107,14 +107,14 @@ TEST(io, json_axis_payload) {
 /// This tests the json io for a grid bin
 TEST(io, json_bin_payload) {
 
-    detray::grid_bin_payload b;
+    detray::grid_bin_payload<> b;
     b.loc_index = std::vector<unsigned int>{1u, 0u, 2u};
     b.content = std::vector<std::size_t>{0u, 1u, 2u, 3u};
 
     nlohmann::ordered_json jbin;
     jbin["bin"] = b;
 
-    detray::grid_bin_payload pb = jbin["bin"];
+    detray::grid_bin_payload<> pb = jbin["bin"];
 
     EXPECT_EQ(b.loc_index.size(), pb.loc_index.size());
     EXPECT_EQ(b.content.size(), pb.content.size());
@@ -123,7 +123,7 @@ TEST(io, json_bin_payload) {
 /// This tests the json io for a grid
 TEST(io, json_grid_payload) {
 
-    std::vector<detray::grid_bin_payload> bins = {
+    std::vector<detray::grid_bin_payload<>> bins = {
         {{0u, 1u}, {0u, 2u}}, {{1u, 1u}, {1u, 2u}}, {{2u, 1u}, {2u, 2u}}};
 
     detray::axis_payload a0{
@@ -136,8 +136,8 @@ TEST(io, json_grid_payload) {
         detray::n_axis::binning::e_regular, detray::n_axis::bounds::e_closed,
         detray::n_axis::label::e_r, 2u, std::vector<detray::real_io>{0.f, 2.f}};
 
-    detray::grid_payload g;
-    g.acc_link = {detray::grid_payload::grid_type::polar2_grid, 12u};
+    detray::grid_payload<> g;
+    g.acc_link = {detray::grid_payload<>::grid_type::polar2_grid, 12u};
     g.volume_link = {2u};
     g.axes = {a0, a1};
     g.bins = bins;
@@ -145,7 +145,7 @@ TEST(io, json_grid_payload) {
     nlohmann::ordered_json j;
     j["grid"] = g;
 
-    detray::grid_payload pg = j["grid"];
+    detray::grid_payload<> pg = j["grid"];
 
     EXPECT_EQ(g.acc_link.type, pg.acc_link.type);
     EXPECT_EQ(g.acc_link.index, pg.acc_link.index);
diff --git a/utils/include/detray/detectors/create_toy_geometry.hpp b/utils/include/detray/detectors/create_toy_geometry.hpp
index da4399500..e60887adc 100644
--- a/utils/include/detray/detectors/create_toy_geometry.hpp
+++ b/utils/include/detray/detectors/create_toy_geometry.hpp
@@ -190,11 +190,13 @@ inline void create_barrel_modules(context_t &ctx, volume_type &vol,
 
 /// Helper function that creates a surface grid of rectangular barrel modules.
 ///
-/// @param surfaces_grid the grid to be created with proper axes
-/// @param resource vecmem memory resource
-/// @param cfg config struct for module creation
+/// @param ctx the detector geometry context
+/// @param resource memory resource for the grid data allocations
+/// @param vol the detector volume in which to build the grid
+/// @param det the detector to fill the grid into
+/// @param module_factory factory to create the volume surfaces
 template <
-    typename detector_t, typename config_t, typename factory_t,
+    typename detector_t, typename factory_t,
     std::enable_if_t<
         std::is_invocable_v<factory_t, typename detector_t::geometry_context &,
                             typename detector_t::volume_type &,
@@ -206,8 +208,7 @@ template <
 inline void add_cylinder_grid(const typename detector_t::geometry_context &ctx,
                               vecmem::memory_resource &resource,
                               typename detector_t::volume_type &vol,
-                              detector_t &det, const config_t &cfg,
-                              factory_t &module_factory) {
+                              detector_t &det, factory_t &module_factory) {
     // Get relevant ids
     using geo_obj_ids = typename detector_t::geo_obj_ids;
 
@@ -217,13 +218,25 @@ inline void add_cylinder_grid(const typename detector_t::geometry_context &ctx,
     using cyl_grid_t =
         typename detector_t::surface_container::template get_type<grid_id>;
     auto gbuilder =
-        grid_builder<detector_t, cyl_grid_t, detray::detail::fill_by_pos>{};
+        grid_builder<detector_t, cyl_grid_t, detray::detail::fill_by_pos>{
+            nullptr};
 
-    // The disc portals are at the end of the portal range by construction
-    auto portal_mask_idx = (det.portals(vol).end() - 3u)->mask().index();
-    const auto &cyl_mask =
+    // The portals are at the end of the portal range by construction
+    auto portal_mask_idx = (det.portals(vol).end() - 4u)->mask().index();
+    const auto &inner_cyl_mask =
+        det.mask_store().template get<cyl_id>().at(portal_mask_idx);
+    portal_mask_idx = (det.portals(vol).end() - 3u)->mask().index();
+    const auto &outer_cyl_mask =
         det.mask_store().template get<cyl_id>().at(portal_mask_idx);
 
+    // Correct cylinder radius so that the grid lies in the middle
+    using cyl_mask_t = detail::remove_cvref_t<decltype(outer_cyl_mask)>;
+    typename cyl_mask_t::mask_values mask_values{outer_cyl_mask.values()};
+    mask_values[cylinder2D<>::e_r] =
+        0.5f * (inner_cyl_mask.values()[cylinder2D<>::e_r] +
+                outer_cyl_mask.values()[cylinder2D<>::e_r]);
+    const cyl_mask_t cyl_mask{mask_values, 0u};
+
     // Create the sensitive surfaces
     typename detector_t::surface_container_t surfaces(&resource);
     typename detector_t::mask_container masks(resource);
@@ -249,16 +262,18 @@ inline void add_cylinder_grid(const typename detector_t::geometry_context &ctx,
         det.add_surface_to_lookup(sf_desc);
     }
 
-    // Add new grid to the detector
-    gbuilder.init_grid(cyl_mask, {cfg.m_binning.first, cfg.m_binning.second});
-    gbuilder.fill_grid(detector_volume{det, vol}, surfaces, transforms, masks,
-                       ctx);
-    assert(gbuilder.get().all().size() == surfaces.size());
-
     // Add transforms, masks and material to detector
     det.append_masks(std::move(masks));
     det.append_transforms(std::move(transforms));
     det.append_materials(std::move(materials));
+
+    // Add new grid to the detector
+    gbuilder.init_grid(cyl_mask, {module_factory.cfg.m_binning.first,
+                                  module_factory.cfg.m_binning.second});
+    gbuilder.fill_grid(detector_volume{det, vol}, det.surface_lookup(),
+                       det.transform_store(), det.mask_store(), ctx);
+    assert(gbuilder.get().all().size() == surfaces.size());
+
     det.surface_store().template push_back<grid_id>(gbuilder.get());
     vol.template set_link<geo_obj_ids::e_sensitive>(
         grid_id, det.surface_store().template size<grid_id>() - 1u);
@@ -266,12 +281,13 @@ inline void add_cylinder_grid(const typename detector_t::geometry_context &ctx,
 
 /// Helper function that creates a surface grid of trapezoidal endcap modules.
 ///
-/// @param vol the detector volume that should be equipped with a grid
-/// @param surfaces_grid the grid to be created with proper axes
-/// @param resource vecmem memory resource
-/// @param cfg config struct for module creation
+/// @param ctx the detector geometry context
+/// @param resource memory resource for the grid data allocations
+/// @param vol the detector volume in which to build the grid
+/// @param det the detector to fill the grid into
+/// @param module_factory factory to create the volume surfaces
 template <
-    typename detector_t, typename config_t, typename factory_t,
+    typename detector_t, typename factory_t,
     std::enable_if_t<
         std::is_invocable_v<factory_t, typename detector_t::geometry_context &,
                             typename detector_t::volume_type &,
@@ -283,8 +299,7 @@ template <
 inline void add_disc_grid(const typename detector_t::geometry_context &ctx,
                           vecmem::memory_resource &resource,
                           typename detector_t::volume_type &vol,
-                          detector_t &det, const config_t &cfg,
-                          factory_t &module_factory) {
+                          detector_t &det, factory_t &module_factory) {
     // Get relevant ids
     using geo_obj_ids = typename detector_t::geo_obj_ids;
 
@@ -294,7 +309,8 @@ inline void add_disc_grid(const typename detector_t::geometry_context &ctx,
     using disc_grid_t =
         typename detector_t::surface_container::template get_type<grid_id>;
     auto gbuilder =
-        grid_builder<detector_t, disc_grid_t, detray::detail::fill_by_pos>{};
+        grid_builder<detector_t, disc_grid_t, detray::detail::fill_by_pos>{
+            nullptr};
 
     // The disc portals are at the end of the portal range by construction
     auto portal_mask_idx = det.portals(vol).back().mask().index();
@@ -326,17 +342,18 @@ inline void add_disc_grid(const typename detector_t::geometry_context &ctx,
         det.add_surface_to_lookup(sf_desc);
     }
 
-    // Add new grid to the detector
-    gbuilder.init_grid(disc_mask,
-                       {cfg.disc_binning.size(), cfg.disc_binning.front()});
-    gbuilder.fill_grid(detector_volume{det, vol}, surfaces, transforms, masks,
-                       ctx);
-    assert(gbuilder.get().all().size() == surfaces.size());
-
     // Add transforms, masks and material to detector
     det.append_masks(std::move(masks));
     det.append_transforms(std::move(transforms));
     det.append_materials(std::move(materials));
+
+    // Add new grid to the detector
+    gbuilder.init_grid(disc_mask, {module_factory.cfg.disc_binning.size(),
+                                   module_factory.cfg.disc_binning.back()});
+    gbuilder.fill_grid(detector_volume{det, vol}, det.surface_lookup(),
+                       det.transform_store(), det.mask_store(), ctx);
+    assert(gbuilder.get().all().size() == surfaces.size());
+
     det.surface_store().template push_back<grid_id>(gbuilder.get());
     vol.template set_link<geo_obj_ids::e_sensitive>(
         grid_id, det.surface_store().template size<grid_id>() - 1u);
@@ -810,8 +827,7 @@ inline void add_endcap_detector(
             dindex vol_idx = det.volumes().back().index();
             names[vol_idx + 1u] = "endcap_" + std::to_string(vol_idx);
 
-            add_disc_grid(ctx, resource, det.volumes().back(), det, cfg,
-                          m_factory);
+            add_disc_grid(ctx, resource, det.volumes().back(), det, m_factory);
         }
     }
 }
@@ -896,7 +912,7 @@ inline void add_barrel_detector(
             dindex vol_idx = det.volumes().back().index();
             names[vol_idx + 1u] = "barrel_" + std::to_string(vol_idx);
 
-            add_cylinder_grid(ctx, resource, det.volumes().back(), det, cfg,
+            add_cylinder_grid(ctx, resource, det.volumes().back(), det,
                               m_factory);
         }
     }
diff --git a/utils/include/detray/detectors/create_wire_chamber.hpp b/utils/include/detray/detectors/create_wire_chamber.hpp
index 78b3f7b10..4be903867 100644
--- a/utils/include/detray/detectors/create_wire_chamber.hpp
+++ b/utils/include/detray/detectors/create_wire_chamber.hpp
@@ -247,6 +247,11 @@ auto create_wire_chamber(vecmem::memory_resource &resource,
             det.add_surface_to_lookup(sf_desc);
         }
 
+        // Add transforms, masks and material to detector
+        det.append_masks(std::move(masks));
+        det.append_transforms(std::move(transforms));
+        det.append_materials(std::move(materials));
+
         //
         // Fill Grid
         //
@@ -259,27 +264,34 @@ auto create_wire_chamber(vecmem::memory_resource &resource,
         using cyl_grid_t =
             typename detector_t::surface_container::template get_type<grid_id>;
         auto gbuilder =
-            grid_builder<detector_t, cyl_grid_t, detray::detail::fill_by_pos>{};
+            grid_builder<detector_t, cyl_grid_t, detray::detail::fill_by_pos>{
+                nullptr};
 
-        // The disc portals are at the end of the portal range by construction
-        auto portal_mask_idx = (det.portals(vol).end() - 3u)->mask().index();
-        const auto &cyl_mask =
+        // The portal portals are at the end of the portal range by construction
+        auto portal_mask_idx = (det.portals(vol).end() - 4u)->mask().index();
+        const auto &inner_cyl_mask =
+            det.mask_store().template get<cyl_id>().at(portal_mask_idx);
+        portal_mask_idx = (det.portals(vol).end() - 3u)->mask().index();
+        const auto &outer_cyl_mask =
             det.mask_store().template get<cyl_id>().at(portal_mask_idx);
 
+        // Correct cylinder radius so that the grid lies in the middle
+        using cyl_mask_t = detail::remove_cvref_t<decltype(outer_cyl_mask)>;
+        typename cyl_mask_t::mask_values mask_values{outer_cyl_mask.values()};
+        mask_values[cylinder2D<>::e_r] =
+            0.5f * (inner_cyl_mask.values()[cylinder2D<>::e_r] +
+                    outer_cyl_mask.values()[cylinder2D<>::e_r]);
+        const cyl_mask_t cyl_mask{mask_values, 0u};
+
         // Add new grid to the detector
         gbuilder.init_grid(cyl_mask, {100u, 1u});
-        gbuilder.fill_grid(detector_volume{det, vol}, surfaces, transforms,
-                           masks, ctx0);
+        gbuilder.fill_grid(detector_volume{det, vol}, det.surface_lookup(),
+                           det.transform_store(), det.mask_store(), ctx0);
 
         det.surface_store().template push_back<grid_id>(gbuilder.get());
         vol.template set_link<geo_obj_ids::e_sensitive>(
             grid_id, det.surface_store().template size<grid_id>() - 1u);
 
-        // Add transforms, masks and material to detector
-        det.append_masks(std::move(masks));
-        det.append_transforms(std::move(transforms));
-        det.append_materials(std::move(materials));
-
         // Add volume grid
         // TODO: Fill it
 
diff --git a/utils/include/detray/detectors/detector_helper.hpp b/utils/include/detray/detectors/detector_helper.hpp
index 393dc880f..7061ce96b 100644
--- a/utils/include/detray/detectors/detector_helper.hpp
+++ b/utils/include/detray/detectors/detector_helper.hpp
@@ -8,8 +8,6 @@
 #pragma once
 
 // Project include(s).
-#include <iostream>
-
 #include "detray/core/detail/data_context.hpp"
 #include "detray/definitions/geometry.hpp"
 #include "detray/definitions/indexing.hpp"
@@ -211,4 +209,4 @@ struct detector_helper {
     }
 };
 
-}  // namespace detray::detail
\ No newline at end of file
+}  // namespace detray::detail
diff --git a/utils/include/detray/detectors/toy_metadata.hpp b/utils/include/detray/detectors/toy_metadata.hpp
index bcbc5f5ba..8858536ca 100644
--- a/utils/include/detray/detectors/toy_metadata.hpp
+++ b/utils/include/detray/detectors/toy_metadata.hpp
@@ -47,7 +47,7 @@ struct toy_metadata {
     template <typename grid_shape_t, typename bin_entry_t, typename container_t>
     using surface_grid_t =
         grid<coordinate_axes<grid_shape_t, false, container_t>, bin_entry_t,
-             simple_serializer, regular_attacher<9>>;
+             simple_serializer, regular_attacher<1>>;
 
     // cylindrical grid for the barrel layers
     template <typename bin_entry_t, typename container_t>