src/Geometry/CompositePortalLink.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/Geometry/CompositePortalLink.hpp"
0010
0011 #include "Acts/Geometry/GridPortalLink.hpp"
0012 #include "Acts/Geometry/PortalError.hpp"
0013 #include "Acts/Geometry/TrivialPortalLink.hpp"
0014 #include "Acts/Surfaces/CylinderSurface.hpp"
0015 #include "Acts/Surfaces/DiscSurface.hpp"
0016 #include "Acts/Surfaces/PlaneSurface.hpp"
0017 #include "Acts/Surfaces/RadialBounds.hpp"
0018 #include "Acts/Surfaces/RegularSurface.hpp"
0019 #include "Acts/Utilities/Axis.hpp"
0020
0021 #include <algorithm>
0022 #include <iostream>
0023 #include <iterator>
0024 #include <stdexcept>
0025
0026 #include <boost/algorithm/string/join.hpp>
0027
0028 namespace Acts {
0029
0030 namespace {
0031 std::shared_ptr<RegularSurface> mergedSurface(const Surface& a,
0032                                               const Surface& b,
0033                                               AxisDirection direction) {
0034   if (a.type() != b.type()) {
0035     throw std::invalid_argument{"Cannot merge surfaces of different types"};
0036   }
0037
0038   if (const auto* cylinderA = dynamic_cast<const CylinderSurface*>(&a);
0039       cylinderA != nullptr) {
0040     const auto& cylinderB = dynamic_cast<const CylinderSurface&>(b);
0041
0042     auto [merged, reversed] = cylinderA->mergedWith(cylinderB, direction, true);
0043     return merged;
0044   } else if (const auto* discA = dynamic_cast<const DiscSurface*>(&a);
0045              discA != nullptr) {
0046     const auto& discB = dynamic_cast<const DiscSurface&>(b);
0047     auto [merged, reversed] = discA->mergedWith(discB, direction, true);
0048     return merged;
0049   } else if (const auto* planeA = dynamic_cast<const PlaneSurface*>(&a);
0050              planeA != nullptr) {
0051     throw std::logic_error{"Plane surfaces not implemented yet"};
0052   } else {
0053     throw std::invalid_argument{"Unsupported surface type"};
0054   }
0055 }
0056
0057 std::shared_ptr<RegularSurface> mergePortalLinks(
0058     const std::vector<std::unique_ptr<PortalLinkBase>>& links,
0059     AxisDirection direction) {
0060   assert(std::ranges::all_of(links,
0061                              [](const auto& link) { return link != nullptr; }));
0062   assert(!links.empty());
0063
0064   std::shared_ptr<RegularSurface> result = links.front()->surfacePtr();
0065   for (auto it = std::next(links.begin()); it != links.end(); ++it) {
0066     assert(result != nullptr);
0067     result = mergedSurface(*result, it->get()->surface(), direction);
0068   }
0069
0070   return result;
0071 }
0072 }  // namespace
0073
0074 CompositePortalLink::CompositePortalLink(std::unique_ptr<PortalLinkBase> a,
0075                                          std::unique_ptr<PortalLinkBase> b,
0076                                          AxisDirection direction, bool flatten)
0077     : PortalLinkBase(mergedSurface(a->surface(), b->surface(), direction)),
0078       m_direction{direction} {
0079   if (!flatten) {
0080     m_children.push_back(std::move(a));
0081     m_children.push_back(std::move(b));
0082
0083   } else {
0084     auto handle = [&](std::unique_ptr<PortalLinkBase> link) {
0085       if (auto* composite = dynamic_cast<CompositePortalLink*>(link.get());
0086           composite != nullptr) {
0087         m_children.insert(
0088             m_children.end(),
0089             std::make_move_iterator(composite->m_children.begin()),
0090             std::make_move_iterator(composite->m_children.end()));
0091       } else {
0092         m_children.push_back(std::move(link));
0093       }
0094     };
0095
0096     handle(std::move(a));
0097     handle(std::move(b));
0098   }
0099 }
0100
0101 CompositePortalLink::CompositePortalLink(
0102     std::vector<std::unique_ptr<PortalLinkBase>> links, AxisDirection direction,
0103     bool flatten)
0104     : PortalLinkBase(mergePortalLinks(links, direction)),
0105       m_direction(direction) {
0106   if (!flatten) {
0107     for (auto& child : links) {
0108       m_children.push_back(std::move(child));
0109     }
0110
0111   } else {
0112     auto handle = [&](std::unique_ptr<PortalLinkBase> link) {
0113       if (auto* composite = dynamic_cast<CompositePortalLink*>(link.get());
0114           composite != nullptr) {
0115         m_children.insert(
0116             m_children.end(),
0117             std::make_move_iterator(composite->m_children.begin()),
0118             std::make_move_iterator(composite->m_children.end()));
0119       } else {
0120         m_children.push_back(std::move(link));
0121       }
0122     };
0123
0124     for (auto& child : links) {
0125       handle(std::move(child));
0126     }
0127   }
0128 }
0129
0130 Result<const TrackingVolume*> CompositePortalLink::resolveVolume(
0131     const GeometryContext& gctx, const Vector2& position,
0132     double tolerance) const {
0133   // In this overload, we have to go back to global, because the children have
0134   // their own local coordinate systems
0135   Vector3 global = m_surface->localToGlobal(gctx, position);
0136   auto res = resolveVolume(gctx, global, tolerance);
0137   if (!res.ok()) {
0138     return res.error();
0139   }
0140   return *res;
0141 }
0142
0143 Result<const TrackingVolume*> CompositePortalLink::resolveVolume(
0144     const GeometryContext& gctx, const Vector3& position,
0145     double tolerance) const {
0146   assert(m_surface->isOnSurface(gctx, position, BoundaryTolerance::None(),
0147                                 tolerance));
0148
0149   for (const auto& child : m_children) {
0150     if (child->surface().isOnSurface(gctx, position, BoundaryTolerance::None(),
0151                                      tolerance)) {
0152       return child->resolveVolume(gctx, position);
0153     }
0154   }
0155
0156   return PortalError::PositionNotOnAnyChildPortalLink;
0157 }
0158
0159 std::size_t CompositePortalLink::depth() const {
0160   std::size_t result = 1;
0161
0162   for (const auto& child : m_children) {
0163     if (const auto* composite =
0164             dynamic_cast<const CompositePortalLink*>(child.get())) {
0165       result = std::max(result, composite->depth() + 1);
0166     }
0167   }
0168
0169   return result;
0170 }
0171
0172 std::size_t CompositePortalLink::size() const {
0173   return m_children.size();
0174 }
0175
0176 void CompositePortalLink::toStream(std::ostream& os) const {
0177   os << "CompositePortalLink";
0178 }
0179
0180 std::unique_ptr<GridPortalLink> CompositePortalLink::makeGrid(
0181     const GeometryContext& gctx, const Logger& logger) const {
0182   ACTS_VERBOSE("Attempting to make a grid from a composite portal link");
0183   std::vector<TrivialPortalLink*> trivialLinks;
0184   std::ranges::transform(m_children, std::back_inserter(trivialLinks),
0185                          [](const auto& child) {
0186                            return dynamic_cast<TrivialPortalLink*>(child.get());
0187                          });
0188
0189   if (std::ranges::any_of(trivialLinks,
0190                           [](auto link) { return link == nullptr; })) {
0191     ACTS_VERBOSE(
0192         "Failed to make a grid from a composite portal link -> returning "
0193         "nullptr");
0194     return nullptr;
0195   }
0196
0197   // we already know all children have surfaces that are compatible, we produced
0198   // a merged surface for the overall dimensions.
0199
0200   auto printEdges = [](const auto& edges) -> std::string {
0201     std::vector<std::string> strings;
0202     std::ranges::transform(edges, std::back_inserter(strings),
0203                            [](const auto& v) { return std::to_string(v); });
0204     return boost::algorithm::join(strings, ", ");
0205   };
0206
0207   if (surface().type() == Surface::SurfaceType::Cylinder) {
0208     ACTS_VERBOSE("Combining composite into cylinder grid");
0209
0210     if (m_direction != AxisDirection::AxisZ) {
0211       ACTS_ERROR("Cylinder grid only supports binning in Z direction");
0212       throw std::runtime_error{"Unsupported binning direction"};
0213     }
0214
0215     std::vector<double> edges;
0216     edges.reserve(m_children.size() + 1);
0217
0218     const Transform3& groupTransform = m_surface->transform(gctx);
0219     Transform3 itransform = groupTransform.inverse();
0220
0221     std::ranges::sort(
0222         trivialLinks, [&itransform, &gctx](const auto& a, const auto& b) {
0223           return (itransform * a->surface().transform(gctx)).translation()[eZ] <
0224                  (itransform * b->surface().transform(gctx)).translation()[eZ];
0225         });
0226
0227     for (const auto& [i, child] : enumerate(trivialLinks)) {
0228       const auto& bounds =
0229           dynamic_cast<const CylinderBounds&>(child->surface().bounds());
0230       Transform3 ltransform = itransform * child->surface().transform(gctx);
0231       double hlZ = bounds.get(CylinderBounds::eHalfLengthZ);
0232       double minZ = ltransform.translation()[eZ] - hlZ;
0233       double maxZ = ltransform.translation()[eZ] + hlZ;
0234       if (i == 0) {
0235         edges.push_back(minZ);
0236       }
0237       edges.push_back(maxZ);
0238     }
0239
0240     ACTS_VERBOSE("~> Determined bin edges to be " << printEdges(edges));
0241
0242     Axis axis{AxisBound, edges};
0243
0244     auto grid = GridPortalLink::make(m_surface, m_direction, std::move(axis));
0245     for (const auto& [i, child] : enumerate(trivialLinks)) {
0246       grid->atLocalBins({i + 1}) = &child->volume();
0247     }
0248
0249     return grid;
0250
0251   } else if (surface().type() == Surface::SurfaceType::Disc) {
0252     ACTS_VERBOSE("Combining composite into disc grid");
0253
0254     if (m_direction != AxisDirection::AxisR) {
0255       ACTS_ERROR("Disc grid only supports binning in R direction");
0256       throw std::runtime_error{"Unsupported binning direction"};
0257     }
0258
0259     std::vector<double> edges;
0260     edges.reserve(m_children.size() + 1);
0261
0262     std::ranges::sort(trivialLinks, [](const auto& a, const auto& b) {
0263       const auto& boundsA =
0264           dynamic_cast<const RadialBounds&>(a->surface().bounds());
0265       const auto& boundsB =
0266           dynamic_cast<const RadialBounds&>(b->surface().bounds());
0267       return boundsA.get(RadialBounds::eMinR) <
0268              boundsB.get(RadialBounds::eMinR);
0269     });
0270
0271     for (const auto& [i, child] : enumerate(trivialLinks)) {
0272       const auto& bounds =
0273           dynamic_cast<const RadialBounds&>(child->surface().bounds());
0274
0275       if (i == 0) {
0276         edges.push_back(bounds.get(RadialBounds::eMinR));
0277       }
0278       edges.push_back(bounds.get(RadialBounds::eMaxR));
0279     }
0280
0281     ACTS_VERBOSE("~> Determined bin edges to be " << printEdges(edges));
0282
0283     Axis axis{AxisBound, edges};
0284
0285     auto grid = GridPortalLink::make(m_surface, m_direction, std::move(axis));
0286     for (const auto& [i, child] : enumerate(trivialLinks)) {
0287       grid->atLocalBins({i + 1}) = &child->volume();
0288     }
0289
0290     return grid;
0291   } else if (surface().type() == Surface::SurfaceType::Plane) {
0292     throw std::runtime_error{"Plane surfaces not implemented yet"};
0293   } else {
0294     throw std::invalid_argument{"Unsupported surface type"};
0295   }
0296 }
0297
0298 }  // namespace Acts