Geant4/src/SensitiveSurfaceMapper.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2021 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsExamples/Geant4/SensitiveSurfaceMapper.hpp"
0010
0011 #include "Acts/Definitions/Units.hpp"
0012 #include "Acts/Geometry/GeometryContext.hpp"
0013 #include "Acts/Geometry/GeometryIdentifier.hpp"
0014 #include "Acts/Geometry/Layer.hpp"
0015 #include "Acts/Geometry/TrackingGeometry.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017 #include "Acts/Surfaces/SurfaceArray.hpp"
0018
0019 #include <algorithm>
0020 #include <ostream>
0021 #include <stdexcept>
0022 #include <utility>
0023
0024 #include <G4LogicalVolume.hh>
0025 #include <G4Material.hh>
0026 #include <G4VPhysicalVolume.hh>
0027
0028 ActsExamples::SensitiveSurfaceMapper::SensitiveSurfaceMapper(
0029     const Config& cfg, std::unique_ptr<const Acts::Logger> logger)
0030     : m_cfg(cfg), m_logger(std::move(logger)) {
0031   if (m_cfg.trackingGeometry == nullptr) {
0032     throw std::invalid_argument("No Acts::TrackingGeometry provided.");
0033   }
0034 }
0035
0036 void ActsExamples::SensitiveSurfaceMapper::remapSensitiveNames(
0037     G4VPhysicalVolume* g4PhysicalVolume,
0038     const Acts::Transform3& motherTransform, int& sCounter) const {
0039   constexpr double convertLength = CLHEP::mm / Acts::UnitConstants::mm;
0040
0041   auto g4LogicalVolume = g4PhysicalVolume->GetLogicalVolume();
0042   auto g4SensitiveDetector = g4LogicalVolume->GetSensitiveDetector();
0043
0044   // Get the transform of the G4 object
0045   auto g4Translation = g4PhysicalVolume->GetTranslation();
0046   auto g4Rotation = g4PhysicalVolume->GetRotation();
0047   Acts::Vector3 g4RelPosition(g4Translation[0] * convertLength,
0048                               g4Translation[1] * convertLength,
0049                               g4Translation[2] * convertLength);
0050   Acts::Translation3 translation(g4RelPosition);
0051   Acts::Transform3 transform;
0052   if (g4Rotation == nullptr) {
0053     transform = motherTransform * translation;
0054   } else {
0055     Acts::RotationMatrix3 rotation;
0056     rotation << g4Rotation->xx(), g4Rotation->yx(), g4Rotation->zx(),
0057         g4Rotation->xy(), g4Rotation->yy(), g4Rotation->zy(), g4Rotation->xz(),
0058         g4Rotation->yz(), g4Rotation->zz();
0059     transform = motherTransform * (translation * rotation);
0060   }
0061   Acts::Vector3 g4AbsPosition = transform * Acts::Vector3::Zero();
0062
0063   if (G4int nDaughters = g4LogicalVolume->GetNoDaughters(); nDaughters > 0) {
0064     // Step down to all daughters
0065     for (G4int id = 0; id < nDaughters; ++id) {
0066       remapSensitiveNames(g4LogicalVolume->GetDaughter(id), transform,
0067                           sCounter);
0068     }
0069     return;
0070   }
0071
0072   std::string volumeName = g4LogicalVolume->GetName();
0073   std::string volumeMaterialName = g4LogicalVolume->GetMaterial()->GetName();
0074
0075   bool isSensitive = g4SensitiveDetector != nullptr;
0076   bool isMappedMaterial =
0077       std::find(m_cfg.materialMappings.begin(), m_cfg.materialMappings.end(),
0078                 volumeMaterialName) != m_cfg.materialMappings.end();
0079   bool isMappedVolume =
0080       std::find(m_cfg.volumeMappings.begin(), m_cfg.volumeMappings.end(),
0081                 volumeName) != m_cfg.volumeMappings.end();
0082
0083   if (isSensitive || isMappedMaterial || isMappedVolume) {
0084     // Find the associated ACTS object
0085     auto actsLayer = m_cfg.trackingGeometry->associatedLayer(
0086         Acts::GeometryContext(), g4AbsPosition);
0087
0088     // Prepare the mapped surface
0089     const Acts::Surface* mappedSurface = nullptr;
0090
0091     if (actsLayer != nullptr && actsLayer->surfaceArray() != nullptr) {
0092       auto actsSurfaces = actsLayer->surfaceArray()->at(g4AbsPosition);
0093       if (!actsSurfaces.empty()) {
0094         // Fast matching: search
0095         for (const auto& as : actsSurfaces) {
0096           if (as->center(Acts::GeometryContext()).isApprox(g4AbsPosition)) {
0097             mappedSurface = as;
0098             break;
0099           }
0100         }
0101       }
0102       if (mappedSurface == nullptr) {
0103         // Slow matching: Fallback, loop over all layer surfaces
0104         for (const auto& as : actsLayer->surfaceArray()->surfaces()) {
0105           if (as->center(Acts::GeometryContext()).isApprox(g4AbsPosition)) {
0106             mappedSurface = as;
0107             break;
0108           }
0109         }
0110       }
0111     }
0112     // A mapped surface was found, a new name will be set that
0113     // contains the GeometryID/
0114     if (mappedSurface != nullptr) {
0115       ++sCounter;
0116       std::string mappedVolumeName(SensitiveSurfaceMapper::mappingPrefix);
0117       mappedVolumeName += std::to_string(mappedSurface->geometryId().value());
0118       ACTS_VERBOSE("Found matching surface " << mappedSurface->geometryId()
0119                                              << " at position "
0120                                              << g4RelPosition.transpose());
0121       ACTS_VERBOSE("Remap: " << g4PhysicalVolume->GetName() << " -> "
0122                              << mappedVolumeName);
0123       g4PhysicalVolume->SetName(mappedVolumeName.c_str());
0124     } else {
0125       ACTS_VERBOSE("No mapping found for '"
0126                    << volumeName << "' with material '" << volumeMaterialName
0127                    << "' at position " << g4RelPosition.transpose());
0128     }
0129   } else {
0130     ACTS_VERBOSE("Did not try mapping '"
0131                  << g4PhysicalVolume->GetName() << "' at "
0132                  << g4RelPosition.transpose()
0133                  << " because g4SensitiveDetector (=" << g4SensitiveDetector
0134                  << ") is null and volume name (=" << volumeName
0135                  << ") and material name (=" << volumeMaterialName
0136                  << ") were not found");
0137   }
0138 }