Core/Surfaces/SurfaceArrayTests.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/GeometryContext.hpp"
0013 #include "Acts/Surfaces/PlanarBounds.hpp"
0014 #include "Acts/Surfaces/PlaneSurface.hpp"
0015 #include "Acts/Surfaces/RectangleBounds.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017 #include "Acts/Surfaces/SurfaceArray.hpp"
0018 #include "Acts/Utilities/Axis.hpp"
0019 #include "Acts/Utilities/AxisDefinitions.hpp"
0020 #include "Acts/Utilities/BinningType.hpp"
0021 #include "Acts/Utilities/Helpers.hpp"
0022
0023 #include <cmath>
0024 #include <cstddef>
0025 #include <fstream>
0026 #include <iomanip>
0027 #include <iostream>
0028 #include <memory>
0029 #include <numbers>
0030 #include <string>
0031 #include <tuple>
0032 #include <utility>
0033 #include <vector>
0034
0035 #include <boost/format.hpp>
0036
0037 using Acts::VectorHelpers::phi;
0038
0039 namespace Acts::Test {
0040
0041 // Create a test context
0042 GeometryContext tgContext = GeometryContext();
0043
0044 using SrfVec = std::vector<std::shared_ptr<const Surface>>;
0045 struct SurfaceArrayFixture {
0046   std::vector<std::shared_ptr<const Surface>> m_surfaces;
0047
0048   SurfaceArrayFixture() { BOOST_TEST_MESSAGE("setup fixture"); }
0049   ~SurfaceArrayFixture() { BOOST_TEST_MESSAGE("teardown fixture"); }
0050
0051   SrfVec fullPhiTestSurfacesEC(std::size_t n = 10, double shift = 0,
0052                                double zbase = 0, double r = 10) {
0053     SrfVec res;
0054
0055     double phiStep = 2 * std::numbers::pi / n;
0056     for (std::size_t i = 0; i < n; ++i) {
0057       double z = zbase + ((i % 2 == 0) ? 1 : -1) * 0.2;
0058
0059       Transform3 trans;
0060       trans.setIdentity();
0061       trans.rotate(Eigen::AngleAxisd(i * phiStep + shift, Vector3(0, 0, 1)));
0062       trans.translate(Vector3(r, 0, z));
0063
0064       auto bounds = std::make_shared<const RectangleBounds>(2, 1);
0065       std::shared_ptr<const Surface> srf =
0066           Surface::makeShared<PlaneSurface>(trans, bounds);
0067
0068       res.push_back(srf);
0069       m_surfaces.push_back(
0070           std::move(srf));  // keep shared, will get destroyed at the end
0071     }
0072
0073     return res;
0074   }
0075
0076   SrfVec fullPhiTestSurfacesBRL(int n = 10, double shift = 0, double zbase = 0,
0077                                 double incl = std::numbers::pi / 9.,
0078                                 double w = 2, double h = 1.5) {
0079     SrfVec res;
0080
0081     double phiStep = 2 * std::numbers::pi / n;
0082     for (int i = 0; i < n; ++i) {
0083       double z = zbase;
0084
0085       Transform3 trans;
0086       trans.setIdentity();
0087       trans.rotate(Eigen::AngleAxisd(i * phiStep + shift, Vector3(0, 0, 1)));
0088       trans.translate(Vector3(10, 0, z));
0089       trans.rotate(Eigen::AngleAxisd(incl, Vector3(0, 0, 1)));
0090       trans.rotate(Eigen::AngleAxisd(std::numbers::pi / 2., Vector3(0, 1, 0)));
0091
0092       auto bounds = std::make_shared<const RectangleBounds>(w, h);
0093       std::shared_ptr<const Surface> srf =
0094           Surface::makeShared<PlaneSurface>(trans, bounds);
0095
0096       res.push_back(srf);
0097       m_surfaces.push_back(
0098           std::move(srf));  // keep shared, will get destroyed at the end
0099     }
0100
0101     return res;
0102   }
0103
0104   SrfVec straightLineSurfaces(
0105       std::size_t n = 10., double step = 3, const Vector3& origin = {0, 0, 1.5},
0106       const Transform3& pretrans = Transform3::Identity(),
0107       const Vector3& dir = {0, 0, 1}) {
0108     SrfVec res;
0109     for (std::size_t i = 0; i < n; ++i) {
0110       Transform3 trans;
0111       trans.setIdentity();
0112       trans.translate(origin + dir * step * i);
0113       // trans.rotate(AngleAxis3(std::numbers::pi/9., Vector3(0, 0, 1)));
0114       trans.rotate(AngleAxis3(std::numbers::pi / 2., Vector3(1, 0, 0)));
0115       trans = trans * pretrans;
0116
0117       auto bounds = std::make_shared<const RectangleBounds>(2, 1.5);
0118
0119       std::shared_ptr<const Surface> srf =
0120           Surface::makeShared<PlaneSurface>(trans, bounds);
0121
0122       res.push_back(srf);
0123       m_surfaces.push_back(
0124           std::move(srf));  // keep shared, will get destroyed at the end
0125     }
0126
0127     return res;
0128   }
0129
0130   SrfVec makeBarrel(int nPhi, int nZ, double w, double h) {
0131     double z0 = -(nZ - 1) * w;
0132     SrfVec res;
0133
0134     for (int i = 0; i < nZ; i++) {
0135       double z = i * w * 2 + z0;
0136       SrfVec ring =
0137           fullPhiTestSurfacesBRL(nPhi, 0, z, std::numbers::pi / 9., w, h);
0138       res.insert(res.end(), ring.begin(), ring.end());
0139     }
0140
0141     return res;
0142   }
0143
0144   void draw_surfaces(const SrfVec& surfaces, const std::string& fname) {
0145     std::ofstream os;
0146     os.open(fname);
0147
0148     os << std::fixed << std::setprecision(4);
0149
0150     std::size_t nVtx = 0;
0151     for (const auto& srfx : surfaces) {
0152       std::shared_ptr<const PlaneSurface> srf =
0153           std::dynamic_pointer_cast<const PlaneSurface>(srfx);
0154       const PlanarBounds* bounds =
0155           dynamic_cast<const PlanarBounds*>(&srf->bounds());
0156
0157       for (const auto& vtxloc : bounds->vertices()) {
0158         Vector3 vtx =
0159             srf->transform(tgContext) * Vector3(vtxloc.x(), vtxloc.y(), 0);
0160         os << "v " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
0161       }
0162
0163       // connect them
0164       os << "f";
0165       for (std::size_t i = 1; i <= bounds->vertices().size(); ++i) {
0166         os << " " << nVtx + i;
0167       }
0168       os << "\n";
0169
0170       nVtx += bounds->vertices().size();
0171     }
0172
0173     os.close();
0174   }
0175 };
0176
0177 BOOST_AUTO_TEST_SUITE(Surfaces)
0178
0179 BOOST_FIXTURE_TEST_CASE(SurfaceArray_create, SurfaceArrayFixture) {
0180   GeometryContext tgContext = GeometryContext();
0181
0182   SrfVec brl = makeBarrel(30, 7, 2, 1);
0183   std::vector<const Surface*> brlRaw = unpack_shared_vector(brl);
0184   draw_surfaces(brl, "SurfaceArray_create_BRL_1.obj");
0185
0186   Axis<AxisType::Equidistant, AxisBoundaryType::Closed> phiAxis(
0187       -std::numbers::pi, std::numbers::pi, 30u);
0188   Axis<AxisType::Equidistant, AxisBoundaryType::Bound> zAxis(-14, 14, 7u);
0189
0190   double angleShift = 2 * std::numbers::pi / 30. / 2.;
0191   auto transform = [angleShift](const Vector3& pos) {
0192     return Vector2(phi(pos) + angleShift, pos.z());
0193   };
0194   double R = 10;
0195   auto itransform = [angleShift, R](const Vector2& loc) {
0196     return Vector3(R * std::cos(loc[0] - angleShift),
0197                    R * std::sin(loc[0] - angleShift), loc[1]);
0198   };
0199
0200   auto sl = std::make_unique<
0201       SurfaceArray::SurfaceGridLookup<decltype(phiAxis), decltype(zAxis)>>(
0202       transform, itransform,
0203       std::make_tuple(std::move(phiAxis), std::move(zAxis)));
0204   sl->fill(tgContext, brlRaw);
0205   SurfaceArray sa(std::move(sl), brl);
0206
0207   // let's see if we can access all surfaces
0208   sa.toStream(tgContext, std::cout);
0209
0210   for (const auto& srf : brl) {
0211     Vector3 ctr = srf->referencePosition(tgContext, AxisDirection::AxisR);
0212     std::vector<const Surface*> binContent = sa.at(ctr);
0213
0214     BOOST_CHECK_EQUAL(binContent.size(), 1u);
0215     BOOST_CHECK_EQUAL(srf.get(), binContent.at(0));
0216   }
0217
0218   std::vector<const Surface*> neighbors =
0219       sa.neighbors(itransform(Vector2(0, 0)));
0220   BOOST_CHECK_EQUAL(neighbors.size(), 9u);
0221
0222   auto sl2 = std::make_unique<
0223       SurfaceArray::SurfaceGridLookup<decltype(phiAxis), decltype(zAxis)>>(
0224       transform, itransform,
0225       std::make_tuple(std::move(phiAxis), std::move(zAxis)));
0226   // do NOT fill, only complete binning
0227   sl2->completeBinning(tgContext, brlRaw);
0228   SurfaceArray sa2(std::move(sl2), brl);
0229   sa.toStream(tgContext, std::cout);
0230   for (const auto& srf : brl) {
0231     Vector3 ctr = srf->referencePosition(tgContext, AxisDirection::AxisR);
0232     std::vector<const Surface*> binContent = sa2.at(ctr);
0233
0234     BOOST_CHECK_EQUAL(binContent.size(), 1u);
0235     BOOST_CHECK_EQUAL(srf.get(), binContent.at(0));
0236   }
0237 }
0238
0239 BOOST_AUTO_TEST_CASE(SurfaceArray_singleElement) {
0240   const double w = 3;
0241   const double h = 4;
0242   auto bounds = std::make_shared<const RectangleBounds>(w, h);
0243   auto srf = Surface::makeShared<PlaneSurface>(Transform3::Identity(), bounds);
0244
0245   SurfaceArray sa(srf);
0246
0247   auto binContent = sa.at(Vector3(42, 42, 42));
0248   BOOST_CHECK_EQUAL(binContent.size(), 1u);
0249   BOOST_CHECK_EQUAL(binContent.at(0), srf.get());
0250   BOOST_CHECK_EQUAL(sa.surfaces().size(), 1u);
0251   BOOST_CHECK_EQUAL(sa.surfaces().at(0), srf.get());
0252 }
0253
0254 BOOST_AUTO_TEST_CASE(SurfaceArray_manyElementsSingleLookup) {
0255   const double w = 3;
0256   const double h = 4;
0257   auto bounds = std::make_shared<const RectangleBounds>(w, h);
0258   auto srf0 = Surface::makeShared<PlaneSurface>(Transform3::Identity(), bounds);
0259   auto srf1 = Surface::makeShared<PlaneSurface>(Transform3::Identity(), bounds);
0260
0261   std::vector<const Surface*> sfPointers = {srf0.get(), srf1.get()};
0262   std::vector<std::shared_ptr<const Surface>> surfaces = {srf0, srf1};
0263
0264   auto singleLookUp =
0265       std::make_unique<Acts::SurfaceArray::SingleElementLookup>(sfPointers);
0266
0267   SurfaceArray sa(std::move(singleLookUp), surfaces);
0268
0269   auto binContent = sa.at(Vector3(42, 42, 42));
0270   BOOST_CHECK_EQUAL(binContent.size(), 2u);
0271   BOOST_CHECK_EQUAL(sa.surfaces().size(), 2u);
0272 }
0273
0274 BOOST_AUTO_TEST_SUITE_END()
0275 }  // namespace Acts::Test