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     const NavigationArguments& args, AppendOnlyNavigationStream& stream,
0142     const Logger& logger) const {
0143   using enum CylinderVolumeBounds::Face;
0144
0145   if (!args.wantsPortals) {
0146     return;
0147   }
0148
0149   ACTS_VERBOSE("CylinderNavigationPolicy::initializeCandidates for volume "
0150                << m_volume->volumeName()
0151                << " with gpos: " << args.position.transpose()
0152                << " gdir: " << args.direction.transpose());
0153
0154   Vector3 pos;
0155   Vector3 dir;
0156   if (m_itransform.has_value()) {
0157     pos = *m_itransform * args.position;
0158     dir = m_itransform->linear() * args.direction;
0159   } else {
0160     pos = args.position - m_volume->transform().translation();
0161     dir = args.direction;
0162   }
0163
0164   ACTS_VERBOSE("-> lpos: " << pos.transpose() << " ldir: " << dir.transpose());
0165
0166   auto add = [&](auto face) {
0167     ACTS_VERBOSE("~~> Adding portal candidate " << faceName(face));
0168     stream.addPortalCandidate(*m_portals.at(toUnderlying(face)));
0169   };
0170
0171   bool hitDisk = false;
0172   Vector3 diskIntersection;
0173   double diskDistance3D = std::numeric_limits<double>::max();
0174   double zDisk = (dir[2] > 0) ? m_halfLengthZ : -m_halfLengthZ;
0175   if (std::abs(dir[2]) > s_onSurfaceTolerance) {
0176     ACTS_VERBOSE(
0177         "-> Not parallel to the disc, see if we're inside the disc radii");
0178
0179     double t = (zDisk - pos[2]) / dir[2];
0180
0181     diskIntersection = pos + t * dir;
0182     double r2 = diskIntersection[0] * diskIntersection[0] +
0183                 diskIntersection[1] * diskIntersection[1];
0184     if (r2 < m_rMax2 && r2 > m_rMin2 &&
0185         t > 0) {  // Only consider forward intersections
0186       hitDisk = true;
0187       diskDistance3D = t;  // Parameter t is the distance along the ray
0188     }
0189   } else {
0190     ACTS_VERBOSE("-> Parallel to the disc, see if we're inside the disc radii");
0191   }
0192
0193   double rpos2 = pos[0] * pos[0] + pos[1] * pos[1];
0194   ACTS_VERBOSE("-> rpos: " << std::sqrt(rpos2));
0195
0196   bool hitInner = false;
0197   if (std::abs(rpos2 - m_rMin2) > s_onSurfaceTolerance) {
0198     // Find point of closest approach to the origin
0199
0200     Vector2 dir2 = dir.head<2>().normalized();
0201     double d = -1 * pos.head<2>().dot(dir2);
0202     ACTS_VERBOSE("-> d: " << d);
0203
0204     if (d > 0) {  // Point of closest approach is in the direction of the ray
0205
0206       double clippedD = d;
0207       if (hitDisk) {
0208         // Clip to distance of disk intersection
0209         double diskIntersectionDistanceXy =
0210             (diskIntersection.head<2>() - pos.head<2>()).norm();
0211         clippedD = std::min(d, diskIntersectionDistanceXy);
0212       }
0213
0214       Vector2 poc = pos.head<2>() + clippedD * dir2;
0215       double r2 = poc.dot(poc);
0216       hitInner = r2 < m_rMin2;
0217       if (hitInner) {
0218         add(InnerCylinder);
0219         // If we hit the inner cylinder before reaching the disk intersection,
0220         // we can discard the disk as we'll never reach it
0221         if (hitDisk) {
0222           // Calculate 3D distance to inner cylinder intersection
0223           double innerDistance3D =
0224               d / dir.head<2>().norm();  // Convert 2D distance to 3D parameter
0225           if (innerDistance3D < diskDistance3D) {
0226             hitDisk = false;
0227           }
0228         }
0229       }
0230     }
0231   }
0232
0233   // Add disk candidate if we determined it's reachable (not blocked by inner
0234   // cylinder)
0235   if (hitDisk) {
0236     add(dir[2] > 0 ? PositiveDisc : NegativeDisc);
0237   }
0238
0239   if (!hitInner && !hitDisk) {
0240     add(OuterCylinder);
0241   }
0242 }
0243
0244 void CylinderNavigationPolicy::connect(NavigationDelegate& delegate) const {
0245   connectDefault<CylinderNavigationPolicy>(delegate);
0246 }
0247 }  // namespace Acts