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 std::array<std::string, Surface::SurfaceType::Other>
0023     Surface::s_surfaceTypeNames = {"Cone",  "Cylinder", "Disc",       "Perigee",
0024                                    "Plane", "Straw",    "Curvilinear"};
0025
0026 Surface::Surface(const Transform3& transform)
0027     : GeometryObject(), m_transform(std::make_unique<Transform3>(transform)) {}
0028
0029 Surface::Surface(const DetectorElementBase& detelement)
0030     : GeometryObject(), m_associatedDetElement(&detelement) {}
0031
0032 Surface::Surface(const Surface& other)
0033     : GeometryObject(other),
0034       std::enable_shared_from_this<Surface>(),
0035       m_associatedDetElement(other.m_associatedDetElement),
0036       m_surfaceMaterial(other.m_surfaceMaterial) {
0037   if (other.m_transform) {
0038     m_transform = std::make_unique<Transform3>(*other.m_transform);
0039   }
0040 }
0041
0042 Surface::Surface(const GeometryContext& gctx, const Surface& other,
0043                  const Transform3& shift)
0044     : GeometryObject(),
0045       m_transform(std::make_unique<Transform3>(shift * other.transform(gctx))),
0046       m_surfaceMaterial(other.m_surfaceMaterial) {}
0047
0048 Surface::~Surface() = default;
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 bool Surface::insideBounds(const Vector2& lposition,
0250                            const BoundaryTolerance& boundaryTolerance) const {
0251   return bounds().inside(lposition, boundaryTolerance);
0252 }
0253
0254 RotationMatrix3 Surface::referenceFrame(const GeometryContext& gctx,
0255                                         const Vector3& /*position*/,
0256                                         const Vector3& /*direction*/) const {
0257   return transform(gctx).matrix().block<3, 3>(0, 0);
0258 }
0259
0260 BoundToFreeMatrix Surface::boundToFreeJacobian(const GeometryContext& gctx,
0261                                                const Vector3& position,
0262                                                const Vector3& direction) const {
0263   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0264
0265   // retrieve the reference frame
0266   const auto rframe = referenceFrame(gctx, position, direction);
0267
0268   // Initialize the jacobian from local to global
0269   BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();
0270   // the local error components - given by reference frame
0271   jacToGlobal.topLeftCorner<3, 2>() = rframe.topLeftCorner<3, 2>();
0272   // the time component
0273   jacToGlobal(eFreeTime, eBoundTime) = 1;
0274   // the momentum components
0275   jacToGlobal.block<3, 2>(eFreeDir0, eBoundPhi) =
0276       sphericalToFreeDirectionJacobian(direction);
0277   jacToGlobal(eFreeQOverP, eBoundQOverP) = 1;
0278   return jacToGlobal;
0279 }
0280
0281 FreeToBoundMatrix Surface::freeToBoundJacobian(const GeometryContext& gctx,
0282                                                const Vector3& position,
0283                                                const Vector3& direction) const {
0284   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0285
0286   // The measurement frame of the surface
0287   RotationMatrix3 rframeT =
0288       referenceFrame(gctx, position, direction).transpose();
0289
0290   // Initialize the jacobian from global to local
0291   FreeToBoundMatrix jacToLocal = FreeToBoundMatrix::Zero();
0292   // Local position component given by the reference frame
0293   jacToLocal.block<2, 3>(eBoundLoc0, eFreePos0) = rframeT.block<2, 3>(0, 0);
0294   // Time component
0295   jacToLocal(eBoundTime, eFreeTime) = 1;
0296   // Directional and momentum elements for reference frame surface
0297   jacToLocal.block<2, 3>(eBoundPhi, eFreeDir0) =
0298       freeToSphericalDirectionJacobian(direction);
0299   jacToLocal(eBoundQOverP, eFreeQOverP) = 1;
0300   return jacToLocal;
0301 }
0302
0303 FreeToPathMatrix Surface::freeToPathDerivative(const GeometryContext& gctx,
0304                                                const Vector3& position,
0305                                                const Vector3& direction) const {
0306   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0307
0308   // The measurement frame of the surface
0309   const RotationMatrix3 rframe = referenceFrame(gctx, position, direction);
0310   // The measurement frame z axis
0311   const Vector3 refZAxis = rframe.col(2);
0312   // Cosine of angle between momentum direction and measurement frame z axis
0313   const double dz = refZAxis.dot(direction);
0314   // Initialize the derivative
0315   FreeToPathMatrix freeToPath = FreeToPathMatrix::Zero();
0316   freeToPath.segment<3>(eFreePos0) = -1.0 * refZAxis.transpose() / dz;
0317   return freeToPath;
0318 }
0319
0320 const DetectorElementBase* Surface::associatedDetectorElement() const {
0321   return m_associatedDetElement;
0322 }
0323
0324 const Layer* Surface::associatedLayer() const {
0325   return m_associatedLayer;
0326 }
0327
0328 const ISurfaceMaterial* Surface::surfaceMaterial() const {
0329   return m_surfaceMaterial.get();
0330 }
0331
0332 const std::shared_ptr<const ISurfaceMaterial>&
0333 Surface::surfaceMaterialSharedPtr() const {
0334   return m_surfaceMaterial;
0335 }
0336
0337 void Surface::assignDetectorElement(const DetectorElementBase& detelement) {
0338   m_associatedDetElement = &detelement;
0339   // resetting the transform as it will be handled through the detector element
0340   // now
0341   m_transform.reset();
0342 }
0343
0344 void Surface::assignSurfaceMaterial(
0345     std::shared_ptr<const ISurfaceMaterial> material) {
0346   m_surfaceMaterial = std::move(material);
0347 }
0348
0349 void Surface::associateLayer(const Layer& lay) {
0350   m_associatedLayer = (&lay);
0351 }
0352
0353 void Surface::visualize(IVisualization3D& helper, const GeometryContext& gctx,
0354                         const ViewConfig& viewConfig) const {
0355   Polyhedron polyhedron =
0356       polyhedronRepresentation(gctx, viewConfig.quarterSegments);
0357   polyhedron.visualize(helper, viewConfig);
0358 }
0359
0360 }  // namespace Acts