feat: introduce material accumulation interface (MM3) (acts-project#3020

)

This is the first PR in a series that divides the Material Mapping into
logical, unit testable modules:

1) Finding intersections with surfaces and associations to volumes (MM2)
2) Assigning material interactions to those intersections (MM1)
3) Mapping those onto dedicated Surface / Volume Material Mappers (MM3,
this PR, partly)
4) Steer that by a chained algorithm (MM4)

This PR:

It encapsulates the pure accumulation of the material for Surface based
materials into a new module. The accumulation is now decoupled from the
association, i.e. one can run the mapping with a navigator or with a
trial and error method, and still use the exact same material mapper.
The code is largely transferred from the 'SurfaceMaterialMapper' which
will vanish to exist after this re-organization.

@Corentin-Allaire - the State of this Mapper would have all the access
to run the optimisation, I have not yet re-introduced the track
variance, but that should be a five-line change.

This PR is blocked by acts-project#3015 and acts-project#3016.


All cases are showcased and tested in a set of UnitTests.

Core/include/Acts/Material/BinnedSurfaceMaterialAccumulater.hpp

@@ -0,0 +1,88 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// 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/.
+
+#pragma once
+
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Material/AccumulatedSurfaceMaterial.hpp"
+#include "Acts/Material/interface/ISurfaceMaterialAccumulater.hpp"
+#include "Acts/Utilities/Logger.hpp"
+
+namespace Acts {
+
+/// @brief The binned surface material accumulater
+///
+/// It consumes the assigned material interactions and then accumulates
+/// the material on the surfaces in prepared binned containers for averaging
+
+class BinnedSurfaceMaterialAccumulater final
+    : public ISurfaceMaterialAccumulater {
+ public:
+  /// @brief Nested config struct
+  struct Config {
+    GeometryContext geoContext;
+
+    /// Correct for empty bins (recommended)
+    bool emptyBinCorrection = true;
+
+    /// The surfaces to be used for the accummulation
+    std::vector<const Surface*> materialSurfaces = {};
+  };
+
+  /// @brief Nested state struct
+  struct State final : public ISurfaceMaterialAccumulater::State {
+    /// The accumulated material per geometry ID
+    std::map<GeometryIdentifier, AccumulatedSurfaceMaterial>
+        accumulatedMaterial;
+  };
+
+  /// Constructor
+  ///
+  /// @param cfg the configuration struct
+  /// @param mlogger the logger
+  BinnedSurfaceMaterialAccumulater(
+      const Config& cfg,
+      std::unique_ptr<const Logger> mlogger =
+          getDefaultLogger("BinnedSurfaceMaterialAccumulater", Logging::INFO));
+
+  /// Factory for creating the state
+  std::unique_ptr<ISurfaceMaterialAccumulater::State> createState()
+      const override;
+
+  /// @brief Accumulate the material interaction on the surface
+  ///
+  /// @param state is the state of the accumulater
+  /// @param interactions is the material interactions, with assigned surfaces
+  /// @param surfacesWithoutAssignment are the surfaces without assignment
+  ///
+  /// @note this the track average over the binned material
+  void accumulate(ISurfaceMaterialAccumulater::State& state,
+                  const std::vector<MaterialInteraction>& interactions,
+                  const std::vector<IAssignmentFinder::SurfaceAssignment>&
+                      surfacesWithoutAssignment) const override;
+
+  /// Finalize the surface material maps
+  ///
+  /// @param state the state of the accumulator
+  ///
+  /// @note this does the run average over the (binned) material
+  std::map<GeometryIdentifier, std::shared_ptr<const ISurfaceMaterial>>
+  finalizeMaterial(ISurfaceMaterialAccumulater::State& state) const override;
+
+ private:
+  /// Access method to the logger
+  const Logger& logger() const { return *m_logger; }
+
+  /// The configuration
+  Config m_cfg;
+
+  /// The logger
+  std::unique_ptr<const Logger> m_logger;
+};
+
+}  // namespace Acts

Core/include/Acts/Material/MaterialInteractionAssignment.hpp

@@ -12,6 +12,7 @@
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/Geometry/GeometryHierarchyMap.hpp"
 #include "Acts/Material/MaterialInteraction.hpp"
+#include "Acts/Material/interface/IAssignmentFinder.hpp"
 
 #include <array>
 #include <functional>
@@ -31,6 +32,8 @@ struct Result {
   std::vector<MaterialInteraction> assigned = {};
   /// The unassigned material interactions
   std::vector<MaterialInteraction> unassigned = {};
+  /// Surfaces without assignment (for empty hit correction)
+  std::vector<IAssignmentFinder::SurfaceAssignment> unassignedSurfaces = {};
 };
 
 /// @brief definition of a global veto for assigning material interactions
@@ -55,7 +58,7 @@ using GlobalVeto =
 /// @param suggestedAssignment is the suggested assignment: surface, position, direction
 using LocalVeto = std::function<bool(
     const MaterialInteraction&,
-    const std::tuple<const Surface*, Vector3, Vector3>& suggestedAssignment)>;
+    const IAssignmentFinder::SurfaceAssignment& suggestedAssignment)>;
 
 /// @brief definition of possible re-assignments to next surface, this could e.g.
 /// be used for respecting pre/post mapping directives that are not fully
@@ -72,8 +75,8 @@ using LocalVeto = std::function<bool(
 /// @note this changes the MaterialInteraction if the re-assignment is accepted
 using ReAssignment = std::function<void(
     MaterialInteraction& materialInteraction,
-    const std::tuple<const Surface*, Vector3, Vector3>& suggestedAssignment,
-    const std::tuple<const Surface*, Vector3, Vector3>& suggestedReAssignment)>;
+    const IAssignmentFinder::SurfaceAssignment& suggestedAssignment,
+    const IAssignmentFinder::SurfaceAssignment& suggestedReAssignment)>;
 
 /// @brief Options for the material interaction matcher
 /// The options are used to specify the vetos for the assignment
@@ -101,7 +104,7 @@ struct Options {
 /// @return a pair of vectors of assigned and unassigned material interactions
 Result assign(const GeometryContext& gctx,
               const std::vector<MaterialInteraction>& materialInteractions,
-              const std::vector<std::tuple<const Surface*, Vector3, Vector3>>&
+              const std::vector<IAssignmentFinder::SurfaceAssignment>&
                   intersectedSurfaces,
               const Options& options = Options());
 

Core/include/Acts/Material/interface/IAssignmentFinder.hpp

@@ -6,6 +6,8 @@
 // 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/.
 
+#pragma once
+
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/MagneticField/MagneticFieldContext.hpp"
 #include "Acts/Material/MaterialInteraction.hpp"
@@ -14,8 +16,6 @@
 #include <utility>
 #include <vector>
 
-#pragma once
-
 namespace Acts {
 
 /// @brief Interface for the material mapping that seeks the possible

Core/include/Acts/Material/interface/ISurfaceMaterialAccumulater.hpp

@@ -0,0 +1,62 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 CERN for the benefit of the Acts project
+//
+// 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/.
+
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "Acts/Material/ISurfaceMaterial.hpp"
+#include "Acts/Material/MaterialInteraction.hpp"
+#include "Acts/Material/interface/IAssignmentFinder.hpp"
+
+#include <map>
+#include <memory>
+#include <vector>
+
+namespace Acts {
+
+class Surface;
+
+/// @brief Interface for the material mapping, this is the accumulation step
+class ISurfaceMaterialAccumulater {
+ public:
+  /// The state of the material accumulater, this is used
+  /// to cache information across tracks/events
+  class State {
+   public:
+    virtual ~State() = default;
+  };
+
+  /// Virtual destructor
+  virtual ~ISurfaceMaterialAccumulater() = default;
+
+  /// Factory for creating the state
+  virtual std::unique_ptr<State> createState() const = 0;
+
+  /// @brief Accumulate the material interaction on the surface
+  ///
+  /// @param state is the state of the accumulater
+  /// @param interactions is the material interactions, with assigned surfaces
+  /// @param surfacesWithoutAssignment are the surfaces without assignment
+  ///
+  /// @note this the track average over the binned material
+  virtual void accumulate(
+      State& state, const std::vector<MaterialInteraction>& interactions,
+      const std::vector<IAssignmentFinder::SurfaceAssignment>&
+          surfacesWithoutAssignment) const = 0;
+
+  /// Finalize the surface material maps
+  ///
+  /// @param state the state of the accumulator
+  ///
+  /// @note this does the run average over the (binned) material
+  virtual std::map<GeometryIdentifier, std::shared_ptr<const ISurfaceMaterial>>
+  finalizeMaterial(State& state) const = 0;
+};
+
+}  // namespace Acts

Core/include/Acts/Utilities/BinAdjustment.hpp

@@ -31,8 +31,9 @@ namespace Acts {
 /// @param transform Transform for the adjusted @c BinUtility
 ///
 /// @return new updated BinUtiltiy
-BinUtility adjustBinUtility(const BinUtility& bu, const RadialBounds& rBounds,
-                            const Transform3& transform) {
+static inline BinUtility adjustBinUtility(const BinUtility& bu,
+                                          const RadialBounds& rBounds,
+                                          const Transform3& transform) {
   // Default constructor
   BinUtility uBinUtil(transform);
 
@@ -80,8 +81,9 @@ BinUtility adjustBinUtility(const BinUtility& bu, const RadialBounds& rBounds,
 /// @param transform Transform for the adjusted @c BinUtility
 ///
 /// @return new updated BinUtiltiy
-BinUtility adjustBinUtility(const BinUtility& bu, const CylinderBounds& cBounds,
-                            const Transform3& transform) {
+static inline BinUtility adjustBinUtility(const BinUtility& bu,
+                                          const CylinderBounds& cBounds,
+                                          const Transform3& transform) {
   // Default constructor
   BinUtility uBinUtil(transform);
 
@@ -133,9 +135,9 @@ BinUtility adjustBinUtility(const BinUtility& bu, const CylinderBounds& cBounds,
 /// @param transform Transform for the adjusted @c BinUtility
 ///
 /// @return new updated BinUtiltiy
-BinUtility adjustBinUtility(const BinUtility& bu,
-                            const RectangleBounds& pBounds,
-                            const Transform3& transform) {
+static inline BinUtility adjustBinUtility(const BinUtility& bu,
+                                          const RectangleBounds& pBounds,
+                                          const Transform3& transform) {
   // Default constructor
   BinUtility uBinUtil(transform);
 
@@ -183,9 +185,9 @@ BinUtility adjustBinUtility(const BinUtility& bu,
 /// @param transform Transform for the adjusted @c BinUtility
 ///
 /// @return new updated BinUtiltiy
-BinUtility adjustBinUtility(const BinUtility& bu,
-                            const TrapezoidBounds& pBounds,
-                            const Transform3& transform) {
+static inline BinUtility adjustBinUtility(const BinUtility& bu,
+                                          const TrapezoidBounds& pBounds,
+                                          const Transform3& transform) {
   // Default constructor
   BinUtility uBinUtil(transform);
 
@@ -233,8 +235,9 @@ BinUtility adjustBinUtility(const BinUtility& bu,
 /// @param gctx Geometry context to get the surfaces transform
 ///
 /// @return new updated BinUtiltiy
-BinUtility adjustBinUtility(const BinUtility& bu, const Surface& surface,
-                            const GeometryContext& gctx) {
+static inline BinUtility adjustBinUtility(const BinUtility& bu,
+                                          const Surface& surface,
+                                          const GeometryContext& gctx) {
   // The surface type is a cylinder
   if (surface.type() == Surface::Cylinder) {
     // Cast to Cylinder bounds and return