src/Surfaces/CylinderSurface.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2016-2020 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Surfaces/CylinderSurface.hpp"
0010
0011 #include "Acts/Geometry/GeometryObject.hpp"
0012 #include "Acts/Surfaces/SurfaceError.hpp"
0013 #include "Acts/Surfaces/detail/AlignmentHelper.hpp"
0014 #include "Acts/Surfaces/detail/FacesHelper.hpp"
0015 #include "Acts/Utilities/Helpers.hpp"
0016 #include "Acts/Utilities/Intersection.hpp"
0017 #include "Acts/Utilities/ThrowAssert.hpp"
0018
0019 #include <algorithm>
0020 #include <cassert>
0021 #include <cmath>
0022 #include <utility>
0023 #include <vector>
0024
0025 namespace Acts {
0026 class DetectorElementBase;
0027 }  // namespace Acts
0028
0029 using Acts::VectorHelpers::perp;
0030 using Acts::VectorHelpers::phi;
0031
0032 Acts::CylinderSurface::CylinderSurface(const CylinderSurface& other)
0033     : GeometryObject(), RegularSurface(other), m_bounds(other.m_bounds) {}
0034
0035 Acts::CylinderSurface::CylinderSurface(const GeometryContext& gctx,
0036                                        const CylinderSurface& other,
0037                                        const Transform3& shift)
0038     : GeometryObject(),
0039       RegularSurface(gctx, other, shift),
0040       m_bounds(other.m_bounds) {}
0041
0042 Acts::CylinderSurface::CylinderSurface(const Transform3& transform,
0043                                        double radius, double halfz,
0044                                        double halfphi, double avphi,
0045                                        double bevelMinZ, double bevelMaxZ)
0046     : RegularSurface(transform),
0047       m_bounds(std::make_shared<const CylinderBounds>(
0048           radius, halfz, halfphi, avphi, bevelMinZ, bevelMaxZ)) {}
0049
0050 Acts::CylinderSurface::CylinderSurface(
0051     std::shared_ptr<const CylinderBounds> cbounds,
0052     const DetectorElementBase& detelement)
0053     : RegularSurface(detelement), m_bounds(std::move(cbounds)) {
0054   /// surfaces representing a detector element must have bounds
0055   throw_assert(m_bounds, "CylinderBounds must not be nullptr");
0056 }
0057
0058 Acts::CylinderSurface::CylinderSurface(
0059     const Transform3& transform, std::shared_ptr<const CylinderBounds> cbounds)
0060     : RegularSurface(transform), m_bounds(std::move(cbounds)) {
0061   throw_assert(m_bounds, "CylinderBounds must not be nullptr");
0062 }
0063
0064 Acts::CylinderSurface& Acts::CylinderSurface::operator=(
0065     const CylinderSurface& other) {
0066   if (this != &other) {
0067     Surface::operator=(other);
0068     m_bounds = other.m_bounds;
0069   }
0070   return *this;
0071 }
0072
0073 // return the binning position for ordering in the BinnedArray
0074 Acts::Vector3 Acts::CylinderSurface::binningPosition(
0075     const GeometryContext& gctx, BinningValue bValue) const {
0076   // special binning type for R-type methods
0077   if (bValue == Acts::binR || bValue == Acts::binRPhi) {
0078     double R = bounds().get(CylinderBounds::eR);
0079     double phi = bounds().get(CylinderBounds::eAveragePhi);
0080     return localToGlobal(gctx, Vector2{phi * R, 0}, Vector3{});
0081   }
0082   // give the center as default for all of these binning types
0083   // binX, binY, binZ, binR, binPhi, binRPhi, binH, binEta
0084   return center(gctx);
0085 }
0086
0087 // return the measurement frame: it's the tangential plane
0088 Acts::RotationMatrix3 Acts::CylinderSurface::referenceFrame(
0089     const GeometryContext& gctx, const Vector3& position,
0090     const Vector3& /*direction*/) const {
0091   RotationMatrix3 mFrame;
0092   // construct the measurement frame
0093   // measured Y is the z axis
0094   Vector3 measY = rotSymmetryAxis(gctx);
0095   // measured z is the position normalized transverse (in local)
0096   Vector3 measDepth = normal(gctx, position);
0097   // measured X is what comoes out of it
0098   Vector3 measX(measY.cross(measDepth).normalized());
0099   // assign the columnes
0100   mFrame.col(0) = measX;
0101   mFrame.col(1) = measY;
0102   mFrame.col(2) = measDepth;
0103   // return the rotation matrix
0104   return mFrame;
0105 }
0106
0107 Acts::Surface::SurfaceType Acts::CylinderSurface::type() const {
0108   return Surface::Cylinder;
0109 }
0110
0111 Acts::Vector3 Acts::CylinderSurface::localToGlobal(
0112     const GeometryContext& gctx, const Vector2& lposition) const {
0113   // create the position in the local 3d frame
0114   double r = bounds().get(CylinderBounds::eR);
0115   double phi = lposition[Acts::eBoundLoc0] / r;
0116   Vector3 position(r * cos(phi), r * sin(phi), lposition[Acts::eBoundLoc1]);
0117   return transform(gctx) * position;
0118 }
0119
0120 Acts::Result<Acts::Vector2> Acts::CylinderSurface::globalToLocal(
0121     const GeometryContext& gctx, const Vector3& position,
0122     double tolerance) const {
0123   double inttol = tolerance;
0124   if (tolerance == s_onSurfaceTolerance) {
0125     // transform default value!
0126     // @TODO: check if s_onSurfaceTolerance would do here
0127     inttol = bounds().get(CylinderBounds::eR) * 0.0001;
0128   }
0129   if (inttol < 0.01) {
0130     inttol = 0.01;
0131   }
0132   const Transform3& sfTransform = transform(gctx);
0133   Transform3 inverseTrans(sfTransform.inverse());
0134   Vector3 loc3Dframe(inverseTrans * position);
0135   if (std::abs(perp(loc3Dframe) - bounds().get(CylinderBounds::eR)) > inttol) {
0136     return Result<Vector2>::failure(SurfaceError::GlobalPositionNotOnSurface);
0137   }
0138   return Result<Vector2>::success(
0139       {bounds().get(CylinderBounds::eR) * phi(loc3Dframe), loc3Dframe.z()});
0140 }
0141
0142 std::string Acts::CylinderSurface::name() const {
0143   return "Acts::CylinderSurface";
0144 }
0145
0146 Acts::Vector3 Acts::CylinderSurface::normal(
0147     const GeometryContext& gctx, const Acts::Vector2& lposition) const {
0148   double phi = lposition[Acts::eBoundLoc0] / m_bounds->get(CylinderBounds::eR);
0149   Vector3 localNormal(cos(phi), sin(phi), 0.);
0150   return transform(gctx).linear() * localNormal;
0151 }
0152
0153 Acts::Vector3 Acts::CylinderSurface::normal(
0154     const GeometryContext& gctx, const Acts::Vector3& position) const {
0155   const Transform3& sfTransform = transform(gctx);
0156   // get it into the cylinder frame
0157   Vector3 pos3D = sfTransform.inverse() * position;
0158   // set the z coordinate to 0
0159   pos3D.z() = 0.;
0160   // normalize and rotate back into global
0161   return sfTransform.linear() * pos3D.normalized();
0162 }
0163
0164 double Acts::CylinderSurface::pathCorrection(
0165     const GeometryContext& gctx, const Acts::Vector3& position,
0166     const Acts::Vector3& direction) const {
0167   Vector3 normalT = normal(gctx, position);
0168   double cosAlpha = normalT.dot(direction);
0169   return std::fabs(1. / cosAlpha);
0170 }
0171
0172 const Acts::CylinderBounds& Acts::CylinderSurface::bounds() const {
0173   return (*m_bounds.get());
0174 }
0175
0176 Acts::Polyhedron Acts::CylinderSurface::polyhedronRepresentation(
0177     const GeometryContext& gctx, std::size_t lseg) const {
0178   auto ctrans = transform(gctx);
0179
0180   // Prepare vertices and faces
0181   std::vector<Vector3> vertices = bounds().createCircles(ctrans, lseg);
0182   std::vector<Polyhedron::FaceType> faces;
0183   std::vector<Polyhedron::FaceType> triangularMesh;
0184
0185   bool fullCylinder = bounds().coversFullAzimuth();
0186
0187   auto facesMesh =
0188       detail::FacesHelper::cylindricalFaceMesh(vertices, fullCylinder);
0189   return Polyhedron(vertices, facesMesh.first, facesMesh.second, false);
0190 }
0191
0192 Acts::Vector3 Acts::CylinderSurface::rotSymmetryAxis(
0193     const GeometryContext& gctx) const {
0194   // fast access via transform matrix (and not rotation())
0195   return transform(gctx).matrix().block<3, 1>(0, 2);
0196 }
0197
0198 Acts::detail::RealQuadraticEquation Acts::CylinderSurface::intersectionSolver(
0199     const Transform3& transform, const Vector3& position,
0200     const Vector3& direction) const {
0201   // Solve for radius R
0202   double R = bounds().get(CylinderBounds::eR);
0203
0204   // Get the transformation matrtix
0205   const auto& tMatrix = transform.matrix();
0206   Vector3 caxis = tMatrix.block<3, 1>(0, 2).transpose();
0207   Vector3 ccenter = tMatrix.block<3, 1>(0, 3).transpose();
0208
0209   // Check documentation for explanation
0210   Vector3 pc = position - ccenter;
0211   Vector3 pcXcd = pc.cross(caxis);
0212   Vector3 ldXcd = direction.cross(caxis);
0213   double a = ldXcd.dot(ldXcd);
0214   double b = 2. * (ldXcd.dot(pcXcd));
0215   double c = pcXcd.dot(pcXcd) - (R * R);
0216   // And solve the qaudratic equation
0217   return detail::RealQuadraticEquation(a, b, c);
0218 }
0219
0220 Acts::SurfaceMultiIntersection Acts::CylinderSurface::intersect(
0221     const GeometryContext& gctx, const Vector3& position,
0222     const Vector3& direction, const BoundaryCheck& bcheck,
0223     ActsScalar tolerance) const {
0224   const auto& gctxTransform = transform(gctx);
0225
0226   // Solve the quadratic equation
0227   auto qe = intersectionSolver(gctxTransform, position, direction);
0228
0229   // If no valid solution return a non-valid surfaceIntersection
0230   if (qe.solutions == 0) {
0231     return {{Intersection3D::invalid(), Intersection3D::invalid()}, this};
0232   }
0233
0234   // Check the validity of the first solution
0235   Vector3 solution1 = position + qe.first * direction;
0236   Intersection3D::Status status1 = std::abs(qe.first) < std::abs(tolerance)
0237                                        ? Intersection3D::Status::onSurface
0238                                        : Intersection3D::Status::reachable;
0239
0240   // Helper method for boundary check
0241   auto boundaryCheck =
0242       [&](const Vector3& solution,
0243           Intersection3D::Status status) -> Intersection3D::Status {
0244     // No check to be done, return current status
0245     if (!bcheck.isEnabled()) {
0246       return status;
0247     }
0248     const auto& cBounds = bounds();
0249     if (cBounds.coversFullAzimuth() &&
0250         bcheck.type() == BoundaryCheck::Type::eAbsolute) {
0251       // Project out the current Z value via local z axis
0252       // Built-in local to global for speed reasons
0253       const auto& tMatrix = gctxTransform.matrix();
0254       // Create the reference vector in local
0255       const Vector3 vecLocal(solution - tMatrix.block<3, 1>(0, 3));
0256       double cZ = vecLocal.dot(tMatrix.block<3, 1>(0, 2));
0257       double modifiedTolerance = tolerance + bcheck.tolerance()[eBoundLoc1];
0258       double hZ = cBounds.get(CylinderBounds::eHalfLengthZ) + modifiedTolerance;
0259       return std::abs(cZ) < std::abs(hZ) ? status
0260                                          : Intersection3D::Status::missed;
0261     }
0262     return (isOnSurface(gctx, solution, direction, bcheck)
0263                 ? status
0264                 : Intersection3D::Status::missed);
0265   };
0266   // Check first solution for boundary compatibility
0267   status1 = boundaryCheck(solution1, status1);
0268   // Set the intersection
0269   Intersection3D first(solution1, qe.first, status1);
0270   if (qe.solutions == 1) {
0271     return {{first, first}, this};
0272   }
0273   // Check the validity of the second solution
0274   Vector3 solution2 = position + qe.second * direction;
0275   Intersection3D::Status status2 = std::abs(qe.second) < std::abs(tolerance)
0276                                        ? Intersection3D::Status::onSurface
0277                                        : Intersection3D::Status::reachable;
0278   // Check first solution for boundary compatibility
0279   status2 = boundaryCheck(solution2, status2);
0280   Intersection3D second(solution2, qe.second, status2);
0281   // Order based on path length
0282   if (first.pathLength() <= second.pathLength()) {
0283     return {{first, second}, this};
0284   }
0285   return {{second, first}, this};
0286 }
0287
0288 Acts::AlignmentToPathMatrix Acts::CylinderSurface::alignmentToPathDerivative(
0289     const GeometryContext& gctx, const Vector3& position,
0290     const Vector3& direction) const {
0291   assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
0292
0293   // The vector between position and center
0294   const auto pcRowVec = (position - center(gctx)).transpose().eval();
0295   // The rotation
0296   const auto& rotation = transform(gctx).rotation();
0297   // The local frame x/y/z axis
0298   const auto& localXAxis = rotation.col(0);
0299   const auto& localYAxis = rotation.col(1);
0300   const auto& localZAxis = rotation.col(2);
0301   // The local coordinates
0302   const auto localPos = (rotation.transpose() * position).eval();
0303   const auto dx = direction.dot(localXAxis);
0304   const auto dy = direction.dot(localYAxis);
0305   const auto dz = direction.dot(localZAxis);
0306   // The normalization factor
0307   const auto norm = 1 / (1 - dz * dz);
0308   // The direction transpose
0309   const auto& dirRowVec = direction.transpose();
0310   // The derivative of path w.r.t. the local axes
0311   // @note The following calculations assume that the intersection of the track
0312   // with the cylinder always satisfy: perp(localPos) = R
0313   const auto localXAxisToPath =
0314       (-2 * norm * (dx * pcRowVec + localPos.x() * dirRowVec)).eval();
0315   const auto localYAxisToPath =
0316       (-2 * norm * (dy * pcRowVec + localPos.y() * dirRowVec)).eval();
0317   const auto localZAxisToPath =
0318       (-4 * norm * norm * (dx * localPos.x() + dy * localPos.y()) * dz *
0319        dirRowVec)
0320           .eval();
0321   // Calculate the derivative of local frame axes w.r.t its rotation
0322   const auto [rotToLocalXAxis, rotToLocalYAxis, rotToLocalZAxis] =
0323       detail::rotationToLocalAxesDerivative(rotation);
0324   // Initialize the derivative of propagation path w.r.t. local frame
0325   // translation (origin) and rotation
0326   AlignmentToPathMatrix alignToPath = AlignmentToPathMatrix::Zero();
0327   alignToPath.segment<3>(eAlignmentCenter0) =
0328       2 * norm * (dx * localXAxis.transpose() + dy * localYAxis.transpose());
0329   alignToPath.segment<3>(eAlignmentRotation0) =
0330       localXAxisToPath * rotToLocalXAxis + localYAxisToPath * rotToLocalYAxis +
0331       localZAxisToPath * rotToLocalZAxis;
0332
0333   return alignToPath;
0334 }
0335
0336 Acts::ActsMatrix<2, 3>
0337 Acts::CylinderSurface::localCartesianToBoundLocalDerivative(
0338     const GeometryContext& gctx, const Vector3& position) const {
0339   using VectorHelpers::perp;
0340   using VectorHelpers::phi;
0341   // The local frame transform
0342   const auto& sTransform = transform(gctx);
0343   // calculate the transformation to local coordinates
0344   const Vector3 localPos = sTransform.inverse() * position;
0345   const double lr = perp(localPos);
0346   const double lphi = phi(localPos);
0347   const double lcphi = std::cos(lphi);
0348   const double lsphi = std::sin(lphi);
0349   // Solve for radius R
0350   double R = bounds().get(CylinderBounds::eR);
0351   ActsMatrix<2, 3> loc3DToLocBound = ActsMatrix<2, 3>::Zero();
0352   loc3DToLocBound << -R * lsphi / lr, R * lcphi / lr, 0, 0, 0, 1;
0353
0354   return loc3DToLocBound;
0355 }