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 void draw_surfaces(std::shared_ptr<const Acts::TrackingGeometry> trk_geo, const Acts::GeometryContext geo_ctx, const std::string& fname)
0030 {
0031   using namespace Acts;
0032   std::vector<const Surface*> surfaces;
0033
0034   trk_geo->visitSurfaces([&](const Acts::Surface* surface) {
0035     // for now we just require a valid surface
0036     if (surface == nullptr) {
0037       std::cout << " Not a surface \n";
0038       return;
0039     }
0040     surfaces.push_back(surface);
0041   });
0042   std::ofstream os;
0043   os.open(fname);
0044   os << std::fixed << std::setprecision(6);
0045   size_t nVtx = 0;
0046   for (const auto& srfx : surfaces) {
0047     const auto* srf    = dynamic_cast<const PlaneSurface*>(srfx);
0048     const auto* bounds = dynamic_cast<const PlanarBounds*>(&srf->bounds());
0049     for (const auto& vtxloc : bounds->vertices()) {
0050       Vector3 vtx = srf->transform(geo_ctx) * Vector3(vtxloc.x(), vtxloc.y(), 0);
0051       os << "v " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
0052     }
0053     // connect them
0054     os << "f";
0055     for (size_t i = 1; i <= bounds->vertices().size(); ++i) {
0056       os << " " << nVtx + i;
0057     }
0058     os << "\n";
0059     nVtx += bounds->vertices().size();
0060   }
0061   os.close();
0062 }
0063
0064 using namespace Gaudi;
0065
0066 namespace Jug::Reco {
0067
0068 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
0069 DECLARE_COMPONENT(ActsGeoSvc)
0070
0071 ActsGeoSvc::ActsGeoSvc(const std::string& name, ISvcLocator* svc)
0072     : base_class(name, svc)
0073     , m_log(msgSvc(), name) {}
0074
0075 ActsGeoSvc::~ActsGeoSvc() {
0076   if (m_dd4hepGeo != nullptr) {
0077     try {
0078       m_dd4hepGeo->destroyInstance();
0079       m_dd4hepGeo = nullptr;
0080     } catch (...) {
0081     }
0082   }
0083 }
0084
0085 StatusCode ActsGeoSvc::initialize() {
0086   StatusCode sc = Service::initialize();
0087   if (!sc.isSuccess()) {
0088     return sc;
0089   }
0090   // Turn off TGeo printouts if appropriate for the msg level
0091   if (msgLevel() >= MSG::INFO) {
0092     TGeoManager::SetVerboseLevel(0);
0093   }
0094   uint printoutLevel = msgLevel();
0095   dd4hep::setPrintLevel(dd4hep::PrintLevel(printoutLevel));
0096   // m_incidentSvc->addListener(this, "GeometryFailure");
0097   if (buildDD4HepGeo().isFailure()) {
0098     m_log << MSG::ERROR << "Could not build DD4Hep geometry" << endmsg;
0099   } else {
0100     m_log << MSG::INFO << "DD4Hep geometry SUCCESSFULLY built" << endmsg;
0101   }
0102
0103   // Set ACTS logging level
0104   auto im = s_msgMap.find(msgLevel());
0105   if (im != s_msgMap.end()) {
0106     m_actsLoggingLevel = im->second;
0107   }
0108
0109   // Load ACTS materials maps
0110   if (m_jsonFileName.size() > 0) {
0111     m_log << MSG::INFO << "loading materials map from file:  '" << m_jsonFileName << "'" << endmsg;
0112     // Set up the converter first
0113     Acts::MaterialMapJsonConverter::Config jsonGeoConvConfig;
0114     // Set up the json-based decorator
0115     m_materialDeco = std::make_shared<const Acts::JsonMaterialDecorator>(
0116       jsonGeoConvConfig, m_jsonFileName, m_actsLoggingLevel);
0117   }
0118
0119   // Convert DD4hep geometry to ACTS
0120   auto logger = Acts::getDefaultLogger("CONV", m_actsLoggingLevel);
0121   Acts::BinningType bTypePhi = Acts::equidistant;
0122   Acts::BinningType bTypeR = Acts::equidistant;
0123   Acts::BinningType bTypeZ = Acts::equidistant;
0124   double layerEnvelopeR = Acts::UnitConstants::mm;
0125   double layerEnvelopeZ = Acts::UnitConstants::mm;
0126   double defaultLayerThickness = Acts::UnitConstants::fm;
0127   using Acts::sortDetElementsByID;
0128   m_trackingGeo = Acts::convertDD4hepDetector(
0129       m_dd4hepGeo->world(),
0130       *logger,
0131       bTypePhi,
0132       bTypeR,
0133       bTypeZ,
0134       layerEnvelopeR,
0135       layerEnvelopeZ,
0136       defaultLayerThickness,
0137       sortDetElementsByID,
0138       m_trackingGeoCtx,
0139       m_materialDeco);
0140   // Visit surfaces
0141   if (m_trackingGeo) {
0142     draw_surfaces(m_trackingGeo, m_trackingGeoCtx, "tracking_geometry.obj");
0143     debug() << "visiting all the surfaces  " << endmsg;
0144     m_trackingGeo->visitSurfaces([this](const Acts::Surface* surface) {
0145       // for now we just require a valid surface
0146       if (surface == nullptr) {
0147         info() << "no surface??? " << endmsg;
0148         return;
0149       }
0150       const auto* det_element =
0151         dynamic_cast<const Acts::DD4hepDetectorElement*>(surface->associatedDetectorElement());
0152
0153       if (det_element == nullptr) {
0154         error() << "invalid det_element!!! " << endmsg;
0155         return;
0156       }
0157       // more verbose output is lower enum value
0158       debug() << " det_element->identifier() " << det_element->identifier() << endmsg;
0159       auto volman  = m_dd4hepGeo->volumeManager();
0160       auto* vol_ctx = volman.lookupContext(det_element->identifier());
0161       auto vol_id  = vol_ctx->identifier;
0162
0163       if (msgLevel() <= MSG::DEBUG) {
0164         auto de  = vol_ctx->element;
0165         debug() << de.path() << endmsg;
0166         debug() << de.placementPath() << endmsg;
0167       }
0168       this->m_surfaces.insert_or_assign(vol_id, surface);
0169     });
0170   }
0171
0172   // Load ACTS magnetic field
0173   m_magneticField = std::make_shared<const Jug::BField::DD4hepBField>(m_dd4hepGeo);
0174   Acts::MagneticFieldContext m_fieldctx{Jug::BField::BFieldVariant(m_magneticField)};
0175   auto bCache = m_magneticField->makeCache(m_fieldctx);
0176   for (int z : {0, 1000, 2000, 4000}) {
0177     auto b = m_magneticField->getField({0.0, 0.0, double(z)}, bCache).value();
0178     debug() << "B(z=" << z << " mm) = " << b.transpose()  << " T" << endmsg;
0179   }
0180
0181   return StatusCode::SUCCESS;
0182 }
0183
0184 StatusCode ActsGeoSvc::finalize() { return StatusCode::SUCCESS; }
0185
0186 StatusCode ActsGeoSvc::buildDD4HepGeo() {
0187   // we retrieve the the static instance of the DD4HEP::Geometry
0188   m_dd4hepGeo = &(dd4hep::Detector::getInstance());
0189   m_dd4hepGeo->addExtension<IActsGeoSvc>(this);
0190
0191   // load geometry
0192   for (auto& filename : m_xmlFileNames) {
0193     m_log << MSG::INFO << "loading geometry from file:  '" << filename << "'" << endmsg;
0194     m_dd4hepGeo->fromCompact(filename);
0195   }
0196   m_dd4hepGeo->volumeManager();
0197   m_dd4hepGeo->apply("DD4hepVolumeManager", 0, nullptr);
0198   return StatusCode::SUCCESS;
0199 }
0200
0201 }