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.get());
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()}, this};
0164   }
0165
0166   double u = (mab.dot(ea) - mab.dot(eb) * eaTeb) / denom;
0167   // Check if we are on the surface already
0168   IntersectionStatus status = std::abs(u) > std::abs(tolerance)
0169                                   ? IntersectionStatus::reachable
0170                                   : IntersectionStatus::onSurface;
0171   Vector3 result = ma + u * ea;
0172   // Evaluate the boundary check if requested
0173   // m_bounds == nullptr prevents unnecessary calculations for PerigeeSurface
0174   if (m_bounds && !boundaryTolerance.isInfinite()) {
0175     Vector3 vecLocal = result - mb;
0176     double cZ = vecLocal.dot(eb);
0177     double cR = (vecLocal - cZ * eb).norm();
0178     if (!m_bounds->inside({cR, cZ}, boundaryTolerance)) {
0179       status = IntersectionStatus::unreachable;
0180     }
0181   }
0182
0183   return {{Intersection3D(result, u, status), Intersection3D::invalid()}, this};
0184 }
0185
0186 BoundToFreeMatrix LineSurface::boundToFreeJacobian(
0187     const GeometryContext& gctx, const Vector3& position,
0188     const Vector3& direction) const {
0189   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0190
0191   // retrieve the reference frame
0192   auto rframe = referenceFrame(gctx, position, direction);
0193
0194   Vector2 local = *globalToLocal(gctx, position, direction,
0195                                  std::numeric_limits<double>::max());
0196
0197   // For the derivative of global position with bound angles, refer the
0198   // following white paper:
0199   // https://acts.readthedocs.io/en/latest/white_papers/line-surface-jacobian.html
0200
0201   BoundToFreeMatrix jacToGlobal =
0202       Surface::boundToFreeJacobian(gctx, position, direction);
0203
0204   // the projection of direction onto ref frame normal
0205   double ipdn = 1. / direction.dot(rframe.col(2));
0206   // build the cross product of d(D)/d(eBoundPhi) components with y axis
0207   Vector3 dDPhiY = rframe.block<3, 1>(0, 1).cross(
0208       jacToGlobal.block<3, 1>(eFreeDir0, eBoundPhi));
0209   // and the same for the d(D)/d(eTheta) components
0210   Vector3 dDThetaY = rframe.block<3, 1>(0, 1).cross(
0211       jacToGlobal.block<3, 1>(eFreeDir0, eBoundTheta));
0212   // and correct for the x axis components
0213   dDPhiY -= rframe.block<3, 1>(0, 0) * (rframe.block<3, 1>(0, 0).dot(dDPhiY));
0214   dDThetaY -=
0215       rframe.block<3, 1>(0, 0) * (rframe.block<3, 1>(0, 0).dot(dDThetaY));
0216   // set the jacobian components for global d/ phi/Theta
0217   jacToGlobal.block<3, 1>(eFreePos0, eBoundPhi) = dDPhiY * local.x() * ipdn;
0218   jacToGlobal.block<3, 1>(eFreePos0, eBoundTheta) = dDThetaY * local.x() * ipdn;
0219
0220   return jacToGlobal;
0221 }
0222
0223 FreeToPathMatrix LineSurface::freeToPathDerivative(
0224     const GeometryContext& gctx, const Vector3& position,
0225     const Vector3& direction) const {
0226   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0227
0228   // The vector between position and center
0229   Vector3 pcRowVec = position - center(gctx);
0230   // The local frame z axis
0231   Vector3 localZAxis = lineDirection(gctx);
0232   // The local z coordinate
0233   double pz = pcRowVec.dot(localZAxis);
0234   // Cosine of angle between momentum direction and local frame z axis
0235   double dz = localZAxis.dot(direction);
0236   double norm = 1 / (1 - dz * dz);
0237
0238   // Initialize the derivative of propagation path w.r.t. free parameter
0239   FreeToPathMatrix freeToPath = FreeToPathMatrix::Zero();
0240
0241   // The derivative of path w.r.t. position
0242   freeToPath.segment<3>(eFreePos0) =
0243       norm * (dz * localZAxis.transpose() - direction.transpose());
0244
0245   // The derivative of path w.r.t. direction
0246   freeToPath.segment<3>(eFreeDir0) =
0247       norm * (pz * localZAxis.transpose() - pcRowVec.transpose());
0248
0249   return freeToPath;
0250 }
0251
0252 AlignmentToPathMatrix LineSurface::alignmentToPathDerivative(
0253     const GeometryContext& gctx, const Vector3& position,
0254     const Vector3& direction) const {
0255   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0256
0257   // The vector between position and center
0258   Vector3 pcRowVec = position - center(gctx);
0259   // The local frame z axis
0260   Vector3 localZAxis = lineDirection(gctx);
0261   // The local z coordinate
0262   double pz = pcRowVec.dot(localZAxis);
0263   // Cosine of angle between momentum direction and local frame z axis
0264   double dz = localZAxis.dot(direction);
0265   double norm = 1 / (1 - dz * dz);
0266   // Calculate the derivative of local frame axes w.r.t its rotation
0267   auto [rotToLocalXAxis, rotToLocalYAxis, rotToLocalZAxis] =
0268       detail::rotationToLocalAxesDerivative(transform(gctx).rotation());
0269
0270   // Initialize the derivative of propagation path w.r.t. local frame
0271   // translation (origin) and rotation
0272   AlignmentToPathMatrix alignToPath = AlignmentToPathMatrix::Zero();
0273   alignToPath.segment<3>(eAlignmentCenter0) =
0274       norm * (direction.transpose() - dz * localZAxis.transpose());
0275   alignToPath.segment<3>(eAlignmentRotation0) =
0276       norm * (dz * pcRowVec.transpose() + pz * direction.transpose()) *
0277       rotToLocalZAxis;
0278
0279   return alignToPath;
0280 }
0281
0282 ActsMatrix<2, 3> LineSurface::localCartesianToBoundLocalDerivative(
0283     const GeometryContext& gctx, const Vector3& position) const {
0284   // calculate the transformation to local coordinates
0285   Vector3 localPosition = transform(gctx).inverse() * position;
0286   double localPhi = VectorHelpers::phi(localPosition);
0287
0288   ActsMatrix<2, 3> loc3DToLocBound = ActsMatrix<2, 3>::Zero();
0289   loc3DToLocBound << std::cos(localPhi), std::sin(localPhi), 0, 0, 0, 1;
0290
0291   return loc3DToLocBound;
0292 }
0293
0294 Vector3 LineSurface::lineDirection(const GeometryContext& gctx) const {
0295   return transform(gctx).linear().col(2);
0296 }
0297
0298 }  // namespace Acts