Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2026-07-16 08:07:49 

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/GridPortalLink.hpp"
0010 
0011 #include "Acts/Geometry/TrivialPortalLink.hpp"
0012 #include "Acts/Surfaces/PlaneSurface.hpp"
0013 #include "Acts/Surfaces/RadialBounds.hpp"
0014 #include "Acts/Surfaces/RectangleBounds.hpp"
0015 #include "Acts/Utilities/AnyGridView.hpp"
0016 #include "Acts/Utilities/AxisDefinitions.hpp"
0017 #include "Acts/Utilities/Grid.hpp"
0018 
0019 #include <iostream>
0020 #include <numbers>
0021 
0022 namespace Acts {
0023 
0024 GridPortalLink::GridPortalLink(std::shared_ptr<RegularSurface> surface,
0025                                AxisDirection direction)
0026     : PortalLinkBase(std::move(surface)), m_direction(direction) {}
0027 
0028 GridPortalLink::~GridPortalLink() = default;
0029 
0030 std::unique_ptr<GridPortalLink> GridPortalLink::make(
0031     const std::shared_ptr<RegularSurface>& surface, TrackingVolume& volume,
0032     AxisDirection direction) {
0033   std::unique_ptr<GridPortalLink> grid;
0034 
0035   if (const auto* cylinder =
0036           dynamic_cast<const CylinderSurface*>(surface.get());
0037       cylinder != nullptr) {
0038     if (direction == AxisDirection::AxisRPhi) {
0039       double r = cylinder->bounds().get(CylinderBounds::eR);
0040       if (cylinder->bounds().coversFullAzimuth()) {
0041         grid = GridPortalLink::make(
0042             surface, direction,
0043             Axis{AxisClosed, -std::numbers::pi * r, std::numbers::pi * r, 1});
0044       } else {
0045         double hlPhi = cylinder->bounds().get(CylinderBounds::eHalfPhiSector);
0046 
0047         grid = GridPortalLink::make(surface, direction,
0048                                     Axis{AxisBound, -hlPhi * r, hlPhi * r, 1});
0049       }
0050     } else if (direction == AxisDirection::AxisZ) {
0051       double hlZ = cylinder->bounds().get(CylinderBounds::eHalfLengthZ);
0052       grid = GridPortalLink::make(surface, direction,
0053                                   Axis{AxisBound, -hlZ, hlZ, 1});
0054     } else {
0055       throw std::invalid_argument{"Invalid binning direction"};
0056     }
0057   } else if (const auto* disc = dynamic_cast<const DiscSurface*>(surface.get());
0058              disc != nullptr) {
0059     const auto& bounds = dynamic_cast<const RadialBounds&>(disc->bounds());
0060     if (direction == AxisDirection::AxisR) {
0061       double minR = bounds.get(RadialBounds::eMinR);
0062       double maxR = bounds.get(RadialBounds::eMaxR);
0063       grid = GridPortalLink::make(surface, direction,
0064                                   Axis{AxisBound, minR, maxR, 1});
0065     } else if (direction == AxisDirection::AxisPhi) {
0066       if (bounds.coversFullAzimuth()) {
0067         grid = GridPortalLink::make(
0068             surface, direction,
0069             Axis{AxisClosed, -std::numbers::pi, std::numbers::pi, 1});
0070       } else {
0071         double hlPhi = bounds.get(RadialBounds::eHalfPhiSector);
0072         grid = GridPortalLink::make(surface, direction,
0073                                     Axis{AxisBound, -hlPhi, hlPhi, 1});
0074       }
0075     } else {
0076       throw std::invalid_argument{"Invalid binning direction"};
0077     }
0078   } else if (const auto* plane =
0079                  dynamic_cast<const PlaneSurface*>(surface.get());
0080              plane != nullptr) {
0081     const auto& bounds = dynamic_cast<const RectangleBounds&>(plane->bounds());
0082     if (direction != AxisDirection::AxisX &&
0083         direction != AxisDirection::AxisY) {
0084       throw std::invalid_argument{"Invalid binning direction"};
0085     }
0086     double min = (direction == AxisDirection::AxisX)
0087                      ? bounds.get(RectangleBounds::eMinX)
0088                      : bounds.get(RectangleBounds::eMinY);
0089     double max = (direction == AxisDirection::AxisX)
0090                      ? bounds.get(RectangleBounds::eMaxX)
0091                      : bounds.get(RectangleBounds::eMaxY);
0092     grid =
0093         GridPortalLink::make(surface, direction, Axis{AxisBound, min, max, 1});
0094   } else {
0095     throw std::invalid_argument{"Surface type is not supported"};
0096   }
0097 
0098   assert(grid != nullptr);
0099   grid->setVolume(&volume);
0100 
0101   return grid;
0102 }
0103 
0104 void GridPortalLink::checkConsistency(const IGrid& grid,
0105                                       AxisDirection direction,
0106                                       const CylinderSurface& cyl) {
0107   if (cyl.bounds().get(CylinderBounds::eAveragePhi) != 0) {
0108     throw std::invalid_argument(
0109         "GridPortalLink: CylinderBounds: only average phi == 0 is "
0110         "supported. Rotate the cylinder surface.");
0111   }
0112 
0113   constexpr auto tolerance = s_onSurfaceTolerance;
0114   auto same = [](auto a, auto b) { return std::abs(a - b) < tolerance; };
0115 
0116   auto checkZ = [&cyl, same](const IAxis& axis) {
0117     double hlZ = cyl.bounds().get(CylinderBounds::eHalfLengthZ);
0118     if (!same(axis.getMin(), -hlZ) || !same(axis.getMax(), hlZ)) {
0119       throw std::invalid_argument(
0120           "GridPortalLink: CylinderBounds: invalid length setup: " +
0121           std::to_string(axis.getMin()) + " != " + std::to_string(-hlZ) +
0122           " or " + std::to_string(axis.getMax()) +
0123           " != " + std::to_string(hlZ));
0124     }
0125   };
0126   auto checkRPhi = [&cyl, same](const IAxis& axis) {
0127     double hlPhi = cyl.bounds().get(CylinderBounds::eHalfPhiSector);
0128     double r = cyl.bounds().get(CylinderBounds::eR);
0129     if (double hlRPhi = r * hlPhi;
0130         !same(axis.getMin(), -hlRPhi) || !same(axis.getMax(), hlRPhi)) {
0131       throw std::invalid_argument(
0132           "GridPortalLink: CylinderBounds: invalid phi sector setup: axes "
0133           "don't match bounds");
0134     }
0135 
0136     // If full cylinder, make sure axis wraps around
0137     if (same(hlPhi, std::numbers::pi)) {
0138       if (axis.getBoundaryType() != AxisBoundaryType::Closed) {
0139         throw std::invalid_argument(
0140             "GridPortalLink: CylinderBounds: invalid phi sector setup: "
0141             "axis is not closed.");
0142       }
0143     } else {
0144       if (axis.getBoundaryType() == AxisBoundaryType::Closed) {
0145         throw std::invalid_argument(
0146             "GridPortalLink: CylinderBounds: invalid phi sector setup: "
0147             "axis is closed.");
0148       }
0149     }
0150   };
0151 
0152   if (grid.axes().size() == 1) {
0153     const IAxis& axisLoc0 = *grid.axes().front();
0154     if (direction == AxisDirection::AxisRPhi) {
0155       checkRPhi(axisLoc0);
0156     } else {
0157       checkZ(axisLoc0);
0158     }
0159   } else {  // DIM == 2
0160     const auto& axisLoc0 = *grid.axes().front();
0161     const auto& axisLoc1 = *grid.axes().back();
0162     checkRPhi(axisLoc0);
0163     checkZ(axisLoc1);
0164   }
0165 }
0166 
0167 void GridPortalLink::checkConsistency(const IGrid& grid,
0168                                       AxisDirection direction,
0169                                       const DiscSurface& disc) {
0170   constexpr auto tolerance = s_onSurfaceTolerance;
0171   auto same = [](auto a, auto b) { return std::abs(a - b) < tolerance; };
0172 
0173   const auto* bounds = dynamic_cast<const RadialBounds*>(&disc.bounds());
0174   if (bounds == nullptr) {
0175     throw std::invalid_argument(
0176         "GridPortalLink: DiscBounds: invalid bounds type.");
0177   }
0178 
0179   if (bounds->get(RadialBounds::eAveragePhi) != 0) {
0180     throw std::invalid_argument(
0181         "GridPortalLink: DiscBounds: only average phi == 0 is supported. "
0182         "Rotate the disc surface.");
0183   }
0184 
0185   auto checkR = [&bounds, same](const IAxis& axis) {
0186     double minR = bounds->get(RadialBounds::eMinR);
0187     double maxR = bounds->get(RadialBounds::eMaxR);
0188     if (!same(axis.getMin(), minR) || !same(axis.getMax(), maxR)) {
0189       throw std::invalid_argument(
0190           "GridPortalLink: DiscBounds: invalid radius setup.");
0191     }
0192   };
0193 
0194   auto checkPhi = [&bounds, same](const IAxis& axis) {
0195     double hlPhi = bounds->get(RadialBounds::eHalfPhiSector);
0196     if (!same(axis.getMin(), -hlPhi) || !same(axis.getMax(), hlPhi)) {
0197       throw std::invalid_argument(
0198           "GridPortalLink: DiscBounds: invalid phi sector setup.");
0199     }
0200     // If full disc, make sure axis wraps around
0201     if (same(hlPhi, std::numbers::pi)) {
0202       if (axis.getBoundaryType() != Acts::AxisBoundaryType::Closed) {
0203         throw std::invalid_argument(
0204             "GridPortalLink: DiscBounds: invalid phi sector setup: axis is "
0205             "not closed.");
0206       }
0207     } else {
0208       if (axis.getBoundaryType() == Acts::AxisBoundaryType::Closed) {
0209         throw std::invalid_argument(
0210             "GridPortalLink: DiscBounds: invalid phi sector setup: axis "
0211             "is closed.");
0212       }
0213     }
0214   };
0215 
0216   if (grid.axes().size() == 1) {
0217     const IAxis& axisLoc0 = *grid.axes().front();
0218     if (direction == AxisDirection::AxisR) {
0219       checkR(axisLoc0);
0220     } else {
0221       checkPhi(axisLoc0);
0222     }
0223   } else {  // DIM == 2
0224     const auto& axisLoc0 = *grid.axes().front();
0225     const auto& axisLoc1 = *grid.axes().back();
0226     checkR(axisLoc0);
0227     checkPhi(axisLoc1);
0228   }
0229 }
0230 
0231 void GridPortalLink::checkConsistency(const IGrid& grid,
0232                                       AxisDirection direction,
0233                                       const PlaneSurface& plane) {
0234   constexpr auto tolerance = s_onSurfaceTolerance;
0235   auto same = [](auto a, auto b) { return std::abs(a - b) < tolerance; };
0236 
0237   const auto* bounds = dynamic_cast<const RectangleBounds*>(&plane.bounds());
0238   if (bounds == nullptr) {
0239     throw std::invalid_argument(
0240         "GridPortalLink: PlaneBounds: invalid bounds type.");
0241   }
0242   auto check = [&bounds, same](const IAxis& axis, AxisDirection dir) {
0243     double min = (dir == AxisDirection::AxisX)
0244                      ? bounds->get(RectangleBounds::eMinX)
0245                      : bounds->get(RectangleBounds::eMinY);
0246     double max = (dir == AxisDirection::AxisX)
0247                      ? bounds->get(RectangleBounds::eMaxX)
0248                      : bounds->get(RectangleBounds::eMaxY);
0249     if (!same(axis.getMin(), min) || !same(axis.getMax(), max)) {
0250       throw std::invalid_argument(
0251           "GridPortalLink: PlaneBounds: invalid setup.");
0252     }
0253   };
0254 
0255   if (grid.axes().size() == 1) {
0256     const IAxis& axisLoc0 = *grid.axes().front();
0257     check(axisLoc0, direction);
0258   } else {  // DIM == 2
0259     const auto& axisLoc0 = *grid.axes().front();
0260     const auto& axisLoc1 = *grid.axes().back();
0261     check(axisLoc0, AxisDirection::AxisX);
0262     check(axisLoc1, AxisDirection::AxisY);
0263   }
0264 }
0265 
0266 void GridPortalLink::printContents(std::ostream& os) const {
0267   std::size_t dim = grid().axes().size();
0268   os << "----- GRID " << dim << "d -----" << std::endl;
0269   os << grid() << " along " << direction() << std::endl;
0270 
0271   auto lpad = [](const std::string& s, std::size_t n) {
0272     return s.size() < n ? std::string(n - s.size(), ' ') + s : s;
0273   };
0274 
0275   auto rpad = [](const std::string& s, std::size_t n) {
0276     return s.size() < n ? s + std::string(n - s.size(), ' ') : s;
0277   };
0278 
0279   std::string loc0;
0280   std::string loc1;
0281 
0282   bool flipped = false;
0283 
0284   if (surface().type() == Surface::Cylinder) {
0285     loc0 = "rPhi";
0286     loc1 = "z";
0287     flipped = direction() != AxisDirection::AxisRPhi;
0288   } else if (surface().type() == Surface::Disc) {
0289     loc0 = "r";
0290     loc1 = "phi";
0291     flipped = direction() != AxisDirection::AxisR;
0292   } else if (surface().type() == Surface::Plane) {
0293     loc0 = "x";
0294     loc1 = "y";
0295     flipped = direction() != AxisDirection::AxisX;
0296   } else {
0297     throw std::invalid_argument{"Unsupported surface type"};
0298   }
0299 
0300   AnyGridConstView<const TrackingVolume*> view(grid());
0301 
0302   if (dim == 1) {
0303     auto loc = grid().numLocalBinsAny();
0304 
0305     if (flipped) {
0306       os << lpad(loc1, 4) << " > " << lpad("i=0", 10) << " ";
0307       for (std::size_t i = 1; i <= loc.at(0) + 1; i++) {
0308         os << lpad("i=" + std::to_string(i), 13) + " ";
0309       }
0310       os << std::endl;
0311 
0312       os << std::string(4, ' ');
0313       for (std::size_t i = 0; i <= loc.at(0) + 1; i++) {
0314         std::string name = "0x0";
0315         if (const auto* v = view.atLocalBins({i}); v != nullptr) {
0316           name = v->volumeName();
0317         }
0318         name = name.substr(0, std::min(name.size(), std::size_t{13}));
0319         os << lpad(name, 13) << " ";
0320       }
0321       os << std::endl;
0322 
0323     } else {
0324       os << "v " << loc0 << std::endl;
0325       for (std::size_t i = 0; i <= loc.at(0) + 1; i++) {
0326         os << "i=" << i << " ";
0327         std::string name = "0x0";
0328         if (const auto* v = view.atLocalBins({i}); v != nullptr) {
0329           name = v->volumeName();
0330         }
0331         name = name.substr(0, std::min(name.size(), std::size_t{13}));
0332         os << lpad(name, 13) << " ";
0333         os << std::endl;
0334       }
0335     }
0336 
0337   } else {
0338     auto loc = grid().numLocalBinsAny();
0339     os << rpad("v " + loc0 + "|" + loc1 + " >", 14) + "j=0 ";
0340     for (std::size_t j = 1; j <= loc.at(1) + 1; j++) {
0341       os << lpad("j=" + std::to_string(j), 13) << " ";
0342     }
0343     os << std::endl;
0344     for (std::size_t i = 0; i <= loc.at(0) + 1; i++) {
0345       os << "i=" << i << " ";
0346       for (std::size_t j = 0; j <= loc.at(1) + 1; j++) {
0347         std::string name = "0x0";
0348         if (const auto* v = view.atLocalBins({i, j}); v != nullptr) {
0349           name = v->volumeName();
0350         }
0351         name = name.substr(0, std::min(name.size(), std::size_t{13}));
0352         os << lpad(name, 13) << " ";
0353       }
0354       os << std::endl;
0355     }
0356   }
0357 }
0358 
0359 void GridPortalLink::fillGrid1dTo2d(FillDirection dir,
0360                                     const GridPortalLink& grid1d,
0361                                     GridPortalLink& grid2d) {
0362   const auto locSource = grid1d.grid().numLocalBinsAny();
0363   const auto locDest = grid2d.grid().numLocalBinsAny();
0364   assert(grid1d.grid().dimensions() == 1);
0365   assert(grid2d.grid().dimensions() == 2);
0366   assert(locSource.size() == 1);
0367   assert(locDest.size() == 2);
0368 
0369   AnyGridConstView<const TrackingVolume*> sourceView(grid1d.grid());
0370   AnyGridView<const TrackingVolume*> destView(grid2d.grid());
0371 
0372   for (std::size_t i = 0; i <= locSource[0] + 1; ++i) {
0373     const TrackingVolume* source = sourceView.atLocalBins({i});
0374 
0375     if (dir == FillDirection::loc1) {
0376       for (std::size_t j = 0; j <= locDest[1] + 1; ++j) {
0377         destView.atLocalBins({i, j}) = source;
0378       }
0379     } else if (dir == FillDirection::loc0) {
0380       for (std::size_t j = 0; j <= locDest[0] + 1; ++j) {
0381         destView.atLocalBins({j, i}) = source;
0382       }
0383     }
0384   }
0385 }
0386 
0387 std::unique_ptr<GridPortalLink> GridPortalLink::extendTo2dImpl(
0388     const std::shared_ptr<CylinderSurface>& surface, const IAxis* other) const {
0389   assert(dim() == 1);
0390   if (direction() == AxisDirection::AxisRPhi) {
0391     const auto& axisRPhi = *grid().axes().front();
0392     // 1D direction is AxisRPhi, so add a Z axis
0393     double hlZ = surface->bounds().get(CylinderBounds::eHalfLengthZ);
0394 
0395     auto grid = axisRPhi.visit([&](const auto& axis0) {
0396       Axis axisZ{AxisBound, -hlZ, hlZ, 1};
0397       const IAxis* axis = other != nullptr ? other : &axisZ;
0398       return axis->visit(
0399           [&surface,
0400            &axis0](const auto& axis1) -> std::unique_ptr<GridPortalLink> {
0401             return GridPortalLink::make(surface, axis0, axis1);
0402           });
0403     });
0404 
0405     fillGrid1dTo2d(FillDirection::loc1, *this, *grid);
0406     return grid;
0407 
0408   } else {
0409     const auto& axisZ = *grid().axes().front();
0410     // 1D direction is AxisZ, so add an rPhi axis
0411     double r = surface->bounds().get(CylinderBounds::eR);
0412     double hlPhi = surface->bounds().get(CylinderBounds::eHalfPhiSector);
0413     double hlRPhi = r * hlPhi;
0414 
0415     auto makeGrid = [&](auto bdt) {
0416       auto grid = axisZ.visit([&](const auto& axis1) {
0417         Axis axisRPhi{bdt, -hlRPhi, hlRPhi, 1};
0418         const IAxis* axis = other != nullptr ? other : &axisRPhi;
0419         return axis->visit(
0420             [&](const auto& axis0) -> std::unique_ptr<GridPortalLink> {
0421               return GridPortalLink::make(surface, axis0, axis1);
0422             });
0423       });
0424 
0425       fillGrid1dTo2d(FillDirection::loc0, *this, *grid);
0426       return grid;
0427     };
0428 
0429     if (surface->bounds().coversFullAzimuth()) {
0430       return makeGrid(AxisClosed);
0431     } else {
0432       return makeGrid(AxisBound);
0433     }
0434   }
0435 }
0436 
0437 std::unique_ptr<GridPortalLink> GridPortalLink::extendTo2dImpl(
0438     const std::shared_ptr<DiscSurface>& surface, const IAxis* other) const {
0439   assert(dim() == 1);
0440 
0441   const auto* bounds = dynamic_cast<const RadialBounds*>(&surface->bounds());
0442   if (bounds == nullptr) {
0443     throw std::invalid_argument(
0444         "GridPortalLink: DiscBounds: invalid bounds type.");
0445   }
0446 
0447   if (direction() == AxisDirection::AxisR) {
0448     const auto& axisR = *grid().axes().front();
0449     // 1D direction is AxisR, so add a phi axis
0450     double hlPhi = bounds->get(RadialBounds::eHalfPhiSector);
0451 
0452     auto makeGrid = [&](auto bdt) {
0453       auto grid = axisR.visit([&](const auto& axis0) {
0454         Axis axisPhi{bdt, -hlPhi, hlPhi, 1};
0455         const IAxis* axis = other != nullptr ? other : &axisPhi;
0456         return axis->visit(
0457             [&](const auto& axis1) -> std::unique_ptr<GridPortalLink> {
0458               return GridPortalLink::make(surface, axis0, axis1);
0459             });
0460       });
0461 
0462       fillGrid1dTo2d(FillDirection::loc1, *this, *grid);
0463       return grid;
0464     };
0465 
0466     if (bounds->coversFullAzimuth()) {
0467       return makeGrid(AxisClosed);
0468     } else {
0469       return makeGrid(AxisBound);
0470     }
0471   } else {
0472     const auto& axisPhi = *grid().axes().front();
0473     // 1D direction is AxisPhi, so add an R axis
0474     double rMin = bounds->get(RadialBounds::eMinR);
0475     double rMax = bounds->get(RadialBounds::eMaxR);
0476 
0477     auto grid = axisPhi.visit([&](const auto& axis1) {
0478       Axis axisR{AxisBound, rMin, rMax, 1};
0479       const IAxis* axis = other != nullptr ? other : &axisR;
0480       return axis->visit(
0481           [&surface,
0482            &axis1](const auto& axis0) -> std::unique_ptr<GridPortalLink> {
0483             return GridPortalLink::make(surface, axis0, axis1);
0484           });
0485     });
0486     fillGrid1dTo2d(FillDirection::loc0, *this, *grid);
0487     return grid;
0488   }
0489 }
0490 
0491 std::unique_ptr<GridPortalLink> GridPortalLink::extendTo2dImpl(
0492     const std::shared_ptr<PlaneSurface>& surface, const IAxis* other) const {
0493   assert(dim() == 1);
0494 
0495   const auto* bounds = dynamic_cast<const RectangleBounds*>(&surface->bounds());
0496   if (bounds == nullptr) {
0497     throw std::invalid_argument(
0498         "GridPortalLink: RectangleBounds: invalid bounds type.");
0499   }
0500 
0501   bool dirX = direction() == AxisDirection::AxisX;
0502   const auto& thisAxis = *grid().axes().front();
0503 
0504   double minExtend = dirX ? bounds->get(RectangleBounds::BoundValues::eMinY)
0505                           : bounds->get(RectangleBounds::BoundValues::eMinX);
0506   double maxExtend = dirX ? bounds->get(RectangleBounds::BoundValues::eMaxY)
0507                           : bounds->get(RectangleBounds::BoundValues::eMaxX);
0508 
0509   FillDirection fillDir = dirX ? FillDirection::loc1 : FillDirection::loc0;
0510 
0511   auto grid = thisAxis.visit([&](const auto& axis0) {
0512     Axis axisExtend{AxisBound, minExtend, maxExtend, 1};
0513     const IAxis* axis = other != nullptr ? other : &axisExtend;
0514     return axis->visit(
0515         [&](const auto& axis1) -> std::unique_ptr<GridPortalLink> {
0516           if (dirX) {
0517             return GridPortalLink::make(surface, axis0, axis1);
0518           } else {
0519             return GridPortalLink::make(surface, axis1, axis0);
0520           }
0521         });
0522   });
0523 
0524   fillGrid1dTo2d(fillDir, *this, *grid);
0525   return grid;
0526 }
0527 
0528 std::span<const TrivialPortalLink> GridPortalLink::artifactPortalLinks() const {
0529   return {m_artifactPortalLinks};
0530 }
0531 
0532 void GridPortalLink::setArtifactPortalLinks(
0533     std::vector<TrivialPortalLink> links) {
0534   m_artifactPortalLinks = std::move(links);
0535 }
0536 
0537 }  // namespace Acts
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!