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   Segment2D clipped(segment);
0084
0085   // Plane surface section -------------------
0086   if (surfaceType == Acts::Surface::Plane ||
0087       surface.bounds().type() == Acts::SurfaceBounds::eDiscTrapezoid) {
0088     Acts::Vector2 localStart =
0089         (surfaceType == Acts::Surface::Plane)
0090             ? segment[0]
0091             : Acts::Vector2(Acts::VectorHelpers::perp(segment[0]),
0092                             Acts::VectorHelpers::phi(segment[0]));
0093
0094     Acts::Vector2 localEnd =
0095         (surfaceType == Acts::Surface::Plane)
0096             ? segment[1]
0097             : Acts::Vector2(Acts::VectorHelpers::perp(segment[1]),
0098                             Acts::VectorHelpers::phi(segment[1]));
0099
0100     bool startInside =
0101         surface.bounds().inside(localStart, Acts::BoundaryTolerance::None());
0102     bool endInside =
0103         surface.bounds().inside(localEnd, Acts::BoundaryTolerance::None());
0104
0105     // Fast exit, both inside
0106     if (startInside && endInside) {
0107       return segment;
0108     }
0109
0110     // It's either planar or disc trapezoid bounds
0111     const Acts::PlanarBounds* planarBounds = nullptr;
0112     const Acts::DiscTrapezoidBounds* dtbBounds = nullptr;
0113     if (surfaceType == Acts::Surface::Plane) {
0114       planarBounds =
0115           static_cast<const Acts::PlanarBounds*>(&(surface.bounds()));
0116       if (planarBounds->type() == Acts::SurfaceBounds::eEllipse) {
0117         return DigitizationError::UndefinedSurface;
0118       }
0119     } else {
0120       dtbBounds =
0121           static_cast<const Acts::DiscTrapezoidBounds*>(&(surface.bounds()));
0122     }
0123     auto vertices = planarBounds != nullptr ? planarBounds->vertices(1)
0124                                             : dtbBounds->vertices(1);
0125
0126     return polygonMask(vertices, segment, startInside);
0127
0128   } else if (surfaceType == Acts::Surface::Disc) {
0129     // Polar coordinates
0130     Acts::Vector2 sPolar(Acts::VectorHelpers::perp(segment[0]),
0131                          Acts::VectorHelpers::phi(segment[0]));
0132     Acts::Vector2 ePolar(Acts::VectorHelpers::perp(segment[1]),
0133                          Acts::VectorHelpers::phi(segment[1]));
0134
0135     bool startInside =
0136         surface.bounds().inside(sPolar, Acts::BoundaryTolerance::None());
0137     bool endInside =
0138         surface.bounds().inside(ePolar, Acts::BoundaryTolerance::None());
0139
0140     // Fast exit for both inside
0141     if (startInside && endInside) {
0142       return segment;
0143     }
0144
0145     auto boundsType = surface.bounds().type();
0146     if (boundsType == Acts::SurfaceBounds::eDisc) {
0147       auto rBounds =
0148           static_cast<const Acts::RadialBounds*>(&(surface.bounds()));
0149       return radialMask(*rBounds, segment, {sPolar, ePolar}, startInside);
0150
0151     } else if (boundsType == Acts::SurfaceBounds::eAnnulus) {
0152       auto aBounds =
0153           static_cast<const Acts::AnnulusBounds*>(&(surface.bounds()));
0154       return annulusMask(*aBounds, segment, startInside);
0155     }
0156   }
0157   return DigitizationError::UndefinedSurface;
0158 }
0159
0160 Acts::Result<ActsFatras::PlanarSurfaceMask::Segment2D>
0161 ActsFatras::PlanarSurfaceMask::polygonMask(
0162     const std::vector<Acts::Vector2>& vertices, const Segment2D& segment,
0163     bool firstInside) const {
0164   std::vector<Acts::Intersection2D> intersections;
0165   Acts::Vector2 sVector(segment[1] - segment[0]);
0166   Acts::Vector2 sDir = sVector.normalized();
0167   double sLength = sVector.norm();
0168
0169   for (std::size_t iv = 0; iv < vertices.size(); ++iv) {
0170     const Acts::Vector2& s0 = vertices[iv];
0171     const Acts::Vector2& s1 =
0172         (iv + 1) < vertices.size() ? vertices[iv + 1] : vertices[0];
0173     checkIntersection(
0174         intersections,
0175         intersector.intersectSegment(s0, s1, segment[0], sDir, true), sLength);
0176   }
0177   return maskAndReturn(intersections, segment, firstInside);
0178 }
0179
0180 Acts::Result<ActsFatras::PlanarSurfaceMask::Segment2D>
0181 ActsFatras::PlanarSurfaceMask::radialMask(const Acts::RadialBounds& rBounds,
0182                                           const Segment2D& segment,
0183                                           const Segment2D& polarSegment,
0184                                           bool firstInside) const {
0185   double rMin = rBounds.get(Acts::RadialBounds::eMinR);
0186   double rMax = rBounds.get(Acts::RadialBounds::eMaxR);
0187   double hPhi = rBounds.get(Acts::RadialBounds::eHalfPhiSector);
0188   double aPhi = rBounds.get(Acts::RadialBounds::eAveragePhi);
0189
0190   std::array<double, 2> radii = {rMin, rMax};
0191   std::array<double, 2> phii = {aPhi - hPhi, aPhi + hPhi};
0192
0193   std::vector<Acts::Intersection2D> intersections;
0194   Acts::Vector2 sVector(segment[1] - segment[0]);
0195   Acts::Vector2 sDir = sVector.normalized();
0196   double sLength = sVector.norm();
0197
0198   double sR = polarSegment[0][Acts::eBoundLoc0];
0199   double eR = polarSegment[1][Acts::eBoundLoc0];
0200   double sPhi = polarSegment[0][Acts::eBoundLoc1];
0201   double ePhi = polarSegment[1][Acts::eBoundLoc1];
0202
0203   // Helper method to intersect phi boundaries
0204   auto intersectPhiLine = [&](double phi) -> void {
0205     Acts::Vector2 s0(rMin * std::cos(phi), rMin * std::sin(phi));
0206     Acts::Vector2 s1(rMax * std::cos(phi), rMax * std::sin(phi));
0207     checkIntersection(
0208         intersections,
0209         intersector.intersectSegment(s0, s1, segment[0], sDir, true), sLength);
0210   };
0211
0212   // Helper method to intersect radial full boundaries
0213   auto intersectCircle = [&](double r) -> void {
0214     auto cIntersections = intersector.intersectCircle(r, segment[0], sDir);
0215     for (const auto& intersection : cIntersections) {
0216       checkIntersection(intersections, intersection, sLength);
0217     }
0218   };
0219
0220   // Intersect phi lines
0221   if ((std::numbers::pi - hPhi) > Acts::s_epsilon) {
0222     if (sPhi < phii[0] || ePhi < phii[0]) {
0223       intersectPhiLine(phii[0]);
0224     }
0225     if (sPhi > phii[1] || ePhi > phii[1]) {
0226       intersectPhiLine(phii[1]);
0227     }
0228     // Intersect radial segments
0229     if (sR < radii[0] || eR < radii[0]) {
0230       checkIntersection(intersections,
0231                         intersector.intersectCircleSegment(
0232                             radii[0], phii[0], phii[1], segment[0], sDir),
0233                         sLength);
0234     }
0235     if (sR > radii[1] || eR > radii[1]) {
0236       checkIntersection(intersections,
0237                         intersector.intersectCircleSegment(
0238                             radii[1], phii[0], phii[1], segment[0], sDir),
0239                         sLength);
0240     }
0241   } else {
0242     // Full radial set
0243     // Intersect radial segments
0244     if (sR < radii[0] || eR < radii[0]) {
0245       intersectCircle(radii[0]);
0246     }
0247     if (sR > radii[1] || eR > radii[1]) {
0248       intersectCircle(radii[1]);
0249     }
0250   }
0251   return maskAndReturn(intersections, segment, firstInside);
0252 }
0253
0254 Acts::Result<ActsFatras::PlanarSurfaceMask::Segment2D>
0255 ActsFatras::PlanarSurfaceMask::annulusMask(const Acts::AnnulusBounds& aBounds,
0256                                            const Segment2D& segment,
0257                                            bool firstInside) const {
0258   auto vertices = aBounds.vertices(0);
0259   Acts::Vector2 moduleOrigin = aBounds.moduleOrigin();
0260
0261   std::array<std::array<unsigned int, 2>, 2> edgeCombos = {
0262       std::array<unsigned int, 2>{0, 3}, std::array<unsigned int, 2>{1, 2}};
0263
0264   std::vector<Acts::Intersection2D> intersections;
0265   Acts::Vector2 sVector(segment[1] - segment[0]);
0266   Acts::Vector2 sDir = sVector.normalized();
0267   double sLength = sVector.norm();
0268   // First the phi edges in strip system
0269   for (const auto& ec : edgeCombos) {
0270     checkIntersection(
0271         intersections,
0272         intersector.intersectSegment(vertices[ec[0]], vertices[ec[1]],
0273                                      segment[0], sDir, true),
0274         sLength);
0275   }
0276
0277   // Shift them to get the module phi and intersect
0278   std::array<unsigned int, 4> phii = {1, 0, 2, 3};
0279   for (unsigned int iarc = 0; iarc < 2; ++iarc) {
0280     Acts::Intersection2D intersection = intersector.intersectCircleSegment(
0281         aBounds.get(static_cast<Acts::AnnulusBounds::BoundValues>(iarc)),
0282         Acts::VectorHelpers::phi(vertices[phii[iarc * 2]] - moduleOrigin),
0283         Acts::VectorHelpers::phi(vertices[phii[iarc * 2 + 1]] - moduleOrigin),
0284         segment[0] - moduleOrigin, sDir);
0285     if (intersection.isValid()) {
0286       checkIntersection(intersections,
0287                         Acts::Intersection2D(
0288                             intersection.position() + moduleOrigin,
0289                             intersection.pathLength(), intersection.status()),
0290                         sLength);
0291     }
0292   }
0293   return maskAndReturn(intersections, segment, firstInside);
0294 }