Geant4/src/SensitiveSurfaceMapper.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 "ActsExamples/Geant4/SensitiveSurfaceMapper.hpp"
0010
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/Surfaces/AnnulusBounds.hpp"
0013 #include "Acts/Surfaces/SurfaceArray.hpp"
0014 #include "Acts/Utilities/Helpers.hpp"
0015 #include "Acts/Visualization/GeometryView3D.hpp"
0016 #include "Acts/Visualization/ObjVisualization3D.hpp"
0017 #include "ActsExamples/Geant4/AlgebraConverters.hpp"
0018
0019 #include <algorithm>
0020 #include <ostream>
0021 #include <type_traits>
0022 #include <utility>
0023
0024 #include <G4LogicalVolume.hh>
0025 #include <G4Material.hh>
0026 #include <G4Polyhedron.hh>
0027 #include <G4VPhysicalVolume.hh>
0028 #include <G4VSolid.hh>
0029 #include <boost/container/flat_set.hpp>
0030 #include <boost/geometry.hpp>
0031
0032 // Add some type traits for boost::geometry, so we can use the machinery
0033 // directly with Acts::Vector2 / Eigen::Matrix
0034 namespace boost::geometry::traits {
0035
0036 template <typename T, int D>
0037 struct tag<Eigen::Matrix<T, D, 1>> {
0038   using type = point_tag;
0039 };
0040 template <typename T, int D>
0041 struct dimension<Eigen::Matrix<T, D, 1>> : std::integral_constant<int, D> {};
0042 template <typename T, int D>
0043 struct coordinate_type<Eigen::Matrix<T, D, 1>> {
0044   using type = T;
0045 };
0046 template <typename T, int D>
0047 struct coordinate_system<Eigen::Matrix<T, D, 1>> {
0048   using type = boost::geometry::cs::cartesian;
0049 };
0050
0051 template <typename T, int D, std::size_t Index>
0052 struct access<Eigen::Matrix<T, D, 1>, Index> {
0053   static_assert(Index < D, "Out of range");
0054   using Point = Eigen::Matrix<T, D, 1>;
0055   using CoordinateType = typename coordinate_type<Point>::type;
0056   static inline CoordinateType get(Point const& p) { return p[Index]; }
0057   static inline void set(Point& p, CoordinateType const& value) {
0058     p[Index] = value;
0059   }
0060 };
0061
0062 }  // namespace boost::geometry::traits
0063
0064 namespace {
0065
0066 void writeG4Polyhedron(
0067     Acts::IVisualization3D& visualizer, const G4Polyhedron& polyhedron,
0068     const Acts::Transform3& trafo = Acts::Transform3::Identity(),
0069     Acts::Color color = {0, 0, 0}) {
0070   constexpr double convertLength = CLHEP::mm / Acts::UnitConstants::mm;
0071
0072   for (int i = 1; i <= polyhedron.GetNoFacets(); ++i) {
0073     // This is a bit ugly but I didn't find an easy way to compute this
0074     // beforehand.
0075     constexpr std::size_t maxPoints = 1000;
0076     G4Point3D points[maxPoints];
0077     int nPoints = 0;
0078     polyhedron.GetFacet(i, nPoints, points);
0079     assert(static_cast<std::size_t>(nPoints) < maxPoints);
0080
0081     std::vector<Acts::Vector3> faces;
0082     for (int j = 0; j < nPoints; ++j) {
0083       faces.emplace_back(points[j][0] * convertLength,
0084                          points[j][1] * convertLength,
0085                          points[j][2] * convertLength);
0086       faces.back() = trafo * faces.back();
0087     }
0088
0089     visualizer.face(faces, color);
0090   }
0091 }
0092
0093 }  // namespace
0094
0095 namespace ActsExamples::Geant4 {
0096
0097 SensitiveCandidates::SensitiveCandidates(
0098     const std::shared_ptr<const Acts::TrackingGeometry>& trackingGeometry,
0099     std::unique_ptr<const Acts::Logger> _logger)
0100     : m_trackingGeo{trackingGeometry}, m_logger{std::move(_logger)} {}
0101 std::vector<const Acts::Surface*> SensitiveCandidates::queryPosition(
0102     const Acts::GeometryContext& gctx, const Acts::Vector3& position) const {
0103   std::vector<const Acts::Surface*> surfaces{};
0104   ACTS_VERBOSE("Try to fetch the surfaces close to " << position.transpose());
0105
0106   switch (m_trackingGeo->geometryVersion()) {
0107     using enum Acts::TrackingGeometry::GeometryVersion;
0108     case Gen1: {
0109       // In case we do not find a layer at this position for whatever reason
0110       const auto layer = m_trackingGeo->associatedLayer(gctx, position);
0111       if (layer == nullptr) {
0112         return surfaces;
0113       }
0114
0115       const auto surfaceArray = layer->surfaceArray();
0116       if (surfaceArray == nullptr) {
0117         return surfaces;
0118       }
0119
0120       for (const auto& surface : surfaceArray->surfaces()) {
0121         if (surface->isSensitive()) {
0122           surfaces.push_back(surface);
0123         }
0124       }
0125       break;
0126     }
0127     case Gen3: {
0128       const auto* refVolume =
0129           m_trackingGeo->lowestTrackingVolume(gctx, position);
0130       if (refVolume != nullptr) {
0131         constexpr bool restrictToSensitives = true;
0132         refVolume->visitSurfaces(
0133             [&](const Acts::Surface* surface) { surfaces.push_back(surface); },
0134             restrictToSensitives);
0135       }
0136       break;
0137     }
0138   }
0139   return surfaces;
0140 }
0141 std::vector<const Acts::Surface*> SensitiveCandidates::queryAll() const {
0142   std::vector<const Acts::Surface*> surfaces;
0143
0144   constexpr bool restrictToSensitives = true;
0145   m_trackingGeo->visitSurfaces(
0146       [&](auto surface) { surfaces.push_back(surface); }, restrictToSensitives);
0147
0148   return surfaces;
0149 }
0150
0151 SensitiveSurfaceMapper::SensitiveSurfaceMapper(
0152     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0153     : m_cfg(cfg), m_logger(std::move(logger)) {}
0154
0155 void SensitiveSurfaceMapper::remapSensitiveNames(
0156     State& state, const Acts::GeometryContext& gctx,
0157     G4VPhysicalVolume* g4PhysicalVolume,
0158     const Acts::Transform3& motherTransform) const {
0159   // Make sure the unit conversion is correct
0160
0161   auto g4LogicalVolume = g4PhysicalVolume->GetLogicalVolume();
0162   auto g4SensitiveDetector = g4LogicalVolume->GetSensitiveDetector();
0163
0164   // Get the transform of the G4 object
0165   Acts::Transform3 localG4ToGlobal{Acts::Transform3::Identity()};
0166   {
0167     auto g4Translation = g4PhysicalVolume->GetTranslation();
0168     auto g4Rotation = g4PhysicalVolume->GetRotation();
0169     Acts::Vector3 g4RelPosition = convertPosition(g4Translation);
0170     Acts::Translation3 translation(g4RelPosition);
0171     if (g4Rotation == nullptr) {
0172       localG4ToGlobal = motherTransform * translation;
0173     } else {
0174       Acts::RotationMatrix3 rotation;
0175       rotation << g4Rotation->xx(), g4Rotation->yx(), g4Rotation->zx(),
0176           g4Rotation->xy(), g4Rotation->yy(), g4Rotation->zy(),
0177           g4Rotation->xz(), g4Rotation->yz(), g4Rotation->zz();
0178       localG4ToGlobal = motherTransform * (translation * rotation);
0179     }
0180   }
0181
0182   const Acts::Vector3 g4AbsPosition = localG4ToGlobal.translation();
0183
0184   if (G4int nDaughters = g4LogicalVolume->GetNoDaughters(); nDaughters > 0) {
0185     // Step down to all daughters
0186     for (G4int id = 0; id < nDaughters; ++id) {
0187       remapSensitiveNames(state, gctx, g4LogicalVolume->GetDaughter(id),
0188                           localG4ToGlobal);
0189     }
0190   }
0191
0192   const std::string& volumeName{g4LogicalVolume->GetName()};
0193   const std::string& volumeMaterialName{
0194       g4LogicalVolume->GetMaterial()->GetName()};
0195
0196   const bool isSensitive = g4SensitiveDetector != nullptr;
0197   const bool isMappedMaterial =
0198       Acts::rangeContainsValue(m_cfg.materialMappings, volumeMaterialName);
0199   const bool isMappedVolume =
0200       Acts::rangeContainsSubstring(m_cfg.volumeMappings, volumeName);
0201
0202   if (!(isSensitive || isMappedMaterial || isMappedVolume)) {
0203     ACTS_VERBOSE("Did not try mapping '"
0204                  << g4PhysicalVolume->GetName() << "' at "
0205                  << g4AbsPosition.transpose()
0206                  << " because g4SensitiveDetector (=" << g4SensitiveDetector
0207                  << ") is null and volume name (=" << volumeName
0208                  << ") and material name (=" << volumeMaterialName
0209                  << ") were not found");
0210     return;
0211   }
0212   ACTS_VERBOSE("Attempt to map " << g4PhysicalVolume->GetName() << "' at "
0213                                  << g4AbsPosition.transpose()
0214                                  << " to the tracking geometry");
0215
0216   // Prepare the mapped surface
0217   const Acts::Surface* mappedSurface = nullptr;
0218
0219   std::vector<const Acts::Surface*> candidateSurfaces;
0220   const auto g4Polyhedron = g4LogicalVolume->GetSolid()->GetPolyhedron();
0221   for (int i = 1; i < g4Polyhedron->GetNoVertices(); ++i) {
0222     auto vtx = convertPosition(g4Polyhedron->GetVertex(i));
0223     auto vtxGlobal = localG4ToGlobal * vtx;
0224
0225     candidateSurfaces = m_cfg.candidateSurfaces->queryPosition(gctx, vtxGlobal);
0226
0227     if (!candidateSurfaces.empty()) {
0228       break;
0229     }
0230   }
0231
0232   // Fall back to query all surfaces
0233   if (candidateSurfaces.empty()) {
0234     ACTS_DEBUG("No candidate surfaces for volume '" << volumeName << "' at "
0235                                                     << g4AbsPosition.transpose()
0236                                                     << ", query all surfaces");
0237     candidateSurfaces = m_cfg.candidateSurfaces->queryAll();
0238   }
0239
0240   ACTS_VERBOSE("Found " << candidateSurfaces.size()
0241                         << " candidate surfaces for " << volumeName);
0242
0243   Acts::detail::TransformComparator trfSorter{};
0244   for (const auto& candidateSurface : candidateSurfaces) {
0245     if (trfSorter.compare<3>(candidateSurface->center(gctx), g4AbsPosition) ==
0246         0) {
0247       ACTS_DEBUG("Successful match with center: "
0248                  << candidateSurface->center(gctx).transpose()
0249                  << ", G4-position: " << g4AbsPosition.transpose());
0250       mappedSurface = candidateSurface;
0251       break;
0252     } else if (candidateSurface->bounds().type() ==
0253                Acts::SurfaceBounds::eAnnulus) {
0254       const auto& bounds =
0255           *static_cast<const Acts::AnnulusBounds*>(&candidateSurface->bounds());
0256
0257       const auto vertices = bounds.vertices(0);
0258
0259       constexpr bool clockwise = false;
0260       constexpr bool closed = false;
0261       using Polygon =
0262           boost::geometry::model::polygon<Acts::Vector2, clockwise, closed>;
0263
0264       Polygon poly;
0265       boost::geometry::assign_points(poly, vertices);
0266
0267       Acts::Vector2 boundsCentroidSurfaceFrame = Acts::Vector2::Zero();
0268       boost::geometry::centroid(poly, boundsCentroidSurfaceFrame);
0269
0270       Acts::Vector3 boundsCentroidGlobal{boundsCentroidSurfaceFrame[0],
0271                                          boundsCentroidSurfaceFrame[1], 0.0};
0272       boundsCentroidGlobal =
0273           candidateSurface->localToGlobalTransform(gctx) * boundsCentroidGlobal;
0274
0275       const auto boundsCentroidG4Frame =
0276           localG4ToGlobal.inverse() * boundsCentroidGlobal;
0277
0278       if (g4LogicalVolume->GetSolid()->Inside(
0279               convertPosition(boundsCentroidG4Frame)) != EInside::kOutside) {
0280         ACTS_VERBOSE("Successful match with centroid matching");
0281         mappedSurface = candidateSurface;
0282         break;
0283       }
0284     }
0285   }
0286
0287   if (mappedSurface == nullptr) {
0288     ACTS_DEBUG("No mapping found for '"
0289                << volumeName << "' with material '" << volumeMaterialName
0290                << "' at position " << g4AbsPosition.transpose());
0291     state.missingVolumes.emplace_back(g4PhysicalVolume, localG4ToGlobal);
0292     return;
0293   }
0294
0295   // A mapped surface was found, a new name will be set that G4PhysVolume
0296   ACTS_DEBUG("Matched " << volumeName << " to " << mappedSurface->geometryId()
0297                         << " at position " << g4AbsPosition.transpose());
0298   // Check if the prefix is not yet assigned
0299   if (volumeName.find(mappingPrefix) == std::string::npos) {
0300     // Set the new name
0301     std::string mappedName = std::string(mappingPrefix) + volumeName;
0302     g4PhysicalVolume->SetName(mappedName);
0303   }
0304   if (state.g4VolumeToSurfaces.find(g4PhysicalVolume) ==
0305       state.g4VolumeToSurfaces.end()) {
0306     state.g4VolumeToSurfaces.insert(
0307         std::make_pair(g4PhysicalVolume, SurfacePosMap_t{trfSorter}));
0308   }
0309   // Insert into the multi-map
0310   if (!state.g4VolumeToSurfaces[g4PhysicalVolume]
0311            .insert(std::make_pair(g4AbsPosition, mappedSurface))
0312            .second) {
0313     ACTS_WARNING("Duplicate surface found for " << volumeName << " @ "
0314                                                 << g4AbsPosition.transpose());
0315   }
0316 }
0317
0318 bool SensitiveSurfaceMapper::checkMapping(
0319     const State& state, const Acts::GeometryContext& gctx,
0320     bool writeMissingG4VolsAsObj, bool writeMissingSurfacesAsObj) const {
0321   auto allSurfaces = m_cfg.candidateSurfaces->queryAll();
0322   std::ranges::sort(allSurfaces);
0323
0324   std::vector<const Acts::Surface*> found;
0325   for (const auto& [_, surfaceMap] : state.g4VolumeToSurfaces) {
0326     for (const auto& [__, surfacePtr] : surfaceMap) {
0327       found.push_back(surfacePtr);
0328     }
0329   }
0330   std::ranges::sort(found);
0331   auto newEnd = std::unique(found.begin(), found.end());
0332   found.erase(newEnd, found.end());
0333
0334   std::vector<const Acts::Surface*> missing;
0335   std::set_difference(allSurfaces.begin(), allSurfaces.end(), found.begin(),
0336                       found.end(), std::back_inserter(missing));
0337
0338   ACTS_INFO("Number of overall sensitive surfaces: " << allSurfaces.size());
0339   ACTS_INFO("Number of mapped volume->surface mappings: " << found.size());
0340   ACTS_INFO(
0341       "Number of sensitive surfaces that are not mapped: " << missing.size());
0342   ACTS_INFO("Number of G4 volumes without a matching Surface: "
0343             << state.missingVolumes.size());
0344
0345   if (writeMissingG4VolsAsObj) {
0346     Acts::ObjVisualization3D visualizer;
0347     for (const auto& [g4vol, trafo] : state.missingVolumes) {
0348       auto polyhedron = g4vol->GetLogicalVolume()->GetSolid()->GetPolyhedron();
0349       writeG4Polyhedron(visualizer, *polyhedron, trafo);
0350     }
0351
0352     std::ofstream os("missing_g4_volumes.obj");
0353     visualizer.write(os);
0354   }
0355
0356   if (writeMissingSurfacesAsObj) {
0357     Acts::ObjVisualization3D visualizer;
0358     Acts::ViewConfig vcfg;
0359     vcfg.quarterSegments = 720;
0360     for (auto srf : missing) {
0361       Acts::GeometryView3D::drawSurface(visualizer, *srf, gctx,
0362                                         Acts::Transform3::Identity(), vcfg);
0363     }
0364
0365     std::ofstream os("missing_acts_surfaces.obj");
0366     visualizer.write(os);
0367   }
0368
0369   return missing.empty();
0370 }
0371
0372 }  // namespace ActsExamples::Geant4