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