src/Surfaces/VerticesHelper.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 "Acts/Surfaces/detail/VerticesHelper.hpp"
0010
0011 #include <cmath>
0012 #include <cstddef>
0013
0014 std::vector<Acts::ActsScalar> Acts::detail::VerticesHelper::phiSegments(
0015     ActsScalar phiMin, ActsScalar phiMax,
0016     const std::vector<ActsScalar>& phiRefs, ActsScalar phiTolerance) {
0017   // This is to ensure that the extrema are built regardless of number
0018   // of segments
0019   std::vector<ActsScalar> phiSegments;
0020   std::vector<ActsScalar> quarters = {-M_PI, -0.5 * M_PI, 0., 0.5 * M_PI, M_PI};
0021   // It does not cover the full azimuth
0022   if (phiMin != -M_PI || phiMax != M_PI) {
0023     phiSegments.push_back(phiMin);
0024     for (unsigned int iq = 1; iq < 4; ++iq) {
0025       if (phiMin < quarters[iq] && phiMax > quarters[iq]) {
0026         phiSegments.push_back(quarters[iq]);
0027       }
0028     }
0029     phiSegments.push_back(phiMax);
0030   } else {
0031     phiSegments = quarters;
0032   }
0033   // Insert the reference phis if
0034   if (!phiRefs.empty()) {
0035     for (const auto& phiRef : phiRefs) {
0036       // Trying to find the right patch
0037       auto match = std::find_if(
0038           phiSegments.begin(), phiSegments.end(), [&](ActsScalar phiSeg) {
0039             return std::abs(phiSeg - phiRef) < phiTolerance;
0040           });
0041       if (match == phiSegments.end()) {
0042         phiSegments.push_back(phiRef);
0043       }
0044     }
0045     std::sort(phiSegments.begin(), phiSegments.end());
0046   }
0047   return phiSegments;
0048 }
0049
0050 std::vector<Acts::Vector2> Acts::detail::VerticesHelper::ellipsoidVertices(
0051     ActsScalar innerRx, ActsScalar innerRy, ActsScalar outerRx,
0052     ActsScalar outerRy, ActsScalar avgPhi, ActsScalar halfPhi,
0053     unsigned int lseg) {
0054   // List of vertices counter-clockwise starting at smallest phi w.r.t center,
0055   // for both inner/outer ring/segment
0056   std::vector<Vector2> rvertices;  // return vertices
0057   std::vector<Vector2> ivertices;  // inner vertices
0058   std::vector<Vector2> overtices;  // outer verices
0059
0060   bool innerExists = (innerRx > 0. && innerRy > 0.);
0061   bool closed = std::abs(halfPhi - M_PI) < s_onSurfaceTolerance;
0062
0063   // Get the phi segments from the helper method
0064   auto phiSegs = detail::VerticesHelper::phiSegments(
0065       avgPhi - halfPhi, avgPhi + halfPhi, {avgPhi});
0066
0067   // The inner (if exists) and outer bow
0068   for (unsigned int iseg = 0; iseg < phiSegs.size() - 1; ++iseg) {
0069     int addon = (iseg == phiSegs.size() - 2 && !closed) ? 1 : 0;
0070     if (innerExists) {
0071       createSegment<Vector2, Transform2>(ivertices, {innerRx, innerRy},
0072                                          phiSegs[iseg], phiSegs[iseg + 1], lseg,
0073                                          addon);
0074     }
0075     createSegment<Vector2, Transform2>(overtices, {outerRx, outerRy},
0076                                        phiSegs[iseg], phiSegs[iseg + 1], lseg,
0077                                        addon);
0078   }
0079
0080   // We want to keep the same counter-clockwise orientation for displaying
0081   if (!innerExists) {
0082     if (!closed) {
0083       // Add the center case we have a sector
0084       rvertices.push_back(Vector2(0., 0.));
0085     }
0086     rvertices.insert(rvertices.end(), overtices.begin(), overtices.end());
0087   } else if (!closed) {
0088     rvertices.insert(rvertices.end(), overtices.begin(), overtices.end());
0089     rvertices.insert(rvertices.end(), ivertices.rbegin(), ivertices.rend());
0090   } else {
0091     rvertices.insert(rvertices.end(), overtices.begin(), overtices.end());
0092     rvertices.insert(rvertices.end(), ivertices.begin(), ivertices.end());
0093   }
0094   return rvertices;
0095 }
0096
0097 std::vector<Acts::Vector2> Acts::detail::VerticesHelper::circularVertices(
0098     ActsScalar innerR, ActsScalar outerR, ActsScalar avgPhi, ActsScalar halfPhi,
0099     unsigned int lseg) {
0100   return ellipsoidVertices(innerR, innerR, outerR, outerR, avgPhi, halfPhi,
0101                            lseg);
0102 }
0103
0104 bool Acts::detail::VerticesHelper::onHyperPlane(
0105     const std::vector<Acts::Vector3>& vertices, ActsScalar tolerance) {
0106   // Obvious always on one surface
0107   if (vertices.size() < 4) {
0108     return true;
0109   }
0110   // Create the hyperplane
0111   auto hyperPlane = Eigen::Hyperplane<ActsScalar, 3>::Through(
0112       vertices[0], vertices[1], vertices[2]);
0113   for (std::size_t ip = 3; ip < vertices.size(); ++ip) {
0114     if (hyperPlane.absDistance(vertices[ip]) > tolerance) {
0115       return false;
0116     }
0117   }
0118   return true;
0119 }