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