From 13b4fcb60340f0d2714203683a489750830cd1e2 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Thu, 6 Aug 2026 17:05:03 +0200 Subject: [PATCH] fix: Repair the Geant4 telescope detector geometry The Geant4 geometry did not match the tracking geometry: - Of the two identically sized envelopes, the inner one was placed at the center of the layer stack and therefore protruded out of its mother, leaving Geant4 navigation undefined over most of the telescope. Both are replaced by a single envelope placed directly in the world. - The silicon planes treated `bounds` as full lengths, while `buildTelescopeDetector` uses them as half lengths, so the Geant4 sensors were half the size of the ACTS surfaces in each direction. - The world size assumed sorted and positive positions. Also turns a segfault on a `TelescopeDetector` without positions into an exception, and makes the Geant4 sources `PRIVATE` so that consumers of `ActsExamplesDetectorTelescope` stop recompiling them. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01KUNabCsitrfbbWnxvZzgJc --- .../TelescopeDetector/CMakeLists.txt | 4 +- .../src/TelescopeDetector.cpp | 8 + .../src/TelescopeG4DetectorConstruction.cpp | 153 +++++++------- Python/Examples/tests/test_detectors.py | 35 +++- Tests/UnitTests/Examples/CMakeLists.txt | 1 + .../Examples/Detectors/CMakeLists.txt | 1 + .../TelescopeDetector/CMakeLists.txt | 6 + .../TelescopeG4DetectorConstructionTests.cpp | 187 ++++++++++++++++++ 8 files changed, 318 insertions(+), 77 deletions(-) create mode 100644 Tests/UnitTests/Examples/Detectors/CMakeLists.txt create mode 100644 Tests/UnitTests/Examples/Detectors/TelescopeDetector/CMakeLists.txt create mode 100644 Tests/UnitTests/Examples/Detectors/TelescopeDetector/TelescopeG4DetectorConstructionTests.cpp diff --git a/Examples/Detectors/TelescopeDetector/CMakeLists.txt b/Examples/Detectors/TelescopeDetector/CMakeLists.txt index 4216ee5a598..521498b89e0 100644 --- a/Examples/Detectors/TelescopeDetector/CMakeLists.txt +++ b/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() diff --git a/Examples/Detectors/TelescopeDetector/src/TelescopeDetector.cpp b/Examples/Detectors/TelescopeDetector/src/TelescopeDetector.cpp index c954f0068cc..3670f9ba34d 100644 --- a/Examples/Detectors/TelescopeDetector/src/TelescopeDetector.cpp +++ b/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 " diff --git a/Examples/Detectors/TelescopeDetector/src/TelescopeG4DetectorConstruction.cpp b/Examples/Detectors/TelescopeDetector/src/TelescopeG4DetectorConstruction.cpp index 157204be9df..bac2e73d5e8 100644 --- a/Examples/Detectors/TelescopeDetector/src/TelescopeG4DetectorConstruction.cpp +++ b/Examples/Detectors/TelescopeDetector/src/TelescopeG4DetectorConstruction.cpp @@ -13,8 +13,12 @@ #include "ActsExamples/TelescopeDetector/BuildTelescopeDetector.hpp" #include "ActsExamples/TelescopeDetector/TelescopeDetector.hpp" +#include +#include + #include "G4Box.hh" #include "G4LogicalVolume.hh" +#include "G4Material.hh" #include "G4NistManager.hh" #include "G4PVPlacement.hh" #include "G4RunManager.hh" @@ -27,7 +31,7 @@ TelescopeG4DetectorConstruction::TelescopeG4DetectorConstruction( const Geant4ConstructionOptions& options) : m_cfg(cfg), m_options(options) { throw_assert(cfg.surfaceType == static_cast(TelescopeSurfaceType::Plane), - "only plan is supported right now"); + "only plane is supported right now"); } G4VPhysicalVolume* TelescopeG4DetectorConstruction::Construct() { @@ -35,38 +39,43 @@ G4VPhysicalVolume* TelescopeG4DetectorConstruction::Construct() { 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(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,50 +116,31 @@ 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 @@ -143,14 +148,14 @@ G4VPhysicalVolume* TelescopeG4DetectorConstruction::Construct() { 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(i), // copy number + checkOverlaps); // overlaps checking } // Create regions diff --git a/Python/Examples/tests/test_detectors.py b/Python/Examples/tests/test_detectors.py index 68149e37ea4..8b382214a62 100644 --- a/Python/Examples/tests/test_detectors.py +++ b/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 diff --git a/Tests/UnitTests/Examples/CMakeLists.txt b/Tests/UnitTests/Examples/CMakeLists.txt index 2a4b693c239..2230f0aa4ba 100644 --- a/Tests/UnitTests/Examples/CMakeLists.txt +++ b/Tests/UnitTests/Examples/CMakeLists.txt @@ -1,4 +1,5 @@ add_subdirectory(Algorithms) +add_subdirectory(Detectors) add_subdirectory(EventData) add_subdirectory(Io) add_subdirectory(Framework) diff --git a/Tests/UnitTests/Examples/Detectors/CMakeLists.txt b/Tests/UnitTests/Examples/Detectors/CMakeLists.txt new file mode 100644 index 00000000000..09d065ed65a --- /dev/null +++ b/Tests/UnitTests/Examples/Detectors/CMakeLists.txt @@ -0,0 +1 @@ +add_subdirectory_if(TelescopeDetector ACTS_BUILD_EXAMPLES_GEANT4) diff --git a/Tests/UnitTests/Examples/Detectors/TelescopeDetector/CMakeLists.txt b/Tests/UnitTests/Examples/Detectors/TelescopeDetector/CMakeLists.txt new file mode 100644 index 00000000000..e0144ad600d --- /dev/null +++ b/Tests/UnitTests/Examples/Detectors/TelescopeDetector/CMakeLists.txt @@ -0,0 +1,6 @@ +set(unittest_extra_libraries ActsExamplesDetectorTelescope ActsExamplesGeant4) + +add_unittest( + TelescopeG4DetectorConstruction + TelescopeG4DetectorConstructionTests.cpp +) diff --git a/Tests/UnitTests/Examples/Detectors/TelescopeDetector/TelescopeG4DetectorConstructionTests.cpp b/Tests/UnitTests/Examples/Detectors/TelescopeDetector/TelescopeG4DetectorConstructionTests.cpp new file mode 100644 index 00000000000..6a9bcf572d6 --- /dev/null +++ b/Tests/UnitTests/Examples/Detectors/TelescopeDetector/TelescopeG4DetectorConstructionTests.cpp @@ -0,0 +1,187 @@ +// This file is part of the ACTS project. +// +// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/. + +#include +#include + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Geometry/GeometryContext.hpp" +#include "Acts/Geometry/TrackingGeometry.hpp" +#include "Acts/Surfaces/RectangleBounds.hpp" +#include "Acts/Surfaces/Surface.hpp" +#include "ActsExamples/Geant4/AlgebraConverters.hpp" +#include "ActsExamples/Geant4/Geant4ConstructionOptions.hpp" +#include "ActsExamples/TelescopeDetector/TelescopeDetector.hpp" +#include "ActsExamples/TelescopeDetector/TelescopeG4DetectorConstruction.hpp" + +#include +#include +#include + +#include +#include +#include +#include +#include + +using namespace ActsExamples; + +namespace { + +/// A sensitive Geant4 volume in global coordinates +struct G4Sensitive { + Acts::Vector3 center; + Acts::Vector2 halfLengths; +}; + +/// Collect all silicon volumes of the tree, composing transforms the way +/// `SensitiveSurfaceMapper` does +void collectSensitives(const G4VPhysicalVolume& physicalVolume, + const Acts::Transform3& motherTransform, + std::vector& sensitives) { + const G4LogicalVolume& logicalVolume = *physicalVolume.GetLogicalVolume(); + + Acts::Transform3 localToGlobal = + motherTransform * Acts::Translation3(Geant4::convertPosition( + physicalVolume.GetTranslation())); + if (const G4RotationMatrix* g4Rotation = physicalVolume.GetRotation(); + g4Rotation != nullptr) { + // Geant4 stores the rotation of the mother relative to the daughter + // frame, hence the transpose + Acts::RotationMatrix3 rotation; + rotation << g4Rotation->xx(), g4Rotation->yx(), g4Rotation->zx(), + g4Rotation->xy(), g4Rotation->yy(), g4Rotation->zy(), g4Rotation->xz(), + g4Rotation->yz(), g4Rotation->zz(); + localToGlobal.rotate(rotation); + } + + if (logicalVolume.GetMaterial()->GetName() == "Silicon") { + const auto& box = dynamic_cast(*logicalVolume.GetSolid()); + constexpr double convertLength = CLHEP::mm / Acts::UnitConstants::mm; + sensitives.push_back({localToGlobal.translation(), + {box.GetXHalfLength() * convertLength, + box.GetYHalfLength() * convertLength}}); + } + + for (std::size_t i = 0; i < logicalVolume.GetNoDaughters(); ++i) { + collectSensitives(*logicalVolume.GetDaughter(i), localToGlobal, sensitives); + } +} + +/// Check the whole tree for overlaps +bool hasOverlaps(G4VPhysicalVolume& physicalVolume) { + constexpr G4int resolution = 1000; + constexpr G4double tolerance = 0.; + constexpr G4bool verbose = false; + bool overlaps = physicalVolume.CheckOverlaps(resolution, tolerance, verbose); + + G4LogicalVolume& logicalVolume = *physicalVolume.GetLogicalVolume(); + for (std::size_t i = 0; i < logicalVolume.GetNoDaughters(); ++i) { + overlaps |= hasOverlaps(*logicalVolume.GetDaughter(i)); + } + return overlaps; +} + +TelescopeDetector::Config makeConfig(int binValue) { + TelescopeDetector::Config cfg; + cfg.positions = {30, 60, 90, 120, 150, 180}; + cfg.stereos = {0, 0, 0, 0, 0, 0}; + cfg.offsets = {10, -20}; + cfg.bounds = {25, 100}; + cfg.binValue = binValue; + return cfg; +} + +} // namespace + +BOOST_AUTO_TEST_SUITE(TelescopeG4DetectorConstructionTests) + +BOOST_AUTO_TEST_CASE(ConstructSmoke) { + TelescopeDetector detector{TelescopeDetector::Config{}}; + + auto construction = + detector.buildGeant4DetectorConstruction(Geant4ConstructionOptions{}); + BOOST_REQUIRE(construction != nullptr); + + G4VPhysicalVolume* world = construction->Construct(); + BOOST_REQUIRE(world != nullptr); + // The world is built once and cached + BOOST_CHECK_EQUAL(construction->Construct(), world); + + BOOST_CHECK(!hasOverlaps(*world)); +} + +// The Geant4 sensors have to match the tracking geometry surfaces on all axes +BOOST_DATA_TEST_CASE(MatchesTrackingGeometry, + boost::unit_test::data::make({0, 1, 2}), binValue) { + const auto cfg = makeConfig(binValue); + TelescopeDetector detector{cfg}; + + auto construction = + detector.buildGeant4DetectorConstruction(Geant4ConstructionOptions{}); + G4VPhysicalVolume* world = construction->Construct(); + BOOST_REQUIRE(world != nullptr); + + BOOST_CHECK(!hasOverlaps(*world)); + + std::vector sensitives; + collectSensitives(*world, Acts::Transform3::Identity(), sensitives); + BOOST_REQUIRE_EQUAL(sensitives.size(), cfg.positions.size()); + + const auto& gctx = detector.nominalGeometryContext(); + std::vector surfaces; + detector.trackingGeometry()->visitSurfaces( + [&](const Acts::Surface* surface) { surfaces.push_back(surface); }, true); + BOOST_REQUIRE_EQUAL(surfaces.size(), cfg.positions.size()); + + for (const auto* surface : surfaces) { + const Acts::Vector3 center = surface->center(gctx); + + auto it = std::ranges::find_if(sensitives, [&](const G4Sensitive& g4) { + return (g4.center - center).norm() < 1e-6; + }); + BOOST_REQUIRE_MESSAGE( + it != sensitives.end(), + "no Geant4 volume at surface center " << center.transpose()); + + const auto& bounds = + dynamic_cast(surface->bounds()); + BOOST_CHECK_CLOSE(it->halfLengths[0], bounds.halfLengthX(), 1e-6); + BOOST_CHECK_CLOSE(it->halfLengths[1], bounds.halfLengthY(), 1e-6); + } +} + +// Unsorted positions and a stack that is not centered on the origin +BOOST_AUTO_TEST_CASE(UnsortedAndShiftedPositions) { + TelescopeDetector::Config cfg; + cfg.positions = {-100, 200, 50}; + cfg.stereos = {0, 0, 0}; + cfg.offsets = {40, 40}; + + TelescopeDetector detector{cfg}; + auto construction = + detector.buildGeant4DetectorConstruction(Geant4ConstructionOptions{}); + G4VPhysicalVolume* world = construction->Construct(); + BOOST_REQUIRE(world != nullptr); + + BOOST_CHECK(!hasOverlaps(*world)); + + std::vector sensitives; + collectSensitives(*world, Acts::Transform3::Identity(), sensitives); + BOOST_CHECK_EQUAL(sensitives.size(), cfg.positions.size()); +} + +BOOST_AUTO_TEST_CASE(EmptyPositionsThrows) { + TelescopeDetector::Config cfg; + cfg.positions = {}; + cfg.stereos = {}; + + BOOST_CHECK_THROW(TelescopeDetector{cfg}, std::invalid_argument); +} + +BOOST_AUTO_TEST_SUITE_END()