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.transform().linear().isApprox(SquareMatrix3::Identity())) {
0084     m_itransform = volume.transform().inverse();
0085   }
0086
0087   // Since the volume does not store the shell assignment, we have to recover
0088   // this from the raw portals
0089   for (const auto& portal : m_volume->portals()) {
0090     if (const auto* cylBounds =
0091             dynamic_cast<const CylinderBounds*>(&portal.surface().bounds());
0092         cylBounds != nullptr) {
0093       if (std::abs(cylBounds->get(CylinderBounds::eR) - rMin) <
0094           s_onSurfaceTolerance) {
0095         m_portals.at(toUnderlying(InnerCylinder)) = &portal;
0096         continue;
0097       }
0098
0099       if (std::abs(cylBounds->get(CylinderBounds::eR) - rMax) <
0100           s_onSurfaceTolerance) {
0101         m_portals.at(toUnderlying(OuterCylinder)) = &portal;
0102         continue;
0103       }
0104     }
0105
0106     if (const auto* radBounds =
0107             dynamic_cast<const RadialBounds*>(&portal.surface().bounds());
0108         radBounds != nullptr) {
0109       Transform3 localTransform =
0110           m_volume->transform().inverse() * portal.surface().transform(gctx);
0111       Vector3 localPosition = localTransform.translation();
0112       double localZ = localPosition.z();
0113       if (std::abs(localZ - m_halfLengthZ) < s_onSurfaceTolerance) {
0114         m_portals.at(toUnderlying(PositiveDisc)) = &portal;
0115         continue;
0116       }
0117
0118       if (std::abs(localZ + m_halfLengthZ) < s_onSurfaceTolerance) {
0119         m_portals.at(toUnderlying(NegativeDisc)) = &portal;
0120         continue;
0121       }
0122     }
0123   }
0124
0125   ACTS_VERBOSE("Portal assignment:");
0126   for (const auto& [i, portal] : enumerate(m_portals)) {
0127     auto face = static_cast<CylinderVolumeBounds::Face>(i);
0128
0129     if (portal == nullptr) {
0130       ACTS_ERROR("Have no portal for " << faceName(face));
0131       throw std::invalid_argument("Have no portal for " +
0132                                   std::string{faceName(face)});
0133     }
0134
0135     ACTS_VERBOSE("  " << faceName(face) << " -> "
0136                       << portal->surface().toStream(gctx));
0137   }
0138 }
0139
0140 void CylinderNavigationPolicy::initializeCandidates(
0141     [[maybe_unused]] const GeometryContext& gctx,
0142     const NavigationArguments& args, AppendOnlyNavigationStream& stream,
0143     const Logger& logger) const {
0144   using enum CylinderVolumeBounds::Face;
0145
0146   ACTS_VERBOSE("CylinderNavigationPolicy::initializeCandidates for volume "
0147                << m_volume->volumeName()
0148                << " with gpos: " << args.position.transpose()
0149                << " gdir: " << args.direction.transpose());
0150
0151   Vector3 pos;
0152   Vector3 dir;
0153   if (m_itransform.has_value()) {
0154     pos = *m_itransform * args.position;
0155     dir = m_itransform->linear() * args.direction;
0156   } else {
0157     pos = args.position - m_volume->transform().translation();
0158     dir = args.direction;
0159   }
0160
0161   ACTS_VERBOSE("-> lpos: " << pos.transpose() << " ldir: " << dir.transpose());
0162
0163   auto add = [&](auto face) {
0164     ACTS_VERBOSE("~~> Adding portal candidate " << faceName(face));
0165     stream.addPortalCandidate(*m_portals.at(toUnderlying(face)));
0166   };
0167
0168   bool hitDisk = false;
0169   Vector3 diskIntersection;
0170   double diskDistance3D = std::numeric_limits<double>::max();
0171   double zDisk = (dir[2] > 0) ? m_halfLengthZ : -m_halfLengthZ;
0172   if (std::abs(dir[2]) > s_onSurfaceTolerance) {
0173     ACTS_VERBOSE(
0174         "-> Not parallel to the disc, see if we're inside the disc radii");
0175
0176     double t = (zDisk - pos[2]) / dir[2];
0177
0178     diskIntersection = pos + t * dir;
0179     double r2 = diskIntersection[0] * diskIntersection[0] +
0180                 diskIntersection[1] * diskIntersection[1];
0181     if (r2 < m_rMax2 && r2 > m_rMin2 &&
0182         t > 0) {  // Only consider forward intersections
0183       hitDisk = true;
0184       diskDistance3D = t;  // Parameter t is the distance along the ray
0185     }
0186   } else {
0187     ACTS_VERBOSE("-> Parallel to the disc, see if we're inside the disc radii");
0188   }
0189
0190   double rpos2 = pos[0] * pos[0] + pos[1] * pos[1];
0191   ACTS_VERBOSE("-> rpos: " << std::sqrt(rpos2));
0192
0193   bool hitInner = false;
0194   if (std::abs(rpos2 - m_rMin2) > s_onSurfaceTolerance) {
0195     // Find point of closest approach to the origin
0196
0197     Vector2 dir2 = dir.head<2>().normalized();
0198     double d = -1 * pos.head<2>().dot(dir2);
0199     ACTS_VERBOSE("-> d: " << d);
0200
0201     if (d > 0) {  // Point of closest approach is in the direction of the ray
0202
0203       double clippedD = d;
0204       if (hitDisk) {
0205         // Clip to distance of disk intersection
0206         double diskIntersectionDistanceXy =
0207             (diskIntersection.head<2>() - pos.head<2>()).norm();
0208         clippedD = std::min(d, diskIntersectionDistanceXy);
0209       }
0210
0211       Vector2 poc = pos.head<2>() + clippedD * dir2;
0212       double r2 = poc.dot(poc);
0213       hitInner = r2 < m_rMin2;
0214       if (hitInner) {
0215         add(InnerCylinder);
0216         // If we hit the inner cylinder before reaching the disk intersection,
0217         // we can discard the disk as we'll never reach it
0218         if (hitDisk) {
0219           // Calculate 3D distance to inner cylinder intersection
0220           double innerDistance3D =
0221               d / dir.head<2>().norm();  // Convert 2D distance to 3D parameter
0222           if (innerDistance3D < diskDistance3D) {
0223             hitDisk = false;
0224           }
0225         }
0226       }
0227     }
0228   }
0229
0230   // Add disk candidate if we determined it's reachable (not blocked by inner
0231   // cylinder)
0232   if (hitDisk) {
0233     add(dir[2] > 0 ? PositiveDisc : NegativeDisc);
0234   }
0235
0236   if (!hitInner && !hitDisk) {
0237     add(OuterCylinder);
0238   }
0239 }
0240
0241 void CylinderNavigationPolicy::connect(NavigationDelegate& delegate) const {
0242   connectDefault<CylinderNavigationPolicy>(delegate);
0243 }
0244 }  // namespace Acts