File indexing completed on 2025-07-09 07:49:57
0001
0002
0003
0004
0005
0006
0007
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 }
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
0140
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
0206
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->transform(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().transform(gctx)).translation()[eZ] <
0232 (itransform * b.surface().transform(gctx)).translation()[eZ];
0233 });
0234
0235 for (const auto& [i, child] : enumerate(trivialLinks)) {
0236 const auto& bounds =
0237 dynamic_cast<const CylinderBounds&>(child.surface().bounds());
0238 Transform3 ltransform = itransform * child.surface().transform(gctx);
0239 double hlZ = bounds.get(CylinderBounds::eHalfLengthZ);
0240 double minZ = ltransform.translation()[eZ] - hlZ;
0241 double maxZ = ltransform.translation()[eZ] + hlZ;
0242 if (i == 0) {
0243 edges.push_back(minZ);
0244 }
0245 edges.push_back(maxZ);
0246 }
0247
0248 ACTS_VERBOSE("~> Determined bin edges to be " << printEdges(edges));
0249
0250 Axis axis{AxisBound, edges};
0251
0252 auto gridPortalLink =
0253 GridPortalLink::make(m_surface, m_direction, std::move(axis));
0254 for (const auto& [i, child] : enumerate(trivialLinks)) {
0255 gridPortalLink->grid().atLocalBins({i + 1}) = &child.volume();
0256 }
0257
0258 return gridPortalLink;
0259
0260 } else if (surface().type() == Surface::SurfaceType::Disc) {
0261 ACTS_VERBOSE("Combining composite into disc grid");
0262
0263 if (m_direction != AxisDirection::AxisR) {
0264 ACTS_ERROR("Disc grid only supports binning in R direction");
0265 throw std::runtime_error{"Unsupported binning direction"};
0266 }
0267
0268 std::vector<double> edges;
0269 edges.reserve(m_children.size() + 1);
0270
0271 std::ranges::sort(trivialLinks, [](const auto& a, const auto& b) {
0272 const auto& boundsA =
0273 dynamic_cast<const RadialBounds&>(a.surface().bounds());
0274 const auto& boundsB =
0275 dynamic_cast<const RadialBounds&>(b.surface().bounds());
0276 return boundsA.get(RadialBounds::eMinR) <
0277 boundsB.get(RadialBounds::eMinR);
0278 });
0279
0280 for (const auto& [i, child] : enumerate(trivialLinks)) {
0281 const auto& bounds =
0282 dynamic_cast<const RadialBounds&>(child.surface().bounds());
0283
0284 if (i == 0) {
0285 edges.push_back(bounds.get(RadialBounds::eMinR));
0286 }
0287 edges.push_back(bounds.get(RadialBounds::eMaxR));
0288 }
0289
0290 ACTS_VERBOSE("~> Determined bin edges to be " << printEdges(edges));
0291
0292 Axis axis{AxisBound, edges};
0293
0294 auto grid = GridPortalLink::make(m_surface, m_direction, std::move(axis));
0295 for (const auto& [i, child] : enumerate(trivialLinks)) {
0296 grid->grid().atLocalBins({i + 1}) = &child.volume();
0297 }
0298
0299 return grid;
0300 } else if (surface().type() == Surface::SurfaceType::Plane) {
0301 ACTS_VERBOSE("Combining composite into plane grid");
0302
0303 if (m_direction != AxisDirection::AxisX &&
0304 m_direction != AxisDirection::AxisY) {
0305 ACTS_ERROR("Plane grid only supports binning in x/y direction");
0306 throw std::runtime_error{"Unsupported binning direction"};
0307 }
0308
0309 bool dirX = m_direction == AxisDirection::AxisX;
0310
0311 std::vector<double> edges;
0312 edges.reserve(m_children.size() + 1);
0313
0314 const Transform3& groupTransform = m_surface->transform(gctx);
0315 Transform3 itransform = groupTransform.inverse();
0316
0317 std::size_t sortingDir = dirX ? eX : eY;
0318 std::ranges::sort(trivialLinks, [&itransform, &gctx, sortingDir](
0319 const auto& a, const auto& b) {
0320 return (itransform * a.surface().transform(gctx))
0321 .translation()[sortingDir] <
0322 (itransform * b.surface().transform(gctx))
0323 .translation()[sortingDir];
0324 });
0325
0326 for (const auto& [i, child] : enumerate(trivialLinks)) {
0327 const auto& bounds =
0328 dynamic_cast<const RectangleBounds&>(child.surface().bounds());
0329 Transform3 ltransform = itransform * child.surface().transform(gctx);
0330 double half = dirX ? bounds.halfLengthX() : bounds.halfLengthY();
0331 double min = ltransform.translation()[sortingDir] - half;
0332 double max = ltransform.translation()[sortingDir] + half;
0333 if (i == 0) {
0334 edges.push_back(min);
0335 }
0336 edges.push_back(max);
0337 }
0338
0339 ACTS_VERBOSE("~> Determined bin edges to be " << printEdges(edges));
0340
0341 Axis axis{AxisBound, edges};
0342
0343 auto grid = GridPortalLink::make(m_surface, m_direction, std::move(axis));
0344 for (const auto& [i, child] : enumerate(trivialLinks)) {
0345 grid->grid().atLocalBins({i + 1}) = &child.volume();
0346 }
0347
0348 grid->setArtifactPortalLinks(std::move(trivialLinks));
0349
0350 return grid;
0351
0352 } else {
0353 throw std::invalid_argument{"Unsupported surface type"};
0354 }
0355 }
0356
0357 }