ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SurfaceArrayCreatorTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SurfaceArrayCreatorTests.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2017-2018 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include <boost/format.hpp>
10 #include <boost/test/data/test_case.hpp>
11 #include <boost/test/unit_test.hpp>
12 
13 #include <fstream>
14 #include <random>
15 
16 #include "Acts/Geometry/GeometryContext.hpp"
17 #include "Acts/Geometry/SurfaceArrayCreator.hpp"
18 #include "Acts/Surfaces/PlaneSurface.hpp"
19 #include "Acts/Surfaces/RectangleBounds.hpp"
20 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
21 #include "Acts/Utilities/BinningType.hpp"
22 #include "Acts/Utilities/Definitions.hpp"
23 
24 using Acts::VectorHelpers::perp;
25 using Acts::VectorHelpers::phi;
26 
27 namespace bdata = boost::unit_test::data;
28 namespace tt = boost::test_tools;
29 
30 namespace Acts {
31 
32 namespace Test {
33 
34 // Create a test context
35 GeometryContext tgContext = GeometryContext();
36 
37 #define CHECK_ROTATION_ANGLE(t, a, tolerance) \
38  { \
39  Vector3D v = (*t) * Vector3D(1, 0, 0); \
40  CHECK_CLOSE_ABS(phi(v), (a), tolerance); \
41  }
42 
43 using SrfVec = std::vector<std::shared_ptr<const Surface>>;
44 
45 struct SurfaceArrayCreatorFixture {
46  SurfaceArrayCreator m_SAC;
47  std::vector<std::shared_ptr<const Surface>> m_surfaces;
48 
49  SurfaceArrayCreatorFixture()
50  : m_SAC(SurfaceArrayCreator::Config(),
51  Acts::getDefaultLogger("SurfaceArrayCreator",
52  Acts::Logging::VERBOSE)) {
53  BOOST_TEST_MESSAGE("setup fixture");
54  }
55  ~SurfaceArrayCreatorFixture() { BOOST_TEST_MESSAGE("teardown fixture"); }
56 
57  template <typename... Args>
58  SurfaceArrayCreator::ProtoAxis createEquidistantAxis(Args&&... args) {
59  return m_SAC.createEquidistantAxis(std::forward<Args>(args)...);
60  }
61 
62  template <typename... Args>
63  SurfaceArrayCreator::ProtoAxis createVariableAxis(Args&&... args) {
64  return m_SAC.createVariableAxis(std::forward<Args>(args)...);
65  }
66 
67  template <detail::AxisBoundaryType bdtA, detail::AxisBoundaryType bdtB,
68  typename... Args>
69  std::unique_ptr<SurfaceArray::ISurfaceGridLookup> makeSurfaceGridLookup2D(
70  Args&&... args) {
71  return m_SAC.makeSurfaceGridLookup2D<bdtA, bdtB>(
72  std::forward<Args>(args)...);
73  }
74 
75  SrfVec fullPhiTestSurfacesEC(size_t n = 10, double shift = 0,
76  double zbase = 0, double r = 10, double w = 2,
77  double h = 1) {
78  SrfVec res;
79  // TODO: The test is extremely numerically unstable in the face of upward
80  // rounding in this multiplication and division. Find out why.
81  double phiStep = 2 * M_PI / n;
82  for (size_t i = 0; i < n; ++i) {
83  double z = zbase + ((i % 2 == 0) ? 1 : -1) * 0.2;
84  double phi = std::fma(i, phiStep, shift);
85 
86  Transform3D trans;
87  trans.setIdentity();
88  trans.rotate(Eigen::AngleAxisd(phi, Vector3D(0, 0, 1)));
89  trans.translate(Vector3D(r, 0, z));
90 
91  auto bounds = std::make_shared<const RectangleBounds>(w, h);
92 
93  auto transptr = std::make_shared<const Transform3D>(trans);
94  std::shared_ptr<Surface> srf =
95  Surface::makeShared<PlaneSurface>(transptr, bounds);
96 
97  res.push_back(srf);
98  m_surfaces.push_back(
99  std::move(srf)); // keep shared, will get destroyed at the end
100  }
101 
102  return res;
103  }
104 
105  SrfVec fullPhiTestSurfacesBRL(size_t n = 10, double shift = 0,
106  double zbase = 0, double incl = M_PI / 9.,
107  double w = 2, double h = 1.5) {
108  SrfVec res;
109  // TODO: The test is extremely numerically unstable in the face of upward
110  // rounding in this multiplication and division. Find out why.
111  double phiStep = 2 * M_PI / n;
112  for (size_t i = 0; i < n; ++i) {
113  double z = zbase;
114  double phi = std::fma(i, phiStep, shift);
115 
116  Transform3D trans;
117  trans.setIdentity();
118  trans.rotate(Eigen::AngleAxisd(phi, Vector3D(0, 0, 1)));
119  trans.translate(Vector3D(10, 0, z));
120  trans.rotate(Eigen::AngleAxisd(incl, Vector3D(0, 0, 1)));
121  trans.rotate(Eigen::AngleAxisd(M_PI / 2., Vector3D(0, 1, 0)));
122 
123  auto bounds = std::make_shared<const RectangleBounds>(w, h);
124 
125  auto transptr = std::make_shared<const Transform3D>(trans);
126  std::shared_ptr<Surface> srf =
127  Surface::makeShared<PlaneSurface>(transptr, bounds);
128 
129  res.push_back(srf);
130  m_surfaces.push_back(
131  std::move(srf)); // keep shared, will get destroyed at the end
132  }
133 
134  return res;
135  }
136 
137  SrfVec straightLineSurfaces(
138  size_t n = 10., double step = 3, const Vector3D& origin = {0, 0, 1.5},
139  const Transform3D& pretrans = Transform3D::Identity(),
140  const Vector3D& dir = {0, 0, 1}) {
141  SrfVec res;
142  for (size_t i = 0; i < n; ++i) {
143  Transform3D trans;
144  trans.setIdentity();
145  trans.translate(origin + dir * step * i);
146  // trans.rotate(AngleAxis3D(M_PI/9., Vector3D(0, 0, 1)));
147  trans.rotate(AngleAxis3D(M_PI / 2., Vector3D(1, 0, 0)));
148  trans = trans * pretrans;
149 
150  auto bounds = std::make_shared<const RectangleBounds>(2, 1.5);
151 
152  auto transptr = std::make_shared<const Transform3D>(trans);
153  std::shared_ptr<Surface> srf =
154  Surface::makeShared<PlaneSurface>(transptr, bounds);
155 
156  res.push_back(srf);
157  m_surfaces.push_back(
158  std::move(srf)); // keep shared, will get destroyed at the end
159  }
160 
161  return res;
162  }
163 
164  SrfVec makeBarrel(int nPhi, int nZ, double w, double h) {
165  double z0 = -(nZ - 1) * w;
166  SrfVec res;
167 
168  for (int i = 0; i < nZ; i++) {
169  double z = i * w * 2 + z0;
170  // std::cout << "z=" << z << std::endl;
171  SrfVec ring = fullPhiTestSurfacesBRL(nPhi, 0, z, M_PI / 9., w, h);
172  res.insert(res.end(), ring.begin(), ring.end());
173  }
174 
175  return res;
176  }
177 
178  std::pair<SrfVec, std::vector<std::pair<const Surface*, const Surface*>>>
179  makeBarrelStagger(int nPhi, int nZ, double shift = 0, double incl = M_PI / 9.,
180  double w = 2, double h = 1.5) {
181  double z0 = -(nZ - 1) * w;
182  SrfVec res;
183  std::vector<std::pair<const Surface*, const Surface*>> pairs;
184  // TODO: The test is extremely numerically unstable in the face of upward
185  // rounding in this multiplication and division. Find out why.
186  double phiStep = 2 * M_PI / nPhi;
187  for (int i = 0; i < nZ; i++) {
188  double z = i * w * 2 + z0;
189  for (int j = 0; j < nPhi; ++j) {
190  double phi = std::fma(j, phiStep, shift);
191  Transform3D trans;
192  trans.setIdentity();
193  trans.rotate(Eigen::AngleAxisd(phi, Vector3D(0, 0, 1)));
194  trans.translate(Vector3D(10, 0, z));
195  trans.rotate(Eigen::AngleAxisd(incl, Vector3D(0, 0, 1)));
196  trans.rotate(Eigen::AngleAxisd(M_PI / 2., Vector3D(0, 1, 0)));
197 
198  auto bounds = std::make_shared<const RectangleBounds>(w, h);
199 
200  auto transAptr = std::make_shared<const Transform3D>(trans);
201 
202  std::shared_ptr<Surface> srfA =
203  Surface::makeShared<PlaneSurface>(transAptr, bounds);
204 
205  Vector3D nrm = srfA->normal(tgContext);
206  Transform3D transB = trans;
207  transB.pretranslate(nrm * 0.1);
208  auto transBptr = std::make_shared<const Transform3D>(transB);
209  std::shared_ptr<Surface> srfB =
210  Surface::makeShared<PlaneSurface>(transBptr, bounds);
211 
212  pairs.push_back(std::make_pair(srfA.get(), srfB.get()));
213 
214  res.push_back(srfA);
215  res.push_back(srfB);
216  m_surfaces.push_back(std::move(srfA));
217  m_surfaces.push_back(std::move(srfB));
218  }
219  }
220 
221  return std::make_pair(res, pairs);
222  }
223 };
224 
225 void draw_surfaces(SrfVec surfaces, const std::string& fname) {
226  std::ofstream os;
227  os.open(fname);
228 
229  os << std::fixed << std::setprecision(4);
230 
231  size_t nVtx = 0;
232  for (const auto& srfx : surfaces) {
233  std::shared_ptr<const PlaneSurface> srf =
234  std::dynamic_pointer_cast<const PlaneSurface>(srfx);
235  const PlanarBounds* bounds =
236  dynamic_cast<const PlanarBounds*>(&srf->bounds());
237 
238  for (const auto& vtxloc : bounds->vertices()) {
239  Vector3D vtx =
240  srf->transform(tgContext) * Vector3D(vtxloc.x(), vtxloc.y(), 0);
241  os << "v " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
242  }
243 
244  // connect them
245  os << "f";
246  for (size_t i = 1; i <= bounds->vertices().size(); ++i) {
247  os << " " << nVtx + i;
248  }
249  os << "\n";
250 
251  nVtx += bounds->vertices().size();
252  }
253 
254  os.close();
255 }
256 
257 BOOST_AUTO_TEST_SUITE(Tools)
258 
259 BOOST_FIXTURE_TEST_CASE(SurfaceArrayCreator_createEquidistantAxis_Phi,
260  SurfaceArrayCreatorFixture) {
261  // fail on empty srf vector
262  std::vector<const Surface*> emptyRaw;
263  ProtoLayer pl(tgContext, emptyRaw);
264  auto tr = Transform3D::Identity();
265  BOOST_CHECK_THROW(
266  createEquidistantAxis(tgContext, emptyRaw, BinningValue::binPhi, pl, tr),
267  std::logic_error);
268 
269  std::vector<float> bdExp = {
270  -3.14159, -2.93215, -2.72271, -2.51327, -2.30383, -2.0944, -1.88496,
271  -1.67552, -1.46608, -1.25664, -1.0472, -0.837758, -0.628319, -0.418879,
272  -0.20944, 0, 0.20944, 0.418879, 0.628319, 0.837758, 1.0472,
273  1.25664, 1.46608, 1.67552, 1.88496, 2.09439, 2.30383, 2.51327,
274  2.72271, 2.93215, 3.14159};
275 
276  double step = 2 * M_PI / 30.;
277 
278  // endcap style modules
279 
280  for (int i = -1; i <= 2; i += 2) {
281  double z = 10 * i;
282  // case 1: one module sits at pi / -pi
283  double angleShift = step / 2.;
284  auto surfaces = fullPhiTestSurfacesEC(30, angleShift, z);
285  std::vector<const Surface*> surfacesRaw = unpack_shared_vector(surfaces);
286  pl = ProtoLayer(tgContext, surfacesRaw);
287  tr = Transform3D::Identity();
288  auto axis = createEquidistantAxis(tgContext, surfacesRaw,
289  BinningValue::binPhi, pl, tr);
290 
291  BOOST_CHECK_EQUAL(axis.nBins, 30u);
292  CHECK_CLOSE_REL(axis.max, M_PI, 1e-6);
293  CHECK_CLOSE_REL(axis.min, -M_PI, 1e-6);
294  BOOST_CHECK_EQUAL(axis.bType, equidistant);
295  CHECK_SMALL(phi(tr * Vector3D::UnitX()), 1e-6);
296 
297  // case 2: two modules sit symmetrically around pi / -pi
298  angleShift = 0.;
299  surfaces = fullPhiTestSurfacesEC(30, angleShift, z);
300  surfacesRaw = unpack_shared_vector(surfaces);
301  pl = ProtoLayer(tgContext, surfacesRaw);
302  tr = Transform3D::Identity();
303  axis = createEquidistantAxis(tgContext, surfacesRaw, BinningValue::binPhi,
304  pl, tr);
305  draw_surfaces(surfaces,
306  "SurfaceArrayCreator_createEquidistantAxis_EC_2.obj");
307  BOOST_CHECK_EQUAL(axis.nBins, 30u);
308  CHECK_CLOSE_REL(axis.max, M_PI, 1e-6);
309  CHECK_CLOSE_REL(axis.min, -M_PI, 1e-6);
310  BOOST_CHECK_EQUAL(axis.bType, equidistant);
311  // CHECK_CLOSE_REL(bdExp, axis.binEdges, 0.001);
312  CHECK_CLOSE_REL(phi(tr * Vector3D::UnitX()), -0.5 * step, 1e-3);
313  // case 3: two modules sit asymmetrically around pi / -pi shifted up
314  angleShift = step / -4.;
315  surfaces = fullPhiTestSurfacesEC(30, angleShift, z);
316  surfacesRaw = unpack_shared_vector(surfaces);
317  pl = ProtoLayer(tgContext, surfacesRaw);
318  tr = Transform3D::Identity();
319  axis = createEquidistantAxis(tgContext, surfacesRaw, BinningValue::binPhi,
320  pl, tr);
321  draw_surfaces(surfaces,
322  "SurfaceArrayCreator_createEquidistantAxis_EC_3.obj");
323  BOOST_CHECK_EQUAL(axis.nBins, 30u);
324  CHECK_CLOSE_REL(axis.max, M_PI, 1e-6);
325  CHECK_CLOSE_REL(axis.min, -M_PI, 1e-6);
326  BOOST_CHECK_EQUAL(axis.bType, equidistant);
327  CHECK_CLOSE_REL(phi(tr * Vector3D::UnitX()), step / -4., 1e-3);
328 
329  // case 4: two modules sit asymmetrically around pi / -pi shifted down
330  angleShift = step / 4.;
331  surfaces = fullPhiTestSurfacesEC(30, angleShift, z);
332  surfacesRaw = unpack_shared_vector(surfaces);
333  pl = ProtoLayer(tgContext, surfaces);
334  surfacesRaw = unpack_shared_vector(surfaces);
335  tr = Transform3D::Identity();
336  axis = createEquidistantAxis(tgContext, surfacesRaw, BinningValue::binPhi,
337  pl, tr);
338  surfacesRaw = unpack_shared_vector(surfaces);
339  draw_surfaces(surfaces,
340  "SurfaceArrayCreator_createEquidistantAxis_EC_4.obj");
341  BOOST_CHECK_EQUAL(axis.nBins, 30u);
342  CHECK_CLOSE_REL(axis.max, M_PI, 1e-6);
343  CHECK_CLOSE_REL(axis.min, -M_PI, 1e-6);
344  BOOST_CHECK_EQUAL(axis.bType, equidistant);
345  CHECK_CLOSE_REL(phi(tr * Vector3D::UnitX()), step / 4., 1e-3);
346  }
347 
348  for (int i = -1; i <= 2; i += 2) {
349  double z = 10 * i;
350  // case 1: one module sits at pi / -pi
351  double angleShift = step / 2.;
352  auto surfaces = fullPhiTestSurfacesBRL(30, angleShift, z);
353  auto surfacesRaw = unpack_shared_vector(surfaces);
354  pl = ProtoLayer(tgContext, surfacesRaw);
355  tr = Transform3D::Identity();
356  auto axis = createEquidistantAxis(tgContext, surfacesRaw,
357  BinningValue::binPhi, pl, tr);
358  draw_surfaces(surfaces,
359  "SurfaceArrayCreator_createEquidistantAxis_BRL_1.obj");
360  BOOST_CHECK_EQUAL(axis.nBins, 30u);
361  CHECK_CLOSE_REL(axis.max, M_PI, 1e-6);
362  CHECK_CLOSE_REL(axis.min, -M_PI, 1e-6);
363  BOOST_CHECK_EQUAL(axis.bType, equidistant);
364  CHECK_SMALL(phi(tr * Vector3D::UnitX()), 1e-6);
365 
366  // case 2: two modules sit symmetrically around pi / -pi
367  angleShift = 0.;
368  surfaces = fullPhiTestSurfacesBRL(30, angleShift, z);
369  surfacesRaw = unpack_shared_vector(surfaces);
370  pl = ProtoLayer(tgContext, surfacesRaw);
371  tr = Transform3D::Identity();
372  axis = createEquidistantAxis(tgContext, surfacesRaw, BinningValue::binPhi,
373  pl, tr);
374  draw_surfaces(surfaces,
375  "SurfaceArrayCreator_createEquidistantAxis_BRL_2.obj");
376  BOOST_CHECK_EQUAL(axis.nBins, 30u);
377  CHECK_CLOSE_REL(axis.max, M_PI, 1e-6);
378  CHECK_CLOSE_REL(axis.min, -M_PI, 1e-6);
379  BOOST_CHECK_EQUAL(axis.bType, equidistant);
380  // CHECK_CLOSE_REL(bdExp, axis.binEdges, 0.001);
381  CHECK_CLOSE_REL(phi(tr * Vector3D::UnitX()), -0.5 * step, 1e-3);
382 
383  // case 3: two modules sit asymmetrically around pi / -pi shifted up
384  angleShift = step / -4.;
385  surfaces = fullPhiTestSurfacesBRL(30, angleShift, z);
386  surfacesRaw = unpack_shared_vector(surfaces);
387  pl = ProtoLayer(tgContext, surfacesRaw);
388  tr = Transform3D::Identity();
389  axis = createEquidistantAxis(tgContext, surfacesRaw, BinningValue::binPhi,
390  pl, tr);
391  draw_surfaces(surfaces,
392  "SurfaceArrayCreator_createEquidistantAxis_BRL_3.obj");
393  BOOST_CHECK_EQUAL(axis.nBins, 30u);
394  CHECK_CLOSE_REL(axis.max, M_PI, 1e-6);
395  CHECK_CLOSE_REL(axis.min, -M_PI, 1e-6);
396  BOOST_CHECK_EQUAL(axis.bType, equidistant);
397  // CHECK_CLOSE_REL(bdExp, axis.binEdges, 0.001);
398  CHECK_CLOSE_REL(phi(tr * Vector3D::UnitX()), step / -4., 1e-3);
399 
400  // case 4: two modules sit asymmetrically around pi / -pi shifted down
401  angleShift = step / 4.;
402  surfaces = fullPhiTestSurfacesBRL(30, angleShift, z);
403  surfacesRaw = unpack_shared_vector(surfaces);
404  pl = ProtoLayer(tgContext, surfacesRaw);
405  tr = Transform3D::Identity();
406  axis = createEquidistantAxis(tgContext, surfacesRaw, BinningValue::binPhi,
407  pl, tr);
408  draw_surfaces(surfaces,
409  "SurfaceArrayCreator_createEquidistantAxis_BRL_4.obj");
410  BOOST_CHECK_EQUAL(axis.nBins, 30u);
411  CHECK_CLOSE_REL(axis.max, M_PI, 1e-6);
412  CHECK_CLOSE_REL(axis.min, -M_PI, 1e-6);
413  BOOST_CHECK_EQUAL(axis.bType, equidistant);
414  // CHECK_CLOSE_REL(bdExp, axis.binEdges, 0.001);
415  CHECK_CLOSE_REL(phi(tr * Vector3D::UnitX()), step / 4., 1e-3);
416  }
417 
418  SrfVec surfaces;
419 
420  // single element in phi
421  surfaces = fullPhiTestSurfacesEC(1);
422  auto surfacesRaw = unpack_shared_vector(surfaces);
423  draw_surfaces(surfaces,
424  "SurfaceArrayCreator_createEquidistantAxis_EC_Single.obj");
425 
426  pl = ProtoLayer(tgContext, surfacesRaw);
427  tr = Transform3D::Identity();
428  auto axis = createEquidistantAxis(tgContext, surfacesRaw,
429  BinningValue::binPhi, pl, tr);
430  BOOST_CHECK_EQUAL(axis.nBins, 1u);
431 
432  CHECK_CLOSE_ABS(axis.max, phi(Vector3D(8, 1, 0)), 1e-3);
433  CHECK_CLOSE_ABS(axis.min, phi(Vector3D(8, -1, 0)), 1e-3);
434  BOOST_CHECK_EQUAL(axis.bType, equidistant);
435 }
436 
437 BOOST_FIXTURE_TEST_CASE(SurfaceArrayCreator_createEquidistantAxis_Z,
438  SurfaceArrayCreatorFixture) {
439  // single element in z
440  auto surfaces = straightLineSurfaces(1);
441  auto surfacesRaw = unpack_shared_vector(surfaces);
442  ProtoLayer pl = ProtoLayer(tgContext, surfacesRaw);
443  auto trf = Transform3D::Identity();
444  auto axis = createEquidistantAxis(tgContext, surfacesRaw, BinningValue::binZ,
445  pl, trf);
446  draw_surfaces(surfaces, "SurfaceArrayCreator_createEquidistantAxis_Z_1.obj");
447  BOOST_CHECK_EQUAL(axis.nBins, 1u);
448  CHECK_CLOSE_ABS(axis.max, 3, 1e-6);
449  CHECK_CLOSE_ABS(axis.min, 0, 1e-6);
450  BOOST_CHECK_EQUAL(axis.bType, equidistant);
451 
452  // z rows with varying starting point
453  for (size_t i = 0; i <= 20; i++) {
454  double z0 = -10 + 1. * i;
455  surfaces = straightLineSurfaces(10, 3, Vector3D(0, 0, z0 + 1.5));
456  surfacesRaw = unpack_shared_vector(surfaces);
457  pl = ProtoLayer(tgContext, surfacesRaw);
458  trf = Transform3D::Identity();
459  axis = createEquidistantAxis(tgContext, surfacesRaw, BinningValue::binZ, pl,
460  trf);
461  draw_surfaces(
462  surfaces,
463  (boost::format(
464  "SurfaceArrayCreator_createEquidistantAxis_Z_2_%1%.obj") %
465  i)
466  .str());
467  BOOST_CHECK_EQUAL(axis.nBins, 10u);
468  CHECK_CLOSE_ABS(axis.max, 30 + z0, 1e-6);
469  CHECK_CLOSE_ABS(axis.min, z0, 1e-6);
470  BOOST_CHECK_EQUAL(axis.bType, equidistant);
471  }
472 
473  // z row where elements are rotated around y
474  Transform3D tr = Transform3D::Identity();
475  tr.rotate(AngleAxis3D(M_PI / 4., Vector3D(0, 0, 1)));
476  surfaces = straightLineSurfaces(10, 3, Vector3D(0, 0, 0 + 1.5), tr);
477  surfacesRaw = unpack_shared_vector(surfaces);
478  pl = ProtoLayer(tgContext, surfacesRaw);
479  trf = Transform3D::Identity();
480  axis = createEquidistantAxis(tgContext, surfacesRaw, BinningValue::binZ, pl,
481  trf);
482  draw_surfaces(surfaces, "SurfaceArrayCreator_createEquidistantAxis_Z_3.obj");
483  BOOST_CHECK_EQUAL(axis.nBins, 10u);
484  CHECK_CLOSE_ABS(axis.max, 30.9749, 1e-3);
485  CHECK_CLOSE_ABS(axis.min, -0.974873, 1e-3);
486  BOOST_CHECK_EQUAL(axis.bType, equidistant);
487 }
488 
489 BOOST_FIXTURE_TEST_CASE(SurfaceArrayCreator_createEquidistantAxis_R,
490  SurfaceArrayCreatorFixture) {
491  // single element in r
492  auto surfaces = fullPhiTestSurfacesEC(1, 0, 0, 15);
493  auto surfacesRaw = unpack_shared_vector(surfaces);
494  draw_surfaces(surfaces, "SurfaceArrayCreator_createEquidistantAxis_R_1.obj");
495  auto trf = Transform3D::Identity();
496  ProtoLayer pl = ProtoLayer(tgContext, surfacesRaw);
497  auto axis = createEquidistantAxis(tgContext, surfacesRaw, BinningValue::binR,
498  pl, trf);
499  BOOST_CHECK_EQUAL(axis.nBins, 1u);
500  CHECK_CLOSE_ABS(axis.max, perp(Vector3D(17, 1, 0)), 1e-3);
501  CHECK_CLOSE_ABS(axis.min, 13, 1e-3);
502  BOOST_CHECK_EQUAL(axis.bType, equidistant);
503 
504  // multiple rings
505  surfaces.resize(0);
506  auto ringa = fullPhiTestSurfacesEC(30, 0, 0, 10);
507  surfaces.insert(surfaces.end(), ringa.begin(), ringa.end());
508  auto ringb = fullPhiTestSurfacesEC(30, 0, 0, 15);
509  surfaces.insert(surfaces.end(), ringb.begin(), ringb.end());
510  auto ringc = fullPhiTestSurfacesEC(30, 0, 0, 20);
511  surfaces.insert(surfaces.end(), ringc.begin(), ringc.end());
512  draw_surfaces(surfaces, "SurfaceArrayCreator_createEquidistantAxis_R_2.obj");
513 
514  surfacesRaw = unpack_shared_vector(surfaces);
515  pl = ProtoLayer(tgContext, surfacesRaw);
516  trf = Transform3D::Identity();
517  axis = createEquidistantAxis(tgContext, surfacesRaw, BinningValue::binR, pl,
518  trf);
519 
520  BOOST_CHECK_EQUAL(axis.nBins, 3u);
521  CHECK_CLOSE_REL(axis.max, perp(Vector3D(20 + 2, 1, 0)), 1e-3);
522  CHECK_CLOSE_ABS(axis.min, 8, 1e-3);
523  BOOST_CHECK_EQUAL(axis.bType, equidistant);
524 }
525 
526 // if there are concentring disc or barrel modules, the bin count might be off
527 // we want to create _as few bins_ as possible, meaning the r-ring with
528 // the lowest number of surfaces should be used for the bin count or
529 // as basis for the variable edge procedure
530 // double filling will make sure no surfaces are dropped
531 BOOST_FIXTURE_TEST_CASE(SurfaceArrayCreator_dependentBinCounts,
532  SurfaceArrayCreatorFixture) {
533  auto ringA = fullPhiTestSurfacesEC(10, 0, 0, 10, 2, 3);
534  auto ringB = fullPhiTestSurfacesEC(15, 0, 0, 15, 2, 3.5);
535  auto ringC = fullPhiTestSurfacesEC(20, 0, 0, 20, 2, 3.8);
536 
537  std::vector<std::shared_ptr<const Surface>> surfaces;
538  std::copy(ringA.begin(), ringA.end(), std::back_inserter(surfaces));
539  std::copy(ringB.begin(), ringB.end(), std::back_inserter(surfaces));
540  std::copy(ringC.begin(), ringC.end(), std::back_inserter(surfaces));
541  draw_surfaces(surfaces, "SurfaceArrayCreator_dependentBinCounts.obj");
542 
543  std::unique_ptr<SurfaceArray> sArray =
544  m_SAC.surfaceArrayOnDisc(tgContext, surfaces, equidistant, equidistant);
545  auto axes = sArray->getAxes();
546  BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), 3u);
547  BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), 10u);
548 }
549 
550 BOOST_FIXTURE_TEST_CASE(SurfaceArrayCreator_completeBinning,
551  SurfaceArrayCreatorFixture) {
552  SrfVec brl = makeBarrel(30, 7, 2, 1);
553  std::vector<const Surface*> brlRaw = unpack_shared_vector(brl);
554  draw_surfaces(brl, "SurfaceArrayCreator_completeBinning_BRL.obj");
555 
556  detail::Axis<detail::AxisType::Equidistant, detail::AxisBoundaryType::Closed>
557  phiAxis(-M_PI, M_PI, 30u);
558  detail::Axis<detail::AxisType::Equidistant, detail::AxisBoundaryType::Bound>
559  zAxis(-14, 14, 7u);
560 
561  double R = 10.;
562  auto globalToLocal = [](const Vector3D& pos) {
563  return Vector2D(phi(pos) + 2 * M_PI / 30 / 2, pos.z());
564  };
565  auto localToGlobal = [R](const Vector2D& loc) {
566  double phi = loc[0] - 2 * M_PI / 30 / 2;
567  return Vector3D(R * std::cos(phi), R * std::sin(phi), loc[1]);
568  };
569 
570  auto sl = std::make_unique<
571  SurfaceArray::SurfaceGridLookup<decltype(phiAxis), decltype(zAxis)>>(
572  globalToLocal, localToGlobal,
573  std::make_tuple(std::move(phiAxis), std::move(zAxis)));
574  sl->fill(tgContext, brlRaw);
575  SurfaceArray sa(std::move(sl), brl);
576 
577  // actually filled SA
578  for (const auto& srf : brl) {
579  Vector3D ctr = srf->binningPosition(tgContext, binR);
580  auto binContent = sa.at(ctr);
581 
582  BOOST_CHECK_EQUAL(binContent.size(), 1u);
583  BOOST_CHECK_EQUAL(srf.get(), binContent.at(0));
584  }
585 }
586 
587 BOOST_FIXTURE_TEST_CASE(SurfaceArrayCreator_barrelStagger,
588  SurfaceArrayCreatorFixture) {
589  auto barrel = makeBarrelStagger(30, 7, 0, M_PI / 9.);
590  auto brl = barrel.first;
591  std::vector<const Surface*> brlRaw = unpack_shared_vector(brl);
592  draw_surfaces(brl, "SurfaceArrayCreator_barrelStagger.obj");
593 
594  ProtoLayer pl(tgContext, brl);
595 
596  // EQUIDISTANT
597  Transform3D tr = Transform3D::Identity();
598 
599  auto pAxisPhi =
600  createEquidistantAxis(tgContext, brlRaw, BinningValue::binPhi, pl, tr);
601  auto pAxisZ =
602  createEquidistantAxis(tgContext, brlRaw, BinningValue::binZ, pl, tr);
603 
604  double R = 10.;
605  Transform3D itr = tr.inverse();
606 
607  auto globalToLocal = [tr](const Vector3D& pos) {
608  Vector3D rot = tr * pos;
609  return Vector2D(phi(rot), rot.z());
610  };
611  auto localToGlobal = [R, itr](const Vector2D& loc) {
612  return itr * Vector3D(R * std::cos(loc[0]), R * std::sin(loc[0]), loc[1]);
613  };
614 
615  auto sl = makeSurfaceGridLookup2D<detail::AxisBoundaryType::Closed,
616  detail::AxisBoundaryType::Bound>(
617  globalToLocal, localToGlobal, pAxisPhi, pAxisZ);
618 
619  sl->fill(tgContext, brlRaw);
620  SurfaceArray sa(std::move(sl), brl);
621  auto axes = sa.getAxes();
622  BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), 30u);
623  BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), 7u);
624 
625  for (const auto& pr : barrel.second) {
626  auto A = pr.first;
627  auto B = pr.second;
628 
629  Vector3D ctr = A->binningPosition(tgContext, binR);
630  auto binContent = sa.at(ctr);
631  BOOST_CHECK_EQUAL(binContent.size(), 2u);
632  std::set<const Surface*> act;
633  act.insert(binContent[0]);
634  act.insert(binContent[1]);
635 
636  std::set<const Surface*> exp;
637  exp.insert(A);
638  exp.insert(B);
639 
640  BOOST_CHECK(act == exp);
641  }
642 
643  // VARIABLE
644  BOOST_TEST_CONTEXT("Barrel Stagger Variable binning") {
645  tr = Transform3D::Identity();
646 
647  auto pAxisPhiVar =
648  createVariableAxis(tgContext, brlRaw, BinningValue::binPhi, pl, tr);
649  auto pAxisZVar =
650  createVariableAxis(tgContext, brlRaw, BinningValue::binZ, pl, tr);
651 
652  itr = tr.inverse();
653 
654  auto globalToLocalVar = [tr](const Vector3D& pos) {
655  Vector3D rot = tr * pos;
656  return Vector2D(phi(rot), rot.z());
657  };
658  auto localToGlobalVar = [R, itr](const Vector2D& loc) {
659  return itr * Vector3D(R * std::cos(loc[0]), R * std::sin(loc[0]), loc[1]);
660  };
661 
662  auto sl2 = makeSurfaceGridLookup2D<detail::AxisBoundaryType::Closed,
663  detail::AxisBoundaryType::Bound>(
664  globalToLocalVar, localToGlobalVar, pAxisPhiVar, pAxisZVar);
665 
666  sl2->fill(tgContext, brlRaw);
667  SurfaceArray sa2(std::move(sl2), brl);
668  axes = sa2.getAxes();
669  BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), 30u);
670  BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), 7u);
671 
672  // check bin edges
673  std::vector<double> phiEdgesExp = {
674  -3.14159, -2.93215, -2.72271, -2.51327, -2.30383, -2.0944,
675  -1.88496, -1.67552, -1.46608, -1.25664, -1.0472, -0.837758,
676  -0.628319, -0.418879, -0.20944, 4.44089e-16, 0.20944, 0.418879,
677  0.628319, 0.837758, 1.0472, 1.25664, 1.46608, 1.67552,
678  1.88496, 2.0944, 2.30383, 2.51327, 2.72271, 3.00831,
679  3.14159};
680  std::vector<double> zEdgesExp = {-14, -10, -6, -2, 2, 6, 10, 14};
681  size_t i = 0;
682  for (const auto& edge : axes.at(0)->getBinEdges()) {
683  BOOST_TEST_INFO("phi edge index " << i);
684  auto phiEdge = phiEdgesExp.at(i);
685  CHECK_CLOSE_ABS(edge, phiEdge, 1e-5 * M_PI);
686  i++;
687  }
688  i = 0;
689  for (const auto& edge : axes.at(1)->getBinEdges()) {
690  BOOST_TEST_INFO("z edge index " << i);
691  CHECK_CLOSE_ABS(edge, zEdgesExp.at(i), 1e-6);
692  i++;
693  }
694 
695  for (const auto& pr : barrel.second) {
696  auto A = pr.first;
697  auto B = pr.second;
698 
699  Vector3D ctr = A->binningPosition(tgContext, binR);
700  auto binContent = sa2.at(ctr);
701  BOOST_CHECK_EQUAL(binContent.size(), 2u);
702  std::set<const Surface*> act;
703  act.insert(binContent[0]);
704  act.insert(binContent[1]);
705 
706  std::set<const Surface*> exp;
707  exp.insert(A);
708  exp.insert(B);
709 
710  BOOST_CHECK(act == exp);
711  }
712  }
713 }
714 
715 BOOST_AUTO_TEST_SUITE_END()
716 } // namespace Test
717 
718 } // namespace Acts