src/Surfaces/LineSurface.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/LineSurface.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/GeometryObject.hpp"
0013 #include "Acts/Surfaces/InfiniteBounds.hpp"
0014 #include "Acts/Surfaces/LineBounds.hpp"
0015 #include "Acts/Surfaces/SurfaceBounds.hpp"
0016 #include "Acts/Surfaces/SurfaceError.hpp"
0017 #include "Acts/Surfaces/detail/AlignmentHelper.hpp"
0018 #include "Acts/Utilities/AlgebraHelpers.hpp"
0019 #include "Acts/Utilities/Diagnostics.hpp"
0020 #include "Acts/Utilities/Intersection.hpp"
0021 #include "Acts/Utilities/ThrowAssert.hpp"
0022
0023 #include <cmath>
0024 #include <limits>
0025 #include <utility>
0026
0027 namespace Acts {
0028
0029 LineSurface::LineSurface(const Transform3& transform, double radius,
0030                          double halez)
0031     : Surface(transform),
0032       m_bounds(std::make_shared<const LineBounds>(radius, halez)) {}
0033
0034 LineSurface::LineSurface(const Transform3& transform,
0035                          std::shared_ptr<const LineBounds> lbounds)
0036     : Surface(transform), m_bounds(std::move(lbounds)) {}
0037
0038 LineSurface::LineSurface(std::shared_ptr<const LineBounds> lbounds,
0039                          const SurfacePlacementBase& placement)
0040     : Surface{placement}, m_bounds(std::move(lbounds)) {
0041   throw_assert(m_bounds, "LineBounds must not be nullptr");
0042 }
0043
0044 LineSurface::LineSurface(const LineSurface& other)
0045     : GeometryObject{}, Surface(other), m_bounds(other.m_bounds) {}
0046
0047 LineSurface::LineSurface(const GeometryContext& gctx, const LineSurface& other,
0048                          const Transform3& shift)
0049     : Surface(gctx, other, shift), m_bounds(other.m_bounds) {}
0050
0051 LineSurface& LineSurface::operator=(const LineSurface& other) {
0052   if (this != &other) {
0053     ACTS_PUSH_IGNORE_DEPRECATED()
0054     Surface::operator=(other);
0055     ACTS_POP_IGNORE_DEPRECATED()
0056     m_bounds = other.m_bounds;
0057   }
0058   return *this;
0059 }
0060
0061 Vector3 LineSurface::localToGlobal(const GeometryContext& gctx,
0062                                    const Vector2& lposition,
0063                                    const Vector3& direction) const {
0064   Vector3 unitZ0 = lineDirection(gctx);
0065
0066   // get the vector perpendicular to the momentum direction and the straw axis
0067   Vector3 radiusAxisGlobal = unitZ0.cross(direction);
0068   Vector3 locZinGlobal =
0069       localToGlobalTransform(gctx) * Vector3(0., 0., lposition[1]);
0070   // add loc0 * radiusAxis
0071   return Vector3(locZinGlobal + lposition[0] * radiusAxisGlobal.normalized());
0072 }
0073
0074 Result<Vector2> LineSurface::globalToLocal(const GeometryContext& gctx,
0075                                            const Vector3& position,
0076                                            const Vector3& direction,
0077                                            double tolerance) const {
0078   using VectorHelpers::perp;
0079
0080   // Bring the global position into the local frame. First remove the
0081   // translation then the rotation.
0082   Vector3 localPosition =
0083       referenceFrame(gctx, position, direction).inverse() *
0084       (position - localToGlobalTransform(gctx).translation());
0085
0086   // `localPosition.z()` is not the distance to the PCA but the smallest
0087   // distance between `position` and the imaginary plane surface defined by the
0088   // local x,y axes in the global frame and the position of the line surface.
0089   //
0090   // This check is also done for the `PlaneSurface` so I aligned the
0091   // `LineSurface` to do the same thing.
0092   if (std::abs(localPosition.z()) > std::abs(tolerance)) {
0093     return Result<Vector2>::failure(SurfaceError::GlobalPositionNotOnSurface);
0094   }
0095
0096   // Construct result from local x,y
0097   Vector2 localXY = localPosition.head<2>();
0098
0099   return Result<Vector2>::success(localXY);
0100 }
0101
0102 std::string LineSurface::name() const {
0103   return "Acts::LineSurface";
0104 }
0105
0106 RotationMatrix3 LineSurface::referenceFrame(const GeometryContext& gctx,
0107                                             const Vector3& /*position*/,
0108                                             const Vector3& direction) const {
0109   Vector3 unitZ0 = lineDirection(gctx);
0110   Vector3 unitD0 = unitZ0.cross(direction).normalized();
0111   Vector3 unitDistance = unitD0.cross(unitZ0);
0112
0113   RotationMatrix3 mFrame;
0114   mFrame.col(0) = unitD0;
0115   mFrame.col(1) = unitZ0;
0116   mFrame.col(2) = unitDistance;
0117
0118   return mFrame;
0119 }
0120
0121 double LineSurface::pathCorrection(const GeometryContext& /*gctx*/,
0122                                    const Vector3& /*pos*/,
0123                                    const Vector3& /*mom*/) const {
0124   return 1.;
0125 }
0126
0127 Vector3 LineSurface::referencePosition(const GeometryContext& gctx,
0128                                        AxisDirection /*aDir*/) const {
0129   return center(gctx);
0130 }
0131
0132 Vector3 LineSurface::normal(const GeometryContext& gctx, const Vector3& pos,
0133                             const Vector3& direction) const {
0134   auto ref = referenceFrame(gctx, pos, direction);
0135   return ref.col(2);
0136 }
0137
0138 const SurfaceBounds& LineSurface::bounds() const {
0139   if (m_bounds) {
0140     return *m_bounds;
0141   }
0142   return s_noBounds;
0143 }
0144
0145 MultiIntersection3D LineSurface::intersect(
0146     const GeometryContext& gctx, const Vector3& position,
0147     const Vector3& direction, const BoundaryTolerance& boundaryTolerance,
0148     double tolerance) const {
0149   // The nomenclature is following the header file and doxygen documentation
0150
0151   const Vector3& ma = position;
0152   const Vector3& ea = direction;
0153
0154   // Origin of the line surface
0155   Vector3 mb = localToGlobalTransform(gctx).translation();
0156   // Line surface axis
0157   Vector3 eb = lineDirection(gctx);
0158
0159   // Now go ahead and solve for the closest approach
0160   Vector3 mab = mb - ma;
0161   double eaTeb = ea.dot(eb);
0162   double denom = 1 - eaTeb * eaTeb;
0163
0164   // `tolerance` does not really have a meaning here it is just a sufficiently
0165   // small number so `u` does not explode
0166   if (std::abs(denom) < std::abs(tolerance)) {
0167     // return a false intersection
0168     return MultiIntersection3D(Intersection3D::Invalid());
0169   }
0170
0171   double u = (mab.dot(ea) - mab.dot(eb) * eaTeb) / denom;
0172   // Check if we are on the surface already
0173   IntersectionStatus status = std::abs(u) > std::abs(tolerance)
0174                                   ? IntersectionStatus::reachable
0175                                   : IntersectionStatus::onSurface;
0176   Vector3 result = ma + u * ea;
0177   // Evaluate the boundary check if requested
0178   // m_bounds == nullptr prevents unnecessary calculations for PerigeeSurface
0179   if (m_bounds && !boundaryTolerance.isInfinite()) {
0180     Vector3 vecLocal = result - mb;
0181     double cZ = vecLocal.dot(eb);
0182     double cR = (vecLocal - cZ * eb).norm();
0183     if (!m_bounds->inside({cR, cZ}, boundaryTolerance)) {
0184       status = IntersectionStatus::unreachable;
0185     }
0186   }
0187
0188   return MultiIntersection3D(Intersection3D(result, u, status));
0189 }
0190
0191 BoundToFreeMatrix LineSurface::boundToFreeJacobian(
0192     const GeometryContext& gctx, const Vector3& position,
0193     const Vector3& direction) const {
0194   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0195
0196   // retrieve the reference frame
0197   auto rframe = referenceFrame(gctx, position, direction);
0198
0199   Vector2 local = *globalToLocal(gctx, position, direction,
0200                                  std::numeric_limits<double>::max());
0201
0202   // For the derivative of global position with bound angles, refer the
0203   // following white paper:
0204   // https://acts.readthedocs.io/en/latest/white_papers/line-surface-jacobian.html
0205
0206   BoundToFreeMatrix jacToGlobal =
0207       Surface::boundToFreeJacobian(gctx, position, direction);
0208
0209   // the projection of direction onto ref frame normal
0210   double ipdn = 1. / direction.dot(rframe.col(2));
0211   // build the cross product of d(D)/d(eBoundPhi) components with y axis
0212   Vector3 dDPhiY = rframe.block<3, 1>(0, 1).cross(
0213       jacToGlobal.block<3, 1>(eFreeDir0, eBoundPhi));
0214   // and the same for the d(D)/d(eTheta) components
0215   Vector3 dDThetaY = rframe.block<3, 1>(0, 1).cross(
0216       jacToGlobal.block<3, 1>(eFreeDir0, eBoundTheta));
0217   // and correct for the x axis components
0218   dDPhiY -= rframe.block<3, 1>(0, 0) * (rframe.block<3, 1>(0, 0).dot(dDPhiY));
0219   dDThetaY -=
0220       rframe.block<3, 1>(0, 0) * (rframe.block<3, 1>(0, 0).dot(dDThetaY));
0221   // set the jacobian components for global d/ phi/Theta
0222   jacToGlobal.block<3, 1>(eFreePos0, eBoundPhi) = dDPhiY * local.x() * ipdn;
0223   jacToGlobal.block<3, 1>(eFreePos0, eBoundTheta) = dDThetaY * local.x() * ipdn;
0224
0225   return jacToGlobal;
0226 }
0227
0228 FreeToPathMatrix LineSurface::freeToPathDerivative(
0229     const GeometryContext& gctx, const Vector3& position,
0230     const Vector3& direction) const {
0231   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0232
0233   // The vector between position and center
0234   Vector3 pcRowVec = position - center(gctx);
0235   // The local frame z axis
0236   Vector3 localZAxis = lineDirection(gctx);
0237   // The local z coordinate
0238   double pz = pcRowVec.dot(localZAxis);
0239   // Cosine of angle between momentum direction and local frame z axis
0240   double dz = localZAxis.dot(direction);
0241   double norm = 1 / (1 - dz * dz);
0242
0243   // Initialize the derivative of propagation path w.r.t. free parameter
0244   FreeToPathMatrix freeToPath = FreeToPathMatrix::Zero();
0245
0246   // The derivative of path w.r.t. position
0247   freeToPath.segment<3>(eFreePos0) =
0248       norm * (dz * localZAxis.transpose() - direction.transpose());
0249
0250   // The derivative of path w.r.t. direction
0251   freeToPath.segment<3>(eFreeDir0) =
0252       norm * (pz * localZAxis.transpose() - pcRowVec.transpose());
0253
0254   return freeToPath;
0255 }
0256
0257 AlignmentToPathMatrix LineSurface::alignmentToPathDerivative(
0258     const GeometryContext& gctx, const Vector3& position,
0259     const Vector3& direction) const {
0260   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0261
0262   // The vector between position and center
0263   Vector3 pcRowVec = position - center(gctx);
0264   // The local frame z axis
0265   Vector3 localZAxis = lineDirection(gctx);
0266   // The local z coordinate
0267   double pz = pcRowVec.dot(localZAxis);
0268   // Cosine of angle between momentum direction and local frame z axis
0269   double dz = localZAxis.dot(direction);
0270   double norm = 1 / (1 - dz * dz);
0271   // Calculate the derivative of local frame axes w.r.t its rotation
0272   auto [rotToLocalXAxis, rotToLocalYAxis, rotToLocalZAxis] =
0273       detail::rotationToLocalAxesDerivative(
0274           localToGlobalTransform(gctx).rotation());
0275
0276   // Initialize the derivative of propagation path w.r.t. local frame
0277   // translation (origin) and rotation
0278   AlignmentToPathMatrix alignToPath = AlignmentToPathMatrix::Zero();
0279   alignToPath.segment<3>(eAlignmentCenter0) =
0280       norm * (direction.transpose() - dz * localZAxis.transpose());
0281   alignToPath.segment<3>(eAlignmentRotation0) =
0282       norm * (dz * pcRowVec.transpose() + pz * direction.transpose()) *
0283       rotToLocalZAxis;
0284
0285   return alignToPath;
0286 }
0287
0288 Matrix<2, 3> LineSurface::localCartesianToBoundLocalDerivative(
0289     const GeometryContext& gctx, const Vector3& position) const {
0290   // calculate the transformation to local coordinates
0291   Vector3 localPosition = localToGlobalTransform(gctx).inverse() * position;
0292   double localPhi = VectorHelpers::phi(localPosition);
0293
0294   Matrix<2, 3> loc3DToLocBound = Matrix<2, 3>::Zero();
0295   loc3DToLocBound << std::cos(localPhi), std::sin(localPhi), 0, 0, 0, 1;
0296
0297   return loc3DToLocBound;
0298 }
0299
0300 Vector3 LineSurface::lineDirection(const GeometryContext& gctx) const {
0301   return localToGlobalTransform(gctx).linear().col(2);
0302 }
0303
0304 const std::shared_ptr<const LineBounds>& LineSurface::boundsPtr() const {
0305   return m_bounds;
0306 }
0307
0308 void LineSurface::assignSurfaceBounds(
0309     std::shared_ptr<const LineBounds> newBounds) {
0310   m_bounds = std::move(newBounds);
0311 }
0312
0313 }  // namespace Acts