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/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/GeometryObject.hpp"
0013 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0014 #include "Acts/Surfaces/CurvilinearSurface.hpp"
0015 #include "Acts/Surfaces/EllipseBounds.hpp"
0016 #include "Acts/Surfaces/InfiniteBounds.hpp"
0017 #include "Acts/Surfaces/PlanarBounds.hpp"
0018 #include "Acts/Surfaces/RectangleBounds.hpp"
0019 #include "Acts/Surfaces/SurfaceBounds.hpp"
0020 #include "Acts/Surfaces/SurfaceError.hpp"
0021 #include "Acts/Surfaces/SurfaceMergingException.hpp"
0022 #include "Acts/Surfaces/detail/FacesHelper.hpp"
0023 #include "Acts/Surfaces/detail/PlanarHelper.hpp"
0024 #include "Acts/Utilities/AlgebraHelpers.hpp"
0025 #include "Acts/Utilities/Intersection.hpp"
0026 #include "Acts/Utilities/ThrowAssert.hpp"
0027
0028 #include <cmath>
0029 #include <numbers>
0030 #include <numeric>
0031 #include <stdexcept>
0032 #include <utility>
0033 #include <vector>
0034
0035 namespace Acts {
0036
0037 PlaneSurface::PlaneSurface(const PlaneSurface& other)
0038     : GeometryObject{}, RegularSurface(other), m_bounds(other.m_bounds) {}
0039
0040 PlaneSurface::PlaneSurface(const GeometryContext& gctx,
0041                            const PlaneSurface& other,
0042                            const Transform3& transform)
0043     : RegularSurface(gctx, other, transform), m_bounds(other.m_bounds) {}
0044
0045 PlaneSurface::PlaneSurface(std::shared_ptr<const PlanarBounds> pbounds,
0046                            const SurfacePlacementBase& placement)
0047     : RegularSurface{placement}, m_bounds(std::move(pbounds)) {
0048   // surfaces representing a detector element must have bounds
0049   throw_assert(m_bounds, "PlaneBounds must not be nullptr");
0050 }
0051
0052 PlaneSurface::PlaneSurface(const Transform3& transform,
0053                            std::shared_ptr<const PlanarBounds> pbounds)
0054     : RegularSurface(transform), m_bounds(std::move(pbounds)) {}
0055
0056 PlaneSurface& PlaneSurface::operator=(const PlaneSurface& other) {
0057   if (this != &other) {
0058     Surface::operator=(other);
0059     m_bounds = other.m_bounds;
0060   }
0061   return *this;
0062 }
0063
0064 RotationMatrix3 PlaneSurface::referenceFrame(
0065     const GeometryContext& gctx) const {
0066   return localToGlobalTransform(gctx).matrix().block<3, 3>(0, 0);
0067 }
0068
0069 Surface::SurfaceType PlaneSurface::type() const {
0070   return Surface::Plane;
0071 }
0072
0073 Vector3 PlaneSurface::localToGlobal(const GeometryContext& gctx,
0074                                     const Vector2& lposition) const {
0075   return localToGlobalTransform(gctx) * Vector3(lposition[0], lposition[1], 0.);
0076 }
0077
0078 Result<Vector2> PlaneSurface::globalToLocal(const GeometryContext& gctx,
0079                                             const Vector3& position,
0080                                             double tolerance) const {
0081   Vector3 loc3Dframe = localToGlobalTransform(gctx).inverse() * position;
0082   if (std::abs(loc3Dframe.z()) > std::abs(tolerance)) {
0083     return Result<Vector2>::failure(SurfaceError::GlobalPositionNotOnSurface);
0084   }
0085   return Result<Vector2>::success({loc3Dframe.x(), loc3Dframe.y()});
0086 }
0087
0088 std::string PlaneSurface::name() const {
0089   return "Acts::PlaneSurface";
0090 }
0091
0092 const SurfaceBounds& PlaneSurface::bounds() const {
0093   if (m_bounds) {
0094     return *m_bounds;
0095   }
0096   return s_noBounds;
0097 }
0098
0099 Polyhedron PlaneSurface::polyhedronRepresentation(
0100     const GeometryContext& gctx, unsigned int quarterSegments) const {
0101   // Prepare vertices and faces
0102   std::vector<Vector3> vertices;
0103   bool exactPolyhedron = true;
0104
0105   // If you have bounds you can create a polyhedron representation
0106   if (m_bounds) {
0107     auto vertices2D = m_bounds->vertices(quarterSegments);
0108     vertices.reserve(vertices2D.size() + 1);
0109     for (const auto& v2D : vertices2D) {
0110       vertices.push_back(localToGlobalTransform(gctx) *
0111                          Vector3(v2D.x(), v2D.y(), 0.));
0112     }
0113     bool isEllipse = bounds().type() == SurfaceBounds::eEllipse;
0114     bool innerExists = false, coversFull = false;
0115     if (isEllipse) {
0116       exactPolyhedron = false;
0117       auto vStore = bounds().values();
0118       innerExists = vStore[EllipseBounds::eInnerRx] > s_epsilon &&
0119                     vStore[EllipseBounds::eInnerRy] > s_epsilon;
0120       coversFull = std::abs(vStore[EllipseBounds::eHalfPhiSector] -
0121                             std::numbers::pi) < s_epsilon;
0122     }
0123     // All of those can be described as convex
0124     // @todo same as for Discs: coversFull is not the right criterium
0125     // for triangulation
0126     if (!isEllipse || !innerExists || !coversFull) {
0127       auto [faces, triangularMesh] =
0128           detail::FacesHelper::convexFaceMesh(vertices);
0129       return Polyhedron(vertices, faces, triangularMesh, exactPolyhedron);
0130     } else {
0131       // Two concentric rings, we use the pure concentric method momentarily,
0132       // but that creates too  many unneccesarry faces, when only two
0133       // are needed to describe the mesh, @todo investigate merging flag
0134       auto [faces, triangularMesh] =
0135           detail::FacesHelper::cylindricalFaceMesh(vertices);
0136       return Polyhedron(vertices, faces, triangularMesh, exactPolyhedron);
0137     }
0138   }
0139   throw std::domain_error(
0140       "Polyhedron representation of boundless surface not possible.");
0141 }
0142
0143 Vector3 PlaneSurface::normal(const GeometryContext& gctx,
0144                              const Vector2& /*lpos*/) const {
0145   return normal(gctx);
0146 }
0147
0148 Vector3 PlaneSurface::normal(const GeometryContext& gctx,
0149                              const Vector3& /*pos*/) const {
0150   return normal(gctx);
0151 }
0152
0153 Vector3 PlaneSurface::normal(const GeometryContext& gctx) const {
0154   return localToGlobalTransform(gctx).linear().col(2);
0155 }
0156
0157 Vector3 PlaneSurface::referencePosition(const GeometryContext& gctx,
0158                                         AxisDirection /*aDir*/) const {
0159   return center(gctx);
0160 }
0161
0162 double PlaneSurface::pathCorrection(const GeometryContext& gctx,
0163                                     const Vector3& /*position*/,
0164                                     const Vector3& direction) const {
0165   // We can ignore the global position here
0166   return 1. / std::abs(normal(gctx).dot(direction));
0167 }
0168
0169 MultiIntersection3D PlaneSurface::intersect(
0170     const GeometryContext& gctx, const Vector3& position,
0171     const Vector3& direction, const BoundaryTolerance& boundaryTolerance,
0172     double tolerance) const {
0173   // Get the contextual transform
0174   const auto& gctxTransform = localToGlobalTransform(gctx);
0175   // Use the intersection helper for planar surfaces
0176   auto intersection =
0177       PlanarHelper::intersect(gctxTransform, position, direction, tolerance);
0178   auto status = intersection.status();
0179   // Evaluate boundary check if requested (and reachable)
0180   if (intersection.status() != IntersectionStatus::unreachable) {
0181     // Built-in local to global for speed reasons
0182     const auto& tMatrix = gctxTransform.matrix();
0183     // Create the reference vector in local
0184     const Vector3 vecLocal(intersection.position() - tMatrix.block<3, 1>(0, 3));
0185     if (!insideBounds(tMatrix.block<3, 2>(0, 0).transpose() * vecLocal,
0186                       boundaryTolerance)) {
0187       status = IntersectionStatus::unreachable;
0188     }
0189   }
0190   return MultiIntersection3D(Intersection3D(intersection.position(),
0191                                             intersection.pathLength(), status));
0192 }
0193
0194 Matrix<2, 3> PlaneSurface::localCartesianToBoundLocalDerivative(
0195     const GeometryContext& /*gctx*/, const Vector3& /*position*/) const {
0196   const Matrix<2, 3> loc3DToLocBound = Matrix<2, 3>::Identity();
0197   return loc3DToLocBound;
0198 }
0199
0200 std::pair<std::shared_ptr<PlaneSurface>, bool> PlaneSurface::mergedWith(
0201     const PlaneSurface& other, AxisDirection direction,
0202     const Logger& logger) const {
0203   ACTS_VERBOSE("Merging plane surfaces in " << axisDirectionName(direction)
0204                                             << " direction");
0205
0206   if (isAlignable() || other.isAlignable()) {
0207     throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0208                                   "PlaneSurface::merge: surfaces are "
0209                                   "associated with a detector element");
0210   }
0211
0212   assert(m_transform != nullptr && other.m_transform != nullptr);
0213
0214   Transform3 otherLocal = m_transform->inverse() * *other.m_transform;
0215
0216   // TODO: Is it a good tolerance?
0217   constexpr auto tolerance = s_onSurfaceTolerance;
0218
0219   // Surface cannot have any relative rotation
0220   if ((otherLocal.rotation().matrix() - RotationMatrix3::Identity()).norm() >
0221       tolerance) {
0222     ACTS_ERROR("PlaneSurface::merge: surfaces have relative rotation");
0223     throw SurfaceMergingException(
0224         getSharedPtr(), other.getSharedPtr(),
0225         "PlaneSurface::merge: surfaces have relative rotation");
0226   }
0227
0228   const auto* thisBounds = dynamic_cast<const RectangleBounds*>(&bounds());
0229   const auto* otherBounds =
0230       dynamic_cast<const RectangleBounds*>(&other.bounds());
0231
0232   if (thisBounds == nullptr || otherBounds == nullptr) {
0233     throw SurfaceMergingException(
0234         getSharedPtr(), other.getSharedPtr(),
0235         "PlaneSurface::merge: only Rectangle Bounds are supported");
0236   }
0237
0238   if (direction != AxisDirection::AxisX && direction != AxisDirection::AxisY) {
0239     throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0240                                   "PlaneSurface::merge: invalid direction " +
0241                                       axisDirectionName(direction));
0242   }
0243
0244   bool mergeX = direction == AxisDirection::AxisX;
0245
0246   double thisHalfMerge =
0247       mergeX ? thisBounds->halfLengthX() : thisBounds->halfLengthY();
0248   double otherHalfMerge =
0249       mergeX ? otherBounds->halfLengthX() : otherBounds->halfLengthY();
0250
0251   double thisHalfNonMerge =
0252       mergeX ? thisBounds->halfLengthY() : thisBounds->halfLengthX();
0253   double otherHalfNonMerge =
0254       mergeX ? otherBounds->halfLengthY() : otherBounds->halfLengthX();
0255
0256   if (std::abs(thisHalfNonMerge - otherHalfNonMerge) > tolerance) {
0257     ACTS_ERROR(
0258         "PlaneSurface::merge: surfaces have different non-merging lengths");
0259     throw SurfaceMergingException(
0260         getSharedPtr(), other.getSharedPtr(),
0261         "PlaneSurface::merge: surfaces have different non-merging lengths");
0262   }
0263   Vector3 otherTranslation = otherLocal.translation();
0264
0265   // No translation in non-merging direction/z is allowed
0266   double nonMergeShift = mergeX ? otherTranslation.y() : otherTranslation.x();
0267
0268   if (std::abs(nonMergeShift) > tolerance ||
0269       std::abs(otherTranslation.z()) > tolerance) {
0270     ACTS_ERROR(
0271         "PlaneSurface::merge: surfaces have relative translation in y/z");
0272     throw SurfaceMergingException(
0273         getSharedPtr(), other.getSharedPtr(),
0274         "PlaneSurface::merge: surfaces have relative translation in y/z");
0275   }
0276
0277   double mergeShift = mergeX ? otherTranslation.x() : otherTranslation.y();
0278
0279   double thisMinMerge = -thisHalfMerge;
0280   double thisMaxMerge = thisHalfMerge;
0281
0282   double otherMinMerge = mergeShift - otherHalfMerge;
0283   double otherMaxMerge = mergeShift + otherHalfMerge;
0284
0285   // Surfaces have to "touch" along merging direction
0286   if (std::abs(thisMaxMerge - otherMinMerge) > tolerance &&
0287       std::abs(thisMinMerge - otherMaxMerge) > tolerance) {
0288     ACTS_ERROR(
0289         "PlaneSurface::merge: surfaces have incompatible merge bound location");
0290     throw SurfaceMergingException(
0291         getSharedPtr(), other.getSharedPtr(),
0292         "PlaneSurface::merge: surfaces have incompatible merge bound location");
0293   }
0294
0295   double newMaxMerge = std::max(thisMaxMerge, otherMaxMerge);
0296   double newMinMerge = std::min(thisMinMerge, otherMinMerge);
0297
0298   double newHalfMerge = std::midpoint(newMaxMerge, -newMinMerge);
0299   double newMidMerge = std::midpoint(newMaxMerge, newMinMerge);
0300
0301   auto newBounds =
0302       mergeX
0303           ? std::make_shared<RectangleBounds>(newHalfMerge, thisHalfNonMerge)
0304           : std::make_shared<RectangleBounds>(thisHalfNonMerge, newHalfMerge);
0305
0306   Vector3 unitDir = mergeX ? Vector3::UnitX() : Vector3::UnitY();
0307   Transform3 newTransform = *m_transform * Translation3{unitDir * newMidMerge};
0308   return {Surface::makeShared<PlaneSurface>(newTransform, newBounds),
0309           mergeShift < 0};
0310 }
0311 const std::shared_ptr<const PlanarBounds>& PlaneSurface::boundsPtr() const {
0312   return m_bounds;
0313 }
0314
0315 void PlaneSurface::assignSurfaceBounds(
0316     std::shared_ptr<const PlanarBounds> newBounds) {
0317   m_bounds = std::move(newBounds);
0318 }
0319
0320 }  // namespace Acts