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