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