From 7f498fbb0280e78b11c35eb33fc9f173666981b3 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 material grids as material maps and add the relevant
 types to the detector metadata. The material maps of different shapes are
 tested in a dedicated unittest.

Mapeterial maps are added to the toy detector and wire chamber via
the detector helper. The material is a silicon and aluminium mixture
that is mapped in configurable binning onto the disc and cylinder
portals.
---
 core/include/detray/core/detector.hpp         |  11 +-
 .../include/detray/core/detector_metadata.hpp |  63 ++++-
 core/include/detray/geometry/surface.hpp      |  11 +-
 core/include/detray/masks/masks.hpp           |   4 +
 core/include/detray/materials/material.hpp    |  57 +++-
 .../include/detray/materials/material_map.hpp |  29 +++
 .../include/detray/materials/material_rod.hpp |  10 +
 .../detray/materials/material_slab.hpp        |  12 +
 .../actors/pointwise_material_interactor.hpp  | 100 ++++---
 core/include/detray/propagator/navigator.hpp  |   8 +
 .../surface_finders/brute_force_finder.hpp    |   9 +
 .../detray/surface_finders/grid/grid.hpp      |  16 +-
 .../surface_finders/grid/grid_collection.hpp  |  38 +++
 core/include/detray/tools/volume_builder.hpp  |  13 +-
 .../detray/io/common/detail/type_traits.hpp   |  19 +-
 .../detray/io/common/detector_writer.hpp      |   2 +-
 .../io/common/homogeneous_material_writer.hpp |  32 ++-
 .../tests/common/test_toy_detector.hpp        |  24 +-
 tests/unit_tests/cpu/CMakeLists.txt           |   1 +
 tests/unit_tests/cpu/material_maps.cpp        | 237 +++++++++++++++++
 .../detray/tutorial/detector_metadata.hpp     |   7 +-
 .../detray/detectors/create_toy_geometry.hpp  | 245 ++++++++++++------
 .../detray/detectors/create_wire_chamber.hpp  |  52 +++-
 .../detray/detectors/detector_helper.hpp      | 186 ++++++++-----
 .../include/detray/detectors/itk_metadata.hpp |  44 +++-
 .../detray/detectors/telescope_metadata.hpp   |  10 +-
 .../include/detray/detectors/toy_metadata.hpp |  36 ++-
 27 files changed, 998 insertions(+), 278 deletions(-)
 create mode 100644 core/include/detray/materials/material_map.hpp
 create mode 100644 tests/unit_tests/cpu/material_maps.cpp

diff --git a/core/include/detray/core/detector.hpp b/core/include/detray/core/detector.hpp
index 2c217bf41c..9f2d8b3f90 100644
--- a/core/include/detray/core/detector.hpp
+++ b/core/include/detray/core/detector.hpp
@@ -97,7 +97,7 @@ class detector {
 
     /// Forward mask types that are present in this detector
     using material_container =
-        typename metadata::template material_store<tuple_type, vector_type>;
+        typename metadata::template material_store<tuple_type, container_t>;
     using materials = typename material_container::value_types;
     using material_link = typename material_container::single_link;
 
@@ -314,15 +314,16 @@ class detector {
     /// @returns the portals of a given volume @param v - const
     /// @note Depending on the detector type, this can also container other
     /// surfaces
-    DETRAY_HOST_DEVICE constexpr auto portals(const volume_type &v) const {
+    DETRAY_HOST_DEVICE constexpr auto portals(const volume_type &v) {
 
         // Index of the portal collection in the surface store
         const dindex pt_coll_idx{
             v.template link<geo_obj_ids::e_portal>().index()};
 
-        const auto &pt_coll =
-            _surfaces
-                .template get<sf_finders::id::e_brute_force>()[pt_coll_idx];
+        auto &bf_finder =
+            _surfaces.template get<sf_finders::id::e_brute_force>();
+
+        auto pt_coll = bf_finder[pt_coll_idx];
 
         return pt_coll.all();
     }
diff --git a/core/include/detray/core/detector_metadata.hpp b/core/include/detray/core/detector_metadata.hpp
index da5140f80a..e209e0e0d4 100644
--- a/core/include/detray/core/detector_metadata.hpp
+++ b/core/include/detray/core/detector_metadata.hpp
@@ -19,6 +19,7 @@
 #include "detray/masks/masks.hpp"
 #include "detray/masks/unbounded.hpp"
 #include "detray/masks/unmasked.hpp"
+#include "detray/materials/material_map.hpp"
 #include "detray/materials/material_rod.hpp"
 #include "detray/materials/material_slab.hpp"
 #include "detray/surface_finders/accelerator_grid.hpp"
@@ -61,9 +62,34 @@ struct default_metadata {
     using slab = material_slab<detray::scalar>;
     using rod = material_rod<detray::scalar>;
 
-    /// @TODO material grids
+    /// Material grid types (default: closed bounds, regular binning)
     /// @{
-    /// ...
+
+    // Disc material grid
+    template <typename container_t>
+    using disc_map_t =
+        material_map<ring2D<>::axes<>, detray::scalar, container_t>;
+
+    // Cylindrical material grid
+    template <typename container_t>
+    using cylinder2_map_t =
+        material_map<cylinder2D<>::axes<>, detray::scalar, container_t>;
+
+    // Rectangular material grid
+    template <typename container_t>
+    using rectangular_map_t =
+        material_map<rectangle2D<>::axes<>, detray::scalar, container_t>;
+
+    // Cuboid volume material grid
+    template <typename container_t>
+    using cuboid_map_t =
+        material_map<cuboid3D<>::axes<>, detray::scalar, container_t>;
+
+    // Cylindrical volume material grid
+    template <typename container_t>
+    using cylinder3_map_t =
+        material_map<cylinder3D::axes<>, detray::scalar, container_t>;
+
     /// @}
 
     /// surface grid types (default boundaries: closed binning)
@@ -162,19 +188,36 @@ unbounded_cell, unmasked_plane*/>;
 
     /// Give your material types a name (needs to be consecutive and has to
     /// match the types position in the mask store!)
-    /// @TODO: Add the material grid types for every surface shape
     enum class material_ids {
-        e_slab = 0,
-        e_rod = 1,
-        // ... material map types
-        e_none = 2,
+        // Material texture (grid) shapes
+        e_disc2_map = 0u,
+        e_cylinder2_map = 1u,
+        e_rectangle2_map = 2u,
+        e_trapezoid2_map = 2u,
+        e_annulus2_map = 0u,
+        e_cell_wire_map = 6u,
+        e_straw_wire_map = 6u,
+        // Volume material
+        e_cuboid3_map = 3u,
+        e_cylinder3_map = 4u,
+        // Homogeneous mapetrial
+        e_slab = 5u,
+        e_rod = 6u,
+        e_none = 7u,
     };
 
     /// How to store materials
     template <template <typename...> class tuple_t = dtuple,
-              template <typename...> class vector_t = dvector>
-    using material_store = regular_multi_store<material_ids, empty_context,
-                                               tuple_t, vector_t, slab, rod>;
+              typename container_t = host_container_types>
+    using material_store =
+        multi_store<material_ids, empty_context, tuple_t,
+                    grid_collection<disc_map_t<container_t>>,
+                    grid_collection<cylinder2_map_t<container_t>>,
+                    grid_collection<rectangular_map_t<container_t>>,
+                    grid_collection<cuboid_map_t<container_t>>,
+                    grid_collection<cylinder3_map_t<container_t>>,
+                    typename container_t::template vector_type<slab>,
+                    typename container_t::template vector_type<rod>>;
 
     /// How to link to the entries in the data stores
     using transform_link = typename transform_store<>::link_type;
diff --git a/core/include/detray/geometry/surface.hpp b/core/include/detray/geometry/surface.hpp
index a3f66cdd0f..6cee92fb0d 100644
--- a/core/include/detray/geometry/surface.hpp
+++ b/core/include/detray/geometry/surface.hpp
@@ -269,7 +269,7 @@ class surface {
     /// @tparam functor_t the prescription to be applied to the mask
     /// @tparam Args      types of additional arguments to the functor
     template <typename functor_t, typename... Args>
-    DETRAY_HOST_DEVICE constexpr auto visit_mask(Args &&... args) const {
+    DETRAY_HOST_DEVICE constexpr auto visit_mask(Args &&...args) const {
         const auto &masks = m_detector.mask_store();
 
         return masks.template visit<functor_t>(m_desc.mask(),
@@ -281,7 +281,7 @@ class surface {
     /// @tparam functor_t the prescription to be applied to the mask
     /// @tparam Args      types of additional arguments to the functor
     template <typename functor_t, typename... Args>
-    DETRAY_HOST_DEVICE constexpr auto visit_material(Args &&... args) const {
+    DETRAY_HOST_DEVICE constexpr auto visit_material(Args &&...args) const {
         const auto &materials = m_detector.material_store();
 
         return materials.template visit<functor_t>(m_desc.material(),
@@ -325,7 +325,8 @@ class surface {
         }
         // Only check, if there is material in the detector
         if (not m_detector.material_store().all_empty()) {
-            if (is_invalid_value(m_desc.material())) {
+            if (m_desc.material().id() != detector_t::materials::id::e_none and
+                m_desc.material().is_invalid_index()) {
                 os << "ERROR: Surface does not have valid material:\n"
                    << *this << std::endl;
                 return false;
@@ -370,10 +371,10 @@ class surface {
 
 template <typename detector_t, typename descr_t>
 DETRAY_HOST_DEVICE surface(const detector_t &, const descr_t &)
-    ->surface<detector_t>;
+    -> surface<detector_t>;
 
 template <typename detector_t>
 DETRAY_HOST_DEVICE surface(const detector_t &, const geometry::barcode)
-    ->surface<detector_t>;
+    -> surface<detector_t>;
 
 }  // namespace detray
diff --git a/core/include/detray/masks/masks.hpp b/core/include/detray/masks/masks.hpp
index 778cc07f0b..9bd1340004 100644
--- a/core/include/detray/masks/masks.hpp
+++ b/core/include/detray/masks/masks.hpp
@@ -181,6 +181,10 @@ class mask {
     DETRAY_HOST_DEVICE
     auto values() const -> const mask_values& { return _values; }
 
+    /// @returns the boundary values
+    DETRAY_HOST_DEVICE
+    auto values() -> mask_values& { return _values; }
+
     /// @returns the volume link - const reference
     DETRAY_HOST_DEVICE
     auto volume_link() const -> const links_type& { return _volume_link; }
diff --git a/core/include/detray/materials/material.hpp b/core/include/detray/materials/material.hpp
index 5de706e236..6acd3087f5 100644
--- a/core/include/detray/materials/material.hpp
+++ b/core/include/detray/materials/material.hpp
@@ -11,10 +11,11 @@
 #include "detray/definitions/math.hpp"
 #include "detray/definitions/qualifiers.hpp"
 #include "detray/definitions/units.hpp"
+#include "detray/utils/invalid_values.hpp"
 
 // System include(s)
-#include <limits>
 #include <ratio>
+#include <sstream>
 
 namespace detray {
 
@@ -26,7 +27,12 @@ struct homogeneous_material_tag {};
 }  // namespace detail
 
 /// Material State
-enum class material_state { e_solid = 0, e_liquid = 1, e_gas = 2 };
+enum class material_state {
+    e_solid = 0,
+    e_liquid = 1,
+    e_gas = 2,
+    e_unknown = 3
+};
 
 template <typename scalar_t, typename R = std::ratio<1, 1>>
 struct material {
@@ -53,7 +59,7 @@ struct material {
     ///
     /// @param rhs is the right hand side to be compared to
     DETRAY_HOST_DEVICE
-    constexpr bool operator==(const material<scalar_t>& rhs) const {
+    constexpr bool operator==(const material<scalar_t> &rhs) const {
         return (m_x0 == rhs.X0() && m_l0 == rhs.L0() && m_ar == rhs.Ar() &&
                 m_z == rhs.Z());
     }
@@ -105,6 +111,45 @@ struct material {
         }
     }
 
+    /// @returns a string that contains the material details
+    DETRAY_HOST
+    std::string to_string() const {
+        std::stringstream strm;
+        strm << "material: ";
+        strm << " X0: " << m_x0;
+        strm << " | L0: " << m_l0;
+        strm << " | Z: " << m_z;
+
+        strm << " | state: ";
+        switch (m_state) {
+            case material_state::e_solid: {
+                strm << "solid";
+                break;
+            }
+            case material_state::e_liquid: {
+                strm << "liquid";
+                break;
+            }
+            case material_state::e_gas: {
+                strm << "gaseous";
+                break;
+            }
+            default: {
+                strm << "unknown";
+                break;
+            }
+        };
+
+        return strm.str();
+    }
+
+    /// @returns a string stream that prints the material details
+    DETRAY_HOST
+    friend std::ostream &operator<<(std::ostream &os, const material &mat) {
+        os << mat.to_string();
+        return os;
+    }
+
     protected:
     DETRAY_HOST_DEVICE
     constexpr scalar_type mass_to_molar_density(double ar, double mass_rho) {
@@ -119,13 +164,13 @@ struct material {
     }
 
     // Material properties
-    scalar_type m_x0 = std::numeric_limits<scalar_type>::infinity();
-    scalar_type m_l0 = std::numeric_limits<scalar_type>::infinity();
+    scalar_type m_x0 = detail::invalid_value<scalar_type>();
+    scalar_type m_l0 = detail::invalid_value<scalar_type>();
     scalar_type m_ar = 0.f;
     scalar_type m_z = 0.f;
     scalar_type m_mass_rho = 0.f;
     scalar_type m_molar_rho = 0.f;
-    material_state m_state = material_state::e_solid;
+    material_state m_state = material_state::e_unknown;
 };
 
 // Macro for declaring the predefined materials (with Density effect data)
diff --git a/core/include/detray/materials/material_map.hpp b/core/include/detray/materials/material_map.hpp
new file mode 100644
index 0000000000..b0524e2e24
--- /dev/null
+++ b/core/include/detray/materials/material_map.hpp
@@ -0,0 +1,29 @@
+/** 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
+
+// Detray include(s)
+#include "detray/definitions/containers.hpp"
+#include "detray/materials/material_slab.hpp"
+#include "detray/surface_finders/grid/grid.hpp"
+#include "detray/tools/grid_factory.hpp"
+
+namespace detray {
+
+/// Definition of binned material
+template <typename axes_shape, typename scalar_t,
+          typename container_t = host_container_types, bool owning = false>
+using material_map = grid<coordinate_axes<axes_shape, owning, container_t>,
+                          material_slab<scalar_t>, simple_serializer, replacer>;
+
+/// How to build material maps of various shapes
+template <typename scalar_t = detray::scalar>
+using material_map_factory =
+    grid_factory<material_slab<scalar_t>, simple_serializer, replacer>;
+
+}  // namespace detray
diff --git a/core/include/detray/materials/material_rod.hpp b/core/include/detray/materials/material_rod.hpp
index 20a124c5a4..83200b4c55 100644
--- a/core/include/detray/materials/material_rod.hpp
+++ b/core/include/detray/materials/material_rod.hpp
@@ -91,6 +91,16 @@ struct material_rod : public detail::homogeneous_material_tag {
         return this->path_segment(cos_inc_angle, approach) / m_material.L0();
     }
 
+    /// @returns a string stream that prints the material details
+    DETRAY_HOST
+    friend std::ostream& operator<<(std::ostream& os, const material_rod& mat) {
+        os << "rod: ";
+        os << mat.get_material().to_string();
+        os << " | radius: " << mat.radius() << "mm";
+
+        return os;
+    }
+
     private:
     material_type m_material = {};
     scalar_type m_radius = std::numeric_limits<scalar>::epsilon();
diff --git a/core/include/detray/materials/material_slab.hpp b/core/include/detray/materials/material_slab.hpp
index eba8a03d3f..5285269167 100644
--- a/core/include/detray/materials/material_slab.hpp
+++ b/core/include/detray/materials/material_slab.hpp
@@ -15,6 +15,7 @@
 
 // System include(s)
 #include <limits>
+#include <ostream>
 
 namespace detray {
 
@@ -93,6 +94,17 @@ struct material_slab : public detail::homogeneous_material_tag {
         return m_thickness_in_L0 / cos_inc_angle;
     }
 
+    /// @returns a string stream that prints the material details
+    DETRAY_HOST
+    friend std::ostream& operator<<(std::ostream& os,
+                                    const material_slab& mat) {
+        os << "slab: ";
+        os << mat.get_material().to_string();
+        os << " | thickness: " << mat.thickness() << "mm";
+
+        return os;
+    }
+
     private:
     material_type m_material = {};
     scalar_type m_thickness = std::numeric_limits<scalar>::epsilon();
diff --git a/core/include/detray/propagator/actors/pointwise_material_interactor.hpp b/core/include/detray/propagator/actors/pointwise_material_interactor.hpp
index cec6deee59..eed84494a7 100644
--- a/core/include/detray/propagator/actors/pointwise_material_interactor.hpp
+++ b/core/include/detray/propagator/actors/pointwise_material_interactor.hpp
@@ -39,7 +39,6 @@ struct pointwise_material_interactor : actor {
         scalar_type mass{105.7f * unit<scalar_type>::MeV};
         /// The particle pdg
         int pdg = 13;  // default muon
-
         /// Evaluated energy loss
         scalar_type e_loss{0.f};
         /// Evaluated projected scattering angle
@@ -67,54 +66,77 @@ struct pointwise_material_interactor : actor {
 
         template <typename material_group_t, typename index_t>
         DETRAY_HOST_DEVICE inline bool operator()(
-            const material_group_t &material_group,
-            const index_t &material_range, state &s,
+            const material_group_t &material_group, const index_t &mat_index,
+            state &s,
             const bound_track_parameters<transform3_type> &bound_params,
             const scalar_type cos_inc_angle, const scalar_type approach) const {
 
-            const scalar qop = bound_params.qop();
-            const scalar charge = bound_params.charge();
+            const auto &mat = this->get_material(material_group, mat_index,
+                                                 bound_params.bound_local());
 
-            bool success = false;
+            // return early in case of vacuum or zero thickness
+            if (not mat) {
+                return false;
+            }
 
-            for (const auto &mat :
-                 detray::ranges::subrange(material_group, material_range)) {
+            const scalar qop = bound_params.qop();
+            const scalar charge = bound_params.charge();
 
-                // return early in case of vacuum or zero thickness
-                if (not mat) {
-                    continue;
-                }
-                success |= true;
+            const scalar_type path_segment{
+                mat.path_segment(cos_inc_angle, approach)};
 
-                const scalar_type path_segment{
-                    mat.path_segment(cos_inc_angle, approach)};
+            // Energy Loss
+            if (s.do_energy_loss) {
+                s.e_loss = interaction_type().compute_energy_loss_bethe(
+                    path_segment, mat, s.pdg, s.mass, qop, charge);
+            }
 
-                // Energy Loss
-                if (s.do_energy_loss) {
-                    s.e_loss = interaction_type().compute_energy_loss_bethe(
+            // @todo: include the radiative loss (Bremsstrahlung)
+            if (s.do_energy_loss && s.do_covariance_transport) {
+                s.sigma_qop =
+                    interaction_type().compute_energy_loss_landau_sigma_QOverP(
                         path_segment, mat, s.pdg, s.mass, qop, charge);
-                }
-
-                // @todo: include the radiative loss (Bremsstrahlung)
-                if (s.do_energy_loss && s.do_covariance_transport) {
-                    s.sigma_qop =
-                        interaction_type()
-                            .compute_energy_loss_landau_sigma_QOverP(
-                                path_segment, mat, s.pdg, s.mass, qop, charge);
-                }
+            }
 
-                // Covariance update
-                if (s.do_multiple_scattering) {
-                    // @todo: use momentum before or after energy loss in
-                    // backward mode?
-                    s.projected_scattering_angle =
-                        interaction_type().compute_multiple_scattering_theta0(
-                            mat.path_segment_in_X0(cos_inc_angle, approach),
-                            s.pdg, s.mass, qop, charge);
-                }
+            // Covariance update
+            if (s.do_multiple_scattering) {
+                // @todo: use momentum before or after energy loss in
+                // backward mode?
+                s.projected_scattering_angle =
+                    interaction_type().compute_multiple_scattering_theta0(
+                        mat.path_segment_in_X0(cos_inc_angle, approach), s.pdg,
+                        s.mass, qop, charge);
             }
 
-            return success;
+            return true;
+        }
+
+        private:
+        /// Access to material slabs or rods in a homogeneous material
+        /// description
+        template <class material_coll_t, class point_t,
+                  std::enable_if_t<not detail::is_grid_v<
+                                       typename material_coll_t::value_type>,
+                                   bool> = true>
+        inline constexpr decltype(auto) get_material(
+            const material_coll_t &material_coll, const dindex idx,
+            const point_t &) const noexcept {
+            return material_coll[idx];
+        }
+
+        /// Access to material slabs in a material map
+        template <class material_coll_t, class point_t,
+                  std::enable_if_t<
+                      detail::is_grid_v<typename material_coll_t::value_type>,
+                      bool> = true>
+        inline constexpr decltype(auto) get_material(
+            const material_coll_t &material_coll, const dindex idx,
+            const point_t &loc_point) const noexcept {
+            // Get the material grid instance
+            const auto &mat_grid = material_coll[idx];
+
+            // Find the material slab (only one entry per bin)
+            return *mat_grid.search(loc_point);
         }
     };
 
@@ -126,8 +148,8 @@ struct pointwise_material_interactor : actor {
 
         const auto &navigation = prop_state._navigation;
 
-        // Do material interaction when the track is on surface
-        if (navigation.is_on_module()) {
+        // Do material interaction when the track is on material surface
+        if (navigation.encountered_material()) {
 
             auto &stepping = prop_state._stepping;
 
diff --git a/core/include/detray/propagator/navigator.hpp b/core/include/detray/propagator/navigator.hpp
index 244c3981e4..daa63caaaa 100644
--- a/core/include/detray/propagator/navigator.hpp
+++ b/core/include/detray/propagator/navigator.hpp
@@ -344,6 +344,14 @@ class navigator {
             return m_status == navigation::status::e_on_portal;
         }
 
+        /// Helper method to check the track has encountered material
+        DETRAY_HOST_DEVICE
+        inline auto encountered_material() const -> bool {
+            return (is_on_module() or is_on_portal()) and
+                   (current()->sf_desc.material().id() !=
+                    detector_t::materials::id::e_none);
+        }
+
         /// Helper method to check if a kernel is exhausted - const
         DETRAY_HOST_DEVICE
         inline auto is_exhausted() const -> bool {
diff --git a/core/include/detray/surface_finders/brute_force_finder.hpp b/core/include/detray/surface_finders/brute_force_finder.hpp
index 50ad35b0e9..642e74c357 100644
--- a/core/include/detray/surface_finders/brute_force_finder.hpp
+++ b/core/include/detray/surface_finders/brute_force_finder.hpp
@@ -75,6 +75,9 @@ class brute_force_collection {
         /// @returns an iterator over all surfaces in the data structure
         DETRAY_HOST_DEVICE constexpr auto all() const { return *this; }
 
+        /// @returns an iterator over all surfaces in the data structure
+        DETRAY_HOST_DEVICE auto all() { return *this; }
+
         /// @return the maximum number of surface candidates during a
         /// neighborhood lookup
         DETRAY_HOST_DEVICE constexpr auto n_max_candidates() const
@@ -152,6 +155,12 @@ class brute_force_collection {
         return {m_surfaces, dindex_range{m_offsets[i], m_offsets[i + 1u]}};
     }
 
+    /// Create brute force surface finder from surface container - const
+    DETRAY_HOST_DEVICE
+    auto operator[](const size_type i) -> value_type {
+        return {m_surfaces, dindex_range{m_offsets[i], m_offsets[i + 1u]}};
+    }
+
     /// Add a new surface collection
     template <
         typename sf_container_t,
diff --git a/core/include/detray/surface_finders/grid/grid.hpp b/core/include/detray/surface_finders/grid/grid.hpp
index 17939a75d8..0d3c0b2c62 100644
--- a/core/include/detray/surface_finders/grid/grid.hpp
+++ b/core/include/detray/surface_finders/grid/grid.hpp
@@ -284,11 +284,21 @@ class grid {
 
     /// Find the value of a single bin - const
     ///
-    /// @param p is point in the local frame
+    /// @param p is point in the local (bound) frame
     ///
     /// @return the iterable view of the bin content
-    DETRAY_HOST_DEVICE auto search(
-        const typename local_frame_type::point3 &p) const {
+    template <typename point_t>
+    DETRAY_HOST_DEVICE auto search(const point_t &p) const {
+        return bin(m_axes.bins(p));
+    }
+
+    /// Find the value of a single bin
+    ///
+    /// @param p is point in the local (bound) frame
+    ///
+    /// @return the iterable view of the bin content
+    template <typename point_t>
+    DETRAY_HOST_DEVICE auto search(const point_t &p) {
         return bin(m_axes.bins(p));
     }
 
diff --git a/core/include/detray/surface_finders/grid/grid_collection.hpp b/core/include/detray/surface_finders/grid/grid_collection.hpp
index 9e86f901f5..26ebfef400 100644
--- a/core/include/detray/surface_finders/grid/grid_collection.hpp
+++ b/core/include/detray/surface_finders/grid/grid_collection.hpp
@@ -115,10 +115,48 @@ class grid_collection<
         return static_cast<dindex>(m_offsets.size());
     }
 
+    /// @returns an iterator that points to the first grid
+    /// @note Not implemented!
+    DETRAY_HOST_DEVICE
+    constexpr auto begin() noexcept -> bool { return true; }
+
+    /// @returns an iterator that points to the coll. end
+    /// @note Not implemented!
+    DETRAY_HOST_DEVICE
+    constexpr auto end() noexcept -> bool { return false; }
+
     /// @returns the number of grids in the collection - const
     DETRAY_HOST_DEVICE
     constexpr auto empty() const noexcept -> bool { return m_offsets.empty(); }
 
+    /// @brief Resize the underlying containers
+    /// @note Not defined! The amount of memory can differ for every grid
+    DETRAY_HOST_DEVICE
+    constexpr void resize(std::size_t) noexcept { /*Not defined*/
+    }
+
+    /// @brief Reserve memory
+    /// @note Not defined! The amount of memory can differ for every grid
+    DETRAY_HOST_DEVICE
+    constexpr void reserve(std::size_t) noexcept { /*Not defined*/
+    }
+
+    /// Removes all data from the grid collection containers
+    DETRAY_HOST_DEVICE
+    constexpr void clear() noexcept {
+        m_offsets.clear();
+        m_bins.clear();
+        m_axes_data.clear();
+        m_bin_edges.clear();
+    }
+
+    /// Insert a number of grids
+    /// @note Not defined! There is no grid iterator implementation
+    template <typename... Args>
+    DETRAY_HOST_DEVICE constexpr void insert(Args &&...) noexcept {
+        /*Not defined*/
+    }
+
     /// @returns the offsets for the grids in the bin storage - const
     DETRAY_HOST_DEVICE
     constexpr auto offsets() const -> const vector_type<size_type> & {
diff --git a/core/include/detray/tools/volume_builder.hpp b/core/include/detray/tools/volume_builder.hpp
index 39d6610e47..278ee3a0f8 100644
--- a/core/include/detray/tools/volume_builder.hpp
+++ b/core/include/detray/tools/volume_builder.hpp
@@ -11,6 +11,8 @@
 #include "detray/definitions/geometry.hpp"
 #include "detray/geometry/surface.hpp"
 #include "detray/tools/volume_builder_interface.hpp"
+// TODO: Remove once 'is_grid_v' can be removed
+#include "detray/surface_finders/accelerator_grid.hpp"
 
 // System include(s)
 #include <memory>
@@ -202,10 +204,13 @@ struct mask_index_update {
 struct material_index_update {
 
     template <typename group_t, typename index_t, typename surface_t>
-    DETRAY_HOST inline void operator()(const group_t& group,
-                                       const index_t& /*index*/,
-                                       surface_t& sf) const {
-        sf.update_material(static_cast<dindex>(group.size()));
+    DETRAY_HOST inline void operator()(
+        [[maybe_unused]] const group_t& group,
+        [[maybe_unused]] const index_t& /*index*/,
+        [[maybe_unused]] surface_t& sf) const {
+        if constexpr (!detail::is_grid_v<typename group_t::value_type>) {
+            sf.update_material(static_cast<dindex>(group.size()));
+        }
     }
 };
 
diff --git a/io/include/detray/io/common/detail/type_traits.hpp b/io/include/detray/io/common/detail/type_traits.hpp
index a895f2a55e..179428e6f6 100644
--- a/io/include/detray/io/common/detail/type_traits.hpp
+++ b/io/include/detray/io/common/detail/type_traits.hpp
@@ -157,22 +157,21 @@ struct mask_info<io::detail::mask_shape::trapezoid2, detector_t,
 /// @}
 
 template <class detector_t, typename = void>
-struct is_homogeneous_material : public std::false_type {};
+struct has_homogeneous_material : public std::false_type {};
 
 /// Is the value type in the detector material store a simple material or is it
 /// wrapped in another class (e.g. grids for material maps)
 template <class detector_t>
-struct is_homogeneous_material<
+struct has_homogeneous_material<
     detector_t,
-    std::enable_if_t<
-        std::is_base_of_v<detail::homogeneous_material_tag,
-                          typename detail::tuple_element_t<
-                              0, typename detector_t::material_container::
-                                     tuple_type>::value_type>,
-        void>> : public std::true_type {};
+    std::enable_if_t<detector_t::materials::template is_defined<
+                         material_slab<typename detector_t::scalar_type>>() or
+                         detector_t::materials::template is_defined<
+                             material_rod<typename detector_t::scalar_type>>(),
+                     void>> : public std::true_type {};
 
 template <typename T>
-inline constexpr bool is_homogeneous_material_v =
-    is_homogeneous_material<T>::value;
+inline constexpr bool has_homogeneous_material_v =
+    has_homogeneous_material<T>::value;
 
 }  // namespace detray::detail
diff --git a/io/include/detray/io/common/detector_writer.hpp b/io/include/detray/io/common/detector_writer.hpp
index 97b82b7d3a..3c0b865ee6 100644
--- a/io/include/detray/io/common/detector_writer.hpp
+++ b/io/include/detray/io/common/detector_writer.hpp
@@ -84,7 +84,7 @@ detail::detector_component_writers<detector_t> assemble_writer(
         // Find other writers, depending on the detector type
         if (cfg.write_material()) {
             // Simple material
-            if constexpr (detail::is_homogeneous_material_v<detector_t>) {
+            if constexpr (detail::has_homogeneous_material_v<detector_t>) {
                 writers.template add<json_homogeneous_material_writer>();
             }
             // Material maps
diff --git a/io/include/detray/io/common/homogeneous_material_writer.hpp b/io/include/detray/io/common/homogeneous_material_writer.hpp
index b16bd5ef25..902163e139 100644
--- a/io/include/detray/io/common/homogeneous_material_writer.hpp
+++ b/io/include/detray/io/common/homogeneous_material_writer.hpp
@@ -39,7 +39,6 @@ class homogeneous_material_writer : public writer_interface<detector_t> {
     /// Serialize the header information into its payload
     static homogeneous_material_header_payload write_header(
         const detector_t& det, const std::string_view det_name) {
-        using mat_types = typename detector_t::material_container::value_types;
 
         homogeneous_material_header_payload header_data;
 
@@ -49,15 +48,16 @@ class homogeneous_material_writer : public writer_interface<detector_t> {
 
         header_data.sub_header.emplace();
         auto& mat_sub_header = header_data.sub_header.value();
-        mat_sub_header.n_slabs =
-            materials.template size<mat_types::to_id(0u)>();
+        if constexpr (detector_t::materials::template is_defined<
+                          material_slab<scalar_t>>()) {
+            mat_sub_header.n_slabs =
+                materials.template size<detector_t::materials::id::e_slab>();
+        }
         mat_sub_header.n_rods = 0u;
-        if constexpr (mat_types::n_types == 2u) {
-            // The compiler looks at this code, even if the number of material
-            // types is one. Therefore, "mat_types::n_types - 1" is safer to use
+        if constexpr (detector_t::materials::template is_defined<
+                          material_rod<scalar_t>>()) {
             mat_sub_header.n_rods =
-                materials
-                    .template size<mat_types::to_id(mat_types::n_types - 1)>();
+                materials.template size<detector_t::materials::id::e_rod>();
         }
 
         return header_data;
@@ -161,8 +161,20 @@ class homogeneous_material_writer : public writer_interface<detector_t> {
         template <typename material_group_t, typename index_t>
         inline auto operator()(const material_group_t& material_group,
                                const index_t& index) const {
-            return homogeneous_material_writer<detector_t>::serialize(
-                material_group[index], index);
+            using material_t = typename material_group_t::value_type;
+
+            constexpr bool is_slab =
+                std::is_same_v<material_t, material_slab<scalar_t>>;
+            constexpr bool is_rod =
+                std::is_same_v<material_t, material_rod<scalar_t>>;
+
+            if constexpr (is_slab or is_rod) {
+
+                return homogeneous_material_writer<detector_t>::serialize(
+                    material_group[index], index);
+            } else {
+                return material_slab_payload{};
+            }
         }
     };
 };
diff --git a/tests/common/include/tests/common/test_toy_detector.hpp b/tests/common/include/tests/common/test_toy_detector.hpp
index 790d7b643e..885251b2c7 100644
--- a/tests/common/include/tests/common/test_toy_detector.hpp
+++ b/tests/common/include/tests/common/test_toy_detector.hpp
@@ -97,6 +97,8 @@ inline bool test_toy_detector(
         (sf_finders.template size<sf_finder_ids::e_disc_grid>() != 0u);
     const bool has_material =
         (materials.template size<material_ids::e_slab>() != 0);
+    const bool has_material_maps =
+        (materials.template size<material_ids::e_disc2_map>() != 0);
 
     // Check number of geomtery objects
     EXPECT_EQ(volumes.size(), 20u);
@@ -112,8 +114,13 @@ inline bool test_toy_detector(
                   4);
         EXPECT_EQ(sf_finders.template size<sf_finder_ids::e_disc_grid>(), 6);
     }
-    if (has_material) {
+    if (has_material and !has_material_maps) {
         EXPECT_EQ(materials.template size<material_ids::e_slab>(), 3244u);
+    } else if (has_material and has_material_maps) {
+        EXPECT_EQ(materials.template size<material_ids::e_slab>(), 3141u);
+        EXPECT_EQ(materials.template size<material_ids::e_cylinder2_map>(),
+                  51u);
+        EXPECT_EQ(materials.template size<material_ids::e_disc2_map>(), 52u);
     }
 
     /** Test the links of a volume.
@@ -163,18 +170,18 @@ inline bool test_toy_detector(
                 const auto volume_link =
                     masks.template visit<volume_link_getter>(sf_itr->mask());
                 EXPECT_EQ(volume_link, volume_links[pti - range[0]]);
-                if (has_material) {
+                /*if (has_material) {
                     EXPECT_EQ(sf_itr->material(), material_index);
                     EXPECT_EQ(
                         materials.template get<
                             material_ids::e_slab>()[sf_itr->material().index()],
                         mat);
-                }
+                }*/
 
                 ++sf_itr;
                 ++trf_index;
                 ++mask_link;
-                ++material_index;
+                //++material_index;
             }
         };
 
@@ -205,18 +212,19 @@ inline bool test_toy_detector(
                 const auto volume_link =
                     masks.template visit<volume_link_getter>(sf_itr->mask());
                 EXPECT_EQ(volume_link, volume_links[0]);
-                if (has_material) {
+                /*if (has_material) {
                     EXPECT_EQ(sf_itr->material(), material_index);
                     EXPECT_EQ(
                         materials.template get<
                             material_ids::e_slab>()[sf_itr->material().index()],
-                        mat);
-                }
+                        mat)
+                        << sf_itr->material();
+                }*/
 
                 ++sf_itr;
                 ++trf_index;
                 ++mask_index;
-                ++material_index;
+                //++material_index;
             }
         };
 
diff --git a/tests/unit_tests/cpu/CMakeLists.txt b/tests/unit_tests/cpu/CMakeLists.txt
index 246d60beeb..141198a7db 100644
--- a/tests/unit_tests/cpu/CMakeLists.txt
+++ b/tests/unit_tests/cpu/CMakeLists.txt
@@ -47,6 +47,7 @@ macro( detray_add_cpu_test algebra )
       "masks_unbounded.cpp"
       "masks_unmasked.cpp"
       "material_interaction.cpp"
+      "material_maps.cpp"
       "materials.cpp"
       "scattering.cpp"
       "sf_finder_brute_force.cpp"
diff --git a/tests/unit_tests/cpu/material_maps.cpp b/tests/unit_tests/cpu/material_maps.cpp
new file mode 100644
index 0000000000..384f58b9fe
--- /dev/null
+++ b/tests/unit_tests/cpu/material_maps.cpp
@@ -0,0 +1,237 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Detray include(s)
+#include "detray/definitions/indexing.hpp"
+#include "detray/masks/masks.hpp"
+#include "detray/materials/material_map.hpp"
+
+// GTest include(s)
+#include <gtest/gtest.h>
+
+using namespace detray;
+using namespace detray::n_axis;
+
+using material_t = typename material_map_factory<scalar>::value_type;
+
+namespace {
+
+material_map_factory<scalar> mat_map_factory{};
+
+}  // anonymous namespace
+
+/// Unittest: Test the construction of an annulus shaped material map
+GTEST_TEST(detray_material, annulus_map) {
+
+    constexpr scalar minR{7.2f * unit<scalar>::mm};
+    constexpr scalar maxR{12.0f * unit<scalar>::mm};
+    constexpr scalar minPhi{0.74195f};
+    constexpr scalar maxPhi{1.33970f};
+
+    mask<annulus2D<>> ann2{0u, minR, maxR, minPhi, maxPhi, 0.f, -2.f, 2.f};
+
+    auto annulus_map = mat_map_factory.template new_grid<>(ann2, {10u, 20u});
+
+    EXPECT_EQ(annulus_map.Dim, 2u);
+    EXPECT_EQ(annulus_map.nbins(), 200u);
+
+    auto r_axis = annulus_map.get_axis<label::e_r>();
+    EXPECT_EQ(r_axis.nbins(), 10u);
+    EXPECT_EQ(r_axis.min(), minR);
+    EXPECT_EQ(r_axis.max(), maxR);
+
+    auto phi_axis = annulus_map.get_axis<label::e_phi>();
+    EXPECT_EQ(phi_axis.nbins(), 20u);
+    EXPECT_EQ(phi_axis.min(), -minPhi);
+    EXPECT_EQ(phi_axis.max(), maxPhi);
+
+    // Add some material entries
+    scalar thickness = 2.f * unit<scalar>::mm;
+    for (dindex gbin = 0; gbin < annulus_map.nbins(); ++gbin) {
+        annulus_map.populate(gbin,
+                             material_t(silicon_tml<scalar>{}, thickness));
+        thickness += 1.f * unit<scalar>::mm;
+    }
+
+    EXPECT_EQ(annulus_map.at(0, 0),
+              material_t(silicon_tml<scalar>{}, 2.f * unit<scalar>::mm));
+    EXPECT_EQ(annulus_map.at(1, 0),
+              material_t(silicon_tml<scalar>{}, 3.f * unit<scalar>::mm));
+    EXPECT_EQ(annulus_map.at(22, 0),
+              material_t(silicon_tml<scalar>{}, 24.f * unit<scalar>::mm));
+    EXPECT_EQ(annulus_map.at(199, 0),
+              material_t(silicon_tml<scalar>{}, 201.f * unit<scalar>::mm));
+    EXPECT_FALSE(annulus_map.at(199, 0) ==
+                 material_t(gold<scalar>{}, 201.f * unit<scalar>::mm));
+}
+
+/// Unittest: Test the construction of a cylindrical material map
+GTEST_TEST(detray_material, cylinder_map) {
+
+    constexpr scalar r{3.f * unit<scalar>::mm};
+    constexpr scalar hz{4.f * unit<scalar>::mm};
+
+    mask<cylinder2D<>> cyl{0u, r, -hz, hz};
+
+    auto cylinder_map = mat_map_factory.template new_grid<>(cyl, {10u, 20u});
+
+    EXPECT_EQ(cylinder_map.Dim, 2u);
+    EXPECT_EQ(cylinder_map.nbins(), 200u);
+
+    auto rphi_axis = cylinder_map.get_axis<label::e_rphi>();
+    EXPECT_EQ(rphi_axis.nbins(), 10u);
+    EXPECT_EQ(rphi_axis.min(), -constant<scalar>::pi * r);
+    EXPECT_EQ(rphi_axis.max(), constant<scalar>::pi * r);
+
+    auto z_axis = cylinder_map.get_axis<label::e_cyl_z>();
+    EXPECT_EQ(z_axis.nbins(), 20u);
+    EXPECT_EQ(z_axis.min(), -hz);
+    EXPECT_EQ(z_axis.max(), hz);
+
+    scalar thickness = 2.f * unit<scalar>::mm;
+    for (dindex gbin = 0; gbin < cylinder_map.nbins(); ++gbin) {
+        cylinder_map.populate(gbin, material_t(vacuum<scalar>{}, thickness));
+        thickness += 1.f * unit<scalar>::mm;
+    }
+
+    EXPECT_EQ(cylinder_map.at(0, 0),
+              material_t(vacuum<scalar>{}, 2.f * unit<scalar>::mm));
+    EXPECT_EQ(cylinder_map.at(1, 0),
+              material_t(vacuum<scalar>{}, 3.f * unit<scalar>::mm));
+    EXPECT_EQ(cylinder_map.at(22, 0),
+              material_t(vacuum<scalar>{}, 24.f * unit<scalar>::mm));
+    EXPECT_EQ(cylinder_map.at(199, 0),
+              material_t(vacuum<scalar>{}, 201.f * unit<scalar>::mm));
+    EXPECT_FALSE(cylinder_map.at(199, 0) ==
+                 material_t(gold<scalar>{}, 201.f * unit<scalar>::mm));
+}
+
+/// Unittest: Test the construction of a rectangular material map
+GTEST_TEST(detray_material, rectangle_map) {
+
+    constexpr scalar hx{1.f * unit<scalar>::mm};
+    constexpr scalar hy{9.3f * unit<scalar>::mm};
+
+    mask<rectangle2D<>> r2{0u, hx, hy};
+
+    auto rectangle_map = mat_map_factory.template new_grid<>(r2, {10u, 20u});
+
+    EXPECT_EQ(rectangle_map.Dim, 2u);
+    EXPECT_EQ(rectangle_map.nbins(), 200u);
+
+    auto x_axis = rectangle_map.get_axis<label::e_x>();
+    EXPECT_EQ(x_axis.nbins(), 10u);
+    EXPECT_EQ(x_axis.min(), -hx);
+    EXPECT_EQ(x_axis.max(), hx);
+
+    auto y_axis = rectangle_map.get_axis<label::e_y>();
+    EXPECT_EQ(y_axis.nbins(), 20u);
+    EXPECT_EQ(y_axis.min(), -hy);
+    EXPECT_EQ(y_axis.max(), hy);
+
+    scalar thickness = 2.f * unit<scalar>::mm;
+    for (dindex gbin = 0; gbin < rectangle_map.nbins(); ++gbin) {
+        rectangle_map.populate(gbin,
+                               material_t(oxygen_gas<scalar>{}, thickness));
+        thickness += 1.f * unit<scalar>::mm;
+    }
+
+    EXPECT_EQ(rectangle_map.at(0, 0),
+              material_t(oxygen_gas<scalar>{}, 2.f * unit<scalar>::mm));
+    EXPECT_EQ(rectangle_map.at(1, 0),
+              material_t(oxygen_gas<scalar>{}, 3.f * unit<scalar>::mm));
+    EXPECT_EQ(rectangle_map.at(22, 0),
+              material_t(oxygen_gas<scalar>{}, 24.f * unit<scalar>::mm));
+    EXPECT_EQ(rectangle_map.at(199, 0),
+              material_t(oxygen_gas<scalar>{}, 201.f * unit<scalar>::mm));
+    EXPECT_FALSE(rectangle_map.at(199, 0) ==
+                 material_t(gold<scalar>{}, 201.f * unit<scalar>::mm));
+}
+
+/// Unittest: Test the construction of a ring shaped material map
+GTEST_TEST(detray_material, disc_map) {
+
+    constexpr scalar inner_r{0.f * unit<scalar>::mm};
+    constexpr scalar outer_r{3.5f * unit<scalar>::mm};
+
+    mask<ring2D<>> r2{0u, inner_r, outer_r};
+
+    auto disc_map = mat_map_factory.template new_grid<>(r2, {10u, 20u});
+
+    EXPECT_EQ(disc_map.Dim, 2u);
+    EXPECT_EQ(disc_map.nbins(), 200u);
+
+    auto r_axis = disc_map.get_axis<label::e_r>();
+    EXPECT_EQ(r_axis.nbins(), 10u);
+    EXPECT_EQ(r_axis.min(), inner_r);
+    EXPECT_EQ(r_axis.max(), outer_r);
+
+    auto phi_axis = disc_map.get_axis<label::e_phi>();
+    EXPECT_EQ(phi_axis.nbins(), 20u);
+    EXPECT_EQ(phi_axis.min(), -constant<scalar>::pi);
+    EXPECT_EQ(phi_axis.max(), constant<scalar>::pi);
+
+    scalar thickness = 2.f * unit<scalar>::mm;
+    for (dindex gbin = 0; gbin < disc_map.nbins(); ++gbin) {
+        disc_map.populate(gbin, material_t(aluminium<scalar>{}, thickness));
+        thickness += 1.f * unit<scalar>::mm;
+    }
+
+    EXPECT_EQ(disc_map.at(0, 0),
+              material_t(aluminium<scalar>{}, 2.f * unit<scalar>::mm));
+    EXPECT_EQ(disc_map.at(1, 0),
+              material_t(aluminium<scalar>{}, 3.f * unit<scalar>::mm));
+    EXPECT_EQ(disc_map.at(22, 0),
+              material_t(aluminium<scalar>{}, 24.f * unit<scalar>::mm));
+    EXPECT_EQ(disc_map.at(199, 0),
+              material_t(aluminium<scalar>{}, 201.f * unit<scalar>::mm));
+    EXPECT_FALSE(disc_map.at(199, 0) ==
+                 material_t(gold<scalar>{}, 201.f * unit<scalar>::mm));
+}
+
+/// Unittest: Test the construction of a trapezoid material map
+GTEST_TEST(detray_material, trapezoid_map) {
+
+    constexpr scalar hx_miny{1.f * unit<scalar>::mm};
+    constexpr scalar hx_maxy{3.f * unit<scalar>::mm};
+    constexpr scalar hy{2.f * unit<scalar>::mm};
+    constexpr scalar divisor{1.f / (2.f * hy)};
+
+    mask<trapezoid2D<>> t2{0u, hx_miny, hx_maxy, hy, divisor};
+
+    auto trapezoid_map = mat_map_factory.template new_grid<>(t2, {10u, 20u});
+
+    EXPECT_EQ(trapezoid_map.Dim, 2u);
+    EXPECT_EQ(trapezoid_map.nbins(), 200u);
+
+    auto x_axis = trapezoid_map.get_axis<label::e_x>();
+    EXPECT_EQ(x_axis.nbins(), 10u);
+    EXPECT_EQ(x_axis.min(), -hx_maxy);
+    EXPECT_EQ(x_axis.max(), hx_maxy);
+
+    auto y_axis = trapezoid_map.get_axis<label::e_y>();
+    EXPECT_EQ(y_axis.nbins(), 20u);
+    EXPECT_EQ(y_axis.min(), -hy);
+    EXPECT_EQ(y_axis.max(), hy);
+
+    scalar thickness = 2.f * unit<scalar>::mm;
+    for (dindex gbin = 0; gbin < trapezoid_map.nbins(); ++gbin) {
+        trapezoid_map.populate(gbin, material_t(gold<scalar>{}, thickness));
+        thickness += 1.f * unit<scalar>::mm;
+    }
+
+    EXPECT_EQ(trapezoid_map.at(0, 0),
+              material_t(gold<scalar>{}, 2.f * unit<scalar>::mm));
+    EXPECT_EQ(trapezoid_map.at(1, 0),
+              material_t(gold<scalar>{}, 3.f * unit<scalar>::mm));
+    EXPECT_EQ(trapezoid_map.at(22, 0),
+              material_t(gold<scalar>{}, 24.f * unit<scalar>::mm));
+    EXPECT_EQ(trapezoid_map.at(199, 0),
+              material_t(gold<scalar>{}, 201.f * unit<scalar>::mm));
+    EXPECT_FALSE(trapezoid_map.at(199, 0) ==
+                 material_t(aluminium<scalar>{}, 201.f * unit<scalar>::mm));
+}
diff --git a/tutorials/include/detray/tutorial/detector_metadata.hpp b/tutorials/include/detray/tutorial/detector_metadata.hpp
index 13aee548c8..45f7d174dc 100644
--- a/tutorials/include/detray/tutorial/detector_metadata.hpp
+++ b/tutorials/include/detray/tutorial/detector_metadata.hpp
@@ -113,9 +113,10 @@ struct my_metadata {
     /// How to store and link materials. The material does not make use of
     /// conditions data ( @c empty_context )
     template <template <typename...> class tuple_t = dtuple,
-              template <typename...> class vector_t = dvector>
-    using material_store = regular_multi_store<material_ids, empty_context,
-                                               tuple_t, vector_t, slab>;
+              typename container_t = host_container_types>
+    using material_store =
+        multi_store<material_ids, empty_context, tuple_t,
+                    typename container_t::template vector_type<slab>>;
 
     /// Surface descriptor type used for sensitives, passives and portals
     /// It holds the indices to the surface data in the detector data stores
diff --git a/utils/include/detray/detectors/create_toy_geometry.hpp b/utils/include/detray/detectors/create_toy_geometry.hpp
index e60887adc9..591c8542ec 100644
--- a/utils/include/detray/detectors/create_toy_geometry.hpp
+++ b/utils/include/detray/detectors/create_toy_geometry.hpp
@@ -17,6 +17,7 @@
 #include "detray/geometry/detector_volume.hpp"
 #include "detray/geometry/surface.hpp"
 #include "detray/io/common/detail/file_handle.hpp"
+#include "detray/materials/mixture.hpp"
 #include "detray/materials/predefined_materials.hpp"
 #include "detray/tools/grid_builder.hpp"
 #include "detray/tools/volume_builder.hpp"
@@ -53,6 +54,15 @@ struct toy_det_config {
     std::string m_bfield_file{!std::getenv("DETRAY_BFIELD_FILE")
                                   ? ""
                                   : std::getenv("DETRAY_BFIELD_FILE")};
+    /// Do material maps on portals
+    bool m_use_material_maps{false};
+    /// Number of bins for material maps
+    std::array<std::size_t, 2> m_cyl_map_bins{20u, 20u};
+    std::array<std::size_t, 2> m_disc_map_bins{3u, 20u};
+    /// Material to be filled into the maps
+    material<scalar> m_material =
+        mixture<scalar, silicon_tml<scalar, std::ratio<9, 10>>,
+                aluminium<scalar, std::ratio<1, 10>>>{};
 
     /// Setters
     /// @{
@@ -76,6 +86,24 @@ struct toy_det_config {
         m_bfield_file = file_name;
         return *this;
     }
+    constexpr toy_det_config &use_material_maps(const bool b) {
+        m_use_material_maps = b;
+        return *this;
+    }
+    constexpr toy_det_config &cyl_map_bins(const std::size_t n_rphi,
+                                           const std::size_t n_z) {
+        m_cyl_map_bins = {n_rphi, n_z};
+        return *this;
+    }
+    constexpr toy_det_config &disc_map_bins(const std::size_t n_r,
+                                            const std::size_t n_phi) {
+        m_disc_map_bins = {n_r, n_phi};
+        return *this;
+    }
+    constexpr toy_det_config &mapped_material(const material<scalar> &mat) {
+        m_material = mat;
+        return *this;
+    }
     /// @}
 
     /// Getters
@@ -84,6 +112,23 @@ struct toy_det_config {
     constexpr unsigned int n_edc_layers() const { return m_n_edc_layers; }
     constexpr const vector3 &bfield_vec() const { return m_bfield_vec; }
     const std::string &bfield_file() const { return m_bfield_file; }
+    constexpr bool use_material_maps() const { return m_use_material_maps; }
+    constexpr const std::array<std::size_t, 2> &cyl_map_bins() const {
+        return m_cyl_map_bins;
+    }
+    constexpr const std::array<std::size_t, 2> &disc_map_bins() const {
+        return m_disc_map_bins;
+    }
+    scalar thickness() const { return 1.f * unit<scalar>::mm; }
+    std::vector<scalar> thickness(
+        const std::array<std::size_t, 2> &nbins) const {
+        return std::vector<scalar>(nbins[0] * nbins[1], 1.f * unit<scalar>::mm);
+    }
+    material<scalar> mapped_material() const { return m_material; }
+    std::vector<material<scalar>> mapped_material(
+        const std::array<std::size_t, 2> &nbins) const {
+        return std::vector<material<scalar>>(nbins[0] * nbins[1], m_material);
+    }
     /// @}
 };
 
@@ -544,7 +589,8 @@ inline void create_endcap_modules(context_t &ctx, volume_type &vol,
  */
 template <typename detector_t>
 inline void add_beampipe(
-    detector_t &det, vecmem::memory_resource &resource,
+    const toy_det_config &cfg, detector_t &det,
+    vecmem::memory_resource &resource,
     typename detector_t::geometry_context &ctx,
     typename detector_t::name_map &names, const unsigned int n_edc_layers,
     const unsigned int n_brl_layers,
@@ -554,10 +600,14 @@ inline void add_beampipe(
 
     // Get the object types handled by the volume
     using object_id = typename detector_t::volume_type::object_id;
+    using material_id = typename detector_t::materials::id;
+    using material_link_t = typename detector_t::surface_type::material_link;
     using nav_link_t = typename detector_t::surface_type::navigation_link;
     constexpr auto leaving_world{detail::invalid_value<nav_link_t>()};
     constexpr auto cyl_id = detector_t::masks::id::e_portal_cylinder2;
 
+    const detail::detector_helper<transform3_t> det_helper{};
+
     scalar max_z{n_edc_layers == 0u ? brl_half_z
                                     : edc_lay_sizes[n_edc_layers - 1u].second};
     scalar min_z{-max_z};
@@ -593,49 +643,61 @@ inline void add_beampipe(
 
     // negative endcap portals
     dindex link = beampipe_idx;
+    auto &material_coll =
+        cfg.use_material_maps() ? det.material_store() : materials;
     for (int i = static_cast<int>(vol_sizes.size()) - 1; i >= 0; --i) {
         volume_link = ++link;
-        detail::detector_helper<transform3_t>().add_cylinder_surface<cyl_id>(
-            beampipe_idx, ctx, surfaces, masks, materials, transforms,
-            edc_inner_r, -vol_sizes[static_cast<unsigned int>(i)].second,
-            -vol_sizes[static_cast<unsigned int>(i)].first, volume_link,
-            vacuum<scalar>(), 0.f * unit<scalar>::mm);
+        auto [neg_cyl, neg_cyl_mask] = det_helper.add_cylinder_surface<cyl_id>(
+            beampipe_idx, ctx, surfaces, masks, transforms, edc_inner_r,
+            -vol_sizes[static_cast<unsigned int>(i)].second,
+            -vol_sizes[static_cast<unsigned int>(i)].first, volume_link);
+
+        det_helper.create_material(cfg, neg_cyl, neg_cyl_mask, material_coll);
     }
 
     // barrel portals
     volume_link = n_brl_layers <= 0u ? leaving_world : link + 1u;
-    detail::detector_helper<transform3_t>().add_cylinder_surface<cyl_id>(
-        beampipe_idx, ctx, surfaces, masks, materials, transforms, edc_inner_r,
-        -brl_half_z, brl_half_z, volume_link, vacuum<scalar>(),
-        0.f * unit<scalar>::mm);
+    auto [barr_cyl, barr_cyl_mask] = det_helper.add_cylinder_surface<cyl_id>(
+        beampipe_idx, ctx, surfaces, masks, transforms, edc_inner_r,
+        -brl_half_z, brl_half_z, volume_link);
+
+    det_helper.create_material(cfg, barr_cyl, barr_cyl_mask, material_coll);
 
     // positive endcap portals
     link += 7u;
     for (unsigned int i = 0u; i < vol_sizes.size(); ++i) {
         volume_link = ++link;
-        detail::detector_helper<transform3_t>().add_cylinder_surface<cyl_id>(
-            beampipe_idx, ctx, surfaces, masks, materials, transforms,
-            edc_inner_r, vol_sizes[i].second, vol_sizes[i].first, volume_link,
-            vacuum<scalar>(), 0.f * unit<scalar>::mm);
+        auto [pos_cyl, pos_cyl_mask] = det_helper.add_cylinder_surface<cyl_id>(
+            beampipe_idx, ctx, surfaces, masks, transforms, edc_inner_r,
+            vol_sizes[i].second, vol_sizes[i].first, volume_link);
+
+        det_helper.create_material(cfg, pos_cyl, pos_cyl_mask, material_coll);
     }
 
     // disc portals
     volume_link = leaving_world;
-    detail::detector_helper<transform3_t>().add_disc_surface(
-        beampipe_idx, ctx, surfaces, masks, materials, transforms,
-        beampipe_vol_size.first, beampipe_vol_size.second, min_z, volume_link,
-        vacuum<scalar>(), 0.f * unit<scalar>::mm);
-    detail::detector_helper<transform3_t>().add_disc_surface(
-        beampipe_idx, ctx, surfaces, masks, materials, transforms,
-        beampipe_vol_size.first, beampipe_vol_size.second, max_z, volume_link,
-        vacuum<scalar>(), 0.f * unit<scalar>::mm);
+    auto [neg_disc, neg_disc_mask] = det_helper.add_disc_surface(
+        beampipe_idx, ctx, surfaces, masks, transforms, beampipe_vol_size.first,
+        beampipe_vol_size.second, min_z, volume_link);
+    det_helper.create_material(cfg, neg_disc, neg_disc_mask, material_coll);
+
+    auto [pos_disc, pos_disc_mask] = det_helper.add_disc_surface(
+        beampipe_idx, ctx, surfaces, masks, transforms, beampipe_vol_size.first,
+        beampipe_vol_size.second, max_z, volume_link);
+    det_helper.create_material(cfg, pos_disc, pos_disc_mask, material_coll);
 
     // This is the beampipe surface
-    detail::detector_helper<transform3_t>()
-        .add_cylinder_surface<detector_t::masks::id::e_cylinder2>(
-            beampipe_idx, ctx, surfaces, masks, materials, transforms,
-            beampipe_r, min_z, max_z, beampipe_idx, beryllium_tml<scalar>(),
-            0.8f * unit<scalar>::mm);
+    auto [bp, bp_mask] =
+        det_helper.add_cylinder_surface<detector_t::masks::id::e_cylinder2>(
+            beampipe_idx, ctx, surfaces, masks, transforms, beampipe_r, min_z,
+            max_z, beampipe_idx);
+
+    materials.template emplace_back<material_id::e_slab>(
+        {}, beryllium_tml<scalar>(), 0.8f * unit<scalar>::mm);
+
+    bp.material() =
+        material_link_t{material_id::e_slab,
+                        materials.template size<material_id::e_slab>() - 1u};
 
     det.add_objects_per_volume(ctx, beampipe, surfaces, masks, transforms,
                                materials);
@@ -659,7 +721,8 @@ inline void add_beampipe(
  */
 template <typename detector_t>
 inline void add_endcap_barrel_connection(
-    detector_t &det, vecmem::memory_resource &resource,
+    const toy_det_config &cfg, detector_t &det,
+    vecmem::memory_resource &resource,
     typename detector_t::geometry_context &ctx,
     typename detector_t::name_map &names, const int side,
     const unsigned int n_brl_layers, const dindex beampipe_idx,
@@ -672,6 +735,8 @@ inline void add_endcap_barrel_connection(
     constexpr auto leaving_world{detail::invalid_value<nav_link_t>()};
     constexpr auto cyl_id = detector_t::masks::id::e_portal_cylinder2;
 
+    const detail::detector_helper<transform3_t> det_helper{};
+
     const scalar sign{static_cast<scalar>(side)};
     const scalar min_z{std::min(sign * gap_lower_z, sign * gap_upper_z)};
     const scalar max_z{std::max(sign * gap_lower_z, sign * gap_upper_z)};
@@ -695,20 +760,24 @@ inline void add_endcap_barrel_connection(
     det.transform_store().emplace_back(ctx, t);
 
     dindex volume_link{beampipe_idx};
-    detail::detector_helper<transform3_t>().add_cylinder_surface<cyl_id>(
-        connector_gap_idx, ctx, surfaces, masks, materials, transforms,
-        edc_inner_r, min_z, max_z, volume_link, vacuum<scalar>(),
-        0.f * unit<scalar>::mm);
+    auto &material_coll =
+        cfg.use_material_maps() ? det.material_store() : materials;
+    auto [inner_cyl, inner_cyl_mask] = det_helper.add_cylinder_surface<cyl_id>(
+        connector_gap_idx, ctx, surfaces, masks, transforms, edc_inner_r, min_z,
+        max_z, volume_link);
+    det_helper.create_material(cfg, inner_cyl, inner_cyl_mask, material_coll);
+
     volume_link = leaving_world;
-    detail::detector_helper<transform3_t>().add_cylinder_surface<cyl_id>(
-        connector_gap_idx, ctx, surfaces, masks, materials, transforms,
-        edc_outer_r, min_z, max_z, volume_link, vacuum<scalar>(),
-        0.f * unit<scalar>::mm);
+    auto [outer_cyl, outer_cyl_mask] = det_helper.add_cylinder_surface<cyl_id>(
+        connector_gap_idx, ctx, surfaces, masks, transforms, edc_outer_r, min_z,
+        max_z, volume_link);
+    det_helper.create_material(cfg, outer_cyl, outer_cyl_mask, material_coll);
+
     volume_link = edc_vol_idx;
-    detail::detector_helper<transform3_t>().add_disc_surface(
-        connector_gap_idx, ctx, surfaces, masks, materials, transforms,
-        edc_inner_r, edc_outer_r, edc_disc_z, volume_link, vacuum<scalar>(),
-        0.f * unit<scalar>::mm);
+    auto [inner_disc, inner_disc_mask] = det_helper.add_disc_surface(
+        connector_gap_idx, ctx, surfaces, masks, transforms, edc_inner_r,
+        edc_outer_r, edc_disc_z, volume_link);
+    det_helper.create_material(cfg, inner_disc, inner_disc_mask, material_coll);
 
     // Get vol sizes in z also for gap volumes
     std::vector<std::pair<scalar, scalar>> vol_sizes;
@@ -721,10 +790,11 @@ inline void add_endcap_barrel_connection(
     volume_link = brl_vol_idx;
     for (unsigned int i = 0u; i < 2u * n_brl_layers - 1u; ++i) {
         volume_link = brl_vol_idx++;
-        detail::detector_helper<transform3_t>().add_disc_surface(
-            connector_gap_idx, ctx, surfaces, masks, materials, transforms,
-            vol_sizes[i].first, vol_sizes[i].second, brl_disc_z, volume_link,
-            vacuum<scalar>(), 0.f * unit<scalar>::mm);
+        auto [outer_disc, outer_disc_mask] = det_helper.add_disc_surface(
+            connector_gap_idx, ctx, surfaces, masks, transforms,
+            vol_sizes[i].first, vol_sizes[i].second, brl_disc_z, volume_link);
+        det_helper.create_material(cfg, outer_disc, outer_disc_mask,
+                                   material_coll);
     }
 
     det.add_objects_per_volume(ctx, connector_gap, surfaces, masks, transforms,
@@ -744,14 +814,17 @@ inline void add_endcap_barrel_connection(
  */
 template <typename edc_module_factory, typename detector_t, typename config_t>
 inline void add_endcap_detector(
-    detector_t &det, vecmem::memory_resource &resource,
+    const toy_det_config &cfg, detector_t &det,
+    vecmem::memory_resource &resource,
     typename detector_t::geometry_context &ctx,
     typename detector_t::name_map &names, dindex n_layers, dindex beampipe_idx,
     const std::vector<std::pair<scalar, scalar>> &lay_sizes,
-    const std::vector<scalar> &lay_positions, config_t cfg) {
+    const std::vector<scalar> &lay_positions, config_t factory_cfg) {
     using nav_link_t = typename detector_t::surface_type::navigation_link;
     constexpr auto leaving_world{detail::invalid_value<nav_link_t>()};
 
+    const detail::detector_helper<transform3_t> det_helper{};
+
     // Generate consecutive linking between volumes (all volume_links for every
     // vol.)
     std::vector<std::vector<dindex>> volume_links_vec;
@@ -765,7 +838,7 @@ inline void add_endcap_detector(
             {beampipe_idx, leaving_world, ++prev_vol_idx, ++next_vol_idx});
     }
     // Edge of the world is flipped
-    if (cfg.side < 0) {
+    if (factory_cfg.side < 0) {
         volume_links_vec.insert(
             volume_links_vec.begin(),
             {beampipe_idx, leaving_world, leaving_world, first_vol_idx + 1u});
@@ -792,12 +865,12 @@ inline void add_endcap_detector(
         vol_sizes.emplace_back(lay_sizes[i].first, lay_sizes[i].second);
     }
 
-    edc_module_factory m_factory{cfg};
+    edc_module_factory m_factory{factory_cfg};
 
     auto vol_size_itr = vol_sizes.begin();
     auto pos_itr = lay_positions.begin();
     // Reverse iteration for negative endcap
-    if (cfg.side < 0) {
+    if (factory_cfg.side < 0) {
         std::advance(vol_size_itr, 2u * n_layers - 2u);
         std::advance(pos_itr, n_layers - 1u);
     }
@@ -805,23 +878,25 @@ inline void add_endcap_detector(
     for (int i = 0; i < 2 * static_cast<int>(n_layers) - 1; ++i) {
         // Every second layer is a gap volume
         is_gap = !is_gap;
-        const scalar sign{static_cast<scalar>(cfg.side)};
+        const scalar sign{static_cast<scalar>(factory_cfg.side)};
         if (is_gap) {
-            detail::detector_helper<transform3_t>().create_cyl_volume(
-                det, resource, ctx, cfg.inner_r, cfg.outer_r,
-                sign * (vol_size_itr + cfg.side * i)->first,
-                sign * (vol_size_itr + cfg.side * i)->second,
+            det_helper.create_cyl_volume(
+                cfg, det, resource, ctx, factory_cfg.inner_r,
+                factory_cfg.outer_r,
+                sign * (vol_size_itr + factory_cfg.side * i)->first,
+                sign * (vol_size_itr + factory_cfg.side * i)->second,
                 volume_links_vec[static_cast<unsigned int>(i)]);
 
             dindex vol_idx = det.volumes().back().index();
             names[vol_idx + 1u] = "gap_" + std::to_string(vol_idx);
 
         } else {
-            m_factory.cfg.edc_position = *(pos_itr + cfg.side * i / 2);
-            detail::detector_helper<transform3_t>().create_cyl_volume(
-                det, resource, ctx, cfg.inner_r, cfg.outer_r,
-                sign * (vol_size_itr + cfg.side * i)->first,
-                sign * (vol_size_itr + cfg.side * i)->second,
+            m_factory.cfg.edc_position = *(pos_itr + factory_cfg.side * i / 2);
+            det_helper.create_cyl_volume(
+                cfg, det, resource, ctx, factory_cfg.inner_r,
+                factory_cfg.outer_r,
+                sign * (vol_size_itr + factory_cfg.side * i)->first,
+                sign * (vol_size_itr + factory_cfg.side * i)->second,
                 volume_links_vec[static_cast<unsigned int>(i)]);
 
             dindex vol_idx = det.volumes().back().index();
@@ -846,17 +921,20 @@ inline void add_endcap_detector(
  */
 template <typename brl_module_factory, typename detector_t, typename config_t>
 inline void add_barrel_detector(
-    detector_t &det, vecmem::memory_resource &resource,
+    const toy_det_config &cfg, detector_t &det,
+    vecmem::memory_resource &resource,
     typename detector_t::geometry_context &ctx,
     typename detector_t::name_map &names, const unsigned int n_layers,
     dindex beampipe_idx, const scalar brl_half_z,
     const std::vector<std::pair<scalar, scalar>> &lay_sizes,
     const std::vector<scalar> &lay_positions,
-    const std::vector<std::pair<int, int>> &m_binning, config_t cfg) {
+    const std::vector<std::pair<int, int>> &m_binning, config_t factory_cfg) {
 
     using nav_link_t = typename detector_t::surface_type::navigation_link;
     constexpr auto leaving_world{detail::invalid_value<nav_link_t>()};
 
+    const detail::detector_helper<transform3_t> det_helper{};
+
     // Generate consecutive linking between volumes
     dindex first_vol_idx = det.volumes().back().index();
     dindex last_vol_idx = first_vol_idx + 2u * n_layers;
@@ -889,25 +967,27 @@ inline void add_barrel_detector(
         vol_sizes.emplace_back(lay_sizes[i].first, lay_sizes[i].second);
     }
 
-    brl_module_factory m_factory{cfg};
+    brl_module_factory m_factory{factory_cfg};
     bool is_gap = true;
     for (unsigned int i = 0u; i < 2u * n_layers - 1u; ++i) {
         unsigned int j = (i + 2u) / 2u;
         // Every second layer is a gap volume
         is_gap = !is_gap;
         if (is_gap) {
-            detail::detector_helper<transform3_t>().create_cyl_volume(
-                det, resource, ctx, vol_sizes[i].first, vol_sizes[i].second,
-                -brl_half_z, brl_half_z, volume_links_vec[i]);
+            det_helper.create_cyl_volume(cfg, det, resource, ctx,
+                                         vol_sizes[i].first,
+                                         vol_sizes[i].second, -brl_half_z,
+                                         brl_half_z, volume_links_vec[i]);
 
             dindex vol_idx = det.volumes().back().index();
             names[vol_idx + 1u] = "gap_" + std::to_string(vol_idx);
         } else {
             m_factory.cfg.m_binning = m_binning[j];
             m_factory.cfg.layer_r = lay_positions[j];
-            detail::detector_helper<transform3_t>().create_cyl_volume(
-                det, resource, ctx, vol_sizes[i].first, vol_sizes[i].second,
-                -brl_half_z, brl_half_z, volume_links_vec[i]);
+            det_helper.create_cyl_volume(cfg, det, resource, ctx,
+                                         vol_sizes[i].first,
+                                         vol_sizes[i].second, -brl_half_z,
+                                         brl_half_z, volume_links_vec[i]);
 
             dindex vol_idx = det.volumes().back().index();
             names[vol_idx + 1u] = "barrel_" + std::to_string(vol_idx);
@@ -1076,7 +1156,7 @@ inline auto create_toy_geometry(vecmem::memory_resource &resource,
 
     // beampipe
     const dindex beampipe_idx{0u};
-    add_beampipe(det, resource, ctx0, name_map, cfg.n_edc_layers(),
+    add_beampipe(cfg, det, resource, ctx0, name_map, cfg.n_edc_layers(),
                  cfg.n_brl_layers(), edc_lay_sizes, brl_lay_sizes[0],
                  brl_positions[0], brl_half_z, edc_config.inner_r);
 
@@ -1084,8 +1164,8 @@ inline auto create_toy_geometry(vecmem::memory_resource &resource,
         edc_config.side = -1;
         // negative endcap layers
         add_endcap_detector<edc_module_factory>(
-            det, resource, ctx0, name_map, cfg.n_edc_layers(), beampipe_idx,
-            edc_lay_sizes, edc_positions, edc_config);
+            cfg, det, resource, ctx0, name_map, cfg.n_edc_layers(),
+            beampipe_idx, edc_lay_sizes, edc_positions, edc_config);
 
         // gap volume that connects barrel and neg. endcap
         dindex prev_vol_idx = det.volumes().back().index();
@@ -1095,15 +1175,17 @@ inline auto create_toy_geometry(vecmem::memory_resource &resource,
                                   : det.volumes().back().index() + 2u;
 
         add_endcap_barrel_connection(
-            det, resource, ctx0, name_map, edc_config.side, cfg.n_brl_layers(),
-            beampipe_idx, brl_lay_sizes, edc_config.inner_r, edc_config.outer_r,
-            edc_lay_sizes[0].first, brl_half_z, next_vol_idx, prev_vol_idx);
+            cfg, det, resource, ctx0, name_map, edc_config.side,
+            cfg.n_brl_layers(), beampipe_idx, brl_lay_sizes, edc_config.inner_r,
+            edc_config.outer_r, edc_lay_sizes[0].first, brl_half_z,
+            next_vol_idx, prev_vol_idx);
     }
     if (cfg.n_brl_layers() > 0u) {
         // barrel
         add_barrel_detector<brl_module_factory>(
-            det, resource, ctx0, name_map, cfg.n_brl_layers(), beampipe_idx,
-            brl_half_z, brl_lay_sizes, brl_positions, brl_binning, brl_config);
+            cfg, det, resource, ctx0, name_map, cfg.n_brl_layers(),
+            beampipe_idx, brl_half_z, brl_lay_sizes, brl_positions, brl_binning,
+            brl_config);
     }
     if (cfg.n_edc_layers() > 0u) {
         // gap layer that connects barrel and pos. endcap
@@ -1117,14 +1199,15 @@ inline auto create_toy_geometry(vecmem::memory_resource &resource,
                                   : det.volumes().back().index() + 2u;
 
         add_endcap_barrel_connection(
-            det, resource, ctx0, name_map, edc_config.side, cfg.n_brl_layers(),
-            beampipe_idx, brl_lay_sizes, edc_config.inner_r, edc_config.outer_r,
-            brl_half_z, edc_lay_sizes[0].first, prev_vol_idx, next_vol_idx);
+            cfg, det, resource, ctx0, name_map, edc_config.side,
+            cfg.n_brl_layers(), beampipe_idx, brl_lay_sizes, edc_config.inner_r,
+            edc_config.outer_r, brl_half_z, edc_lay_sizes[0].first,
+            prev_vol_idx, next_vol_idx);
 
         // positive endcap layers
         add_endcap_detector<edc_module_factory>(
-            det, resource, ctx0, name_map, cfg.n_edc_layers(), beampipe_idx,
-            edc_lay_sizes, edc_positions, edc_config);
+            cfg, det, resource, ctx0, name_map, cfg.n_edc_layers(),
+            beampipe_idx, edc_lay_sizes, edc_positions, edc_config);
     }
 
     // Add volume grid
diff --git a/utils/include/detray/detectors/create_wire_chamber.hpp b/utils/include/detray/detectors/create_wire_chamber.hpp
index 4be9038673..9d809cb884 100644
--- a/utils/include/detray/detectors/create_wire_chamber.hpp
+++ b/utils/include/detray/detectors/create_wire_chamber.hpp
@@ -15,6 +15,7 @@
 #include "detray/detectors/detector_helper.hpp"
 #include "detray/io/common/detail/file_handle.hpp"
 #include "detray/masks/masks.hpp"
+#include "detray/materials/mixture.hpp"
 #include "detray/materials/predefined_materials.hpp"
 #include "detray/tools/bounding_volume.hpp"
 #include "detray/tools/grid_builder.hpp"
@@ -50,6 +51,15 @@ struct wire_chamber_config {
     std::string m_bfield_file{!std::getenv("DETRAY_BFIELD_FILE")
                                   ? ""
                                   : std::getenv("DETRAY_BFIELD_FILE")};
+    /// Do material maps on portals
+    bool m_use_material_maps{false};
+    /// Number of bins for material maps
+    std::array<std::size_t, 2> m_cyl_map_bins{20u, 20u};
+    std::array<std::size_t, 2> m_disc_map_bins{3u, 20u};
+    /// Material to be filled into the maps
+    material<scalar> m_material =
+        mixture<scalar, silicon_tml<scalar, std::ratio<9, 10>>,
+                aluminium<scalar, std::ratio<1, 10>>>{};
 
     constexpr wire_chamber_config &n_layers(const unsigned int n) {
         m_n_layers = n;
@@ -76,11 +86,47 @@ struct wire_chamber_config {
         m_bfield_file = file_name;
         return *this;
     }
+    constexpr wire_chamber_config &use_material_maps(const bool b) {
+        m_use_material_maps = b;
+        return *this;
+    }
+    constexpr wire_chamber_config &cyl_map_bins(const std::size_t n_rphi,
+                                                const std::size_t n_z) {
+        m_cyl_map_bins = {n_rphi, n_z};
+        return *this;
+    }
+    constexpr wire_chamber_config &disc_map_bins(const std::size_t n_r,
+                                                 const std::size_t n_phi) {
+        m_disc_map_bins = {n_r, n_phi};
+        return *this;
+    }
+    constexpr wire_chamber_config &mapped_material(
+        const material<scalar> &mat) {
+        m_material = mat;
+        return *this;
+    }
 
     constexpr unsigned int n_layers() const { return m_n_layers; }
     constexpr scalar half_z() const { return m_half_z; }
     constexpr const vector3 &bfield_vec() const { return m_bfield_vec; }
     const std::string &bfield_file() const { return m_bfield_file; }
+    constexpr bool use_material_maps() const { return m_use_material_maps; }
+    constexpr const std::array<std::size_t, 2> &cyl_map_bins() const {
+        return m_cyl_map_bins;
+    }
+    constexpr const std::array<std::size_t, 2> &disc_map_bins() const {
+        return m_disc_map_bins;
+    }
+    scalar thickness() const { return 1.f * unit<scalar>::mm; }
+    std::vector<scalar> thickness(
+        const std::array<std::size_t, 2> &nbins) const {
+        return std::vector<scalar>(nbins[0] * nbins[1], 1.f * unit<scalar>::mm);
+    }
+    material<scalar> mapped_material() const { return m_material; }
+    std::vector<material<scalar>> mapped_material(
+        const std::array<std::size_t, 2> &nbins) const {
+        return std::vector<material<scalar>>(nbins[0] * nbins[1], m_material);
+    }
 
 };  // wire chamber config
 
@@ -136,7 +182,7 @@ auto create_wire_chamber(vecmem::memory_resource &resource,
 
     // Beam collision volume
     detail::detector_helper<transform3_t>().create_cyl_volume(
-        det, resource, ctx0, 0.f, inner_cyl_rad, -cyl_half_z, cyl_half_z,
+        cfg, det, resource, ctx0, 0.f, inner_cyl_rad, -cyl_half_z, cyl_half_z,
         {leaving_world, 1u, leaving_world, leaving_world});
 
     name_map[1u] = "beam_vol_0";
@@ -153,12 +199,12 @@ auto create_wire_chamber(vecmem::memory_resource &resource,
 
         if (i_lay < n_layers - 1) {
             detail::detector_helper<transform3_t>().create_cyl_volume(
-                det, resource, ctx0, inner_layer_rad, outer_layer_rad,
+                cfg, det, resource, ctx0, inner_layer_rad, outer_layer_rad,
                 -cyl_half_z, cyl_half_z,
                 {i_lay, i_lay + 2, leaving_world, leaving_world});
         } else {
             detail::detector_helper<transform3_t>().create_cyl_volume(
-                det, resource, ctx0, inner_layer_rad, outer_layer_rad,
+                cfg, det, resource, ctx0, inner_layer_rad, outer_layer_rad,
                 -cyl_half_z, cyl_half_z,
                 {i_lay, leaving_world, leaving_world, leaving_world});
         }
diff --git a/utils/include/detray/detectors/detector_helper.hpp b/utils/include/detray/detectors/detector_helper.hpp
index 7061ce96be..aae50c9e3e 100644
--- a/utils/include/detray/detectors/detector_helper.hpp
+++ b/utils/include/detray/detectors/detector_helper.hpp
@@ -11,8 +11,8 @@
 #include "detray/core/detail/data_context.hpp"
 #include "detray/definitions/geometry.hpp"
 #include "detray/definitions/indexing.hpp"
-#include "detray/materials/material.hpp"
-#include "detray/materials/predefined_materials.hpp"
+#include "detray/materials/material_map.hpp"
+#include "detray/materials/material_slab.hpp"
 
 namespace detray::detail {
 
@@ -37,22 +37,20 @@ struct detector_helper {
      * @param volume_link link to next volume for the masks
      */
     template <auto cyl_id, typename context_t, typename surface_container_t,
-              typename mask_container_t, typename material_container_t,
-              typename transform_container_t, typename volume_links>
-    inline void add_cylinder_surface(
-        const dindex volume_idx, context_t &ctx, surface_container_t &surfaces,
-        mask_container_t &masks, material_container_t &materials,
-        transform_container_t &transforms, const scalar_t r,
-        const scalar_t lower_z, const scalar_t upper_z,
-        const volume_links volume_link, const material<scalar_t> &mat,
-        const scalar_t thickness) {
+              typename mask_container_t, typename transform_container_t,
+              typename volume_links>
+    inline auto add_cylinder_surface(const dindex volume_idx, context_t &ctx,
+                                     surface_container_t &surfaces,
+                                     mask_container_t &masks,
+                                     transform_container_t &transforms,
+                                     const scalar_t r, const scalar_t lower_z,
+                                     const scalar_t upper_z,
+                                     const volume_links volume_link) const {
         using surface_type = typename surface_container_t::value_type;
         using mask_link_type = typename surface_type::mask_link;
         using material_id = typename surface_type::material_id;
         using material_link_type = typename surface_type::material_link;
 
-        constexpr auto slab_id = material_id::e_slab;
-
         const scalar_t min_z{std::min(lower_z, upper_z)};
         const scalar_t max_z{std::max(lower_z, upper_z)};
 
@@ -61,23 +59,21 @@ struct detector_helper {
 
         // add transform and masks
         transforms.emplace_back(ctx, tsl);
-        masks.template emplace_back<cyl_id>(empty_context{}, volume_link, r,
-                                            min_z, max_z);
-
-        // Add material slab
-        materials.template emplace_back<slab_id>(empty_context{}, mat,
-                                                 thickness);
+        auto &mask_ref = masks.template emplace_back<cyl_id>(
+            empty_context{}, volume_link, r, min_z, max_z);
 
         // add surface
         mask_link_type mask_link{cyl_id, masks.template size<cyl_id>() - 1u};
-        material_link_type material_link{
-            slab_id, materials.template size<slab_id>() - 1u};
+        material_link_type material_link{material_id::e_none, dindex_invalid};
         const surface_id sf_id = (volume_link != volume_idx)
                                      ? surface_id::e_portal
                                      : surface_id::e_passive;
 
-        surfaces.emplace_back(transforms.size(ctx) - 1u, mask_link,
-                              material_link, volume_idx, dindex_invalid, sf_id);
+        auto &sf_ref = surfaces.emplace_back(transforms.size(ctx) - 1u,
+                                             mask_link, material_link,
+                                             volume_idx, dindex_invalid, sf_id);
+
+        return std::tie(sf_ref, mask_ref);
     }
 
     /** Function that adds a disc portal.
@@ -89,21 +85,21 @@ struct detector_helper {
      * @param surfaces container to add new surface to
      * @param masks container to add new cylinder mask to
      * @param transforms container to add new transform to
-     * @param min_r lower radius of the disc
-     * @param max_r upper radius of the disc
+     * @param inner_r lower radius of the disc
+     * @param outer_r upper radius of the disc
      * @param z z position of the disc
      * @param volume_link link to next volume for the masks
      */
     template <typename context_t, typename surface_container_t,
-              typename mask_container_t, typename material_container_t,
-              typename transform_container_t, typename volume_links>
-    inline void add_disc_surface(
-        const dindex volume_idx, context_t &ctx, surface_container_t &surfaces,
-        mask_container_t &masks, material_container_t &materials,
-        transform_container_t &transforms, const scalar_t inner_r,
-        const scalar_t outer_r, const scalar_t z,
-        const volume_links volume_link, const material<scalar_t> &mat,
-        const scalar_t thickness) {
+              typename mask_container_t, typename transform_container_t,
+              typename volume_links>
+    inline auto add_disc_surface(const dindex volume_idx, context_t &ctx,
+                                 surface_container_t &surfaces,
+                                 mask_container_t &masks,
+                                 transform_container_t &transforms,
+                                 const scalar_t inner_r, const scalar_t outer_r,
+                                 const scalar_t z,
+                                 const volume_links volume_link) const {
         using surface_type = typename surface_container_t::value_type;
         using mask_id = typename surface_type::mask_id;
         using mask_link_type = typename surface_type::mask_link;
@@ -111,7 +107,6 @@ struct detector_helper {
         using material_link_type = typename surface_type::material_link;
 
         constexpr auto disc_id = mask_id::e_portal_ring2;
-        constexpr auto slab_id = material_id::e_slab;
 
         const scalar_t min_r{std::min(inner_r, outer_r)};
         const scalar_t max_r{std::max(inner_r, outer_r)};
@@ -121,22 +116,20 @@ struct detector_helper {
 
         // add transform and mask
         transforms.emplace_back(ctx, tsl);
-        masks.template emplace_back<disc_id>(empty_context{}, volume_link,
-                                             min_r, max_r);
-
-        // Add material slab
-        materials.template emplace_back<slab_id>(empty_context{}, mat,
-                                                 thickness);
+        auto &mask_ref = masks.template emplace_back<disc_id>(
+            empty_context{}, volume_link, min_r, max_r);
 
         // add surface
         mask_link_type mask_link{disc_id, masks.template size<disc_id>() - 1u};
-        material_link_type material_link{
-            slab_id, materials.template size<slab_id>() - 1u};
+        material_link_type material_link{material_id::e_none, dindex_invalid};
         const surface_id sf_id = (volume_link != volume_idx)
                                      ? surface_id::e_portal
                                      : surface_id::e_sensitive;
-        surfaces.emplace_back(transforms.size(ctx) - 1u, mask_link,
-                              material_link, volume_idx, dindex_invalid, sf_id);
+        auto &sf_ref = surfaces.emplace_back(transforms.size(ctx) - 1u,
+                                             mask_link, material_link,
+                                             volume_idx, dindex_invalid, sf_id);
+
+        return std::tie(sf_ref, mask_ref);
     }
 
     /** Function that adds a generic cylinder volume, using a factory for
@@ -145,6 +138,7 @@ struct detector_helper {
      * @tparam detector_t the detector type
      * @tparam factory_t type of module factory. Must be callable on containers.
      *
+     * @param cfg detector building configuration
      * @param det detector the volume should be added to
      * @param resource vecmem memory resource for the temporary containers
      * @param ctx geometry context
@@ -153,15 +147,15 @@ struct detector_helper {
      * @param lay_neg_r lower extend of volume
      * @param lay_pos_r upper extend of volume
      * @param volume_links volume links for the portals of the volume
-     * @param module_factory functor that adds module surfaces to volume
      */
-    template <typename detector_t>
-    void create_cyl_volume(detector_t &det, vecmem::memory_resource &resource,
+    template <typename config_t, typename detector_t>
+    void create_cyl_volume(const config_t &cfg, detector_t &det,
+                           vecmem::memory_resource &resource,
                            typename detector_t::geometry_context &ctx,
                            const scalar_t lay_inner_r,
                            const scalar_t lay_outer_r, const scalar_t lay_neg_z,
                            const scalar_t lay_pos_z,
-                           const std::vector<dindex> &volume_links) {
+                           const std::vector<dindex> &volume_links) const {
         // volume bounds
         const scalar_t inner_r{std::min(lay_inner_r, lay_outer_r)};
         const scalar_t outer_r{std::max(lay_inner_r, lay_outer_r)};
@@ -186,27 +180,91 @@ struct detector_helper {
         // negative and positive, inner and outer portal surface
         constexpr auto cyl_id = detector_t::masks::id::e_portal_cylinder2;
 
-        // If inner radius is 0, skip add the innder cylinder
+        auto &material_coll =
+            cfg.use_material_maps() ? det.material_store() : materials;
+
+        // If inner radius is 0, skip adding the inner cylinder
         if (inner_r > 0.f) {
-            add_cylinder_surface<cyl_id>(
-                cyl_volume.index(), ctx, surfaces, masks, materials, transforms,
-                inner_r, lower_z, upper_z, volume_links[0], vacuum<scalar_t>(),
-                0.f * unit<scalar_t>::mm);
+            auto [inner_cyl, inner_cyl_mask] = add_cylinder_surface<cyl_id>(
+                cyl_volume.index(), ctx, surfaces, masks, transforms, inner_r,
+                lower_z, upper_z, volume_links[0]);
+            create_material(cfg, inner_cyl, inner_cyl_mask, material_coll);
         }
-        add_cylinder_surface<cyl_id>(
-            cyl_volume.index(), ctx, surfaces, masks, materials, transforms,
-            outer_r, lower_z, upper_z, volume_links[1], vacuum<scalar_t>(),
-            0.f * unit<scalar_t>::mm);
-        add_disc_surface(cyl_volume.index(), ctx, surfaces, masks, materials,
-                         transforms, inner_r, outer_r, lower_z, volume_links[2],
-                         vacuum<scalar_t>(), 0.f * unit<scalar_t>::mm);
-        add_disc_surface(cyl_volume.index(), ctx, surfaces, masks, materials,
-                         transforms, inner_r, outer_r, upper_z, volume_links[3],
-                         vacuum<scalar_t>(), 0.f * unit<scalar_t>::mm);
+
+        auto [outer_cyl, outer_cyl_mask] = add_cylinder_surface<cyl_id>(
+            cyl_volume.index(), ctx, surfaces, masks, transforms, outer_r,
+            lower_z, upper_z, volume_links[1]);
+        create_material(cfg, outer_cyl, outer_cyl_mask, material_coll);
+
+        auto [neg_disc, neg_disc_mask] = add_disc_surface(
+            cyl_volume.index(), ctx, surfaces, masks, transforms, inner_r,
+            outer_r, lower_z, volume_links[2]);
+        create_material(cfg, neg_disc, neg_disc_mask, material_coll);
+
+        auto [pos_disc, pos_disc_mask] = add_disc_surface(
+            cyl_volume.index(), ctx, surfaces, masks, transforms, inner_r,
+            outer_r, upper_z, volume_links[3]);
+        create_material(cfg, pos_disc, pos_disc_mask, material_coll);
 
         det.add_objects_per_volume(ctx, cyl_volume, surfaces, masks, transforms,
                                    materials);
     }
+
+    template <typename config_t, typename surface_desc_t, typename mask_t,
+              typename material_container_t>
+    inline void create_material(const config_t &cfg, surface_desc_t &sf,
+                                const mask_t &sf_mask,
+                                material_container_t &materials) const {
+        using material_id = typename surface_desc_t::material_id;
+        using material_link_type = typename surface_desc_t::material_link;
+
+        if (cfg.use_material_maps()) {
+            material_map_factory<scalar_t> mat_map_factory{};
+
+            constexpr auto is_disc_map{
+                std::is_same_v<typename mask_t::shape, ring2D<>>};
+
+            const auto &nbins =
+                is_disc_map ? cfg.disc_map_bins() : cfg.cyl_map_bins();
+            const auto &mat = cfg.mapped_material(nbins);
+            const auto &t = cfg.thickness(nbins);
+
+            auto material_map =
+                mat_map_factory.template new_grid<>(sf_mask, nbins);
+
+            assert(mat.size() == material_map.nbins());
+            assert(t.size() == material_map.nbins());
+
+            for (dindex gbin = 0; gbin < material_map.nbins(); ++gbin) {
+                material_map.populate(
+                    gbin, material_slab<scalar_t>(mat[gbin], t[gbin]));
+            }
+
+            if constexpr (is_disc_map) {
+                materials.template push_back<material_id::e_disc2_map>(
+                    material_map);
+
+                sf.material() = material_link_type{
+                    material_id::e_disc2_map,
+                    materials.template size<material_id::e_disc2_map>() - 1u};
+            } else {
+                materials.template push_back<material_id::e_cylinder2_map>(
+                    material_map);
+
+                sf.material() = material_link_type{
+                    material_id::e_cylinder2_map,
+                    materials.template size<material_id::e_cylinder2_map>() -
+                        1u};
+            }
+        } else {
+            materials.template emplace_back<material_id::e_slab>(
+                {}, cfg.mapped_material(), cfg.thickness());
+
+            sf.material() = material_link_type{
+                material_id::e_slab,
+                materials.template size<material_id::e_slab>() - 1u};
+        }
+    }
 };
 
 }  // namespace detray::detail
diff --git a/utils/include/detray/detectors/itk_metadata.hpp b/utils/include/detray/detectors/itk_metadata.hpp
index 47e4987fce..62ca84745a 100644
--- a/utils/include/detray/detectors/itk_metadata.hpp
+++ b/utils/include/detray/detectors/itk_metadata.hpp
@@ -16,7 +16,7 @@
 #include "detray/intersection/cylinder_portal_intersector.hpp"
 #include "detray/io/common/detail/type_traits.hpp"  // mask_info
 #include "detray/masks/masks.hpp"
-#include "detray/materials/material_slab.hpp"
+#include "detray/materials/material_map.hpp"
 #include "detray/surface_finders/accelerator_grid.hpp"
 #include "detray/surface_finders/brute_force_finder.hpp"
 
@@ -55,9 +55,19 @@ struct itk_metadata {
     // Material Description
     //
 
-    /// The material types to be mapped onto the surfaces: Here homogeneous
-    /// material
-    using slab = material_slab<detray::scalar>;
+    // Cylindrical material map
+    template <typename container_t>
+    using cylinder_map_t =
+        material_map<cylinder2D<>::axes<>, scalar, container_t>;
+
+    // Disc material map
+    template <typename container_t>
+    using disc_map_t = material_map<ring2D<>::axes<>, scalar, container_t>;
+
+    // Rectangular material map
+    template <typename container_t>
+    using rectangular_map_t =
+        material_map<rectangle2D<>::axes<>, scalar, container_t>;
 
     /// How to store and link transforms. The geometry context allows to resolve
     /// the conditions data for e.g. module alignment
@@ -70,10 +80,10 @@ struct itk_metadata {
     /// order to minimize the depth of the 'unrolling' before a mask is found
     /// in the tuple
     enum class mask_ids {
-        e_rectangle2 = 0,
-        e_annulus2 = 1,
-        e_portal_cylinder2 = 2,
-        e_portal_ring2 = 3,
+        e_rectangle2 = 0u,
+        e_annulus2 = 1u,
+        e_portal_cylinder2 = 2u,
+        e_portal_ring2 = 3u,
     };
 
     /// This is the mask collections tuple (in the detector called 'mask store')
@@ -87,16 +97,22 @@ struct itk_metadata {
 
     /// Similar to the mask store, there is a material store, which
     enum class material_ids {
-        e_slab = 0,
-        e_none = 1,
+        e_disc2_map = 0u,
+        e_annulus2_map = 0u,
+        e_cylinder2_map = 1u,
+        e_reactangle2_map = 2u,
+        e_none = 3u,
     };
 
     /// How to store and link materials. The material does not make use of
     /// conditions data ( @c empty_context )
     template <template <typename...> class tuple_t = dtuple,
-              template <typename...> class vector_t = dvector>
-    using material_store = regular_multi_store<material_ids, empty_context,
-                                               tuple_t, vector_t, slab>;
+              typename container_t = host_container_types>
+    using material_store =
+        multi_store<material_ids, empty_context, tuple_t,
+                    grid_collection<disc_map_t<container_t>>,
+                    grid_collection<cylinder_map_t<container_t>>,
+                    grid_collection<rectangular_map_t<container_t>>>;
 
     /// Surface descriptor type used for sensitives, passives and portals
     /// It holds the indices to the surface data in the detector data stores
@@ -123,7 +139,7 @@ struct itk_metadata {
     /// The acceleration data structures live in another tuple that needs to be
     /// indexed correctly:
     enum class sf_finder_ids {
-        e_brute_force = 0,  //< test all surfaces in a volume (brute force)
+        e_brute_force = 0u,  //< test all surfaces in a volume (brute force)
         e_default = e_brute_force,
     };
 
diff --git a/utils/include/detray/detectors/telescope_metadata.hpp b/utils/include/detray/detectors/telescope_metadata.hpp
index c761842735..2269b0e311 100644
--- a/utils/include/detray/detectors/telescope_metadata.hpp
+++ b/utils/include/detray/detectors/telescope_metadata.hpp
@@ -86,14 +86,14 @@ struct telescope_metadata {
 
     /// How to store materials
     template <template <typename...> class tuple_t = dtuple,
-              template <typename...> class vector_t = dvector>
+              typename container_t = host_container_types>
     using material_store = std::conditional_t<
         std::is_same_v<mask<mask_shape_t, nav_link>, cell_wire> |
             std::is_same_v<mask<mask_shape_t, nav_link>, straw_wire>,
-        regular_multi_store<material_ids, empty_context, tuple_t, vector_t,
-                            slab, rod>,
-        regular_multi_store<material_ids, empty_context, tuple_t, vector_t,
-                            slab>>;
+        regular_multi_store<material_ids, empty_context, tuple_t,
+                            container_t::template vector_type, slab, rod>,
+        regular_multi_store<material_ids, empty_context, tuple_t,
+                            container_t::template vector_type, slab>>;
 
     /// How to link to the entries in the data stores
     using transform_link = typename transform_store<>::link_type;
diff --git a/utils/include/detray/detectors/toy_metadata.hpp b/utils/include/detray/detectors/toy_metadata.hpp
index 8858536ca2..a73f3b6300 100644
--- a/utils/include/detray/detectors/toy_metadata.hpp
+++ b/utils/include/detray/detectors/toy_metadata.hpp
@@ -17,6 +17,7 @@
 #include "detray/intersection/cylinder_portal_intersector.hpp"
 #include "detray/intersection/plane_intersector.hpp"
 #include "detray/masks/masks.hpp"
+#include "detray/materials/material_map.hpp"
 #include "detray/materials/material_slab.hpp"
 #include "detray/surface_finders/accelerator_grid.hpp"
 #include "detray/surface_finders/brute_force_finder.hpp"
@@ -40,10 +41,24 @@ struct toy_metadata {
     /// Material types
     using slab = material_slab<detray::scalar>;
 
-    /// surface grid types (regular, open binning)
+    // Cylindrical material grid
+    template <typename container_t>
+    using cylinder_map_t =
+        material_map<cylinder2D<>::axes<>, scalar, container_t>;
+
+    // Disc material grid
+    template <typename container_t>
+    using disc_map_t = material_map<ring2D<>::axes<>, scalar, container_t>;
+
+    // Rectangular material grid
+    template <typename container_t>
+    using rectangular_map_t =
+        material_map<rectangle2D<>::axes<>, scalar, container_t>;
+
+    /// Surface grid types (regular, open binning)
     /// @{
 
-    // surface grid definition: bin-content: std::array<dindex, 9>
+    // Surface grid definition: bin-content: std::array<sf_descriptor, 1>
     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,
@@ -84,15 +99,22 @@ struct toy_metadata {
 
     /// Material type ids
     enum class material_ids {
-        e_slab = 0,
-        e_none = 1,
+        e_disc2_map = 0u,
+        e_cylinder2_map = 1u,
+        e_reactangle2_map = 2u,
+        e_slab = 3u,
+        e_none = 4u,
     };
 
     /// How to store materials
     template <template <typename...> class tuple_t = dtuple,
-              template <typename...> class vector_t = dvector>
-    using material_store = regular_multi_store<material_ids, empty_context,
-                                               tuple_t, vector_t, slab>;
+              typename container_t = host_container_types>
+    using material_store =
+        multi_store<material_ids, empty_context, tuple_t,
+                    grid_collection<disc_map_t<container_t>>,
+                    grid_collection<cylinder_map_t<container_t>>,
+                    grid_collection<rectangular_map_t<container_t>>,
+                    typename container_t::template vector_type<slab>>;
 
     /// How to link to the entries in the data stores
     using transform_link = typename transform_store<>::link_type;