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/detail/periodic.hpp"
0014
0015 #include <cmath>
0016 #include <iomanip>
0017 #include <iostream>
0018 #include <stdexcept>
0019
0020 namespace Acts {
0021
0022 DiscTrapezoidBounds::DiscTrapezoidBounds(double halfXminR, double halfXmaxR,
0023                                          double minR, double maxR,
0024                                          double avgPhi,
0025                                          double stereo) noexcept(false)
0026     : m_values({halfXminR, halfXmaxR, minR, maxR, avgPhi, stereo}) {
0027   checkConsistency();
0028   m_ymax = std::sqrt(get(eMaxR) * get(eMaxR) -
0029                      get(eHalfLengthXmaxR) * 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   Vector2 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   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   Vector2 cAxis(std::cos(get(eAveragePhi)), std::sin(get(eAveragePhi)));
0096   Vector2 nAxis(cAxis.y(), -cAxis.x());
0097   auto ymin = std::sqrt(get(eMinR) * get(eMinR) -
0098                         get(eHalfLengthXminR) * get(eHalfLengthXminR));
0099   auto 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   double rCentroid = 0.5 * (get(eMinR) + get(eMaxR));
0110   double phiCentroid = get(eAveragePhi);
0111   return Vector2(rCentroid, phiCentroid);
0112 }
0113
0114 std::ostream& DiscTrapezoidBounds::toStream(std::ostream& sl) const {
0115   sl << std::setiosflags(std::ios::fixed);
0116   sl << 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   sl << std::setprecision(-1);
0126   return sl;
0127 }
0128
0129 }  // namespace Acts