Core/Surfaces/SurfaceArrayTests.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 <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/GeometryContext.hpp"
0013 #include "Acts/Surfaces/CylinderSurface.hpp"
0014 #include "Acts/Surfaces/DiscSurface.hpp"
0015 #include "Acts/Surfaces/PlanarBounds.hpp"
0016 #include "Acts/Surfaces/PlaneSurface.hpp"
0017 #include "Acts/Surfaces/RectangleBounds.hpp"
0018 #include "Acts/Surfaces/Surface.hpp"
0019 #include "Acts/Surfaces/SurfaceArray.hpp"
0020 #include "Acts/Utilities/Axis.hpp"
0021 #include "Acts/Utilities/AxisDefinitions.hpp"
0022 #include "Acts/Utilities/Diagnostics.hpp"
0023 #include "Acts/Utilities/Helpers.hpp"
0024
0025 #include <algorithm>
0026 #include <array>
0027 #include <cfenv>
0028 #include <cmath>
0029 #include <cstddef>
0030 #include <cstdint>
0031 #include <fstream>
0032 #include <iomanip>
0033 #include <iostream>
0034 #include <memory>
0035 #include <numbers>
0036 #include <ranges>
0037 #include <sstream>
0038 #include <string>
0039 #include <tuple>
0040 #include <utility>
0041 #include <vector>
0042
0043 #include <boost/format.hpp>
0044
0045 using Acts::VectorHelpers::phi;
0046
0047 using namespace Acts;
0048
0049 namespace ActsTests {
0050
0051 // Create a test context
0052 GeometryContext tgContext = GeometryContext::dangerouslyDefaultConstruct();
0053
0054 using SrfVec = std::vector<std::shared_ptr<const Surface>>;
0055 struct SurfaceArrayFixture {
0056   std::vector<std::shared_ptr<const Surface>> m_surfaces;
0057
0058   SurfaceArrayFixture() { BOOST_TEST_MESSAGE("setup fixture"); }
0059   ~SurfaceArrayFixture() { BOOST_TEST_MESSAGE("teardown fixture"); }
0060
0061   SrfVec fullPhiTestSurfacesEC(std::size_t n = 10, double shift = 0,
0062                                double zbase = 0, double r = 10) {
0063     SrfVec res;
0064
0065     double phiStep = 2 * std::numbers::pi / n;
0066     for (std::size_t i = 0; i < n; ++i) {
0067       double z = zbase + ((i % 2 == 0) ? 1 : -1) * 0.2;
0068
0069       Transform3 trans;
0070       trans.setIdentity();
0071       trans.rotate(Eigen::AngleAxisd(i * phiStep + shift, Vector3(0, 0, 1)));
0072       trans.translate(Vector3(r, 0, z));
0073
0074       auto bounds = std::make_shared<const RectangleBounds>(2, 1);
0075       std::shared_ptr<const Surface> srf =
0076           Surface::makeShared<PlaneSurface>(trans, bounds);
0077
0078       res.push_back(srf);
0079       m_surfaces.push_back(
0080           std::move(srf));  // keep shared, will get destroyed at the end
0081     }
0082
0083     return res;
0084   }
0085
0086   SrfVec fullPhiTestSurfacesBRL(int n = 10, double shift = 0, double zbase = 0,
0087                                 double incl = std::numbers::pi / 9.,
0088                                 double w = 2, double h = 1.5) {
0089     SrfVec res;
0090
0091     double phiStep = 2 * std::numbers::pi / n;
0092     for (int i = 0; i < n; ++i) {
0093       double z = zbase;
0094
0095       Transform3 trans;
0096       trans.setIdentity();
0097       trans.rotate(Eigen::AngleAxisd(i * phiStep + shift, Vector3(0, 0, 1)));
0098       trans.translate(Vector3(10, 0, z));
0099       trans.rotate(Eigen::AngleAxisd(incl, Vector3(0, 0, 1)));
0100       trans.rotate(Eigen::AngleAxisd(std::numbers::pi / 2., Vector3(0, 1, 0)));
0101
0102       auto bounds = std::make_shared<const RectangleBounds>(w, h);
0103       std::shared_ptr<const Surface> srf =
0104           Surface::makeShared<PlaneSurface>(trans, bounds);
0105
0106       res.push_back(srf);
0107       m_surfaces.push_back(
0108           std::move(srf));  // keep shared, will get destroyed at the end
0109     }
0110
0111     return res;
0112   }
0113
0114   SrfVec straightLineSurfaces(
0115       std::size_t n = 10., double step = 3, const Vector3& origin = {0, 0, 1.5},
0116       const Transform3& pretrans = Transform3::Identity(),
0117       const Vector3& dir = {0, 0, 1}) {
0118     SrfVec res;
0119     for (std::size_t i = 0; i < n; ++i) {
0120       Transform3 trans;
0121       trans.setIdentity();
0122       trans.translate(origin + dir * step * i);
0123       // trans.rotate(AngleAxis3(std::numbers::pi/9., Vector3(0, 0, 1)));
0124       trans.rotate(AngleAxis3(std::numbers::pi / 2., Vector3(1, 0, 0)));
0125       trans = trans * pretrans;
0126
0127       auto bounds = std::make_shared<const RectangleBounds>(2, 1.5);
0128
0129       std::shared_ptr<const Surface> srf =
0130           Surface::makeShared<PlaneSurface>(trans, bounds);
0131
0132       res.push_back(srf);
0133       m_surfaces.push_back(
0134           std::move(srf));  // keep shared, will get destroyed at the end
0135     }
0136
0137     return res;
0138   }
0139
0140   SrfVec makeBarrel(int nPhi, int nZ, double w, double h) {
0141     double z0 = -(nZ - 1) * w;
0142     SrfVec res;
0143
0144     for (int i = 0; i < nZ; i++) {
0145       double z = i * w * 2 + z0;
0146       SrfVec ring =
0147           fullPhiTestSurfacesBRL(nPhi, 0, z, std::numbers::pi / 9., w, h);
0148       res.insert(res.end(), ring.begin(), ring.end());
0149     }
0150
0151     return res;
0152   }
0153
0154   void draw_surfaces(const SrfVec& surfaces, const std::string& fname) {
0155     std::ofstream os;
0156     os.open(fname);
0157
0158     os << std::fixed << std::setprecision(4);
0159
0160     std::size_t nVtx = 0;
0161     for (const auto& srfx : surfaces) {
0162       std::shared_ptr<const PlaneSurface> srf =
0163           std::dynamic_pointer_cast<const PlaneSurface>(srfx);
0164       const PlanarBounds* bounds =
0165           dynamic_cast<const PlanarBounds*>(&srf->bounds());
0166
0167       for (const auto& vtxloc : bounds->vertices()) {
0168         Vector3 vtx = srf->localToGlobalTransform(tgContext) *
0169                       Vector3(vtxloc.x(), vtxloc.y(), 0);
0170         os << "v " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
0171       }
0172
0173       // connect them
0174       os << "f";
0175       for (std::size_t i = 1; i <= bounds->vertices().size(); ++i) {
0176         os << " " << nVtx + i;
0177       }
0178       os << "\n";
0179
0180       nVtx += bounds->vertices().size();
0181     }
0182
0183     os.close();
0184   }
0185 };
0186
0187 BOOST_AUTO_TEST_SUITE(SurfacesSuite)
0188
0189 BOOST_FIXTURE_TEST_CASE(SurfaceArray_create, SurfaceArrayFixture) {
0190   GeometryContext tgContext = GeometryContext::dangerouslyDefaultConstruct();
0191
0192   SrfVec brl = makeBarrel(30, 7, 2, 1);
0193   std::vector<const Surface*> brlRaw = unpackSmartPointers(brl);
0194   draw_surfaces(brl, "SurfaceArray_create_BRL_1.obj");
0195
0196   Axis<AxisType::Equidistant, AxisBoundaryType::Closed> phiAxis(
0197       -std::numbers::pi, std::numbers::pi, 30u);
0198   Axis<AxisType::Equidistant, AxisBoundaryType::Bound> zAxis(-14, 14, 7u);
0199
0200   double R = 10;
0201   auto itransform = [R](const Vector2& loc) {
0202     return Vector3(R * std::cos(loc[0]), R * std::sin(loc[0]), loc[1]);
0203   };
0204
0205   auto cylinder =
0206       Surface::makeShared<CylinderSurface>(Transform3::Identity(), R, 10);
0207   SurfaceArray sa(tgContext, brl, cylinder, 1., std::tuple{phiAxis, zAxis});
0208
0209   // let's see if we can access all surfaces
0210   sa.toStream(tgContext, std::cout);
0211
0212   // the bin a surface's own centre falls into has to hold it
0213   for (const auto& srf : brl) {
0214     const Vector3 ctr = srf->referencePosition(tgContext, AxisDirection::AxisR);
0215     const Vector3 normal = srf->normal(tgContext, ctr, Vector3::UnitZ());
0216     const auto binContent = sa.at(tgContext, ctr, normal);
0217
0218     BOOST_CHECK(std::ranges::find(binContent, srf.get()) != binContent.end());
0219   }
0220
0221   const Vector3 crossing = itransform(Vector2(0, 0));
0222   const auto axes = sa.getAxes();
0223   const std::array<std::size_t, 2> crossingBins{axes.at(0)->getBin(0.),
0224                                                 axes.at(1)->getBin(0.)};
0225
0226   // at normal incidence the track does not slide, so the lookup is the bin
0227   BOOST_CHECK(std::ranges::equal(
0228       sa.neighbors(tgContext, crossing, crossing.normalized()),
0229       sa.neighbors(crossingBins, {0, 0})));
0230
0231   // an inclined track slides along z only, so only z widens. the z bins are 4
0232   // wide against a tolerance of 1: a slope of 3 reaches one bin, 7 two.
0233   for (const auto& [slope, distance] :
0234        std::vector<std::pair<double, std::uint8_t>>{{3., 1}, {7., 2}}) {
0235     const Vector3 direction = Vector3(1, 0, slope).normalized();
0236     BOOST_CHECK(std::ranges::equal(sa.neighbors(tgContext, crossing, direction),
0237                                    sa.neighbors(crossingBins, {0, distance})));
0238     BOOST_CHECK(
0239         !std::ranges::equal(sa.neighbors(tgContext, crossing, direction),
0240                             sa.neighbors(crossingBins, {distance, distance})));
0241   }
0242
0243   // the floor is served regardless of the crossing angle
0244   SurfaceArray floored(tgContext, brl, cylinder, 1., std::tuple{phiAxis, zAxis},
0245                        {{1, 1}, {1, 2}});
0246   BOOST_CHECK(std::ranges::equal(
0247       floored.neighbors(tgContext, crossing, crossing.normalized()),
0248       floored.neighbors(crossingBins, {1, 1})));
0249   BOOST_CHECK_THROW(SurfaceArray(tgContext, brl, cylinder, 1.,
0250                                  std::tuple{phiAxis, zAxis}, {{2, 2}, {1, 2}}),
0251                     std::invalid_argument);
0252
0253   // nothing is cached past the bound, so asking for it is an error
0254   BOOST_CHECK_THROW(floored.neighbors(crossingBins, {2, 0}), std::out_of_range);
0255   BOOST_CHECK_THROW(floored.neighbors(crossingBins, {0, 3}), std::out_of_range);
0256   BOOST_CHECK_THROW(floored.neighbors({10000, 0}, {0, 0}), std::out_of_range);
0257
0258   // a scalar bound is the isotropic window it used to describe
0259   ACTS_PUSH_IGNORE_DEPRECATED()
0260   const SurfaceArray scalarBound(tgContext, brl, cylinder, 1.,
0261                                  std::tuple{phiAxis, zAxis}, std::uint8_t{1});
0262   BOOST_CHECK_EQUAL(scalarBound.maxNeighborDistance(), 1u);
0263   ACTS_POP_IGNORE_DEPRECATED()
0264   const SurfaceArray::NeighborWindow scalarWindow =
0265       scalarBound.neighborWindow();
0266   BOOST_CHECK((scalarWindow.min == std::array<std::uint8_t, 2>{1, 1}));
0267   BOOST_CHECK((scalarWindow.max == std::array<std::uint8_t, 2>{1, 1}));
0268   BOOST_CHECK(std::ranges::equal(
0269       scalarBound.neighbors(tgContext, crossing, crossing.normalized()),
0270       scalarBound.neighbors(crossingBins, {1, 1})));
0271 }
0272
0273 BOOST_AUTO_TEST_CASE(SurfaceArray_overfill) {
0274   const auto representative = Surface::makeShared<PlaneSurface>(
0275       Transform3::Identity(), std::make_shared<RectangleBounds>(3.5, 3.5));
0276   // The footprint covers three cells, so expansion has overlapping
0277   // contributions.
0278   const auto module = Surface::makeShared<PlaneSurface>(
0279       Transform3::Identity(), std::make_shared<RectangleBounds>(0.6, 0.2));
0280   const auto check = [&]<AxisBoundaryType boundary>() {
0281     const Axis<AxisType::Equidistant, boundary> axis(-3.5, 3.5, 7);
0282     for (const std::uint8_t radius : {0, 1, 2, 3, 255}) {
0283       const SurfaceArray array(tgContext, {module}, representative, 0.,
0284                                {axis, axis}, {{0, 0}, {0, 0}}, radius);
0285       for (int x = -3; x <= 3; ++x) {
0286         for (int y = -3; y <= 3; ++y) {
0287           const auto content =
0288               array.at(tgContext, Vector3(x, y, 0), Vector3::UnitZ());
0289           const bool expected =
0290               std::abs(x) <= radius + 1 && std::abs(y) <= radius;
0291           BOOST_CHECK_EQUAL(content.size(), expected ? 1u : 0u);
0292           if (expected) {
0293             BOOST_CHECK_EQUAL(content.front(), module.get());
0294           }
0295         }
0296       }
0297       for (std::size_t bin = 0; bin < array.size(); ++bin) {
0298         if (!array.isValidBin(bin)) {
0299           BOOST_CHECK(array.at(bin).empty());
0300         }
0301       }
0302       BOOST_CHECK(!array.isValidBin(array.size() + 10));
0303     }
0304   };
0305   check.template operator()<AxisBoundaryType::Bound>();
0306   check.template operator()<AxisBoundaryType::Open>();
0307 }
0308
0309 BOOST_FIXTURE_TEST_CASE(SurfaceArray_overfillPhiSeam, SurfaceArrayFixture) {
0310   const auto modules = makeBarrel(12, 3, 2, 0.2);
0311   const auto cylinder =
0312       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 10., 10.);
0313   const Axis<AxisType::Variable, AxisBoundaryType::Closed> phiAxis(
0314       {-std::numbers::pi, -3., -2., -1., 0., 1., 2., 3., std::numbers::pi});
0315   const Axis<AxisType::Equidistant, AxisBoundaryType::Bound> zAxis(-6., 6., 3);
0316   const SurfaceArray nominal(tgContext, modules, cylinder, 1., {phiAxis, zAxis},
0317                              {{0, 0}, {3, 3}});
0318   for (const std::uint8_t radius : {0, 1, 2, 3}) {
0319     const SurfaceArray expanded(tgContext, modules, cylinder, 1.,
0320                                 {phiAxis, zAxis}, {{0, 0}, {0, 0}}, radius);
0321     for (std::size_t phiBin = 1; phiBin <= phiAxis.getNBins(); ++phiBin) {
0322       for (std::size_t zBin = 1; zBin <= zAxis.getNBins(); ++zBin) {
0323         BOOST_CHECK(
0324             std::ranges::equal(expanded.neighbors({phiBin, zBin}, 0),
0325                                nominal.neighbors({phiBin, zBin}, radius)));
0326       }
0327     }
0328   }
0329 }
0330
0331 BOOST_AUTO_TEST_CASE(SurfaceArray_maximumWindow) {
0332   const auto plane = Surface::makeShared<PlaneSurface>(
0333       Transform3::Identity(), std::make_shared<RectangleBounds>(1., 1.));
0334   const Axis<AxisType::Equidistant, AxisBoundaryType::Bound> axis(-1., 1., 1);
0335   // Both counters must terminate at the largest representable distance.
0336   const SurfaceArray array(tgContext, {plane}, plane, 0., {axis, axis},
0337                            {{0, 0}, {255, 255}});
0338   BOOST_CHECK_EQUAL(array.neighbors({1, 1}, {255, 255}).size(), 1u);
0339 }
0340
0341 BOOST_FIXTURE_TEST_CASE(SurfaceArray_periodicWindow, SurfaceArrayFixture) {
0342   const auto modules = makeBarrel(30, 1, 2, 0.2);
0343   const auto cylinder =
0344       Surface::makeShared<CylinderSurface>(Transform3::Identity(), 10., 10.);
0345   const Axis<AxisType::Equidistant, AxisBoundaryType::Closed> phiAxis(
0346       -std::numbers::pi, std::numbers::pi, 30);
0347   const Axis<AxisType::Equidistant, AxisBoundaryType::Bound> zAxis(-3., 3., 1);
0348   const SurfaceArray array(tgContext, modules, cylinder, 1., {phiAxis, zAxis},
0349                            {{0, 0}, {2, 0}});
0350   const double phi = 0.1;
0351   const Vector3 crossing(10. * std::cos(phi), 10. * std::sin(phi), 0.);
0352   for (const double slide : {0., 2. * std::numbers::pi - 0.01,
0353                              2. * std::numbers::pi, 4. * std::numbers::pi}) {
0354     const Vector3 direction =
0355         Vector3(std::cos(phi) - 10. * slide * std::sin(phi),
0356                 std::sin(phi) + 10. * slide * std::cos(phi), 0.)
0357             .normalized();
0358     const auto expected =
0359         array.neighbors({phiAxis.getBin(phi), 1},
0360                         {static_cast<std::uint8_t>(slide == 0. ? 0 : 2), 0});
0361     BOOST_CHECK(std::ranges::equal(
0362         array.neighbors(tgContext, crossing, direction), expected));
0363     BOOST_CHECK(std::ranges::equal(
0364         array.neighbors(tgContext, crossing, -direction), expected));
0365   }
0366 }
0367
0368 BOOST_FIXTURE_TEST_CASE(SurfaceArray_variablePeriodicWindow,
0369                         SurfaceArrayFixture) {
0370   const auto modules = fullPhiTestSurfacesEC(30);
0371   const auto disc =
0372       Surface::makeShared<DiscSurface>(Transform3::Identity(), 0., 20.);
0373   const Axis<AxisType::Equidistant, AxisBoundaryType::Bound> rAxis(0., 20., 1);
0374   const Axis<AxisType::Variable, AxisBoundaryType::Closed> phiAxis(
0375       {-std::numbers::pi, -3., -2.9, -2.8, -2.7, -1., 0., 1., 2., 3.,
0376        std::numbers::pi});
0377   const SurfaceArray array(tgContext, modules, disc, 1., {rAxis, phiAxis},
0378                            {{0, 0}, {0, 4}});
0379   const double phi = 3.05;
0380   const Vector3 crossing(10. * std::cos(phi), 10. * std::sin(phi), 0.);
0381   const Vector3 direction =
0382       Vector3(-8. * std::sin(phi), 8. * std::cos(phi), 1.).normalized();
0383   BOOST_CHECK(
0384       std::ranges::equal(array.neighbors(tgContext, crossing, direction),
0385                          array.neighbors({1, phiAxis.getBin(phi)}, {0, 4})));
0386   // At the polar singularity, an undefined derivative must not enter getBin.
0387   std::feclearexcept(FE_INVALID | FE_DIVBYZERO);
0388   BOOST_CHECK(std::ranges::equal(
0389       array.neighbors(tgContext, Vector3::Zero(), Vector3::UnitZ()),
0390       array.neighbors({1, phiAxis.getBin(0.)}, {0, 4})));
0391   BOOST_CHECK(
0392       std::ranges::equal(array.neighbors(tgContext, Vector3::Zero(), direction),
0393                          array.neighbors({1, phiAxis.getBin(0.)}, {0, 4})));
0394   BOOST_CHECK_EQUAL(std::fetestexcept(FE_INVALID | FE_DIVBYZERO), 0);
0395 }
0396
0397 BOOST_AUTO_TEST_CASE(SurfaceArray_singleElement) {
0398   const double w = 3;
0399   const double h = 4;
0400   const auto bounds = std::make_shared<const RectangleBounds>(w, h);
0401   auto srf = Surface::makeShared<PlaneSurface>(Transform3::Identity(), bounds);
0402
0403   SurfaceArray sa(srf);
0404
0405   const auto binContent =
0406       sa.at(tgContext, Vector3(42, 42, 42), Vector3::UnitX());
0407   BOOST_CHECK_EQUAL(binContent.size(), 1u);
0408   BOOST_CHECK_EQUAL(binContent[0], srf.get());
0409   BOOST_CHECK_EQUAL(sa.surfaces().size(), 1u);
0410   BOOST_CHECK_EQUAL(sa.surfaces().at(0), srf.get());
0411 }
0412
0413 BOOST_AUTO_TEST_CASE(SurfaceArrayToStreamPreservesStreamState) {
0414   const auto bounds = std::make_shared<const RectangleBounds>(3., 4.);
0415   auto surface =
0416       Surface::makeShared<PlaneSurface>(Transform3::Identity(), bounds);
0417   SurfaceArray surfaceArray(surface);
0418
0419   std::ostringstream stream;
0420   stream << std::scientific << std::showpos << std::setfill('#')
0421          << std::setprecision(3);
0422   stream.width(17);
0423
0424   const auto flags = stream.flags();
0425   const auto precision = stream.precision();
0426   const auto width = stream.width();
0427   const auto fill = stream.fill();
0428
0429   surfaceArray.toStream(tgContext, stream);
0430
0431   BOOST_CHECK(stream.flags() == flags);
0432   BOOST_CHECK_EQUAL(stream.precision(), precision);
0433   BOOST_CHECK_EQUAL(stream.width(), width);
0434   BOOST_CHECK_EQUAL(stream.fill(), fill);
0435 }
0436
0437 BOOST_AUTO_TEST_SUITE_END()
0438
0439 }  // namespace ActsTests