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