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/Surfaces/PlaneSurface.hpp"
0012 #include "Acts/Surfaces/RadialBounds.hpp"
0013
0014 #include <numbers>
0015
0016 namespace Acts {
0017
0018 std::unique_ptr<GridPortalLink> GridPortalLink::make(
0019     const std::shared_ptr<RegularSurface>& surface, TrackingVolume& volume,
0020     AxisDirection direction) {
0021   std::unique_ptr<GridPortalLink> grid;
0022
0023   if (const auto* cylinder =
0024           dynamic_cast<const CylinderSurface*>(surface.get());
0025       cylinder != nullptr) {
0026     if (direction == AxisDirection::AxisRPhi) {
0027       double r = cylinder->bounds().get(CylinderBounds::eR);
0028       if (cylinder->bounds().coversFullAzimuth()) {
0029         grid = GridPortalLink::make(
0030             surface, direction,
0031             Axis{AxisClosed, -std::numbers::pi * r, std::numbers::pi * r, 1});
0032       } else {
0033         double hlPhi = cylinder->bounds().get(CylinderBounds::eHalfPhiSector);
0034
0035         grid = GridPortalLink::make(surface, direction,
0036                                     Axis{AxisBound, -hlPhi * r, hlPhi * r, 1});
0037       }
0038     } else if (direction == AxisDirection::AxisZ) {
0039       double hlZ = cylinder->bounds().get(CylinderBounds::eHalfLengthZ);
0040       grid = GridPortalLink::make(surface, direction,
0041                                   Axis{AxisBound, -hlZ, hlZ, 1});
0042     } else {
0043       throw std::invalid_argument{"Invalid binning direction"};
0044     }
0045   } else if (const auto* disc = dynamic_cast<const DiscSurface*>(surface.get());
0046              disc != nullptr) {
0047     const auto& bounds = dynamic_cast<const RadialBounds&>(disc->bounds());
0048     if (direction == AxisDirection::AxisR) {
0049       double minR = bounds.get(RadialBounds::eMinR);
0050       double maxR = bounds.get(RadialBounds::eMaxR);
0051       grid = GridPortalLink::make(surface, direction,
0052                                   Axis{AxisBound, minR, maxR, 1});
0053     } else if (direction == AxisDirection::AxisPhi) {
0054       if (bounds.coversFullAzimuth()) {
0055         grid = GridPortalLink::make(
0056             surface, direction,
0057             Axis{AxisClosed, -std::numbers::pi, std::numbers::pi, 1});
0058       } else {
0059         double hlPhi = bounds.get(RadialBounds::eHalfPhiSector);
0060         grid = GridPortalLink::make(surface, direction,
0061                                     Axis{AxisBound, -hlPhi, hlPhi, 1});
0062       }
0063     } else {
0064       throw std::invalid_argument{"Invalid binning direction"};
0065     }
0066   } else if (const auto* plane =
0067                  dynamic_cast<const PlaneSurface*>(surface.get());
0068              plane != nullptr) {
0069     throw std::invalid_argument{"Plane surface is not implemented yet"};
0070
0071   } else {
0072     throw std::invalid_argument{"Surface type is not supported"};
0073   }
0074
0075   assert(grid != nullptr);
0076   grid->setVolume(&volume);
0077
0078   return grid;
0079 }
0080
0081 void GridPortalLink::checkConsistency(const CylinderSurface& cyl) const {
0082   if (cyl.bounds().get(CylinderBounds::eAveragePhi) != 0) {
0083     throw std::invalid_argument(
0084         "GridPortalLink: CylinderBounds: only average phi == 0 is "
0085         "supported. Rotate the cylinder surface.");
0086   }
0087
0088   constexpr auto tolerance = s_onSurfaceTolerance;
0089   auto same = [](auto a, auto b) { return std::abs(a - b) < tolerance; };
0090
0091   auto checkZ = [&cyl, same](const IAxis& axis) {
0092     double hlZ = cyl.bounds().get(CylinderBounds::eHalfLengthZ);
0093     if (!same(axis.getMin(), -hlZ) || !same(axis.getMax(), hlZ)) {
0094       throw std::invalid_argument(
0095           "GridPortalLink: CylinderBounds: invalid length setup: " +
0096           std::to_string(axis.getMin()) + " != " + std::to_string(-hlZ) +
0097           " or " + std::to_string(axis.getMax()) +
0098           " != " + std::to_string(hlZ));
0099     }
0100   };
0101   auto checkRPhi = [&cyl, same](const IAxis& axis) {
0102     double hlPhi = cyl.bounds().get(CylinderBounds::eHalfPhiSector);
0103     double r = cyl.bounds().get(CylinderBounds::eR);
0104     if (double hlRPhi = r * hlPhi;
0105         !same(axis.getMin(), -hlRPhi) || !same(axis.getMax(), hlRPhi)) {
0106       throw std::invalid_argument(
0107           "GridPortalLink: CylinderBounds: invalid phi sector setup: axes "
0108           "don't match bounds");
0109     }
0110
0111     // If full cylinder, make sure axis wraps around
0112     if (same(hlPhi, std::numbers::pi)) {
0113       if (axis.getBoundaryType() != AxisBoundaryType::Closed) {
0114         throw std::invalid_argument(
0115             "GridPortalLink: CylinderBounds: invalid phi sector setup: "
0116             "axis is not closed.");
0117       }
0118     } else {
0119       if (axis.getBoundaryType() == AxisBoundaryType::Closed) {
0120         throw std::invalid_argument(
0121             "GridPortalLink: CylinderBounds: invalid phi sector setup: "
0122             "axis is closed.");
0123       }
0124     }
0125   };
0126
0127   if (dim() == 1) {
0128     const IAxis& axisLoc0 = *grid().axes().front();
0129     if (direction() == AxisDirection::AxisRPhi) {
0130       checkRPhi(axisLoc0);
0131     } else {
0132       checkZ(axisLoc0);
0133     }
0134   } else {  // DIM == 2
0135     const auto& axisLoc0 = *grid().axes().front();
0136     const auto& axisLoc1 = *grid().axes().back();
0137     checkRPhi(axisLoc0);
0138     checkZ(axisLoc1);
0139   }
0140 }
0141
0142 void GridPortalLink::checkConsistency(const DiscSurface& disc) const {
0143   constexpr auto tolerance = s_onSurfaceTolerance;
0144   auto same = [](auto a, auto b) { return std::abs(a - b) < tolerance; };
0145
0146   const auto* bounds = dynamic_cast<const RadialBounds*>(&disc.bounds());
0147   if (bounds == nullptr) {
0148     throw std::invalid_argument(
0149         "GridPortalLink: DiscBounds: invalid bounds type.");
0150   }
0151
0152   if (bounds->get(RadialBounds::eAveragePhi) != 0) {
0153     throw std::invalid_argument(
0154         "GridPortalLink: DiscBounds: only average phi == 0 is supported. "
0155         "Rotate the disc surface.");
0156   }
0157
0158   auto checkR = [&bounds, same](const IAxis& axis) {
0159     double minR = bounds->get(RadialBounds::eMinR);
0160     double maxR = bounds->get(RadialBounds::eMaxR);
0161     if (!same(axis.getMin(), minR) || !same(axis.getMax(), maxR)) {
0162       throw std::invalid_argument(
0163           "GridPortalLink: DiscBounds: invalid radius setup.");
0164     }
0165   };
0166
0167   auto checkPhi = [&bounds, same](const IAxis& axis) {
0168     double hlPhi = bounds->get(RadialBounds::eHalfPhiSector);
0169     if (!same(axis.getMin(), -hlPhi) || !same(axis.getMax(), hlPhi)) {
0170       throw std::invalid_argument(
0171           "GridPortalLink: DiscBounds: invalid phi sector setup.");
0172     }
0173     // If full disc, make sure axis wraps around
0174     if (same(hlPhi, std::numbers::pi)) {
0175       if (axis.getBoundaryType() != Acts::AxisBoundaryType::Closed) {
0176         throw std::invalid_argument(
0177             "GridPortalLink: DiscBounds: invalid phi sector setup: axis is "
0178             "not closed.");
0179       }
0180     } else {
0181       if (axis.getBoundaryType() == Acts::AxisBoundaryType::Closed) {
0182         throw std::invalid_argument(
0183             "GridPortalLink: DiscBounds: invalid phi sector setup: axis "
0184             "is closed.");
0185       }
0186     }
0187   };
0188
0189   if (dim() == 1) {
0190     const IAxis& axisLoc0 = *grid().axes().front();
0191     if (direction() == AxisDirection::AxisR) {
0192       checkR(axisLoc0);
0193     } else {
0194       checkPhi(axisLoc0);
0195     }
0196   } else {  // DIM == 2
0197     const auto& axisLoc0 = *grid().axes().front();
0198     const auto& axisLoc1 = *grid().axes().back();
0199     checkR(axisLoc0);
0200     checkPhi(axisLoc1);
0201   }
0202 }
0203
0204 void GridPortalLink::printContents(std::ostream& os) const {
0205   std::size_t dim = grid().axes().size();
0206   os << "----- GRID " << dim << "d -----" << std::endl;
0207   os << grid() << " along " << direction() << std::endl;
0208
0209   auto lpad = [](const std::string& s, std::size_t n) {
0210     return s.size() < n ? std::string(n - s.size(), ' ') + s : s;
0211   };
0212
0213   auto rpad = [](const std::string& s, std::size_t n) {
0214     return s.size() < n ? s + std::string(n - s.size(), ' ') : s;
0215   };
0216
0217   std::string loc0;
0218   std::string loc1;
0219
0220   bool flipped = false;
0221
0222   if (surface().type() == Surface::Cylinder) {
0223     loc0 = "rPhi";
0224     loc1 = "z";
0225     flipped = direction() != AxisDirection::AxisRPhi;
0226   } else if (surface().type() == Surface::Disc) {
0227     loc0 = "r";
0228     loc1 = "phi";
0229     flipped = direction() != AxisDirection::AxisR;
0230   } else if (surface().type() == Surface::Plane) {
0231     loc0 = "x";
0232     loc1 = "y";
0233     flipped = direction() != AxisDirection::AxisX;
0234   } else {
0235     throw std::invalid_argument{"Unsupported surface type"};
0236   }
0237
0238   if (dim == 1) {
0239     auto loc = numLocalBins();
0240
0241     if (flipped) {
0242       os << lpad(loc1, 4) << " > " << lpad("i=0", 10) << " ";
0243       for (std::size_t i = 1; i <= loc.at(0) + 1; i++) {
0244         os << lpad("i=" + std::to_string(i), 13) + " ";
0245       }
0246       os << std::endl;
0247
0248       os << std::string(4, ' ');
0249       for (std::size_t i = 0; i <= loc.at(0) + 1; i++) {
0250         std::string name = "0x0";
0251         if (const auto* v = atLocalBins({i}); v != nullptr) {
0252           name = v->volumeName();
0253         }
0254         name = name.substr(0, std::min(name.size(), std::size_t{13}));
0255         os << lpad(name, 13) << " ";
0256       }
0257       os << std::endl;
0258
0259     } else {
0260       os << "v " << loc0 << std::endl;
0261       for (std::size_t i = 0; i <= loc.at(0) + 1; i++) {
0262         os << "i=" << i << " ";
0263         std::string name = "0x0";
0264         if (const auto* v = atLocalBins({i}); v != nullptr) {
0265           name = v->volumeName();
0266         }
0267         name = name.substr(0, std::min(name.size(), std::size_t{13}));
0268         os << lpad(name, 13) << " ";
0269         os << std::endl;
0270       }
0271     }
0272
0273   } else {
0274     auto loc = numLocalBins();
0275     os << rpad("v " + loc0 + "|" + loc1 + " >", 14) + "j=0 ";
0276     for (std::size_t j = 1; j <= loc.at(1) + 1; j++) {
0277       os << lpad("j=" + std::to_string(j), 13) << " ";
0278     }
0279     os << std::endl;
0280     for (std::size_t i = 0; i <= loc.at(0) + 1; i++) {
0281       os << "i=" << i << " ";
0282       for (std::size_t j = 0; j <= loc.at(1) + 1; j++) {
0283         std::string name = "0x0";
0284         if (const auto* v = atLocalBins({i, j}); v != nullptr) {
0285           name = v->volumeName();
0286         }
0287         name = name.substr(0, std::min(name.size(), std::size_t{13}));
0288         os << lpad(name, 13) << " ";
0289       }
0290       os << std::endl;
0291     }
0292   }
0293 }
0294
0295 void GridPortalLink::fillGrid1dTo2d(FillDirection dir,
0296                                     const GridPortalLink& grid1d,
0297                                     GridPortalLink& grid2d) {
0298   const auto locSource = grid1d.numLocalBins();
0299   const auto locDest = grid2d.numLocalBins();
0300   assert(locSource.size() == 1);
0301   assert(locDest.size() == 2);
0302
0303   for (std::size_t i = 0; i <= locSource[0] + 1; ++i) {
0304     const TrackingVolume* source = grid1d.atLocalBins({i});
0305
0306     if (dir == FillDirection::loc1) {
0307       for (std::size_t j = 0; j <= locDest[1] + 1; ++j) {
0308         grid2d.atLocalBins({i, j}) = source;
0309       }
0310     } else if (dir == FillDirection::loc0) {
0311       for (std::size_t j = 0; j <= locDest[0] + 1; ++j) {
0312         grid2d.atLocalBins({j, i}) = source;
0313       }
0314     }
0315   }
0316 }
0317
0318 std::unique_ptr<GridPortalLink> GridPortalLink::extendTo2dImpl(
0319     const std::shared_ptr<CylinderSurface>& surface, const IAxis* other) const {
0320   assert(dim() == 1);
0321   if (direction() == AxisDirection::AxisRPhi) {
0322     const auto& axisRPhi = *grid().axes().front();
0323     // 1D direction is AxisRPhi, so add a Z axis
0324     double hlZ = surface->bounds().get(CylinderBounds::eHalfLengthZ);
0325
0326     auto grid = axisRPhi.visit([&](const auto& axis0) {
0327       Axis axisZ{AxisBound, -hlZ, hlZ, 1};
0328       const IAxis* axis = other != nullptr ? other : &axisZ;
0329       return axis->visit(
0330           [&surface,
0331            &axis0](const auto& axis1) -> std::unique_ptr<GridPortalLink> {
0332             return GridPortalLink::make(surface, axis0, axis1);
0333           });
0334     });
0335
0336     fillGrid1dTo2d(FillDirection::loc1, *this, *grid);
0337     return grid;
0338
0339   } else {
0340     const auto& axisZ = *grid().axes().front();
0341     // 1D direction is AxisZ, so add an rPhi axis
0342     double r = surface->bounds().get(CylinderBounds::eR);
0343     double hlPhi = surface->bounds().get(CylinderBounds::eHalfPhiSector);
0344     double hlRPhi = r * hlPhi;
0345
0346     auto makeGrid = [&](auto bdt) {
0347       auto grid = axisZ.visit([&](const auto& axis1) {
0348         Axis axisRPhi{bdt, -hlRPhi, hlRPhi, 1};
0349         const IAxis* axis = other != nullptr ? other : &axisRPhi;
0350         return axis->visit(
0351             [&](const auto& axis0) -> std::unique_ptr<GridPortalLink> {
0352               return GridPortalLink::make(surface, axis0, axis1);
0353             });
0354       });
0355
0356       fillGrid1dTo2d(FillDirection::loc0, *this, *grid);
0357       return grid;
0358     };
0359
0360     if (surface->bounds().coversFullAzimuth()) {
0361       return makeGrid(AxisClosed);
0362     } else {
0363       return makeGrid(AxisBound);
0364     }
0365   }
0366 }
0367
0368 std::unique_ptr<GridPortalLink> GridPortalLink::extendTo2dImpl(
0369     const std::shared_ptr<DiscSurface>& surface, const IAxis* other) const {
0370   assert(dim() == 1);
0371
0372   const auto* bounds = dynamic_cast<const RadialBounds*>(&surface->bounds());
0373   if (bounds == nullptr) {
0374     throw std::invalid_argument(
0375         "GridPortalLink: DiscBounds: invalid bounds type.");
0376   }
0377
0378   if (direction() == AxisDirection::AxisR) {
0379     const auto& axisR = *grid().axes().front();
0380     // 1D direction is AxisR, so add a phi axis
0381     double hlPhi = bounds->get(RadialBounds::eHalfPhiSector);
0382
0383     auto makeGrid = [&](auto bdt) {
0384       auto grid = axisR.visit([&](const auto& axis0) {
0385         Axis axisPhi{bdt, -hlPhi, hlPhi, 1};
0386         const IAxis* axis = other != nullptr ? other : &axisPhi;
0387         return axis->visit(
0388             [&](const auto& axis1) -> std::unique_ptr<GridPortalLink> {
0389               return GridPortalLink::make(surface, axis0, axis1);
0390             });
0391       });
0392
0393       fillGrid1dTo2d(FillDirection::loc1, *this, *grid);
0394       return grid;
0395     };
0396
0397     if (bounds->coversFullAzimuth()) {
0398       return makeGrid(AxisClosed);
0399     } else {
0400       return makeGrid(AxisBound);
0401     }
0402   } else {
0403     const auto& axisPhi = *grid().axes().front();
0404     // 1D direction is AxisPhi, so add an R axis
0405     double rMin = bounds->get(RadialBounds::eMinR);
0406     double rMax = bounds->get(RadialBounds::eMaxR);
0407
0408     auto grid = axisPhi.visit([&](const auto& axis1) {
0409       Axis axisR{AxisBound, rMin, rMax, 1};
0410       const IAxis* axis = other != nullptr ? other : &axisR;
0411       return axis->visit(
0412           [&surface,
0413            &axis1](const auto& axis0) -> std::unique_ptr<GridPortalLink> {
0414             return GridPortalLink::make(surface, axis0, axis1);
0415           });
0416     });
0417     fillGrid1dTo2d(FillDirection::loc0, *this, *grid);
0418     return grid;
0419   }
0420 }
0421
0422 }  // namespace Acts