feat: make Geant4 surface mapper configurable (acts-project#3157)

This PR breaks the direct dependency of `ActsExamples::Geant4Simulation` from `TrackingGeometry` and allows to set a customised SensitiveSurface mapper.

The current SensitiveSurface mapper would fail for e.g. ITK geometry because the center comparison would probably not be valid. Also, it allows to test with Gen2 geometry, which will be needed (for the moment) for the `traccc` input simulation.

Examples/Algorithms/Digitization/include/ActsExamples/Digitization/DigitizationAlgorithm.hpp

@@ -41,11 +41,6 @@ namespace ActsFatras {
 class Barcode;
 }  // namespace ActsFatras
 
-namespace Acts {
-class Surface;
-class TrackingGeometry;
-}  // namespace Acts
-
 namespace ActsExamples {
 struct AlgorithmContext;
 

Examples/Algorithms/Geant4/include/ActsExamples/Geant4/Geant4Simulation.hpp

@@ -16,6 +16,7 @@
 #include "ActsExamples/Framework/IAlgorithm.hpp"
 #include "ActsExamples/Framework/ProcessCode.hpp"
 #include "ActsExamples/Framework/RandomNumbers.hpp"
+#include "ActsExamples/Geant4/SensitiveSurfaceMapper.hpp"
 
 #include <cstddef>
 #include <memory>
@@ -34,7 +35,6 @@ class G4VUserPhysicsList;
 class G4FieldManager;
 
 namespace Acts {
-class TrackingGeometry;
 class MagneticFieldProvider;
 class Volume;
 }  // namespace Acts
@@ -111,6 +111,9 @@ class Geant4SimulationBase : public IAlgorithm {
 /// Algorithm to run Geant4 simulation in the ActsExamples framework
 class Geant4Simulation final : public Geant4SimulationBase {
  public:
+  using SensitiveCandidates = std::function<std::vector<const Acts::Surface*>(
+      const Acts::GeometryContext&, const Acts::Vector3&)>;
+
   struct Config : public Geant4SimulationBase::Config {
     /// Name of the output collection : hits
     std::string outputSimHits = "simhits";
@@ -121,11 +124,12 @@ class Geant4Simulation final : public Geant4SimulationBase {
     /// Name of the output collection : final particles
     std::string outputParticlesFinal = "particles_final";
 
-    /// The ACTS tracking geometry
-    std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry;
+    /// The ACTS sensitive surfaces in a mapper, used for hit creation
+    std::shared_ptr<const SensitiveSurfaceMapper> sensitiveSurfaceMapper =
+        nullptr;
 
     /// The ACTS Magnetic field provider
-    std::shared_ptr<const Acts::MagneticFieldProvider> magneticField;
+    std::shared_ptr<const Acts::MagneticFieldProvider> magneticField = nullptr;
 
     /// If a physics list has to be instantiated this one is chosen.
     std::string physicsList = "FTFP_BERT";

Examples/Algorithms/Geant4/include/ActsExamples/Geant4/SensitiveSurfaceMapper.hpp

@@ -9,6 +9,9 @@
 #pragma once
 
 #include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Geometry/TrackingGeometry.hpp"
+#include "Acts/Utilities/GridAccessHelpers.hpp"
 #include "Acts/Utilities/Logger.hpp"
 
 #include <memory>
@@ -18,11 +21,41 @@
 class G4VPhysicalVolume;
 
 namespace Acts {
-class TrackingGeometry;
+class Surface;
 }
 
 namespace ActsExamples {
 
+// Helper struct to find the sensitive surface candidates
+struct SensitiveCandidates {
+  std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry = nullptr;
+  /// Find the sensitive surfaces for a given position
+  ///
+  /// This fulfills the concept of a SensitiveCandidates
+  ///
+  /// @param gctx the geometry context
+  /// @param position the position to look for sensitive surfaces
+  ///
+  /// @return a vector of sensitive surfaces
+  std::vector<const Acts::Surface*> operator()(
+      const Acts::GeometryContext& gctx, const Acts::Vector3& position) const {
+    std::vector<const Acts::Surface*> surfaces;
+
+    if (trackingGeometry != nullptr) {
+      auto layer = trackingGeometry->associatedLayer(gctx, position);
+
+      if (layer->surfaceArray() != nullptr) {
+        for (const auto& surface : layer->surfaceArray()->surfaces()) {
+          if (surface->associatedDetectorElement() != nullptr) {
+            surfaces.push_back(surface);
+          }
+        }
+      }
+    }
+    return surfaces;
+  }
+};
+
 /// This Mapper takes a (non-const) Geant4 geometry and maps
 /// it such that name will be containing the mapping prefix
 /// and the Acts::GeometryIdentifier of the surface.
@@ -34,6 +67,9 @@ namespace ActsExamples {
 /// elements of the Acts::TrackingGeoemtry w/o map lookup.
 class SensitiveSurfaceMapper {
  public:
+  using SensitiveCandidates = std::function<std::vector<const Acts::Surface*>(
+      const Acts::GeometryContext&, const Acts::Vector3&)>;
+
   constexpr static std::string_view mappingPrefix = "ActsGeoID#";
 
   /// Configuration struct for the surface mapper
@@ -44,8 +80,8 @@ class SensitiveSurfaceMapper {
     /// For which G4 volume names we try to find a mapping
     std::vector<std::string> volumeMappings;
 
-    /// The tracking geometry we try to map
-    std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry;
+    /// The sensitive surfaces that are being mapped to
+    SensitiveCandidates candidateSurfaces;
   };
 
   /// Constructor with:
@@ -62,10 +98,12 @@ class SensitiveSurfaceMapper {
   /// hierarchy and applies name remapping to the Physical volumes
   /// of the Geant4 geometry.
   ///
+  /// @param gctx the geometry context
   /// @param g4PhysicalVolume the current physical volume in process
   /// @param motherPosition the absolute position of the mother
   /// @param sCounter  a counter of how many volumes have been remapped
-  void remapSensitiveNames(G4VPhysicalVolume* g4PhysicalVolume,
+  void remapSensitiveNames(const Acts::GeometryContext& gctx,
+                           G4VPhysicalVolume* g4PhysicalVolume,
                            const Acts::Transform3& motherTransform,
                            int& sCounter) const;
 

Examples/Algorithms/Geant4/src/Geant4Simulation.cpp

@@ -269,22 +269,14 @@ ActsExamples::Geant4Simulation::Geant4Simulation(const Config& cfg,
     g4World->GetLogicalVolume()->SetFieldManager(m_fieldManager.get(), true);
   }
 
-  // An ACTS TrackingGeometry is provided, so simulation for sensitive
-  // detectors is turned on - they need to get matched first
-  if (cfg.trackingGeometry) {
+  // ACTS sensitive surfaces are provided, so hit creation is turned on
+  if (cfg.sensitiveSurfaceMapper != nullptr) {
     ACTS_INFO(
         "Remapping selected volumes from Geant4 to Acts::Surface::GeometryID");
-
-    SensitiveSurfaceMapper::Config ssmCfg;
-    ssmCfg.trackingGeometry = cfg.trackingGeometry;
-    ssmCfg.volumeMappings = cfg.volumeMappings;
-    ssmCfg.materialMappings = cfg.materialMappings;
-
-    SensitiveSurfaceMapper sensitiveSurfaceMapper(
-        ssmCfg, m_logger->cloneWithSuffix("SensitiveSurfaceMapper"));
     int sCounter = 0;
-    sensitiveSurfaceMapper.remapSensitiveNames(
-        g4World, Acts::Transform3::Identity(), sCounter);
+    cfg.sensitiveSurfaceMapper->remapSensitiveNames(
+        Acts::GeometryContext{}, g4World, Acts::Transform3::Identity(),
+        sCounter);
 
     ACTS_INFO("Remapping successful for " << sCounter << " selected volumes.");
   }

Examples/Algorithms/Geant4/src/SensitiveSurfaceMapper.cpp

@@ -9,12 +9,6 @@
 #include "ActsExamples/Geant4/SensitiveSurfaceMapper.hpp"
 
 #include "Acts/Definitions/Units.hpp"
-#include "Acts/Geometry/GeometryContext.hpp"
-#include "Acts/Geometry/GeometryIdentifier.hpp"
-#include "Acts/Geometry/Layer.hpp"
-#include "Acts/Geometry/TrackingGeometry.hpp"
-#include "Acts/Surfaces/Surface.hpp"
-#include "Acts/Surfaces/SurfaceArray.hpp"
 
 #include <algorithm>
 #include <ostream>
@@ -27,14 +21,10 @@
 
 ActsExamples::SensitiveSurfaceMapper::SensitiveSurfaceMapper(
     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
-    : m_cfg(cfg), m_logger(std::move(logger)) {
-  if (m_cfg.trackingGeometry == nullptr) {
-    throw std::invalid_argument("No Acts::TrackingGeometry provided.");
-  }
-}
+    : m_cfg(cfg), m_logger(std::move(logger)) {}
 
 void ActsExamples::SensitiveSurfaceMapper::remapSensitiveNames(
-    G4VPhysicalVolume* g4PhysicalVolume,
+    const Acts::GeometryContext& gctx, G4VPhysicalVolume* g4PhysicalVolume,
     const Acts::Transform3& motherTransform, int& sCounter) const {
   constexpr double convertLength = CLHEP::mm / Acts::UnitConstants::mm;
 
@@ -63,7 +53,7 @@ void ActsExamples::SensitiveSurfaceMapper::remapSensitiveNames(
   if (G4int nDaughters = g4LogicalVolume->GetNoDaughters(); nDaughters > 0) {
     // Step down to all daughters
     for (G4int id = 0; id < nDaughters; ++id) {
-      remapSensitiveNames(g4LogicalVolume->GetDaughter(id), transform,
+      remapSensitiveNames(gctx, g4LogicalVolume->GetDaughter(id), transform,
                           sCounter);
     }
     return;
@@ -81,36 +71,20 @@ void ActsExamples::SensitiveSurfaceMapper::remapSensitiveNames(
                 volumeName) != m_cfg.volumeMappings.end();
 
   if (isSensitive || isMappedMaterial || isMappedVolume) {
-    // Find the associated ACTS object
-    auto actsLayer = m_cfg.trackingGeometry->associatedLayer(
-        Acts::GeometryContext(), g4AbsPosition);
-
     // Prepare the mapped surface
     const Acts::Surface* mappedSurface = nullptr;
 
-    if (actsLayer != nullptr && actsLayer->surfaceArray() != nullptr) {
-      auto actsSurfaces = actsLayer->surfaceArray()->at(g4AbsPosition);
-      if (!actsSurfaces.empty()) {
-        // Fast matching: search
-        for (const auto& as : actsSurfaces) {
-          if (as->center(Acts::GeometryContext()).isApprox(g4AbsPosition)) {
-            mappedSurface = as;
-            break;
-          }
-        }
-      }
-      if (mappedSurface == nullptr) {
-        // Slow matching: Fallback, loop over all layer surfaces
-        for (const auto& as : actsLayer->surfaceArray()->surfaces()) {
-          if (as->center(Acts::GeometryContext()).isApprox(g4AbsPosition)) {
-            mappedSurface = as;
-            break;
-          }
-        }
+    // Get the candidate surfaces
+    auto candidateSurfaces = m_cfg.candidateSurfaces(gctx, g4AbsPosition);
+    for (const auto& candidateSurface : candidateSurfaces) {
+      // Check for center matching at the moment (needs to be extended)
+      if (candidateSurface->center(gctx).isApprox(g4AbsPosition)) {
+        mappedSurface = candidateSurface;
+        break;
       }
     }
-    // A mapped surface was found, a new name will be set that
-    // contains the GeometryID/
+
+    // A mapped surface was found, a new name will be set that G4PhysVolume
     if (mappedSurface != nullptr) {
       ++sCounter;
       std::string mappedVolumeName(SensitiveSurfaceMapper::mappingPrefix);

Examples/Python/python/acts/examples/simulation.py

@@ -616,7 +616,7 @@ def getG4DetectorConstructionFactory(
 def addGeant4(
     s: acts.examples.Sequencer,
     detector: Optional[Any],
-    trackingGeometry: acts.TrackingGeometry,
+    trackingGeometry: Union[acts.TrackingGeometry, acts.Detector],
     field: acts.MagneticFieldProvider,
     rnd: acts.examples.RandomNumbers,
     g4DetectorConstructionFactory: Optional[Any] = None,
@@ -667,7 +667,7 @@ def addGeant4(
         if given, secondary particles are removed from simulation
     """
 
-    from acts.examples.geant4 import Geant4Simulation
+    from acts.examples.geant4 import Geant4Simulation, SensitiveSurfaceMapper
 
     customLogLevel = acts.examples.defaultLogging(s, logLevel)
 
@@ -692,6 +692,13 @@ def addGeant4(
 
     global __geant4Handle
 
+    smmConfig = SensitiveSurfaceMapper.Config()
+    smmConfig.volumeMappings = volumeMappings
+    smmConfig.materialMappings = materialMappings
+    sensitiveMapper = SensitiveSurfaceMapper.create(
+        smmConfig, acts.logging.INFO, trackingGeometry
+    )
+
     # Simulation
     alg = Geant4Simulation(
         level=customLogLevel(),
@@ -702,7 +709,7 @@ def addGeant4(
         outputParticlesInitial=outputParticlesInitial,
         outputParticlesFinal=outputParticlesFinal,
         outputSimHits=outputSimHits,
-        trackingGeometry=trackingGeometry,
+        sensitiveSurfaceMapper=sensitiveMapper,
         magneticField=field,
         physicsList=physicsList,
         volumeMappings=volumeMappings,

Examples/Python/src/Geant4Component.cpp

@@ -22,6 +22,7 @@
 #include "ActsExamples/Geant4/Geant4Manager.hpp"
 #include "ActsExamples/Geant4/Geant4Simulation.hpp"
 #include "ActsExamples/Geant4/RegionCreator.hpp"
+#include "ActsExamples/Geant4/SensitiveSurfaceMapper.hpp"
 #include "ActsExamples/Geant4Detector/Geant4Detector.hpp"
 #include "ActsExamples/MuonSpectrometerMockupDetector/MockupSectorBuilder.hpp"
 #include "ActsExamples/TelescopeDetector/TelescopeDetector.hpp"
@@ -97,6 +98,66 @@ PYBIND11_MODULE(ActsPythonBindingsGeant4, mod) {
     ACTS_PYTHON_STRUCT_END();
   }
 
+  {
+    using Config = SensitiveSurfaceMapper::Config;
+    auto sm =
+        py::class_<SensitiveSurfaceMapper,
+                   std::shared_ptr<SensitiveSurfaceMapper>>(
+            mod, "SensitiveSurfaceMapper")
+            .def(py::init([](const Config& cfg, Acts::Logging::Level level) {
+              return std::make_shared<SensitiveSurfaceMapper>(
+                  cfg, getDefaultLogger("SensitiveSurfaceMapper", level));
+            }));
+
+    auto c = py::class_<Config>(sm, "Config").def(py::init<>());
+    ACTS_PYTHON_STRUCT_BEGIN(c, Config);
+    ACTS_PYTHON_MEMBER(materialMappings);
+    ACTS_PYTHON_MEMBER(volumeMappings);
+    ACTS_PYTHON_MEMBER(candidateSurfaces);
+    ACTS_PYTHON_STRUCT_END();
+
+    sm.def(
+        "create", [](const Config& cfg, Acts::Logging::Level level,
+                     const std::shared_ptr<const TrackingGeometry> tGeometry) {
+          // Set a new surface finder
+          Config ccfg = cfg;
+          ccfg.candidateSurfaces = ActsExamples::SensitiveCandidates{tGeometry};
+          return std::make_shared<SensitiveSurfaceMapper>(
+              ccfg, getDefaultLogger("SensitiveSurfaceMapper", level));
+        });
+
+    sm.def("create",
+           [](const Config& cfg, Acts::Logging::Level level,
+              const std::shared_ptr<const Experimental::Detector>& detector) {
+             // Helper struct to find the sensitive surface candidates
+             struct SensitiveCandidates {
+               std::shared_ptr<const Experimental::Detector> detector;
+
+               /// Find the sensitive surfaces for a given position
+               std::vector<const Acts::Surface*> operator()(
+                   const Acts::GeometryContext& gctx,
+                   const Acts::Vector3& position) const {
+                 std::vector<const Acts::Surface*> surfaces;
+                 // Here's the detector volume
+                 auto volume = detector->findDetectorVolume(gctx, position);
+                 if (volume != nullptr) {
+                   for (const auto& surface : volume->surfaces()) {
+                     if (surface->associatedDetectorElement() != nullptr) {
+                       surfaces.push_back(surface);
+                     }
+                   }
+                 }
+                 return surfaces;
+               }
+             };
+             // Set a new surface finder
+             Config ccfg = cfg;
+             ccfg.candidateSurfaces = SensitiveCandidates{detector};
+             return std::make_shared<SensitiveSurfaceMapper>(
+                 ccfg, getDefaultLogger("SensitiveSurfaceMapper", level));
+           });
+  }
+
   {
     using Algorithm = Geant4Simulation;
     using Config = Algorithm::Config;
@@ -114,7 +175,7 @@ PYBIND11_MODULE(ActsPythonBindingsGeant4, mod) {
     ACTS_PYTHON_MEMBER(outputSimHits);
     ACTS_PYTHON_MEMBER(outputParticlesInitial);
     ACTS_PYTHON_MEMBER(outputParticlesFinal);
-    ACTS_PYTHON_MEMBER(trackingGeometry);
+    ACTS_PYTHON_MEMBER(sensitiveSurfaceMapper);
     ACTS_PYTHON_MEMBER(magneticField);
     ACTS_PYTHON_MEMBER(physicsList);
     ACTS_PYTHON_MEMBER(volumeMappings);