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