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