DD4hep/src/DD4hepLayerBuilder.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsPlugins/DD4hep/DD4hepLayerBuilder.hpp"
0010
0011 #include "Acts/Definitions/Common.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Geometry/ApproachDescriptor.hpp"
0014 #include "Acts/Geometry/CylinderLayer.hpp"
0015 #include "Acts/Geometry/DiscLayer.hpp"
0016 #include "Acts/Geometry/Extent.hpp"
0017 #include "Acts/Geometry/Layer.hpp"
0018 #include "Acts/Geometry/LayerCreator.hpp"
0019 #include "Acts/Geometry/ProtoLayer.hpp"
0020 #include "Acts/Surfaces/CylinderBounds.hpp"
0021 #include "Acts/Surfaces/RadialBounds.hpp"
0022 #include "Acts/Surfaces/Surface.hpp"
0023 #include "Acts/Surfaces/SurfaceArray.hpp"
0024 #include "Acts/Utilities/Logger.hpp"
0025 #include "ActsPlugins/DD4hep/DD4hepConversionHelpers.hpp"
0026 #include "ActsPlugins/DD4hep/DD4hepMaterialHelpers.hpp"
0027 #include "ActsPlugins/Root/TGeoPrimitivesHelper.hpp"
0028
0029 #include <algorithm>
0030 #include <cmath>
0031 #include <cstddef>
0032 #include <iterator>
0033 #include <ostream>
0034 #include <stdexcept>
0035 #include <utility>
0036
0037 #include <DD4hep/Alignments.h>
0038 #include <DD4hep/DetElement.h>
0039 #include <DD4hep/Volumes.h>
0040 #include <DDRec/DetectorData.h>
0041 #include <RtypesCore.h>
0042 #include <TGeoManager.h>
0043 #include <TGeoMatrix.h>
0044 #include <boost/algorithm/string.hpp>
0045
0046 using namespace Acts;
0047
0048 namespace ActsPlugins {
0049
0050 DD4hepLayerBuilder::DD4hepLayerBuilder(const DD4hepLayerBuilder::Config& config,
0051                                        std::unique_ptr<const Logger> logger)
0052     : m_cfg(), m_logger(std::move(logger)) {
0053   setConfiguration(config);
0054 }
0055
0056 DD4hepLayerBuilder::~DD4hepLayerBuilder() = default;
0057
0058 void DD4hepLayerBuilder::setConfiguration(
0059     const DD4hepLayerBuilder::Config& config) {
0060   m_cfg = config;
0061 }
0062
0063 const LayerVector DD4hepLayerBuilder::endcapLayers(
0064     const GeometryContext& gctx,
0065     const std::vector<dd4hep::DetElement>& dendcapLayers,
0066     const std::string& side) const {
0067   LayerVector layers;
0068   if (dendcapLayers.empty()) {
0069     ACTS_VERBOSE(" No layers handed over for " << side << " volume!");
0070   } else {
0071     ACTS_VERBOSE(" Received layers for " << side
0072                                          << " volume -> creating "
0073                                             "disc layers");
0074     // go through layers
0075     for (auto& detElement : dendcapLayers) {
0076       ACTS_VERBOSE("=> Translating layer from: " << detElement.name());
0077       // prepare the layer surfaces
0078       std::vector<std::shared_ptr<const Surface>> layerSurfaces;
0079       // access the extension of the layer
0080       // at this stage all layer detElements have extension (checked in
0081       // ConvertDD4hepDetector)
0082       auto& params = getParams(detElement);
0083       // collect the sensitive detector elements possibly contained by the layer
0084       resolveSensitive(detElement, layerSurfaces);
0085       // access the global transformation matrix of the layer
0086       auto transform =
0087           convertTransform(&(detElement.nominal().worldTransformation()));
0088       auto itransform = transform.inverse();
0089       // get the shape of the layer
0090       TGeoShape* geoShape =
0091           detElement.placement().ptr()->GetVolume()->GetShape();
0092       // create the proto layer
0093       ProtoLayer pl(gctx, layerSurfaces, itransform);
0094       if (logger().doPrint(Logging::VERBOSE)) {
0095         std::stringstream ss;
0096         pl.toStream(ss);
0097         ACTS_VERBOSE("Extent from surfaces: " << ss.str());
0098
0099         std::vector<double> rvalues;
0100         std::transform(layerSurfaces.begin(), layerSurfaces.end(),
0101                        std::back_inserter(rvalues), [&](const auto& surface) {
0102                          return VectorHelpers::perp(surface->center(gctx));
0103                        });
0104         std::ranges::sort(rvalues);
0105         std::vector<std::string> locs;
0106         std::transform(rvalues.begin(),
0107                        std::unique(rvalues.begin(), rvalues.end()),
0108                        std::back_inserter(locs),
0109                        [](const auto& v) { return std::to_string(v); });
0110         ACTS_VERBOSE(
0111             "-> unique r locations: " << boost::algorithm::join(locs, ", "));
0112       }
0113
0114       if (params.contains("envelope_r_min") &&
0115           params.contains("envelope_r_max") &&
0116           params.contains("envelope_z_min") &&
0117           params.contains("envelope_z_max")) {
0118         // set the values of the proto layer in case enevelopes are handed
0119         // over
0120         pl.envelope[AxisDirection::AxisR] = {
0121             params.get<double>("envelope_r_min"),
0122             params.get<double>("envelope_r_max")};
0123         pl.envelope[AxisDirection::AxisZ] = {
0124             params.get<double>("envelope_z_min"),
0125             params.get<double>("envelope_z_max")};
0126       } else if (geoShape != nullptr) {
0127         TGeoTubeSeg* tube = dynamic_cast<TGeoTubeSeg*>(geoShape);
0128         if (tube == nullptr) {
0129           ACTS_ERROR(" Disc layer has wrong shape - needs to be TGeoTubeSeg!");
0130           throw std::logic_error{
0131               "Disc layer has wrong shape - needs to be TGeoTubeSeg!"};
0132         }
0133         // extract the boundaries
0134         double rMin = tube->GetRmin() * UnitConstants::cm;
0135         double rMax = tube->GetRmax() * UnitConstants::cm;
0136
0137         // For disc layers, since ProtoLayer uses local coordinates,
0138         // we can simply use ±dz directly in local coordinates
0139         double dz = tube->GetDz() * UnitConstants::cm;
0140         double zMin = -dz;
0141         double zMax = +dz;
0142
0143         // check if layer has surfaces
0144         if (layerSurfaces.empty()) {
0145           ACTS_VERBOSE(" Disc layer has no sensitive surfaces.");
0146           // in case no surfaces are handed over the layer thickness will be
0147           // set to a default value to allow attaching material layers
0148           double eiz = (transform.translation().z() != 0.)
0149                            ? -m_cfg.defaultThickness
0150                            : 0.;
0151           double eoz = (transform.translation().z() != 0.)
0152                            ? +m_cfg.defaultThickness
0153                            : 0.;
0154           pl.extent.range(AxisDirection::AxisZ).set(eiz, eoz);
0155           pl.extent.range(AxisDirection::AxisR).set(rMin, rMax);
0156           pl.envelope[AxisDirection::AxisR] = {0., 0.};
0157           pl.envelope[AxisDirection::AxisZ] = {0., 0.};
0158         } else {
0159           ACTS_VERBOSE(" Disc layer has " << layerSurfaces.size()
0160                                           << " sensitive surfaces.");
0161           // set the values of the proto layer in case dimensions are given by
0162           // geometry
0163           pl.envelope[AxisDirection::AxisZ] = {
0164               std::abs(zMin - pl.min(AxisDirection::AxisZ)),
0165               std::abs(zMax - pl.max(AxisDirection::AxisZ))};
0166           pl.envelope[AxisDirection::AxisR] = {
0167               std::abs(rMin - pl.min(AxisDirection::AxisR)),
0168               std::abs(rMax - pl.max(AxisDirection::AxisR))};
0169           pl.extent.range(AxisDirection::AxisR).set(rMin, rMax);
0170         }
0171       } else {
0172         throw std::logic_error(
0173             std::string("Layer DetElement: ") + detElement.name() +
0174             std::string(" has neither a shape nor tolerances for envelopes "
0175                         "added to its extension. Please check your detector "
0176                         "constructor!"));
0177       }
0178
0179       std::shared_ptr<Layer> endcapLayer = nullptr;
0180
0181       // Check if DD4hep pre-defines the surface binning
0182       bool hasSurfaceBinning =
0183           getParamOr<bool>("surface_binning", detElement, true);
0184       std::size_t nPhi = 1;
0185       std::size_t nR = 1;
0186       if (hasSurfaceBinning) {
0187         if (params.contains("surface_binning_n_phi")) {
0188           nPhi = static_cast<std::size_t>(
0189               params.get<int>("surface_binning_n_phi"));
0190         }
0191         if (params.contains("surface_binning_n_r")) {
0192           nR = static_cast<std::size_t>(params.get<int>("surface_binning_n_r"));
0193         }
0194         hasSurfaceBinning = nR * nPhi > 1;
0195       }
0196
0197       // In case the layer is sensitive
0198       if (detElement.volume().isSensitive()) {
0199         // Create the sensitive surface
0200         auto sensitiveSurf = createSensitiveSurface(detElement, true);
0201         // Create the surfaceArray
0202         auto sArray = std::make_unique<SurfaceArray>(sensitiveSurf);
0203
0204         // create the share disc bounds
0205         auto dBounds = std::make_shared<const RadialBounds>(
0206             pl.min(AxisDirection::AxisR), pl.max(AxisDirection::AxisR));
0207         double thickness = std::abs(pl.max(AxisDirection::AxisZ) -
0208                                     pl.min(AxisDirection::AxisZ));
0209         // Create the layer containing the sensitive surface
0210         endcapLayer = DiscLayer::create(transform, dBounds, std::move(sArray),
0211                                         thickness, nullptr, active);
0212
0213       } else if (hasSurfaceBinning) {
0214         // This method uses the binning from DD4hep/xml
0215         endcapLayer = m_cfg.layerCreator->discLayer(
0216             gctx, layerSurfaces, nR, nPhi, pl, transform, nullptr);
0217       } else {
0218         // This method determines the binning automatically
0219         endcapLayer = m_cfg.layerCreator->discLayer(
0220             gctx, layerSurfaces, m_cfg.bTypeR, m_cfg.bTypePhi, pl, transform,
0221             nullptr);
0222       }
0223       // Add the ProtoMaterial if present
0224       addDiscLayerProtoMaterial(detElement, *endcapLayer, logger());
0225       // push back created layer
0226       layers.push_back(endcapLayer);
0227     }
0228   }
0229   return layers;
0230 }
0231
0232 const LayerVector DD4hepLayerBuilder::negativeLayers(
0233     const GeometryContext& gctx) const {
0234   return endcapLayers(gctx, m_cfg.negativeLayers, "negative");
0235 }
0236
0237 const LayerVector DD4hepLayerBuilder::centralLayers(
0238     const GeometryContext& gctx) const {
0239   LayerVector layers;
0240   if (m_cfg.centralLayers.empty()) {
0241     ACTS_VERBOSE(" No layers handed over for central volume!");
0242   } else {
0243     ACTS_VERBOSE(
0244         " Received layers for central volume -> creating "
0245         "cylindrical layers");
0246     // go through layers
0247     for (auto& detElement : m_cfg.centralLayers) {
0248       ACTS_VERBOSE("=> Translating layer from: " << detElement.name());
0249       // prepare the layer surfaces
0250       std::vector<std::shared_ptr<const Surface>> layerSurfaces;
0251       // access the extension of the layer
0252       // at this stage all layer detElements have extension (checked in
0253       // ConvertDD4hepDetector)
0254       auto& params = getParams(detElement);
0255       // collect the sensitive detector elements possibly contained by the layer
0256       resolveSensitive(detElement, layerSurfaces);
0257       // access the global transformation matrix of the layer
0258       auto transform =
0259           convertTransform(&(detElement.nominal().worldTransformation()));
0260       // get the shape of the layer
0261       TGeoShape* geoShape =
0262           detElement.placement().ptr()->GetVolume()->GetShape();
0263       // create the proto layer
0264       ProtoLayer pl(gctx, layerSurfaces);
0265       if (logger().doPrint(Logging::VERBOSE)) {
0266         std::stringstream ss;
0267         pl.toStream(ss);
0268         ACTS_VERBOSE("Extent from surfaces: " << ss.str());
0269         std::vector<double> zvalues;
0270         std::transform(layerSurfaces.begin(), layerSurfaces.end(),
0271                        std::back_inserter(zvalues), [&](const auto& surface) {
0272                          return surface->center(gctx)[eZ];
0273                        });
0274         std::ranges::sort(zvalues);
0275         std::vector<std::string> locs;
0276         std::transform(zvalues.begin(),
0277                        std::unique(zvalues.begin(), zvalues.end()),
0278                        std::back_inserter(locs),
0279                        [](const auto& v) { return std::to_string(v); });
0280         ACTS_VERBOSE(
0281             "-> unique z locations: " << boost::algorithm::join(locs, ", "));
0282       }
0283
0284       if (params.contains("envelope_r_min") &&
0285           params.contains("envelope_r_max") &&
0286           params.contains("envelope_z_min") &&
0287           params.contains("envelope_z_max")) {
0288         // set the values of the proto layer in case enevelopes are handed over
0289         pl.envelope[AxisDirection::AxisR] = {
0290             params.get<double>("envelope_r_min"),
0291             params.get<double>("envelope_r_max")};
0292         pl.envelope[AxisDirection::AxisZ] = {
0293             params.get<double>("envelope_z_min"),
0294             params.get<double>("envelope_z_max")};
0295       } else if (geoShape != nullptr) {
0296         TGeoTubeSeg* tube = dynamic_cast<TGeoTubeSeg*>(geoShape);
0297         if (tube == nullptr) {
0298           ACTS_ERROR(
0299               " Cylinder layer has wrong shape - needs to be TGeoTubeSeg!");
0300           throw std::logic_error{
0301               " Cylinder layer has wrong shape - needs to be TGeoTubeSeg!"};
0302         }
0303
0304         // extract the boundaries
0305         double rMin = tube->GetRmin() * UnitConstants::cm;
0306         double rMax = tube->GetRmax() * UnitConstants::cm;
0307         double dz = tube->GetDz() * UnitConstants::cm;
0308         // check if layer has surfaces
0309         if (layerSurfaces.empty()) {
0310           // in case no surfaces are handed over the layer thickness will be set
0311           // to a default value to allow attaching material layers
0312           double r = (rMin + rMax) * 0.5;
0313           // create layer without surfaces
0314           // manually create a proto layer
0315           double eir = (r != 0.) ? r - m_cfg.defaultThickness : 0.;
0316           double eor = (r != 0.) ? r + m_cfg.defaultThickness : 0.;
0317           pl.extent.range(AxisDirection::AxisR).set(eir, eor);
0318           pl.extent.range(AxisDirection::AxisZ).set(-dz, dz);
0319           pl.envelope[AxisDirection::AxisR] = {0., 0.};
0320           pl.envelope[AxisDirection::AxisZ] = {0., 0.};
0321         } else {
0322           // set the values of the proto layer in case dimensions are given by
0323           // geometry
0324           pl.envelope[AxisDirection::AxisZ] = {
0325               std::abs(-dz - pl.min(AxisDirection::AxisZ)),
0326               std::abs(dz - pl.max(AxisDirection::AxisZ))};
0327           pl.envelope[AxisDirection::AxisR] = {
0328               std::abs(rMin - pl.min(AxisDirection::AxisR)),
0329               std::abs(rMax - pl.max(AxisDirection::AxisR))};
0330         }
0331       } else {
0332         throw std::logic_error(
0333             std::string("Layer DetElement: ") + detElement.name() +
0334             std::string(" has neither a shape nor tolerances for envelopes "
0335                         "added to it¥s extension. Please check your detector "
0336                         "constructor!"));
0337       }
0338
0339       double halfZ =
0340           (pl.min(AxisDirection::AxisZ) - pl.max(AxisDirection::AxisZ)) * 0.5;
0341
0342       std::shared_ptr<Layer> centralLayer = nullptr;
0343       // In case the layer is sensitive
0344       if (detElement.volume().isSensitive()) {
0345         // Create the sensitive surface
0346         auto sensitiveSurf = createSensitiveSurface(detElement);
0347         // Create the surfaceArray
0348         std::unique_ptr<SurfaceArray> sArray =
0349             std::make_unique<SurfaceArray>(sensitiveSurf);
0350
0351         // create the layer
0352         double layerR =
0353             (pl.min(AxisDirection::AxisR) + pl.max(AxisDirection::AxisR)) * 0.5;
0354         double thickness = std::abs(pl.max(AxisDirection::AxisR) -
0355                                     pl.min(AxisDirection::AxisR));
0356         auto cBounds = std::make_shared<CylinderBounds>(layerR, halfZ);
0357         // Create the layer containing the sensitive surface
0358         centralLayer = CylinderLayer::create(
0359             transform, cBounds, std::move(sArray), thickness, nullptr, active);
0360
0361       } else {
0362         centralLayer = m_cfg.layerCreator->cylinderLayer(
0363             gctx, layerSurfaces, m_cfg.bTypePhi, m_cfg.bTypeZ, pl, transform,
0364             nullptr);
0365       }
0366       // Add the ProtoMaterial if present
0367       addCylinderLayerProtoMaterial(detElement, *centralLayer, logger());
0368       // push back created layer
0369       layers.push_back(centralLayer);
0370     }
0371   }
0372   return layers;
0373 }
0374
0375 const LayerVector DD4hepLayerBuilder::positiveLayers(
0376     const GeometryContext& gctx) const {
0377   return endcapLayers(gctx, m_cfg.positiveLayers, "positive");
0378 }
0379
0380 void DD4hepLayerBuilder::resolveSensitive(
0381     const dd4hep::DetElement& detElement,
0382     std::vector<std::shared_ptr<const Surface>>& surfaces) const {
0383   const dd4hep::DetElement::Children& children = detElement.children();
0384   if (!children.empty()) {
0385     for (auto& child : children) {
0386       dd4hep::DetElement childDetElement = child.second;
0387       if (childDetElement.volume().isSensitive()) {
0388         // create the surface
0389         surfaces.push_back(createSensitiveSurface(childDetElement, false));
0390       }
0391       resolveSensitive(childDetElement, surfaces);
0392     }
0393   }
0394 }
0395
0396 std::shared_ptr<const Surface> DD4hepLayerBuilder::createSensitiveSurface(
0397     const dd4hep::DetElement& detElement, bool isDisc) const {
0398   std::string detAxis =
0399       getParamOr<std::string>("axis_definitions", detElement, "XYZ");
0400   // Create the corresponding detector element !- memory leak --!
0401   auto dd4hepDetElement = m_cfg.detectorElementFactory(
0402       detElement, detAxis, UnitConstants::cm, isDisc, nullptr);
0403
0404   detElement.addExtension<DD4hepDetectorElementExtension>(
0405       new dd4hep::rec::StructExtension(
0406           DD4hepDetectorElementExtension(dd4hepDetElement)));
0407
0408   // return the surface
0409   return dd4hepDetElement->surface().getSharedPtr();
0410 }
0411
0412 Transform3 DD4hepLayerBuilder::convertTransform(
0413     const TGeoMatrix* tGeoTrans) const {
0414   // get the placement and orientation in respect to its mother
0415   const Double_t* rotation = tGeoTrans->GetRotationMatrix();
0416   const Double_t* translation = tGeoTrans->GetTranslation();
0417   return TGeoPrimitivesHelper::makeTransform(
0418       Vector3(rotation[0], rotation[3], rotation[6]),
0419       Vector3(rotation[1], rotation[4], rotation[7]),
0420       Vector3(rotation[2], rotation[5], rotation[8]),
0421       Vector3(translation[0] * UnitConstants::cm,
0422               translation[1] * UnitConstants::cm,
0423               translation[2] * UnitConstants::cm));
0424 }
0425
0426 std::shared_ptr<DD4hepDetectorElement>
0427 DD4hepLayerBuilder::defaultDetectorElementFactory(
0428     const dd4hep::DetElement& detElement, const std::string& detAxis,
0429     double thickness, bool isDisc,
0430     std::shared_ptr<const ISurfaceMaterial> surfaceMaterial) {
0431   return std::make_shared<DD4hepDetectorElement>(
0432       detElement, detAxis, thickness, isDisc, std::move(surfaceMaterial));
0433 }
0434
0435 }  // namespace ActsPlugins