feat: Robust surface matching for Geant4

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

@@ -16,8 +16,13 @@
 
 #include <G4UserSteppingAction.hh>
 
+class G4VPhysicalVolume;
 class G4Step;
 
+namespace Acts {
+class Surface;
+}
+
 namespace ActsExamples {
 
 /// The G4SteppingAction that is called for every step in
@@ -52,6 +57,15 @@ class SensitiveSteppingAction : public G4UserSteppingAction {
   /// @param step is the Geant4 step of the particle
   void UserSteppingAction(const G4Step* step) override;
 
+  /// Set the multimap that correlates G4VPhysicalVolumes to Acts::Surfaces
+  ///
+  /// @param surfaceMapping the multimap of physical volumes to surfaces
+  void assignSurfaceMapping(
+      const std::multimap<const G4VPhysicalVolume*, const Acts::Surface*>&
+          surfaceMapping) {
+    m_surfaceMapping = std::move(surfaceMapping);
+  }
+
  protected:
   Config m_cfg;
 
@@ -64,6 +78,9 @@ class SensitiveSteppingAction : public G4UserSteppingAction {
 
   /// The looging instance
   std::unique_ptr<const Acts::Logger> m_logger;
+
+  std::multimap<const G4VPhysicalVolume*, const Acts::Surface*>
+      m_surfaceMapping;
 };
 
 }  // namespace ActsExamples

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

@@ -10,10 +10,12 @@
 
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Geometry/GeometryIdentifier.hpp"
 #include "Acts/Geometry/TrackingGeometry.hpp"
 #include "Acts/Utilities/GridAccessHelpers.hpp"
 #include "Acts/Utilities/Logger.hpp"
 
+#include <map>
 #include <memory>
 #include <string>
 #include <vector>
@@ -67,11 +69,12 @@ struct SensitiveCandidates {
 /// elements of the Acts::TrackingGeoemtry w/o map lookup.
 class SensitiveSurfaceMapper {
  public:
+  /// This prefix is used to indicate a sensitive volume that is matched
+  constexpr static std::string_view mappingPrefix = "ActsSensitive#";
+
   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
   struct Config {
     /// For which G4 material names we try to find a mapping
@@ -84,6 +87,15 @@ class SensitiveSurfaceMapper {
     SensitiveCandidates candidateSurfaces;
   };
 
+  /// State object that coutns the assignments and makes
+  /// a replica save copy association map
+  struct State {
+    /// The map of G4 physical volumes to the mapped surfaces (can be many as
+    /// there can be replicas)
+    std::multimap<const G4VPhysicalVolume*, const Acts::Surface*>
+        g4VolumeToSurfaces;
+  };
+
   /// Constructor with:
   ///
   /// @param cfg the configuration struct
@@ -98,14 +110,13 @@ class SensitiveSurfaceMapper {
   /// hierarchy and applies name remapping to the Physical volumes
   /// of the Geant4 geometry.
   ///
+  /// @param state is the state object (for caching)
   /// @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(const Acts::GeometryContext& gctx,
+  void remapSensitiveNames(State& state, const Acts::GeometryContext& gctx,
                            G4VPhysicalVolume* g4PhysicalVolume,
-                           const Acts::Transform3& motherTransform,
-                           int& sCounter) const;
+                           const Acts::Transform3& motherTransform) const;
 
  protected:
   /// Configuration object

Examples/Algorithms/Geant4/src/Geant4Simulation.cpp

@@ -215,6 +215,7 @@ ActsExamples::Geant4Simulation::Geant4Simulation(const Config& cfg,
   }
 
   // Stepping actions
+  SensitiveSteppingAction* sensitiveSteppingActionAccess = nullptr;
   {
     // Clear stepping action if it exists
     if (runManager().GetUserSteppingAction() != nullptr) {
@@ -237,8 +238,12 @@ ActsExamples::Geant4Simulation::Geant4Simulation(const Config& cfg,
     SteppingActionList::Config steppingCfg;
     steppingCfg.actions.push_back(std::make_unique<ParticleKillAction>(
         particleKillCfg, m_logger->cloneWithSuffix("Killer")));
-    steppingCfg.actions.push_back(std::make_unique<SensitiveSteppingAction>(
-        stepCfg, m_logger->cloneWithSuffix("SensitiveStepping")));
+
+    auto sensitiveSteppingAction = std::make_unique<SensitiveSteppingAction>(
+        stepCfg, m_logger->cloneWithSuffix("SensitiveStepping"));
+    sensitiveSteppingActionAccess = sensitiveSteppingAction.get();
+
+    steppingCfg.actions.push_back(std::move(sensitiveSteppingAction));
 
     // G4RunManager will take care of deletion
     auto steppingAction = new SteppingActionList(steppingCfg);
@@ -271,14 +276,16 @@ ActsExamples::Geant4Simulation::Geant4Simulation(const Config& cfg,
 
   // ACTS sensitive surfaces are provided, so hit creation is turned on
   if (cfg.sensitiveSurfaceMapper != nullptr) {
+    SensitiveSurfaceMapper::State sState;
     ACTS_INFO(
         "Remapping selected volumes from Geant4 to Acts::Surface::GeometryID");
-    int sCounter = 0;
     cfg.sensitiveSurfaceMapper->remapSensitiveNames(
-        Acts::GeometryContext{}, g4World, Acts::Transform3::Identity(),
-        sCounter);
+        sState, Acts::GeometryContext{}, g4World, Acts::Transform3::Identity());
+    ACTS_INFO("Remapping successful for " << sState.g4VolumeToSurfaces.size()
+                                          << " selected volumes.");
 
-    ACTS_INFO("Remapping successful for " << sCounter << " selected volumes.");
+    sensitiveSteppingActionAccess->assignSurfaceMapping(
+        sState.g4VolumeToSurfaces);
   }
 
   m_inputParticles.initialize(cfg.inputParticles);

Examples/Algorithms/Geant4/src/SensitiveSteppingAction.cpp

@@ -11,6 +11,7 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Definitions/Units.hpp"
 #include "Acts/Geometry/GeometryIdentifier.hpp"
+#include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/MultiIndex.hpp"
 #include "ActsExamples/EventData/SimHit.hpp"
 #include "ActsExamples/Geant4/EventStore.hpp"
@@ -28,6 +29,7 @@
 #include <G4Track.hh>
 #include <G4UnitsTable.hh>
 #include <G4VPhysicalVolume.hh>
+#include <G4VTouchable.hh>
 
 class G4PrimaryParticle;
 
@@ -97,6 +99,7 @@ ActsExamples::SensitiveSteppingAction::SensitiveSteppingAction(
 void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
     const G4Step* step) {
   // Unit conversions G4->::ACTS
+  static constexpr double convertLength = Acts::UnitConstants::mm / CLHEP::mm;
 
   // The particle after the step
   G4Track* track = step->GetTrack();
@@ -125,19 +128,53 @@ void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
   }
 
   // Get the physical volume & check if it has the sensitive string name
-  std::string_view volumeName = track->GetVolume()->GetName();
+  auto volume = track->GetVolume();
+  std::string volumeName = volume->GetName();
 
   if (volumeName.find(SensitiveSurfaceMapper::mappingPrefix) ==
       std::string_view::npos) {
     return;
   }
 
-  // Cast out the GeometryIdentifier
-  volumeName = volumeName.substr(SensitiveSurfaceMapper::mappingPrefix.size(),
-                                 std::string_view::npos);
-  char* end = nullptr;
-  const Acts::GeometryIdentifier geoId(
-      std::strtoul(volumeName.data(), &end, 10));
+  // Get PreStepPoint and PostStepPoint
+  const G4StepPoint* preStepPoint = step->GetPreStepPoint();
+  const G4StepPoint* postStepPoint = step->GetPostStepPoint();
+
+  // The G4Touchable for the matching
+  const G4VTouchable* touchable = track->GetTouchable();
+
+  Acts::GeometryIdentifier geoId{};
+
+  // Find the range of candidate surfaces for the current position in the
+  // mapping multimap
+  auto [bsf, esf] = m_surfaceMapping.equal_range(volume);
+  std::size_t nSurfaces = std::distance(bsf, esf);
+
+  ACTS_VERBOSE("Found " << nSurfaces << " candidate surfaces for volume "
+                        << volume->GetName());
+
+  if (nSurfaces == 0) {
+    ACTS_ERROR("No candidate surfaces found for volume " << volume->GetName());
+    return;
+  } else if (nSurfaces == 1u) {
+    geoId = bsf->second->geometryId();
+    ACTS_VERBOSE("Unique assignment successful -> to surface " << geoId);
+  } else {
+    // Find the closest surface to the current position
+    Acts::GeometryContext gctx;
+    for (; bsf != esf; ++bsf) {
+      auto surface = bsf->second;
+      auto translation = touchable->GetTranslation();
+      Acts::Vector3 g4VolumePosition(convertLength * translation.x(),
+                                     convertLength * translation.y(),
+                                     convertLength * translation.z());
+      if (surface->center(gctx).isApprox(g4VolumePosition)) {
+        geoId = surface->geometryId();
+        break;
+      }
+    }
+    ACTS_VERBOSE("Replica assignment successful -> to surface " << geoId);
+  }
 
   // This is not the case if we have a particle-ID collision
   if (eventStore().trackIdMapping.find(track->GetTrackID()) ==
@@ -149,10 +186,6 @@ void ActsExamples::SensitiveSteppingAction::UserSteppingAction(
 
   ACTS_VERBOSE("Step of " << particleId << " in sensitive volume " << geoId);
 
-  // Get PreStepPoint and PostStepPoint
-  const G4StepPoint* preStepPoint = step->GetPreStepPoint();
-  const G4StepPoint* postStepPoint = step->GetPostStepPoint();
-
   // Set particle hit count to zero, so we have this entry in the map later
   if (eventStore().particleHitCount.find(particleId) ==
       eventStore().particleHitCount.end()) {

Examples/Algorithms/Geant4/src/SensitiveSurfaceMapper.cpp

@@ -24,8 +24,10 @@ ActsExamples::SensitiveSurfaceMapper::SensitiveSurfaceMapper(
     : m_cfg(cfg), m_logger(std::move(logger)) {}
 
 void ActsExamples::SensitiveSurfaceMapper::remapSensitiveNames(
-    const Acts::GeometryContext& gctx, G4VPhysicalVolume* g4PhysicalVolume,
-    const Acts::Transform3& motherTransform, int& sCounter) const {
+    State& state, const Acts::GeometryContext& gctx,
+    G4VPhysicalVolume* g4PhysicalVolume,
+    const Acts::Transform3& motherTransform) const {
+  // Make sure the unit conversion is correct
   constexpr double convertLength = CLHEP::mm / Acts::UnitConstants::mm;
 
   auto g4LogicalVolume = g4PhysicalVolume->GetLogicalVolume();
@@ -53,8 +55,8 @@ void ActsExamples::SensitiveSurfaceMapper::remapSensitiveNames(
   if (G4int nDaughters = g4LogicalVolume->GetNoDaughters(); nDaughters > 0) {
     // Step down to all daughters
     for (G4int id = 0; id < nDaughters; ++id) {
-      remapSensitiveNames(gctx, g4LogicalVolume->GetDaughter(id), transform,
-                          sCounter);
+      remapSensitiveNames(state, gctx, g4LogicalVolume->GetDaughter(id),
+                          transform);
     }
     return;
   }
@@ -86,15 +88,18 @@ void ActsExamples::SensitiveSurfaceMapper::remapSensitiveNames(
 
     // A mapped surface was found, a new name will be set that G4PhysVolume
     if (mappedSurface != nullptr) {
-      ++sCounter;
-      std::string mappedVolumeName(SensitiveSurfaceMapper::mappingPrefix);
-      mappedVolumeName += std::to_string(mappedSurface->geometryId().value());
       ACTS_VERBOSE("Found matching surface " << mappedSurface->geometryId()
                                              << " at position "
                                              << g4RelPosition.transpose());
-      ACTS_VERBOSE("Remap: " << g4PhysicalVolume->GetName() << " -> "
-                             << mappedVolumeName);
-      g4PhysicalVolume->SetName(mappedVolumeName.c_str());
+      // Check if the prefix is not yet assigned
+      std::string prefixName = std::string(mappingPrefix);
+      if (volumeName.find(prefixName) == std::string::npos) {
+        // Set the new name
+        std::string mappedName = prefixName + volumeName;
+        g4PhysicalVolume->SetName(mappedName);
+      }
+      // Insert into the multi-map
+      state.g4VolumeToSurfaces.insert({g4PhysicalVolume, mappedSurface});
     } else {
       ACTS_VERBOSE("No mapping found for '"
                    << volumeName << "' with material '" << volumeMaterialName