feat: Robust surface matching for Geant4 (acts-project#3169)

This PR fixes:

acts-project#3165

It first creates a `std::multimap`  that relates the `G4VPhysVolume` to the possible `Acts::Surface` objects (there can be more due to the replica mechanism).

In the Stepping action, the right one will then be picked.

The output performance change in geant_hits is expected.

Digitisation configuration file:
[odd-digi-geometric-config.json](https://github.com/acts-project/acts/files/15188041/odd-digi-geometric-config.json)

Example how to run full simulation with this:
```sh
python <path_to_your_acts_source_dir>Examples/Scripts/Python/sim_digi_odd.py --digi-config odd-digi-geometric-config.json --geant4 --gun-multiplicity 10 --gun-pt-range 100 100 -n100 -o g4_mu_100GeV_m10
```

Options are:
```sh
usage: sim_digi_odd.py [-h] [--output OUTPUT] [--events EVENTS] [--skip SKIP] [--edm4hep EDM4HEP] [--geant4] [--ttbar] [--ttbar-pu TTBAR_PU] [--gun-multiplicity GUN_MULTIPLICITY] [--gun-eta-range GUN_ETA_RANGE [GUN_ETA_RANGE ...]] [--gun-pt-range GUN_PT_RANGE [GUN_PT_RANGE ...]] [--rnd-seed RND_SEED]
                       [--digi-config DIGI_CONFIG]

Sim-digi chain with the OpenDataDetector

options:
  -h, --help            show this help message and exit
  --output OUTPUT, -o OUTPUT
                        Output directory
  --events EVENTS, -n EVENTS
                        Number of events
  --skip SKIP, -s SKIP  Number of events
  --edm4hep EDM4HEP     Use edm4hep inputs
  --geant4              Use Geant4 instead of fatras
  --ttbar               Use Pythia8 (ttbar, pile-up 200) instead of particle gun
  --ttbar-pu TTBAR_PU   Number of pile-up events for ttbar
  --gun-multiplicity GUN_MULTIPLICITY
                        Multiplicity of the particle gun
  --gun-eta-range GUN_ETA_RANGE [GUN_ETA_RANGE ...]
                        Eta range of the particle gun
  --gun-pt-range GUN_PT_RANGE [GUN_PT_RANGE ...]
                        Pt range of the particle gun (GeV)
  --rnd-seed RND_SEED   Random seed
  --digi-config DIGI_CONFIG
                        Digitization configuration file
```

Author: asalzburger

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 = 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();
+  const G4VPhysicalVolume* 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 "
+                        << volumeName);
+
+  if (nSurfaces == 0) {
+    ACTS_ERROR("No candidate surfaces found for volume " << volumeName);
+    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) {
+      const Acts::Surface* surface = bsf->second;
+      const G4ThreeVector& 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,17 @@ 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
+      if (volumeName.find(mappingPrefix) == std::string::npos) {
+        // Set the new name
+        std::string mappedName = std::string(mappingPrefix) + 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

Examples/Python/tests/root_file_hashes.txt

@@ -2,7 +2,7 @@ test_pythia8__pythia8_particles.root: 49b89c3458a51aa9407f887be50e6bbcd9a2a0c897
 test_fatras__particles_simulation.root: bc970873fef0c2efd86ed5413623802353d2cd04abea72de14e8cdfc0e40076f
 test_fatras__hits.root: 821d9f3ce364bfea6c08df86b0204d7257cc1810f5b0bdfda45950e697ef5c52
 test_geant4__particles_simulation.root: a51f15555abefec7e9d9de31f9d856a3026af1b991db73fefaec5ae160a4af2f
-test_geant4__hits.root: 794fbae710b8ec68089efb53c5722025bd2cfd961f2f1dc7102c55e65e75dc5f
+test_geant4__hits.root: adf5dcdf000a580412dc5089e17460897d6535c978eafa021584ba4281d0a1ac
 test_seeding__estimatedparams.root: 1fa2e879142059344e32ab62e3763dce14b5138b591c6e919073cb985782075f
 test_seeding__performance_seeding.root: 992f9c611d30dde0d3f3ab676bab19ada61ab6a4442828e27b65ec5e5b7a2880
 test_seeding__particles.root: 7855b021f39ad238bca098e4282667be0666f2d1630e5bcb9d51d3b5ee39fa14