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