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 <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 namespace Acts {
0047
0048 DD4hepLayerBuilder::DD4hepLayerBuilder(const DD4hepLayerBuilder::Config& config,
0049                                        std::unique_ptr<const Logger> logger)
0050     : m_cfg(), m_logger(std::move(logger)) {
0051   setConfiguration(config);
0052 }
0053
0054 DD4hepLayerBuilder::~DD4hepLayerBuilder() = default;
0055
0056 void DD4hepLayerBuilder::setConfiguration(
0057     const DD4hepLayerBuilder::Config& config) {
0058   m_cfg = config;
0059 }
0060
0061 const LayerVector 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       auto itransform = transform.inverse();
0087       // get the shape of the layer
0088       TGeoShape* geoShape =
0089           detElement.placement().ptr()->GetVolume()->GetShape();
0090       // create the proto layer
0091       ProtoLayer pl(gctx, layerSurfaces, itransform);
0092       if (logger().doPrint(Logging::VERBOSE)) {
0093         std::stringstream ss;
0094         pl.toStream(ss);
0095         ACTS_VERBOSE("Extent from surfaces: " << ss.str());
0096
0097         std::vector<double> rvalues;
0098         std::transform(layerSurfaces.begin(), layerSurfaces.end(),
0099                        std::back_inserter(rvalues), [&](const auto& surface) {
0100                          return VectorHelpers::perp(surface->center(gctx));
0101                        });
0102         std::ranges::sort(rvalues);
0103         std::vector<std::string> locs;
0104         std::transform(rvalues.begin(),
0105                        std::unique(rvalues.begin(), rvalues.end()),
0106                        std::back_inserter(locs),
0107                        [](const auto& v) { return std::to_string(v); });
0108         ACTS_VERBOSE(
0109             "-> unique r locations: " << boost::algorithm::join(locs, ", "));
0110       }
0111
0112       if (params.contains("envelope_r_min") &&
0113           params.contains("envelope_r_max") &&
0114           params.contains("envelope_z_min") &&
0115           params.contains("envelope_z_max")) {
0116         // set the values of the proto layer in case enevelopes are handed
0117         // over
0118         pl.envelope[AxisDirection::AxisR] = {
0119             params.get<double>("envelope_r_min"),
0120             params.get<double>("envelope_r_max")};
0121         pl.envelope[AxisDirection::AxisZ] = {
0122             params.get<double>("envelope_z_min"),
0123             params.get<double>("envelope_z_max")};
0124       } else if (geoShape != nullptr) {
0125         TGeoTubeSeg* tube = dynamic_cast<TGeoTubeSeg*>(geoShape);
0126         if (tube == nullptr) {
0127           ACTS_ERROR(" Disc layer has wrong shape - needs to be TGeoTubeSeg!");
0128           throw std::logic_error{
0129               "Disc layer has wrong shape - needs to be TGeoTubeSeg!"};
0130         }
0131         // extract the boundaries
0132         double rMin = tube->GetRmin() * UnitConstants::cm;
0133         double rMax = tube->GetRmax() * UnitConstants::cm;
0134         double zMin =
0135             (transform.translation() -
0136              transform.rotation().col(2) * tube->GetDz() * UnitConstants::cm)
0137                 .z();
0138         double zMax =
0139             (transform.translation() +
0140              transform.rotation().col(2) * tube->GetDz() * UnitConstants::cm)
0141                 .z();
0142         if (zMin > zMax) {
0143           std::swap(zMin, zMax);
0144         }
0145         // check if layer has surfaces
0146         if (layerSurfaces.empty()) {
0147           ACTS_VERBOSE(" Disc layer has no sensitive surfaces.");
0148           // in case no surfaces are handed over the layer thickness will be
0149           // set to a default value to allow attaching material layers
0150           double z = (zMin + zMax) * 0.5;
0151           // create layer without surfaces
0152           // manually create a proto layer
0153           double eiz = (z != 0.) ? z - m_cfg.defaultThickness : 0.;
0154           double eoz = (z != 0.) ? z + m_cfg.defaultThickness : 0.;
0155           pl.extent.range(AxisDirection::AxisZ).set(eiz, eoz);
0156           pl.extent.range(AxisDirection::AxisR).set(rMin, rMax);
0157           pl.envelope[AxisDirection::AxisR] = {0., 0.};
0158           pl.envelope[AxisDirection::AxisZ] = {0., 0.};
0159         } else {
0160           ACTS_VERBOSE(" Disc layer has " << layerSurfaces.size()
0161                                           << " sensitive surfaces.");
0162           // set the values of the proto layer in case dimensions are given by
0163           // geometry
0164           pl.envelope[AxisDirection::AxisZ] = {
0165               std::abs(zMin - pl.min(AxisDirection::AxisZ)),
0166               std::abs(zMax - pl.max(AxisDirection::AxisZ))};
0167           pl.envelope[AxisDirection::AxisR] = {
0168               std::abs(rMin - pl.min(AxisDirection::AxisR)),
0169               std::abs(rMax - pl.max(AxisDirection::AxisR))};
0170           pl.extent.range(AxisDirection::AxisR).set(rMin, rMax);
0171         }
0172       } else {
0173         throw std::logic_error(
0174             std::string("Layer DetElement: ") + detElement.name() +
0175             std::string(" has neither a shape nor tolerances for envelopes "
0176                         "added to its extension. Please check your detector "
0177                         "constructor!"));
0178       }
0179
0180       std::shared_ptr<Layer> endcapLayer = nullptr;
0181
0182       // Check if DD4hep pre-defines the surface binning
0183       bool hasSurfaceBinning =
0184           getParamOr<bool>("surface_binning", detElement, true);
0185       std::size_t nPhi = 1;
0186       std::size_t nR = 1;
0187       if (hasSurfaceBinning) {
0188         if (params.contains("surface_binning_n_phi")) {
0189           nPhi = static_cast<std::size_t>(
0190               params.get<int>("surface_binning_n_phi"));
0191         }
0192         if (params.contains("surface_binning_n_r")) {
0193           nR = static_cast<std::size_t>(params.get<int>("surface_binning_n_r"));
0194         }
0195         hasSurfaceBinning = nR * nPhi > 1;
0196       }
0197
0198       // In case the layer is sensitive
0199       if (detElement.volume().isSensitive()) {
0200         // Create the sensitive surface
0201         auto sensitiveSurf = createSensitiveSurface(detElement, true);
0202         // Create the surfaceArray
0203         auto sArray = std::make_unique<SurfaceArray>(sensitiveSurf);
0204
0205         // create the share disc bounds
0206         auto dBounds = std::make_shared<const RadialBounds>(
0207             pl.min(AxisDirection::AxisR), pl.max(AxisDirection::AxisR));
0208         double thickness = std::abs(pl.max(AxisDirection::AxisZ) -
0209                                     pl.min(AxisDirection::AxisZ));
0210         // Create the layer containing the sensitive surface
0211         endcapLayer = DiscLayer::create(transform, dBounds, std::move(sArray),
0212                                         thickness, nullptr, 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 LayerVector DD4hepLayerBuilder::negativeLayers(
0234     const GeometryContext& gctx) const {
0235   return endcapLayers(gctx, m_cfg.negativeLayers, "negative");
0236 }
0237
0238 const LayerVector 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[AxisDirection::AxisR] = {
0291             params.get<double>("envelope_r_min"),
0292             params.get<double>("envelope_r_max")};
0293         pl.envelope[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(AxisDirection::AxisR).set(eir, eor);
0319           pl.extent.range(AxisDirection::AxisZ).set(-dz, dz);
0320           pl.envelope[AxisDirection::AxisR] = {0., 0.};
0321           pl.envelope[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[AxisDirection::AxisZ] = {
0326               std::abs(-dz - pl.min(AxisDirection::AxisZ)),
0327               std::abs(dz - pl.max(AxisDirection::AxisZ))};
0328           pl.envelope[AxisDirection::AxisR] = {
0329               std::abs(rMin - pl.min(AxisDirection::AxisR)),
0330               std::abs(rMax - pl.max(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 =
0341           (pl.min(AxisDirection::AxisZ) - pl.max(AxisDirection::AxisZ)) * 0.5;
0342
0343       std::shared_ptr<Layer> centralLayer = nullptr;
0344       // In case the layer is sensitive
0345       if (detElement.volume().isSensitive()) {
0346         // Create the sensitive surface
0347         auto sensitiveSurf = createSensitiveSurface(detElement);
0348         // Create the surfaceArray
0349         std::unique_ptr<SurfaceArray> sArray =
0350             std::make_unique<SurfaceArray>(sensitiveSurf);
0351
0352         // create the layer
0353         double layerR =
0354             (pl.min(AxisDirection::AxisR) + pl.max(AxisDirection::AxisR)) * 0.5;
0355         double thickness = std::abs(pl.max(AxisDirection::AxisR) -
0356                                     pl.min(AxisDirection::AxisR));
0357         auto cBounds = std::make_shared<CylinderBounds>(layerR, halfZ);
0358         // Create the layer containing the sensitive surface
0359         centralLayer = CylinderLayer::create(
0360             transform, cBounds, std::move(sArray), thickness, nullptr, active);
0361
0362       } else {
0363         centralLayer = m_cfg.layerCreator->cylinderLayer(
0364             gctx, layerSurfaces, m_cfg.bTypePhi, m_cfg.bTypeZ, pl, transform,
0365             nullptr);
0366       }
0367       // Add the ProtoMaterial if present
0368       addCylinderLayerProtoMaterial(detElement, *centralLayer, logger());
0369       // push back created layer
0370       layers.push_back(centralLayer);
0371     }
0372   }
0373   return layers;
0374 }
0375
0376 const LayerVector DD4hepLayerBuilder::positiveLayers(
0377     const GeometryContext& gctx) const {
0378   return endcapLayers(gctx, m_cfg.positiveLayers, "positive");
0379 }
0380
0381 void DD4hepLayerBuilder::resolveSensitive(
0382     const dd4hep::DetElement& detElement,
0383     std::vector<std::shared_ptr<const Surface>>& surfaces) const {
0384   const dd4hep::DetElement::Children& children = detElement.children();
0385   if (!children.empty()) {
0386     for (auto& child : children) {
0387       dd4hep::DetElement childDetElement = child.second;
0388       if (childDetElement.volume().isSensitive()) {
0389         // create the surface
0390         surfaces.push_back(createSensitiveSurface(childDetElement, false));
0391       }
0392       resolveSensitive(childDetElement, surfaces);
0393     }
0394   }
0395 }
0396
0397 std::shared_ptr<const Surface> DD4hepLayerBuilder::createSensitiveSurface(
0398     const dd4hep::DetElement& detElement, bool isDisc) const {
0399   std::string detAxis =
0400       getParamOr<std::string>("axis_definitions", detElement, "XYZ");
0401   // Create the corresponding detector element !- memory leak --!
0402   auto dd4hepDetElement = m_cfg.detectorElementFactory(
0403       detElement, detAxis, UnitConstants::cm, isDisc, nullptr);
0404
0405   detElement.addExtension<DD4hepDetectorElementExtension>(
0406       new dd4hep::rec::StructExtension(
0407           DD4hepDetectorElementExtension(dd4hepDetElement)));
0408
0409   // return the surface
0410   return dd4hepDetElement->surface().getSharedPtr();
0411 }
0412
0413 Transform3 DD4hepLayerBuilder::convertTransform(
0414     const TGeoMatrix* tGeoTrans) const {
0415   // get the placement and orientation in respect to its mother
0416   const Double_t* rotation = tGeoTrans->GetRotationMatrix();
0417   const Double_t* translation = tGeoTrans->GetTranslation();
0418   return TGeoPrimitivesHelper::makeTransform(
0419       Vector3(rotation[0], rotation[3], rotation[6]),
0420       Vector3(rotation[1], rotation[4], rotation[7]),
0421       Vector3(rotation[2], rotation[5], rotation[8]),
0422       Vector3(translation[0] * UnitConstants::cm,
0423               translation[1] * UnitConstants::cm,
0424               translation[2] * UnitConstants::cm));
0425 }
0426
0427 std::shared_ptr<DD4hepDetectorElement>
0428 DD4hepLayerBuilder::defaultDetectorElementFactory(
0429     const dd4hep::DetElement& detElement, const std::string& detAxis,
0430     double thickness, bool isDisc,
0431     std::shared_ptr<const ISurfaceMaterial> surfaceMaterial) {
0432   return std::make_shared<DD4hepDetectorElement>(
0433       detElement, detAxis, thickness, isDisc, std::move(surfaceMaterial));
0434 }
0435
0436 }  // namespace Acts