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