src/Surfaces/DiscTrapezoidBounds.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/DiscTrapezoidBounds.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Surfaces/detail/VerticesHelper.hpp"
0013 #include "Acts/Utilities/MathHelpers.hpp"
0014 #include "Acts/Utilities/detail/OstreamStateGuard.hpp"
0015 #include "Acts/Utilities/detail/periodic.hpp"
0016
0017 #include <cmath>
0018 #include <iomanip>
0019 #include <iostream>
0020 #include <stdexcept>
0021
0022 namespace Acts {
0023
0024 DiscTrapezoidBounds::DiscTrapezoidBounds(double halfXminR, double halfXmaxR,
0025                                          double minR, double maxR,
0026                                          double avgPhi,
0027                                          double stereo) noexcept(false)
0028     : m_values({halfXminR, halfXmaxR, minR, maxR, avgPhi, stereo}) {
0029   checkConsistency();
0030   m_ymax = fastCathetus(get(eMaxR), get(eHalfLengthXmaxR));
0031 }
0032
0033 std::vector<double> DiscTrapezoidBounds::values() const {
0034   return {m_values.begin(), m_values.end()};
0035 }
0036
0037 void DiscTrapezoidBounds::checkConsistency() noexcept(false) {
0038   if (get(eMinR) < 0. || get(eMaxR) <= 0. || get(eMinR) > get(eMaxR)) {
0039     throw std::invalid_argument("DiscTrapezoidBounds: invalid radial setup.");
0040   }
0041   if (get(eHalfLengthXminR) < 0. || get(eHalfLengthXmaxR) <= 0.) {
0042     throw std::invalid_argument("DiscTrapezoidBounds: negative length given.");
0043   }
0044   if (get(eAveragePhi) != detail::radian_sym(get(eAveragePhi))) {
0045     throw std::invalid_argument(
0046         "DiscTrapezoidBounds: invalid phi positioning.");
0047   }
0048 }
0049
0050 Vector2 DiscTrapezoidBounds::toLocalCartesian(const Vector2& lposition) const {
0051   return {lposition[0] * std::sin(lposition[1] - get(eAveragePhi)),
0052           lposition[0] * std::cos(lposition[1] - get(eAveragePhi))};
0053 }
0054
0055 SquareMatrix2 DiscTrapezoidBounds::boundToCartesianJacobian(
0056     const Vector2& lposition) const {
0057   SquareMatrix2 j;
0058   j(0, 0) = std::cos(lposition[1] - get(eAveragePhi));
0059   j(0, 1) = -lposition[0] * std::sin(lposition[1] - get(eAveragePhi));
0060   j(1, 0) = std::sin(lposition[1] - get(eAveragePhi));
0061   j(1, 1) = lposition[0] * std::cos(lposition[1] - get(eAveragePhi));
0062   return j;
0063 }
0064
0065 SquareMatrix2 DiscTrapezoidBounds::boundToCartesianMetric(
0066     const Vector2& lposition) const {
0067   SquareMatrix2 m;
0068   m(0, 0) = 1;
0069   m(0, 1) = 0;
0070   m(1, 0) = 0;
0071   m(1, 1) = lposition[0] * lposition[0];
0072   return m;
0073 }
0074
0075 bool DiscTrapezoidBounds::inside(const Vector2& lposition) const {
0076   const std::array<Vector2, 4> vertices{{{get(eHalfLengthXminR), get(eMinR)},
0077                                          {get(eHalfLengthXmaxR), m_ymax},
0078                                          {-get(eHalfLengthXmaxR), m_ymax},
0079                                          {-get(eHalfLengthXminR), get(eMinR)}}};
0080   return detail::VerticesHelper::isInsidePolygon(toLocalCartesian(lposition),
0081                                                  vertices);
0082 }
0083
0084 Vector2 DiscTrapezoidBounds::closestPoint(const Vector2& lposition,
0085                                           const SquareMatrix2& metric) const {
0086   const std::array<Vector2, 4> vertices{{{get(eHalfLengthXminR), get(eMinR)},
0087                                          {get(eHalfLengthXmaxR), m_ymax},
0088                                          {-get(eHalfLengthXmaxR), m_ymax},
0089                                          {-get(eHalfLengthXminR), get(eMinR)}}};
0090   return detail::VerticesHelper::computeClosestPointOnPolygon(
0091       toLocalCartesian(lposition), vertices, metric);
0092 }
0093
0094 std::vector<Vector2> DiscTrapezoidBounds::vertices(
0095     unsigned int /*ignoredSegments*/) const {
0096   const Vector2 cAxis(std::cos(get(eAveragePhi)), std::sin(get(eAveragePhi)));
0097   const Vector2 nAxis(cAxis.y(), -cAxis.x());
0098   const double ymin = fastCathetus(get(eMinR), get(eHalfLengthXminR));
0099   const double halfY = (m_ymax - ymin) / 2;
0100   return {-halfY * cAxis - get(eHalfLengthXminR) * nAxis,
0101           -halfY * cAxis + get(eHalfLengthXminR) * nAxis,
0102           halfY * cAxis + get(eHalfLengthXmaxR) * nAxis,
0103           halfY * cAxis - get(eHalfLengthXmaxR) * nAxis};
0104 }
0105
0106 Vector2 DiscTrapezoidBounds::center() const {
0107   // For disc trapezoid bounds in polar coordinates (r, phi),
0108   // centroid is at the middle radius and average phi
0109   const double rCentroid = 0.5 * (get(eMinR) + get(eMaxR));
0110   const double phiCentroid = get(eAveragePhi);
0111   return Vector2(rCentroid, phiCentroid);
0112 }
0113
0114 std::ostream& DiscTrapezoidBounds::toStream(std::ostream& sl) const {
0115   detail::OstreamStateGuard guard{sl};
0116   sl << std::fixed << std::setprecision(7);
0117   sl << "Acts::DiscTrapezoidBounds: (innerRadius, outerRadius, "
0118         "halfLengthXminR, "
0119         "halfLengthXmaxR, halfLengthY, halfPhiSector, averagePhi, rCenter, "
0120         "stereo) = ";
0121   sl << "(" << get(eMinR) << ", " << get(eMaxR) << ", " << get(eHalfLengthXminR)
0122      << ", " << get(eHalfLengthXmaxR) << ", " << halfLengthY() << ", "
0123      << halfPhiSector() << ", " << get(eAveragePhi) << ", " << rCenter() << ", "
0124      << stereo() << ")";
0125   return sl;
0126 }
0127
0128 }  // namespace Acts