src/components/ActsGeoSvc.cpp

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 Whitney Armstrong, Wouter Deconinck
0003
0004 #include "ActsGeoSvc.h"
0005 #include "GaudiKernel/Service.h"
0006 #include "TGeoManager.h"
0007
0008 #include "DD4hep/Printout.h"
0009
0010 #include "JugTrack/ACTSLogger.h"
0011
0012 #include "Acts/Geometry/TrackingGeometry.hpp"
0013 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0014 #include "Acts/Surfaces/PlaneSurface.hpp"
0015
0016 #if Acts_VERSION_MAJOR >= 44
0017 #include "ActsPlugins/DD4hep/ConvertDD4hepDetector.hpp"
0018 #include "ActsPlugins/Json/JsonMaterialDecorator.hpp"
0019 #include "ActsPlugins/Json/MaterialMapJsonConverter.hpp"
0020 using ActsPlugins::sortDetElementsByID;
0021 using ActsPlugins::convertDD4hepDetector;
0022 using ActsPlugins::DD4hepDetectorElement;
0023 #else
0024 #include "Acts/Plugins/DD4hep/ConvertDD4hepDetector.hpp"
0025 #include "Acts/Plugins/Json/JsonMaterialDecorator.hpp"
0026 #include "Acts/Plugins/Json/MaterialMapJsonConverter.hpp"
0027 using Acts::sortDetElementsByID;
0028 using Acts::convertDD4hepDetector;
0029 using Acts::DD4hepDetectorElement;
0030 #endif
0031
0032 static const std::map<int, Acts::Logging::Level> s_msgMap = {
0033     {MSG::DEBUG, Acts::Logging::DEBUG},
0034     {MSG::VERBOSE, Acts::Logging::VERBOSE},
0035     {MSG::INFO, Acts::Logging::INFO},
0036     {MSG::WARNING, Acts::Logging::WARNING},
0037     {MSG::ERROR, Acts::Logging::ERROR},
0038     {MSG::FATAL, Acts::Logging::FATAL},
0039     {MSG::ALWAYS, Acts::Logging::MAX},
0040 };
0041
0042 using namespace Gaudi;
0043
0044 namespace Jug::Reco {
0045
0046 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0047 DECLARE_COMPONENT(ActsGeoSvc)
0048
0049 ActsGeoSvc::ActsGeoSvc(const std::string& name, ISvcLocator* svc)
0050     : base_class(name, svc)
0051     , m_log(msgSvc(), name) {}
0052
0053 ActsGeoSvc::~ActsGeoSvc() {
0054   if (m_dd4hepGeo != nullptr) {
0055     try {
0056       m_dd4hepGeo->destroyInstance();
0057       m_dd4hepGeo = nullptr;
0058     } catch (...) {
0059     }
0060   }
0061 }
0062
0063 StatusCode ActsGeoSvc::initialize() {
0064   StatusCode sc = Service::initialize();
0065   if (!sc.isSuccess()) {
0066     return sc;
0067   }
0068   // Turn off TGeo printouts if appropriate for the msg level
0069   if (msgLevel() >= MSG::INFO) {
0070     TGeoManager::SetVerboseLevel(0);
0071   }
0072   uint printoutLevel = msgLevel();
0073   dd4hep::setPrintLevel(dd4hep::PrintLevel(printoutLevel));
0074   // m_incidentSvc->addListener(this, "GeometryFailure");
0075   if (buildDD4HepGeo().isFailure()) {
0076     m_log << MSG::ERROR << "Could not build DD4Hep geometry" << endmsg;
0077   } else {
0078     m_log << MSG::INFO << "DD4Hep geometry SUCCESSFULLY built" << endmsg;
0079   }
0080
0081   // Set ACTS logging level
0082   auto im = s_msgMap.find(msgLevel());
0083   if (im != s_msgMap.end()) {
0084     m_actsLoggingLevel = im->second;
0085   }
0086
0087   // Load ACTS materials maps
0088   if (m_jsonFileName.size() > 0) {
0089     m_log << MSG::INFO << "loading materials map from file:  '" << m_jsonFileName << "'" << endmsg;
0090     // Set up the converter first
0091     Acts::MaterialMapJsonConverter::Config jsonGeoConvConfig;
0092     // Set up the json-based decorator
0093     m_materialDeco = std::make_shared<const Acts::JsonMaterialDecorator>(
0094       jsonGeoConvConfig, m_jsonFileName, m_actsLoggingLevel);
0095   }
0096
0097   // Convert DD4hep geometry to ACTS
0098   auto logger = Acts::getDefaultLogger("CONV", m_actsLoggingLevel);
0099   Acts::BinningType bTypePhi = Acts::equidistant;
0100   Acts::BinningType bTypeR = Acts::equidistant;
0101   Acts::BinningType bTypeZ = Acts::equidistant;
0102   double layerEnvelopeR = Acts::UnitConstants::mm;
0103   double layerEnvelopeZ = Acts::UnitConstants::mm;
0104   double defaultLayerThickness = Acts::UnitConstants::fm;
0105   m_trackingGeo = convertDD4hepDetector(
0106       m_dd4hepGeo->world(),
0107       *logger,
0108       bTypePhi,
0109       bTypeR,
0110       bTypeZ,
0111       layerEnvelopeR,
0112       layerEnvelopeZ,
0113       defaultLayerThickness,
0114       sortDetElementsByID,
0115       m_trackingGeoCtx,
0116       m_materialDeco);
0117   // Visit surfaces
0118   if (m_trackingGeo) {
0119     debug() << "visiting all the surfaces  " << endmsg;
0120     m_trackingGeo->visitSurfaces([this](const Acts::Surface* surface) {
0121       // for now we just require a valid surface
0122       if (surface == nullptr) {
0123         info() << "no surface??? " << endmsg;
0124         return;
0125       }
0126 #if Acts_VERSION_MAJOR >= 45
0127       const auto* placement   = surface->surfacePlacement();
0128       const auto* det_element = dynamic_cast<const DD4hepDetectorElement*>(placement);
0129 #else
0130       const auto* det_element =
0131         dynamic_cast<const DD4hepDetectorElement*>(surface->associatedDetectorElement());
0132 #endif
0133
0134       if (det_element == nullptr) {
0135         error() << "invalid det_element!!! " << endmsg;
0136         return;
0137       }
0138       // more verbose output is lower enum value
0139       debug() << " det_element->identifier() " << det_element->identifier() << endmsg;
0140       auto volman  = m_dd4hepGeo->volumeManager();
0141       auto* vol_ctx = volman.lookupContext(det_element->identifier());
0142       auto vol_id  = vol_ctx->identifier;
0143
0144       if (msgLevel() <= MSG::DEBUG) {
0145         auto de  = vol_ctx->element;
0146         debug() << de.path() << endmsg;
0147         debug() << de.placementPath() << endmsg;
0148       }
0149       this->m_surfaces.insert_or_assign(vol_id, surface);
0150     });
0151   }
0152
0153   // Load ACTS magnetic field
0154   m_magneticField = std::make_shared<const Jug::BField::DD4hepBField>(m_dd4hepGeo);
0155   Acts::MagneticFieldContext m_fieldctx{Jug::BField::BFieldVariant(m_magneticField)};
0156   auto bCache = m_magneticField->makeCache(m_fieldctx);
0157   for (int z : {0, 1000, 2000, 4000}) {
0158     auto b = m_magneticField->getField({0.0, 0.0, double(z)}, bCache).value();
0159     debug() << "B(z=" << z << " mm) = " << b.transpose()  << " T" << endmsg;
0160   }
0161
0162   return StatusCode::SUCCESS;
0163 }
0164
0165 StatusCode ActsGeoSvc::finalize() { return StatusCode::SUCCESS; }
0166
0167 StatusCode ActsGeoSvc::buildDD4HepGeo() {
0168   // we retrieve the the static instance of the DD4HEP::Geometry
0169   m_dd4hepGeo = &(dd4hep::Detector::getInstance());
0170   m_dd4hepGeo->addExtension<IActsGeoSvc>(this);
0171
0172   // load geometry
0173   for (auto& filename : m_xmlFileNames) {
0174     m_log << MSG::INFO << "loading geometry from file:  '" << filename << "'" << endmsg;
0175     m_dd4hepGeo->fromCompact(filename);
0176   }
0177   m_dd4hepGeo->volumeManager();
0178   m_dd4hepGeo->apply("DD4hepVolumeManager", 0, nullptr);
0179   return StatusCode::SUCCESS;
0180 }
0181
0182 }