algorithms/tracking/ActsGeometryProvider.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2024 Whitney Armstrong, Wouter Deconinck, Dmitry Romanov
0003
0004 #if Acts_VERSION_MAJOR < 45
0005 #include <Acts/Geometry/DetectorElementBase.hpp>
0006 #endif
0007 #include <Acts/Geometry/GeometryIdentifier.hpp>
0008 #include <Acts/Geometry/TrackingGeometry.hpp>
0009 #include <Acts/Geometry/TrackingVolume.hpp>
0010 #include <Acts/MagneticField/MagneticFieldContext.hpp>
0011 #include <Acts/Material/IMaterialDecorator.hpp>
0012 #include <Acts/Utilities/Logger.hpp>
0013 #include <boost/container/detail/std_fwd.hpp>
0014 #include <fmt/format.h>
0015 #if __has_include(<ActsPlugins/DD4hep/ConvertDD4hepDetector.hpp>)
0016 #include <ActsPlugins/DD4hep/ConvertDD4hepDetector.hpp>
0017 #include <ActsPlugins/DD4hep/DD4hepDetectorElement.hpp>
0018 #include <ActsPlugins/DD4hep/DD4hepFieldAdapter.hpp>
0019 #include <ActsPlugins/Json/JsonMaterialDecorator.hpp>
0020 #include <ActsPlugins/Json/MaterialMapJsonConverter.hpp>
0021 #else
0022 #include <Acts/Plugins/DD4hep/ConvertDD4hepDetector.hpp>
0023 #include <Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp>
0024 #include <Acts/Plugins/DD4hep/DD4hepFieldAdapter.hpp>
0025 #include <Acts/Plugins/Json/JsonMaterialDecorator.hpp>
0026 #include <Acts/Plugins/Json/MaterialMapJsonConverter.hpp>
0027 #endif
0028 #include <Acts/Surfaces/Surface.hpp>
0029 #include <Acts/Utilities/BinningType.hpp>
0030 #include <Acts/Utilities/Result.hpp>
0031 #include <Acts/Visualization/GeometryView3D.hpp>
0032 #include <Acts/Visualization/ObjVisualization3D.hpp>
0033 #include <Acts/Visualization/PlyVisualization3D.hpp>
0034 #include <DD4hep/DetElement.h>
0035 #include <DD4hep/VolumeManager.h>
0036 #include <TGeoManager.h>
0037 #include <fmt/ostream.h>
0038 #include <spdlog/common.h>
0039 // Formatter for Eigen matrices
0040 #include <Eigen/Core>
0041 #include <exception>
0042 #include <filesystem>
0043 #include <functional>
0044 #include <initializer_list>
0045 #include <set>
0046 #include <type_traits>
0047 #include <utility>
0048
0049 #include "ActsGeometryProvider.h"
0050 #include "extensions/spdlog/SpdlogToActs.h"
0051
0052 template <typename T>
0053 struct fmt::formatter<T, std::enable_if_t<std::is_base_of_v<Eigen::MatrixBase<T>, T>, char>>
0054     : fmt::ostream_formatter {};
0055
0056 // Ensure ActsPlugins namespace is used when present
0057 #if __has_include(<ActsPlugins/DD4hep/ConvertDD4hepDetector.hpp>)
0058 // Acts_MAJOR_VERSION >= 44
0059 using DD4hepDetectorElement = ActsPlugins::DD4hepDetectorElement;
0060 using ActsPlugins::convertDD4hepDetector;
0061 using ActsPlugins::DD4hepFieldAdapter;
0062 using ActsPlugins::sortDetElementsByID;
0063 #else
0064 // Acts_MAJOR_VERSION < 44
0065 using DD4hepDetectorElement = Acts::DD4hepDetectorElement;
0066 using Acts::convertDD4hepDetector;
0067 using Acts::DD4hepFieldAdapter;
0068 using Acts::sortDetElementsByID;
0069 #endif
0070
0071 /// @brief Material decorator wrapper that tracks per-layer material assignment
0072 ///
0073 /// Wraps Acts::JsonMaterialDecorator and, for each decorate() call, records
0074 /// whether material was assigned to any approach surface in each
0075 /// (volume, layer) pair. After geometry conversion, call check() to emit
0076 /// critical log messages for every layer that was visited but never had
0077 /// material assigned to any of its approach surfaces.
0078 class EpicJsonMaterialDecorator : public Acts::IMaterialDecorator {
0079 public:
0080   /// Key identifying a layer: (volume id, layer id)
0081   using LayerKey = std::pair<Acts::GeometryIdentifier::Value, Acts::GeometryIdentifier::Value>;
0082
0083   EpicJsonMaterialDecorator(const Acts::MaterialMapJsonConverter::Config& rConfig,
0084                             const std::string& jFileName, Acts::Logging::Level level,
0085                             std::shared_ptr<spdlog::logger> logger)
0086       : m_inner(rConfig, jFileName, level), m_log(std::move(logger)) {}
0087
0088   void decorate(Acts::Surface& surface) const override {
0089     m_inner.decorate(surface);
0090     const auto id = surface.geometryId();
0091     m_log->trace("{} assigned to surface with geometryId=(volume={}, boundary={}, layer={}, "
0092                  "approach={}, sensitive={}, extra={})",
0093                  (surface.surfaceMaterial() != nullptr) ? "Material" : "No material", id.volume(),
0094                  id.boundary(), id.layer(), id.approach(), id.sensitive(), id.extra());
0095     // Only consider approach surfaces
0096     if (id.approach() == 0) {
0097       return;
0098     }
0099     LayerKey key{id.volume(), id.layer()};
0100     // Record that this layer was visited
0101     m_decoratedLayers.insert(key);
0102     // If material was assigned, record that
0103     if (surface.surfaceMaterial() != nullptr) {
0104       m_layersWithMaterial.insert(key);
0105     }
0106   }
0107
0108   void decorate(Acts::TrackingVolume& volume) const override { m_inner.decorate(volume); }
0109
0110   /// Report every decorated layer that never received material on any approach surface.
0111   void check() const {
0112     for (const auto& key : m_decoratedLayers) {
0113       if (m_layersWithMaterial.find(key) == m_layersWithMaterial.end()) {
0114         m_log->critical(
0115             "No material assigned to any approach surface in layer (volume={}, layer={})",
0116             key.first, key.second);
0117       }
0118     }
0119   }
0120
0121 private:
0122   Acts::JsonMaterialDecorator m_inner;
0123   std::shared_ptr<spdlog::logger> m_log;
0124   /// All (volume, layer) pairs seen during decoration
0125   mutable std::set<LayerKey> m_decoratedLayers;
0126   /// Subset of decorated layers that had at least one surface with material
0127   mutable std::set<LayerKey> m_layersWithMaterial;
0128 };
0129
0130 void ActsGeometryProvider::initialize(const dd4hep::Detector* dd4hep_geo, std::string material_file,
0131                                       std::shared_ptr<spdlog::logger> log,
0132                                       std::shared_ptr<spdlog::logger> init_log) {
0133   // LOGGING
0134   m_log      = log;
0135   m_init_log = init_log;
0136
0137   m_init_log->debug("ActsGeometryProvider initializing...");
0138
0139   m_init_log->debug("Set TGeoManager and acts_init_log_level log levels");
0140   // Turn off TGeo printouts if appropriate for the msg level
0141   if (m_log->level() >= (int)spdlog::level::info) {
0142     TGeoManager::SetVerboseLevel(0);
0143   }
0144
0145   // Set ACTS logging level
0146   auto acts_init_log_level = eicrecon::SpdlogToActsLevel(m_init_log->level());
0147
0148   m_dd4hepDetector = dd4hep_geo;
0149
0150   // Load ACTS materials maps
0151   std::shared_ptr<const Acts::IMaterialDecorator> materialDeco{nullptr};
0152   if (!material_file.empty()) {
0153     m_init_log->info("loading materials map from file: '{}'", material_file);
0154     // Set up the converter first
0155     Acts::MaterialMapJsonConverter::Config jsonGeoConvConfig;
0156     // Set up the json-based decorator, wrapped to report undecorated layers
0157     materialDeco = std::make_shared<const EpicJsonMaterialDecorator>(
0158         jsonGeoConvConfig, material_file, acts_init_log_level, m_init_log);
0159   }
0160
0161   // Geometry identifier hook to write detector ID to extra field
0162   class ConvertDD4hepDetectorGeometryIdentifierHook : public Acts::GeometryIdentifierHook {
0163     Acts::GeometryIdentifier decorateIdentifier(Acts::GeometryIdentifier identifier,
0164                                                 const Acts::Surface& surface) const override {
0165 #if Acts_VERSION_MAJOR >= 45
0166       const auto* placement          = surface.surfacePlacement();
0167       const auto* dd4hep_det_element = dynamic_cast<const DD4hepDetectorElement*>(placement);
0168 #else
0169       const auto* dd4hep_det_element =
0170           dynamic_cast<const DD4hepDetectorElement*>(surface.associatedDetectorElement());
0171 #endif
0172       if (dd4hep_det_element == nullptr) {
0173         return identifier;
0174       }
0175       // set 8-bit extra field to 8-bit DD4hep detector ID
0176 #if Acts_VERSION_MAJOR >= 40
0177       return identifier.withExtra(0xff & dd4hep_det_element->identifier());
0178 #else
0179       return identifier.setExtra(0xff & dd4hep_det_element->identifier());
0180 #endif
0181     };
0182   };
0183   auto geometryIdHook = std::make_shared<ConvertDD4hepDetectorGeometryIdentifierHook>();
0184
0185   // Convert DD4hep geometry to ACTS
0186   m_init_log->info("Converting DD4Hep geometry to ACTS...");
0187   auto logger                  = eicrecon::getSpdlogLogger("CONV", m_log);
0188   Acts::BinningType bTypePhi   = Acts::equidistant;
0189   Acts::BinningType bTypeR     = Acts::equidistant;
0190   Acts::BinningType bTypeZ     = Acts::equidistant;
0191   double layerEnvelopeR        = Acts::UnitConstants::mm;
0192   double layerEnvelopeZ        = Acts::UnitConstants::mm;
0193   double defaultLayerThickness = Acts::UnitConstants::fm;
0194
0195   try {
0196     m_trackingGeo =
0197         convertDD4hepDetector(m_dd4hepDetector->world(), *logger, bTypePhi, bTypeR, bTypeZ,
0198                               layerEnvelopeR, layerEnvelopeZ, defaultLayerThickness,
0199                               sortDetElementsByID, m_trackingGeoCtx, materialDeco, geometryIdHook);
0200   } catch (std::exception& ex) {
0201     m_init_log->error("Error during DD4Hep -> ACTS geometry conversion: {}", ex.what());
0202     m_init_log->info("Set parameter acts::InitLogLevel=trace to see conversion info and possibly "
0203                      "identify failing geometry");
0204     throw;
0205   }
0206
0207   m_init_log->info("DD4Hep geometry converted!");
0208
0209   // Report layers that were visited but never had material assigned
0210   if (auto epicDeco = std::dynamic_pointer_cast<const EpicJsonMaterialDecorator>(materialDeco)) {
0211     epicDeco->check();
0212   }
0213
0214   // Visit surfaces
0215   m_init_log->info("Checking surfaces...");
0216   if (m_trackingGeo) {
0217     // Write tracking geometry to collection of obj or ply files
0218     const Acts::TrackingVolume* world = m_trackingGeo->highestTrackingVolume();
0219     if (m_objWriteIt) {
0220       m_init_log->info("Writing obj files to {}...", m_outputDir);
0221       Acts::ObjVisualization3D objVis;
0222       Acts::GeometryView3D::drawTrackingVolume(objVis, *world, m_trackingGeoCtx, m_containerView,
0223                                                m_volumeView, m_passiveView, m_sensitiveView,
0224                                                m_gridView, m_objWriteIt, m_outputTag, m_outputDir);
0225     }
0226     if (m_plyWriteIt) {
0227       m_init_log->info("Writing ply files to {}...", m_outputDir);
0228       Acts::PlyVisualization3D plyVis;
0229       Acts::GeometryView3D::drawTrackingVolume(plyVis, *world, m_trackingGeoCtx, m_containerView,
0230                                                m_volumeView, m_passiveView, m_sensitiveView,
0231                                                m_gridView, m_plyWriteIt, m_outputTag, m_outputDir);
0232     }
0233
0234     m_init_log->debug("visiting all the surfaces  ");
0235     m_trackingGeo->visitSurfaces([this](const Acts::Surface* surface) {
0236       // for now we just require a valid surface
0237       if (surface == nullptr) {
0238         m_init_log->info("no surface??? ");
0239         return;
0240       }
0241 #if Acts_VERSION_MAJOR >= 45
0242       const auto* placement   = surface->surfacePlacement();
0243       const auto* det_element = dynamic_cast<const DD4hepDetectorElement*>(placement);
0244 #else
0245       const auto* det_element =
0246           dynamic_cast<const DD4hepDetectorElement*>(surface->associatedDetectorElement());
0247 #endif
0248
0249       if (det_element == nullptr) {
0250         m_init_log->error("invalid det_element!!! det_element == nullptr ");
0251         return;
0252       }
0253
0254       // more verbose output is lower enum value
0255       m_init_log->debug(" det_element->identifier() = {} ", det_element->identifier());
0256       auto volman   = m_dd4hepDetector->volumeManager();
0257       auto* vol_ctx = volman.lookupContext(det_element->identifier());
0258       auto vol_id   = vol_ctx->identifier;
0259
0260       if (m_init_log->level() <= spdlog::level::debug) {
0261         auto de = vol_ctx->element;
0262         m_init_log->debug("  de.path          = {}", de.path());
0263         m_init_log->debug("  de.placementPath = {}", de.placementPath());
0264       }
0265
0266       this->m_surfaces.insert_or_assign(vol_id, surface);
0267     });
0268   } else {
0269     m_init_log->error("m_trackingGeo==null why am I still alive???");
0270   }
0271
0272   // Load ACTS magnetic field
0273   m_init_log->info("Loading magnetic field...");
0274   m_magneticField = std::make_shared<DD4hepFieldAdapter>(m_dd4hepDetector->field());
0275   auto bCache     = m_magneticField->makeCache(Acts::MagneticFieldContext{});
0276   for (int z : {0, 500, 1000, 1500, 2000, 3000, 4000}) {
0277     auto b = m_magneticField->getField({0.0, 0.0, double(z)}, bCache).value();
0278     m_init_log->debug("B(z = {:>5} [mm]) = {} T", z, b.transpose() / Acts::UnitConstants::T);
0279   }
0280
0281   m_init_log->info("ActsGeometryProvider initialization complete");
0282 }
0283
0284 std::shared_ptr<const Acts::MagneticFieldProvider> ActsGeometryProvider::getFieldProvider() const {
0285   return m_magneticField;
0286 }