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