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