src/Surfaces/DiscSurface.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2016-2021 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/DiscSurface.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/GeometryObject.hpp"
0013 #include "Acts/Surfaces/BoundaryCheck.hpp"
0014 #include "Acts/Surfaces/DiscBounds.hpp"
0015 #include "Acts/Surfaces/DiscTrapezoidBounds.hpp"
0016 #include "Acts/Surfaces/InfiniteBounds.hpp"
0017 #include "Acts/Surfaces/RadialBounds.hpp"
0018 #include "Acts/Surfaces/SurfaceBounds.hpp"
0019 #include "Acts/Surfaces/SurfaceError.hpp"
0020 #include "Acts/Surfaces/detail/FacesHelper.hpp"
0021 #include "Acts/Surfaces/detail/PlanarHelper.hpp"
0022 #include "Acts/Utilities/Intersection.hpp"
0023 #include "Acts/Utilities/JacobianHelpers.hpp"
0024 #include "Acts/Utilities/ThrowAssert.hpp"
0025
0026 #include <cmath>
0027 #include <stdexcept>
0028 #include <utility>
0029 #include <vector>
0030
0031 using Acts::VectorHelpers::perp;
0032 using Acts::VectorHelpers::phi;
0033
0034 Acts::DiscSurface::DiscSurface(const DiscSurface& other)
0035     : GeometryObject(), RegularSurface(other), m_bounds(other.m_bounds) {}
0036
0037 Acts::DiscSurface::DiscSurface(const GeometryContext& gctx,
0038                                const DiscSurface& other,
0039                                const Transform3& shift)
0040     : GeometryObject(),
0041       RegularSurface(gctx, other, shift),
0042       m_bounds(other.m_bounds) {}
0043
0044 Acts::DiscSurface::DiscSurface(const Transform3& transform, double rmin,
0045                                double rmax, double hphisec)
0046     : GeometryObject(),
0047       RegularSurface(transform),
0048       m_bounds(std::make_shared<const RadialBounds>(rmin, rmax, hphisec)) {}
0049
0050 Acts::DiscSurface::DiscSurface(const Transform3& transform, double minhalfx,
0051                                double maxhalfx, double minR, double maxR,
0052                                double avephi, double stereo)
0053     : GeometryObject(),
0054       RegularSurface(transform),
0055       m_bounds(std::make_shared<const DiscTrapezoidBounds>(
0056           minhalfx, maxhalfx, minR, maxR, avephi, stereo)) {}
0057
0058 Acts::DiscSurface::DiscSurface(const Transform3& transform,
0059                                std::shared_ptr<const DiscBounds> dbounds)
0060     : GeometryObject(),
0061       RegularSurface(transform),
0062       m_bounds(std::move(dbounds)) {}
0063
0064 Acts::DiscSurface::DiscSurface(std::shared_ptr<const DiscBounds> dbounds,
0065                                const DetectorElementBase& detelement)
0066     : GeometryObject(),
0067       RegularSurface(detelement),
0068       m_bounds(std::move(dbounds)) {
0069   throw_assert(m_bounds, "nullptr as DiscBounds");
0070 }
0071
0072 Acts::DiscSurface& Acts::DiscSurface::operator=(const DiscSurface& other) {
0073   if (this != &other) {
0074     Acts::Surface::operator=(other);
0075     m_bounds = other.m_bounds;
0076   }
0077   return *this;
0078 }
0079
0080 Acts::Surface::SurfaceType Acts::DiscSurface::type() const {
0081   return Surface::Disc;
0082 }
0083
0084 Acts::Vector3 Acts::DiscSurface::localToGlobal(const GeometryContext& gctx,
0085                                                const Vector2& lposition) const {
0086   // create the position in the local 3d frame
0087   Vector3 loc3Dframe(
0088       lposition[Acts::eBoundLoc0] * cos(lposition[Acts::eBoundLoc1]),
0089       lposition[Acts::eBoundLoc0] * sin(lposition[Acts::eBoundLoc1]), 0.);
0090   // transform to globalframe
0091   return transform(gctx) * loc3Dframe;
0092 }
0093
0094 Acts::Result<Acts::Vector2> Acts::DiscSurface::globalToLocal(
0095     const GeometryContext& gctx, const Vector3& position,
0096     double tolerance) const {
0097   // transport it to the globalframe
0098   Vector3 loc3Dframe = (transform(gctx).inverse()) * position;
0099   if (std::abs(loc3Dframe.z()) > std::abs(tolerance)) {
0100     return Result<Vector2>::failure(SurfaceError::GlobalPositionNotOnSurface);
0101   }
0102   return Result<Acts::Vector2>::success({perp(loc3Dframe), phi(loc3Dframe)});
0103 }
0104
0105 Acts::Vector2 Acts::DiscSurface::localPolarToLocalCartesian(
0106     const Vector2& locpol) const {
0107   const DiscTrapezoidBounds* dtbo =
0108       dynamic_cast<const Acts::DiscTrapezoidBounds*>(&(bounds()));
0109   if (dtbo != nullptr) {
0110     double rMedium = dtbo->rCenter();
0111     double phi = dtbo->get(DiscTrapezoidBounds::eAveragePhi);
0112
0113     Vector2 polarCenter(rMedium, phi);
0114     Vector2 cartCenter = localPolarToCartesian(polarCenter);
0115     Vector2 cartPos = localPolarToCartesian(locpol);
0116     Vector2 Pos = cartPos - cartCenter;
0117
0118     Acts::Vector2 locPos(
0119         Pos[Acts::eBoundLoc0] * sin(phi) - Pos[Acts::eBoundLoc1] * cos(phi),
0120         Pos[Acts::eBoundLoc1] * sin(phi) + Pos[Acts::eBoundLoc0] * cos(phi));
0121     return Vector2(locPos[Acts::eBoundLoc0], locPos[Acts::eBoundLoc1]);
0122   }
0123   return Vector2(locpol[Acts::eBoundLoc0] * cos(locpol[Acts::eBoundLoc1]),
0124                  locpol[Acts::eBoundLoc0] * sin(locpol[Acts::eBoundLoc1]));
0125 }
0126
0127 Acts::Vector3 Acts::DiscSurface::localCartesianToGlobal(
0128     const GeometryContext& gctx, const Vector2& lposition) const {
0129   Vector3 loc3Dframe(lposition[Acts::eBoundLoc0], lposition[Acts::eBoundLoc1],
0130                      0.);
0131   return Vector3(transform(gctx) * loc3Dframe);
0132 }
0133
0134 Acts::Vector2 Acts::DiscSurface::globalToLocalCartesian(
0135     const GeometryContext& gctx, const Vector3& position,
0136     double /*direction*/) const {
0137   Vector3 loc3Dframe = (transform(gctx).inverse()) * position;
0138   return Vector2(loc3Dframe.x(), loc3Dframe.y());
0139 }
0140
0141 std::string Acts::DiscSurface::name() const {
0142   return "Acts::DiscSurface";
0143 }
0144
0145 const Acts::SurfaceBounds& Acts::DiscSurface::bounds() const {
0146   if (m_bounds) {
0147     return (*(m_bounds.get()));
0148   }
0149   return s_noBounds;
0150 }
0151
0152 Acts::Polyhedron Acts::DiscSurface::polyhedronRepresentation(
0153     const GeometryContext& gctx, std::size_t lseg) const {
0154   // Prepare vertices and faces
0155   std::vector<Vector3> vertices;
0156   std::vector<Polyhedron::FaceType> faces;
0157   std::vector<Polyhedron::FaceType> triangularMesh;
0158
0159   // Understand the disc
0160   bool fullDisc = m_bounds->coversFullAzimuth();
0161   bool toCenter = m_bounds->rMin() < s_onSurfaceTolerance;
0162   // If you have bounds you can create a polyhedron representation
0163   bool exactPolyhedron = (m_bounds->type() == SurfaceBounds::eDiscTrapezoid);
0164   bool addCentreFromConvexFace = (m_bounds->type() != SurfaceBounds::eAnnulus);
0165   if (m_bounds) {
0166     auto vertices2D = m_bounds->vertices(lseg);
0167     vertices.reserve(vertices2D.size() + 1);
0168     Vector3 wCenter(0., 0., 0);
0169     for (const auto& v2D : vertices2D) {
0170       vertices.push_back(transform(gctx) * Vector3(v2D.x(), v2D.y(), 0.));
0171       wCenter += (*vertices.rbegin());
0172     }
0173     // These are convex shapes, use the helper method
0174     // For rings there's a sweet spot when this stops working
0175     if (exactPolyhedron || toCenter || !fullDisc) {
0176       // Transform them into the vertex frame
0177       wCenter *= 1. / vertices.size();
0178       if (addCentreFromConvexFace) {
0179         vertices.push_back(wCenter);
0180       }
0181       auto facesMesh = detail::FacesHelper::convexFaceMesh(vertices, true);
0182       faces = facesMesh.first;
0183       triangularMesh = facesMesh.second;
0184     } else {
0185       // Two concentric rings, we use the pure concentric method momentarily,
0186       // but that creates too  many unneccesarry faces, when only two
0187       // are needed to describe the mesh, @todo investigate merging flag
0188       auto facesMesh = detail::FacesHelper::cylindricalFaceMesh(vertices, true);
0189       faces = facesMesh.first;
0190       triangularMesh = facesMesh.second;
0191     }
0192   } else {
0193     throw std::domain_error(
0194         "Polyhedron repr of boundless surface not possible.");
0195   }
0196   return Polyhedron(vertices, faces, triangularMesh, exactPolyhedron);
0197 }
0198
0199 Acts::Vector2 Acts::DiscSurface::localPolarToCartesian(
0200     const Vector2& lpolar) const {
0201   return Vector2(lpolar[eBoundLoc0] * cos(lpolar[eBoundLoc1]),
0202                  lpolar[eBoundLoc0] * sin(lpolar[eBoundLoc1]));
0203 }
0204
0205 Acts::Vector2 Acts::DiscSurface::localCartesianToPolar(
0206     const Vector2& lcart) const {
0207   return Vector2(std::hypot(lcart[eBoundLoc0], lcart[eBoundLoc1]),
0208                  std::atan2(lcart[eBoundLoc1], lcart[eBoundLoc0]));
0209 }
0210
0211 Acts::BoundToFreeMatrix Acts::DiscSurface::boundToFreeJacobian(
0212     const GeometryContext& gctx, const Vector3& position,
0213     const Vector3& direction) const {
0214   assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
0215
0216   // The measurement frame of the surface
0217   RotationMatrix3 rframeT =
0218       referenceFrame(gctx, position, direction).transpose();
0219
0220   // calculate the transformation to local coordinates
0221   const Vector3 posLoc = transform(gctx).inverse() * position;
0222   const double lr = perp(posLoc);
0223   const double lphi = phi(posLoc);
0224   const double lcphi = std::cos(lphi);
0225   const double lsphi = std::sin(lphi);
0226   // rotate into the polar coorindates
0227   auto lx = rframeT.block<1, 3>(0, 0);
0228   auto ly = rframeT.block<1, 3>(1, 0);
0229   // Initialize the jacobian from local to global
0230   BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();
0231   // the local error components - rotated from reference frame
0232   jacToGlobal.block<3, 1>(eFreePos0, eBoundLoc0) = lcphi * lx + lsphi * ly;
0233   jacToGlobal.block<3, 1>(eFreePos0, eBoundLoc1) =
0234       lr * (lcphi * ly - lsphi * lx);
0235   // the time component
0236   jacToGlobal(eFreeTime, eBoundTime) = 1;
0237   // the momentum components
0238   jacToGlobal.block<3, 2>(eFreeDir0, eBoundPhi) =
0239       sphericalToFreeDirectionJacobian(direction);
0240   jacToGlobal(eFreeQOverP, eBoundQOverP) = 1;
0241   return jacToGlobal;
0242 }
0243
0244 Acts::FreeToBoundMatrix Acts::DiscSurface::freeToBoundJacobian(
0245     const GeometryContext& gctx, const Vector3& position,
0246     const Vector3& direction) const {
0247   using VectorHelpers::perp;
0248   using VectorHelpers::phi;
0249
0250   assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
0251
0252   // The measurement frame of the surface
0253   RotationMatrix3 rframeT =
0254       referenceFrame(gctx, position, direction).transpose();
0255
0256   // calculate the transformation to local coordinates
0257   const Vector3 posLoc = transform(gctx).inverse() * position;
0258   const double lr = perp(posLoc);
0259   const double lphi = phi(posLoc);
0260   const double lcphi = std::cos(lphi);
0261   const double lsphi = std::sin(lphi);
0262   // rotate into the polar coorindates
0263   auto lx = rframeT.block<1, 3>(0, 0);
0264   auto ly = rframeT.block<1, 3>(1, 0);
0265   // Initialize the jacobian from global to local
0266   FreeToBoundMatrix jacToLocal = FreeToBoundMatrix::Zero();
0267   // Local position component
0268   jacToLocal.block<1, 3>(eBoundLoc0, eFreePos0) = lcphi * lx + lsphi * ly;
0269   jacToLocal.block<1, 3>(eBoundLoc1, eFreePos0) =
0270       (lcphi * ly - lsphi * lx) / lr;
0271   // Time element
0272   jacToLocal(eBoundTime, eFreeTime) = 1;
0273   // Directional and momentum elements for reference frame surface
0274   jacToLocal.block<2, 3>(eBoundPhi, eFreeDir0) =
0275       freeToSphericalDirectionJacobian(direction);
0276   jacToLocal(eBoundQOverP, eFreeQOverP) = 1;
0277   return jacToLocal;
0278 }
0279
0280 Acts::SurfaceMultiIntersection Acts::DiscSurface::intersect(
0281     const GeometryContext& gctx, const Vector3& position,
0282     const Vector3& direction, const BoundaryCheck& bcheck,
0283     ActsScalar tolerance) const {
0284   // Get the contextual transform
0285   auto gctxTransform = transform(gctx);
0286   // Use the intersection helper for planar surfaces
0287   auto intersection =
0288       PlanarHelper::intersect(gctxTransform, position, direction, tolerance);
0289   auto status = intersection.status();
0290   // Evaluate boundary check if requested (and reachable)
0291   if (intersection.status() != Intersection3D::Status::unreachable &&
0292       bcheck.isEnabled() && m_bounds != nullptr) {
0293     // Built-in local to global for speed reasons
0294     const auto& tMatrix = gctxTransform.matrix();
0295     const Vector3 vecLocal(intersection.position() - tMatrix.block<3, 1>(0, 3));
0296     const Vector2 lcartesian = tMatrix.block<3, 2>(0, 0).transpose() * vecLocal;
0297     if (bcheck.type() == BoundaryCheck::Type::eAbsolute &&
0298         m_bounds->coversFullAzimuth()) {
0299       double modifiedTolerance = tolerance + bcheck.tolerance()[eBoundLoc0];
0300       if (!m_bounds->insideRadialBounds(VectorHelpers::perp(lcartesian),
0301                                         modifiedTolerance)) {
0302         status = Intersection3D::Status::missed;
0303       }
0304     } else if (!insideBounds(localCartesianToPolar(lcartesian), bcheck)) {
0305       status = Intersection3D::Status::missed;
0306     }
0307   }
0308   return {{Intersection3D(intersection.position(), intersection.pathLength(),
0309                           status),
0310            Intersection3D::invalid()},
0311           this};
0312 }
0313
0314 Acts::ActsMatrix<2, 3> Acts::DiscSurface::localCartesianToBoundLocalDerivative(
0315     const GeometryContext& gctx, const Vector3& position) const {
0316   using VectorHelpers::perp;
0317   using VectorHelpers::phi;
0318   // The local frame transform
0319   const auto& sTransform = transform(gctx);
0320   // calculate the transformation to local coordinates
0321   const Vector3 localPos = sTransform.inverse() * position;
0322   const double lr = perp(localPos);
0323   const double lphi = phi(localPos);
0324   const double lcphi = std::cos(lphi);
0325   const double lsphi = std::sin(lphi);
0326   ActsMatrix<2, 3> loc3DToLocBound = ActsMatrix<2, 3>::Zero();
0327   loc3DToLocBound << lcphi, lsphi, 0, -lsphi / lr, lcphi / lr, 0;
0328
0329   return loc3DToLocBound;
0330 }
0331
0332 Acts::Vector3 Acts::DiscSurface::normal(const GeometryContext& gctx,
0333                                         const Vector2& /*lposition*/) const {
0334   return normal(gctx);
0335 }
0336
0337 Acts::Vector3 Acts::DiscSurface::normal(const GeometryContext& gctx,
0338                                         const Vector3& /*position*/) const {
0339   return normal(gctx);
0340 }
0341
0342 Acts::Vector3 Acts::DiscSurface::normal(const GeometryContext& gctx) const {
0343   // fast access via transform matrix (and not rotation())
0344   const auto& tMatrix = transform(gctx).matrix();
0345   return Vector3(tMatrix(0, 2), tMatrix(1, 2), tMatrix(2, 2));
0346 }
0347
0348 Acts::Vector3 Acts::DiscSurface::binningPosition(const GeometryContext& gctx,
0349                                                  BinningValue bValue) const {
0350   if (bValue == binR || bValue == binPhi) {
0351     double r = m_bounds->binningValueR();
0352     double phi = m_bounds->binningValuePhi();
0353     return localToGlobal(gctx, Vector2{r, phi}, Vector3{});
0354   }
0355   return center(gctx);
0356 }
0357
0358 double Acts::DiscSurface::binningPositionValue(const GeometryContext& gctx,
0359                                                BinningValue bValue) const {
0360   if (bValue == binR) {
0361     return VectorHelpers::perp(binningPosition(gctx, bValue));
0362   }
0363   if (bValue == binPhi) {
0364     return VectorHelpers::phi(binningPosition(gctx, bValue));
0365   }
0366
0367   return GeometryObject::binningPositionValue(gctx, bValue);
0368 }
0369
0370 double Acts::DiscSurface::pathCorrection(const GeometryContext& gctx,
0371                                          const Vector3& /*position*/,
0372                                          const Vector3& direction) const {
0373   /// we can ignore the global position here
0374   return 1. / std::abs(normal(gctx).dot(direction));
0375 }