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