DD4hep/src/DD4hepLayerBuilder.cpp

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