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