src/Digitization/PlanarSurfaceMask.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020 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 http://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/SurfaceBounds.hpp"
0014 #include "Acts/Utilities/Intersection.hpp"
0015 #include "ActsFatras/Digitization/DigitizationError.hpp"
0016 #include <Acts/Surfaces/AnnulusBounds.hpp>
0017 #include <Acts/Surfaces/DiscTrapezoidBounds.hpp>
0018 #include <Acts/Surfaces/PlanarBounds.hpp>
0019 #include <Acts/Surfaces/RadialBounds.hpp>
0020 #include <Acts/Surfaces/Surface.hpp>
0021 #include <Acts/Utilities/Helpers.hpp>
0022
0023 #include <algorithm>
0024 #include <cmath>
0025 #include <cstddef>
0026 #include <memory>
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 && 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::sort(intersections.begin(), intersections.end(),
0061             Acts::Intersection2D::pathLengthOrder);
0062   if (intersections.size() >= 2) {
0063     return ActsFatras::PlanarSurfaceMask::Segment2D{
0064         intersections[0].position(), intersections[1].position()};
0065   } else if (intersections.size() == 1) {
0066     return (!firstInside
0067                 ? ActsFatras::PlanarSurfaceMask::Segment2D{intersections[0]
0068                                                                .position(),
0069                                                            segment[1]}
0070                 : ActsFatras::PlanarSurfaceMask::Segment2D{
0071                       segment[0], intersections[0].position()});
0072   }
0073   return ActsFatras::DigitizationError::MaskingError;
0074 }
0075
0076 }  // anonymous namespace
0077
0078 Acts::Result<ActsFatras::PlanarSurfaceMask::Segment2D>
0079 ActsFatras::PlanarSurfaceMask::apply(const Acts::Surface& surface,
0080                                      const Segment2D& segment) const {
0081   auto surfaceType = surface.type();
0082   Segment2D clipped(segment);
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::BoundaryCheck(true));
0101     bool endInside =
0102         surface.bounds().inside(localEnd, Acts::BoundaryCheck(true));
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::BoundaryCheck(true));
0136     bool endInside = surface.bounds().inside(ePolar, Acts::BoundaryCheck(true));
0137
0138     // Fast exit for both inside
0139     if (startInside && endInside) {
0140       return segment;
0141     }
0142
0143     auto boundsType = surface.bounds().type();
0144     if (boundsType == Acts::SurfaceBounds::eDisc) {
0145       auto rBounds =
0146           static_cast<const Acts::RadialBounds*>(&(surface.bounds()));
0147       return radialMask(*rBounds, segment, {sPolar, ePolar}, startInside);
0148
0149     } else if (boundsType == Acts::SurfaceBounds::eAnnulus) {
0150       auto aBounds =
0151           static_cast<const Acts::AnnulusBounds*>(&(surface.bounds()));
0152       return annulusMask(*aBounds, segment, startInside);
0153     }
0154   }
0155   return DigitizationError::UndefinedSurface;
0156 }
0157
0158 Acts::Result<ActsFatras::PlanarSurfaceMask::Segment2D>
0159 ActsFatras::PlanarSurfaceMask::polygonMask(
0160     const std::vector<Acts::Vector2>& vertices, const Segment2D& segment,
0161     bool firstInside) const {
0162   std::vector<Acts::Intersection2D> intersections;
0163   Acts::Vector2 sVector(segment[1] - segment[0]);
0164   Acts::Vector2 sDir = sVector.normalized();
0165   double sLength = sVector.norm();
0166
0167   for (std::size_t iv = 0; iv < vertices.size(); ++iv) {
0168     const Acts::Vector2& s0 = vertices[iv];
0169     const Acts::Vector2& s1 =
0170         (iv + 1) < vertices.size() ? vertices[iv + 1] : vertices[0];
0171     checkIntersection(
0172         intersections,
0173         intersector.intersectSegment(s0, s1, segment[0], sDir, true), sLength);
0174   }
0175   return maskAndReturn(intersections, segment, firstInside);
0176 }
0177
0178 Acts::Result<ActsFatras::PlanarSurfaceMask::Segment2D>
0179 ActsFatras::PlanarSurfaceMask::radialMask(const Acts::RadialBounds& rBounds,
0180                                           const Segment2D& segment,
0181                                           const Segment2D& polarSegment,
0182                                           bool firstInside) const {
0183   double rMin = rBounds.get(Acts::RadialBounds::eMinR);
0184   double rMax = rBounds.get(Acts::RadialBounds::eMaxR);
0185   double hPhi = rBounds.get(Acts::RadialBounds::eHalfPhiSector);
0186   double aPhi = rBounds.get(Acts::RadialBounds::eAveragePhi);
0187
0188   std::array<double, 2> radii = {rMin, rMax};
0189   std::array<double, 2> phii = {aPhi - hPhi, aPhi + hPhi};
0190
0191   std::vector<Acts::Intersection2D> intersections;
0192   Acts::Vector2 sVector(segment[1] - segment[0]);
0193   Acts::Vector2 sDir = sVector.normalized();
0194   double sLength = sVector.norm();
0195
0196   double sR = polarSegment[0][Acts::eBoundLoc0];
0197   double eR = polarSegment[1][Acts::eBoundLoc0];
0198   double sPhi = polarSegment[0][Acts::eBoundLoc1];
0199   double ePhi = polarSegment[1][Acts::eBoundLoc1];
0200
0201   // Helper method to intersect phi boundaries
0202   auto intersectPhiLine = [&](double phi) -> void {
0203     Acts::Vector2 s0(rMin * std::cos(phi), rMin * std::sin(phi));
0204     Acts::Vector2 s1(rMax * std::cos(phi), rMax * std::sin(phi));
0205     checkIntersection(
0206         intersections,
0207         intersector.intersectSegment(s0, s1, segment[0], sDir, true), sLength);
0208   };
0209
0210   // Helper method to intersect radial full boundaries
0211   auto intersectCircle = [&](double r) -> void {
0212     auto cIntersections = intersector.intersectCircle(r, segment[0], sDir);
0213     for (const auto& intersection : cIntersections) {
0214       checkIntersection(intersections, intersection, sLength);
0215     }
0216   };
0217
0218   // Intersect phi lines
0219   if ((M_PI - hPhi) > Acts::s_epsilon) {
0220     if (sPhi < phii[0] || ePhi < phii[0]) {
0221       intersectPhiLine(phii[0]);
0222     }
0223     if (sPhi > phii[1] || ePhi > phii[1]) {
0224       intersectPhiLine(phii[1]);
0225     }
0226     // Intersect radial segments
0227     if (sR < radii[0] || eR < radii[0]) {
0228       checkIntersection(intersections,
0229                         intersector.intersectCircleSegment(
0230                             radii[0], phii[0], phii[1], segment[0], sDir),
0231                         sLength);
0232     }
0233     if (sR > radii[1] || eR > radii[1]) {
0234       checkIntersection(intersections,
0235                         intersector.intersectCircleSegment(
0236                             radii[1], phii[0], phii[1], segment[0], sDir),
0237                         sLength);
0238     }
0239   } else {
0240     // Full radial set
0241     // Intersect radial segments
0242     if (sR < radii[0] || eR < radii[0]) {
0243       intersectCircle(radii[0]);
0244     }
0245     if (sR > radii[1] || eR > radii[1]) {
0246       intersectCircle(radii[1]);
0247     }
0248   }
0249   return maskAndReturn(intersections, segment, firstInside);
0250 }
0251
0252 Acts::Result<ActsFatras::PlanarSurfaceMask::Segment2D>
0253 ActsFatras::PlanarSurfaceMask::annulusMask(const Acts::AnnulusBounds& aBounds,
0254                                            const Segment2D& segment,
0255                                            bool firstInside) const {
0256   auto vertices = aBounds.vertices(0);
0257   Acts::Vector2 moduleOrigin = aBounds.moduleOrigin();
0258
0259   std::array<std::array<unsigned int, 2>, 2> edgeCombos = {
0260       std::array<unsigned int, 2>{0, 3}, std::array<unsigned int, 2>{1, 2}};
0261
0262   std::vector<Acts::Intersection2D> intersections;
0263   Acts::Vector2 sVector(segment[1] - segment[0]);
0264   Acts::Vector2 sDir = sVector.normalized();
0265   double sLength = sVector.norm();
0266   // First the phi edges in strip system
0267   for (const auto& ec : edgeCombos) {
0268     checkIntersection(
0269         intersections,
0270         intersector.intersectSegment(vertices[ec[0]], vertices[ec[1]],
0271                                      segment[0], sDir, true),
0272         sLength);
0273   }
0274
0275   // Shift them to get the module phi and intersect
0276   std::array<unsigned int, 4> phii = {1, 0, 2, 3};
0277   for (unsigned int iarc = 0; iarc < 2; ++iarc) {
0278     Acts::Intersection2D intersection = intersector.intersectCircleSegment(
0279         aBounds.get(static_cast<Acts::AnnulusBounds::BoundValues>(iarc)),
0280         Acts::VectorHelpers::phi(vertices[phii[iarc * 2]] - moduleOrigin),
0281         Acts::VectorHelpers::phi(vertices[phii[iarc * 2 + 1]] - moduleOrigin),
0282         segment[0] - moduleOrigin, sDir);
0283     if (intersection) {
0284       checkIntersection(intersections,
0285                         Acts::Intersection2D(
0286                             intersection.position() + moduleOrigin,
0287                             intersection.pathLength(), intersection.status()),
0288                         sLength);
0289     }
0290   }
0291   return maskAndReturn(intersections, segment, firstInside);
0292 }