src/Navigation/CylinderNavigationPolicy.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/Navigation/CylinderNavigationPolicy.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/CylinderVolumeBounds.hpp"
0013 #include "Acts/Geometry/TrackingVolume.hpp"
0014 #include "Acts/Geometry/VolumeBounds.hpp"
0015 #include "Acts/Surfaces/CylinderBounds.hpp"
0016 #include "Acts/Surfaces/RadialBounds.hpp"
0017
0018 namespace Acts {
0019
0020 namespace {
0021 std::string_view faceName(CylinderVolumeBounds::Face face) {
0022   using enum CylinderVolumeBounds::Face;
0023   switch (face) {
0024     case PositiveDisc:
0025       return "PositiveDisc";
0026     case NegativeDisc:
0027       return "NegativeDisc";
0028     case OuterCylinder:
0029       return "OuterCylinder";
0030     case InnerCylinder:
0031       return "InnerCylinder";
0032     default:
0033       return "Unknown";
0034   }
0035 }
0036 }  // namespace
0037
0038 CylinderNavigationPolicy::CylinderNavigationPolicy(const GeometryContext& gctx,
0039                                                    const TrackingVolume& volume,
0040                                                    const Logger& logger)
0041     : m_volume(&volume) {
0042   ACTS_VERBOSE("CylinderNavigationPolicy constructor for volume "
0043                << volume.volumeName());
0044   using enum CylinderVolumeBounds::Face;
0045
0046   if (m_volume->volumeBounds().type() != VolumeBounds::eCylinder) {
0047     ACTS_ERROR("CylinderNavigationPolicy can only be used with "
0048                << "cylinder volumes");
0049     throw std::invalid_argument(
0050         "CylinderNavigationPolicy can only be used with "
0051         "cylinder volumes");
0052   }
0053
0054   auto& bounds =
0055       dynamic_cast<const CylinderVolumeBounds&>(m_volume->volumeBounds());
0056
0057   double rMin = bounds.get(CylinderVolumeBounds::eMinR);
0058   m_rMin2 = rMin * rMin;
0059   double rMax = bounds.get(CylinderVolumeBounds::eMaxR);
0060   m_rMax2 = rMax * rMax;
0061   m_halfLengthZ = bounds.get(CylinderVolumeBounds::eHalfLengthZ);
0062
0063   if (rMin == 0) {
0064     ACTS_ERROR("CylinderNavigationPolicy can only be used with "
0065                << "non-zero inner radius, which " << m_volume->volumeName()
0066                << " has not");
0067     throw std::invalid_argument(
0068         "CylinderNavigationPolicy can only be used with "
0069         "non-zero inner radius");
0070   }
0071
0072   if (m_volume->portals().size() != 4) {
0073     ACTS_ERROR("CylinderNavigationPolicy can only be used with "
0074                << "volumes with 4 portals");
0075     throw std::invalid_argument(
0076         "CylinderNavigationPolicy can only be used with "
0077         "volumes with 4 portals");
0078   }
0079
0080   ACTS_VERBOSE("CylinderNavigationPolicy created for volume "
0081                << volume.volumeName());
0082
0083   if (!volume.localToGlobalTransform(gctx).linear().isApprox(
0084           SquareMatrix3::Identity())) {
0085     m_itransform = volume.globalToLocalTransform(gctx);
0086   }
0087
0088   // Since the volume does not store the shell assignment, we have to recover
0089   // this from the raw portals
0090   for (const auto& portal : m_volume->portals()) {
0091     if (const auto* cylBounds =
0092             dynamic_cast<const CylinderBounds*>(&portal.surface().bounds());
0093         cylBounds != nullptr) {
0094       if (std::abs(cylBounds->get(CylinderBounds::eR) - rMin) <
0095           s_onSurfaceTolerance) {
0096         m_portals.at(toUnderlying(InnerCylinder)) = &portal;
0097         continue;
0098       }
0099
0100       if (std::abs(cylBounds->get(CylinderBounds::eR) - rMax) <
0101           s_onSurfaceTolerance) {
0102         m_portals.at(toUnderlying(OuterCylinder)) = &portal;
0103         continue;
0104       }
0105     }
0106
0107     if (const auto* radBounds =
0108             dynamic_cast<const RadialBounds*>(&portal.surface().bounds());
0109         radBounds != nullptr) {
0110       Transform3 localTransform = m_volume->globalToLocalTransform(gctx) *
0111                                   portal.surface().localToGlobalTransform(gctx);
0112       Vector3 localPosition = localTransform.translation();
0113       double localZ = localPosition.z();
0114       if (std::abs(localZ - m_halfLengthZ) < s_onSurfaceTolerance) {
0115         m_portals.at(toUnderlying(PositiveDisc)) = &portal;
0116         continue;
0117       }
0118
0119       if (std::abs(localZ + m_halfLengthZ) < s_onSurfaceTolerance) {
0120         m_portals.at(toUnderlying(NegativeDisc)) = &portal;
0121         continue;
0122       }
0123     }
0124   }
0125
0126   ACTS_VERBOSE("Portal assignment:");
0127   for (const auto& [i, portal] : enumerate(m_portals)) {
0128     auto face = static_cast<CylinderVolumeBounds::Face>(i);
0129
0130     if (portal == nullptr) {
0131       ACTS_ERROR("Have no portal for " << faceName(face));
0132       throw std::invalid_argument("Have no portal for " +
0133                                   std::string{faceName(face)});
0134     }
0135
0136     ACTS_VERBOSE("  " << faceName(face) << " -> "
0137                       << portal->surface().toStream(gctx));
0138   }
0139 }
0140
0141 void CylinderNavigationPolicy::initializeCandidates(
0142     [[maybe_unused]] const GeometryContext& gctx,
0143     const NavigationArguments& args, AppendOnlyNavigationStream& stream,
0144     const Logger& logger) const {
0145   using enum CylinderVolumeBounds::Face;
0146
0147   ACTS_VERBOSE("CylinderNavigationPolicy::initializeCandidates for volume "
0148                << m_volume->volumeName()
0149                << " with gpos: " << args.position.transpose()
0150                << " gdir: " << args.direction.transpose());
0151
0152   Vector3 pos;
0153   Vector3 dir;
0154   if (m_itransform.has_value()) {
0155     pos = *m_itransform * args.position;
0156     dir = m_itransform->linear() * args.direction;
0157   } else {
0158     pos = args.position - m_volume->center(gctx);
0159     dir = args.direction;
0160   }
0161
0162   ACTS_VERBOSE("-> lpos: " << pos.transpose() << " ldir: " << dir.transpose());
0163
0164   auto add = [&](auto face) {
0165     ACTS_VERBOSE("~~> Adding portal candidate " << faceName(face));
0166     stream.addPortalCandidate(*m_portals.at(toUnderlying(face)));
0167   };
0168
0169   bool hitDisk = false;
0170   Vector3 diskIntersection;
0171   double diskDistance3D = std::numeric_limits<double>::max();
0172   double zDisk = (dir[2] > 0) ? m_halfLengthZ : -m_halfLengthZ;
0173   if (std::abs(dir[2]) > s_onSurfaceTolerance) {
0174     ACTS_VERBOSE(
0175         "-> Not parallel to the disc, see if we're inside the disc radii");
0176
0177     double t = (zDisk - pos[2]) / dir[2];
0178
0179     diskIntersection = pos + t * dir;
0180     double r2 = diskIntersection[0] * diskIntersection[0] +
0181                 diskIntersection[1] * diskIntersection[1];
0182     if (r2 < m_rMax2 && r2 > m_rMin2 &&
0183         t > 0) {  // Only consider forward intersections
0184       hitDisk = true;
0185       diskDistance3D = t;  // Parameter t is the distance along the ray
0186     }
0187   } else {
0188     ACTS_VERBOSE("-> Parallel to the disc, see if we're inside the disc radii");
0189   }
0190
0191   double rpos2 = pos[0] * pos[0] + pos[1] * pos[1];
0192   ACTS_VERBOSE("-> rpos: " << std::sqrt(rpos2));
0193
0194   bool hitInner = false;
0195   if (std::abs(rpos2 - m_rMin2) > s_onSurfaceTolerance) {
0196     // Find point of closest approach to the origin
0197
0198     Vector2 dir2 = dir.head<2>().normalized();
0199     double d = -1 * pos.head<2>().dot(dir2);
0200     ACTS_VERBOSE("-> d: " << d);
0201
0202     if (d > 0) {  // Point of closest approach is in the direction of the ray
0203
0204       double clippedD = d;
0205       if (hitDisk) {
0206         // Clip to distance of disk intersection
0207         double diskIntersectionDistanceXy =
0208             (diskIntersection.head<2>() - pos.head<2>()).norm();
0209         clippedD = std::min(d, diskIntersectionDistanceXy);
0210       }
0211
0212       Vector2 poc = pos.head<2>() + clippedD * dir2;
0213       double r2 = poc.dot(poc);
0214       hitInner = r2 < m_rMin2;
0215       if (hitInner) {
0216         add(InnerCylinder);
0217         // If we hit the inner cylinder before reaching the disk intersection,
0218         // we can discard the disk as we'll never reach it
0219         if (hitDisk) {
0220           // Calculate 3D distance to inner cylinder intersection
0221           double innerDistance3D =
0222               d / dir.head<2>().norm();  // Convert 2D distance to 3D parameter
0223           if (innerDistance3D < diskDistance3D) {
0224             hitDisk = false;
0225           }
0226         }
0227       }
0228     }
0229   }
0230
0231   // Add disk candidate if we determined it's reachable (not blocked by inner
0232   // cylinder)
0233   if (hitDisk) {
0234     add(dir[2] > 0 ? PositiveDisc : NegativeDisc);
0235   }
0236
0237   if (!hitInner && !hitDisk) {
0238     add(OuterCylinder);
0239   }
0240 }
0241
0242 void CylinderNavigationPolicy::connect(NavigationDelegate& delegate) const {
0243   connectDefault<CylinderNavigationPolicy>(delegate);
0244 }
0245 }  // namespace Acts