Core/Geometry/SurfaceArrayCreatorTests.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/Geometry/ProtoLayer.hpp"
0014 #include "Acts/Geometry/SurfaceArrayCreator.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/Surfaces/SurfaceBounds.hpp"
0021 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0022 #include "Acts/Utilities/Axis.hpp"
0023 #include "Acts/Utilities/AxisDefinitions.hpp"
0024 #include "Acts/Utilities/BinningType.hpp"
0025 #include "Acts/Utilities/Grid.hpp"
0026 #include "Acts/Utilities/Helpers.hpp"
0027 #include "Acts/Utilities/IAxis.hpp"
0028 #include "Acts/Utilities/Logger.hpp"
0029 #include "Acts/Utilities/detail/grid_helper.hpp"
0030 #include "Acts/Visualization/GeometryView3D.hpp"
0031 #include "Acts/Visualization/ObjVisualization3D.hpp"
0032
0033 #include <algorithm>
0034 #include <cmath>
0035 #include <cstddef>
0036 #include <fstream>
0037 #include <iomanip>
0038 #include <iterator>
0039 #include <memory>
0040 #include <numbers>
0041 #include <set>
0042 #include <stdexcept>
0043 #include <string>
0044 #include <tuple>
0045 #include <utility>
0046 #include <vector>
0047
0048 #include <boost/format.hpp>
0049
0050 using Acts::VectorHelpers::perp;
0051 using Acts::VectorHelpers::phi;
0052
0053 namespace Acts::Test {
0054
0055 // Create a test context
0056 GeometryContext tgContext = GeometryContext();
0057
0058 using SrfVec = std::vector<std::shared_ptr<const Surface>>;
0059
0060 struct SurfaceArrayCreatorFixture {
0061   SurfaceArrayCreator m_SAC;
0062   std::vector<std::shared_ptr<const Surface>> m_surfaces;
0063
0064   SurfaceArrayCreatorFixture()
0065       : m_SAC(SurfaceArrayCreator::Config(),
0066               Acts::getDefaultLogger("SurfaceArrayCreator",
0067                                      Acts::Logging::VERBOSE)) {
0068     BOOST_TEST_MESSAGE("setup fixture");
0069   }
0070   ~SurfaceArrayCreatorFixture() { BOOST_TEST_MESSAGE("teardown fixture"); }
0071
0072   template <typename... Args>
0073   SurfaceArrayCreator::ProtoAxis createEquidistantAxis(Args&&... args) {
0074     return m_SAC.createEquidistantAxis(std::forward<Args>(args)...);
0075   }
0076
0077   template <typename... Args>
0078   SurfaceArrayCreator::ProtoAxis createVariableAxis(Args&&... args) {
0079     return m_SAC.createVariableAxis(std::forward<Args>(args)...);
0080   }
0081
0082   template <AxisBoundaryType bdtA, AxisBoundaryType bdtB, typename... Args>
0083   std::unique_ptr<SurfaceArray::ISurfaceGridLookup> makeSurfaceGridLookup2D(
0084       Args&&... args) {
0085     return m_SAC.makeSurfaceGridLookup2D<bdtA, bdtB>(
0086         std::forward<Args>(args)...);
0087   }
0088
0089   SrfVec fullPhiTestSurfacesEC(std::size_t n = 10, double shift = 0,
0090                                double zbase = 0, double r = 10, double w = 2,
0091                                double h = 1) {
0092     SrfVec res;
0093     // TODO: The test is extremely numerically unstable in the face of upward
0094     //       rounding in this multiplication and division. Find out why.
0095     double phiStep = 2 * std::numbers::pi / n;
0096     for (std::size_t i = 0; i < n; ++i) {
0097       double z = zbase + ((i % 2 == 0) ? 1 : -1) * 0.2;
0098       double phi = std::fma(i, phiStep, shift);
0099
0100       Transform3 trans;
0101       trans.setIdentity();
0102       trans.rotate(Eigen::AngleAxisd(phi, Vector3(0, 0, 1)));
0103       trans.translate(Vector3(r, 0, z));
0104
0105       auto bounds = std::make_shared<const RectangleBounds>(w, h);
0106
0107       std::shared_ptr<Surface> srf =
0108           Surface::makeShared<PlaneSurface>(trans, bounds);
0109
0110       res.push_back(srf);
0111       m_surfaces.push_back(
0112           std::move(srf));  // keep shared, will get destroyed at the end
0113     }
0114
0115     return res;
0116   }
0117
0118   SrfVec fullPhiTestSurfacesBRL(std::size_t n = 10, double shift = 0,
0119                                 double zbase = 0,
0120                                 double incl = std::numbers::pi / 9.,
0121                                 double w = 2, double h = 1.5) {
0122     SrfVec res;
0123     // TODO: The test is extremely numerically unstable in the face of upward
0124     //       rounding in this multiplication and division. Find out why.
0125     double phiStep = 2 * std::numbers::pi / n;
0126     for (std::size_t i = 0; i < n; ++i) {
0127       double z = zbase;
0128       double phi = std::fma(i, phiStep, shift);
0129
0130       Transform3 trans;
0131       trans.setIdentity();
0132       trans.rotate(Eigen::AngleAxisd(phi, Vector3(0, 0, 1)));
0133       trans.translate(Vector3(10, 0, z));
0134       trans.rotate(Eigen::AngleAxisd(incl, Vector3(0, 0, 1)));
0135       trans.rotate(Eigen::AngleAxisd(std::numbers::pi / 2., Vector3(0, 1, 0)));
0136
0137       auto bounds = std::make_shared<const RectangleBounds>(w, h);
0138       std::shared_ptr<Surface> srf =
0139           Surface::makeShared<PlaneSurface>(trans, bounds);
0140
0141       res.push_back(srf);
0142       m_surfaces.push_back(
0143           std::move(srf));  // keep shared, will get destroyed at the end
0144     }
0145
0146     return res;
0147   }
0148
0149   SrfVec straightLineSurfaces(
0150       std::size_t n = 10., double step = 3, const Vector3& origin = {0, 0, 1.5},
0151       const Transform3& pretrans = Transform3::Identity(),
0152       const Vector3& dir = {0, 0, 1}) {
0153     SrfVec res;
0154     for (std::size_t i = 0; i < n; ++i) {
0155       Transform3 trans;
0156       trans.setIdentity();
0157       trans.translate(origin + dir * step * i);
0158       // trans.rotate(AngleAxis3(std::numbers::pi / 9., Vector3(0, 0, 1)));
0159       trans.rotate(AngleAxis3(std::numbers::pi / 2., Vector3(1, 0, 0)));
0160       trans = trans * pretrans;
0161
0162       auto bounds = std::make_shared<const RectangleBounds>(2, 1.5);
0163
0164       std::shared_ptr<Surface> srf =
0165           Surface::makeShared<PlaneSurface>(trans, bounds);
0166
0167       res.push_back(srf);
0168       m_surfaces.push_back(
0169           std::move(srf));  // keep shared, will get destroyed at the end
0170     }
0171
0172     return res;
0173   }
0174
0175   SrfVec makeBarrel(int nPhi, int nZ, double w, double h) {
0176     double z0 = -(nZ - 1) * w;
0177     SrfVec res;
0178
0179     for (int i = 0; i < nZ; i++) {
0180       double z = i * w * 2 + z0;
0181       // std::cout << "z=" << z << std::endl;
0182       SrfVec ring =
0183           fullPhiTestSurfacesBRL(nPhi, 0, z, std::numbers::pi / 9., w, h);
0184       res.insert(res.end(), ring.begin(), ring.end());
0185     }
0186
0187     return res;
0188   }
0189
0190   std::pair<SrfVec, std::vector<std::pair<const Surface*, const Surface*>>>
0191   makeBarrelStagger(int nPhi, int nZ, double shift = 0,
0192                     double incl = std::numbers::pi / 9., double w = 2,
0193                     double h = 1.5) {
0194     double z0 = -(nZ - 1) * w;
0195     SrfVec res;
0196     std::vector<std::pair<const Surface*, const Surface*>> pairs;
0197     // TODO: The test is extremely numerically unstable in the face of upward
0198     //       rounding in this multiplication and division. Find out why.
0199     double phiStep = 2 * std::numbers::pi / nPhi;
0200     for (int i = 0; i < nZ; i++) {
0201       double z = i * w * 2 + z0;
0202       for (int j = 0; j < nPhi; ++j) {
0203         double phi = std::fma(j, phiStep, shift);
0204         Transform3 trans;
0205         trans.setIdentity();
0206         trans.rotate(Eigen::AngleAxisd(phi, Vector3(0, 0, 1)));
0207         trans.translate(Vector3(10, 0, z));
0208         trans.rotate(Eigen::AngleAxisd(incl, Vector3(0, 0, 1)));
0209         trans.rotate(
0210             Eigen::AngleAxisd(std::numbers::pi / 2., Vector3(0, 1, 0)));
0211
0212         auto bounds = std::make_shared<const RectangleBounds>(w, h);
0213         std::shared_ptr<PlaneSurface> srfA =
0214             Surface::makeShared<PlaneSurface>(trans, bounds);
0215
0216         Vector3 nrm = srfA->normal(tgContext);
0217         Transform3 transB = trans;
0218         transB.pretranslate(nrm * 0.1);
0219         std::shared_ptr<Surface> srfB =
0220             Surface::makeShared<PlaneSurface>(transB, bounds);
0221
0222         pairs.push_back(std::make_pair(srfA.get(), srfB.get()));
0223
0224         res.push_back(srfA);
0225         res.push_back(srfB);
0226         m_surfaces.push_back(std::move(srfA));
0227         m_surfaces.push_back(std::move(srfB));
0228       }
0229     }
0230
0231     return {res, pairs};
0232   }
0233 };
0234
0235 void draw_surfaces(const SrfVec& surfaces, const std::string& fname) {
0236   std::ofstream os;
0237   os.open(fname);
0238
0239   os << std::fixed << std::setprecision(4);
0240
0241   std::size_t nVtx = 0;
0242   for (const auto& srfx : surfaces) {
0243     std::shared_ptr<const PlaneSurface> srf =
0244         std::dynamic_pointer_cast<const PlaneSurface>(srfx);
0245     const PlanarBounds* bounds =
0246         dynamic_cast<const PlanarBounds*>(&srf->bounds());
0247
0248     for (const auto& vtxloc : bounds->vertices()) {
0249       Vector3 vtx =
0250           srf->transform(tgContext) * Vector3(vtxloc.x(), vtxloc.y(), 0);
0251       os << "v " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
0252     }
0253
0254     // connect them
0255     os << "f";
0256     for (std::size_t i = 1; i <= bounds->vertices().size(); ++i) {
0257       os << " " << nVtx + i;
0258     }
0259     os << "\n";
0260
0261     nVtx += bounds->vertices().size();
0262   }
0263
0264   os.close();
0265 }
0266
0267 BOOST_AUTO_TEST_SUITE(Tools)
0268
0269 BOOST_FIXTURE_TEST_CASE(SurfaceArrayCreator_createEquidistantAxis_Phi,
0270                         SurfaceArrayCreatorFixture) {
0271   // fail on empty srf vector
0272   std::vector<const Surface*> emptyRaw;
0273   ProtoLayer pl(tgContext, emptyRaw);
0274   auto tr = Transform3::Identity();
0275   BOOST_CHECK_THROW(createEquidistantAxis(tgContext, emptyRaw,
0276                                           AxisDirection::AxisPhi, pl, tr),
0277                     std::logic_error);
0278
0279   std::vector<float> bdExp = {
0280       -3.14159, -2.93215, -2.72271, -2.51327, -2.30383,  -2.0944,   -1.88496,
0281       -1.67552, -1.46608, -1.25664, -1.0472,  -0.837758, -0.628319, -0.418879,
0282       -0.20944, 0,        0.20944,  0.418879, 0.628319,  0.837758,  1.0472,
0283       1.25664,  1.46608,  1.67552,  1.88496,  2.09439,   2.30383,   2.51327,
0284       2.72271,  2.93215,  3.14159};
0285
0286   double step = 2 * std::numbers::pi / 30.;
0287
0288   // endcap style modules
0289
0290   for (int i = -1; i <= 2; i += 2) {
0291     double z = 10 * i;
0292     // case 1: one module sits at pi / -pi
0293     double angleShift = step / 2.;
0294     auto surfaces = fullPhiTestSurfacesEC(30, angleShift, z);
0295     std::vector<const Surface*> surfacesRaw = unpack_shared_vector(surfaces);
0296     pl = ProtoLayer(tgContext, surfacesRaw);
0297     tr = Transform3::Identity();
0298     auto axis = createEquidistantAxis(tgContext, surfacesRaw,
0299                                       AxisDirection::AxisPhi, pl, tr);
0300
0301     BOOST_CHECK_EQUAL(axis.nBins, 30u);
0302     CHECK_CLOSE_REL(axis.max, std::numbers::pi, 1e-6);
0303     CHECK_CLOSE_REL(axis.min, -std::numbers::pi, 1e-6);
0304     BOOST_CHECK_EQUAL(axis.bType, equidistant);
0305     CHECK_SMALL(phi(tr * Vector3::UnitX()), 1e-6);
0306
0307     // case 2: two modules sit symmetrically around pi / -pi
0308     angleShift = 0.;
0309     surfaces = fullPhiTestSurfacesEC(30, angleShift, z);
0310     surfacesRaw = unpack_shared_vector(surfaces);
0311     pl = ProtoLayer(tgContext, surfacesRaw);
0312     tr = Transform3::Identity();
0313     axis = createEquidistantAxis(tgContext, surfacesRaw, AxisDirection::AxisPhi,
0314                                  pl, tr);
0315     draw_surfaces(surfaces,
0316                   "SurfaceArrayCreator_createEquidistantAxis_EC_2.obj");
0317     BOOST_CHECK_EQUAL(axis.nBins, 30u);
0318     CHECK_CLOSE_REL(axis.max, std::numbers::pi, 1e-6);
0319     CHECK_CLOSE_REL(axis.min, -std::numbers::pi, 1e-6);
0320     BOOST_CHECK_EQUAL(axis.bType, equidistant);
0321     // CHECK_CLOSE_REL(bdExp, axis.binEdges, 0.001);
0322     CHECK_CLOSE_REL(phi(tr * Vector3::UnitX()), -0.5 * step, 1e-3);
0323     // case 3: two modules sit asymmetrically around pi / -pi shifted up
0324     angleShift = step / -4.;
0325     surfaces = fullPhiTestSurfacesEC(30, angleShift, z);
0326     surfacesRaw = unpack_shared_vector(surfaces);
0327     pl = ProtoLayer(tgContext, surfacesRaw);
0328     tr = Transform3::Identity();
0329     axis = createEquidistantAxis(tgContext, surfacesRaw, AxisDirection::AxisPhi,
0330                                  pl, tr);
0331     draw_surfaces(surfaces,
0332                   "SurfaceArrayCreator_createEquidistantAxis_EC_3.obj");
0333     BOOST_CHECK_EQUAL(axis.nBins, 30u);
0334     CHECK_CLOSE_REL(axis.max, std::numbers::pi, 1e-6);
0335     CHECK_CLOSE_REL(axis.min, -std::numbers::pi, 1e-6);
0336     BOOST_CHECK_EQUAL(axis.bType, equidistant);
0337     CHECK_CLOSE_REL(phi(tr * Vector3::UnitX()), step / -4., 1e-3);
0338
0339     // case 4: two modules sit asymmetrically around pi / -pi shifted down
0340     angleShift = step / 4.;
0341     surfaces = fullPhiTestSurfacesEC(30, angleShift, z);
0342     surfacesRaw = unpack_shared_vector(surfaces);
0343     pl = ProtoLayer(tgContext, surfaces);
0344     surfacesRaw = unpack_shared_vector(surfaces);
0345     tr = Transform3::Identity();
0346     axis = createEquidistantAxis(tgContext, surfacesRaw, AxisDirection::AxisPhi,
0347                                  pl, tr);
0348     surfacesRaw = unpack_shared_vector(surfaces);
0349     draw_surfaces(surfaces,
0350                   "SurfaceArrayCreator_createEquidistantAxis_EC_4.obj");
0351     BOOST_CHECK_EQUAL(axis.nBins, 30u);
0352     CHECK_CLOSE_REL(axis.max, std::numbers::pi, 1e-6);
0353     CHECK_CLOSE_REL(axis.min, -std::numbers::pi, 1e-6);
0354     BOOST_CHECK_EQUAL(axis.bType, equidistant);
0355     CHECK_CLOSE_REL(phi(tr * Vector3::UnitX()), step / 4., 1e-3);
0356   }
0357
0358   for (int i = -1; i <= 2; i += 2) {
0359     double z = 10 * i;
0360     // case 1: one module sits at pi / -pi
0361     double angleShift = step / 2.;
0362     auto surfaces = fullPhiTestSurfacesBRL(30, angleShift, z);
0363     auto surfacesRaw = unpack_shared_vector(surfaces);
0364     pl = ProtoLayer(tgContext, surfacesRaw);
0365     tr = Transform3::Identity();
0366     auto axis = createEquidistantAxis(tgContext, surfacesRaw,
0367                                       AxisDirection::AxisPhi, pl, tr);
0368     draw_surfaces(surfaces,
0369                   "SurfaceArrayCreator_createEquidistantAxis_BRL_1.obj");
0370     BOOST_CHECK_EQUAL(axis.nBins, 30u);
0371     CHECK_CLOSE_REL(axis.max, std::numbers::pi, 1e-6);
0372     CHECK_CLOSE_REL(axis.min, -std::numbers::pi, 1e-6);
0373     BOOST_CHECK_EQUAL(axis.bType, equidistant);
0374     CHECK_SMALL(phi(tr * Vector3::UnitX()), 1e-6);
0375
0376     // case 2: two modules sit symmetrically around pi / -pi
0377     angleShift = 0.;
0378     surfaces = fullPhiTestSurfacesBRL(30, angleShift, z);
0379     surfacesRaw = unpack_shared_vector(surfaces);
0380     pl = ProtoLayer(tgContext, surfacesRaw);
0381     tr = Transform3::Identity();
0382     axis = createEquidistantAxis(tgContext, surfacesRaw, AxisDirection::AxisPhi,
0383                                  pl, tr);
0384     draw_surfaces(surfaces,
0385                   "SurfaceArrayCreator_createEquidistantAxis_BRL_2.obj");
0386     BOOST_CHECK_EQUAL(axis.nBins, 30u);
0387     CHECK_CLOSE_REL(axis.max, std::numbers::pi, 1e-6);
0388     CHECK_CLOSE_REL(axis.min, -std::numbers::pi, 1e-6);
0389     BOOST_CHECK_EQUAL(axis.bType, equidistant);
0390     // CHECK_CLOSE_REL(bdExp, axis.binEdges, 0.001);
0391     CHECK_CLOSE_REL(phi(tr * Vector3::UnitX()), -0.5 * step, 1e-3);
0392
0393     // case 3: two modules sit asymmetrically around pi / -pi shifted up
0394     angleShift = step / -4.;
0395     surfaces = fullPhiTestSurfacesBRL(30, angleShift, z);
0396     surfacesRaw = unpack_shared_vector(surfaces);
0397     pl = ProtoLayer(tgContext, surfacesRaw);
0398     tr = Transform3::Identity();
0399     axis = createEquidistantAxis(tgContext, surfacesRaw, AxisDirection::AxisPhi,
0400                                  pl, tr);
0401     draw_surfaces(surfaces,
0402                   "SurfaceArrayCreator_createEquidistantAxis_BRL_3.obj");
0403     BOOST_CHECK_EQUAL(axis.nBins, 30u);
0404     CHECK_CLOSE_REL(axis.max, std::numbers::pi, 1e-6);
0405     CHECK_CLOSE_REL(axis.min, -std::numbers::pi, 1e-6);
0406     BOOST_CHECK_EQUAL(axis.bType, equidistant);
0407     // CHECK_CLOSE_REL(bdExp, axis.binEdges, 0.001);
0408     CHECK_CLOSE_REL(phi(tr * Vector3::UnitX()), step / -4., 1e-3);
0409
0410     // case 4: two modules sit asymmetrically around pi / -pi shifted down
0411     angleShift = step / 4.;
0412     surfaces = fullPhiTestSurfacesBRL(30, angleShift, z);
0413     surfacesRaw = unpack_shared_vector(surfaces);
0414     pl = ProtoLayer(tgContext, surfacesRaw);
0415     tr = Transform3::Identity();
0416     axis = createEquidistantAxis(tgContext, surfacesRaw, AxisDirection::AxisPhi,
0417                                  pl, tr);
0418     draw_surfaces(surfaces,
0419                   "SurfaceArrayCreator_createEquidistantAxis_BRL_4.obj");
0420     BOOST_CHECK_EQUAL(axis.nBins, 30u);
0421     CHECK_CLOSE_REL(axis.max, std::numbers::pi, 1e-6);
0422     CHECK_CLOSE_REL(axis.min, -std::numbers::pi, 1e-6);
0423     BOOST_CHECK_EQUAL(axis.bType, equidistant);
0424     // CHECK_CLOSE_REL(bdExp, axis.binEdges, 0.001);
0425     CHECK_CLOSE_REL(phi(tr * Vector3::UnitX()), step / 4., 1e-3);
0426   }
0427
0428   SrfVec surfaces;
0429
0430   // single element in phi
0431   surfaces = fullPhiTestSurfacesEC(1);
0432   auto surfacesRaw = unpack_shared_vector(surfaces);
0433   draw_surfaces(surfaces,
0434                 "SurfaceArrayCreator_createEquidistantAxis_EC_Single.obj");
0435
0436   pl = ProtoLayer(tgContext, surfacesRaw);
0437   tr = Transform3::Identity();
0438   auto axis = createEquidistantAxis(tgContext, surfacesRaw,
0439                                     AxisDirection::AxisPhi, pl, tr);
0440   BOOST_CHECK_EQUAL(axis.nBins, 1u);
0441
0442   CHECK_CLOSE_ABS(axis.max, phi(Vector3(8, 1, 0)), 1e-3);
0443   CHECK_CLOSE_ABS(axis.min, phi(Vector3(8, -1, 0)), 1e-3);
0444   BOOST_CHECK_EQUAL(axis.bType, equidistant);
0445 }
0446
0447 BOOST_FIXTURE_TEST_CASE(SurfaceArrayCreator_createEquidistantAxis_Z,
0448                         SurfaceArrayCreatorFixture) {
0449   // single element in z
0450   auto surfaces = straightLineSurfaces(1);
0451   auto surfacesRaw = unpack_shared_vector(surfaces);
0452   ProtoLayer pl = ProtoLayer(tgContext, surfacesRaw);
0453   auto trf = Transform3::Identity();
0454   auto axis = createEquidistantAxis(tgContext, surfacesRaw,
0455                                     AxisDirection::AxisZ, pl, trf);
0456   draw_surfaces(surfaces, "SurfaceArrayCreator_createEquidistantAxis_Z_1.obj");
0457   BOOST_CHECK_EQUAL(axis.nBins, 1u);
0458   CHECK_CLOSE_ABS(axis.max, 3, 1e-6);
0459   CHECK_CLOSE_ABS(axis.min, 0, 1e-6);
0460   BOOST_CHECK_EQUAL(axis.bType, equidistant);
0461
0462   // z rows with varying starting point
0463   for (std::size_t i = 0; i <= 20; i++) {
0464     double z0 = -10 + 1. * i;
0465     surfaces = straightLineSurfaces(10, 3, Vector3(0, 0, z0 + 1.5));
0466     surfacesRaw = unpack_shared_vector(surfaces);
0467     pl = ProtoLayer(tgContext, surfacesRaw);
0468     trf = Transform3::Identity();
0469     axis = createEquidistantAxis(tgContext, surfacesRaw, AxisDirection::AxisZ,
0470                                  pl, trf);
0471     draw_surfaces(
0472         surfaces,
0473         (boost::format(
0474              "SurfaceArrayCreator_createEquidistantAxis_Z_2_%1%.obj") %
0475          i)
0476             .str());
0477     BOOST_CHECK_EQUAL(axis.nBins, 10u);
0478     CHECK_CLOSE_ABS(axis.max, 30 + z0, 1e-6);
0479     CHECK_CLOSE_ABS(axis.min, z0, 1e-6);
0480     BOOST_CHECK_EQUAL(axis.bType, equidistant);
0481   }
0482
0483   // z row where elements are rotated around y
0484   Transform3 tr = Transform3::Identity();
0485   tr.rotate(AngleAxis3(std::numbers::pi / 4., Vector3(0, 0, 1)));
0486   surfaces = straightLineSurfaces(10, 3, Vector3(0, 0, 0 + 1.5), tr);
0487   surfacesRaw = unpack_shared_vector(surfaces);
0488   pl = ProtoLayer(tgContext, surfacesRaw);
0489   trf = Transform3::Identity();
0490   axis = createEquidistantAxis(tgContext, surfacesRaw, AxisDirection::AxisZ, pl,
0491                                trf);
0492   draw_surfaces(surfaces, "SurfaceArrayCreator_createEquidistantAxis_Z_3.obj");
0493   BOOST_CHECK_EQUAL(axis.nBins, 10u);
0494   CHECK_CLOSE_ABS(axis.max, 30.9749, 1e-3);
0495   CHECK_CLOSE_ABS(axis.min, -0.974873, 1e-3);
0496   BOOST_CHECK_EQUAL(axis.bType, equidistant);
0497 }
0498
0499 BOOST_FIXTURE_TEST_CASE(SurfaceArrayCreator_createEquidistantAxis_R,
0500                         SurfaceArrayCreatorFixture) {
0501   // single element in r
0502   auto surfaces = fullPhiTestSurfacesEC(1, 0, 0, 15);
0503   auto surfacesRaw = unpack_shared_vector(surfaces);
0504   draw_surfaces(surfaces, "SurfaceArrayCreator_createEquidistantAxis_R_1.obj");
0505   auto trf = Transform3::Identity();
0506   ProtoLayer pl = ProtoLayer(tgContext, surfacesRaw);
0507   auto axis = createEquidistantAxis(tgContext, surfacesRaw,
0508                                     AxisDirection::AxisR, pl, trf);
0509   BOOST_CHECK_EQUAL(axis.nBins, 1u);
0510   CHECK_CLOSE_ABS(axis.max, perp(Vector3(17, 1, 0)), 1e-3);
0511   CHECK_CLOSE_ABS(axis.min, 13, 1e-3);
0512   BOOST_CHECK_EQUAL(axis.bType, equidistant);
0513
0514   // multiple rings
0515   surfaces.resize(0);
0516   auto ringa = fullPhiTestSurfacesEC(30, 0, 0, 10);
0517   surfaces.insert(surfaces.end(), ringa.begin(), ringa.end());
0518   auto ringb = fullPhiTestSurfacesEC(30, 0, 0, 15);
0519   surfaces.insert(surfaces.end(), ringb.begin(), ringb.end());
0520   auto ringc = fullPhiTestSurfacesEC(30, 0, 0, 20);
0521   surfaces.insert(surfaces.end(), ringc.begin(), ringc.end());
0522   draw_surfaces(surfaces, "SurfaceArrayCreator_createEquidistantAxis_R_2.obj");
0523
0524   surfacesRaw = unpack_shared_vector(surfaces);
0525   pl = ProtoLayer(tgContext, surfacesRaw);
0526   trf = Transform3::Identity();
0527   axis = createEquidistantAxis(tgContext, surfacesRaw, AxisDirection::AxisR, pl,
0528                                trf);
0529
0530   BOOST_CHECK_EQUAL(axis.nBins, 3u);
0531   CHECK_CLOSE_REL(axis.max, perp(Vector3(20 + 2, 1, 0)), 1e-3);
0532   CHECK_CLOSE_ABS(axis.min, 8, 1e-3);
0533   BOOST_CHECK_EQUAL(axis.bType, equidistant);
0534 }
0535
0536 // if there are concentring disc or barrel modules, the bin count might be off
0537 // we want to create _as few bins_ as possible, meaning the r-ring with
0538 // the lowest number of surfaces should be used for the bin count or
0539 // as basis for the variable edge procedure
0540 // double filling will make sure no surfaces are dropped
0541 BOOST_FIXTURE_TEST_CASE(SurfaceArrayCreator_dependentBinCounts,
0542                         SurfaceArrayCreatorFixture) {
0543   auto ringA = fullPhiTestSurfacesEC(10, 0, 0, 10, 2, 3);
0544   auto ringB = fullPhiTestSurfacesEC(15, 0, 0, 15, 2, 3.5);
0545   auto ringC = fullPhiTestSurfacesEC(20, 0, 0, 20, 2, 3.8);
0546
0547   std::vector<std::shared_ptr<const Surface>> surfaces;
0548   std::copy(ringA.begin(), ringA.end(), std::back_inserter(surfaces));
0549   std::copy(ringB.begin(), ringB.end(), std::back_inserter(surfaces));
0550   std::copy(ringC.begin(), ringC.end(), std::back_inserter(surfaces));
0551   draw_surfaces(surfaces, "SurfaceArrayCreator_dependentBinCounts.obj");
0552
0553   std::unique_ptr<SurfaceArray> sArray =
0554       m_SAC.surfaceArrayOnDisc(tgContext, surfaces, equidistant, equidistant);
0555   auto axes = sArray->getAxes();
0556   BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), 3u);
0557   BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), 10u);
0558
0559   // Write the surrace array with grid
0560   ObjVisualization3D objVis;
0561   GeometryView3D::drawSurfaceArray(objVis, *sArray, tgContext);
0562   objVis.write("SurfaceArrayCreator_EndcapGrid");
0563 }
0564
0565 BOOST_FIXTURE_TEST_CASE(SurfaceArrayCreator_completeBinning,
0566                         SurfaceArrayCreatorFixture) {
0567   SrfVec brl = makeBarrel(30, 7, 2, 1);
0568   std::vector<const Surface*> brlRaw = unpack_shared_vector(brl);
0569   draw_surfaces(brl, "SurfaceArrayCreator_completeBinning_BRL.obj");
0570
0571   Axis<AxisType::Equidistant, AxisBoundaryType::Closed> phiAxis(
0572       -std::numbers::pi, std::numbers::pi, 30u);
0573   Axis<AxisType::Equidistant, AxisBoundaryType::Bound> zAxis(-14, 14, 7u);
0574
0575   double R = 10.;
0576   auto globalToLocal = [](const Vector3& pos) {
0577     return Vector2(phi(pos) + 2 * std::numbers::pi / 30 / 2, pos.z());
0578   };
0579   auto localToGlobal = [R](const Vector2& loc) {
0580     double phi = loc[0] - 2 * std::numbers::pi / 30 / 2;
0581     return Vector3(R * std::cos(phi), R * std::sin(phi), loc[1]);
0582   };
0583
0584   auto sl = std::make_unique<
0585       SurfaceArray::SurfaceGridLookup<decltype(phiAxis), decltype(zAxis)>>(
0586       globalToLocal, localToGlobal,
0587       std::make_tuple(std::move(phiAxis), std::move(zAxis)));
0588   sl->fill(tgContext, brlRaw);
0589   SurfaceArray sa(std::move(sl), brl);
0590
0591   // Write the surrace array with grid
0592   ObjVisualization3D objVis;
0593   GeometryView3D::drawSurfaceArray(objVis, sa, tgContext);
0594   objVis.write("SurfaceArrayCreator_BarrelGrid");
0595
0596   // actually filled SA
0597   for (const auto& srf : brl) {
0598     Vector3 ctr = srf->referencePosition(tgContext, AxisDirection::AxisR);
0599     auto binContent = sa.at(ctr);
0600
0601     BOOST_CHECK_EQUAL(binContent.size(), 1u);
0602     BOOST_CHECK_EQUAL(srf.get(), binContent.at(0));
0603   }
0604 }
0605
0606 BOOST_FIXTURE_TEST_CASE(SurfaceArrayCreator_barrelStagger,
0607                         SurfaceArrayCreatorFixture) {
0608   auto barrel = makeBarrelStagger(30, 7, 0, std::numbers::pi / 9.);
0609   auto brl = barrel.first;
0610   std::vector<const Surface*> brlRaw = unpack_shared_vector(brl);
0611   draw_surfaces(brl, "SurfaceArrayCreator_barrelStagger.obj");
0612
0613   ProtoLayer pl(tgContext, brl);
0614
0615   // EQUIDISTANT
0616   Transform3 tr = Transform3::Identity();
0617
0618   auto pAxisPhi =
0619       createEquidistantAxis(tgContext, brlRaw, AxisDirection::AxisPhi, pl, tr);
0620   auto pAxisZ =
0621       createEquidistantAxis(tgContext, brlRaw, AxisDirection::AxisZ, pl, tr);
0622
0623   double R = 10.;
0624   Transform3 itr = tr.inverse();
0625
0626   auto globalToLocal = [tr](const Vector3& pos) {
0627     Vector3 rot = tr * pos;
0628     return Vector2(phi(rot), rot.z());
0629   };
0630   auto localToGlobal = [R, itr](const Vector2& loc) {
0631     return itr * Vector3(R * std::cos(loc[0]), R * std::sin(loc[0]), loc[1]);
0632   };
0633
0634   auto sl = makeSurfaceGridLookup2D<AxisBoundaryType::Closed,
0635                                     AxisBoundaryType::Bound>(
0636       globalToLocal, localToGlobal, pAxisPhi, pAxisZ);
0637
0638   sl->fill(tgContext, brlRaw);
0639   SurfaceArray sa(std::move(sl), brl);
0640   auto axes = sa.getAxes();
0641   BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), 30u);
0642   BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), 7u);
0643
0644   for (const auto& pr : barrel.second) {
0645     auto A = pr.first;
0646     auto B = pr.second;
0647
0648     Vector3 ctr = A->referencePosition(tgContext, AxisDirection::AxisR);
0649     auto binContent = sa.at(ctr);
0650     BOOST_CHECK_EQUAL(binContent.size(), 2u);
0651     std::set<const Surface*> act;
0652     act.insert(binContent[0]);
0653     act.insert(binContent[1]);
0654
0655     std::set<const Surface*> exp;
0656     exp.insert(A);
0657     exp.insert(B);
0658
0659     BOOST_CHECK(act == exp);
0660   }
0661
0662   // VARIABLE
0663   BOOST_TEST_CONTEXT("Barrel Stagger Variable binning") {
0664     tr = Transform3::Identity();
0665
0666     auto pAxisPhiVar =
0667         createVariableAxis(tgContext, brlRaw, AxisDirection::AxisPhi, pl, tr);
0668     auto pAxisZVar =
0669         createVariableAxis(tgContext, brlRaw, AxisDirection::AxisZ, pl, tr);
0670
0671     itr = tr.inverse();
0672
0673     auto globalToLocalVar = [tr](const Vector3& pos) {
0674       Vector3 rot = tr * pos;
0675       return Vector2(phi(rot), rot.z());
0676     };
0677     auto localToGlobalVar = [R, itr](const Vector2& loc) {
0678       return itr * Vector3(R * std::cos(loc[0]), R * std::sin(loc[0]), loc[1]);
0679     };
0680
0681     auto sl2 = makeSurfaceGridLookup2D<AxisBoundaryType::Closed,
0682                                        AxisBoundaryType::Bound>(
0683         globalToLocalVar, localToGlobalVar, pAxisPhiVar, pAxisZVar);
0684
0685     sl2->fill(tgContext, brlRaw);
0686     SurfaceArray sa2(std::move(sl2), brl);
0687     axes = sa2.getAxes();
0688     BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), 30u);
0689     BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), 7u);
0690
0691     // check bin edges
0692     std::vector<double> phiEdgesExp = {
0693         -3.14159,  -2.93215,  -2.72271, -2.51327,    -2.30383, -2.0944,
0694         -1.88496,  -1.67552,  -1.46608, -1.25664,    -1.0472,  -0.837758,
0695         -0.628319, -0.418879, -0.20944, 4.44089e-16, 0.20944,  0.418879,
0696         0.628319,  0.837758,  1.0472,   1.25664,     1.46608,  1.67552,
0697         1.88496,   2.0944,    2.30383,  2.51327,     2.72271,  3.00831,
0698         3.14159};
0699     std::vector<double> zEdgesExp = {-14, -10, -6, -2, 2, 6, 10, 14};
0700     std::size_t i = 0;
0701     for (const auto& edge : axes.at(0)->getBinEdges()) {
0702       BOOST_TEST_INFO("phi edge index " << i);
0703       auto phiEdge = phiEdgesExp.at(i);
0704       CHECK_CLOSE_ABS(edge, phiEdge, 1e-5 * std::numbers::pi);
0705       i++;
0706     }
0707     i = 0;
0708     for (const auto& edge : axes.at(1)->getBinEdges()) {
0709       BOOST_TEST_INFO("z edge index " << i);
0710       CHECK_CLOSE_ABS(edge, zEdgesExp.at(i), 1e-6);
0711       i++;
0712     }
0713
0714     for (const auto& pr : barrel.second) {
0715       auto A = pr.first;
0716       auto B = pr.second;
0717
0718       Vector3 ctr = A->referencePosition(tgContext, AxisDirection::AxisR);
0719       auto binContent = sa2.at(ctr);
0720       BOOST_CHECK_EQUAL(binContent.size(), 2u);
0721       std::set<const Surface*> act;
0722       act.insert(binContent[0]);
0723       act.insert(binContent[1]);
0724
0725       std::set<const Surface*> exp;
0726       exp.insert(A);
0727       exp.insert(B);
0728
0729       BOOST_CHECK(act == exp);
0730     }
0731   }
0732 }
0733
0734 BOOST_AUTO_TEST_SUITE_END()
0735 }  // namespace Acts::Test