src/Surfaces/ConeSurface.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/ConeSurface.hpp"
0010
0011 #include "Acts/Geometry/GeometryObject.hpp"
0012 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0013 #include "Acts/Surfaces/SurfaceError.hpp"
0014 #include "Acts/Surfaces/detail/AlignmentHelper.hpp"
0015 #include "Acts/Surfaces/detail/FacesHelper.hpp"
0016 #include "Acts/Surfaces/detail/VerticesHelper.hpp"
0017 #include "Acts/Utilities/AlgebraHelpers.hpp"
0018 #include "Acts/Utilities/Intersection.hpp"
0019 #include "Acts/Utilities/ThrowAssert.hpp"
0020 #include "Acts/Utilities/detail/RealQuadraticEquation.hpp"
0021
0022 #include <algorithm>
0023 #include <cmath>
0024 #include <limits>
0025 #include <numbers>
0026 #include <stdexcept>
0027 #include <utility>
0028 #include <vector>
0029
0030 namespace Acts {
0031
0032 using VectorHelpers::perp;
0033 using VectorHelpers::phi;
0034
0035 ConeSurface::ConeSurface(const ConeSurface& other)
0036     : GeometryObject{}, RegularSurface(other), m_bounds(other.m_bounds) {}
0037
0038 ConeSurface::ConeSurface(const GeometryContext& gctx, const ConeSurface& other,
0039                          const Transform3& shift)
0040     : RegularSurface(gctx, other, shift), m_bounds(other.m_bounds) {}
0041
0042 ConeSurface::ConeSurface(const Transform3& transform, double alpha,
0043                          bool symmetric)
0044     : RegularSurface(transform),
0045       m_bounds(std::make_shared<const ConeBounds>(alpha, symmetric)) {}
0046
0047 ConeSurface::ConeSurface(const Transform3& transform, double alpha, double zmin,
0048                          double zmax, double halfPhi)
0049     : RegularSurface(transform),
0050       m_bounds(std::make_shared<const ConeBounds>(alpha, zmin, zmax, halfPhi)) {
0051 }
0052
0053 ConeSurface::ConeSurface(const Transform3& transform,
0054                          std::shared_ptr<const ConeBounds> cbounds)
0055     : RegularSurface(transform), m_bounds(std::move(cbounds)) {
0056   throw_assert(m_bounds, "ConeBounds must not be nullptr");
0057 }
0058
0059 Vector3 ConeSurface::referencePosition(const GeometryContext& gctx,
0060                                        AxisDirection aDir) const {
0061   const Vector3& sfCenter = center(gctx);
0062
0063   // special binning type for R-type methods
0064   if (aDir == AxisDirection::AxisR || aDir == AxisDirection::AxisRPhi) {
0065     return Vector3(sfCenter.x() + bounds().r(sfCenter.z()), sfCenter.y(),
0066                    sfCenter.z());
0067   }
0068   // give the center as default for all of these binning types
0069   // AxisDirection::AxisX, AxisDirection::AxisY, AxisDirection::AxisZ,
0070   // AxisDirection::AxisR, AxisDirection::AxisPhi, AxisDirection::AxisRPhi,
0071   // AxisDirection::AxisTheta, AxisDirection::AxisEta
0072   return sfCenter;
0073 }
0074
0075 Surface::SurfaceType ConeSurface::type() const {
0076   return Surface::Cone;
0077 }
0078
0079 ConeSurface& ConeSurface::operator=(const ConeSurface& other) {
0080   if (this != &other) {
0081     Surface::operator=(other);
0082     m_bounds = other.m_bounds;
0083   }
0084   return *this;
0085 }
0086
0087 Vector3 ConeSurface::rotSymmetryAxis(const GeometryContext& gctx) const {
0088   return localToGlobalTransform(gctx).matrix().block<3, 1>(0, 2);
0089 }
0090
0091 RotationMatrix3 ConeSurface::referenceFrame(
0092     const GeometryContext& gctx, const Vector3& position,
0093     const Vector3& /*direction*/) const {
0094   RotationMatrix3 mFrame;
0095   // construct the measurement frame
0096   // measured Y is the local z axis
0097   Vector3 measY = rotSymmetryAxis(gctx);
0098   // measured z is the position transverse normalized
0099   Vector3 measDepth = Vector3(position.x(), position.y(), 0.).normalized();
0100   // measured X is what comoes out of it
0101   Vector3 measX(measY.cross(measDepth).normalized());
0102   // the columnes
0103   mFrame.col(0) = measX;
0104   mFrame.col(1) = measY;
0105   mFrame.col(2) = measDepth;
0106   // return the rotation matrix
0107   //!< @todo fold in alpha
0108   // return it
0109   return mFrame;
0110 }
0111
0112 Vector3 ConeSurface::localToGlobal(const GeometryContext& gctx,
0113                                    const Vector2& lposition) const {
0114   // create the position in the local 3d frame
0115   double r = lposition[1] * bounds().tanAlpha();
0116   double phi = lposition[0] / r;
0117   Vector3 loc3Dframe(r * std::cos(phi), r * std::sin(phi), lposition[1]);
0118   return localToGlobalTransform(gctx) * loc3Dframe;
0119 }
0120
0121 Result<Vector2> ConeSurface::globalToLocal(const GeometryContext& gctx,
0122                                            const Vector3& position,
0123                                            double tolerance) const {
0124   Vector3 loc3Dframe =
0125       inverseTransform(localToGlobalTransform(gctx)) * position;
0126   double r = loc3Dframe.z() * bounds().tanAlpha();
0127   if (std::abs(perp(loc3Dframe) - r) > tolerance) {
0128     return Result<Vector2>::failure(SurfaceError::GlobalPositionNotOnSurface);
0129   }
0130   return Result<Vector2>::success(
0131       Vector2(r * std::atan2(loc3Dframe.y(), loc3Dframe.x()), loc3Dframe.z()));
0132 }
0133
0134 double ConeSurface::pathCorrection(const GeometryContext& gctx,
0135                                    const Vector3& position,
0136                                    const Vector3& direction) const {
0137   // (cos phi cos alpha, sin phi cos alpha, sgn z sin alpha)
0138   Vector3 posLocal = inverseTransform(localToGlobalTransform(gctx)) * position;
0139   double phi = VectorHelpers::phi(posLocal);
0140   double sgn = -std::copysign(1., posLocal.z());
0141   double cosAlpha = std::cos(bounds().get(ConeBounds::eAlpha));
0142   double sinAlpha = std::sin(bounds().get(ConeBounds::eAlpha));
0143   Vector3 normalC(std::cos(phi) * cosAlpha, std::sin(phi) * cosAlpha,
0144                   sgn * sinAlpha);
0145   normalC = localToGlobalTransform(gctx).linear() * normalC;
0146   // Back to the global frame
0147   double cAlpha = normalC.dot(direction);
0148   return std::abs(1. / cAlpha);
0149 }
0150
0151 std::string ConeSurface::name() const {
0152   return "Acts::ConeSurface";
0153 }
0154
0155 Vector3 ConeSurface::normal(const GeometryContext& gctx,
0156                             const Vector2& lposition) const {
0157   // (cos phi cos alpha, sin phi cos alpha, sgn z sin alpha)
0158   double phi = lposition[0] / (bounds().r(lposition[1])),
0159          sgn = -std::copysign(1., lposition[1]);
0160   double cosAlpha = std::cos(bounds().get(ConeBounds::eAlpha));
0161   double sinAlpha = std::sin(bounds().get(ConeBounds::eAlpha));
0162   Vector3 localNormal(std::cos(phi) * cosAlpha, std::sin(phi) * cosAlpha,
0163                       sgn * sinAlpha);
0164   return Vector3(localToGlobalTransform(gctx).linear() * localNormal);
0165 }
0166
0167 Vector3 ConeSurface::normal(const GeometryContext& gctx,
0168                             const Vector3& position) const {
0169   // get it into the cylinder frame if needed
0170   // @todo respect opening angle
0171   Vector3 pos3D = inverseTransform(localToGlobalTransform(gctx)) * position;
0172   pos3D.z() = 0;
0173   return pos3D.normalized();
0174 }
0175
0176 const ConeBounds& ConeSurface::bounds() const {
0177   // is safe because no constructor w/o bounds exists
0178   return *m_bounds;
0179 }
0180
0181 Polyhedron ConeSurface::polyhedronRepresentation(
0182     const GeometryContext& gctx, unsigned int quarterSegments) const {
0183   // Prepare vertices and faces
0184   std::vector<Vector3> vertices;
0185   std::vector<Polyhedron::FaceType> faces;
0186   std::vector<Polyhedron::FaceType> triangularMesh;
0187   double minZ = bounds().get(ConeBounds::eMinZ);
0188   double maxZ = bounds().get(ConeBounds::eMaxZ);
0189
0190   if (minZ == -std::numeric_limits<double>::infinity() ||
0191       maxZ == std::numeric_limits<double>::infinity()) {
0192     throw std::domain_error(
0193         "Polyhedron representation of boundless surface is not possible");
0194   }
0195
0196   auto ctransform = localToGlobalTransform(gctx);
0197
0198   // The tip - created only once and only, if it is not a cut-off cone
0199   bool tipExists = false;
0200   if (minZ * maxZ <= s_onSurfaceTolerance) {
0201     vertices.push_back(ctransform * Vector3(0., 0., 0.));
0202     tipExists = true;
0203   }
0204
0205   // Cone parameters
0206   double hPhiSec = bounds().get(ConeBounds::eHalfPhiSector);
0207   double avgPhi = bounds().get(ConeBounds::eAveragePhi);
0208   std::vector<double> refPhi = {};
0209   if (bool fullCone =
0210           std::abs(hPhiSec - std::numbers::pi) < s_fullAzimuthTolerance;
0211       !fullCone) {
0212     refPhi = {avgPhi};
0213   }
0214
0215   // Add the cone sizes
0216   std::vector<double> coneSides;
0217   if (std::abs(minZ) > s_onSurfaceTolerance) {
0218     coneSides.push_back(minZ);
0219   }
0220   if (std::abs(maxZ) > s_onSurfaceTolerance) {
0221     coneSides.push_back(maxZ);
0222   }
0223
0224   for (auto& z : coneSides) {
0225     std::size_t firstIv = vertices.size();
0226     // Radius and z offset
0227     double r = std::abs(z) * bounds().tanAlpha();
0228     Vector3 zoffset(0., 0., z);
0229     auto svertices = detail::VerticesHelper::segmentVertices(
0230         {r, r}, avgPhi - hPhiSec, avgPhi + hPhiSec, refPhi, quarterSegments,
0231         zoffset, ctransform);
0232     vertices.insert(vertices.end(), svertices.begin(), svertices.end());
0233     // If the tip exists, the faces need to be triangular
0234     if (tipExists) {
0235       for (std::size_t iv = firstIv + 1; iv < svertices.size() + firstIv;
0236            ++iv) {
0237         std::size_t one = 0, two = iv, three = iv - 1;
0238         if (z < 0.) {
0239           std::swap(two, three);
0240         }
0241         faces.push_back({one, two, three});
0242       }
0243     }
0244   }
0245
0246   // if no tip exists, connect the two bows
0247   if (tipExists) {
0248     triangularMesh = faces;
0249   } else {
0250     auto facesMesh = detail::FacesHelper::cylindricalFaceMesh(vertices);
0251     faces = facesMesh.first;
0252     triangularMesh = facesMesh.second;
0253   }
0254
0255   return Polyhedron(vertices, faces, triangularMesh, false);
0256 }
0257
0258 detail::RealQuadraticEquation ConeSurface::intersectionSolver(
0259     const GeometryContext& gctx, const Vector3& position,
0260     const Vector3& direction) const {
0261   // Transform into the local frame
0262   Transform3 invTrans = inverseTransform(localToGlobalTransform(gctx));
0263   Vector3 point1 = invTrans * position;
0264   Vector3 dir1 = invTrans.linear() * direction;
0265
0266   // See file header for the formula derivation
0267   double tan2Alpha = bounds().tanAlpha() * bounds().tanAlpha(),
0268          A = dir1.x() * dir1.x() + dir1.y() * dir1.y() -
0269              tan2Alpha * dir1.z() * dir1.z(),
0270          B = 2 * (dir1.x() * point1.x() + dir1.y() * point1.y() -
0271                   tan2Alpha * dir1.z() * point1.z()),
0272          C = point1.x() * point1.x() + point1.y() * point1.y() -
0273              tan2Alpha * point1.z() * point1.z();
0274   if (A == 0.) {
0275     A += 1e-16;  // avoid division by zero
0276   }
0277
0278   return detail::RealQuadraticEquation(A, B, C);
0279 }
0280
0281 MultiIntersection3D ConeSurface::intersect(
0282     const GeometryContext& gctx, const Vector3& position,
0283     const Vector3& direction, const BoundaryTolerance& boundaryTolerance,
0284     double tolerance) const {
0285   // Solve the quadratic equation
0286   auto qe = intersectionSolver(gctx, position, direction);
0287
0288   // If no valid solution return a non-valid surfaceIntersection
0289   if (qe.solutions == 0) {
0290     return MultiIntersection3D(Intersection3D::Invalid(),
0291                                Intersection3D::Invalid());
0292   }
0293
0294   // Check the validity of the first solution
0295   Vector3 solution1 = position + qe.first * direction;
0296   IntersectionStatus status1 = std::abs(qe.first) < std::abs(tolerance)
0297                                    ? IntersectionStatus::onSurface
0298                                    : IntersectionStatus::reachable;
0299
0300   if (!boundaryTolerance.isInfinite() &&
0301       !isOnSurface(gctx, solution1, direction, boundaryTolerance)) {
0302     status1 = IntersectionStatus::unreachable;
0303   }
0304
0305   // Check the validity of the second solution
0306   Vector3 solution2 = position + qe.first * direction;
0307   IntersectionStatus status2 = std::abs(qe.second) < std::abs(tolerance)
0308                                    ? IntersectionStatus::onSurface
0309                                    : IntersectionStatus::reachable;
0310   if (!boundaryTolerance.isInfinite() &&
0311       !isOnSurface(gctx, solution2, direction, boundaryTolerance)) {
0312     status2 = IntersectionStatus::unreachable;
0313   }
0314
0315   const auto& tf = localToGlobalTransform(gctx);
0316   // Set the intersection
0317   Intersection3D first(tf * solution1, qe.first, status1);
0318   Intersection3D second(tf * solution2, qe.second, status2);
0319   // Order based on path length
0320   if (first.pathLength() <= second.pathLength()) {
0321     return MultiIntersection3D(first, second);
0322   }
0323   return MultiIntersection3D(second, first);
0324 }
0325
0326 AlignmentToPathMatrix ConeSurface::alignmentToPathDerivative(
0327     const GeometryContext& gctx, const Vector3& position,
0328     const Vector3& direction) const {
0329   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0330
0331   // The vector between position and center
0332   const auto pcRowVec = (position - center(gctx)).transpose().eval();
0333   // The rotation
0334   const auto& rotation = localToGlobalTransform(gctx).rotation();
0335   // The local frame x/y/z axis
0336   const auto& localXAxis = rotation.col(0);
0337   const auto& localYAxis = rotation.col(1);
0338   const auto& localZAxis = rotation.col(2);
0339   // The local coordinates
0340   const auto localPos = (rotation.transpose() * position).eval();
0341   const auto dx = direction.dot(localXAxis);
0342   const auto dy = direction.dot(localYAxis);
0343   const auto dz = direction.dot(localZAxis);
0344   // The normalization factor
0345   const auto tanAlpha2 = bounds().tanAlpha() * bounds().tanAlpha();
0346   const auto norm = 1 / (1 - dz * dz * (1 + tanAlpha2));
0347   // The direction transpose
0348   const auto& dirRowVec = direction.transpose();
0349   // The derivative of path w.r.t. the local axes
0350   // @note The following calculations assume that the intersection of the track
0351   // with the cone always satisfy: localPos.z()*tanAlpha =perp(localPos)
0352   const auto localXAxisToPath =
0353       (-2 * norm * (dx * pcRowVec + localPos.x() * dirRowVec)).eval();
0354   const auto localYAxisToPath =
0355       (-2 * norm * (dy * pcRowVec + localPos.y() * dirRowVec)).eval();
0356   const auto localZAxisToPath =
0357       (2 * norm * tanAlpha2 * (dz * pcRowVec + localPos.z() * dirRowVec) -
0358        4 * norm * norm * (1 + tanAlpha2) *
0359            (dx * localPos.x() + dy * localPos.y() -
0360             dz * localPos.z() * tanAlpha2) *
0361            dz * dirRowVec)
0362           .eval();
0363   // Calculate the derivative of local frame axes w.r.t its rotation
0364   const auto [rotToLocalXAxis, rotToLocalYAxis, rotToLocalZAxis] =
0365       detail::rotationToLocalAxesDerivative(rotation);
0366   // Initialize the derivative of propagation path w.r.t. local frame
0367   // translation (origin) and rotation
0368   AlignmentToPathMatrix alignToPath = AlignmentToPathMatrix::Zero();
0369   alignToPath.segment<3>(eAlignmentCenter0) =
0370       2 * norm * (dx * localXAxis.transpose() + dy * localYAxis.transpose());
0371   alignToPath.segment<3>(eAlignmentRotation0) =
0372       localXAxisToPath * rotToLocalXAxis + localYAxisToPath * rotToLocalYAxis +
0373       localZAxisToPath * rotToLocalZAxis;
0374
0375   return alignToPath;
0376 }
0377
0378 Matrix<2, 3> ConeSurface::localCartesianToBoundLocalDerivative(
0379     const GeometryContext& gctx, const Vector3& position) const {
0380   using VectorHelpers::perp;
0381   using VectorHelpers::phi;
0382   // The local frame transform
0383   const auto& sTransform = localToGlobalTransform(gctx);
0384   // calculate the transformation to local coordinates
0385   const Vector3 localPos = inverseTransform(sTransform) * position;
0386   const double lr = perp(localPos);
0387   const double lphi = phi(localPos);
0388   const double lcphi = std::cos(lphi);
0389   const double lsphi = std::sin(lphi);
0390   // Solve for radius R
0391   const double R = localPos.z() * bounds().tanAlpha();
0392   Matrix<2, 3> loc3DToLocBound = Matrix<2, 3>::Zero();
0393   loc3DToLocBound << -R * lsphi / lr, R * lcphi / lr,
0394       lphi * bounds().tanAlpha(), 0, 0, 1;
0395
0396   return loc3DToLocBound;
0397 }
0398
0399 const std::shared_ptr<const ConeBounds>& ConeSurface::boundsPtr() const {
0400   return m_bounds;
0401 }
0402
0403 void ConeSurface::assignSurfaceBounds(
0404     std::shared_ptr<const ConeBounds> newBounds) {
0405   m_bounds = std::move(newBounds);
0406 }
0407
0408 }  // namespace Acts