Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-09 07:49:57

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->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 }  // namespace Acts