src/Surfaces/Surface.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/Surface.hpp"
0010
0011 #include "Acts/Definitions/Common.hpp"
0012 #include "Acts/Surfaces/SurfaceBounds.hpp"
0013 #include "Acts/Surfaces/detail/AlignmentHelper.hpp"
0014 #include "Acts/Utilities/JacobianHelpers.hpp"
0015 #include "Acts/Utilities/VectorHelpers.hpp"
0016
0017 #include <iomanip>
0018 #include <utility>
0019
0020 std::array<std::string, Acts::Surface::SurfaceType::Other>
0021     Acts::Surface::s_surfaceTypeNames = {
0022         "Cone", "Cylinder", "Disc", "Perigee", "Plane", "Straw", "Curvilinear"};
0023
0024 Acts::Surface::Surface(const Transform3& transform)
0025     : GeometryObject(), m_transform(transform) {}
0026
0027 Acts::Surface::Surface(const DetectorElementBase& detelement)
0028     : GeometryObject(), m_associatedDetElement(&detelement) {}
0029
0030 Acts::Surface::Surface(const Surface& other)
0031     : GeometryObject(other),
0032       std::enable_shared_from_this<Surface>(),
0033       m_transform(other.m_transform),
0034       m_surfaceMaterial(other.m_surfaceMaterial) {}
0035
0036 Acts::Surface::Surface(const GeometryContext& gctx, const Surface& other,
0037                        const Transform3& shift)
0038     : GeometryObject(),
0039       m_transform(shift * other.transform(gctx)),
0040       m_surfaceMaterial(other.m_surfaceMaterial) {}
0041
0042 Acts::Surface::~Surface() = default;
0043
0044 bool Acts::Surface::isOnSurface(const GeometryContext& gctx,
0045                                 const Vector3& position,
0046                                 const Vector3& direction,
0047                                 const BoundaryCheck& bcheck) const {
0048   // global to local transformation
0049   auto lpResult = globalToLocal(gctx, position, direction);
0050   if (lpResult.ok()) {
0051     return bcheck.isEnabled() ? bounds().inside(lpResult.value(), bcheck)
0052                               : true;
0053   }
0054   return false;
0055 }
0056
0057 Acts::AlignmentToBoundMatrix Acts::Surface::alignmentToBoundDerivative(
0058     const GeometryContext& gctx, const Vector3& position,
0059     const Vector3& direction, const FreeVector& pathDerivative) const {
0060   assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
0061
0062   // 1) Calculate the derivative of bound parameter local position w.r.t.
0063   // alignment parameters without path length correction
0064   const auto alignToBoundWithoutCorrection =
0065       alignmentToBoundDerivativeWithoutCorrection(gctx, position, direction);
0066   // 2) Calculate the derivative of path length w.r.t. alignment parameters
0067   const auto alignToPath = alignmentToPathDerivative(gctx, position, direction);
0068   // 3) Calculate the jacobian from free parameters to bound parameters
0069   FreeToBoundMatrix jacToLocal = freeToBoundJacobian(gctx, position, direction);
0070   // 4) The derivative of bound parameters w.r.t. alignment
0071   // parameters is alignToBoundWithoutCorrection +
0072   // jacToLocal*pathDerivative*alignToPath
0073   AlignmentToBoundMatrix alignToBound =
0074       alignToBoundWithoutCorrection + jacToLocal * pathDerivative * alignToPath;
0075
0076   return alignToBound;
0077 }
0078
0079 Acts::AlignmentToBoundMatrix
0080 Acts::Surface::alignmentToBoundDerivativeWithoutCorrection(
0081     const GeometryContext& gctx, const Vector3& position,
0082     const Vector3& direction) const {
0083   (void)direction;
0084   assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
0085
0086   // The vector between position and center
0087   const auto pcRowVec = (position - center(gctx)).transpose().eval();
0088   // The local frame rotation
0089   const auto& rotation = transform(gctx).rotation();
0090   // The axes of local frame
0091   const auto& localXAxis = rotation.col(0);
0092   const auto& localYAxis = rotation.col(1);
0093   const auto& localZAxis = rotation.col(2);
0094   // Calculate the derivative of local frame axes w.r.t its rotation
0095   const auto [rotToLocalXAxis, rotToLocalYAxis, rotToLocalZAxis] =
0096       detail::rotationToLocalAxesDerivative(rotation);
0097   // Calculate the derivative of local 3D Cartesian coordinates w.r.t.
0098   // alignment parameters (without path correction)
0099   AlignmentToPositionMatrix alignToLoc3D = AlignmentToPositionMatrix::Zero();
0100   alignToLoc3D.block<1, 3>(eX, eAlignmentCenter0) = -localXAxis.transpose();
0101   alignToLoc3D.block<1, 3>(eY, eAlignmentCenter0) = -localYAxis.transpose();
0102   alignToLoc3D.block<1, 3>(eZ, eAlignmentCenter0) = -localZAxis.transpose();
0103   alignToLoc3D.block<1, 3>(eX, eAlignmentRotation0) =
0104       pcRowVec * rotToLocalXAxis;
0105   alignToLoc3D.block<1, 3>(eY, eAlignmentRotation0) =
0106       pcRowVec * rotToLocalYAxis;
0107   alignToLoc3D.block<1, 3>(eZ, eAlignmentRotation0) =
0108       pcRowVec * rotToLocalZAxis;
0109   // The derivative of bound local w.r.t. local 3D Cartesian coordinates
0110   ActsMatrix<2, 3> loc3DToBoundLoc =
0111       localCartesianToBoundLocalDerivative(gctx, position);
0112   // Initialize the derivative of bound parameters w.r.t. alignment
0113   // parameters without path correction
0114   AlignmentToBoundMatrix alignToBound = AlignmentToBoundMatrix::Zero();
0115   // It's only relevant with the bound local position without path correction
0116   alignToBound.block<2, eAlignmentSize>(eBoundLoc0, eAlignmentCenter0) =
0117       loc3DToBoundLoc * alignToLoc3D;
0118   return alignToBound;
0119 }
0120
0121 Acts::AlignmentToPathMatrix Acts::Surface::alignmentToPathDerivative(
0122     const GeometryContext& gctx, const Vector3& position,
0123     const Vector3& direction) const {
0124   assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
0125
0126   // The vector between position and center
0127   const auto pcRowVec = (position - center(gctx)).transpose().eval();
0128   // The local frame rotation
0129   const auto& rotation = transform(gctx).rotation();
0130   // The local frame z axis
0131   const auto& localZAxis = rotation.col(2);
0132   // Cosine of angle between momentum direction and local frame z axis
0133   const auto dz = localZAxis.dot(direction);
0134   // Calculate the derivative of local frame axes w.r.t its rotation
0135   const auto [rotToLocalXAxis, rotToLocalYAxis, rotToLocalZAxis] =
0136       detail::rotationToLocalAxesDerivative(rotation);
0137   // Initialize the derivative of propagation path w.r.t. local frame
0138   // translation (origin) and rotation
0139   AlignmentToPathMatrix alignToPath = AlignmentToPathMatrix::Zero();
0140   alignToPath.segment<3>(eAlignmentCenter0) = localZAxis.transpose() / dz;
0141   alignToPath.segment<3>(eAlignmentRotation0) =
0142       -pcRowVec * rotToLocalZAxis / dz;
0143
0144   return alignToPath;
0145 }
0146
0147 std::shared_ptr<Acts::Surface> Acts::Surface::getSharedPtr() {
0148   return shared_from_this();
0149 }
0150
0151 std::shared_ptr<const Acts::Surface> Acts::Surface::getSharedPtr() const {
0152   return shared_from_this();
0153 }
0154
0155 Acts::Surface& Acts::Surface::operator=(const Surface& other) {
0156   if (&other != this) {
0157     GeometryObject::operator=(other);
0158     // detector element, identifier & layer association are unique
0159     m_transform = other.m_transform;
0160     m_associatedLayer = other.m_associatedLayer;
0161     m_surfaceMaterial = other.m_surfaceMaterial;
0162     m_associatedDetElement = other.m_associatedDetElement;
0163   }
0164   return *this;
0165 }
0166
0167 bool Acts::Surface::operator==(const Surface& other) const {
0168   // (a) fast exit for pointer comparison
0169   if (&other == this) {
0170     return true;
0171   }
0172   // (b) fast exit for type
0173   if (other.type() != type()) {
0174     return false;
0175   }
0176   // (c) fast exit for bounds
0177   if (other.bounds() != bounds()) {
0178     return false;
0179   }
0180   // (d) compare  detector elements
0181   if (m_associatedDetElement != other.m_associatedDetElement) {
0182     return false;
0183   }
0184   // (e) compare transform values
0185   if (!m_transform.isApprox(other.m_transform, 1e-9)) {
0186     return false;
0187   }
0188   // (f) compare material
0189   if (m_surfaceMaterial != other.m_surfaceMaterial) {
0190     return false;
0191   }
0192
0193   // we should be good
0194   return true;
0195 }
0196
0197 // overload dump for stream operator
0198 std::ostream& Acts::Surface::toStream(const GeometryContext& gctx,
0199                                       std::ostream& sl) const {
0200   sl << std::setiosflags(std::ios::fixed);
0201   sl << std::setprecision(4);
0202   sl << name() << std::endl;
0203   const Vector3& sfcenter = center(gctx);
0204   sl << "     Center position  (x, y, z) = (" << sfcenter.x() << ", "
0205      << sfcenter.y() << ", " << sfcenter.z() << ")" << std::endl;
0206   Acts::RotationMatrix3 rot(transform(gctx).matrix().block<3, 3>(0, 0));
0207   Acts::Vector3 rotX(rot.col(0));
0208   Acts::Vector3 rotY(rot.col(1));
0209   Acts::Vector3 rotZ(rot.col(2));
0210   sl << std::setprecision(6);
0211   sl << "     Rotation:             colX = (" << rotX(0) << ", " << rotX(1)
0212      << ", " << rotX(2) << ")" << std::endl;
0213   sl << "                           colY = (" << rotY(0) << ", " << rotY(1)
0214      << ", " << rotY(2) << ")" << std::endl;
0215   sl << "                           colZ = (" << rotZ(0) << ", " << rotZ(1)
0216      << ", " << rotZ(2) << ")" << std::endl;
0217   sl << "     Bounds  : " << bounds();
0218   sl << std::setprecision(-1);
0219   return sl;
0220 }
0221
0222 std::string Acts::Surface::toString(const GeometryContext& gctx) const {
0223   std::stringstream ss;
0224   toStream(gctx, ss);
0225   return ss.str();
0226 }
0227
0228 bool Acts::Surface::operator!=(const Acts::Surface& sf) const {
0229   return !(operator==(sf));
0230 }
0231
0232 Acts::Vector3 Acts::Surface::center(const GeometryContext& gctx) const {
0233   // fast access via transform matrix (and not translation())
0234   auto tMatrix = transform(gctx).matrix();
0235   return Vector3(tMatrix(0, 3), tMatrix(1, 3), tMatrix(2, 3));
0236 }
0237
0238 const Acts::Transform3& Acts::Surface::transform(
0239     const GeometryContext& gctx) const {
0240   if (m_associatedDetElement != nullptr) {
0241     return m_associatedDetElement->transform(gctx);
0242   }
0243   return m_transform;
0244 }
0245
0246 bool Acts::Surface::insideBounds(const Vector2& lposition,
0247                                  const BoundaryCheck& bcheck) const {
0248   return bounds().inside(lposition, bcheck);
0249 }
0250
0251 Acts::RotationMatrix3 Acts::Surface::referenceFrame(
0252     const GeometryContext& gctx, const Vector3& /*position*/,
0253     const Vector3& /*direction*/) const {
0254   return transform(gctx).matrix().block<3, 3>(0, 0);
0255 }
0256
0257 Acts::BoundToFreeMatrix Acts::Surface::boundToFreeJacobian(
0258     const GeometryContext& gctx, const Vector3& position,
0259     const Vector3& direction) const {
0260   assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
0261
0262   // retrieve the reference frame
0263   const auto rframe = referenceFrame(gctx, position, direction);
0264
0265   // Initialize the jacobian from local to global
0266   BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();
0267   // the local error components - given by reference frame
0268   jacToGlobal.topLeftCorner<3, 2>() = rframe.topLeftCorner<3, 2>();
0269   // the time component
0270   jacToGlobal(eFreeTime, eBoundTime) = 1;
0271   // the momentum components
0272   jacToGlobal.block<3, 2>(eFreeDir0, eBoundPhi) =
0273       sphericalToFreeDirectionJacobian(direction);
0274   jacToGlobal(eFreeQOverP, eBoundQOverP) = 1;
0275   return jacToGlobal;
0276 }
0277
0278 Acts::FreeToBoundMatrix Acts::Surface::freeToBoundJacobian(
0279     const GeometryContext& gctx, const Vector3& position,
0280     const Vector3& direction) const {
0281   assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
0282
0283   // The measurement frame of the surface
0284   RotationMatrix3 rframeT =
0285       referenceFrame(gctx, position, direction).transpose();
0286
0287   // Initialize the jacobian from global to local
0288   FreeToBoundMatrix jacToLocal = FreeToBoundMatrix::Zero();
0289   // Local position component given by the reference frame
0290   jacToLocal.block<2, 3>(eBoundLoc0, eFreePos0) = rframeT.block<2, 3>(0, 0);
0291   // Time component
0292   jacToLocal(eBoundTime, eFreeTime) = 1;
0293   // Directional and momentum elements for reference frame surface
0294   jacToLocal.block<2, 3>(eBoundPhi, eFreeDir0) =
0295       freeToSphericalDirectionJacobian(direction);
0296   jacToLocal(eBoundQOverP, eFreeQOverP) = 1;
0297   return jacToLocal;
0298 }
0299
0300 Acts::FreeToPathMatrix Acts::Surface::freeToPathDerivative(
0301     const GeometryContext& gctx, const Vector3& position,
0302     const Vector3& direction) const {
0303   assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
0304
0305   // The measurement frame of the surface
0306   const RotationMatrix3 rframe = referenceFrame(gctx, position, direction);
0307   // The measurement frame z axis
0308   const Vector3 refZAxis = rframe.col(2);
0309   // Cosine of angle between momentum direction and measurement frame z axis
0310   const ActsScalar dz = refZAxis.dot(direction);
0311   // Initialize the derivative
0312   FreeToPathMatrix freeToPath = FreeToPathMatrix::Zero();
0313   freeToPath.segment<3>(eFreePos0) = -1.0 * refZAxis.transpose() / dz;
0314   return freeToPath;
0315 }
0316
0317 const Acts::DetectorElementBase* Acts::Surface::associatedDetectorElement()
0318     const {
0319   return m_associatedDetElement;
0320 }
0321
0322 const Acts::Layer* Acts::Surface::associatedLayer() const {
0323   return m_associatedLayer;
0324 }
0325
0326 const Acts::ISurfaceMaterial* Acts::Surface::surfaceMaterial() const {
0327   return m_surfaceMaterial.get();
0328 }
0329
0330 const std::shared_ptr<const Acts::ISurfaceMaterial>&
0331 Acts::Surface::surfaceMaterialSharedPtr() const {
0332   return m_surfaceMaterial;
0333 }
0334
0335 void Acts::Surface::assignDetectorElement(
0336     const DetectorElementBase& detelement) {
0337   m_associatedDetElement = &detelement;
0338   // resetting the transform as it will be handled through the detector element
0339   // now
0340   m_transform = Transform3::Identity();
0341 }
0342
0343 void Acts::Surface::assignSurfaceMaterial(
0344     std::shared_ptr<const Acts::ISurfaceMaterial> material) {
0345   m_surfaceMaterial = std::move(material);
0346 }
0347
0348 void Acts::Surface::associateLayer(const Acts::Layer& lay) {
0349   m_associatedLayer = (&lay);
0350 }