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