Core/Surfaces/SurfaceArrayTests.cpp

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