fix: Update volume material mapping (#1108)

This PR update the volume material mapping to solve 2 issue with it.

-First it had the same flaw as the one from the surface mapping (#1088)
-Second all the material hits where kept until the end of the mapping process, this meant that if the detector was too big and too many material track were used it could run out of memory. The are now progressively added to the grid during the mapping and averaged at the end.

Those change shouldn't affect the result mapping itself. (it was tested again with the ATLAS detector and appear to be fine)

Core/include/Acts/Material/MaterialGridHelper.hpp

@@ -94,30 +94,20 @@ Grid3D createGrid3D(
     const BinUtility& bins,
     std::function<Acts::Vector3(Acts::Vector3)>& transfoGlobalToLocal);
 
-/// @brief Concatenate a set of material at arbitrary space points on a set of
-/// grid points and produces a grid containing the averaged material values.
+/// @brief Average the material collected in a 2D grid and use it to create a 2D material grid
 ///
 /// @param [in] grid The material collecting grid
-/// @param [in] mPoints The set of material at the space points
-/// @param [in] transfoGlobalToLocal tranformation from local to local
 /// coordinate
 ///
 /// @return The average material grid decomposed into classification numbers
-MaterialGrid2D mapMaterialPoints(
-    Grid2D& grid, const Acts::RecordedMaterialVolumePoint& mPoints,
-    std::function<Acts::Vector2(Acts::Vector3)>& transfoGlobalToLocal);
+MaterialGrid2D mapMaterialPoints(Grid2D& grid);
 
-/// @brief Concatenate a set of material at arbitrary space points on a set of
-/// grid points and produces a grid containing the averaged material values.
+/// @brief Average the material collected in a 3D grid and use it to create a 3D material grid
 ///
 /// @param [in] grid The material collecting grid
-/// @param [in] mPoints The set of material at the space points
-/// @param [in] transfoGlobalToLocal tranformation from local to local
 /// coordinate
 ///
 /// @return The average material grid decomposed into classification numbers
-MaterialGrid3D mapMaterialPoints(
-    Grid3D& grid, const Acts::RecordedMaterialVolumePoint& mPoints,
-    std::function<Acts::Vector3(Acts::Vector3)>& transfoGlobalToLocal);
+MaterialGrid3D mapMaterialPoints(Grid3D& grid);
 
 }  // namespace Acts

Core/include/Acts/Material/VolumeMaterialMapper.hpp

@@ -16,6 +16,8 @@
 #include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "Acts/Geometry/TrackingVolume.hpp"
 #include "Acts/MagneticField/MagneticFieldContext.hpp"
+#include "Acts/Material/AccumulatedVolumeMaterial.hpp"
+#include "Acts/Material/MaterialGridHelper.hpp"
 #include "Acts/Material/MaterialSlab.hpp"
 #include "Acts/Propagator/MaterialInteractor.hpp"
 #include "Acts/Propagator/Navigator.hpp"
@@ -81,10 +83,28 @@ class VolumeMaterialMapper {
         : geoContext(gctx), magFieldContext(mctx) {}
 
     /// The recorded material per geometry ID
-    std::map<GeometryIdentifier, RecordedMaterialVolumePoint> recordedMaterial;
+    std::map<const GeometryIdentifier, Acts::AccumulatedVolumeMaterial>
+        homogeneousGrid;
 
-    /// The binning per geometry ID
-    std::map<GeometryIdentifier, BinUtility> materialBin;
+    /// The recorded 2D transform associated the grid for each geometry ID
+    std::map<const GeometryIdentifier,
+             std::function<Acts::Vector2(Acts::Vector3)>>
+        transform2D;
+
+    /// The 2D material grid for each geometry ID
+    std::map<const GeometryIdentifier, Grid2D> grid2D;
+
+    /// The recorded 3D transform associated the material grid for each geometry
+    /// ID
+    std::map<const GeometryIdentifier,
+             std::function<Acts::Vector3(Acts::Vector3)>>
+        transform3D;
+
+    /// The 3D material grid for each geometry ID
+    std::map<const GeometryIdentifier, Grid3D> grid3D;
+
+    /// The binning for each geometry ID
+    std::map<const GeometryIdentifier, BinUtility> materialBin;
 
     /// The surface material of the input tracking geometry
     std::map<GeometryIdentifier, std::shared_ptr<const ISurfaceMaterial>>
@@ -129,8 +149,8 @@ class VolumeMaterialMapper {
   /// @brief Method to finalize the maps
   ///
   /// It calls the final run averaging and then transforms
-  /// the AccumulatedVolume material class to a surface material
-  /// class type
+  /// the Homogeneous material into HomogeneousVolumeMaterial and
+  /// the 2D and 3D grid into a InterpolatedMaterialMap
   ///
   /// @param mState
   void finalizeMaps(State& mState) const;
@@ -179,15 +199,17 @@ class VolumeMaterialMapper {
   void collectMaterialSurfaces(State& mState,
                                const TrackingVolume& tVolume) const;
 
-  /// Create extra material point for the mapping
+  /// Create extra material point for the mapping and add them to the grid
   ///
-  /// @param matPoint RecordedMaterialVolumePoint where the extra hit are stored
+  /// @param mState The state to be filled
+  /// @param currentBinning a pair containing the current geometry ID and the current binning
   /// @param properties material properties of the original hit
   /// @param position position of the original hit
   /// @param direction direction of the track
-  void createExtraHits(RecordedMaterialVolumePoint& matPoint,
-                       Acts::MaterialSlab properties, Vector3 position,
-                       Vector3 direction) const;
+  void createExtraHits(
+      State& mState,
+      std::pair<const GeometryIdentifier, BinUtility>& currentBinning,
+      Acts::MaterialSlab properties, Vector3 position, Vector3 direction) const;
 
   /// Standard logger method
   const Logger& logger() const { return *m_logger; }

Core/src/Material/MaterialGridHelper.cpp

@@ -220,20 +220,7 @@ Acts::Grid3D Acts::createGrid3D(
                            std::move(gridAxis3)));
 }
 
-Acts::MaterialGrid2D Acts::mapMaterialPoints(
-    Acts::Grid2D& grid, const Acts::RecordedMaterialVolumePoint& mPoints,
-    std::function<Acts::Vector2(Acts::Vector3)>& transfoGlobalToLocal) {
-  // Walk over each properties
-  for (const auto& rm : mPoints) {
-    // Walk over each point associated with the properties
-    for (const auto& point : rm.second) {
-      // Search for fitting grid point and accumulate
-      Acts::Grid2D::index_t index =
-          grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(point));
-      grid.atLocalBins(index).accumulate(rm.first);
-    }
-  }
-
+Acts::MaterialGrid2D Acts::mapMaterialPoints(Acts::Grid2D& grid) {
   // Build material grid
   // Re-build the axes
   Acts::Grid2D::point_t min = grid.minPosition();
@@ -243,7 +230,7 @@ Acts::MaterialGrid2D Acts::mapMaterialPoints(
   Acts::EAxis axis1(min[0], max[0], nBins[0]);
   Acts::EAxis axis2(min[1], max[1], nBins[1]);
 
-  // Build the grid and fill it with data
+  // Fill the material Grid by averaging the material in the 2D grid
   Acts::MaterialGrid2D mGrid(std::make_tuple(axis1, axis2));
   for (size_t index = 0; index < grid.size(); index++) {
     mGrid.at(index) = grid.at(index).average().parameters();
@@ -252,20 +239,7 @@ Acts::MaterialGrid2D Acts::mapMaterialPoints(
   return mGrid;
 }
 
-Acts::MaterialGrid3D Acts::mapMaterialPoints(
-    Acts::Grid3D& grid, const Acts::RecordedMaterialVolumePoint& mPoints,
-    std::function<Acts::Vector3(Acts::Vector3)>& transfoGlobalToLocal) {
-  // Walk over each properties
-  for (const auto& rm : mPoints) {
-    // Walk over each point associated with the properties
-    for (const auto& point : rm.second) {
-      // Search for fitting grid point and accumulate
-      Acts::Grid3D::index_t index =
-          grid.localBinsFromLowerLeftEdge(transfoGlobalToLocal(point));
-      grid.atLocalBins(index).accumulate(rm.first);
-    }
-  }
-
+Acts::MaterialGrid3D Acts::mapMaterialPoints(Acts::Grid3D& grid) {
   // Build material grid
   // Re-build the axes
   Acts::Grid3D::point_t min = grid.minPosition();
@@ -276,7 +250,7 @@ Acts::MaterialGrid3D Acts::mapMaterialPoints(
   Acts::EAxis axis2(min[1], max[1], nBins[1]);
   Acts::EAxis axis3(min[2], max[2], nBins[2]);
 
-  // Build the grid and fill it with data
+  // Fill the material Grid by averaging the material in the 3D grid
   Acts::MaterialGrid3D mGrid(std::make_tuple(axis1, axis2, axis3));
   for (size_t index = 0; index < grid.size(); index++) {
     mGrid.at(index) = grid.at(index).average().parameters();