feat: example material mapping using new core mapper (MM5)

Core/include/Acts/Material/IntersectionMaterialAssigner.hpp

@@ -44,7 +44,14 @@ class IntersectionMaterialAssigner final : public IAssignmentFinder {
   };
 
   /// @brief Construct with the configuration
-  IntersectionMaterialAssigner(const Config& cfg) : m_cfg(cfg) {}
+  ///
+  /// @param cfg is the configuration struct
+  /// @param mlogger is the logger
+  IntersectionMaterialAssigner(
+      const Config& cfg,
+      std::unique_ptr<const Logger> mlogger =
+          getDefaultLogger("IntersectionMaterialAssigner", Logging::INFO))
+      : m_cfg(cfg), m_logger(std::move(mlogger)) {}
 
   /// @brief Method for generating assignment candidates for the
   /// material interaction assignment to surfaces or volumes
@@ -63,8 +70,14 @@ class IntersectionMaterialAssigner final : public IAssignmentFinder {
                        const Vector3& direction) const final;
 
  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/src/Material/IntersectionMaterialAssigner.cpp

@@ -10,6 +10,7 @@
 
 #include "Acts/Surfaces/BoundaryCheck.hpp"
 #include "Acts/Surfaces/Surface.hpp"
+#include "Acts/Utilities/StringHelpers.hpp"
 
 namespace {
 
@@ -51,6 +52,9 @@ Acts::IntersectionMaterialAssigner::assignmentCandidates(
             std::vector<Acts::IAssignmentFinder::VolumeAssignment>>
       candidates;
 
+  ACTS_DEBUG("Finding material assignment from position "
+             << toString(position) << " and direction " << toString(direction));
+
   // Try the surfaces first
   auto sIntersections =
       forwardOrderedIntersections(gctx, position, direction, m_cfg.surfaces);
@@ -102,6 +106,9 @@ Acts::IntersectionMaterialAssigner::assignmentCandidates(
     }
   }
 
+  ACTS_DEBUG("Found " << candidates.first.size() << " surface candidates and "
+                      << candidates.second.size() << " volume candidates");
+
   // Return the result
   return candidates;
 }

Examples/Algorithms/MaterialMapping/CMakeLists.txt

@@ -1,6 +1,7 @@
 add_library(
   ActsExamplesMaterialMapping SHARED
-  src/MaterialMapping.cpp)
+  src/MaterialMapping.cpp
+  src/CoreMaterialMapping.cpp)
 target_include_directories(
   ActsExamplesMaterialMapping
   PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>)

Examples/Algorithms/MaterialMapping/include/ActsExamples/MaterialMapping/CoreMaterialMapping.hpp

@@ -0,0 +1,89 @@
+// 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/MaterialMapper.hpp"
+#include "Acts/Utilities/Logger.hpp"
+#include "ActsExamples/Framework/DataHandle.hpp"
+#include "ActsExamples/Framework/IAlgorithm.hpp"
+#include "ActsExamples/Framework/ProcessCode.hpp"
+#include "ActsExamples/MaterialMapping/IMaterialWriter.hpp"
+
+namespace ActsExamples {
+
+/// @class CoreMaterialMapping
+///
+/// @brief Initiates and executes material mapping using the MaterialMapper
+/// from the core component of ACTS
+///
+/// By construction, the material mapping needs inter-event information
+/// to build the material maps of accumulated single particle views.
+/// However, running it in one single event, puts enormous pressure onto
+/// the I/O structure.
+///
+/// It therefore saves the mapping state/cache as a private member variable
+/// and is designed to be executed in a single threaded mode.
+class CoreMaterialMapping : public IAlgorithm {
+ public:
+  /// @class nested Config class
+  /// of the MaterialMapping algorithm
+  struct Config {
+    /// Input collection
+    std::string inputMaterialTracks = "material_tracks";
+
+    /// The actually mapped material tracks
+    std::string mappedMaterialTracks = "mapped_material_tracks";
+
+    /// Theunmapped part of the material tracks
+    std::string unmappedMaterialTracks = "unmapped_material_tracks";
+
+    /// The ACTS material mapper from the core component
+    std::shared_ptr<Acts::MaterialMapper> materialMapper = nullptr;
+
+    /// The writer of the material
+    std::vector<std::shared_ptr<IMaterialWriter>> materiaMaplWriters{};
+  };
+
+  /// Constructor
+  ///
+  /// @param cfg The configuration struct carrying the used tools
+  /// @param level The output logging level
+  CoreMaterialMapping(const Config& cfg,
+                      Acts::Logging::Level level = Acts::Logging::INFO);
+
+  /// Destructor
+  /// - it also writes out the file
+  ~CoreMaterialMapping() override;
+
+  /// Framework execute method
+  ///
+  /// @param context The algorithm context for event consistency
+  ActsExamples::ProcessCode execute(
+      const AlgorithmContext& context) const override;
+
+  /// Readonly access to the config
+  const Config& config() const { return m_cfg; }
+
+ private:
+  Config m_cfg;  //!< internal config object
+
+  std::unique_ptr<Acts::MaterialMapper::State> m_mappingState{nullptr};
+
+  ReadDataHandle<std::unordered_map<std::size_t, Acts::RecordedMaterialTrack>>
+      m_inputMaterialTracks{this, "InputMaterialTracks"};
+
+  WriteDataHandle<std::unordered_map<std::size_t, Acts::RecordedMaterialTrack>>
+      m_outputMappedMaterialTracks{this, "OutputMappedMaterialTracks"};
+
+  WriteDataHandle<std::unordered_map<std::size_t, Acts::RecordedMaterialTrack>>
+      m_outputUnmappedMaterialTracks{this, "OutputUnmappedMaterialTracks"};
+};
+
+}  // namespace ActsExamples

Examples/Algorithms/MaterialMapping/src/CoreMaterialMapping.cpp

@@ -0,0 +1,81 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2017-2020 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/.
+
+#include "ActsExamples/MaterialMapping/CoreMaterialMapping.hpp"
+
+#include "Acts/Material/AccumulatedMaterialSlab.hpp"
+#include "Acts/Material/AccumulatedSurfaceMaterial.hpp"
+#include "ActsExamples/MaterialMapping/IMaterialWriter.hpp"
+
+#include <stdexcept>
+#include <unordered_map>
+
+namespace ActsExamples {
+
+CoreMaterialMapping::CoreMaterialMapping(const CoreMaterialMapping::Config& cfg,
+                                         Acts::Logging::Level level)
+    : IAlgorithm("CoreMaterialMapping", level), m_cfg(cfg) {
+  // Prepare the I/O collections
+  m_inputMaterialTracks.initialize(m_cfg.inputMaterialTracks);
+  m_outputMappedMaterialTracks.initialize(m_cfg.mappedMaterialTracks);
+  m_outputUnmappedMaterialTracks.initialize(m_cfg.unmappedMaterialTracks);
+
+  ACTS_INFO("This algorithm requires inter-event information, "
+            << "run in single-threaded mode!");
+
+  if (m_cfg.materialMapper == nullptr) {
+    throw std::invalid_argument("Missing material mapper");
+  }
+  // Create the state object
+  m_mappingState = m_cfg.materialMapper->createState();
+}
+
+CoreMaterialMapping::~CoreMaterialMapping() {
+  Acts::DetectorMaterialMaps detectorMaterial =
+      m_cfg.materialMapper->finalizeMaps(*m_mappingState);
+  // Loop over the available writers and write the maps
+  for (auto& imw : m_cfg.materiaMaplWriters) {
+    imw->writeMaterial(detectorMaterial);
+  }
+}
+
+ProcessCode CoreMaterialMapping::execute(
+    const AlgorithmContext& context) const {
+  // Take the collection from the EventStore: input collection
+  std::unordered_map<std::size_t, Acts::RecordedMaterialTrack>
+      mtrackCollection = m_inputMaterialTracks(context);
+
+  // Write the output collections to the Event store : mapped and unmapped
+  std::unordered_map<std::size_t, Acts::RecordedMaterialTrack>
+      mappedTrackCollection;
+
+  std::unordered_map<std::size_t, Acts::RecordedMaterialTrack>
+      unmappedTrackCollection;
+
+  // To make it work with the framework needs a lock guard
+  auto mappingState =
+      const_cast<Acts::MaterialMapper::State*>(m_mappingState.get());
+
+  for (auto& [idTrack, mTrack] : mtrackCollection) {
+    auto [mapped, unmapped] = m_cfg.materialMapper->mapMaterial(
+        *mappingState, context.geoContext, context.magFieldContext, mTrack);
+
+    mappedTrackCollection.emplace_hint(mappedTrackCollection.end(), idTrack,
+                                       mapped);
+    unmappedTrackCollection.emplace_hint(unmappedTrackCollection.end(), idTrack,
+                                         unmapped);
+  }
+
+  // Write the mapped and unmapped material tracks to the output
+  m_outputMappedMaterialTracks(context, std::move(mappedTrackCollection));
+  m_outputUnmappedMaterialTracks(context, std::move(unmappedTrackCollection));
+
+  return ProcessCode::SUCCESS;
+}
+
+}  // namespace ActsExamples

Examples/Python/src/Geometry.cpp

@@ -59,6 +59,21 @@ struct GeometryIdentifierHookBinding : public Acts::GeometryIdentifierHook {
         .cast<Acts::GeometryIdentifier>();
   }
 };
+
+struct MaterialSurfaceSelector {
+  std::vector<const Acts::Surface*> surfaces = {};
+
+  /// @param surface is the test surface
+  void operator()(const Acts::Surface* surface) {
+    if (surface->surfaceMaterial() != nullptr) {
+      if (std::find(surfaces.begin(), surfaces.end(), surface) ==
+          surfaces.end()) {
+        surfaces.push_back(surface);
+      }
+    }
+  }
+};
+
 }  // namespace
 
 namespace Acts::Python {
@@ -126,6 +141,12 @@ void addGeometry(Context& ctx) {
              [](Acts::TrackingGeometry& self, py::function& func) {
                self.visitSurfaces(func);
              })
+        .def("extractMaterialSurfaces",
+             [](Acts::TrackingGeometry& self) {
+               MaterialSurfaceSelector selector;
+               self.visitSurfaces(selector, false);
+               return selector.surfaces;
+             })
         .def_property_readonly(
             "worldVolume",
             &Acts::TrackingGeometry::highestTrackingVolumeShared);
@@ -190,8 +211,13 @@ void addExperimentalGeometry(Context& ctx) {
 
   // Detector definition
   py::class_<Detector, std::shared_ptr<Detector>>(m, "Detector")
-      .def("number_volumes",
-           [](Detector& self) { return self.volumes().size(); });
+      .def("numberVolumes",
+           [](Detector& self) { return self.volumes().size(); })
+      .def("extractMaterialSurfaces", [](Detector& self) {
+        MaterialSurfaceSelector selector;
+        self.visitSurfaces(selector);
+        return selector.surfaces;
+      });
 
   // Portal definition
   py::class_<Portal, std::shared_ptr<Portal>>(m, "Portal");