Algorithms/Geant4/SensitiveSurfaceMapperTests.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/data/test_case.hpp>
0010 #include <boost/test/unit_test.hpp>
0011
0012 #include "Acts/Definitions/Algebra.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/Geometry/GeometryContext.hpp"
0015 #include "Acts/Surfaces/AnnulusBounds.hpp"
0016 #include "Acts/Surfaces/CylinderBounds.hpp"
0017 #include "Acts/Surfaces/CylinderSurface.hpp"
0018 #include "Acts/Surfaces/DiscSurface.hpp"
0019 #include "Acts/Surfaces/PlaneSurface.hpp"
0020 #include "Acts/Surfaces/RectangleBounds.hpp"
0021 #include "Acts/Surfaces/Surface.hpp"
0022 #include "Acts/Utilities/Logger.hpp"
0023 #include "ActsExamples/Geant4/SensitiveSurfaceMapper.hpp"
0024
0025 #include <algorithm>
0026 #include <cmath>
0027 #include <limits>
0028 #include <memory>
0029 #include <numbers>
0030 #include <string>
0031 #include <utility>
0032 #include <vector>
0033
0034 #include <G4Box.hh>
0035 #include <G4LogicalVolume.hh>
0036 #include <G4NistManager.hh>
0037 #include <G4PVPlacement.hh>
0038 #include <G4Tubs.hh>
0039 #include <G4VPhysicalVolume.hh>
0040
0041 using namespace Acts::UnitLiterals;
0042 using namespace ActsExamples::Geant4;
0043
0044 namespace {
0045
0046 auto gctx = Acts::GeometryContext::dangerouslyDefaultConstruct();
0047
0048 /// Translation tolerance with which the mapper compares surface centers and
0049 /// Geant4 positions
0050 constexpr double kTolerance = 0.1_um;
0051
0052 /// All surfaces are candidates for every position. With `cacheByRegion` they
0053 /// form a single region, so the mapper builds its lookup once; without it the
0054 /// mapper queries the candidates for every Geant4 volume.
0055 struct ListCandidates : public SensitiveCandidatesBase {
0056   std::vector<const Acts::Surface*> surfaces;
0057   bool cacheByRegion = false;
0058
0059   std::vector<const Acts::Surface*> queryPosition(
0060       const Acts::GeometryContext& /*gctx*/,
0061       const Acts::Vector3& /*position*/) const override {
0062     return surfaces;
0063   }
0064
0065   std::vector<const Acts::Surface*> queryAll() const override {
0066     return surfaces;
0067   }
0068
0069   const Acts::GeometryObject* queryRegion(
0070       const Acts::GeometryContext& /*gctx*/,
0071       const Acts::Vector3& /*position*/) const override {
0072     // The first surface serves as the key of the single region
0073     if (!cacheByRegion || surfaces.empty()) {
0074       return nullptr;
0075     }
0076     return surfaces.front();
0077   }
0078 };
0079
0080 std::shared_ptr<Acts::Surface> makePlane(const Acts::Vector3& center,
0081                                          double halfLength = 0.5_mm) {
0082   return Acts::Surface::makeShared<Acts::PlaneSurface>(
0083       Acts::Transform3(Acts::Translation3(center)),
0084       std::make_shared<Acts::RectangleBounds>(halfLength, halfLength));
0085 }
0086
0087 /// A Geant4 world with daughters placed without rotation
0088 struct World {
0089   G4Material* material =
0090       G4NistManager::Instance()->FindOrBuildMaterial("G4_Si");
0091   G4LogicalVolume* logical = new G4LogicalVolume(
0092       new G4Box("World", 10_m, 10_m, 10_m), material, "World");
0093   G4VPhysicalVolume* physical = new G4PVPlacement(
0094       nullptr, G4ThreeVector(), logical, "World", nullptr, false, 0);
0095
0096   const G4VPhysicalVolume* place(G4VSolid* solid, const std::string& name,
0097                                  const Acts::Vector3& position) {
0098     auto* volume = new G4LogicalVolume(solid, material, name);
0099     return new G4PVPlacement(
0100         nullptr, G4ThreeVector(position.x(), position.y(), position.z()),
0101         volume, name, logical, false, 0);
0102   }
0103
0104   const G4VPhysicalVolume* placeBox(const std::string& name,
0105                                     const Acts::Vector3& position,
0106                                     double halfLength = 0.5_mm) {
0107     return place(new G4Box(name, halfLength, halfLength, 10_um), name,
0108                  position);
0109   }
0110 };
0111
0112 struct MappingResult {
0113   SensitiveSurfaceMapper::State state;
0114   std::shared_ptr<SensitiveSurfaceMapper> mapper;
0115
0116   /// The surface a Geant4 volume is mapped to, nullptr if none
0117   const Acts::Surface* surfaceOf(const G4VPhysicalVolume* volume) const {
0118     const auto it = state.g4VolumeToSurfaces.find(volume);
0119     if (it == state.g4VolumeToSurfaces.end()) {
0120       return nullptr;
0121     }
0122     BOOST_REQUIRE_EQUAL(it->second.size(), 1u);
0123     return it->second.begin()->second;
0124   }
0125 };
0126
0127 MappingResult runMapper(World& world,
0128                         const std::vector<const Acts::Surface*>& surfaces,
0129                         bool cacheByRegion) {
0130   auto candidates = std::make_shared<ListCandidates>();
0131   candidates->surfaces = surfaces;
0132   candidates->cacheByRegion = cacheByRegion;
0133
0134   SensitiveSurfaceMapper::Config cfg;
0135   cfg.volumeMappings = {"Sensor"};
0136   cfg.candidateSurfaces = candidates;
0137
0138   MappingResult result;
0139   result.mapper = std::make_shared<SensitiveSurfaceMapper>(
0140       cfg,
0141       Acts::getDefaultLogger("SensitiveSurfaceMapper", Acts::Logging::WARNING));
0142   result.mapper->remapSensitiveNames(result.state, gctx, world.physical,
0143                                      Acts::Transform3::Identity());
0144   return result;
0145 }
0146
0147 }  // namespace
0148
0149 BOOST_AUTO_TEST_SUITE(Geant4SensitiveSurfaceMapperSuite)
0150
0151 BOOST_DATA_TEST_CASE(CenterMatchingTolerance,
0152                      boost::unit_test::data::make({false, true}),
0153                      cacheByRegion) {
0154   World world;
0155   // Centers on and next to multiples of 1 mm (and of 10 um)
0156   const std::vector<Acts::Vector3> centers = {
0157       {0., 0., 0.}, {1., 2., -3.}, {1.5, 2.5, -3.5}, {-7.99997, 4., 12.}};
0158   std::vector<std::shared_ptr<Acts::Surface>> planes;
0159   std::vector<const Acts::Surface*> surfaces;
0160   for (const auto& center : centers) {
0161     planes.push_back(makePlane(center));
0162     surfaces.push_back(planes.back().get());
0163   }
0164
0165   std::vector<const G4VPhysicalVolume*> inside;
0166   std::vector<const G4VPhysicalVolume*> outside;
0167   for (const auto& center : centers) {
0168     // Within the tolerance in every component. For the second and last
0169     // center this crosses a multiple of 1 mm
0170     inside.push_back(world.placeBox(
0171         "Sensor", center + Acts::Vector3(-0.5, 0.5, -0.9) * kTolerance));
0172     outside.push_back(world.placeBox(
0173         "Sensor", center + Acts::Vector3(0., 0., 2. * kTolerance)));
0174   }
0175
0176   const auto result = runMapper(world, surfaces, cacheByRegion);
0177   for (std::size_t i = 0; i < centers.size(); ++i) {
0178     BOOST_CHECK_EQUAL(result.surfaceOf(inside[i]), surfaces[i]);
0179     BOOST_CHECK_EQUAL(result.surfaceOf(outside[i]), nullptr);
0180   }
0181   BOOST_CHECK_EQUAL(result.state.missingVolumes.size(), centers.size());
0182 }
0183
0184 BOOST_DATA_TEST_CASE(FarAndNonFiniteCenters,
0185                      boost::unit_test::data::make({false, true}),
0186                      cacheByRegion) {
0187   // Centers far apart, centers just within and too far away to be indexed
0188   // by the mapper (1000 km), and a center with a NaN component, which
0189   // compares equal to everything in that component
0190   World world;
0191   const double nan = std::numeric_limits<double>::quiet_NaN();
0192   const std::vector<Acts::Vector3> centers = {
0193       {0., 0., 0.},  {800_m, -800_m, 800_m}, {1e20, 0., 0.},
0194       {nan, 5., 5.}, {999_km, -1_m, 1_m},    {5e12, 0., 0.}};
0195   std::vector<std::shared_ptr<Acts::Surface>> planes;
0196   std::vector<const Acts::Surface*> surfaces;
0197   for (const auto& center : centers) {
0198     planes.push_back(makePlane(center));
0199     surfaces.push_back(planes.back().get());
0200   }
0201
0202   // Positions next to the centers and far away from them, with the surface
0203   // they are expected to be mapped to
0204   const std::vector<std::pair<Acts::Vector3, const Acts::Surface*>> expected = {
0205       {Acts::Vector3(-0.9, 0.5, -0.5) * kTolerance, surfaces[0]},
0206       {centers[1] + Acts::Vector3(0.9, -0.9, 0.9) * kTolerance, surfaces[1]},
0207       {centers[1] + Acts::Vector3(2., 0., 0.) * kTolerance, nullptr},
0208       {centers[2], surfaces[2]},
0209       {{-30_m, 5., 5. + 0.5 * kTolerance}, surfaces[3]},
0210       {{-30_m, 5., 5. + 2. * kTolerance}, nullptr},
0211       {centers[4] + Acts::Vector3(-0.9, 0.9, -0.9) * kTolerance, surfaces[4]},
0212       {centers[4] + Acts::Vector3(0., 2., 0.) * kTolerance, nullptr},
0213       {centers[5], surfaces[5]}};
0214   std::vector<const G4VPhysicalVolume*> boxes;
0215   for (const auto& entry : expected) {
0216     boxes.push_back(world.placeBox("Sensor", entry.first));
0217   }
0218
0219   const auto result = runMapper(world, surfaces, cacheByRegion);
0220   for (std::size_t i = 0; i < expected.size(); ++i) {
0221     BOOST_CHECK_EQUAL(result.surfaceOf(boxes[i]), expected[i].second);
0222   }
0223 }
0224
0225 BOOST_DATA_TEST_CASE(SharedCenterConcentricCylinders,
0226                      boost::unit_test::data::make({false, true}) *
0227                          boost::unit_test::data::make({false, true}),
0228                      cacheByRegion, reverseOrder) {
0229   // Concentric cylinders around the beam line, like a vertex detector, all
0230   // with the center of their Geant4 tubes
0231   World world;
0232   const std::vector<double> radii = {5_mm, 12_mm, 25_mm};
0233   std::vector<std::shared_ptr<Acts::Surface>> cylinders;
0234   std::vector<const G4VPhysicalVolume*> tubes;
0235   for (const double r : radii) {
0236     cylinders.push_back(Acts::Surface::makeShared<Acts::CylinderSurface>(
0237         Acts::Transform3::Identity(),
0238         std::make_shared<Acts::CylinderBounds>(r + 10_um, 250_mm)));
0239     tubes.push_back(world.place(
0240         new G4Tubs("Tube", r, r + 20_um, 250_mm, 0., 2. * std::numbers::pi),
0241         "Sensor", Acts::Vector3::Zero()));
0242   }
0243   std::vector<const Acts::Surface*> surfaces;
0244   for (const auto& cylinder : cylinders) {
0245     surfaces.push_back(cylinder.get());
0246   }
0247   if (reverseOrder) {
0248     std::ranges::reverse(surfaces);
0249   }
0250
0251   const auto result = runMapper(world, surfaces, cacheByRegion);
0252   for (std::size_t i = 0; i < radii.size(); ++i) {
0253     BOOST_CHECK_EQUAL(result.surfaceOf(tubes[i]), cylinders[i].get());
0254   }
0255 }
0256
0257 BOOST_DATA_TEST_CASE(IdenticalSurfacesKeepCandidateOrder,
0258                      boost::unit_test::data::make({false, true}),
0259                      cacheByRegion) {
0260   World world;
0261   const Acts::Vector3 center{10., -20., 30.};
0262   const auto first = makePlane(center);
0263   const auto second = makePlane(center);
0264   const auto* box = world.placeBox("Sensor", center);
0265
0266   const auto result =
0267       runMapper(world, {first.get(), second.get()}, cacheByRegion);
0268   BOOST_CHECK_EQUAL(result.surfaceOf(box), first.get());
0269 }
0270
0271 BOOST_DATA_TEST_CASE(ManyIdenticalAndCoplanarCenters,
0272                      boost::unit_test::data::make({false, true}),
0273                      cacheByRegion) {
0274   // Many surfaces with an identical center, whose components have an odd
0275   // last mantissa bit, and many centers in one plane at a z value shared by
0276   // hundreds of ALICE 3 sensors
0277   World world;
0278   const Acts::Vector3 shared{std::nextafter(37., 38.), std::nextafter(-41., 0.),
0279                              std::nextafter(43., 44.)};
0280   std::vector<std::shared_ptr<Acts::Surface>> planes;
0281   for (int i = 0; i < 200; ++i) {
0282     planes.push_back(makePlane(shared));
0283   }
0284   const auto* sharedBox = world.placeBox("Sensor", shared);
0285
0286   std::vector<const G4VPhysicalVolume*> inside;
0287   std::vector<const G4VPhysicalVolume*> outside;
0288   for (int ix = 0; ix < 30; ++ix) {
0289     for (int iy = 0; iy < 30; ++iy) {
0290       const Acts::Vector3 center{-100. + 7. * ix, -100. + 7. * iy,
0291                                  1801.0500000000002};
0292       planes.push_back(makePlane(center));
0293       inside.push_back(world.placeBox(
0294           "Sensor", center + Acts::Vector3(0.9, -0.9, 0.9) * kTolerance));
0295       outside.push_back(world.placeBox(
0296           "Sensor", center + Acts::Vector3(0., 0., -2.) * kTolerance));
0297     }
0298   }
0299   std::vector<const Acts::Surface*> surfaces;
0300   for (const auto& plane : planes) {
0301     surfaces.push_back(plane.get());
0302   }
0303
0304   const auto result = runMapper(world, surfaces, cacheByRegion);
0305   BOOST_CHECK_EQUAL(result.surfaceOf(sharedBox), surfaces.front());
0306   for (std::size_t i = 0; i < inside.size(); ++i) {
0307     BOOST_CHECK_EQUAL(result.surfaceOf(inside[i]), surfaces[200 + i]);
0308     BOOST_CHECK_EQUAL(result.surfaceOf(outside[i]), nullptr);
0309   }
0310 }
0311
0312 BOOST_DATA_TEST_CASE(AnnulusCentroidAndCenterMatchFollowCandidateOrder,
0313                      boost::unit_test::data::make({false, true}) *
0314                          boost::unit_test::data::make({false, true}),
0315                      cacheByRegion, annulusFirst) {
0316   // An annulus whose bounds centroid lies inside a large Geant4 box, and a
0317   // plane centered on the same box: whichever comes first in the candidate
0318   // order is the match
0319   World world;
0320   const Acts::Vector3 discCenter{0., 0., 100.};
0321   const auto annulus = Acts::Surface::makeShared<Acts::DiscSurface>(
0322       Acts::Transform3(Acts::Translation3(discCenter)),
0323       std::make_shared<Acts::AnnulusBounds>(380_mm, 500_mm, -0.07, 0.07,
0324                                             Acts::Vector2{2_mm, -1_mm}));
0325   const Acts::Vector3 boxCenter = discCenter + Acts::Vector3(440., 0., 0.);
0326   const auto plane = makePlane(boxCenter);
0327   const auto* box =
0328       world.place(new G4Box("BigBox", 1_m, 1_m, 1_m), "Sensor", boxCenter);
0329
0330   std::vector<const Acts::Surface*> surfaces = {plane.get(), annulus.get()};
0331   if (annulusFirst) {
0332     std::ranges::reverse(surfaces);
0333   }
0334   const auto result = runMapper(world, surfaces, cacheByRegion);
0335   BOOST_CHECK_EQUAL(result.surfaceOf(box), surfaces.front());
0336 }
0337
0338 BOOST_AUTO_TEST_CASE(CheckMappingRequiresOnlySurfacesFlaggedSensitive) {
0339   World world;
0340   const auto sensor = makePlane({0., 0., 0.});
0341   sensor->assignIsSensitive(true);
0342   // E.g. a passive surface with a sensitive geometry id, without Geant4
0343   // counterpart
0344   const auto passive = makePlane({0., 0., 50.});
0345   passive->assignIsSensitive(false);
0346   world.placeBox("Sensor", {0., 0., 0.});
0347
0348   const auto result = runMapper(world, {sensor.get(), passive.get()}, true);
0349   BOOST_CHECK(result.mapper->checkMapping(result.state, gctx));
0350
0351   // A sensor without Geant4 counterpart is still reported
0352   passive->assignIsSensitive(true);
0353   BOOST_CHECK(!result.mapper->checkMapping(result.state, gctx));
0354 }
0355
0356 BOOST_AUTO_TEST_SUITE_END()