src/Surfaces/PlaneSurface.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/PlaneSurface.hpp"
0010
0011 #include "Acts/Geometry/GeometryObject.hpp"
0012 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0013 #include "Acts/Surfaces/CurvilinearSurface.hpp"
0014 #include "Acts/Surfaces/EllipseBounds.hpp"
0015 #include "Acts/Surfaces/InfiniteBounds.hpp"
0016 #include "Acts/Surfaces/PlanarBounds.hpp"
0017 #include "Acts/Surfaces/SurfaceBounds.hpp"
0018 #include "Acts/Surfaces/SurfaceError.hpp"
0019 #include "Acts/Surfaces/detail/FacesHelper.hpp"
0020 #include "Acts/Surfaces/detail/PlanarHelper.hpp"
0021 #include "Acts/Utilities/Intersection.hpp"
0022 #include "Acts/Utilities/ThrowAssert.hpp"
0023
0024 #include <cmath>
0025 #include <numbers>
0026 #include <stdexcept>
0027 #include <utility>
0028 #include <vector>
0029
0030 namespace Acts {
0031
0032 PlaneSurface::PlaneSurface(const PlaneSurface& other)
0033     : GeometryObject(), RegularSurface(other), m_bounds(other.m_bounds) {}
0034
0035 PlaneSurface::PlaneSurface(const GeometryContext& gctx,
0036                            const PlaneSurface& other,
0037                            const Transform3& transform)
0038     : GeometryObject(),
0039       RegularSurface(gctx, other, transform),
0040       m_bounds(other.m_bounds) {}
0041
0042 PlaneSurface::PlaneSurface(std::shared_ptr<const PlanarBounds> pbounds,
0043                            const DetectorElementBase& detelement)
0044     : RegularSurface(detelement), m_bounds(std::move(pbounds)) {
0045   // surfaces representing a detector element must have bounds
0046   throw_assert(m_bounds, "PlaneBounds must not be nullptr");
0047 }
0048
0049 PlaneSurface::PlaneSurface(const Transform3& transform,
0050                            std::shared_ptr<const PlanarBounds> pbounds)
0051     : RegularSurface(transform), m_bounds(std::move(pbounds)) {}
0052
0053 PlaneSurface& PlaneSurface::operator=(const PlaneSurface& other) {
0054   if (this != &other) {
0055     Surface::operator=(other);
0056     m_bounds = other.m_bounds;
0057   }
0058   return *this;
0059 }
0060
0061 Surface::SurfaceType PlaneSurface::type() const {
0062   return Surface::Plane;
0063 }
0064
0065 Vector3 PlaneSurface::localToGlobal(const GeometryContext& gctx,
0066                                     const Vector2& lposition) const {
0067   return transform(gctx) * Vector3(lposition[0], lposition[1], 0.);
0068 }
0069
0070 Result<Vector2> PlaneSurface::globalToLocal(const GeometryContext& gctx,
0071                                             const Vector3& position,
0072                                             double tolerance) const {
0073   Vector3 loc3Dframe = transform(gctx).inverse() * position;
0074   if (std::abs(loc3Dframe.z()) > std::abs(tolerance)) {
0075     return Result<Vector2>::failure(SurfaceError::GlobalPositionNotOnSurface);
0076   }
0077   return Result<Vector2>::success({loc3Dframe.x(), loc3Dframe.y()});
0078 }
0079
0080 std::string PlaneSurface::name() const {
0081   return "Acts::PlaneSurface";
0082 }
0083
0084 const SurfaceBounds& PlaneSurface::bounds() const {
0085   if (m_bounds) {
0086     return (*m_bounds.get());
0087   }
0088   return s_noBounds;
0089 }
0090
0091 Polyhedron PlaneSurface::polyhedronRepresentation(
0092     const GeometryContext& gctx, unsigned int quarterSegments) const {
0093   // Prepare vertices and faces
0094   std::vector<Vector3> vertices;
0095   bool exactPolyhedron = true;
0096
0097   // If you have bounds you can create a polyhedron representation
0098   if (m_bounds) {
0099     auto vertices2D = m_bounds->vertices(quarterSegments);
0100     vertices.reserve(vertices2D.size() + 1);
0101     for (const auto& v2D : vertices2D) {
0102       vertices.push_back(transform(gctx) * Vector3(v2D.x(), v2D.y(), 0.));
0103     }
0104     bool isEllipse = bounds().type() == SurfaceBounds::eEllipse;
0105     bool innerExists = false, coversFull = false;
0106     if (isEllipse) {
0107       exactPolyhedron = false;
0108       auto vStore = bounds().values();
0109       innerExists = vStore[EllipseBounds::eInnerRx] > s_epsilon &&
0110                     vStore[EllipseBounds::eInnerRy] > s_epsilon;
0111       coversFull = std::abs(vStore[EllipseBounds::eHalfPhiSector] -
0112                             std::numbers::pi) < s_epsilon;
0113     }
0114     // All of those can be described as convex
0115     // @todo same as for Discs: coversFull is not the right criterium
0116     // for triangulation
0117     if (!isEllipse || !innerExists || !coversFull) {
0118       auto [faces, triangularMesh] =
0119           detail::FacesHelper::convexFaceMesh(vertices);
0120       return Polyhedron(vertices, faces, triangularMesh, exactPolyhedron);
0121     } else {
0122       // Two concentric rings, we use the pure concentric method momentarily,
0123       // but that creates too  many unneccesarry faces, when only two
0124       // are needed to describe the mesh, @todo investigate merging flag
0125       auto [faces, triangularMesh] =
0126           detail::FacesHelper::cylindricalFaceMesh(vertices);
0127       return Polyhedron(vertices, faces, triangularMesh, exactPolyhedron);
0128     }
0129   }
0130   throw std::domain_error(
0131       "Polyhedron representation of boundless surface not possible.");
0132 }
0133
0134 Vector3 PlaneSurface::normal(const GeometryContext& gctx,
0135                              const Vector2& /*lpos*/) const {
0136   return normal(gctx);
0137 }
0138
0139 Vector3 PlaneSurface::normal(const GeometryContext& gctx,
0140                              const Vector3& /*pos*/) const {
0141   return normal(gctx);
0142 }
0143
0144 Vector3 PlaneSurface::normal(const GeometryContext& gctx) const {
0145   return transform(gctx).linear().col(2);
0146 }
0147
0148 Vector3 PlaneSurface::referencePosition(const GeometryContext& gctx,
0149                                         AxisDirection /*aDir*/) const {
0150   return center(gctx);
0151 }
0152
0153 double PlaneSurface::pathCorrection(const GeometryContext& gctx,
0154                                     const Vector3& /*position*/,
0155                                     const Vector3& direction) const {
0156   // We can ignore the global position here
0157   return 1. / std::abs(normal(gctx).dot(direction));
0158 }
0159
0160 SurfaceMultiIntersection PlaneSurface::intersect(
0161     const GeometryContext& gctx, const Vector3& position,
0162     const Vector3& direction, const BoundaryTolerance& boundaryTolerance,
0163     double tolerance) const {
0164   // Get the contextual transform
0165   const auto& gctxTransform = transform(gctx);
0166   // Use the intersection helper for planar surfaces
0167   auto intersection =
0168       PlanarHelper::intersect(gctxTransform, position, direction, tolerance);
0169   auto status = intersection.status();
0170   // Evaluate boundary check if requested (and reachable)
0171   if (intersection.status() != IntersectionStatus::unreachable) {
0172     // Built-in local to global for speed reasons
0173     const auto& tMatrix = gctxTransform.matrix();
0174     // Create the reference vector in local
0175     const Vector3 vecLocal(intersection.position() - tMatrix.block<3, 1>(0, 3));
0176     if (!insideBounds(tMatrix.block<3, 2>(0, 0).transpose() * vecLocal,
0177                       boundaryTolerance)) {
0178       status = IntersectionStatus::unreachable;
0179     }
0180   }
0181   return {{Intersection3D(intersection.position(), intersection.pathLength(),
0182                           status),
0183            Intersection3D::invalid()},
0184           this};
0185 }
0186
0187 ActsMatrix<2, 3> PlaneSurface::localCartesianToBoundLocalDerivative(
0188     const GeometryContext& /*gctx*/, const Vector3& /*position*/) const {
0189   const ActsMatrix<2, 3> loc3DToLocBound = ActsMatrix<2, 3>::Identity();
0190   return loc3DToLocBound;
0191 }
0192
0193 }  // namespace Acts