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{Acts::Transform3::Identity()};
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   const Acts::Vector3 g4AbsPosition = localG4ToGlobal.translation();
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   }
0186
0187   const std::string& volumeName{g4LogicalVolume->GetName()};
0188   const std::string& volumeMaterialName{
0189       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   ACTS_VERBOSE("Attempt to map " << g4PhysicalVolume->GetName() << "' at "
0208                                  << g4AbsPosition.transpose()
0209                                  << " to the tracking geometry");
0210
0211   // Prepare the mapped surface
0212   const Acts::Surface* mappedSurface = nullptr;
0213
0214   std::vector<const Acts::Surface*> candidateSurfaces;
0215   const auto g4Polyhedron = g4LogicalVolume->GetSolid()->GetPolyhedron();
0216   for (int i = 1; i < g4Polyhedron->GetNoVertices(); ++i) {
0217     auto vtx = g4ToActsVector(g4Polyhedron->GetVertex(i));
0218     auto vtxGlobal = localG4ToGlobal * vtx;
0219
0220     candidateSurfaces = m_cfg.candidateSurfaces->queryPosition(gctx, vtxGlobal);
0221
0222     if (!candidateSurfaces.empty()) {
0223       break;
0224     }
0225   }
0226
0227   // Fall back to query all surfaces
0228   if (candidateSurfaces.empty()) {
0229     ACTS_DEBUG("No candidate surfaces for volume '" << volumeName << "' at "
0230                                                     << g4AbsPosition.transpose()
0231                                                     << ", query all surfaces");
0232     candidateSurfaces = m_cfg.candidateSurfaces->queryAll();
0233   }
0234
0235   ACTS_VERBOSE("Found " << candidateSurfaces.size()
0236                         << " candidate surfaces for " << volumeName);
0237
0238   for (const auto& candidateSurface : candidateSurfaces) {
0239     if (candidateSurface->center(gctx).isApprox(g4AbsPosition)) {
0240       ACTS_VERBOSE("Successful match with center matching");
0241       mappedSurface = candidateSurface;
0242       break;
0243     } else if (candidateSurface->bounds().type() ==
0244                Acts::SurfaceBounds::eAnnulus) {
0245       const auto& bounds =
0246           *static_cast<const Acts::AnnulusBounds*>(&candidateSurface->bounds());
0247
0248       const auto vertices = bounds.vertices(0);
0249
0250       constexpr bool clockwise = false;
0251       constexpr bool closed = false;
0252       using Polygon =
0253           boost::geometry::model::polygon<Acts::Vector2, clockwise, closed>;
0254
0255       Polygon poly;
0256       boost::geometry::assign_points(poly, vertices);
0257
0258       Acts::Vector2 boundsCentroidSurfaceFrame = Acts::Vector2::Zero();
0259       boost::geometry::centroid(poly, boundsCentroidSurfaceFrame);
0260
0261       Acts::Vector3 boundsCentroidGlobal{boundsCentroidSurfaceFrame[0],
0262                                          boundsCentroidSurfaceFrame[1], 0.0};
0263       boundsCentroidGlobal =
0264           candidateSurface->transform(gctx) * boundsCentroidGlobal;
0265
0266       const auto boundsCentroidG4Frame =
0267           localG4ToGlobal.inverse() * boundsCentroidGlobal;
0268
0269       if (g4LogicalVolume->GetSolid()->Inside(
0270               actsToG4Vec(boundsCentroidG4Frame)) != EInside::kOutside) {
0271         ACTS_VERBOSE("Successful match with centroid matching");
0272         mappedSurface = candidateSurface;
0273         break;
0274       }
0275     }
0276   }
0277
0278   if (mappedSurface == nullptr) {
0279     ACTS_DEBUG("No mapping found for '"
0280                << volumeName << "' with material '" << volumeMaterialName
0281                << "' at position " << g4AbsPosition.transpose());
0282     state.missingVolumes.emplace_back(g4PhysicalVolume, localG4ToGlobal);
0283     return;
0284   }
0285
0286   // A mapped surface was found, a new name will be set that G4PhysVolume
0287   ACTS_DEBUG("Matched " << volumeName << " to " << mappedSurface->geometryId()
0288                         << " at position " << g4AbsPosition.transpose());
0289   // Check if the prefix is not yet assigned
0290   if (volumeName.find(mappingPrefix) == std::string::npos) {
0291     // Set the new name
0292     std::string mappedName = std::string(mappingPrefix) + volumeName;
0293     g4PhysicalVolume->SetName(mappedName);
0294   }
0295   // Insert into the multi-map
0296   state.g4VolumeToSurfaces.insert({g4PhysicalVolume, mappedSurface});
0297 }
0298
0299 bool SensitiveSurfaceMapper::checkMapping(
0300     const State& state, const Acts::GeometryContext& gctx,
0301     bool writeMissingG4VolsAsObj, bool writeMissingSurfacesAsObj) const {
0302   auto allSurfaces = m_cfg.candidateSurfaces->queryAll();
0303   std::ranges::sort(allSurfaces);
0304
0305   std::vector<const Acts::Surface*> found;
0306   for (const auto [_, surfacePtr] : state.g4VolumeToSurfaces) {
0307     found.push_back(surfacePtr);
0308   }
0309   std::ranges::sort(found);
0310   auto newEnd = std::unique(found.begin(), found.end());
0311   found.erase(newEnd, found.end());
0312
0313   std::vector<const Acts::Surface*> missing;
0314   std::set_difference(allSurfaces.begin(), allSurfaces.end(), found.begin(),
0315                       found.end(), std::back_inserter(missing));
0316
0317   ACTS_INFO("Number of overall sensitive surfaces: " << allSurfaces.size());
0318   ACTS_INFO("Number of mapped volume->surface mappings: "
0319             << state.g4VolumeToSurfaces.size());
0320   ACTS_INFO(
0321       "Number of sensitive surfaces that are not mapped: " << missing.size());
0322   ACTS_INFO("Number of G4 volumes without a matching Surface: "
0323             << state.missingVolumes.size());
0324
0325   if (writeMissingG4VolsAsObj) {
0326     Acts::ObjVisualization3D visualizer;
0327     for (const auto& [g4vol, trafo] : state.missingVolumes) {
0328       auto polyhedron = g4vol->GetLogicalVolume()->GetSolid()->GetPolyhedron();
0329       writeG4Polyhedron(visualizer, *polyhedron, trafo);
0330     }
0331
0332     std::ofstream os("missing_g4_volumes.obj");
0333     visualizer.write(os);
0334   }
0335
0336   if (writeMissingSurfacesAsObj) {
0337     Acts::ObjVisualization3D visualizer;
0338     Acts::ViewConfig vcfg;
0339     vcfg.quarterSegments = 720;
0340     for (auto srf : missing) {
0341       Acts::GeometryView3D::drawSurface(visualizer, *srf, gctx,
0342                                         Acts::Transform3::Identity(), vcfg);
0343     }
0344
0345     std::ofstream os("missing_acts_surfaces.obj");
0346     visualizer.write(os);
0347   }
0348
0349   return missing.empty();
0350 }
0351
0352 }  // namespace ActsExamples::Geant4