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