Surfaces/detail/VerticesHelper.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 "Acts/Surfaces/detail/VerticesHelper.hpp"
0010
0011 #include <algorithm>
0012 #include <cmath>
0013 #include <cstddef>
0014 #include <numbers>
0015
0016 namespace Acts {
0017
0018 std::vector<double> detail::VerticesHelper::phiSegments(
0019     double phiMin, double phiMax, const std::vector<double>& phiRefs,
0020     unsigned int quarterSegments) {
0021   // Check that the phi range is valid
0022   if (phiMin > phiMax) {
0023     throw std::invalid_argument(
0024         "VerticesHelper::phiSegments ... Minimum phi must be smaller than "
0025         "maximum phi");
0026   }
0027
0028   // First check that no reference phi is outside the range
0029   for (double phiRef : phiRefs) {
0030     if (phiRef < phiMin || phiRef > phiMax) {
0031       throw std::invalid_argument(
0032           "VerticesHelper::phiSegments ... Reference phi is outside the range "
0033           "of the segment");
0034     }
0035   }
0036   if (quarterSegments == 0u) {
0037     throw std::invalid_argument(
0038         "VerticesHelper::phiSegments ... Number of segments must be larger "
0039         "than 0.");
0040   }
0041   std::vector<double> phiSegments = {phiMin, phiMax};
0042   // Minimum approximation for a circle need
0043   // - if the circle is closed the last point is given twice
0044   for (unsigned int i = 0; i < 4 * quarterSegments + 1; ++i) {
0045     double phiExt =
0046         -std::numbers::pi + i * 2 * std::numbers::pi / (4 * quarterSegments);
0047     if (phiExt > phiMin && phiExt < phiMax &&
0048         std::ranges::none_of(phiSegments, [&phiExt](double phi) {
0049           return std::abs(phi - phiExt) <
0050                  std::numeric_limits<double>::epsilon();
0051         })) {
0052       phiSegments.push_back(phiExt);
0053     }
0054   }
0055   // Add the reference phis
0056   for (const auto& phiRef : phiRefs) {
0057     if (phiRef > phiMin && phiRef < phiMax) {
0058       if (std::ranges::none_of(phiSegments, [&phiRef](double phi) {
0059             return std::abs(phi - phiRef) <
0060                    std::numeric_limits<double>::epsilon();
0061           })) {
0062         phiSegments.push_back(phiRef);
0063       }
0064     }
0065   }
0066
0067   // Sort the phis
0068   std::ranges::sort(phiSegments);
0069   return phiSegments;
0070 }
0071
0072 std::vector<Vector2> detail::VerticesHelper::ellipsoidVertices(
0073     double innerRx, double innerRy, double outerRx, double outerRy,
0074     double avgPhi, double halfPhi, unsigned int quarterSegments) {
0075   // List of vertices counter-clockwise starting at smallest phi w.r.t center,
0076   // for both inner/outer ring/segment
0077   std::vector<Vector2> rvertices;  // return vertices
0078   std::vector<Vector2> ivertices;  // inner vertices
0079   std::vector<Vector2> overtices;  // outer verices
0080
0081   bool innerExists = (innerRx > 0. && innerRy > 0.);
0082   bool closed = std::abs(halfPhi - std::numbers::pi) < s_onSurfaceTolerance;
0083
0084   std::vector<double> refPhi = {};
0085   if (avgPhi != 0.) {
0086     refPhi.push_back(avgPhi);
0087   }
0088
0089   // The inner (if exists) and outer bow
0090   if (innerExists) {
0091     ivertices = segmentVertices<Vector2, Transform2>(
0092         {innerRx, innerRy}, avgPhi - halfPhi, avgPhi + halfPhi, refPhi,
0093         quarterSegments);
0094   }
0095   overtices = segmentVertices<Vector2, Transform2>(
0096       {outerRx, outerRy}, avgPhi - halfPhi, avgPhi + halfPhi, refPhi,
0097       quarterSegments);
0098
0099   // We want to keep the same counter-clockwise orientation for displaying
0100   if (!innerExists) {
0101     if (!closed) {
0102       // Add the center case we have a sector
0103       rvertices.push_back(Vector2(0., 0.));
0104     }
0105     rvertices.insert(rvertices.end(), overtices.begin(), overtices.end());
0106   } else if (!closed) {
0107     rvertices.insert(rvertices.end(), overtices.begin(), overtices.end());
0108     rvertices.insert(rvertices.end(), ivertices.rbegin(), ivertices.rend());
0109   } else {
0110     rvertices.insert(rvertices.end(), overtices.begin(), overtices.end());
0111     rvertices.insert(rvertices.end(), ivertices.begin(), ivertices.end());
0112   }
0113   return rvertices;
0114 }
0115
0116 std::vector<Vector2> detail::VerticesHelper::circularVertices(
0117     double innerR, double outerR, double avgPhi, double halfPhi,
0118     unsigned int quarterSegments) {
0119   return ellipsoidVertices(innerR, innerR, outerR, outerR, avgPhi, halfPhi,
0120                            quarterSegments);
0121 }
0122
0123 bool detail::VerticesHelper::onHyperPlane(const std::vector<Vector3>& vertices,
0124                                           double tolerance) {
0125   // Obvious always on one surface
0126   if (vertices.size() < 4) {
0127     return true;
0128   }
0129   // Create the hyperplane
0130   auto hyperPlane = Eigen::Hyperplane<double, 3>::Through(
0131       vertices[0], vertices[1], vertices[2]);
0132   for (std::size_t ip = 3; ip < vertices.size(); ++ip) {
0133     if (hyperPlane.absDistance(vertices[ip]) > tolerance) {
0134       return false;
0135     }
0136   }
0137   return true;
0138 }
0139
0140 }  // namespace Acts