Core/include/Acts/Detector/TrackingGeometry.hpp

@@ -13,6 +13,7 @@
 #pragma once
 #include "Acts/Utilities/Definitions.hpp"
 #include "Acts/Utilities/GeometryContext.hpp"
+#include "Acts/Utilities/GeometryID.hpp"
 #include "Acts/Utilities/GeometrySignature.hpp"
 
 #include <functional>
@@ -26,9 +27,12 @@ class TrackingVolume;
 class Layer;
 class Surface;
 class PerigeeSurface;
+class SurfaceMaterial;
 
 using TrackingVolumePtr        = std::shared_ptr<const TrackingVolume>;
 using MutableTrackingVolumePtr = std::shared_ptr<TrackingVolume>;
+using SurfaceMaterialMap
+    = std::map<GeometryID, std::shared_ptr<const SurfaceMaterial>>;
 
 ///  @class TrackingGeometry
 ///
@@ -48,7 +52,11 @@ public:
   /// Constructor
   ///
   /// @param highestVolume is the world volume
-  TrackingGeometry(const MutableTrackingVolumePtr& highestVolume);
+  /// @param surfaceMaterialMap is a map that holds optionally the
+  ///        the surface material to be assigned by GeometryID
+  TrackingGeometry(const MutableTrackingVolumePtr& highestVolume,
+                   const SurfaceMaterialMap&       surfaceMaterialMap
+                   = SurfaceMaterialMap());
 
   /// Destructor
   ~TrackingGeometry();

Core/include/Acts/Detector/TrackingVolume.hpp

@@ -21,6 +21,7 @@
 #include "Acts/Utilities/BinnedArray.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 #include "Acts/Utilities/GeometryContext.hpp"
+#include "Acts/Utilities/GeometryID.hpp"
 #include "Acts/Utilities/GeometrySignature.hpp"
 #include "Acts/Volumes/BoundarySurfaceT.hpp"
 #include "Acts/Volumes/Volume.hpp"
@@ -30,6 +31,7 @@ namespace Acts {
 class GlueVolumesDescriptor;
 class VolumeBounds;
 class Material;
+class SurfaceMaterial;
 
 using TrackingVolumeBoundaryPtr
     = std::shared_ptr<const BoundarySurfaceT<TrackingVolume>>;
@@ -47,6 +49,10 @@ using LayerVector          = std::vector<LayerPtr>;
 // full intersection with Layer
 using LayerIntersection = FullIntersection<Layer, Surface>;
 
+// SurfaceMaterial assignment map
+using SurfaceMaterialMap
+    = std::map<GeometryID, std::shared_ptr<const SurfaceMaterial>>;
+
 // full intersection with surface
 using BoundaryIntersection
     = FullIntersection<BoundarySurfaceT<TrackingVolume>, Surface>;
@@ -411,15 +417,18 @@ private:
   void
   synchronizeLayers(double envelope = 1.) const;
 
-  /// close the Geometry, i.e. set the TDD_ID
+  /// close the Geometry, i.e. set the GeometryID and assign material
   ///
+  /// @param surfaceMaterialMap is a map that contains
+  ///        the surface material mapped iwth the GeometryID
   /// @param volumeMap is a map to find the a volume
   ///        by a given name
   /// @param vol is the geometry id of the volume
   ///        as calculated by the TrackingGeometry
   ///
   void
-  closeGeometry(std::map<std::string, const TrackingVolume*>& volumeMap,
+  closeGeometry(const SurfaceMaterialMap& surfaceMaterialMap,
+                std::map<std::string, const TrackingVolume*>& volumeMap,
                 size_t& vol);
 
   /// interlink the layers in this TrackingVolume

Core/include/Acts/Extrapolator/MaterialInteractor.hpp

@@ -321,4 +321,13 @@ private:
   }
 };
 
+/// Using some short hands for Recorded Material
+using RecordedMaterial = MaterialInteractor::result_type;
+
+/// And recorded material track
+/// - this is start:  position, start momentum
+///   and the Recorded material
+using RecordedMaterialTrack
+    = std::pair<std::pair<Acts::Vector3D, Acts::Vector3D>, RecordedMaterial>;
+
 }  // end of namespace Acts

Core/include/Acts/Layers/Layer.hpp

@@ -21,6 +21,7 @@
 #include "Acts/Utilities/BinnedArray.hpp"
 #include "Acts/Utilities/Definitions.hpp"
 #include "Acts/Utilities/GeometryContext.hpp"
+#include "Acts/Utilities/GeometryID.hpp"
 #include "Acts/Utilities/GeometryObject.hpp"
 #include "Acts/Utilities/GeometryStatics.hpp"
 #include "Acts/Utilities/Intersection.hpp"
@@ -29,6 +30,7 @@
 namespace Acts {
 
 class Surface;
+class SurfaceMaterial;
 class BinUtility;
 class Volume;
 class VolumeBounds;
@@ -38,6 +40,10 @@ class ApproachDescriptor;
 // Simple surface intersection
 using SurfaceIntersection = ObjectIntersection<Surface>;
 
+// SurfaceMaterial assignment
+using SurfaceMaterialMap
+    = std::map<GeometryID, std::shared_ptr<const SurfaceMaterial>>;
+
 // master typedef
 class Layer;
 using LayerPtr        = std::shared_ptr<const Layer>;
@@ -370,12 +376,17 @@ protected:
 
 private:
   /// Private helper method to close the geometry
+  /// - it will assign material to the surfaces if needed
   /// - it will set the layer geometry ID for a unique identification
   /// - it will also register the internal sub strucutre
+  ///
+  /// @param surfaceMaterialMap is a map that contains the surface
+  ///        the surface material
   /// @param layerID is the geometry id of the volume
   ///                as calculated by the TrackingGeometry
   void
-  closeGeometry(const GeometryID& layerID);
+  closeGeometry(const SurfaceMaterialMap& surfaceMaterialMap,
+                const GeometryID&         layerID);
 };
 
 /// Layers are constructedd with shared_ptr factories, hence the layer array is

Core/include/Acts/MagneticField/InterpolatedBFieldMap.hpp

@@ -160,10 +160,9 @@ public:
   getFieldCell(const Vector3D& position) const
   {
     const auto& gridPosition = m_transformPos(position);
-    size_t      bin          = m_grid.getGlobalBinIndex(gridPosition);
-    const auto& indices      = m_grid.getLocalBinIndices(bin);
-    const auto& lowerLeft    = m_grid.getLowerLeftBinEdge(indices);
-    const auto& upperRight   = m_grid.getUpperRightBinEdge(indices);
+    const auto& indices      = m_grid.localBinsFromPosition(gridPosition);
+    const auto& lowerLeft    = m_grid.lowerLeftBinEdge(indices);
+    const auto& upperRight   = m_grid.upperRightBinEdge(indices);
 
     // loop through all corner points
     constexpr size_t nCorners = 1 << DIM_POS;
@@ -185,7 +184,7 @@ public:
   std::vector<size_t>
   getNBins() const
   {
-    auto nBinsArray = m_grid.getNBins();
+    auto nBinsArray = m_grid.numLocalBins();
     return std::vector<size_t>(nBinsArray.begin(), nBinsArray.end());
   }
 
@@ -195,7 +194,7 @@ public:
   std::vector<double>
   getMin() const
   {
-    auto minArray = m_grid.getMin();
+    auto minArray = m_grid.minPosition();
     return std::vector<double>(minArray.begin(), minArray.end());
   }
 
@@ -205,7 +204,7 @@ public:
   std::vector<double>
   getMax() const
   {
-    auto maxArray = m_grid.getMax();
+    auto maxArray = m_grid.maxPosition();
     return std::vector<double>(maxArray.begin(), maxArray.end());
   }
 

Core/include/Acts/Surfaces/Surface.hpp

@@ -277,7 +277,7 @@ public:
   /// @param material Material description this given and stored as a shared
   /// pointer
   void
-  setAssociatedMaterial(std::shared_ptr<const SurfaceMaterial> material);
+  assignSurfaceMaterial(std::shared_ptr<const SurfaceMaterial> material);
 
   /// The templated Parameters onSurface method
   /// In order to avoid unneccessary geometrical operations, it checks on the

Core/include/Acts/Surfaces/SurfaceArray.hpp

@@ -231,7 +231,7 @@ public:
     SurfaceVector&
     lookup(const Vector3D& pos) override
     {
-      return m_grid.at(m_globalToLocal(pos));
+      return m_grid.atPosition(m_globalToLocal(pos));
     }
 
     /// @brief Performs lookup at @c pos and returns bin content as const
@@ -241,7 +241,7 @@ public:
     const SurfaceVector&
     lookup(const Vector3D& pos) const override
     {
-      return m_grid.at(m_globalToLocal(pos));
+      return m_grid.atPosition(m_globalToLocal(pos));
     }
 
     /// @brief Performs lookup at global bin and returns bin content as
@@ -272,7 +272,7 @@ public:
     neighbors(const Vector3D& pos) const override
     {
       auto loc = m_globalToLocal(pos);
-      return m_neighborMap.at(m_grid.getGlobalBinIndex(loc));
+      return m_neighborMap.at(m_grid.globalBinFromPosition(loc));
     }
 
     /// @brief Returns the total size of the grid (including under/overflow
@@ -299,7 +299,7 @@ public:
     std::vector<const IAxis*>
     getAxes() const override
     {
-      auto arr = m_grid.getAxes();
+      auto arr = m_grid.axes();
       return std::vector<const IAxis*>(arr.begin(), arr.end());
     }
 
@@ -319,8 +319,8 @@ public:
     bool
     isValidBin(size_t bin) const override
     {
-      std::array<size_t, DIM> indices = m_grid.getLocalBinIndices(bin);
-      std::array<size_t, DIM> nBins   = m_grid.getNBins();
+      std::array<size_t, DIM> indices = m_grid.localBinsFromGlobalBin(bin);
+      std::array<size_t, DIM> nBins   = m_grid.numLocalBins();
       for (size_t i = 0; i < indices.size(); ++i) {
         size_t idx = indices.at(i);
         if (idx <= 0 || idx >= nBins.at(i) + 1) {
@@ -340,8 +340,8 @@ public:
         if (!isValidBin(i)) {
           continue;
         }
-        typename Grid_t::index_t loc  = m_grid.getLocalBinIndices(i);
-        std::set<size_t> neighborIdxs = m_grid.neighborHoodIndices(loc, 1u);
+        typename Grid_t::index_t loc = m_grid.localBinsFromGlobalBin(i);
+        auto neighborIdxs            = m_grid.neighborHoodIndices(loc, 1u);
         std::vector<const Surface*>& neighbors = m_neighborMap.at(i);
         neighbors.clear();
 
@@ -369,7 +369,7 @@ public:
     getBinCenterImpl(size_t bin) const
     {
       return m_localToGlobal(ActsVectorD<DIM>(
-          m_grid.getBinCenter(m_grid.getLocalBinIndices(bin)).data()));
+          m_grid.binCenter(m_grid.localBinsFromGlobalBin(bin)).data()));
     }
 
     /// Internal method, see above.
@@ -378,7 +378,7 @@ public:
     Vector3D
     getBinCenterImpl(size_t bin) const
     {
-      point_t pos = m_grid.getBinCenter(m_grid.getLocalBinIndices(bin));
+      point_t pos = m_grid.binCenter(m_grid.localBinsFromGlobalBin(bin));
       return m_localToGlobal(pos);
     }
 

Core/include/Acts/Surfaces/detail/Surface.ipp

@@ -159,7 +159,7 @@ Surface::associatedMaterial() const
 }
 
 inline void
-Surface::setAssociatedMaterial(std::shared_ptr<const SurfaceMaterial> material)
+Surface::assignSurfaceMaterial(std::shared_ptr<const SurfaceMaterial> material)
 {
   m_associatedMaterial = std::move(material);
 }

Core/include/Acts/Utilities/BinAdjustment.hpp

+// This file is part of the Acts project.
+//
+// Copyright (C) 2019 Acts project team
+//
+// This Source Code Form is subject to the terms of the Mozilla Public
+// License, v. 2.0. If a copy of the MPL was not distributed with this
+// file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+///////////////////////////////////////////////////////////////////
+// BinAdjustment.hpp, Acts project
+///////////////////////////////////////////////////////////////////
+
+#pragma once
+
+#include <stdexcept>
+
+#include "Acts/Surfaces/CylinderBounds.hpp"
+#include "Acts/Surfaces/RadialBounds.hpp"
+#include "Acts/Surfaces/Surface.hpp"
+#include "Acts/Utilities/BinUtility.hpp"
+
+namespace Acts {
+
+/// @brief adjust the BinUtility bu to the dimensions of radial bounds
+///
+/// @param bu BinUtility at source
+/// @param rBounds the Radial bounds to adjust to
+///
+/// @return new updated BinUtiltiy
+BinUtility
+adjustBinUtility(const BinUtility& bu, const RadialBounds& rBounds)
+{
+  // Default constructor
+  BinUtility uBinUtil;
+  // The parameters from the cylinder bounds
+  double minR   = rBounds.rMin();
+  double maxR   = rBounds.rMax();
+  double minPhi = rBounds.averagePhi() - rBounds.halfPhiSector();
+  double maxPhi = rBounds.averagePhi() + rBounds.halfPhiSector();
+  // Retrieve the binning data
+  const std::vector<BinningData>& bData = bu.binningData();
+  // Loop over the binning data and adjust the dimensions
+  for (auto& bd : bData) {
+    // The binning value
+    BinningValue bval = bd.binvalue;
+    // Throw exceptions is stuff doesn't make sense:
+    // - not the right binning value
+    // - not equidistant
+    if (bd.type == arbitrary) {
+      throw std::invalid_argument("Arbirary binning can not be adjusted.");
+    } else if (bval != binR and bval != binPhi) {
+      throw std::invalid_argument("Disc binning must be: phi, r");
+    }
+    float min, max = 0.;
+    // Perform the value adjustment
+    if (bval == binPhi) {
+      min = minPhi;
+      max = maxPhi;
+    } else {
+      min = minR;
+      max = maxR;
+    }
+    // Create the updated BinningData
+    BinningData uBinData(bd.option, bval, bd.bins(), min, max);
+    uBinUtil += BinUtility(uBinData);
+  }
+  return uBinUtil;
+}
+
+/// @brief adjust the BinUtility bu to the dimensions of cylinder bounds
+///
+/// @param bu BinUtility at source
+/// @param cBounds the Cylinder bounds to adjust to
+///
+/// @return new updated BinUtiltiy
+BinUtility
+adjustBinUtility(const BinUtility& bu, const CylinderBounds& cBounds)
+{
+  // Default constructor
+  BinUtility uBinUtil;
+  // The parameters from the cylinder bounds
+  double cR     = cBounds.r();
+  double cHz    = cBounds.halflengthZ();
+  double minPhi = cBounds.averagePhi() - cBounds.halfPhiSector();
+  double maxPhi = cBounds.averagePhi() + cBounds.halfPhiSector();
+  // Retrieve the binning data
+  const std::vector<BinningData>& bData = bu.binningData();
+  // Loop over the binning data and adjust the dimensions
+  for (auto& bd : bData) {
+    // The binning value
+    BinningValue bval = bd.binvalue;
+    // Throw exceptions if stuff doesn't make sense:
+    // - not the right binning value
+    // - not equidistant
+    if (bd.type == arbitrary) {
+      throw std::invalid_argument("Arbitrary binning can not be adjusted.");
+    } else if (bval != binRPhi and bval != binPhi and bval != binZ) {
+      throw std::invalid_argument("Cylinder binning must be: rphi, phi, z");
+    }
+    float min, max = 0.;
+    // Perform the value adjustment
+    if (bval == binPhi) {
+      min = minPhi;
+      max = maxPhi;
+    } else if (bval == binRPhi) {
+      min = cR * minPhi;
+      max = cR * maxPhi;
+    } else {
+      min = -cHz;
+      max = cHz;
+    }
+    // Create the updated BinningData
+    BinningData uBinData(bd.option, bval, bd.bins(), min, max);
+    uBinUtil += BinUtility(uBinData);
+  }
+  return uBinUtil;
+}
+
+/// @brief adjust the BinUtility bu to a surface
+///
+/// @param bu BinUtility at source
+/// @param Surface to which the adjustment is being done
+///
+/// @return new updated BinUtiltiy
+BinUtility
+adjustBinUtility(const BinUtility& bu, const Surface& surface)
+{
+  // The surface type is a cylinder
+  if (surface.type() == Surface::Cylinder) {
+    // Cast to Cylinder bounds and return
+    auto cBounds = dynamic_cast<const CylinderBounds*>(&(surface.bounds()));
+    // Return specific adjustment
+    return adjustBinUtility(bu, *cBounds);
+
+  } else if (surface.type() == Surface::Disc) {
+    // Cast to Cylinder bounds and return
+    auto rBounds = dynamic_cast<const RadialBounds*>(&(surface.bounds()));
+    // Return specific adjustment
+    return adjustBinUtility(bu, *rBounds);
+  }
+
+  throw std::invalid_argument(
+      "Bin adjustment not implemented for this surface yet!");
+
+  return BinUtility();
+}
+
+}  // namespace Acts

Core/include/Acts/Utilities/GeometryID.hpp

@@ -104,6 +104,15 @@ public:
     return (m_value == tddID.value());
   }
 
+  /// Non-equality operator
+  ///
+  /// @param tddID is the geometry ID that will be compared on equality
+  bool
+  operator!=(const GeometryID& tddID) const
+  {
+    return !operator==(tddID);
+  }
+
   /// Add some stuff
   ///
   /// @param type_id which identifier do you wanna add

Core/include/Acts/Utilities/detail/Axis.hpp

@@ -10,7 +10,6 @@
 
 #include <algorithm>
 #include <cmath>
-#include <set>
 #include <vector>
 #include "Acts/Utilities/IAxis.hpp"
 
@@ -31,6 +30,111 @@ namespace detail {
   /// Enum which determines the binning type of the axis
   enum class AxisType { Equidistant, Variable };
 
+  // This object can be iterated to produce up to two sequences of integer
+  // indices, corresponding to the half-open integer ranges [begin1, end1[ and
+  // [begin2, end2[.
+  //
+  // The goal is to emulate the effect of enumerating a range of neighbor
+  // indices on an axis (which may go out of bounds and wrap around since we
+  // have AxisBoundaryType::Closed), inserting them into an std::vector, and
+  // discarding duplicates, without paying the price of duplicate removal
+  // and dynamic memory allocation in hot magnetic field interpolation code.
+  //
+  class NeighborHoodIndices
+  {
+  public:
+    NeighborHoodIndices() = default;
+
+    NeighborHoodIndices(size_t begin, size_t end)
+      : m_begin1(begin), m_end1(end), m_begin2(end), m_end2(end)
+    {
+    }
+
+    NeighborHoodIndices(size_t begin1, size_t end1, size_t begin2, size_t end2)
+      : m_begin1(begin1), m_end1(end1), m_begin2(begin2), m_end2(end2)
+    {
+    }
+
+    class iterator
+    {
+    public:
+      iterator() = default;
+
+      // Specialized constructor for end() iterator
+      iterator(size_t current) : m_current(current), m_wrapped(true) {}
+
+      iterator(size_t begin1, size_t end1, size_t begin2)
+        : m_current(begin1)
+        , m_end1(end1)
+        , m_begin2(begin2)
+        , m_wrapped(begin1 == begin2)
+      {
+      }
+
+      size_t operator*() const { return m_current; }
+
+      iterator& operator++()
+      {
+        ++m_current;
+        if (m_current == m_end1) {
+          m_current = m_begin2;
+          m_wrapped = true;
+        }
+        return *this;
+      }
+
+      bool
+      operator==(const iterator& it) const
+      {
+        return (m_current == it.m_current) && (m_wrapped == it.m_wrapped);
+      }
+
+      bool
+      operator!=(const iterator& it) const
+      {
+        return !(*this == it);
+      }
+
+    private:
+      size_t m_current, m_end1, m_begin2;
+      bool   m_wrapped;
+    };
+
+    iterator
+    begin() const
+    {
+      return iterator(m_begin1, m_end1, m_begin2);
+    }
+
+    iterator
+    end() const
+    {
+      return iterator(m_end2);
+    }
+
+    // Number of indices that will be produced if this sequence is iterated
+    size_t
+    size() const
+    {
+      return (m_end1 - m_begin1) + (m_end2 - m_begin2);
+    }
+
+    // Collect the sequence of indices into an std::vector
+    std::vector<size_t>
+    collect() const
+    {
+      std::vector<size_t> result;
+      result.reserve(this->size());
+      for (size_t idx : *this) {
+        result.push_back(idx);
+      }
+      return result;
+    }
+
+  private:
+    size_t m_begin1 = 0, m_end1 = 0, m_begin2 = 0, m_end2 = 0;
+  };
+
   /// @brief calculate bin indices from a given binning structure
   ///
   /// This class provides some basic functionality for calculating bin indices
@@ -102,8 +206,8 @@ namespace detail {
     ///
     /// @param [in] idx requested bin index
     /// @param [in] sizes how many neighboring bins (up/down)
-    /// @return std::set of neighboring bin indices (global)
-    std::set<size_t>
+    /// @return Set of neighboring bin indices (global)
+    NeighborHoodIndices
     neighborHoodIndices(size_t idx, size_t size = 1) const
     {
       return neighborHoodIndices(idx, std::make_pair(size, size));
@@ -115,25 +219,21 @@ namespace detail {
     ///
     /// @param [in] idx requested bin index
     /// @param [in] sizes how many neighboring bins (up/down)
-    /// @return std::set of neighboring bin indices (global)
+    /// @return Set of neighboring bin indices (global)
     /// @note Open varies given bin and allows 0 and NBins+1 (underflow,
     /// overflow)
     ///       as neighbors
     template <AxisBoundaryType T = bdt,
               std::enable_if_t<T == AxisBoundaryType::Open, int> = 0>
-    std::set<size_t>
+    NeighborHoodIndices
     neighborHoodIndices(size_t idx,
                         std::pair<size_t, size_t> sizes = {1, 1}) const
     {
-      std::set<size_t> result;
-      constexpr int    min   = 0;
-      const int        max   = getNBins() + 1;
-      const int        itmin = std::max(min, int(idx - sizes.first));
-      const int        itmax = std::min(max, int(idx + sizes.second));
-      for (int i = itmin; i <= itmax; i++) {
-        result.insert(i);
-      }
-      return result;
+      constexpr int min   = 0;
+      const int     max   = getNBins() + 1;
+      const int     itmin = std::max(min, int(idx - sizes.first));
+      const int     itmax = std::min(max, int(idx + sizes.second));
+      return NeighborHoodIndices(itmin, itmax + 1);
     }
 
     /// @brief Get #size bins which neighbor the one given
@@ -142,27 +242,23 @@ namespace detail {
     ///
     /// @param [in] idx requested bin index
     /// @param [in] sizes how many neighboring bins (up/down)
-    /// @return std::set of neighboring bin indices (global)
+    /// @return Set of neighboring bin indices (global)
     /// @note Bound varies given bin and allows 1 and NBins (regular bins)
     ///       as neighbors
     template <AxisBoundaryType T = bdt,
               std::enable_if_t<T == AxisBoundaryType::Bound, int> = 0>
-    std::set<size_t>
+    NeighborHoodIndices
     neighborHoodIndices(size_t idx,
                         std::pair<size_t, size_t> sizes = {1, 1}) const
     {
-      std::set<size_t> result;
       if (idx <= 0 || idx >= (getNBins() + 1)) {
-        return result;
+        return NeighborHoodIndices();
       }
       constexpr int min   = 1;
       const int     max   = getNBins();
       const int     itmin = std::max(min, int(idx - sizes.first));
       const int     itmax = std::min(max, int(idx + sizes.second));
-      for (int i = itmin; i <= itmax; i++) {
-        result.insert(i);
-      }
-      return result;
+      return NeighborHoodIndices(itmin, itmax + 1);
     }
 
     /// @brief Get #size bins which neighbor the one given
@@ -171,25 +267,44 @@ namespace detail {
     ///
     /// @param [in] idx requested bin index
     /// @param [in] sizes how many neighboring bins (up/down)
-    /// @return std::set of neighboring bin indices (global)
+    /// @return Set of neighboring bin indices (global)
     /// @note Closed varies given bin and allows bins on the opposite
     ///       side of the axis as neighbors. (excludes underflow / overflow)
     template <AxisBoundaryType T = bdt,
               std::enable_if_t<T == AxisBoundaryType::Closed, int> = 0>
-    std::set<size_t>
+    NeighborHoodIndices
     neighborHoodIndices(size_t idx,
                         std::pair<size_t, size_t> sizes = {1, 1}) const
     {
-      std::set<size_t> result;
+      // Handle invalid indices
       if (idx <= 0 || idx >= (getNBins() + 1)) {
-        return result;
+        return NeighborHoodIndices();
       }
-      const int itmin = idx - sizes.first;
-      const int itmax = idx + sizes.second;
-      for (int i = itmin; i <= itmax; i++) {
-        result.insert(wrapBin(i));
+
+      // Handle corner case where user requests more neighbours than the number
+      // of bins on the axis. We do not want to double-count bins in that case.
+      sizes.first %= getNBins();
+      sizes.second %= getNBins();
+      if (sizes.first + sizes.second + 1 > getNBins()) {
+        sizes.second -= (sizes.first + sizes.second + 1) - getNBins();
+      }
+
+      // If the entire index range is not covered, we must wrap the range of
+      // targeted neighbor indices into the range of valid bin indices. This may
+      // split the range of neighbor indices in two parts:
+      //
+      // Before wraparound - [        XXXXX]XXX
+      // After wraparound  - [ XXXX   XXXX ]
+      //
+      const int    itmin   = idx - sizes.first;
+      const int    itmax   = idx + sizes.second;
+      const size_t itfirst = wrapBin(itmin);
+      const size_t itlast  = wrapBin(itmax);
+      if (itfirst <= itlast) {
+        return NeighborHoodIndices(itfirst, itlast + 1);
+      } else {
+        return NeighborHoodIndices(itfirst, getNBins() + 1, 1, itlast + 1);
       }
-      return result;
     }
 
     /// @brief Converts bin index into a valid one for this axis.
@@ -418,8 +533,8 @@ namespace detail {
     ///
     /// @param [in] idx requested bin index
     /// @param [in] size how many neighboring bins
-    /// @return std::set of neighboring bin indices (global)
-    std::set<size_t>
+    /// @return Set of neighboring bin indices (global)
+    NeighborHoodIndices
     neighborHoodIndices(size_t idx, size_t size = 1) const
     {
       return neighborHoodIndices(idx, std::make_pair(size, size));
@@ -431,25 +546,21 @@ namespace detail {
     ///
     /// @param [in] idx requested bin index
     /// @param [in] sizes how many neighboring bins (up/down)
-    /// @return std::set of neighboring bin indices (global)
+    /// @return Set of neighboring bin indices (global)
     /// @note Open varies given bin and allows 0 and NBins+1 (underflow,
     /// overflow)
     ///       as neighbors
     template <AxisBoundaryType T = bdt,
               std::enable_if_t<T == AxisBoundaryType::Open, int> = 0>
-    std::set<size_t>
+    NeighborHoodIndices
     neighborHoodIndices(size_t idx,
                         std::pair<size_t, size_t> sizes = {1, 1}) const
     {
-      std::set<size_t> result;
-      constexpr int    min   = 0;
-      const int        max   = getNBins() + 1;
-      const int        itmin = std::max(min, int(idx - sizes.first));
-      const int        itmax = std::min(max, int(idx + sizes.second));
-      for (int i = itmin; i <= itmax; i++) {
-        result.insert(i);
-      }
-      return result;
+      constexpr int min   = 0;
+      const int     max   = getNBins() + 1;
+      const int     itmin = std::max(min, int(idx - sizes.first));
+      const int     itmax = std::min(max, int(idx + sizes.second));
+      return NeighborHoodIndices(itmin, itmax + 1);
     }
 
     /// @brief Get #size bins which neighbor the one given
@@ -458,27 +569,23 @@ namespace detail {
     ///
     /// @param [in] idx requested bin index
     /// @param [in] sizes how many neighboring bins (up/down)
-    /// @return std::set of neighboring bin indices (global)
+    /// @return Set of neighboring bin indices (global)
     /// @note Bound varies given bin and allows 1 and NBins (regular bins)
     ///       as neighbors
     template <AxisBoundaryType T = bdt,
               std::enable_if_t<T == AxisBoundaryType::Bound, int> = 0>
-    std::set<size_t>
+    NeighborHoodIndices
     neighborHoodIndices(size_t idx,
                         std::pair<size_t, size_t> sizes = {1, 1}) const
     {
-      std::set<size_t> result;
       if (idx <= 0 || idx >= (getNBins() + 1)) {
-        return result;
+        return NeighborHoodIndices();
       }
       constexpr int min   = 1;
       const int     max   = getNBins();
       const int     itmin = std::max(min, int(idx - sizes.first));
       const int     itmax = std::min(max, int(idx + sizes.second));
-      for (int i = itmin; i <= itmax; i++) {
-        result.insert(i);
-      }
-      return result;
+      return NeighborHoodIndices(itmin, itmax + 1);
     }
 
     /// @brief Get #size bins which neighbor the one given
@@ -487,25 +594,44 @@ namespace detail {
     ///
     /// @param [in] idx requested bin index
     /// @param [in] sizes how many neighboring bins (up/down)
-    /// @return std::set of neighboring bin indices (global)
+    /// @return Set of neighboring bin indices (global)
     /// @note Closed varies given bin and allows bins on the opposite
     ///       side of the axis as neighbors. (excludes underflow / overflow)
     template <AxisBoundaryType T = bdt,
               std::enable_if_t<T == AxisBoundaryType::Closed, int> = 0>
-    std::set<size_t>
+    NeighborHoodIndices
     neighborHoodIndices(size_t idx,
                         std::pair<size_t, size_t> sizes = {1, 1}) const
     {
-      std::set<size_t> result;
+      // Handle invalid indices
       if (idx <= 0 || idx >= (getNBins() + 1)) {
-        return result;
+        return NeighborHoodIndices();
       }
-      const int itmin = idx - sizes.first;
-      const int itmax = idx + sizes.second;
-      for (int i = itmin; i <= itmax; i++) {
-        result.insert(wrapBin(i));
+
+      // Handle corner case where user requests more neighbours than the number
+      // of bins on the axis. We do not want to double-count bins in that case.
+      sizes.first %= getNBins();
+      sizes.second %= getNBins();
+      if (sizes.first + sizes.second + 1 > getNBins()) {
+        sizes.second -= (sizes.first + sizes.second + 1) - getNBins();
+      }
+
+      // If the entire index range is not covered, we must wrap the range of
+      // targeted neighbor indices into the range of valid bin indices. This may
+      // split the range of neighbor indices in two parts:
+      //
+      // Before wraparound - [        XXXXX]XXX
+      // After wraparound  - [ XXXX   XXXX ]
+      //
+      const int    itmin   = idx - sizes.first;
+      const int    itmax   = idx + sizes.second;
+      const size_t itfirst = wrapBin(itmin);
+      const size_t itlast  = wrapBin(itmax);
+      if (itfirst <= itlast) {
+        return NeighborHoodIndices(itfirst, itlast + 1);
+      } else {
+        return NeighborHoodIndices(itfirst, getNBins() + 1, 1, itlast + 1);
       }
-      return result;
     }
 
     /// @brief Converts bin index into a valid one for this axis.

Core/include/Acts/Utilities/detail/Grid.hpp

@@ -74,9 +74,9 @@ namespace detail {
     //
     template <class Point>
     reference
-    at(const Point& point)
+    atPosition(const Point& point)
     {
-      return m_values.at(getGlobalBinIndex(point));
+      return m_values.at(globalBinFromPosition(point));
     }
 
     /// @brief access value stored in bin for a given point
@@ -95,9 +95,9 @@ namespace detail {
     ///       Therefore, the look-up will never fail.
     template <class Point>
     const_reference
-    at(const Point& point) const
+    atPosition(const Point& point) const
     {
-      return m_values.at(getGlobalBinIndex(point));
+      return m_values.at(globalBinFromPosition(point));
     }
 
     /// @brief access value stored in bin with given global bin number
@@ -131,9 +131,9 @@ namespace detail {
     /// @pre All local bin indices must be a valid index for the corresponding
     ///      axis (including the under-/overflow bin for this axis).
     reference
-    at(const index_t& localBins)
+    atLocalBins(const index_t& localBins)
     {
-      return m_values.at(getGlobalBinIndex(localBins));
+      return m_values.at(globalBinFromLocalBins(localBins));
     }
 
     /// @brief access value stored in bin with given local bin numbers
@@ -145,9 +145,9 @@ namespace detail {
     /// @pre All local bin indices must be a valid index for the corresponding
     ///      axis (including the under-/overflow bin for this axis).
     const_reference
-    at(const index_t& localBins) const
+    atLocalBins(const index_t& localBins) const
     {
-      return m_values.at(getGlobalBinIndex(localBins));
+      return m_values.at(globalBinFromLocalBins(localBins));
     }
 
     /// @brief get global bin indices for closest points on grid
@@ -155,18 +155,17 @@ namespace detail {
     /// @tparam Point any type with point semantics supporting component access
     ///               through @c operator[]
     /// @param [in] position point of interest
-    /// @return set of global bin indices for bins whose lower-left corners are
-    ///         the closest points on the grid to the given point
+    /// @return Iterable thatemits the indices of bins whose lower-left corners
+    ///         are the closest points on the grid to the input.
     ///
     /// @pre The given @c Point type must represent a point in d (or higher)
     ///      dimensions where d is dimensionality of the grid. It must lie
     ///      within the grid range (i.e. not within a under-/overflow bin).
     template <class Point>
-    std::set<size_t>
+    detail::GlobalNeighborHoodIndices<DIM>
     closestPointsIndices(const Point& position) const
     {
-      return grid_helper::closestPointsIndices(getGlobalBinIndex(position),
-                                               m_axes);
+      return rawClosestPointsIndices(localBinsFromPosition(position));
     }
 
     /// @brief dimensionality of grid
@@ -186,7 +185,7 @@ namespace detail {
     /// @pre All local bin indices must be a valid index for the corresponding
     ///      axis (excluding the under-/overflow bins for each axis).
     std::array<double, DIM>
-    getBinCenter(const index_t& localBins) const
+    binCenter(const index_t& localBins) const
     {
       return grid_helper::getBinCenter(localBins, m_axes);
     }
@@ -204,9 +203,9 @@ namespace detail {
     /// @note This could be a under-/overflow bin along one or more axes.
     template <class Point>
     size_t
-    getGlobalBinIndex(const Point& point) const
+    globalBinFromPosition(const Point& point) const
     {
-      return grid_helper::getGlobalBin(point, m_axes);
+      return globalBinFromLocalBins(localBinsFromPosition(point));
     }
 
     /// @brief determine global bin index from local bin indices along each axis
@@ -217,11 +216,30 @@ namespace detail {
     /// @pre All local bin indices must be a valid index for the corresponding
     ///      axis (including the under-/overflow bin for this axis).
     size_t
-    getGlobalBinIndex(const index_t& localBins) const
+    globalBinFromLocalBins(const index_t& localBins) const
     {
       return grid_helper::getGlobalBin(localBins, m_axes);
     }
 
+    /// @brief  determine local bin index for each axis from the given point
+    ///
+    /// @tparam Point any type with point semantics supporting component access
+    ///               through @c operator[]
+    ///
+    /// @param  [in] point point to look up in the grid
+    /// @return array with local bin indices along each axis (in same order as
+    ///         given @c axes object)
+    ///
+    /// @pre The given @c Point type must represent a point in d (or higher)
+    ///      dimensions where d is dimensionality of the grid.
+    /// @note This could be a under-/overflow bin along one or more axes.
+    template <class Point>
+    index_t
+    localBinsFromPosition(const Point& point) const
+    {
+      return grid_helper::getLocalBinIndices(point, m_axes);
+    }
+
     /// @brief determine local bin index for each axis from global bin index
     ///
     /// @param  [in] bin global bin index
@@ -231,7 +249,7 @@ namespace detail {
     /// @note Local bin indices can contain under-/overflow bins along the
     ///       corresponding axis.
     index_t
-    getLocalBinIndices(size_t bin) const
+    localBinsFromGlobalBin(size_t bin) const
     {
       return grid_helper::getLocalBinIndices(bin, m_axes);
     }
@@ -244,18 +262,31 @@ namespace detail {
     /// @pre @c localBins must only contain valid bin indices (excluding
     ///      underflow bins).
     point_t
-    getLowerLeftBinEdge(const index_t& localBins) const
+    lowerLeftBinEdge(const index_t& localBins) const
     {
       return grid_helper::getLowerLeftBinEdge(localBins, m_axes);
     }
 
+    /// @brief retrieve upper-right bin edge from set of local bin indices
+    ///
+    /// @param  [in] localBins local bin indices along each axis
+    /// @return generalized upper-right bin edge position
+    ///
+    /// @pre @c localBins must only contain valid bin indices (excluding
+    ///      overflow bins).
+    point_t
+    upperRightBinEdge(const index_t& localBins) const
+    {
+      return grid_helper::getUpperRightBinEdge(localBins, m_axes);
+    }
+
     /// @brief get number of bins along each specific axis
     ///
     /// @return array giving the number of bins along all axes
     ///
     /// @note Not including under- and overflow bins
     index_t
-    getNBins() const
+    numLocalBins() const
     {
       return grid_helper::getNBins(m_axes);
     }
@@ -264,7 +295,7 @@ namespace detail {
     ///
     /// @return array returning the minima of all given axes
     point_t
-    getMin() const
+    minPosition() const
     {
       return grid_helper::getMin(m_axes);
     }
@@ -273,24 +304,11 @@ namespace detail {
     ///
     /// @return array returning the maxima of all given axes
     point_t
-    getMax() const
+    maxPosition() const
     {
       return grid_helper::getMax(m_axes);
     }
 
-    /// @brief retrieve upper-right bin edge from set of local bin indices
-    ///
-    /// @param  [in] localBins local bin indices along each axis
-    /// @return generalized upper-right bin edge position
-    ///
-    /// @pre @c localBins must only contain valid bin indices (excluding
-    ///      overflow bins).
-    point_t
-    getUpperRightBinEdge(const index_t& localBins) const
-    {
-      return grid_helper::getUpperRightBinEdge(localBins, m_axes);
-    }
-
     /// @brief set all overflow and underflow bins to a certain value
     ///
     /// @param [in] value value to be inserted in every overflow and underflow
@@ -348,25 +366,20 @@ namespace detail {
       // get local indices for current bin
       // value of bin is interpreted as being the field value at its lower left
       // corner
-      const auto& llIndices = getLocalBinIndices(getGlobalBinIndex(point));
-
-      // get local indices for "upper right" bin (i.e. all local indices
-      // incremented by one but avoid overflows)
-      auto urIndices = grid_helper::getUpperRightBinIndices(llIndices, m_axes);
+      const auto& llIndices = localBinsFromPosition(point);
 
       // get global indices for all surrounding corner points
-      const auto& closestIndices = closestPointsIndices(point);
+      const auto& closestIndices = rawClosestPointsIndices(llIndices);
 
       // get values on grid points
       size_t i = 0;
       for (size_t index : closestIndices) {
-        // In case it is the last bin return the last value as neighbour
-        neighbors.at(i++) = (index < size()) ? at(index) : at(size() - 1);
+        neighbors.at(i++) = at(index);
       }
 
       return Acts::interpolate(point,
-                               getLowerLeftBinEdge(llIndices),
-                               getLowerLeftBinEdge(urIndices),
+                               lowerLeftBinEdge(llIndices),
+                               upperRightBinEdge(llIndices),
                                neighbors);
     }
 
@@ -403,26 +416,12 @@ namespace detail {
     ///       Ignoring the truncation of the neighborhood size reaching beyond
     ///       over-/underflow bins, the neighborhood is of size \f$2 \times
     ///       \text{size}+1\f$ along each dimension.
-    std::set<size_t>
+    detail::GlobalNeighborHoodIndices<DIM>
     neighborHoodIndices(const index_t& localBins, size_t size = 1u) const
     {
       return grid_helper::neighborHoodIndices(localBins, size, m_axes);
     }
 
-    /// @brief get global bin indices for neighborhood of bin identified by @p
-    /// pos
-    /// @param pos position around which to look
-    /// @param size how many neighbors (defaul: 1)
-    /// @return set of global bin indices pointing to the neighbors
-    template <class Point>
-    std::set<size_t>
-    neighborHoodIndices(const Point& pos, size_t size = 1u) const
-    {
-      const size_t  bin       = getGlobalBinIndex(pos);
-      const index_t localBins = getLocalBinIndices(bin);
-      return grid_helper::neighborHoodIndices(localBins, size, m_axes);
-    }
-
     /// @brief total number of bins
     ///
     /// @return total number of bins in the grid
@@ -431,7 +430,7 @@ namespace detail {
     size_t
     size() const
     {
-      index_t nBinsArray = getNBins();
+      index_t nBinsArray = numLocalBins();
       // add under-and overflow bins for each axis and multiply all bins
       return std::accumulate(
           nBinsArray.begin(),
@@ -441,7 +440,7 @@ namespace detail {
     }
 
     std::array<const IAxis*, DIM>
-    getAxes() const
+    axes() const
     {
       return grid_helper::getAxes(m_axes);
     }
@@ -451,6 +450,15 @@ namespace detail {
     std::tuple<Axes...> m_axes;
     /// linear value store for each bin
     std::vector<T> m_values;
+
+    // Part of closestPointsIndices that goes after local bins resolution.
+    // Used as an interpolation performance optimization, but not exposed as it
+    // doesn't make that much sense from an API design standpoint.
+    detail::GlobalNeighborHoodIndices<DIM>
+    rawClosestPointsIndices(const index_t& localBins) const
+    {
+      return grid_helper::closestPointsIndices(localBins, m_axes);
+    }
   };
 }  // namespace detail
 

Core/include/Acts/Utilities/detail/grid_helper.hpp

@@ -17,6 +17,144 @@ namespace Acts {
 
 namespace detail {
 
+  // This object can be iterated to produce the (ordered) set of global indices
+  // associated with a neighborhood around a certain point on a grid.
+  //
+  // The goal is to emulate the effect of enumerating the global indices into
+  // an std::set (or into an std::vector that gets subsequently sorted), without
+  // paying the price of dynamic memory allocation in hot magnetic field
+  // interpolation code.
+  //
+  template <size_t DIM>
+  class GlobalNeighborHoodIndices
+  {
+  public:
+    // You can get the local neighbor indices from
+    // grid_helper_impl<DIM>::neighborHoodIndices and the number of bins in
+    // each direction from grid_helper_impl<DIM>::getNBins.
+    GlobalNeighborHoodIndices(
+        std::array<NeighborHoodIndices, DIM>& neighborIndices,
+        const std::array<size_t, DIM>&        nBinsArray)
+      : m_localIndices(neighborIndices)
+    {
+      if (DIM == 1) return;
+      size_t globalStride = 1;
+      for (long i = DIM - 2; i >= 0; --i) {
+        globalStride *= (nBinsArray[i + 1] + 2);
+        m_globalStrides[i] = globalStride;
+      }
+    }
+
+    class iterator
+    {
+    public:
+      iterator() = default;
+
+      iterator(
+          const GlobalNeighborHoodIndices& parent,
+          std::array<NeighborHoodIndices::iterator, DIM>&& localIndicesIter)
+        : m_localIndicesIter(std::move(localIndicesIter)), m_parent(&parent)
+      {
+      }
+
+      size_t operator*() const
+      {
+        size_t globalIndex = *m_localIndicesIter[DIM - 1];
+        if (DIM == 1) return globalIndex;
+        for (size_t i = 0; i < DIM - 1; ++i) {
+          globalIndex
+              += m_parent->m_globalStrides[i] * (*m_localIndicesIter[i]);
+        }
+        return globalIndex;
+      }
+
+      iterator& operator++()
+      {
+        const auto& localIndices = m_parent->m_localIndices;
+
+        // Go to the next global index via a lexicographic increment:
+        // - Start by incrementing the last local index
+        // - If it reaches the end, reset it and increment the previous one...
+        for (long i = DIM - 1; i > 0; --i) {
+          ++m_localIndicesIter[i];
+          if (m_localIndicesIter[i] != localIndices[i].end()) return *this;
+          m_localIndicesIter[i] = localIndices[i].begin();
+        }
+
+        // The first index should stay at the end value when it reaches it, so
+        // that we know when we've reached the end of iteration.
+        ++m_localIndicesIter[0];
+        return *this;
+      }
+
+      bool
+      operator==(const iterator& it)
+      {
+        // We know when we've reached the end, so we don't need an end-iterator.
+        // Sadly, in C++, there has to be one. Therefore, we special-case it
+        // heavily so that it's super-efficient to create and compare to.
+        if (it.m_parent == nullptr) {
+          return m_localIndicesIter[0] == m_parent->m_localIndices[0].end();
+        } else {
+          return m_localIndicesIter == it.m_localIndicesIter;
+        }
+      }
+
+      bool
+      operator!=(const iterator& it)
+      {
+        return !(*this == it);
+      }
+
+    private:
+      std::array<NeighborHoodIndices::iterator, DIM> m_localIndicesIter;
+      const GlobalNeighborHoodIndices* m_parent = nullptr;
+    };
+
+    iterator
+    begin() const
+    {
+      std::array<NeighborHoodIndices::iterator, DIM> localIndicesIter;
+      for (size_t i = 0; i < DIM; ++i) {
+        localIndicesIter[i] = m_localIndices[i].begin();
+      }
+      return iterator(*this, std::move(localIndicesIter));
+    }
+
+    iterator
+    end() const
+    {
+      return iterator();
+    }
+
+    // Number of indices that will be produced if this sequence is iterated
+    size_t
+    size() const
+    {
+      size_t result = m_localIndices[0].size();
+      for (size_t i = 1; i < DIM; ++i) {
+        result *= m_localIndices[i].size();
+      }
+      return result;
+    }
+
+    // Collect the sequence of indices into an std::vector
+    std::vector<size_t>
+    collect() const
+    {
+      std::vector<size_t> result;
+      result.reserve(this->size());
+      for (size_t idx : *this) {
+        result.push_back(idx);
+      }
+      return result;
+    }
+
+  private:
+    std::array<NeighborHoodIndices, DIM> m_localIndices;
+    std::array<size_t, DIM - 1>          m_globalStrides;
+  };
+
   /// @cond
   /// @brief helper struct to calculate number of bins inside a grid
   ///
@@ -37,32 +175,29 @@ namespace detail {
       grid_helper_impl<N - 1>::getBinCenter(center, localIndices, axes);
     }
 
-    template <class Point, class... Axes>
+    template <class... Axes>
     static void
-    getGlobalBin(const Point&               point,
+    getGlobalBin(const std::array<size_t, sizeof...(Axes)>& localBins,
                  const std::tuple<Axes...>& axes,
                  size_t&                    bin,
                  size_t&                    area)
     {
       const auto& thisAxis = std::get<N>(axes);
-      bin += area * thisAxis.getBin(point[N]);
+      bin += area * localBins.at(N);
       // make sure to account for under-/overflow bins
       area *= (thisAxis.getNBins() + 2);
-      grid_helper_impl<N - 1>::getGlobalBin(point, axes, bin, area);
+      grid_helper_impl<N - 1>::getGlobalBin(localBins, axes, bin, area);
     }
 
-    template <class... Axes>
+    template <class Point, class... Axes>
     static void
-    getGlobalBin(const std::array<size_t, sizeof...(Axes)>& localBins,
-                 const std::tuple<Axes...>& axes,
-                 size_t&                    bin,
-                 size_t&                    area)
+    getLocalBinIndices(const Point&               point,
+                       const std::tuple<Axes...>& axes,
+                       std::array<size_t, sizeof...(Axes)>& indices)
     {
       const auto& thisAxis = std::get<N>(axes);
-      bin += area * localBins.at(N);
-      // make sure to account for under-/overflow bins
-      area *= (thisAxis.getNBins() + 2);
-      grid_helper_impl<N - 1>::getGlobalBin(localBins, axes, bin, area);
+      indices.at(N)        = thisAxis.getBin(point[N]);
+      grid_helper_impl<N - 1>::getLocalBinIndices(point, axes, indices);
     }
 
     template <class... Axes>
@@ -171,12 +306,12 @@ namespace detail {
         const std::array<size_t, sizeof...(Axes)>& localIndices,
         std::pair<size_t, size_t>                  sizes,
         const std::tuple<Axes...>& axes,
-        std::array<std::set<size_t>, sizeof...(Axes)>& neighborIndices)
+        std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices)
     {
 
       // ask n-th axis
-      size_t           locIdx = localIndices.at(N);
-      std::set<size_t> locNeighbors
+      size_t              locIdx = localIndices.at(N);
+      NeighborHoodIndices locNeighbors
           = std::get<N>(axes).neighborHoodIndices(locIdx, sizes);
       neighborIndices.at(N) = locNeighbors;
 
@@ -184,24 +319,6 @@ namespace detail {
           localIndices, sizes, axes, neighborIndices);
     }
 
-    template <class... Axes>
-    static void
-    combineNeighborHoodIndices(
-        std::array<size_t, sizeof...(Axes)>&           idx,
-        std::array<std::set<size_t>, sizeof...(Axes)>& neighborIndices,
-        std::set<size_t>&          combinations,
-        const std::tuple<Axes...>& axes)
-    {
-      // iterate over this axis' neighbors
-      std::set<size_t>& neighbors = neighborIndices.at(N);
-      for (const auto& i : neighbors) {
-        idx.at(N) = i;
-        // vary other axes recursively
-        grid_helper_impl<N - 1>::combineNeighborHoodIndices(
-            idx, neighborIndices, combinations, axes);
-      }
-    }
-
     template <class... Axes>
     static void
     exteriorBinIndices(std::array<size_t, sizeof...(Axes)>& idx,
@@ -233,17 +350,6 @@ namespace detail {
       center.at(0u) = std::get<0u>(axes).getBinCenter(localIndices.at(0u));
     }
 
-    template <class Point, class... Axes>
-    static void
-    getGlobalBin(const Point&               point,
-                 const std::tuple<Axes...>& axes,
-                 size_t&                    bin,
-                 size_t&                    area)
-    {
-      const auto& thisAxis = std::get<0u>(axes);
-      bin += area * thisAxis.getBin(point[0u]);
-    }
-
     template <class... Axes>
     static void
     getGlobalBin(const std::array<size_t, sizeof...(Axes)>& localBins,
@@ -254,6 +360,16 @@ namespace detail {
       bin += area * localBins.at(0u);
     }
 
+    template <class Point, class... Axes>
+    static void
+    getLocalBinIndices(const Point&               point,
+                       const std::tuple<Axes...>& axes,
+                       std::array<size_t, sizeof...(Axes)>& indices)
+    {
+      const auto& thisAxis = std::get<0u>(axes);
+      indices.at(0u)       = thisAxis.getBin(point[0u]);
+    }
+
     template <class... Axes>
     static void
     getLocalBinIndices(size_t& bin,
@@ -347,36 +463,16 @@ namespace detail {
         const std::array<size_t, sizeof...(Axes)>& localIndices,
         std::pair<size_t, size_t>                  sizes,
         const std::tuple<Axes...>& axes,
-        std::array<std::set<size_t>, sizeof...(Axes)>& neighborIndices)
+        std::array<NeighborHoodIndices, sizeof...(Axes)>& neighborIndices)
     {
 
       // ask 0-th axis
-      size_t           locIdx = localIndices.at(0u);
-      std::set<size_t> locNeighbors
+      size_t              locIdx = localIndices.at(0u);
+      NeighborHoodIndices locNeighbors
           = std::get<0u>(axes).neighborHoodIndices(locIdx, sizes);
       neighborIndices.at(0u) = locNeighbors;
     }
 
-    template <class... Axes>
-    static void
-    combineNeighborHoodIndices(
-        std::array<size_t, sizeof...(Axes)>&           idx,
-        std::array<std::set<size_t>, sizeof...(Axes)>& neighborIndices,
-        std::set<size_t>&          combinations,
-        const std::tuple<Axes...>& axes)
-    {
-      // iterate over this axis' neighbors
-      std::set<size_t>& neighbors = neighborIndices.at(0u);
-      for (const auto& i : neighbors) {
-        idx.at(0u) = i;
-        // at this point, combinations are complete, so fill
-        size_t bin = 0, area = 1;
-        grid_helper_impl<sizeof...(Axes) - 1>::getGlobalBin(
-            idx, axes, bin, area);
-        combinations.insert(bin);
-      }
-    }
-
     template <class... Axes>
     static void
     exteriorBinIndices(std::array<size_t, sizeof...(Axes)>& idx,
@@ -425,22 +521,20 @@ namespace detail {
     /// @tparam Axes parameter pack of axis types defining the grid
     /// @param  [in] bin  global bin index for bin of interest
     /// @param  [in] axes actual axis objects spanning the grid
-    /// @return set of global bin indices for bins whose lower-left corners are
-    ///         the closest points on the grid to every point in the given bin
+    /// @return Sorted collection of global bin indices for bins whose
+    ///         lower-left corners are the closest points on the grid to every
+    ///         point in the given bin
     ///
     /// @note @c bin must be a valid bin index (excluding under-/overflow bins
     ///       along any axis).
     template <class... Axes>
-    static std::set<size_t>
-    closestPointsIndices(size_t bin, const std::tuple<Axes...>& axes)
+    static GlobalNeighborHoodIndices<sizeof...(Axes)>
+    closestPointsIndices(
+        const std::array<size_t, sizeof...(Axes)>& localIndices,
+        const std::tuple<Axes...>& axes)
     {
-      std::array<size_t, sizeof...(Axes)> localIndices
-          = getLocalBinIndices(bin, axes);
       // get neighboring bins, but only increment.
-      std::array<std::set<size_t>, sizeof...(Axes)> neighborIndices;
-      std::set<size_t> comb
-          = neighborHoodIndices(localIndices, std::make_pair(0, 1), axes);
-      return comb;
+      return neighborHoodIndices(localIndices, std::make_pair(0, 1), axes);
     }
 
     /// @brief retrieve bin center from set of local bin indices
@@ -464,54 +558,54 @@ namespace detail {
       return center;
     }
 
-    /// @brief determine global bin index in grid defined by a set of axes
+    /// @brief determine global bin index from local indices along each axis
     ///
-    /// @tparam Point any type with point semantics supporting component access
-    ///               through @c operator[]
     /// @tparam Axes parameter pack of axis types defining the grid
     ///
-    /// @param  [in] point point to look up in the grid
+    /// @param  [in] localBins local bin indices along each axis
     /// @param  [in] axes  actual axis objects spanning the grid
-    /// @return global index for bin containing the given point
+    /// @return global index for bin defined by the local bin indices
     ///
-    /// @pre The given @c Point type must represent a point in d (or higher)
-    ///      dimensions where d is the number of axis objects in the tuple.
-    /// @note This could be a under-/overflow bin along one or more axes.
-    template <class Point, class... Axes>
+    /// @pre All local bin indices must be a valid index for the corresponding
+    ///      axis (including the under-/overflow bin for this axis).
+    template <class... Axes>
     static size_t
-    getGlobalBin(const Point& point, const std::tuple<Axes...>& axes)
+    getGlobalBin(const std::array<size_t, sizeof...(Axes)>& localBins,
+                 const std::tuple<Axes...>& axes)
     {
       constexpr size_t MAX  = sizeof...(Axes) - 1;
       size_t           area = 1;
       size_t           bin  = 0;
 
-      grid_helper_impl<MAX>::getGlobalBin(point, axes, bin, area);
+      grid_helper_impl<MAX>::getGlobalBin(localBins, axes, bin, area);
 
       return bin;
     }
 
-    /// @brief determine global bin index from local indices along each axis
+    /// @brief determine local bin index for each axis from point
     ///
+    /// @tparam Point any type with point semantics supporting component access
+    ///               through @c operator[]
     /// @tparam Axes parameter pack of axis types defining the grid
     ///
-    /// @param  [in] localBins local bin indices along each axis
+    /// @param  [in] point point to look up in the grid
     /// @param  [in] axes  actual axis objects spanning the grid
-    /// @return global index for bin defined by the local bin indices
+    /// @return array with local bin indices along each axis (in same order as
+    ///         given @c axes object)
     ///
-    /// @pre All local bin indices must be a valid index for the corresponding
-    ///      axis (including the under-/overflow bin for this axis).
-    template <class... Axes>
-    static size_t
-    getGlobalBin(const std::array<size_t, sizeof...(Axes)>& localBins,
-                 const std::tuple<Axes...>& axes)
+    /// @pre The given @c Point type must represent a point in d (or higher)
+    ///      dimensions where d is the number of axis objects in the tuple.
+    /// @note This could be a under-/overflow bin along one or more axes.
+    template <class Point, class... Axes>
+    static std::array<size_t, sizeof...(Axes)>
+    getLocalBinIndices(const Point& point, const std::tuple<Axes...>& axes)
     {
-      constexpr size_t MAX  = sizeof...(Axes) - 1;
-      size_t           area = 1;
-      size_t           bin  = 0;
+      constexpr size_t MAX = sizeof...(Axes) - 1;
+      std::array<size_t, sizeof...(Axes)> indices;
 
-      grid_helper_impl<MAX>::getGlobalBin(localBins, axes, bin, area);
+      grid_helper_impl<MAX>::getLocalBinIndices(point, axes, indices);
 
-      return bin;
+      return indices;
     }
 
     /// @brief determine local bin index for each axis from global bin index
@@ -699,7 +793,8 @@ namespace detail {
     /// @param  [in] size         size of neighborhood determining how many
     ///                           adjacent bins along each axis are considered
     /// @param  [in] axes         actual axis objects spanning the grid
-    /// @return set of global bin indices for all bins in neighborhood
+    /// @return Sorted collection of global bin indices for all bins in
+    ///         the neighborhood
     ///
     /// @note Over-/underflow bins are included in the neighborhood.
     /// @note The @c size parameter sets the range by how many units each local
@@ -711,13 +806,8 @@ namespace detail {
     /// @note The concrete bins which are returned depend on the WrappingTypes
     ///       of the contained axes
     ///
-    /// @todo At the moment, this function relies on local-to-global bin index
-    ///       conversions for each neighboring bin. In principle, all
-    ///       information should be available to work directly on global bin
-    ///       indices. The problematic part is the check when going beyond
-    ///       under-/overflow bins.
     template <class... Axes>
-    static std::set<size_t>
+    static GlobalNeighborHoodIndices<sizeof...(Axes)>
     neighborHoodIndices(const std::array<size_t, sizeof...(Axes)>& localIndices,
                         std::pair<size_t, size_t>                  sizes,
                         const std::tuple<Axes...>& axes)
@@ -725,21 +815,20 @@ namespace detail {
       constexpr size_t MAX = sizeof...(Axes) - 1;
 
       // length N array which contains local neighbors based on size par
-      std::array<std::set<size_t>, sizeof...(Axes)> neighborIndices;
+      std::array<NeighborHoodIndices, sizeof...(Axes)> neighborIndices;
       // get local bin indices for neighboring bins
       grid_helper_impl<MAX>::neighborHoodIndices(
           localIndices, sizes, axes, neighborIndices);
 
-      std::array<size_t, sizeof...(Axes)> idx;
-      std::set<size_t> combinations;
-      grid_helper_impl<MAX>::combineNeighborHoodIndices(
-          idx, neighborIndices, combinations, axes);
+      // Query the number of bins
+      std::array<size_t, sizeof...(Axes)> nBinsArray = getNBins(axes);
 
-      return combinations;
+      // Produce iterator of global indices
+      return GlobalNeighborHoodIndices(neighborIndices, nBinsArray);
     }
 
     template <class... Axes>
-    static std::set<size_t>
+    static GlobalNeighborHoodIndices<sizeof...(Axes)>
     neighborHoodIndices(const std::array<size_t, sizeof...(Axes)>& localIndices,
                         size_t                     size,
                         const std::tuple<Axes...>& axes)

Core/src/Detector/TrackingGeometry.cpp

 // This file is part of the Acts project.
 //
-// Copyright (C) 2016-2018 Acts project team
+// Copyright (C) 2016-2019 Acts project team
 //
 // This Source Code Form is subject to the terms of the Mozilla Public
 // License, v. 2.0. If a copy of the MPL was not distributed with this
@@ -19,13 +19,14 @@
 #include "Acts/Surfaces/Surface.hpp"
 
 Acts::TrackingGeometry::TrackingGeometry(
-    const MutableTrackingVolumePtr& highestVolume)
+    const MutableTrackingVolumePtr& highestVolume,
+    const SurfaceMaterialMap&       surfaceMaterialMap)
   : m_world(highestVolume)
   , m_beam(Surface::makeShared<PerigeeSurface>(s_origin))
 {
   // close up the geometry
   size_t volumeID = 0;
-  highestVolume->closeGeometry(m_trackingVolumes, volumeID);
+  highestVolume->closeGeometry(surfaceMaterialMap, m_trackingVolumes, volumeID);
 }
 
 Acts::TrackingGeometry::~TrackingGeometry() = default;

Core/src/Detector/TrackingVolume.cpp

@@ -308,6 +308,7 @@ Acts::TrackingVolume::interlinkLayers()
 
 void
 Acts::TrackingVolume::closeGeometry(
+    const SurfaceMaterialMap& surfaceMaterialMap,
     std::map<std::string, const TrackingVolume*>& volumeMap,
     size_t& vol)
 {
@@ -321,6 +322,16 @@ Acts::TrackingVolume::closeGeometry(
   auto thisVolume = const_cast<TrackingVolume*>(this);
   thisVolume->assignGeoID(volumeID);
 
+  // This functor checks and assigns the material for a given
+  auto assignSurfaceMaterial
+      = [&surfaceMaterialMap](Surface& sf, GeometryID geoID) -> void {
+    // Try to find the surface in the map
+    auto sMaterial = surfaceMaterialMap.find(geoID);
+    if (sMaterial != surfaceMaterialMap.end()) {
+      sf.assignSurfaceMaterial(sMaterial->second);
+    }
+  };
+
   // loop over the boundary surfaces
   geo_id_value iboundary = 0;
   // loop over the boundary surfaces
@@ -333,6 +344,7 @@ Acts::TrackingVolume::closeGeometry(
     // now assign to the boundary surface
     auto& mutableBSurface = *(const_cast<Surface*>(&bSurface));
     mutableBSurface.assignGeoID(boundaryID);
+    assignSurfaceMaterial(mutableBSurface, boundaryID);
   }
 
   // A) this is NOT a container volume, volumeID is already incremented
@@ -347,7 +359,7 @@ Acts::TrackingVolume::closeGeometry(
         layerID.add(++ilayer, GeometryID::layer_mask);
         // now close the geometry
         auto mutableLayerPtr = std::const_pointer_cast<Layer>(layerPtr);
-        mutableLayerPtr->closeGeometry(layerID);
+        mutableLayerPtr->closeGeometry(surfaceMaterialMap, layerID);
       }
     }
   } else {
@@ -356,7 +368,7 @@ Acts::TrackingVolume::closeGeometry(
     for (auto& volumesIter : m_confinedVolumes->arrayObjects()) {
       auto mutableVolumesIter
           = std::const_pointer_cast<TrackingVolume>(volumesIter);
-      mutableVolumesIter->closeGeometry(volumeMap, vol);
+      mutableVolumesIter->closeGeometry(surfaceMaterialMap, volumeMap, vol);
     }
   }
 }

Core/src/Layers/Layer.cpp

@@ -53,10 +53,25 @@ Acts::Layer::approachDescriptor()
 }
 
 void
-Acts::Layer::closeGeometry(const GeometryID& layerID)
+Acts::Layer::closeGeometry(const SurfaceMaterialMap& surfaceMaterialMap,
+                           const GeometryID&         layerID)
 {
+
+  // This functor checks and assigns the material for a given
+  auto assignSurfaceMaterial
+      = [&surfaceMaterialMap](Surface* sf, GeometryID geoID) -> void {
+    // Try to find the surface in the map
+    auto sMaterial = surfaceMaterialMap.find(geoID);
+    if (sMaterial != surfaceMaterialMap.end()) {
+      sf->assignSurfaceMaterial(sMaterial->second);
+    }
+  };
+
   // set the volumeID of this
   assignGeoID(layerID);
+  // assign to the representing surface
+  Surface* rSurface = const_cast<Surface*>(&surfaceRepresentation());
+  assignSurfaceMaterial(rSurface, layerID);
 
   // also find out how the sub structure is defined
   if (surfaceRepresentation().associatedMaterial() != nullptr) {
@@ -74,6 +89,7 @@ Acts::Layer::closeGeometry(const GeometryID& layerID)
       asurfaceID.add(++iasurface, GeometryID::approach_mask);
       auto mutableASurface = const_cast<Surface*>(aSurface);
       mutableASurface->assignGeoID(asurfaceID);
+      assignSurfaceMaterial(mutableASurface, asurfaceID);
       // if any of the approach surfaces has material
       if (aSurface->associatedMaterial() != nullptr) {
         m_ssApproachSurfaces = 2;
@@ -91,6 +107,7 @@ Acts::Layer::closeGeometry(const GeometryID& layerID)
       ssurfaceID.add(++issurface, GeometryID::sensitive_mask);
       auto mutableSSurface = const_cast<Surface*>(sSurface);
       mutableSSurface->assignGeoID(ssurfaceID);
+      assignSurfaceMaterial(mutableSSurface, ssurfaceID);
       // if any of the sensitive surfaces has material
       if (sSurface->associatedMaterial() != nullptr) {
         m_ssSensitiveSurfaces = 2;

Core/src/Tools/CuboidVolumeBuilder.cpp

@@ -47,7 +47,7 @@ Acts::CuboidVolumeBuilder::buildSurface(
     surface = Surface::makeShared<PlaneSurface>(
         std::make_shared<const Transform3D>(trafo), cfg.rBounds);
   }
-  surface->setAssociatedMaterial(cfg.surMat);
+  surface->assignSurfaceMaterial(cfg.surMat);
   return surface;
 }
 

Core/src/Tools/PassiveLayerBuilder.cpp

@@ -96,7 +96,7 @@ Acts::PassiveLayerBuilder::endcapLayers(const Acts::GeometryContext& /*gctx*/,
         // create homogeneous material
         material = m_cfg.posnegLayerMaterial.at(ipnl);
         // sign it to the surface
-        eLayer->surfaceRepresentation().setAssociatedMaterial(material);
+        eLayer->surfaceRepresentation().assignSurfaceMaterial(material);
       }
       // push it into the layer vector
       eLayers.push_back(eLayer);
@@ -137,7 +137,7 @@ Acts::PassiveLayerBuilder::centralLayers(
         // create homogeneous material
         material = m_cfg.centralLayerMaterial.at(icl);
         // sign it to the surface
-        cLayer->surfaceRepresentation().setAssociatedMaterial(material);
+        cLayer->surfaceRepresentation().assignSurfaceMaterial(material);
       }
       // push it into the layer vector
       cLayers.push_back(cLayer);

Core/src/Utilities/BFieldMapUtils.cpp

@@ -86,14 +86,14 @@ Acts::
         size_t          n = std::abs(int(j) - int(zPos.size()));
         Grid_t::index_t indicesFirstQuadrant = {{i - 1, n}};
 
-        grid.at(indices)
+        grid.atLocalBins(indices)
             = bField.at(localToGlobalBin(indicesFirstQuadrant, nIndices))
             * BFieldUnit;
       } else {
         // std::vectors begin with 0 and we do not want the user needing to
         // take underflow or overflow bins in account this is why we need to
         // subtract by one
-        grid.at(indices)
+        grid.atLocalBins(indices)
             = bField.at(localToGlobalBin({{i - 1, j - 1}}, nIndices))
             * BFieldUnit;
       }
@@ -109,11 +109,20 @@ Acts::
 
   // [4] Create the transformation for the bfield
   // map (Br,Bz) -> (Bx,By,Bz)
-  auto transformBField
-      = [](const Acts::Vector2D& field, const Acts::Vector3D& pos) {
-          return Acts::Vector3D(
-              field.x() * cos(phi(pos)), field.x() * sin(phi(pos)), field.y());
-        };
+  auto transformBField = [](const Acts::Vector2D& field,
+                            const Acts::Vector3D& pos) {
+    double r_sin_theta_2 = pos.x() * pos.x() + pos.y() * pos.y();
+    double cos_phi, sin_phi;
+    if (r_sin_theta_2 > std::numeric_limits<double>::min()) {
+      double inv_r_sin_theta = 1. / sqrt(r_sin_theta_2);
+      cos_phi                = pos.x() * inv_r_sin_theta;
+      sin_phi                = pos.y() * inv_r_sin_theta;
+    } else {
+      cos_phi = 1.;
+      sin_phi = 0.;
+    }
+    return Acts::Vector3D(field.x() * cos_phi, field.x() * sin_phi, field.y());
+  };
 
   // [5] Create the mapper & BField Service
   // create field mapping
@@ -216,7 +225,7 @@ Acts::
           size_t          l = std::abs(int(k) - (int(zPos.size())));
           Grid_t::index_t indicesFirstOctant = {{m, n, l}};
 
-          grid.at(indices)
+          grid.atLocalBins(indices)
               = bField.at(localToGlobalBin(indicesFirstOctant, nIndices))
               * BFieldUnit;
 
@@ -224,7 +233,7 @@ Acts::
           // std::vectors begin with 0 and we do not want the user needing to
           // take underflow or overflow bins in account this is why we need to
           // subtract by one
-          grid.at(indices)
+          grid.atLocalBins(indices)
               = bField.at(localToGlobalBin({{i - 1, j - 1, k - 1}}, nIndices))
               * BFieldUnit;
         }
@@ -288,11 +297,21 @@ Acts::
 
   // Create the transformation for the bfield
   // map (Br,Bz) -> (Bx,By,Bz)
-  auto transformBField = [](const Acts::Vector2D& bfield,
-                            const Acts::Vector3D& pos) {
-    return Acts::Vector3D(
-        bfield.x() * cos(phi(pos)), bfield.x() * sin(phi(pos)), bfield.y());
-  };
+  auto transformBField
+      = [](const Acts::Vector2D& bfield, const Acts::Vector3D& pos) {
+          double r_sin_theta_2 = pos.x() * pos.x() + pos.y() * pos.y();
+          double cos_phi, sin_phi;
+          if (r_sin_theta_2 > std::numeric_limits<double>::min()) {
+            double inv_r_sin_theta = 1. / sqrt(r_sin_theta_2);
+            cos_phi                = pos.x() * inv_r_sin_theta;
+            sin_phi                = pos.y() * inv_r_sin_theta;
+          } else {
+            cos_phi = 1.;
+            sin_phi = 0.;
+          }
+          return Acts::Vector3D(
+              bfield.x() * cos_phi, bfield.x() * sin_phi, bfield.y());
+        };
 
   // iterate over all bins, set their value to the solenoid value
   // at their lower left position
@@ -301,13 +320,13 @@ Acts::
       Grid_t::index_t index({i, j});
       if (i == 0 || j == 0 || i == nBinsR + 1 || j == nBinsZ + 1) {
         // under or overflow bin, set zero
-        grid.at(index) = Grid_t::value_type(0, 0);
+        grid.atLocalBins(index) = Grid_t::value_type(0, 0);
       } else {
         // regular bin, get lower left boundary
-        Grid_t::point_t lowerLeft = grid.getLowerLeftBinEdge(index);
+        Grid_t::point_t lowerLeft = grid.lowerLeftBinEdge(index);
         // do lookup
-        Vector2D B     = field.getField(Vector2D(lowerLeft[0], lowerLeft[1]));
-        grid.at(index) = B;
+        Vector2D B = field.getField(Vector2D(lowerLeft[0], lowerLeft[1]));
+        grid.atLocalBins(index) = B;
       }
     }
   }

Plugins/DD4hep/src/DD4hepLayerBuilder.cpp

@@ -198,15 +198,15 @@ Acts::DD4hepLayerBuilder::negativeLayers(const GeometryContext& gctx) const
 
         // set material surface
         if (layerPos == Acts::LayerMaterialPos::inner) {
-          innerBoundary->setAssociatedMaterial(materialProxy);
+          innerBoundary->assignSurfaceMaterial(materialProxy);
         }
 
         if (layerPos == Acts::LayerMaterialPos::outer) {
-          outerBoundary->setAssociatedMaterial(materialProxy);
+          outerBoundary->assignSurfaceMaterial(materialProxy);
         }
 
         if (layerPos == Acts::LayerMaterialPos::central) {
-          centralSurface->setAssociatedMaterial(materialProxy);
+          centralSurface->assignSurfaceMaterial(materialProxy);
         }
 
         // collect approach surfaces
@@ -268,7 +268,7 @@ Acts::DD4hepLayerBuilder::negativeLayers(const GeometryContext& gctx) const
             materialProperties);
       }
 
-      negativeLayer->surfaceRepresentation().setAssociatedMaterial(
+      negativeLayer->surfaceRepresentation().assignSurfaceMaterial(
           surfMaterial);
       // push back created layer
       layers.push_back(negativeLayer);
@@ -403,15 +403,15 @@ Acts::DD4hepLayerBuilder::centralLayers(const GeometryContext& gctx) const
         // check if the material should be set to the inner or outer boundary
         // and set it in case
         if (layerPos == Acts::LayerMaterialPos::inner) {
-          innerBoundary->setAssociatedMaterial(materialProxy);
+          innerBoundary->assignSurfaceMaterial(materialProxy);
         }
 
         if (layerPos == Acts::LayerMaterialPos::outer) {
-          outerBoundary->setAssociatedMaterial(materialProxy);
+          outerBoundary->assignSurfaceMaterial(materialProxy);
         }
 
         if (layerPos == Acts::LayerMaterialPos::central) {
-          centralSurface->setAssociatedMaterial(materialProxy);
+          centralSurface->assignSurfaceMaterial(materialProxy);
         }
 
         // collect the surfaces
@@ -475,10 +475,10 @@ Acts::DD4hepLayerBuilder::centralLayers(const GeometryContext& gctx) const
         surfMaterial = std::make_shared<const HomogeneousSurfaceMaterial>(
             materialProperties);
 
-        //   innerBoundary->setAssociatedMaterial(surfMaterial);
+        //   innerBoundary->assignSurfaceMaterial(surfMaterial);
       }
 
-      centralLayer->surfaceRepresentation().setAssociatedMaterial(surfMaterial);
+      centralLayer->surfaceRepresentation().assignSurfaceMaterial(surfMaterial);
 
       // push back created layer
       layers.push_back(centralLayer);
@@ -634,15 +634,15 @@ Acts::DD4hepLayerBuilder::positiveLayers(const GeometryContext& gctx) const
 
         // set material surface
         if (layerPos == Acts::LayerMaterialPos::inner) {
-          innerBoundary->setAssociatedMaterial(materialProxy);
+          innerBoundary->assignSurfaceMaterial(materialProxy);
         }
 
         if (layerPos == Acts::LayerMaterialPos::outer) {
-          outerBoundary->setAssociatedMaterial(materialProxy);
+          outerBoundary->assignSurfaceMaterial(materialProxy);
         }
 
         if (layerPos == Acts::LayerMaterialPos::central) {
-          centralSurface->setAssociatedMaterial(materialProxy);
+          centralSurface->assignSurfaceMaterial(materialProxy);
         }
         // collect approach surfaces
         aSurfaces.push_back(innerBoundary);
@@ -702,7 +702,7 @@ Acts::DD4hepLayerBuilder::positiveLayers(const GeometryContext& gctx) const
         surfMaterial = std::make_shared<const HomogeneousSurfaceMaterial>(
             materialProperties);
       }
-      positiveLayer->surfaceRepresentation().setAssociatedMaterial(
+      positiveLayer->surfaceRepresentation().assignSurfaceMaterial(
           surfMaterial);
 
       // push back created layer