TGeo/src/TGeoCylinderDiscSplitter.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2021 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/TGeo/TGeoCylinderDiscSplitter.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Plugins/TGeo/TGeoDetectorElement.hpp"
0013 #include "Acts/Surfaces/CylinderBounds.hpp"
0014 #include "Acts/Surfaces/RadialBounds.hpp"
0015 #include "Acts/Surfaces/RectangleBounds.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017 #include "Acts/Surfaces/SurfaceBounds.hpp"
0018 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0019
0020 #include <cmath>
0021 #include <cstdlib>
0022 #include <utility>
0023
0024 Acts::TGeoCylinderDiscSplitter::TGeoCylinderDiscSplitter(
0025     const TGeoCylinderDiscSplitter::Config& cfg,
0026     std::unique_ptr<const Acts::Logger> logger)
0027     : m_cfg(cfg), m_logger(std::move(logger)) {}
0028
0029 std::vector<std::shared_ptr<const Acts::TGeoDetectorElement>>
0030 Acts::TGeoCylinderDiscSplitter::split(
0031     const GeometryContext& gctx,
0032     std::shared_ptr<const Acts::TGeoDetectorElement> tgde) const {
0033   const Acts::Surface& sf = tgde->surface();
0034   // Thickness
0035   auto tgIdentifier = tgde->identifier();
0036   const TGeoNode& tgNode = tgde->tgeoNode();
0037   ActsScalar tgThickness = tgde->thickness();
0038
0039   // Disc segments are detected, attempt a split
0040   if (m_cfg.discPhiSegments > 0 || m_cfg.discRadialSegments > 0) {
0041     // Splitting for discs detected
0042     if (sf.type() == Acts::Surface::Disc &&
0043         sf.bounds().type() == Acts::SurfaceBounds::eDisc) {
0044       ACTS_DEBUG("- splitting detected for a Disc shaped sensor.");
0045
0046       std::vector<std::shared_ptr<const Acts::TGeoDetectorElement>>
0047           tgDetectorElements = {};
0048       tgDetectorElements.reserve(std::abs(m_cfg.discPhiSegments) *
0049                                  std::abs(m_cfg.discRadialSegments));
0050
0051       const Acts::Vector3 discCenter = sf.center(gctx);
0052
0053       const auto& boundValues = sf.bounds().values();
0054       ActsScalar discMinR = boundValues[Acts::RadialBounds::eMinR];
0055       ActsScalar discMaxR = boundValues[Acts::RadialBounds::eMaxR];
0056
0057       ActsScalar phiStep = 2 * M_PI / m_cfg.discPhiSegments;
0058       ActsScalar cosPhiHalf = std::cos(0.5 * phiStep);
0059       ActsScalar sinPhiHalf = std::sin(0.5 * phiStep);
0060
0061       std::vector<ActsScalar> radialValues = {};
0062       if (m_cfg.discRadialSegments > 1) {
0063         ActsScalar rStep = (discMaxR - discMinR) / m_cfg.discRadialSegments;
0064         radialValues.reserve(m_cfg.discRadialSegments);
0065         for (int ir = 0; ir <= m_cfg.discRadialSegments; ++ir) {
0066           radialValues.push_back(discMinR + ir * rStep);
0067         }
0068       } else {
0069         radialValues = {discMinR, discMaxR};
0070       }
0071
0072       for (std::size_t ir = 1; ir < radialValues.size(); ++ir) {
0073         ActsScalar minR = radialValues[ir - 1];
0074         ActsScalar maxR = radialValues[ir];
0075
0076         ActsScalar maxLocY = maxR * cosPhiHalf;
0077         ActsScalar minLocY = minR * cosPhiHalf;
0078
0079         ActsScalar hR = 0.5 * (maxLocY + minLocY);
0080         ActsScalar hY = 0.5 * (maxLocY - minLocY);
0081         ActsScalar hXminY = minR * sinPhiHalf;
0082         ActsScalar hXmaxY = maxR * sinPhiHalf;
0083
0084         auto tgTrapezoid =
0085             std::make_shared<Acts::TrapezoidBounds>(hXminY, hXmaxY, hY);
0086
0087         for (int im = 0; im < m_cfg.discPhiSegments; ++im) {
0088           // Get the moduleTransform
0089           ActsScalar phi = -M_PI + im * phiStep;
0090           auto tgTransform =
0091               Transform3(Translation3(hR * std::cos(phi), hR * std::sin(phi),
0092                                       discCenter.z()) *
0093                          AngleAxis3(phi - 0.5 * M_PI, Vector3::UnitZ()));
0094
0095           // Create a new detector element per split
0096           auto tgDetectorElement = std::make_shared<Acts::TGeoDetectorElement>(
0097               tgIdentifier, tgNode, tgTransform, tgTrapezoid, tgThickness);
0098
0099           tgDetectorElements.push_back(tgDetectorElement);
0100         }
0101       }
0102
0103       return tgDetectorElements;
0104     }
0105   }
0106
0107   // Cylinder segments are detected, attempt a split
0108   if (m_cfg.cylinderPhiSegments > 0 || m_cfg.cylinderLongitudinalSegments > 0) {
0109     if (sf.type() == Acts::Surface::Cylinder &&
0110         sf.bounds().type() == Acts::SurfaceBounds::eCylinder) {
0111       ACTS_DEBUG("- splitting detected for a Cylinder shaped sensor.");
0112
0113       std::vector<std::shared_ptr<const Acts::TGeoDetectorElement>>
0114           tgDetectorElements = {};
0115       tgDetectorElements.reserve(std::abs(m_cfg.cylinderPhiSegments) *
0116                                  std::abs(m_cfg.cylinderLongitudinalSegments));
0117
0118       const auto& boundValues = sf.bounds().values();
0119       ActsScalar cylinderR = boundValues[Acts::CylinderBounds::eR];
0120       ActsScalar cylinderHalfZ =
0121           boundValues[Acts::CylinderBounds::eHalfLengthZ];
0122
0123       ActsScalar phiStep = 2 * M_PI / m_cfg.cylinderPhiSegments;
0124       ActsScalar cosPhiHalf = std::cos(0.5 * phiStep);
0125       ActsScalar sinPhiHalf = std::sin(0.5 * phiStep);
0126
0127       std::vector<ActsScalar> zValues = {};
0128       if (m_cfg.cylinderLongitudinalSegments > 1) {
0129         ActsScalar zStep =
0130             2 * cylinderHalfZ / m_cfg.cylinderLongitudinalSegments;
0131         zValues.reserve(m_cfg.cylinderLongitudinalSegments);
0132         for (int ir = 0; ir <= m_cfg.cylinderLongitudinalSegments; ++ir) {
0133           zValues.push_back(-cylinderHalfZ + ir * zStep);
0134         }
0135       } else {
0136         zValues = {-cylinderHalfZ, cylinderHalfZ};
0137       }
0138
0139       ActsScalar planeR = cylinderR * cosPhiHalf;
0140       ActsScalar planeHalfX = cylinderR * sinPhiHalf;
0141
0142       for (std::size_t iz = 1; iz < zValues.size(); ++iz) {
0143         ActsScalar minZ = zValues[iz - 1];
0144         ActsScalar maxZ = zValues[iz];
0145
0146         ActsScalar planeZ = 0.5 * (minZ + maxZ);
0147         ActsScalar planeHalfY = 0.5 * (maxZ - minZ);
0148
0149         auto tgRectangle =
0150             std::make_shared<Acts::RectangleBounds>(planeHalfX, planeHalfY);
0151
0152         for (int im = 0; im < m_cfg.cylinderPhiSegments; ++im) {
0153           // Get the moduleTransform
0154           ActsScalar phi = -M_PI + im * phiStep;
0155           ActsScalar cosPhi = std::cos(phi);
0156           ActsScalar sinPhi = std::sin(phi);
0157           ActsScalar planeX = planeR * cosPhi;
0158           ActsScalar planeY = planeR * sinPhi;
0159
0160           Acts::Vector3 planeCenter(planeX, planeY, planeZ);
0161           Acts::Vector3 planeAxisZ = {cosPhi, sinPhi, 0.};
0162           Acts::Vector3 planeAxisY{0., 0., 1.};
0163           Acts::Vector3 planeAxisX = planeAxisY.cross(planeAxisZ);
0164
0165           RotationMatrix3 planeRotation;
0166           planeRotation.col(0) = planeAxisX;
0167           planeRotation.col(1) = planeAxisY;
0168           planeRotation.col(2) = planeAxisZ;
0169
0170           // curvilinear surfaces are boundless
0171           Transform3 planeTransform{planeRotation};
0172           planeTransform.pretranslate(planeCenter);
0173
0174           Transform3 tgTransform = planeTransform;
0175
0176           // Create a new detector element per split
0177           auto tgDetectorElement = std::make_shared<Acts::TGeoDetectorElement>(
0178               tgIdentifier, tgNode, tgTransform, tgRectangle, tgThickness);
0179
0180           tgDetectorElements.push_back(tgDetectorElement);
0181         }
0182       }
0183       return tgDetectorElements;
0184     }
0185   }
0186   return {tgde};
0187 }