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 }
0190
0191 ActsMatrix<2, 3> PlaneSurface::localCartesianToBoundLocalDerivative(
0192     const GeometryContext& /*gctx*/, const Vector3& /*position*/) const {
0193   const ActsMatrix<2, 3> loc3DToLocBound = ActsMatrix<2, 3>::Identity();
0194   return loc3DToLocBound;
0195 }
0196
0197 std::pair<std::shared_ptr<PlaneSurface>, bool> PlaneSurface::mergedWith(
0198     const PlaneSurface& other, AxisDirection direction,
0199     const Logger& logger) const {
0200   ACTS_VERBOSE("Merging plane surfaces in " << axisDirectionName(direction)
0201                                             << " direction");
0202
0203   if (m_associatedDetElement != nullptr ||
0204       other.m_associatedDetElement != nullptr) {
0205     throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0206                                   "PlaneSurface::merge: surfaces are "
0207                                   "associated with a detector element");
0208   }
0209
0210   assert(m_transform != nullptr && other.m_transform != nullptr);
0211
0212   Transform3 otherLocal = m_transform->inverse() * *other.m_transform;
0213
0214   // TODO: Is it a good tolerance?
0215   constexpr auto tolerance = s_onSurfaceTolerance;
0216
0217   // Surface cannot have any relative rotation
0218   if ((otherLocal.rotation().matrix() - RotationMatrix3::Identity()).norm() >
0219       tolerance) {
0220     ACTS_ERROR("PlaneSurface::merge: surfaces have relative rotation");
0221     throw SurfaceMergingException(
0222         getSharedPtr(), other.getSharedPtr(),
0223         "PlaneSurface::merge: surfaces have relative rotation");
0224   }
0225
0226   const auto* thisBounds = dynamic_cast<const RectangleBounds*>(&bounds());
0227   const auto* otherBounds =
0228       dynamic_cast<const RectangleBounds*>(&other.bounds());
0229
0230   if (thisBounds == nullptr || otherBounds == nullptr) {
0231     throw SurfaceMergingException(
0232         getSharedPtr(), other.getSharedPtr(),
0233         "PlaneSurface::merge: only Rectangle Bounds are supported");
0234   }
0235
0236   if (direction != AxisDirection::AxisX && direction != AxisDirection::AxisY) {
0237     throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0238                                   "PlaneSurface::merge: invalid direction " +
0239                                       axisDirectionName(direction));
0240   }
0241
0242   bool mergeX = direction == AxisDirection::AxisX;
0243
0244   double thisHalfMerge =
0245       mergeX ? thisBounds->halfLengthX() : thisBounds->halfLengthY();
0246   double otherHalfMerge =
0247       mergeX ? otherBounds->halfLengthX() : otherBounds->halfLengthY();
0248
0249   double thisHalfNonMerge =
0250       mergeX ? thisBounds->halfLengthY() : thisBounds->halfLengthX();
0251   double otherHalfNonMerge =
0252       mergeX ? otherBounds->halfLengthY() : otherBounds->halfLengthX();
0253
0254   if (std::abs(thisHalfNonMerge - otherHalfNonMerge) > tolerance) {
0255     ACTS_ERROR(
0256         "PlaneSurface::merge: surfaces have different non-merging lengths");
0257     throw SurfaceMergingException(
0258         getSharedPtr(), other.getSharedPtr(),
0259         "PlaneSurface::merge: surfaces have different non-merging lengths");
0260   }
0261   Vector3 otherTranslation = otherLocal.translation();
0262
0263   // No translation in non-merging direction/z is allowed
0264   double nonMergeShift = mergeX ? otherTranslation.y() : otherTranslation.x();
0265
0266   if (std::abs(nonMergeShift) > tolerance ||
0267       std::abs(otherTranslation.z()) > tolerance) {
0268     ACTS_ERROR(
0269         "PlaneSurface::merge: surfaces have relative translation in y/z");
0270     throw SurfaceMergingException(
0271         getSharedPtr(), other.getSharedPtr(),
0272         "PlaneSurface::merge: surfaces have relative translation in y/z");
0273   }
0274
0275   double mergeShift = mergeX ? otherTranslation.x() : otherTranslation.y();
0276
0277   double thisMinMerge = -thisHalfMerge;
0278   double thisMaxMerge = thisHalfMerge;
0279
0280   double otherMinMerge = mergeShift - otherHalfMerge;
0281   double otherMaxMerge = mergeShift + otherHalfMerge;
0282
0283   // Surfaces have to "touch" along merging direction
0284   if (std::abs(thisMaxMerge - otherMinMerge) > tolerance &&
0285       std::abs(thisMinMerge - otherMaxMerge) > tolerance) {
0286     ACTS_ERROR(
0287         "PlaneSurface::merge: surfaces have incompatible merge bound location");
0288     throw SurfaceMergingException(
0289         getSharedPtr(), other.getSharedPtr(),
0290         "PlaneSurface::merge: surfaces have incompatible merge bound location");
0291   }
0292
0293   double newMaxMerge = std::max(thisMaxMerge, otherMaxMerge);
0294   double newMinMerge = std::min(thisMinMerge, otherMinMerge);
0295
0296   double newHalfMerge = std::midpoint(newMaxMerge, -newMinMerge);
0297   double newMidMerge = std::midpoint(newMaxMerge, newMinMerge);
0298
0299   auto newBounds =
0300       mergeX
0301           ? std::make_shared<RectangleBounds>(newHalfMerge, thisHalfNonMerge)
0302           : std::make_shared<RectangleBounds>(thisHalfNonMerge, newHalfMerge);
0303
0304   Vector3 unitDir = mergeX ? Vector3::UnitX() : Vector3::UnitY();
0305   Transform3 newTransform = *m_transform * Translation3{unitDir * newMidMerge};
0306   return {Surface::makeShared<PlaneSurface>(newTransform, newBounds),
0307           mergeShift < 0};
0308 }
0309
0310 }  // namespace Acts