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