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