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