src/Digitization/PlanarSurfaceMask.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 "ActsFatras/Digitization/PlanarSurfaceMask.hpp"
0010
0011 #include "Acts/Definitions/Tolerance.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0014 #include "Acts/Surfaces/SurfaceBounds.hpp"
0015 #include "Acts/Utilities/Intersection.hpp"
0016 #include "ActsFatras/Digitization/DigitizationError.hpp"
0017 #include <Acts/Surfaces/AnnulusBounds.hpp>
0018 #include <Acts/Surfaces/DiscTrapezoidBounds.hpp>
0019 #include <Acts/Surfaces/PlanarBounds.hpp>
0020 #include <Acts/Surfaces/RadialBounds.hpp>
0021 #include <Acts/Surfaces/Surface.hpp>
0022 #include <Acts/Utilities/Helpers.hpp>
0023
0024 #include <algorithm>
0025 #include <cmath>
0026 #include <cstddef>
0027 #include <memory>
0028 #include <numbers>
0029
0030 namespace {
0031
0032 /// Helper method to check if an intersection is good.
0033 ///
0034 /// Good in this context is defined as: along direction,
0035 /// closer than the segment length & reachable
0036 ///
0037 /// @param intersections The confirmed intersections for the segment
0038 /// @param candidate The candidate intersection
0039 /// @param sLength The segment length, maximal allowed length
0040 void checkIntersection(std::vector<Acts::Intersection2D>& intersections,
0041                        const Acts::Intersection2D& candidate, double sLength) {
0042   if (candidate.isValid() && candidate.pathLength() > 0 &&
0043       candidate.pathLength() < sLength) {
0044     intersections.push_back(candidate);
0045   }
0046 }
0047
0048 /// Helper method to apply the mask and return.
0049 ///
0050 /// If two (or more) intersections would be good, apply the first two
0051 /// If only one is available, the boolean tells you which one it is.
0052 /// If no intersection is valid, return an error code for masking.
0053 ///
0054 /// @param intersections All confirmed intersections
0055 /// @param segment The original segment before masking
0056 /// @param firstInside Indicator if the first is inside or not
0057 ///
0058 /// @return a new Segment (clipped) wrapped in a result or error_code
0059 Acts::Result<ActsFatras::PlanarSurfaceMask::Segment2D> maskAndReturn(
0060     std::vector<Acts::Intersection2D>& intersections,
0061     const ActsFatras::PlanarSurfaceMask::Segment2D& segment, bool firstInside) {
0062   std::ranges::sort(intersections, Acts::Intersection2D::pathLengthOrder);
0063   if (intersections.size() >= 2) {
0064     return ActsFatras::PlanarSurfaceMask::Segment2D{
0065         intersections[0].position(), intersections[1].position()};
0066   } else if (intersections.size() == 1) {
0067     return (!firstInside
0068                 ? ActsFatras::PlanarSurfaceMask::Segment2D{intersections[0]
0069                                                                .position(),
0070                                                            segment[1]}
0071                 : ActsFatras::PlanarSurfaceMask::Segment2D{
0072                       segment[0], intersections[0].position()});
0073   }
0074   return ActsFatras::DigitizationError::MaskingError;
0075 }
0076
0077 }  // anonymous namespace
0078
0079 Acts::Result<ActsFatras::PlanarSurfaceMask::Segment2D>
0080 ActsFatras::PlanarSurfaceMask::apply(const Acts::Surface& surface,
0081                                      const Segment2D& segment) const {
0082   auto surfaceType = surface.type();
0083
0084   // Plane surface section -------------------
0085   if (surfaceType == Acts::Surface::Plane ||
0086       surface.bounds().type() == Acts::SurfaceBounds::eDiscTrapezoid) {
0087     Acts::Vector2 localStart =
0088         (surfaceType == Acts::Surface::Plane)
0089             ? segment[0]
0090             : Acts::Vector2(Acts::VectorHelpers::perp(segment[0]),
0091                             Acts::VectorHelpers::phi(segment[0]));
0092
0093     Acts::Vector2 localEnd =
0094         (surfaceType == Acts::Surface::Plane)
0095             ? segment[1]
0096             : Acts::Vector2(Acts::VectorHelpers::perp(segment[1]),
0097                             Acts::VectorHelpers::phi(segment[1]));
0098
0099     bool startInside =
0100         surface.bounds().inside(localStart, Acts::BoundaryTolerance::None());
0101     bool endInside =
0102         surface.bounds().inside(localEnd, Acts::BoundaryTolerance::None());
0103
0104     // Fast exit, both inside
0105     if (startInside && endInside) {
0106       return segment;
0107     }
0108
0109     // It's either planar or disc trapezoid bounds
0110     const Acts::PlanarBounds* planarBounds = nullptr;
0111     const Acts::DiscTrapezoidBounds* dtbBounds = nullptr;
0112     if (surfaceType == Acts::Surface::Plane) {
0113       planarBounds =
0114           static_cast<const Acts::PlanarBounds*>(&(surface.bounds()));
0115       if (planarBounds->type() == Acts::SurfaceBounds::eEllipse) {
0116         return DigitizationError::UndefinedSurface;
0117       }
0118     } else {
0119       dtbBounds =
0120           static_cast<const Acts::DiscTrapezoidBounds*>(&(surface.bounds()));
0121     }
0122     auto vertices = planarBounds != nullptr ? planarBounds->vertices(1)
0123                                             : dtbBounds->vertices(1);
0124
0125     return polygonMask(vertices, segment, startInside);
0126
0127   } else if (surfaceType == Acts::Surface::Disc) {
0128     // Polar coordinates
0129     Acts::Vector2 sPolar(Acts::VectorHelpers::perp(segment[0]),
0130                          Acts::VectorHelpers::phi(segment[0]));
0131     Acts::Vector2 ePolar(Acts::VectorHelpers::perp(segment[1]),
0132                          Acts::VectorHelpers::phi(segment[1]));
0133
0134     bool startInside =
0135         surface.bounds().inside(sPolar, Acts::BoundaryTolerance::None());
0136     bool endInside =
0137         surface.bounds().inside(ePolar, Acts::BoundaryTolerance::None());
0138
0139     // Fast exit for both inside
0140     if (startInside && endInside) {
0141       return segment;
0142     }
0143
0144     auto boundsType = surface.bounds().type();
0145     if (boundsType == Acts::SurfaceBounds::eDisc) {
0146       auto rBounds =
0147           static_cast<const Acts::RadialBounds*>(&(surface.bounds()));
0148       return radialMask(*rBounds, segment, {sPolar, ePolar}, startInside);
0149
0150     } else if (boundsType == Acts::SurfaceBounds::eAnnulus) {
0151       auto aBounds =
0152           static_cast<const Acts::AnnulusBounds*>(&(surface.bounds()));
0153       return annulusMask(*aBounds, segment, startInside);
0154     }
0155   }
0156   return DigitizationError::UndefinedSurface;
0157 }
0158
0159 Acts::Result<ActsFatras::PlanarSurfaceMask::Segment2D>
0160 ActsFatras::PlanarSurfaceMask::polygonMask(
0161     const std::vector<Acts::Vector2>& vertices, const Segment2D& segment,
0162     bool firstInside) const {
0163   std::vector<Acts::Intersection2D> intersections;
0164   Acts::Vector2 sVector(segment[1] - segment[0]);
0165   Acts::Vector2 sDir = sVector.normalized();
0166   double sLength = sVector.norm();
0167
0168   for (std::size_t iv = 0; iv < vertices.size(); ++iv) {
0169     const Acts::Vector2& s0 = vertices[iv];
0170     const Acts::Vector2& s1 =
0171         (iv + 1) < vertices.size() ? vertices[iv + 1] : vertices[0];
0172     checkIntersection(
0173         intersections,
0174         intersector.intersectSegment(s0, s1, segment[0], sDir, true), sLength);
0175   }
0176   return maskAndReturn(intersections, segment, firstInside);
0177 }
0178
0179 Acts::Result<ActsFatras::PlanarSurfaceMask::Segment2D>
0180 ActsFatras::PlanarSurfaceMask::radialMask(const Acts::RadialBounds& rBounds,
0181                                           const Segment2D& segment,
0182                                           const Segment2D& polarSegment,
0183                                           bool firstInside) const {
0184   double rMin = rBounds.get(Acts::RadialBounds::eMinR);
0185   double rMax = rBounds.get(Acts::RadialBounds::eMaxR);
0186   double hPhi = rBounds.get(Acts::RadialBounds::eHalfPhiSector);
0187   double aPhi = rBounds.get(Acts::RadialBounds::eAveragePhi);
0188
0189   std::array<double, 2> radii = {rMin, rMax};
0190   std::array<double, 2> phii = {aPhi - hPhi, aPhi + hPhi};
0191
0192   std::vector<Acts::Intersection2D> intersections;
0193   Acts::Vector2 sVector(segment[1] - segment[0]);
0194   Acts::Vector2 sDir = sVector.normalized();
0195   double sLength = sVector.norm();
0196
0197   double sR = polarSegment[0][Acts::eBoundLoc0];
0198   double eR = polarSegment[1][Acts::eBoundLoc0];
0199   double sPhi = polarSegment[0][Acts::eBoundLoc1];
0200   double ePhi = polarSegment[1][Acts::eBoundLoc1];
0201
0202   // Helper method to intersect phi boundaries
0203   auto intersectPhiLine = [&](double phi) -> void {
0204     Acts::Vector2 s0(rMin * std::cos(phi), rMin * std::sin(phi));
0205     Acts::Vector2 s1(rMax * std::cos(phi), rMax * std::sin(phi));
0206     checkIntersection(
0207         intersections,
0208         intersector.intersectSegment(s0, s1, segment[0], sDir, true), sLength);
0209   };
0210
0211   // Helper method to intersect radial full boundaries
0212   auto intersectCircle = [&](double r) -> void {
0213     auto cIntersections = intersector.intersectCircle(r, segment[0], sDir);
0214     for (const auto& intersection : cIntersections) {
0215       checkIntersection(intersections, intersection, sLength);
0216     }
0217   };
0218
0219   // Intersect phi lines
0220   if ((std::numbers::pi - hPhi) > Acts::s_epsilon) {
0221     if (sPhi < phii[0] || ePhi < phii[0]) {
0222       intersectPhiLine(phii[0]);
0223     }
0224     if (sPhi > phii[1] || ePhi > phii[1]) {
0225       intersectPhiLine(phii[1]);
0226     }
0227     // Intersect radial segments
0228     if (sR < radii[0] || eR < radii[0]) {
0229       checkIntersection(intersections,
0230                         intersector.intersectCircleSegment(
0231                             radii[0], phii[0], phii[1], segment[0], sDir),
0232                         sLength);
0233     }
0234     if (sR > radii[1] || eR > radii[1]) {
0235       checkIntersection(intersections,
0236                         intersector.intersectCircleSegment(
0237                             radii[1], phii[0], phii[1], segment[0], sDir),
0238                         sLength);
0239     }
0240   } else {
0241     // Full radial set
0242     // Intersect radial segments
0243     if (sR < radii[0] || eR < radii[0]) {
0244       intersectCircle(radii[0]);
0245     }
0246     if (sR > radii[1] || eR > radii[1]) {
0247       intersectCircle(radii[1]);
0248     }
0249   }
0250   return maskAndReturn(intersections, segment, firstInside);
0251 }
0252
0253 Acts::Result<ActsFatras::PlanarSurfaceMask::Segment2D>
0254 ActsFatras::PlanarSurfaceMask::annulusMask(const Acts::AnnulusBounds& aBounds,
0255                                            const Segment2D& segment,
0256                                            bool firstInside) const {
0257   auto vertices = aBounds.vertices(0);
0258   Acts::Vector2 moduleOrigin = aBounds.moduleOrigin();
0259
0260   std::array<std::array<unsigned int, 2>, 2> edgeCombos = {
0261       std::array<unsigned int, 2>{0, 3}, std::array<unsigned int, 2>{1, 2}};
0262
0263   std::vector<Acts::Intersection2D> intersections;
0264   Acts::Vector2 sVector(segment[1] - segment[0]);
0265   Acts::Vector2 sDir = sVector.normalized();
0266   double sLength = sVector.norm();
0267   // First the phi edges in strip system
0268   for (const auto& ec : edgeCombos) {
0269     checkIntersection(
0270         intersections,
0271         intersector.intersectSegment(vertices[ec[0]], vertices[ec[1]],
0272                                      segment[0], sDir, true),
0273         sLength);
0274   }
0275
0276   // Shift them to get the module phi and intersect
0277   std::array<unsigned int, 4> phii = {1, 0, 2, 3};
0278   for (unsigned int iarc = 0; iarc < 2; ++iarc) {
0279     Acts::Intersection2D intersection = intersector.intersectCircleSegment(
0280         aBounds.get(static_cast<Acts::AnnulusBounds::BoundValues>(iarc)),
0281         Acts::VectorHelpers::phi(vertices[phii[iarc * 2]] - moduleOrigin),
0282         Acts::VectorHelpers::phi(vertices[phii[iarc * 2 + 1]] - moduleOrigin),
0283         segment[0] - moduleOrigin, sDir);
0284     if (intersection.isValid()) {
0285       checkIntersection(intersections,
0286                         Acts::Intersection2D(
0287                             intersection.position() + moduleOrigin,
0288                             intersection.pathLength(), intersection.status()),
0289                         sLength);
0290     }
0291   }
0292   return maskAndReturn(intersections, segment, firstInside);
0293 }