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