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