Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:11:22 

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