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 }  // namespace
0093
0094 namespace ActsExamples::Geant4 {
0095
0096 SensitiveCandidates::SensitiveCandidates(
0097     const std::shared_ptr<const Acts::TrackingGeometry>& trackingGeometry,
0098     std::unique_ptr<const Acts::Logger> _logger)
0099     : m_trackingGeo{trackingGeometry}, m_logger{std::move(_logger)} {}
0100 std::vector<const Acts::Surface*> SensitiveCandidates::queryPosition(
0101     const Acts::GeometryContext& gctx, const Acts::Vector3& position) const {
0102   std::vector<const Acts::Surface*> surfaces{};
0103   ACTS_VERBOSE("Try to fetch the surfaces close to " << position.transpose());
0104
0105   switch (m_trackingGeo->geometryVersion()) {
0106     using enum Acts::TrackingGeometry::GeometryVersion;
0107     case Gen1: {
0108       // In case we do not find a layer at this position for whatever reason
0109       const auto layer = m_trackingGeo->associatedLayer(gctx, position);
0110       if (layer == nullptr) {
0111         return surfaces;
0112       }
0113
0114       const auto surfaceArray = layer->surfaceArray();
0115       if (surfaceArray == nullptr) {
0116         return surfaces;
0117       }
0118
0119       for (const auto& surface : surfaceArray->surfaces()) {
0120         if (surface->associatedDetectorElement() != nullptr) {
0121           surfaces.push_back(surface);
0122         }
0123       }
0124       break;
0125     }
0126     case Gen3: {
0127       const auto* refVolume =
0128           m_trackingGeo->lowestTrackingVolume(gctx, position);
0129       if (refVolume != nullptr) {
0130         constexpr bool restrictToSensitives = true;
0131         refVolume->visitSurfaces(
0132             [&](const Acts::Surface* surface) { surfaces.push_back(surface); },
0133             restrictToSensitives);
0134       }
0135       break;
0136     }
0137   }
0138   return surfaces;
0139 }
0140 std::vector<const Acts::Surface*> SensitiveCandidates::queryAll() const {
0141   std::vector<const Acts::Surface*> surfaces;
0142
0143   constexpr bool restrictToSensitives = true;
0144   m_trackingGeo->visitSurfaces(
0145       [&](auto surface) { surfaces.push_back(surface); }, restrictToSensitives);
0146
0147   return surfaces;
0148 }
0149
0150 SensitiveSurfaceMapper::SensitiveSurfaceMapper(
0151     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0152     : m_cfg(cfg), m_logger(std::move(logger)) {}
0153
0154 void SensitiveSurfaceMapper::remapSensitiveNames(
0155     State& state, const Acts::GeometryContext& gctx,
0156     G4VPhysicalVolume* g4PhysicalVolume,
0157     const Acts::Transform3& motherTransform) const {
0158   // Make sure the unit conversion is correct
0159
0160   auto g4LogicalVolume = g4PhysicalVolume->GetLogicalVolume();
0161   auto g4SensitiveDetector = g4LogicalVolume->GetSensitiveDetector();
0162
0163   // Get the transform of the G4 object
0164   Acts::Transform3 localG4ToGlobal{Acts::Transform3::Identity()};
0165   {
0166     auto g4Translation = g4PhysicalVolume->GetTranslation();
0167     auto g4Rotation = g4PhysicalVolume->GetRotation();
0168     Acts::Vector3 g4RelPosition = convertPosition(g4Translation);
0169     Acts::Translation3 translation(g4RelPosition);
0170     if (g4Rotation == nullptr) {
0171       localG4ToGlobal = motherTransform * translation;
0172     } else {
0173       Acts::RotationMatrix3 rotation;
0174       rotation << g4Rotation->xx(), g4Rotation->yx(), g4Rotation->zx(),
0175           g4Rotation->xy(), g4Rotation->yy(), g4Rotation->zy(),
0176           g4Rotation->xz(), g4Rotation->yz(), g4Rotation->zz();
0177       localG4ToGlobal = motherTransform * (translation * rotation);
0178     }
0179   }
0180
0181   const Acts::Vector3 g4AbsPosition = localG4ToGlobal.translation();
0182
0183   if (G4int nDaughters = g4LogicalVolume->GetNoDaughters(); nDaughters > 0) {
0184     // Step down to all daughters
0185     for (G4int id = 0; id < nDaughters; ++id) {
0186       remapSensitiveNames(state, gctx, g4LogicalVolume->GetDaughter(id),
0187                           localG4ToGlobal);
0188     }
0189   }
0190
0191   const std::string& volumeName{g4LogicalVolume->GetName()};
0192   const std::string& volumeMaterialName{
0193       g4LogicalVolume->GetMaterial()->GetName()};
0194
0195   const bool isSensitive = g4SensitiveDetector != nullptr;
0196   const bool isMappedMaterial =
0197       Acts::rangeContainsValue(m_cfg.materialMappings, volumeMaterialName);
0198   const bool isMappedVolume =
0199       Acts::rangeContainsValue(m_cfg.volumeMappings, volumeName);
0200
0201   if (!(isSensitive || isMappedMaterial || isMappedVolume)) {
0202     ACTS_VERBOSE("Did not try mapping '"
0203                  << g4PhysicalVolume->GetName() << "' at "
0204                  << g4AbsPosition.transpose()
0205                  << " because g4SensitiveDetector (=" << g4SensitiveDetector
0206                  << ") is null and volume name (=" << volumeName
0207                  << ") and material name (=" << volumeMaterialName
0208                  << ") were not found");
0209     return;
0210   }
0211   ACTS_VERBOSE("Attempt to map " << g4PhysicalVolume->GetName() << "' at "
0212                                  << g4AbsPosition.transpose()
0213                                  << " to the tracking geometry");
0214
0215   // Prepare the mapped surface
0216   const Acts::Surface* mappedSurface = nullptr;
0217
0218   std::vector<const Acts::Surface*> candidateSurfaces;
0219   const auto g4Polyhedron = g4LogicalVolume->GetSolid()->GetPolyhedron();
0220   for (int i = 1; i < g4Polyhedron->GetNoVertices(); ++i) {
0221     auto vtx = convertPosition(g4Polyhedron->GetVertex(i));
0222     auto vtxGlobal = localG4ToGlobal * vtx;
0223
0224     candidateSurfaces = m_cfg.candidateSurfaces->queryPosition(gctx, vtxGlobal);
0225
0226     if (!candidateSurfaces.empty()) {
0227       break;
0228     }
0229   }
0230
0231   // Fall back to query all surfaces
0232   if (candidateSurfaces.empty()) {
0233     ACTS_DEBUG("No candidate surfaces for volume '" << volumeName << "' at "
0234                                                     << g4AbsPosition.transpose()
0235                                                     << ", query all surfaces");
0236     candidateSurfaces = m_cfg.candidateSurfaces->queryAll();
0237   }
0238
0239   ACTS_VERBOSE("Found " << candidateSurfaces.size()
0240                         << " candidate surfaces for " << volumeName);
0241
0242   Acts::detail::TransformComparator trfSorter{};
0243   for (const auto& candidateSurface : candidateSurfaces) {
0244     if (trfSorter.compare<3>(candidateSurface->center(gctx), g4AbsPosition) ==
0245         0) {
0246       ACTS_DEBUG("Successful match with center: "
0247                  << candidateSurface->center(gctx).transpose()
0248                  << ", G4-position: " << g4AbsPosition.transpose());
0249       mappedSurface = candidateSurface;
0250       break;
0251     } else if (candidateSurface->bounds().type() ==
0252                Acts::SurfaceBounds::eAnnulus) {
0253       const auto& bounds =
0254           *static_cast<const Acts::AnnulusBounds*>(&candidateSurface->bounds());
0255
0256       const auto vertices = bounds.vertices(0);
0257
0258       constexpr bool clockwise = false;
0259       constexpr bool closed = false;
0260       using Polygon =
0261           boost::geometry::model::polygon<Acts::Vector2, clockwise, closed>;
0262
0263       Polygon poly;
0264       boost::geometry::assign_points(poly, vertices);
0265
0266       Acts::Vector2 boundsCentroidSurfaceFrame = Acts::Vector2::Zero();
0267       boost::geometry::centroid(poly, boundsCentroidSurfaceFrame);
0268
0269       Acts::Vector3 boundsCentroidGlobal{boundsCentroidSurfaceFrame[0],
0270                                          boundsCentroidSurfaceFrame[1], 0.0};
0271       boundsCentroidGlobal =
0272           candidateSurface->transform(gctx) * boundsCentroidGlobal;
0273
0274       const auto boundsCentroidG4Frame =
0275           localG4ToGlobal.inverse() * boundsCentroidGlobal;
0276
0277       if (g4LogicalVolume->GetSolid()->Inside(
0278               convertPosition(boundsCentroidG4Frame)) != EInside::kOutside) {
0279         ACTS_VERBOSE("Successful match with centroid matching");
0280         mappedSurface = candidateSurface;
0281         break;
0282       }
0283     }
0284   }
0285
0286   if (mappedSurface == nullptr) {
0287     ACTS_DEBUG("No mapping found for '"
0288                << volumeName << "' with material '" << volumeMaterialName
0289                << "' at position " << g4AbsPosition.transpose());
0290     state.missingVolumes.emplace_back(g4PhysicalVolume, localG4ToGlobal);
0291     return;
0292   }
0293
0294   // A mapped surface was found, a new name will be set that G4PhysVolume
0295   ACTS_DEBUG("Matched " << volumeName << " to " << mappedSurface->geometryId()
0296                         << " at position " << g4AbsPosition.transpose());
0297   // Check if the prefix is not yet assigned
0298   if (volumeName.find(mappingPrefix) == std::string::npos) {
0299     // Set the new name
0300     std::string mappedName = std::string(mappingPrefix) + volumeName;
0301     g4PhysicalVolume->SetName(mappedName);
0302   }
0303   if (state.g4VolumeToSurfaces.find(g4PhysicalVolume) ==
0304       state.g4VolumeToSurfaces.end()) {
0305     state.g4VolumeToSurfaces.insert(
0306         std::make_pair(g4PhysicalVolume, SurfacePosMap_t{trfSorter}));
0307   }
0308   // Insert into the multi-map
0309   if (!state.g4VolumeToSurfaces[g4PhysicalVolume]
0310            .insert(std::make_pair(g4AbsPosition, mappedSurface))
0311            .second) {
0312     ACTS_WARNING("Duplicate surface found for " << volumeName << " @ "
0313                                                 << g4AbsPosition.transpose());
0314   }
0315 }
0316
0317 bool SensitiveSurfaceMapper::checkMapping(
0318     const State& state, const Acts::GeometryContext& gctx,
0319     bool writeMissingG4VolsAsObj, bool writeMissingSurfacesAsObj) const {
0320   auto allSurfaces = m_cfg.candidateSurfaces->queryAll();
0321   std::ranges::sort(allSurfaces);
0322
0323   std::vector<const Acts::Surface*> found;
0324   for (const auto& [_, surfaceMap] : state.g4VolumeToSurfaces) {
0325     for (const auto& [__, surfacePtr] : surfaceMap) {
0326       found.push_back(surfacePtr);
0327     }
0328   }
0329   std::ranges::sort(found);
0330   auto newEnd = std::unique(found.begin(), found.end());
0331   found.erase(newEnd, found.end());
0332
0333   std::vector<const Acts::Surface*> missing;
0334   std::set_difference(allSurfaces.begin(), allSurfaces.end(), found.begin(),
0335                       found.end(), std::back_inserter(missing));
0336
0337   ACTS_INFO("Number of overall sensitive surfaces: " << allSurfaces.size());
0338   ACTS_INFO("Number of mapped volume->surface mappings: " << found.size());
0339   ACTS_INFO(
0340       "Number of sensitive surfaces that are not mapped: " << missing.size());
0341   ACTS_INFO("Number of G4 volumes without a matching Surface: "
0342             << state.missingVolumes.size());
0343
0344   if (writeMissingG4VolsAsObj) {
0345     Acts::ObjVisualization3D visualizer;
0346     for (const auto& [g4vol, trafo] : state.missingVolumes) {
0347       auto polyhedron = g4vol->GetLogicalVolume()->GetSolid()->GetPolyhedron();
0348       writeG4Polyhedron(visualizer, *polyhedron, trafo);
0349     }
0350
0351     std::ofstream os("missing_g4_volumes.obj");
0352     visualizer.write(os);
0353   }
0354
0355   if (writeMissingSurfacesAsObj) {
0356     Acts::ObjVisualization3D visualizer;
0357     Acts::ViewConfig vcfg;
0358     vcfg.quarterSegments = 720;
0359     for (auto srf : missing) {
0360       Acts::GeometryView3D::drawSurface(visualizer, *srf, gctx,
0361                                         Acts::Transform3::Identity(), vcfg);
0362     }
0363
0364     std::ofstream os("missing_acts_surfaces.obj");
0365     visualizer.write(os);
0366   }
0367
0368   return missing.empty();
0369 }
0370
0371 }  // namespace ActsExamples::Geant4