Detector/detail/SupportSurfacesHelper.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2023 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/Detector/detail/SupportSurfacesHelper.hpp"
0010
0011 #include "Acts/Definitions/Tolerance.hpp"
0012 #include "Acts/Surfaces/CylinderBounds.hpp"
0013 #include "Acts/Surfaces/CylinderSurface.hpp"
0014 #include "Acts/Surfaces/DiscSurface.hpp"
0015 #include "Acts/Surfaces/PlaneSurface.hpp"
0016 #include "Acts/Surfaces/RadialBounds.hpp"
0017 #include "Acts/Surfaces/RectangleBounds.hpp"
0018 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0019 #include "Acts/Utilities/BinningType.hpp"
0020
0021 #include <algorithm>
0022 #include <cmath>
0023 #include <stdexcept>
0024 #include <utility>
0025
0026 Acts::Experimental::detail::SupportSurfacesHelper::SupportSurfaceComponents
0027 Acts::Experimental::detail::SupportSurfacesHelper::CylindricalSupport::
0028 operator()(const Extent& lExtent) const {
0029   // Bail out if you have no measure of R, Z
0030   if (!lExtent.constrains(binZ) || !lExtent.constrains(binR)) {
0031     throw std::invalid_argument(
0032         "SupportSurfacesHelper::CylindricalSupport::operator() - z or "
0033         "r are not constrained.");
0034   }
0035
0036   // Min / Max z  with clearances adapted
0037   ActsScalar minZ = lExtent.min(binZ) + std::abs(zClearance[0u]);
0038   ActsScalar maxZ = lExtent.max(binZ) - std::abs(zClearance[1u]);
0039
0040   // Phi sector
0041   ActsScalar hPhiSector = M_PI;
0042   ActsScalar avgPhi = 0.;
0043   if (lExtent.constrains(binPhi)) {
0044     // Min / Max phi  with clearances adapted
0045     ActsScalar minPhi = lExtent.min(binPhi) + std::abs(phiClearance[0u]);
0046     ActsScalar maxPhi = lExtent.max(binPhi) - std::abs(phiClearance[1u]);
0047     hPhiSector = 0.5 * (maxPhi - minPhi);
0048     avgPhi = 0.5 * (minPhi + maxPhi);
0049   }
0050
0051   Transform3 transform = Transform3::Identity();
0052   if (std::abs(minZ + maxZ) > s_onSurfaceTolerance) {
0053     transform.pretranslate(Vector3(0., 0., 0.5 * (minZ + maxZ)));
0054   }
0055
0056   // The Radius estimation
0057   ActsScalar r =
0058       rOffset < 0 ? lExtent.min(binR) + rOffset : lExtent.max(binR) + rOffset;
0059   if (rOffset == 0.) {
0060     r = lExtent.medium(binR);
0061   }
0062   // Components are resolved and returned as a tuple
0063   return {
0064       type, {r, 0.5 * (maxZ - minZ), hPhiSector, avgPhi, 0., 0.}, transform};
0065 }
0066
0067 Acts::Experimental::detail::SupportSurfacesHelper::SupportSurfaceComponents
0068 Acts::Experimental::detail::SupportSurfacesHelper::DiscSupport::operator()(
0069     const Extent& lExtent) const {
0070   // Bail out if you have no measure of R, Z
0071   if (!lExtent.constrains(binZ) || !lExtent.constrains(binR)) {
0072     throw std::invalid_argument(
0073         "SupportSurfacesHelper::DiscSupport::operator() - z or "
0074         "r are not constrained.");
0075   }
0076
0077   // Min / Max r  with clearances adapted
0078   ActsScalar minR = lExtent.min(binR) + std::abs(rClearance[0u]);
0079   ActsScalar maxR = lExtent.max(binR) - std::abs(rClearance[1u]);
0080
0081   // Phi sector
0082   ActsScalar hPhiSector = M_PI;
0083   ActsScalar avgPhi = 0.;
0084   if (lExtent.constrains(binPhi)) {
0085     // Min / Max phi  with clearances adapted
0086     ActsScalar minPhi = lExtent.min(binPhi) + std::abs(phiClearance[0u]);
0087     ActsScalar maxPhi = lExtent.max(binPhi) - std::abs(phiClearance[1u]);
0088     hPhiSector = 0.5 * (maxPhi - minPhi);
0089     avgPhi = 0.5 * (minPhi + maxPhi);
0090   }
0091
0092   // The z position estimate
0093   ActsScalar z =
0094       zOffset < 0 ? lExtent.min(binZ) + zOffset : lExtent.max(binZ) + zOffset;
0095   if (zOffset == 0.) {
0096     z = lExtent.medium(binZ);
0097   }
0098
0099   Transform3 transform = Transform3::Identity();
0100   transform.pretranslate(Vector3(0., 0., z));
0101
0102   // Components are resolved and returned as a tuple
0103   return {type, {minR, maxR, hPhiSector, avgPhi}, transform};
0104 }
0105
0106 Acts::Experimental::detail::SupportSurfacesHelper::SupportSurfaceComponents
0107 Acts::Experimental::detail::SupportSurfacesHelper::RectangularSupport::
0108 operator()(const Extent& lExtent) const {
0109   // Bail out if you have no measure of X, Y, Z
0110   if (!(lExtent.constrains(binX) && lExtent.constrains(binY) &&
0111         lExtent.constrains(binZ))) {
0112     throw std::invalid_argument(
0113         "SupportSurfacesHelper::RectangularSupport::operator() - x, y or "
0114         "z are not constrained.");
0115   }
0116
0117   // Set the local coordinates - cyclic permutation
0118   std::array<BinningValue, 2> locals = {binX, binY};
0119   if (pPlacement == binX) {
0120     locals = {binY, binZ};
0121   } else if (pPlacement == binY) {
0122     locals = {binZ, binX};
0123   }
0124
0125   // Make the rectangular shape
0126   ActsScalar minX = lExtent.min(locals[0]) + std::abs(loc0Clearance[0u]);
0127   ActsScalar maxX = lExtent.max(locals[0]) - std::abs(loc0Clearance[1u]);
0128   ActsScalar minY = lExtent.min(locals[1]) + std::abs(loc1Clearance[0u]);
0129   ActsScalar maxY = lExtent.max(locals[1]) - std::abs(loc1Clearance[1u]);
0130
0131   ActsScalar gPlacement = lExtent.medium(pPlacement) + pOffset;
0132   Vector3 placement = Vector3::Zero();
0133   placement[pPlacement] = gPlacement;
0134
0135   Transform3 transform = Transform3::Identity();
0136   transform.pretranslate(placement);
0137
0138   return {type, {minX, minY, maxX, maxY}, transform};
0139 }
0140
0141 std::vector<std::shared_ptr<Acts::Surface>>
0142 Acts::Experimental::detail::SupportSurfacesHelper::cylindricalSupport(
0143     const SupportSurfaceComponents& components, unsigned int splits) {
0144   // Resolve the components
0145   auto [type, values, transform] = components;
0146
0147   // Parameter size check
0148   if (values.size() != 6u) {
0149     throw std::invalid_argument(
0150         "SupportSurfacesHelper::cylindricalSupport(...) - "
0151         "values vector has wrong size, requires 6 parameters.");
0152   }
0153
0154   // Surface type check
0155   if (type != Surface::SurfaceType::Cylinder) {
0156     throw std::invalid_argument(
0157         "SupportSurfacesHelper::cylindricalSupport(...) - "
0158         "surface type is not a cylinder.");
0159   }
0160
0161   std::array<ActsScalar, 6u> bounds = {};
0162   std::copy_n(values.begin(), 6u, bounds.begin());
0163
0164   // Return vector for generated surfaces
0165   std::vector<std::shared_ptr<Acts::Surface>> cSupport;
0166   if (splits == 1u) {
0167     // No splitting is done in this case
0168     cSupport.push_back(Surface::makeShared<CylinderSurface>(
0169         transform, std::make_shared<CylinderBounds>(bounds)));
0170   } else {
0171     // Split into n(splits) planar surfaces, prep work:
0172     ActsScalar r = bounds[0u];
0173     ActsScalar halfZ = bounds[1u];
0174     ActsScalar minPhi = bounds[3u] - bounds[2u];
0175     ActsScalar maxPhi = bounds[3u] + bounds[2u];
0176     ActsScalar dHalfPhi = (maxPhi - minPhi) / (2 * splits);
0177     ActsScalar cosPhiHalf = std::cos(dHalfPhi);
0178     ActsScalar sinPhiHalf = std::sin(dHalfPhi);
0179     ActsScalar planeR = r * cosPhiHalf;
0180     ActsScalar planeHalfX = r * sinPhiHalf;
0181     ActsScalar planeZ = transform.translation().z();
0182
0183     auto sRectangle =
0184         std::make_shared<Acts::RectangleBounds>(planeHalfX, halfZ);
0185     // Now create the Trapezoids
0186     for (unsigned int iphi = 0; iphi < splits; ++iphi) {
0187       // Get the moduleTransform
0188       ActsScalar phi = -M_PI + (iphi + 0.5) * 2 * dHalfPhi;
0189       ActsScalar cosPhi = std::cos(phi);
0190       ActsScalar sinPhi = std::sin(phi);
0191       ActsScalar planeX = planeR * cosPhi;
0192       ActsScalar planeY = planeR * sinPhi;
0193
0194       Acts::Vector3 planeCenter(planeX, planeY, planeZ);
0195       Acts::Vector3 planeAxisZ(cosPhi, sinPhi, 0.);
0196       Acts::Vector3 planeAxisY(0., 0., 1.);
0197       Acts::Vector3 planeAxisX = planeAxisY.cross(planeAxisZ);
0198
0199       RotationMatrix3 planeRotation;
0200       planeRotation.col(0) = planeAxisX;
0201       planeRotation.col(1) = planeAxisY;
0202       planeRotation.col(2) = planeAxisZ;
0203
0204       Transform3 sTransform{planeRotation};
0205       sTransform.pretranslate(planeCenter);
0206       // Place it
0207       cSupport.push_back(
0208           Surface::makeShared<PlaneSurface>(sTransform, sRectangle));
0209     }
0210   }
0211
0212   return cSupport;
0213 }
0214
0215 std::vector<std::shared_ptr<Acts::Surface>>
0216 Acts::Experimental::detail::SupportSurfacesHelper::discSupport(
0217     const SupportSurfaceComponents& components, unsigned int splits) {
0218   // Resolve the components
0219   auto [type, values, transform] = components;
0220
0221   // Parameter size check
0222   if (values.size() != 4u) {
0223     throw std::invalid_argument(
0224         "SupportSurfacesHelper::discSupport(...) - "
0225         "values vector has wrong size, requires 4 parameters.");
0226   }
0227
0228   // Surface type check
0229   if (type != Surface::SurfaceType::Disc) {
0230     throw std::invalid_argument(
0231         "SupportSurfacesHelper::discSupport(...) - "
0232         "surface type is not a disc.");
0233   }
0234
0235   std::array<ActsScalar, 4u> bounds = {};
0236   std::copy_n(values.begin(), 4u, bounds.begin());
0237
0238   // Return vector for generated surfaces
0239   std::vector<std::shared_ptr<Acts::Surface>> dSupport;
0240   if (splits == 1u) {
0241     // No splitting is done in this case
0242     dSupport.push_back(Surface::makeShared<DiscSurface>(
0243         transform, std::make_shared<RadialBounds>(bounds)));
0244   } else {
0245     // Split into n(splits) planar surfaces in phi, prep work:
0246     ActsScalar minR = bounds[0u];
0247     ActsScalar maxR = bounds[1u];
0248     ActsScalar minPhi = bounds[3u] - bounds[2u];
0249     ActsScalar maxPhi = bounds[3u] + bounds[2u];
0250     ActsScalar dHalfPhi = (maxPhi - minPhi) / (2 * splits);
0251     ActsScalar cosPhiHalf = std::cos(dHalfPhi);
0252     ActsScalar sinPhiHalf = std::sin(dHalfPhi);
0253     ActsScalar maxLocY = maxR * cosPhiHalf;
0254     ActsScalar minLocY = minR * cosPhiHalf;
0255     ActsScalar hR = 0.5 * (maxLocY + minLocY);
0256     ActsScalar hY = 0.5 * (maxLocY - minLocY);
0257     ActsScalar hXminY = minR * sinPhiHalf;
0258     ActsScalar hXmaxY = maxR * sinPhiHalf;
0259     // Split trapezoid
0260     auto sTrapezoid =
0261         std::make_shared<Acts::TrapezoidBounds>(hXminY, hXmaxY, hY);
0262     Vector3 zAxis = transform.rotation().col(2);
0263     ActsScalar zPosition = transform.translation().z();
0264     // Now create the Trapezoids
0265     for (unsigned int iphi = 0; iphi < splits; ++iphi) {
0266       // Create the split module transform
0267       ActsScalar phi = -M_PI + (iphi + 0.5) * 2 * dHalfPhi;
0268       auto sTransform = Transform3(
0269           Translation3(hR * std::cos(phi), hR * std::sin(phi), zPosition) *
0270           AngleAxis3(phi - 0.5 * M_PI, zAxis));
0271       // Place it
0272       dSupport.push_back(
0273           Surface::makeShared<PlaneSurface>(sTransform, sTrapezoid));
0274     }
0275   }
0276   return dSupport;
0277 }
0278
0279 std::vector<std::shared_ptr<Acts::Surface>>
0280 Acts::Experimental::detail::SupportSurfacesHelper::rectangularSupport(
0281     const SupportSurfaceComponents& components) {
0282   // Resolve the components
0283   auto [type, values, transform] = components;
0284
0285   // Parameter size check
0286   if (values.size() != 4u) {
0287     throw std::invalid_argument(
0288         "SupportSurfacesHelper::rectangularSupport(...) - "
0289         "values vector has wrong size, requires 4 parameters.");
0290   }
0291
0292   // Surface type check
0293   if (type != Surface::SurfaceType::Plane) {
0294     throw std::invalid_argument(
0295         "SupportSurfacesHelper::rectangularSupport(...) - "
0296         "surface type is not a plane.");
0297   }
0298
0299   std::array<ActsScalar, 4u> bounds = {};
0300   std::copy_n(values.begin(), 4u, bounds.begin());
0301
0302   return {Surface::makeShared<PlaneSurface>(
0303       transform, std::make_shared<RectangleBounds>(bounds))};
0304 }
0305
0306 void Acts::Experimental::detail::SupportSurfacesHelper::addSupport(
0307     std::vector<std::shared_ptr<Surface>>& layerSurfaces,
0308     std::vector<std::size_t>& assignToAll, const Extent& layerExtent,
0309     const SurfaceComponentsCreator& componentCreator,
0310     unsigned int supportSplits) {
0311   // Get the main support surface components
0312   auto supportComponents = componentCreator(layerExtent);
0313   const auto& sType = std::get<0>(supportComponents);
0314
0315   std::vector<std::shared_ptr<Acts::Surface>> supportSurfaces = {};
0316
0317   if (sType == Surface::SurfaceType::Cylinder) {
0318     supportSurfaces = cylindricalSupport(supportComponents, supportSplits);
0319   } else if (sType == Surface::SurfaceType::Disc) {
0320     supportSurfaces = discSupport(supportComponents, supportSplits);
0321   } else if (sType == Surface::SurfaceType::Plane) {
0322     supportSurfaces = rectangularSupport(supportComponents);
0323   } else {
0324     throw std::invalid_argument(
0325         "SupportSurfacesHelper: currently only cylindrical/disc/rectangle "
0326         "support building is possible.");
0327   }
0328
0329   // Remember the surfaces to be assigned to all bins, once the
0330   // support surfaces are split they enter the standard bin assignment
0331   if (supportSplits == 1u && supportSurfaces.size() == 1u) {
0332     assignToAll.push_back(layerSurfaces.size());
0333   }
0334   // Add those to the layer surfaces
0335   layerSurfaces.insert(layerSurfaces.end(), supportSurfaces.begin(),
0336                        supportSurfaces.end());
0337 }