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