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 #include <Acts/Geometry/DetectorElementBase.hpp>
0005 #include <Acts/Geometry/GeometryIdentifier.hpp>
0006 #include <Acts/Geometry/TrackingGeometry.hpp>
0007 #include <Acts/Geometry/TrackingVolume.hpp>
0008 #include <Acts/MagneticField/MagneticFieldContext.hpp>
0009 #include <Acts/Material/IMaterialDecorator.hpp>
0010 #if __has_include(<ActsPlugins/DD4hep/ConvertDD4hepDetector.hpp>)
0011 #include <ActsPlugins/DD4hep/ConvertDD4hepDetector.hpp>
0012 #include <ActsPlugins/DD4hep/DD4hepDetectorElement.hpp>
0013 #include <ActsPlugins/Json/JsonMaterialDecorator.hpp>
0014 #include <ActsPlugins/Json/MaterialMapJsonConverter.hpp>
0015 #else
0016 #include <Acts/Plugins/DD4hep/ConvertDD4hepDetector.hpp>
0017 #include <Acts/Plugins/DD4hep/DD4hepDetectorElement.hpp>
0018 #include <Acts/Plugins/Json/JsonMaterialDecorator.hpp>
0019 #include <Acts/Plugins/Json/MaterialMapJsonConverter.hpp>
0020 #endif
0021 #include <Acts/Surfaces/Surface.hpp>
0022 #include <Acts/Utilities/BinningType.hpp>
0023 #include <Acts/Utilities/Result.hpp>
0024 #include <Acts/Visualization/GeometryView3D.hpp>
0025 #include <Acts/Visualization/ObjVisualization3D.hpp>
0026 #include <Acts/Visualization/PlyVisualization3D.hpp>
0027 #include <DD4hep/DetElement.h>
0028 #include <DD4hep/VolumeManager.h>
0029 #include <JANA/JException.h>
0030 #include <TGeoManager.h>
0031 #include <fmt/core.h>
0032 #include <fmt/format.h>
0033 #include <fmt/ostream.h>
0034 #include <spdlog/common.h>
0035 #include <exception>
0036 #include <filesystem>
0037 #include <functional>
0038 #include <initializer_list>
0039 #include <type_traits>
0040
0041 #include "ActsGeometryProvider.h"
0042 #include "DD4hepBField.h"
0043 #include "extensions/spdlog/SpdlogToActs.h"
0044
0045 // Formatter for Eigen matrices
0046 #if FMT_VERSION >= 90000
0047 #include <Eigen/Core>
0048
0049 template <typename T>
0050 struct fmt::formatter<T, std::enable_if_t<std::is_base_of_v<Eigen::MatrixBase<T>, T>, char>>
0051     : fmt::ostream_formatter {};
0052 #endif // FMT_VERSION >= 90000
0053
0054 // Ensure ActsPlugins namespace is used when present
0055 #if __has_include(<ActsPlugins/DD4hep/ConvertDD4hepDetector.hpp>)
0056 // Acts_MAJOR_VERSION >= 44
0057 using DD4hepDetectorElement = ActsPlugins::DD4hepDetectorElement;
0058 using ActsPlugins::convertDD4hepDetector;
0059 using ActsPlugins::sortDetElementsByID;
0060 #else
0061 // Acts_MAJOR_VERSION < 44
0062 using DD4hepDetectorElement = Acts::DD4hepDetectorElement;
0063 using Acts::convertDD4hepDetector;
0064 using Acts::sortDetElementsByID;
0065 #endif
0066
0067 void ActsGeometryProvider::initialize(const dd4hep::Detector* dd4hep_geo, std::string material_file,
0068                                       std::shared_ptr<spdlog::logger> log,
0069                                       std::shared_ptr<spdlog::logger> init_log) {
0070   // LOGGING
0071   m_log      = log;
0072   m_init_log = init_log;
0073
0074   m_init_log->debug("ActsGeometryProvider initializing...");
0075
0076   m_init_log->debug("Set TGeoManager and acts_init_log_level log levels");
0077   // Turn off TGeo printouts if appropriate for the msg level
0078   if (m_log->level() >= (int)spdlog::level::info) {
0079     TGeoManager::SetVerboseLevel(0);
0080   }
0081
0082   // Set ACTS logging level
0083   auto acts_init_log_level = eicrecon::SpdlogToActsLevel(m_init_log->level());
0084
0085   m_dd4hepDetector = dd4hep_geo;
0086
0087   // Load ACTS materials maps
0088   std::shared_ptr<const Acts::IMaterialDecorator> materialDeco{nullptr};
0089   if (!material_file.empty()) {
0090     m_init_log->info("loading materials map from file: '{}'", material_file);
0091     // Set up the converter first
0092     Acts::MaterialMapJsonConverter::Config jsonGeoConvConfig;
0093     // Set up the json-based decorator
0094     materialDeco = std::make_shared<const Acts::JsonMaterialDecorator>(
0095         jsonGeoConvConfig, material_file, acts_init_log_level);
0096   }
0097
0098   // Geometry identifier hook to write detector ID to extra field
0099   class ConvertDD4hepDetectorGeometryIdentifierHook : public Acts::GeometryIdentifierHook {
0100     Acts::GeometryIdentifier decorateIdentifier(Acts::GeometryIdentifier identifier,
0101                                                 const Acts::Surface& surface) const override {
0102       const auto* dd4hep_det_element =
0103           dynamic_cast<const DD4hepDetectorElement*>(surface.associatedDetectorElement());
0104       if (dd4hep_det_element == nullptr) {
0105         return identifier;
0106       }
0107       // set 8-bit extra field to 8-bit DD4hep detector ID
0108 #if Acts_VERSION_MAJOR >= 40
0109       return identifier.withExtra(0xff & dd4hep_det_element->identifier());
0110 #else
0111       return identifier.setExtra(0xff & dd4hep_det_element->identifier());
0112 #endif
0113     };
0114   };
0115   auto geometryIdHook = std::make_shared<ConvertDD4hepDetectorGeometryIdentifierHook>();
0116
0117   // Convert DD4hep geometry to ACTS
0118   m_init_log->info("Converting DD4Hep geometry to ACTS...");
0119   auto logger                  = eicrecon::getSpdlogLogger("CONV", m_log);
0120   Acts::BinningType bTypePhi   = Acts::equidistant;
0121   Acts::BinningType bTypeR     = Acts::equidistant;
0122   Acts::BinningType bTypeZ     = Acts::equidistant;
0123   double layerEnvelopeR        = Acts::UnitConstants::mm;
0124   double layerEnvelopeZ        = Acts::UnitConstants::mm;
0125   double defaultLayerThickness = Acts::UnitConstants::fm;
0126
0127   try {
0128     m_trackingGeo =
0129         convertDD4hepDetector(m_dd4hepDetector->world(), *logger, bTypePhi, bTypeR, bTypeZ,
0130                               layerEnvelopeR, layerEnvelopeZ, defaultLayerThickness,
0131                               sortDetElementsByID, m_trackingGeoCtx, materialDeco, geometryIdHook);
0132   } catch (std::exception& ex) {
0133     m_init_log->error("Error during DD4Hep -> ACTS geometry conversion: {}", ex.what());
0134     m_init_log->info("Set parameter acts::InitLogLevel=trace to see conversion info and possibly "
0135                      "identify failing geometry");
0136     throw JException(ex.what());
0137   }
0138
0139   m_init_log->info("DD4Hep geometry converted!");
0140
0141   // Visit surfaces
0142   m_init_log->info("Checking surfaces...");
0143   if (m_trackingGeo) {
0144     // Write tracking geometry to collection of obj or ply files
0145     const Acts::TrackingVolume* world = m_trackingGeo->highestTrackingVolume();
0146     if (m_objWriteIt) {
0147       m_init_log->info("Writing obj files to {}...", m_outputDir);
0148       Acts::ObjVisualization3D objVis;
0149       Acts::GeometryView3D::drawTrackingVolume(objVis, *world, m_trackingGeoCtx, m_containerView,
0150                                                m_volumeView, m_passiveView, m_sensitiveView,
0151                                                m_gridView, m_objWriteIt, m_outputTag, m_outputDir);
0152     }
0153     if (m_plyWriteIt) {
0154       m_init_log->info("Writing ply files to {}...", m_outputDir);
0155       Acts::PlyVisualization3D plyVis;
0156       Acts::GeometryView3D::drawTrackingVolume(plyVis, *world, m_trackingGeoCtx, m_containerView,
0157                                                m_volumeView, m_passiveView, m_sensitiveView,
0158                                                m_gridView, m_plyWriteIt, m_outputTag, m_outputDir);
0159     }
0160
0161     m_init_log->debug("visiting all the surfaces  ");
0162     m_trackingGeo->visitSurfaces([this](const Acts::Surface* surface) {
0163       // for now we just require a valid surface
0164       if (surface == nullptr) {
0165         m_init_log->info("no surface??? ");
0166         return;
0167       }
0168       const auto* det_element =
0169           dynamic_cast<const DD4hepDetectorElement*>(surface->associatedDetectorElement());
0170
0171       if (det_element == nullptr) {
0172         m_init_log->error("invalid det_element!!! det_element == nullptr ");
0173         return;
0174       }
0175
0176       // more verbose output is lower enum value
0177       m_init_log->debug(" det_element->identifier() = {} ", det_element->identifier());
0178       auto volman   = m_dd4hepDetector->volumeManager();
0179       auto* vol_ctx = volman.lookupContext(det_element->identifier());
0180       auto vol_id   = vol_ctx->identifier;
0181
0182       if (m_init_log->level() <= spdlog::level::debug) {
0183         auto de = vol_ctx->element;
0184         m_init_log->debug("  de.path          = {}", de.path());
0185         m_init_log->debug("  de.placementPath = {}", de.placementPath());
0186       }
0187
0188       this->m_surfaces.insert_or_assign(vol_id, surface);
0189     });
0190   } else {
0191     m_init_log->error("m_trackingGeo==null why am I still alive???");
0192   }
0193
0194   // Load ACTS magnetic field
0195   m_init_log->info("Loading magnetic field...");
0196   m_magneticField = std::make_shared<const eicrecon::BField::DD4hepBField>(m_dd4hepDetector);
0197   Acts::MagneticFieldContext m_fieldctx{eicrecon::BField::BFieldVariant(m_magneticField)};
0198   auto bCache = m_magneticField->makeCache(m_fieldctx);
0199   for (int z : {0, 500, 1000, 1500, 2000, 3000, 4000}) {
0200     auto b = m_magneticField->getField({0.0, 0.0, double(z)}, bCache).value();
0201     m_init_log->debug("B(z = {:>5} [mm]) = {} T", z, b.transpose() / Acts::UnitConstants::T);
0202   }
0203
0204   m_init_log->info("ActsGeometryProvider initialization complete");
0205 }
0206
0207 std::shared_ptr<const Acts::MagneticFieldProvider> ActsGeometryProvider::getFieldProvider() const {
0208   return m_magneticField;
0209 }