src/Digitization/SurfaceMask.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/SurfaceMask.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/CylinderBounds.hpp"
0016 #include "Acts/Surfaces/DiscTrapezoidBounds.hpp"
0017 #include "Acts/Surfaces/PlanarBounds.hpp"
0018 #include "Acts/Surfaces/RadialBounds.hpp"
0019 #include "Acts/Surfaces/Surface.hpp"
0020 #include "Acts/Surfaces/SurfaceBounds.hpp"
0021 #include "Acts/Utilities/Intersection.hpp"
0022 #include "ActsFatras/Digitization/DigitizationError.hpp"
0023
0024 #include <algorithm>
0025 #include <array>
0026 #include <cmath>
0027 #include <cstddef>
0028 #include <limits>
0029 #include <numbers>
0030
0031 namespace ActsFatras {
0032 namespace {
0033
0034 /// Helper method to check if an intersection is good.
0035 ///
0036 /// Good in this context is defined as: along direction,
0037 /// closer than the segment length & reachable
0038 ///
0039 /// @param intersections The confirmed intersections for the segment
0040 /// @param candidate The candidate intersection
0041 /// @param sLength The segment length, maximal allowed length
0042 void checkIntersection(std::vector<Acts::Intersection2D>& intersections,
0043                        const Acts::Intersection2D& candidate, double sLength) {
0044   if (candidate.isValid() && candidate.pathLength() > 0 &&
0045       candidate.pathLength() < sLength) {
0046     intersections.push_back(candidate);
0047   }
0048 }
0049
0050 /// Helper method to apply the mask and return.
0051 ///
0052 /// If two (or more) intersections would be good, apply the first two
0053 /// If only one is available, the boolean tells you which one it is.
0054 /// If no intersection is valid, return an error code for masking.
0055 ///
0056 /// @param intersections All confirmed intersections
0057 /// @param segment The original segment before masking
0058 /// @param firstInside Indicator if the first is inside or not
0059 ///
0060 /// @return a new Segment (clipped) wrapped in a result or error_code
0061 Acts::Result<SurfaceMask::Segment2D> maskAndReturn(
0062     std::vector<Acts::Intersection2D>& intersections,
0063     const SurfaceMask::Segment2D& segment, bool firstInside) {
0064   std::ranges::sort(intersections, Acts::Intersection2D::pathLengthOrder);
0065   if (intersections.size() >= 2) {
0066     return SurfaceMask::Segment2D{intersections[0].position(),
0067                                   intersections[1].position()};
0068   } else if (intersections.size() == 1) {
0069     return (
0070         !firstInside
0071             ? SurfaceMask::Segment2D{intersections[0].position(), segment[1]}
0072             : SurfaceMask::Segment2D{segment[0], intersections[0].position()});
0073   }
0074   return DigitizationError::MaskingError;
0075 }
0076
0077 /// Liang–Barsky clipping of a segment against an axis-aligned rectangle.
0078 ///
0079 /// @return on success the {tEnter, tExit} parameter pair in [0, 1] describing
0080 /// the clipped segment along (end - start); a DigitizationError::MaskingError
0081 /// if the segment lies completely outside the rectangle.
0082 Acts::Result<std::array<double, 2>> clipLiangBarsky(double x0, double y0,
0083                                                     double x1, double y1,
0084                                                     double xmin, double xmax,
0085                                                     double ymin, double ymax) {
0086   const double dx = x1 - x0;
0087   const double dy = y1 - y0;
0088
0089   double tEnter = 0.0;
0090   double tExit = 1.0;
0091
0092   const std::array<double, 4> p = {-dx, dx, -dy, dy};
0093   const std::array<double, 4> q = {x0 - xmin, xmax - x0, y0 - ymin, ymax - y0};
0094
0095   for (std::size_t i = 0; i < p.size(); ++i) {
0096     if (p[i] == 0.0) {
0097       // Segment is parallel to this clipping edge.
0098       if (q[i] < 0.0) {
0099         return DigitizationError::MaskingError;
0100       }
0101       continue;
0102     }
0103     const double t = q[i] / p[i];
0104     if (p[i] < 0.0) {
0105       tEnter = std::max(tEnter, t);
0106     } else {
0107       tExit = std::min(tExit, t);
0108     }
0109   }
0110
0111   if (tEnter > tExit) {
0112     return DigitizationError::MaskingError;
0113   }
0114   return std::array<double, 2>{tEnter, tExit};
0115 }
0116
0117 }  // anonymous namespace
0118
0119 Acts::Result<SurfaceMask::Segment2D> SurfaceMask::apply(
0120     const Acts::Surface& surface, const Segment2D& segment) const {
0121   auto surfaceType = surface.type();
0122
0123   // Cylinder surface section -------------------
0124   if (surfaceType == Acts::Surface::Cylinder &&
0125       surface.bounds().type() == Acts::SurfaceBounds::eCylinder) {
0126     const auto& cBounds =
0127         static_cast<const Acts::CylinderBounds&>(surface.bounds());
0128     return cylinderMask(cBounds, segment);
0129   }
0130
0131   // Plane surface section -------------------
0132   if (surfaceType == Acts::Surface::Plane ||
0133       surface.bounds().type() == Acts::SurfaceBounds::eDiscTrapezoid) {
0134     Acts::Vector2 localStart =
0135         (surfaceType == Acts::Surface::Plane)
0136             ? segment[0]
0137             : Acts::Vector2(Acts::VectorHelpers::perp(segment[0]),
0138                             Acts::VectorHelpers::phi(segment[0]));
0139
0140     Acts::Vector2 localEnd =
0141         (surfaceType == Acts::Surface::Plane)
0142             ? segment[1]
0143             : Acts::Vector2(Acts::VectorHelpers::perp(segment[1]),
0144                             Acts::VectorHelpers::phi(segment[1]));
0145
0146     bool startInside =
0147         surface.bounds().inside(localStart, Acts::BoundaryTolerance::None());
0148     bool endInside =
0149         surface.bounds().inside(localEnd, Acts::BoundaryTolerance::None());
0150
0151     // Fast exit, both inside
0152     if (startInside && endInside) {
0153       return segment;
0154     }
0155
0156     // It's either planar or disc trapezoid bounds
0157     const Acts::PlanarBounds* planarBounds = nullptr;
0158     const Acts::DiscTrapezoidBounds* dtbBounds = nullptr;
0159     if (surfaceType == Acts::Surface::Plane) {
0160       planarBounds =
0161           static_cast<const Acts::PlanarBounds*>(&(surface.bounds()));
0162       if (planarBounds->type() == Acts::SurfaceBounds::eEllipse) {
0163         return DigitizationError::UndefinedSurface;
0164       }
0165     } else {
0166       dtbBounds =
0167           static_cast<const Acts::DiscTrapezoidBounds*>(&(surface.bounds()));
0168     }
0169     auto vertices = planarBounds != nullptr ? planarBounds->vertices(1)
0170                                             : dtbBounds->vertices(1);
0171
0172     return polygonMask(vertices, segment, startInside);
0173
0174   } else if (surfaceType == Acts::Surface::Disc) {
0175     // Polar coordinates
0176     Acts::Vector2 sPolar(Acts::VectorHelpers::perp(segment[0]),
0177                          Acts::VectorHelpers::phi(segment[0]));
0178     Acts::Vector2 ePolar(Acts::VectorHelpers::perp(segment[1]),
0179                          Acts::VectorHelpers::phi(segment[1]));
0180
0181     bool startInside =
0182         surface.bounds().inside(sPolar, Acts::BoundaryTolerance::None());
0183     bool endInside =
0184         surface.bounds().inside(ePolar, Acts::BoundaryTolerance::None());
0185
0186     // Fast exit for both inside
0187     if (startInside && endInside) {
0188       return segment;
0189     }
0190
0191     auto boundsType = surface.bounds().type();
0192     if (boundsType == Acts::SurfaceBounds::eDisc) {
0193       auto rBounds =
0194           static_cast<const Acts::RadialBounds*>(&(surface.bounds()));
0195       return radialMask(*rBounds, segment, {sPolar, ePolar}, startInside);
0196
0197     } else if (boundsType == Acts::SurfaceBounds::eAnnulus) {
0198       auto aBounds =
0199           static_cast<const Acts::AnnulusBounds*>(&(surface.bounds()));
0200       return annulusMask(*aBounds, segment, startInside);
0201     }
0202   }
0203   return DigitizationError::UndefinedSurface;
0204 }
0205
0206 Acts::Result<SurfaceMask::Segment2D> SurfaceMask::polygonMask(
0207     const std::vector<Acts::Vector2>& vertices, const Segment2D& segment,
0208     bool firstInside) const {
0209   std::vector<Acts::Intersection2D> intersections;
0210   Acts::Vector2 sVector(segment[1] - segment[0]);
0211   Acts::Vector2 sDir = sVector.normalized();
0212   double sLength = sVector.norm();
0213
0214   for (std::size_t iv = 0; iv < vertices.size(); ++iv) {
0215     const Acts::Vector2& s0 = vertices[iv];
0216     const Acts::Vector2& s1 =
0217         (iv + 1) < vertices.size() ? vertices[iv + 1] : vertices[0];
0218     checkIntersection(
0219         intersections,
0220         intersector.intersectSegment(s0, s1, segment[0], sDir, true), sLength);
0221   }
0222   return maskAndReturn(intersections, segment, firstInside);
0223 }
0224
0225 Acts::Result<SurfaceMask::Segment2D> SurfaceMask::radialMask(
0226     const Acts::RadialBounds& rBounds, const Segment2D& segment,
0227     const Segment2D& polarSegment, bool firstInside) const {
0228   double rMin = rBounds.get(Acts::RadialBounds::eMinR);
0229   double rMax = rBounds.get(Acts::RadialBounds::eMaxR);
0230   double hPhi = rBounds.get(Acts::RadialBounds::eHalfPhiSector);
0231   double aPhi = rBounds.get(Acts::RadialBounds::eAveragePhi);
0232
0233   std::array<double, 2> radii = {rMin, rMax};
0234   std::array<double, 2> phii = {aPhi - hPhi, aPhi + hPhi};
0235
0236   std::vector<Acts::Intersection2D> intersections;
0237   Acts::Vector2 sVector(segment[1] - segment[0]);
0238   Acts::Vector2 sDir = sVector.normalized();
0239   double sLength = sVector.norm();
0240
0241   double sR = polarSegment[0][Acts::eBoundLoc0];
0242   double eR = polarSegment[1][Acts::eBoundLoc0];
0243   double sPhi = polarSegment[0][Acts::eBoundLoc1];
0244   double ePhi = polarSegment[1][Acts::eBoundLoc1];
0245
0246   // Helper method to intersect phi boundaries
0247   auto intersectPhiLine = [&](double phi) -> void {
0248     Acts::Vector2 s0(rMin * std::cos(phi), rMin * std::sin(phi));
0249     Acts::Vector2 s1(rMax * std::cos(phi), rMax * std::sin(phi));
0250     checkIntersection(
0251         intersections,
0252         intersector.intersectSegment(s0, s1, segment[0], sDir, true), sLength);
0253   };
0254
0255   // Helper method to intersect radial full boundaries
0256   auto intersectCircle = [&](double r) -> void {
0257     auto cIntersections = intersector.intersectCircle(r, segment[0], sDir);
0258     for (const auto& intersection : cIntersections) {
0259       checkIntersection(intersections, intersection, sLength);
0260     }
0261   };
0262
0263   // Intersect phi lines
0264   if ((std::numbers::pi - hPhi) > Acts::s_epsilon) {
0265     if (sPhi < phii[0] || ePhi < phii[0]) {
0266       intersectPhiLine(phii[0]);
0267     }
0268     if (sPhi > phii[1] || ePhi > phii[1]) {
0269       intersectPhiLine(phii[1]);
0270     }
0271     // Intersect radial segments
0272     if (sR < radii[0] || eR < radii[0]) {
0273       checkIntersection(intersections,
0274                         intersector.intersectCircleSegment(
0275                             radii[0], phii[0], phii[1], segment[0], sDir),
0276                         sLength);
0277     }
0278     if (sR > radii[1] || eR > radii[1]) {
0279       checkIntersection(intersections,
0280                         intersector.intersectCircleSegment(
0281                             radii[1], phii[0], phii[1], segment[0], sDir),
0282                         sLength);
0283     }
0284   } else {
0285     // Full radial set
0286     // Intersect radial segments
0287     if (sR < radii[0] || eR < radii[0]) {
0288       intersectCircle(radii[0]);
0289     }
0290     if (sR > radii[1] || eR > radii[1]) {
0291       intersectCircle(radii[1]);
0292     }
0293   }
0294   return maskAndReturn(intersections, segment, firstInside);
0295 }
0296
0297 Acts::Result<SurfaceMask::Segment2D> SurfaceMask::annulusMask(
0298     const Acts::AnnulusBounds& aBounds, const Segment2D& segment,
0299     bool firstInside) const {
0300   auto vertices = aBounds.vertices(0);
0301   Acts::Vector2 moduleOrigin = aBounds.moduleOrigin();
0302
0303   std::array<std::array<unsigned int, 2>, 2> edgeCombos = {
0304       std::array<unsigned int, 2>{0, 3}, std::array<unsigned int, 2>{1, 2}};
0305
0306   std::vector<Acts::Intersection2D> intersections;
0307   Acts::Vector2 sVector(segment[1] - segment[0]);
0308   Acts::Vector2 sDir = sVector.normalized();
0309   double sLength = sVector.norm();
0310   // First the phi edges in strip system
0311   for (const auto& ec : edgeCombos) {
0312     checkIntersection(
0313         intersections,
0314         intersector.intersectSegment(vertices[ec[0]], vertices[ec[1]],
0315                                      segment[0], sDir, true),
0316         sLength);
0317   }
0318
0319   // Shift them to get the module phi and intersect
0320   std::array<unsigned int, 4> phii = {1, 0, 2, 3};
0321   for (unsigned int iarc = 0; iarc < 2; ++iarc) {
0322     Acts::Intersection2D intersection = intersector.intersectCircleSegment(
0323         aBounds.get(static_cast<Acts::AnnulusBounds::BoundValues>(iarc)),
0324         Acts::VectorHelpers::phi(vertices[phii[iarc * 2]] - moduleOrigin),
0325         Acts::VectorHelpers::phi(vertices[phii[iarc * 2 + 1]] - moduleOrigin),
0326         segment[0] - moduleOrigin, sDir);
0327     if (intersection.isValid()) {
0328       checkIntersection(intersections,
0329                         Acts::Intersection2D(
0330                             intersection.position() + moduleOrigin,
0331                             intersection.pathLength(), intersection.status()),
0332                         sLength);
0333     }
0334   }
0335   return maskAndReturn(intersections, segment, firstInside);
0336 }
0337
0338 Acts::Result<SurfaceMask::Segment2D> SurfaceMask::cylinderMask(
0339     const Acts::CylinderBounds& cBounds, const Segment2D& segment) const {
0340   const double R = cBounds.get(Acts::CylinderBounds::eR);
0341   const double halfZ = cBounds.get(Acts::CylinderBounds::eHalfLengthZ);
0342   const double halfPhi = cBounds.get(Acts::CylinderBounds::eHalfPhiSector);
0343   const double avgPhi = cBounds.get(Acts::CylinderBounds::eAveragePhi);
0344
0345   // (rPhi, z) bounding box of the cylinder bounds.
0346   const double rPhiMin = R * (avgPhi - halfPhi);
0347   const double rPhiMax = R * (avgPhi + halfPhi);
0348   const double zMin = -halfZ;
0349   const double zMax = halfZ;
0350
0351   // Local rPhi may be reported in the (-π R, π R] principal branch by the
0352   // drift step. If the sector is closed (full cylinder) any rPhi is valid and
0353   // we only need to clip in z. Otherwise unwrap by adding/subtracting 2π·R so
0354   // that the segment lies in the same branch as the sector midpoint.
0355   Segment2D seg = segment;
0356   const bool isClosed = std::abs(halfPhi - std::numbers::pi) < Acts::s_epsilon;
0357
0358   if (!isClosed) {
0359     const double midRPhi = R * avgPhi;
0360     const double twoPiR = 2.0 * std::numbers::pi * R;
0361     for (auto& p : seg) {
0362       while (p[0] - midRPhi > std::numbers::pi * R) {
0363         p[0] -= twoPiR;
0364       }
0365       while (midRPhi - p[0] > std::numbers::pi * R) {
0366         p[0] += twoPiR;
0367       }
0368     }
0369   }
0370
0371   // Fast exit: both endpoints inside.
0372   auto inside = [&](const Acts::Vector2& p) {
0373     const bool inZ = (p[1] >= zMin) && (p[1] <= zMax);
0374     if (isClosed) {
0375       return inZ;
0376     }
0377     return inZ && (p[0] >= rPhiMin) && (p[0] <= rPhiMax);
0378   };
0379   if (inside(seg[0]) && inside(seg[1])) {
0380     return seg;
0381   }
0382
0383   // If the cylinder is closed in phi, the only clipping needed is in z.
0384   const double useRPhiMin =
0385       isClosed ? -std::numeric_limits<double>::infinity() : rPhiMin;
0386   const double useRPhiMax =
0387       isClosed ? std::numeric_limits<double>::infinity() : rPhiMax;
0388
0389   const auto clip = clipLiangBarsky(seg[0][0], seg[0][1], seg[1][0], seg[1][1],
0390                                     useRPhiMin, useRPhiMax, zMin, zMax);
0391   if (!clip.ok()) {
0392     return clip.error();
0393   }
0394   const auto& [tEnter, tExit] = *clip;
0395
0396   const Acts::Vector2 d = seg[1] - seg[0];
0397   return Segment2D{seg[0] + tEnter * d, seg[0] + tExit * d};
0398 }
0399 }  // namespace ActsFatras