fix: Repair the Geant4 telescope detector geometry

Examples/Detectors/TelescopeDetector/CMakeLists.txt

@@ -23,13 +23,13 @@ if(ACTS_BUILD_EXAMPLES_GEANT4)
 
     target_sources(
         ActsExamplesDetectorTelescope
-        PUBLIC
+        PRIVATE
             src/TelescopeDetectorGeant4.cpp
             src/TelescopeG4DetectorConstruction.cpp
     )
 else()
     target_sources(
         ActsExamplesDetectorTelescope
-        PUBLIC src/TelescopeDetectorGeant4Stub.cpp
+        PRIVATE src/TelescopeDetectorGeant4Stub.cpp
     )
 endif()

Examples/Detectors/TelescopeDetector/src/TelescopeDetector.cpp

@@ -32,6 +32,10 @@ TelescopeDetector::TelescopeDetector(const Config& cfg)
         "The minR should be smaller than the maxR for disc surface bounds.");
   }
 
+  if (m_cfg.positions.empty()) {
+    throw std::invalid_argument("At least one surface position is required.");
+  }
+
   if (m_cfg.positions.size() != m_cfg.stereos.size()) {
     throw std::invalid_argument(
         "The number of provided positions must match the number of "
@@ -64,6 +68,10 @@ TelescopeDetector::TelescopeDetector(const Config& cfg, NoBuildTag /*unused*/)
         "The minR should be smaller than the maxR for disc surface bounds.");
   }
 
+  if (m_cfg.positions.empty()) {
+    throw std::invalid_argument("At least one surface position is required.");
+  }
+
   if (m_cfg.positions.size() != m_cfg.stereos.size()) {
     throw std::invalid_argument(
         "The number of provided positions must match the number of "

Examples/Detectors/TelescopeDetector/src/TelescopeG4DetectorConstruction.cpp

@@ -13,8 +13,12 @@
 #include "ActsExamples/TelescopeDetector/BuildTelescopeDetector.hpp"
 #include "ActsExamples/TelescopeDetector/TelescopeDetector.hpp"
 
+#include <algorithm>
+#include <cmath>
+
 #include "G4Box.hh"
 #include "G4LogicalVolume.hh"
+#include "G4Material.hh"
 #include "G4NistManager.hh"
 #include "G4PVPlacement.hh"
 #include "G4RunManager.hh"
@@ -27,46 +31,51 @@ TelescopeG4DetectorConstruction::TelescopeG4DetectorConstruction(
     const Geant4ConstructionOptions& options)
     : m_cfg(cfg), m_options(options) {
   throw_assert(cfg.surfaceType == static_cast<int>(TelescopeSurfaceType::Plane),
-               "only plan is supported right now");
+               "only plane is supported right now");
 }
 
 G4VPhysicalVolume* TelescopeG4DetectorConstruction::Construct() {
   if (m_world != nullptr) {
     return m_world;
   }
 
-  G4double center =
-      (m_cfg.positions.back() + m_cfg.positions.front()) * 0.5 * mm;
-  G4double length = (m_cfg.positions.back() - m_cfg.positions.front()) * mm;
+  // Option to switch on/off checking of volumes overlaps
+  constexpr G4bool checkOverlaps = true;
 
-  // Get nist material manager
-  G4NistManager* nist = G4NistManager::Instance();
+  // Keeps nested volumes from sharing surfaces
+  constexpr G4double margin = 1 * mm;
 
-  // World
-  //
-  G4double worldSize =
-      std::max({std::abs(m_cfg.offsets[0]) + m_cfg.bounds[0] * 0.5,
-                std::abs(m_cfg.offsets[1]) + m_cfg.bounds[1] * 0.5,
-                m_cfg.positions.back() + m_cfg.thickness});
+  // The positions need not be sorted
+  const auto [minPosition, maxPosition] = std::ranges::minmax(m_cfg.positions);
 
-  // Envelope parameters
-  //
-  G4double envSizeX = m_cfg.bounds[0] * mm;
-  G4double envSizeY = m_cfg.bounds[1] * mm;
-  G4double envSizeZ = length + m_cfg.thickness * mm;
+  // Extent of the layer stack along the telescope axis
+  const G4double stackCenter = (minPosition + maxPosition) * 0.5 * mm;
+  const G4double stackLength = (maxPosition - minPosition) * mm;
 
-  // Option to switch on/off checking of volumes overlaps
-  //
-  G4bool checkOverlaps = true;
+  // `bounds` are half lengths, matching the `Acts::RectangleBounds` that
+  // `buildTelescopeDetector` gives the sensitive surfaces
+  const G4double layerHalfX = m_cfg.bounds[0] * mm;
+  const G4double layerHalfY = m_cfg.bounds[1] * mm;
+  const G4double layerHalfZ = m_cfg.thickness * 0.5 * mm;
 
-  // Materials
+  const G4double envHalfX = layerHalfX + margin;
+  const G4double envHalfY = layerHalfY + margin;
+  const G4double envHalfZ = stackLength * 0.5 + layerHalfZ + margin;
+
+  // Get nist material manager
+  G4NistManager* nist = G4NistManager::Instance();
+
+  // `SensitiveSurfaceMapper` keys on the silicon material name
   G4Material* galactic = nist->FindOrBuildMaterial("G4_Galactic");
-  G4Material* silicon =
-      new G4Material("Silicon", 14, 28.0855 * g / mole, 2.329 * g / cm3);
+  G4Material* silicon = G4Material::GetMaterial("Silicon", false);
+  if (silicon == nullptr) {
+    silicon =
+        new G4Material("Silicon", 14, 28.0855 * g / mole, 2.329 * g / cm3);
+  }
 
-  // Construct the rotation
-  // This assumes the binValue is AxisDirection::AxisX, AxisDirection::AxisY or
-  // AxisDirection::AxisZ. No reset is necessary in case of AxisDirection::AxisZ
+  // Orientation of the envelope in the world, assuming binValue is AxisX,
+  // AxisY or AxisZ. `G4PVPlacement` takes the rotation of the mother relative
+  // to the daughter frame, i.e. the inverse of that orientation.
   G4RotationMatrix* rotation = nullptr;
   if (static_cast<Acts::AxisDirection>(m_cfg.binValue) ==
       Acts::AxisDirection::AxisX) {
@@ -76,9 +85,24 @@ G4VPhysicalVolume* TelescopeG4DetectorConstruction::Construct() {
     rotation = new G4RotationMatrix({1, 0, 0}, {0, 0, 1}, {0, -1, 0});
   }
 
+  // The envelope center in the world frame
+  G4ThreeVector envCenter(m_cfg.offsets[0] * mm, m_cfg.offsets[1] * mm,
+                          stackCenter);
+  if (rotation != nullptr) {
+    envCenter = rotation->inverse() * envCenter;
+  }
+
+  // The world has to contain the envelope in any orientation
+  const G4double worldHalfSize =
+      envCenter.mag() +
+      std::sqrt(envHalfX * envHalfX + envHalfY * envHalfY +
+                envHalfZ * envHalfZ) +
+      margin;
+
   // World
   //
-  G4Box* solidWorld = new G4Box("World Solid", worldSize, worldSize, worldSize);
+  G4Box* solidWorld =
+      new G4Box("World Solid", worldHalfSize, worldHalfSize, worldHalfSize);
 
   G4LogicalVolume* logicWorld =
       new G4LogicalVolume(solidWorld, galactic, "World Logic");
@@ -92,65 +116,46 @@ G4VPhysicalVolume* TelescopeG4DetectorConstruction::Construct() {
                               0,                // copy number
                               checkOverlaps);   // overlaps checking
 
-  // Envelope 1
-  //
-  G4Box* solidEnv =
-      new G4Box("Envelope Solid",                                 // its name
-                0.5 * envSizeX, 0.5 * envSizeY, 0.5 * envSizeZ);  // its size
-
-  G4LogicalVolume* logicEnv1 =
-      new G4LogicalVolume(solidEnv,              // its solid
-                          galactic,              // its material
-                          "Envelope #1 Logic");  // its name
-
-  G4VPhysicalVolume* physEnv1 =
-      new G4PVPlacement(rotation,            // rotation
-                        G4ThreeVector(),     // at detector center
-                        logicEnv1,           // its logical volume
-                        "Envelope #1 Phys",  // its name
-                        logicWorld,          // its mother volume
-                        false,               // no boolean operation
-                        0,                   // copy number
-                        checkOverlaps);      // overlaps checking
-
-  // Envelope 2
+  // Envelope
   //
-  G4LogicalVolume* logicEnv2 =
-      new G4LogicalVolume(solidEnv,              // its solid
-                          galactic,              // its material
-                          "Envelope #2 Logic");  // its name
-
-  G4VPhysicalVolume* physEnv2 = new G4PVPlacement(
-      nullptr,  // no rotation
-      G4ThreeVector(m_cfg.offsets[0] * mm, m_cfg.offsets[1] * mm,
-                    center),  // at detector center
-      "Envelope #2 Phys",     // its name
-      logicEnv2,              // its logical volume
-      physEnv1,               // its mother volume
-      false,                  // no boolean operation
-      0,                      // copy number
-      checkOverlaps);         // overlaps checking
+  G4Box* solidEnv = new G4Box("Envelope Solid",               // its name
+                              envHalfX, envHalfY, envHalfZ);  // its size
+
+  G4LogicalVolume* logicEnv =
+      new G4LogicalVolume(solidEnv,           // its solid
+                          galactic,           // its material
+                          "Envelope Logic");  // its name
+
+  G4VPhysicalVolume* physEnv =
+      new G4PVPlacement(rotation,         // rotation
+                        envCenter,        // at the center of the stack
+                        logicEnv,         // its logical volume
+                        "Envelope Phys",  // its name
+                        logicWorld,       // its mother volume
+                        false,            // no boolean operation
+                        0,                // copy number
+                        checkOverlaps);   // overlaps checking
 
   // Layer
   //
 
-  G4Box* solidLayer = new G4Box("Layer Solid", 0.5 * m_cfg.bounds[0],
-                                0.5 * m_cfg.bounds[1], 0.5 * m_cfg.thickness);
+  G4Box* solidLayer =
+      new G4Box("Layer Solid", layerHalfX, layerHalfY, layerHalfZ);
 
   G4LogicalVolume* logicLayer = new G4LogicalVolume(solidLayer,  // its solid
                                                     silicon,     // its material
                                                     "Layer Logic");  // its name
 
   for (std::size_t i = 0; i < m_cfg.positions.size(); ++i) {
     new G4PVPlacement(
-        nullptr,                                                // no rotation
-        G4ThreeVector(0, 0, m_cfg.positions[i] * mm - center),  // at position
-        "Layer #" + std::to_string(i) + " Phys",                // its name
-        logicLayer,      // its logical volume
-        physEnv2,        // its mother volume
-        false,           // no boolean operation
-        0,               // copy number
-        checkOverlaps);  // overlaps checking
+        nullptr,  // no rotation
+        G4ThreeVector(0, 0, m_cfg.positions[i] * mm - stackCenter),  // position
+        "Layer #" + std::to_string(i) + " Phys",                     // its name
+        logicLayer,             // its logical volume
+        physEnv,                // its mother volume
+        false,                  // no boolean operation
+        static_cast<G4int>(i),  // copy number
+        checkOverlaps);         // overlaps checking
   }
 
   // Create regions

Python/Examples/tests/test_detectors.py

@@ -1,7 +1,7 @@
 import pytest
 from pathlib import Path
 
-from helpers import dd4hepEnabled
+from helpers import dd4hepEnabled, geant4Enabled
 
 import acts.examples
 
@@ -56,6 +56,39 @@ def test_telescope_geometry():
     assert count_surfaces(trackingGeometry) == n_surfaces
 
 
+@pytest.mark.skipif(not geant4Enabled, reason="Geant4 is not set up")
+@pytest.mark.parametrize("binValue", [0, 1, 2])
+def test_telescope_geant4_geometry(binValue):
+    from acts.examples.geant4 import SensitiveSurfaceMapper
+
+    n_surfaces = 6
+
+    config = acts.examples.TelescopeDetector.Config(
+        bounds=[100, 200],
+        positions=[30 * (i + 1) for i in range(n_surfaces)],
+        stereos=[0] * n_surfaces,
+        offsets=[10, -20],
+        binValue=binValue,
+    )
+    detector = acts.examples.TelescopeDetector(config)
+    trackingGeometry = detector.trackingGeometry()
+    gctx = detector.nominalGeometryContext()
+
+    # every sensitive surface has to be backed by a Geant4 volume in the same place
+    smmConfig = SensitiveSurfaceMapper.Config()
+    smmConfig.materialMappings = ["Silicon"]
+    mapper = SensitiveSurfaceMapper.create(
+        smmConfig, acts.logging.INFO, trackingGeometry
+    )
+
+    state = SensitiveSurfaceMapper.State()
+    mapper.remapSensitiveNames(
+        state, gctx, detector, acts.Transform3(acts.Vector3(0, 0, 0))
+    )
+
+    assert mapper.checkMapping(state, gctx, False, False)
+
+
 @pytest.mark.skipif(not dd4hepEnabled, reason="DD4hep is not set up")
 def test_odd(odd_detector):
     detector = odd_detector

Tests/UnitTests/Examples/CMakeLists.txt

@@ -1,4 +1,5 @@
 add_subdirectory(Algorithms)
+add_subdirectory(Detectors)
 add_subdirectory(EventData)
 add_subdirectory(Io)
 add_subdirectory(Framework)

Tests/UnitTests/Examples/Detectors/CMakeLists.txt

@@ -0,0 +1 @@
+add_subdirectory_if(TelescopeDetector ACTS_BUILD_EXAMPLES_GEANT4)

Tests/UnitTests/Examples/Detectors/TelescopeDetector/CMakeLists.txt

@@ -0,0 +1,6 @@
+set(unittest_extra_libraries ActsExamplesDetectorTelescope ActsExamplesGeant4)
+
+add_unittest(
+    TelescopeG4DetectorConstruction
+    TelescopeG4DetectorConstructionTests.cpp
+)