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