Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-01-06 09:23:56 

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