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