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 * cos(phi), r * 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 * 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 = posLocal.z() > 0. ? -1. : +1.;
0145   double cosAlpha = std::cos(bounds().get(ConeBounds::eAlpha));
0146   double sinAlpha = std::sin(bounds().get(ConeBounds::eAlpha));
0147   Vector3 normalC(cos(phi) * cosAlpha, sin(phi) * cosAlpha, sgn * sinAlpha);
0148   normalC = transform(gctx) * normalC;
0149   // Back to the global frame
0150   double cAlpha = normalC.dot(direction);
0151   return std::abs(1. / cAlpha);
0152 }
0153
0154 std::string ConeSurface::name() const {
0155   return "Acts::ConeSurface";
0156 }
0157
0158 Vector3 ConeSurface::normal(const GeometryContext& gctx,
0159                             const Vector2& lposition) const {
0160   // (cos phi cos alpha, sin phi cos alpha, sgn z sin alpha)
0161   double phi = lposition[0] / (bounds().r(lposition[1])),
0162          sgn = lposition[1] > 0 ? -1. : +1.;
0163   double cosAlpha = std::cos(bounds().get(ConeBounds::eAlpha));
0164   double sinAlpha = std::sin(bounds().get(ConeBounds::eAlpha));
0165   Vector3 localNormal(cos(phi) * cosAlpha, sin(phi) * cosAlpha, sgn * sinAlpha);
0166   return Vector3(transform(gctx).linear() * localNormal);
0167 }
0168
0169 Vector3 ConeSurface::normal(const GeometryContext& gctx,
0170                             const Vector3& position) const {
0171   // get it into the cylinder frame if needed
0172   // @todo respect opening angle
0173   Vector3 pos3D = transform(gctx).inverse() * position;
0174   pos3D.z() = 0;
0175   return pos3D.normalized();
0176 }
0177
0178 const ConeBounds& ConeSurface::bounds() const {
0179   // is safe because no constructor w/o bounds exists
0180   return (*m_bounds.get());
0181 }
0182
0183 Polyhedron ConeSurface::polyhedronRepresentation(
0184     const GeometryContext& gctx, unsigned int quarterSegments) const {
0185   // Prepare vertices and faces
0186   std::vector<Vector3> vertices;
0187   std::vector<Polyhedron::FaceType> faces;
0188   std::vector<Polyhedron::FaceType> triangularMesh;
0189   double minZ = bounds().get(ConeBounds::eMinZ);
0190   double maxZ = bounds().get(ConeBounds::eMaxZ);
0191
0192   if (minZ == -std::numeric_limits<double>::infinity() ||
0193       maxZ == std::numeric_limits<double>::infinity()) {
0194     throw std::domain_error(
0195         "Polyhedron representation of boundless surface is not possible");
0196   }
0197
0198   auto ctransform = transform(gctx);
0199
0200   // The tip - created only once and only, if it is not a cut-off cone
0201   bool tipExists = false;
0202   if (minZ * maxZ <= s_onSurfaceTolerance) {
0203     vertices.push_back(ctransform * Vector3(0., 0., 0.));
0204     tipExists = true;
0205   }
0206
0207   // Cone parameters
0208   double hPhiSec = bounds().get(ConeBounds::eHalfPhiSector);
0209   double avgPhi = bounds().get(ConeBounds::eAveragePhi);
0210   std::vector<double> refPhi = {};
0211   if (bool fullCone = (hPhiSec == std::numbers::pi); !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 = transform(gctx).inverse();
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 SurfaceMultiIntersection 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 {{Intersection3D::invalid(), Intersection3D::invalid()}, this};
0291   }
0292
0293   // Check the validity of the first solution
0294   Vector3 solution1 = position + qe.first * direction;
0295   IntersectionStatus status1 = std::abs(qe.first) < std::abs(tolerance)
0296                                    ? IntersectionStatus::onSurface
0297                                    : IntersectionStatus::reachable;
0298
0299   if (!boundaryTolerance.isInfinite() &&
0300       !isOnSurface(gctx, solution1, direction, boundaryTolerance)) {
0301     status1 = IntersectionStatus::unreachable;
0302   }
0303
0304   // Check the validity of the second solution
0305   Vector3 solution2 = position + qe.first * direction;
0306   IntersectionStatus status2 = std::abs(qe.second) < std::abs(tolerance)
0307                                    ? IntersectionStatus::onSurface
0308                                    : IntersectionStatus::reachable;
0309   if (!boundaryTolerance.isInfinite() &&
0310       !isOnSurface(gctx, solution2, direction, boundaryTolerance)) {
0311     status2 = IntersectionStatus::unreachable;
0312   }
0313
0314   const auto& tf = transform(gctx);
0315   // Set the intersection
0316   Intersection3D first(tf * solution1, qe.first, status1);
0317   Intersection3D second(tf * solution2, qe.second, status2);
0318   // Order based on path length
0319   if (first.pathLength() <= second.pathLength()) {
0320     return {{first, second}, this};
0321   }
0322   return {{second, first}, this};
0323 }
0324
0325 AlignmentToPathMatrix ConeSurface::alignmentToPathDerivative(
0326     const GeometryContext& gctx, const Vector3& position,
0327     const Vector3& direction) const {
0328   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0329
0330   // The vector between position and center
0331   const auto pcRowVec = (position - center(gctx)).transpose().eval();
0332   // The rotation
0333   const auto& rotation = transform(gctx).rotation();
0334   // The local frame x/y/z axis
0335   const auto& localXAxis = rotation.col(0);
0336   const auto& localYAxis = rotation.col(1);
0337   const auto& localZAxis = rotation.col(2);
0338   // The local coordinates
0339   const auto localPos = (rotation.transpose() * position).eval();
0340   const auto dx = direction.dot(localXAxis);
0341   const auto dy = direction.dot(localYAxis);
0342   const auto dz = direction.dot(localZAxis);
0343   // The normalization factor
0344   const auto tanAlpha2 = bounds().tanAlpha() * bounds().tanAlpha();
0345   const auto norm = 1 / (1 - dz * dz * (1 + tanAlpha2));
0346   // The direction transpose
0347   const auto& dirRowVec = direction.transpose();
0348   // The derivative of path w.r.t. the local axes
0349   // @note The following calculations assume that the intersection of the track
0350   // with the cone always satisfy: localPos.z()*tanAlpha =perp(localPos)
0351   const auto localXAxisToPath =
0352       (-2 * norm * (dx * pcRowVec + localPos.x() * dirRowVec)).eval();
0353   const auto localYAxisToPath =
0354       (-2 * norm * (dy * pcRowVec + localPos.y() * dirRowVec)).eval();
0355   const auto localZAxisToPath =
0356       (2 * norm * tanAlpha2 * (dz * pcRowVec + localPos.z() * dirRowVec) -
0357        4 * norm * norm * (1 + tanAlpha2) *
0358            (dx * localPos.x() + dy * localPos.y() -
0359             dz * localPos.z() * tanAlpha2) *
0360            dz * dirRowVec)
0361           .eval();
0362   // Calculate the derivative of local frame axes w.r.t its rotation
0363   const auto [rotToLocalXAxis, rotToLocalYAxis, rotToLocalZAxis] =
0364       detail::rotationToLocalAxesDerivative(rotation);
0365   // Initialize the derivative of propagation path w.r.t. local frame
0366   // translation (origin) and rotation
0367   AlignmentToPathMatrix alignToPath = AlignmentToPathMatrix::Zero();
0368   alignToPath.segment<3>(eAlignmentCenter0) =
0369       2 * norm * (dx * localXAxis.transpose() + dy * localYAxis.transpose());
0370   alignToPath.segment<3>(eAlignmentRotation0) =
0371       localXAxisToPath * rotToLocalXAxis + localYAxisToPath * rotToLocalYAxis +
0372       localZAxisToPath * rotToLocalZAxis;
0373
0374   return alignToPath;
0375 }
0376
0377 ActsMatrix<2, 3> ConeSurface::localCartesianToBoundLocalDerivative(
0378     const GeometryContext& gctx, const Vector3& position) const {
0379   using VectorHelpers::perp;
0380   using VectorHelpers::phi;
0381   // The local frame transform
0382   const auto& sTransform = transform(gctx);
0383   // calculate the transformation to local coordinates
0384   const Vector3 localPos = sTransform.inverse() * position;
0385   const double lr = perp(localPos);
0386   const double lphi = phi(localPos);
0387   const double lcphi = std::cos(lphi);
0388   const double lsphi = std::sin(lphi);
0389   // Solve for radius R
0390   const double R = localPos.z() * bounds().tanAlpha();
0391   ActsMatrix<2, 3> loc3DToLocBound = ActsMatrix<2, 3>::Zero();
0392   loc3DToLocBound << -R * lsphi / lr, R * lcphi / lr,
0393       lphi * bounds().tanAlpha(), 0, 0, 1;
0394
0395   return loc3DToLocBound;
0396 }
0397
0398 }  // namespace Acts