src/Surfaces/CylinderBounds.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2016-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/CylinderBounds.hpp"
0010
0011 #include "Acts/Definitions/TrackParametrization.hpp"
0012 #include "Acts/Surfaces/detail/VerticesHelper.hpp"
0013 #include "Acts/Utilities/VectorHelpers.hpp"
0014
0015 #include <algorithm>
0016 #include <cmath>
0017 #include <iomanip>
0018 #include <iostream>
0019 #include <utility>
0020
0021 using Acts::VectorHelpers::perp;
0022 using Acts::VectorHelpers::phi;
0023
0024 Acts::SurfaceBounds::BoundsType Acts::CylinderBounds::type() const {
0025   return SurfaceBounds::eCylinder;
0026 }
0027
0028 Acts::Vector2 Acts::CylinderBounds::shifted(
0029     const Acts::Vector2& lposition) const {
0030   return {Acts::detail::radian_sym((lposition[Acts::eBoundLoc0] / get(eR)) -
0031                                    get(eAveragePhi)),
0032           lposition[Acts::eBoundLoc1]};
0033 }
0034
0035 Acts::ActsMatrix<2, 2> Acts::CylinderBounds::jacobian() const {
0036   ActsMatrix<2, 2> j;
0037   j(0, eBoundLoc0) = 1 / get(eR);
0038   j(0, eBoundLoc1) = 0;
0039   j(1, eBoundLoc0) = 0;
0040   j(1, eBoundLoc1) = 1;
0041   return j;
0042 }
0043
0044 bool Acts::CylinderBounds::inside(const Vector2& lposition,
0045                                   const BoundaryCheck& bcheck) const {
0046   double bevelMinZ = get(eBevelMinZ);
0047   double bevelMaxZ = get(eBevelMaxZ);
0048
0049   double halfLengthZ = get(eHalfLengthZ);
0050   double halfPhi = get(eHalfPhiSector);
0051   if (bevelMinZ != 0. && bevelMaxZ != 0.) {
0052     double radius = get(eR);
0053     // Beleved sides will unwrap to a trapezoid
0054     ///////////////////////////////////
0055     //  ________
0056     // /| .  . |\ r/phi
0057     // \|______|/ r/phi
0058     // -Z   0  Z
0059     ///////////////////////////////////
0060     float localx =
0061         lposition[0] > radius ? 2 * radius - lposition[0] : lposition[0];
0062     Vector2 shiftedlposition = shifted(lposition);
0063     if ((std::fabs(shiftedlposition[0]) <= halfPhi &&
0064          std::fabs(shiftedlposition[1]) <= halfLengthZ)) {
0065       return true;
0066     } else if ((lposition[1] >=
0067                 -(localx * std::tan(bevelMinZ) + halfLengthZ)) &&
0068                (lposition[1] <= (localx * std::tan(bevelMaxZ) + halfLengthZ))) {
0069       return true;
0070     } else {
0071       // check within tolerance
0072       auto boundaryCheck = bcheck.transformed(jacobian());
0073
0074       Vector2 lowerLeft = {-radius, -halfLengthZ};
0075       Vector2 middleLeft = {0., -(halfLengthZ + radius * std::tan(bevelMinZ))};
0076       Vector2 upperLeft = {radius, -halfLengthZ};
0077       Vector2 upperRight = {radius, halfLengthZ};
0078       Vector2 middleRight = {0., (halfLengthZ + radius * std::tan(bevelMaxZ))};
0079       Vector2 lowerRight = {-radius, halfLengthZ};
0080       Vector2 vertices[] = {lowerLeft,  middleLeft,  upperLeft,
0081                             upperRight, middleRight, lowerRight};
0082       Vector2 closestPoint =
0083           boundaryCheck.computeClosestPointOnPolygon(lposition, vertices);
0084
0085       return boundaryCheck.isTolerated(closestPoint - lposition);
0086     }
0087   } else {
0088     return bcheck.transformed(jacobian())
0089         .isInside(shifted(lposition), Vector2(-halfPhi, -halfLengthZ),
0090                   Vector2(halfPhi, halfLengthZ));
0091   }
0092 }
0093
0094 bool Acts::CylinderBounds::inside3D(const Vector3& position,
0095                                     const BoundaryCheck& bcheck) const {
0096   // additional tolerance from the boundary check if configured
0097   bool checkAbsolute = bcheck.m_type == BoundaryCheck::Type::eAbsolute;
0098
0099   // this fast check only applies to closed cylindrical bounds
0100   double addToleranceR =
0101       (checkAbsolute && m_closed) ? bcheck.m_tolerance[0] : 0.;
0102   double addToleranceZ = checkAbsolute ? bcheck.m_tolerance[1] : 0.;
0103   // check if the position compatible with the radius
0104   if ((s_onSurfaceTolerance + addToleranceR) <=
0105       std::abs(perp(position) - get(eR))) {
0106     return false;
0107   } else if (checkAbsolute && m_closed) {
0108     double bevelMinZ = get(eBevelMinZ);
0109     double bevelMaxZ = get(eBevelMaxZ);
0110
0111     double addedMinZ =
0112         bevelMinZ != 0. ? position.y() * std::sin(bevelMinZ) : 0.;
0113     double addedMaxZ =
0114         bevelMinZ != 0. ? position.y() * std::sin(bevelMaxZ) : 0.;
0115
0116     return ((s_onSurfaceTolerance + addToleranceZ + get(eHalfLengthZ) +
0117              addedMinZ) >= position.z()) &&
0118            ((s_onSurfaceTolerance + addToleranceZ + get(eHalfLengthZ) +
0119              addedMaxZ) <= position.z());
0120   }
0121   // detailed, but slower check
0122   Vector2 lpos(detail::radian_sym(phi(position) - get(eAveragePhi)),
0123                position.z());
0124   return bcheck.transformed(jacobian())
0125       .isInside(lpos, Vector2(-get(eHalfPhiSector), -get(eHalfLengthZ)),
0126                 Vector2(get(eHalfPhiSector), get(eHalfLengthZ)));
0127 }
0128
0129 std::ostream& Acts::CylinderBounds::toStream(std::ostream& sl) const {
0130   sl << std::setiosflags(std::ios::fixed);
0131   sl << std::setprecision(7);
0132   sl << "Acts::CylinderBounds: (radius, halfLengthZ, halfPhiSector, "
0133         "averagePhi, bevelMinZ, bevelMaxZ) = ";
0134   sl << "(" << get(eR) << ", " << get(eHalfLengthZ) << ", ";
0135   sl << get(eHalfPhiSector) << ", " << get(eAveragePhi) << ", ";
0136   sl << get(eBevelMinZ) << ", " << get(eBevelMaxZ) << ")";
0137   sl << std::setprecision(-1);
0138   return sl;
0139 }
0140
0141 std::vector<Acts::Vector3> Acts::CylinderBounds::createCircles(
0142     const Transform3 ctrans, std::size_t lseg) const {
0143   std::vector<Vector3> vertices;
0144
0145   double avgPhi = get(eAveragePhi);
0146   double halfPhi = get(eHalfPhiSector);
0147
0148   bool fullCylinder = coversFullAzimuth();
0149
0150   // Get the phi segments from the helper - ensures extra points
0151   auto phiSegs = fullCylinder ? detail::VerticesHelper::phiSegments()
0152                               : detail::VerticesHelper::phiSegments(
0153                                     avgPhi - halfPhi, avgPhi + halfPhi,
0154                                     {static_cast<ActsScalar>(avgPhi)});
0155
0156   // Write the two bows/circles on either side
0157   std::vector<int> sides = {-1, 1};
0158   for (auto& side : sides) {
0159     for (std::size_t iseg = 0; iseg < phiSegs.size() - 1; ++iseg) {
0160       int addon = (iseg == phiSegs.size() - 2 && !fullCylinder) ? 1 : 0;
0161       /// Helper method to create the segment
0162       detail::VerticesHelper::createSegment(
0163           vertices, {get(eR), get(eR)}, phiSegs[iseg], phiSegs[iseg + 1], lseg,
0164           addon, Vector3(0., 0., side * get(eHalfLengthZ)), ctrans);
0165     }
0166   }
0167
0168   double bevelMinZ = get(eBevelMinZ);
0169   double bevelMaxZ = get(eBevelMaxZ);
0170
0171   // Modify the vertices position if bevel is defined
0172   if ((bevelMinZ != 0. || bevelMaxZ != 0.) && vertices.size() % 2 == 0) {
0173     auto halfWay = vertices.end() - vertices.size() / 2;
0174     double mult{1};
0175     auto invCtrans = ctrans.inverse();
0176     auto func = [&mult, &ctrans, &invCtrans](Vector3& v) {
0177       v = invCtrans * v;
0178       v(2) += v(1) * mult;
0179       v = ctrans * v;
0180     };
0181     if (bevelMinZ != 0.) {
0182       mult = std::tan(-bevelMinZ);
0183       std::for_each(vertices.begin(), halfWay, func);
0184     }
0185     if (bevelMaxZ != 0.) {
0186       mult = std::tan(bevelMaxZ);
0187       std::for_each(halfWay, vertices.end(), func);
0188     }
0189   }
0190   return vertices;
0191 }