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