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