ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LayerCreatorTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LayerCreatorTests.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/CylinderLayer.hpp"
17 #include "Acts/Geometry/DiscLayer.hpp"
18 #include "Acts/Geometry/GeometryContext.hpp"
19 #include "Acts/Geometry/LayerCreator.hpp"
20 #include "Acts/Geometry/ProtoLayer.hpp"
21 #include "Acts/Geometry/SurfaceArrayCreator.hpp"
22 #include "Acts/Surfaces/CylinderBounds.hpp"
23 #include "Acts/Surfaces/PlaneSurface.hpp"
24 #include "Acts/Surfaces/RadialBounds.hpp"
25 #include "Acts/Surfaces/RectangleBounds.hpp"
26 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
27 #include "Acts/Utilities/BinningType.hpp"
28 #include "Acts/Utilities/Definitions.hpp"
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(VectorHelpers::phi(v), (a), tolerance); \
44  }
45 
46 using SrfVec = std::vector<std::shared_ptr<const Surface>>;
47 
48 void draw_surfaces(const SrfVec& surfaces, const std::string& fname) {
49  std::ofstream os;
50  os.open(fname);
51 
52  os << std::fixed << std::setprecision(4);
53 
54  size_t nVtx = 0;
55  for (const auto& srfx : surfaces) {
56  std::shared_ptr<const PlaneSurface> srf =
57  std::dynamic_pointer_cast<const PlaneSurface>(srfx);
58  const PlanarBounds* bounds =
59  dynamic_cast<const PlanarBounds*>(&srf->bounds());
60 
61  for (const auto& vtxloc : bounds->vertices()) {
62  Vector3D vtx =
63  srf->transform(tgContext) * Vector3D(vtxloc.x(), vtxloc.y(), 0);
64  os << "v " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
65  }
66 
67  // connect them
68  os << "f";
69  for (size_t i = 1; i <= bounds->vertices().size(); ++i) {
70  os << " " << nVtx + i;
71  }
72  os << "\n";
73 
74  nVtx += bounds->vertices().size();
75  }
76 
77  os.close();
78 }
79 
80 struct LayerCreatorFixture {
81  std::shared_ptr<const SurfaceArrayCreator> p_SAC;
82  std::shared_ptr<LayerCreator> p_LC;
83  std::vector<std::shared_ptr<const Surface>> m_surfaces;
84 
85  LayerCreatorFixture() {
86  p_SAC = std::make_shared<const SurfaceArrayCreator>(
87  SurfaceArrayCreator::Config(),
88  Acts::getDefaultLogger("SurfaceArrayCreator", Acts::Logging::VERBOSE));
89  LayerCreator::Config cfg;
90  cfg.surfaceArrayCreator = p_SAC;
91  p_LC = std::make_shared<LayerCreator>(
92  cfg, Acts::getDefaultLogger("LayerCreator", Acts::Logging::VERBOSE));
93  }
94 
95  template <typename... Args>
96  bool checkBinning(Args&&... args) {
97  return p_LC->checkBinning(std::forward<Args>(args)...);
98  }
99 
100  bool checkBinContentSize(const SurfaceArray* sArray, size_t n) {
101  size_t nBins = sArray->size();
102  bool result = true;
103  for (size_t i = 0; i < nBins; ++i) {
104  if (!sArray->isValidBin(i)) {
105  continue;
106  }
107  std::vector<const Surface*> binContent = sArray->at(i);
108  BOOST_TEST_INFO("Bin: " << i);
109  BOOST_CHECK_EQUAL(binContent.size(), n);
110  result = result && binContent.size() == n;
111  }
112 
113  return result;
114  }
115 
116  SrfVec fullPhiTestSurfacesEC(size_t n = 10, double shift = 0,
117  double zbase = 0, double r = 10) {
118  SrfVec res;
119 
120  double phiStep = 2 * M_PI / n;
121  for (size_t i = 0; i < n; ++i) {
122  double z = zbase + ((i % 2 == 0) ? 1 : -1) * 0.2;
123 
124  Transform3D trans;
125  trans.setIdentity();
126  trans.rotate(Eigen::AngleAxisd(i * phiStep + shift, Vector3D(0, 0, 1)));
127  trans.translate(Vector3D(r, 0, z));
128 
129  auto bounds = std::make_shared<const RectangleBounds>(2, 1);
130 
131  auto transptr = std::make_shared<const Transform3D>(trans);
132  std::shared_ptr<PlaneSurface> srf =
133  Surface::makeShared<PlaneSurface>(transptr, bounds);
134 
135  res.push_back(srf);
136  m_surfaces.push_back(
137  std::move(srf)); // keep shared, will get destroyed at the end
138  }
139 
140  return res;
141  }
142 
143  SrfVec fullPhiTestSurfacesBRL(int n = 10, double shift = 0, double zbase = 0,
144  double incl = M_PI / 9., double w = 2,
145  double h = 1.5) {
146  SrfVec res;
147 
148  double phiStep = 2 * M_PI / n;
149  for (int i = 0; i < n; ++i) {
150  double z = zbase;
151 
152  Transform3D trans;
153  trans.setIdentity();
154  trans.rotate(Eigen::AngleAxisd(i * phiStep + shift, Vector3D(0, 0, 1)));
155  trans.translate(Vector3D(10, 0, z));
156  trans.rotate(Eigen::AngleAxisd(incl, Vector3D(0, 0, 1)));
157  trans.rotate(Eigen::AngleAxisd(M_PI / 2., Vector3D(0, 1, 0)));
158 
159  auto bounds = std::make_shared<const RectangleBounds>(w, h);
160 
161  auto transptr = std::make_shared<const Transform3D>(trans);
162  std::shared_ptr<PlaneSurface> srf =
163  Surface::makeShared<PlaneSurface>(transptr, bounds);
164 
165  res.push_back(srf);
166  m_surfaces.push_back(
167  std::move(srf)); // keep shared, will get destroyed at the end
168  }
169 
170  return res;
171  }
172 
173  SrfVec makeBarrel(int nPhi, int nZ, double w, double h) {
174  double z0 = -(nZ - 1) * w;
175  SrfVec res;
176 
177  for (int i = 0; i < nZ; i++) {
178  double z = i * w * 2 + z0;
179  std::cout << "z=" << z << std::endl;
180  SrfVec ring = fullPhiTestSurfacesBRL(nPhi, 0, z, M_PI / 9., w, h);
181  res.insert(res.end(), ring.begin(), ring.end());
182  }
183 
184  return res;
185  }
186 
187  std::pair<SrfVec, std::vector<std::pair<const Surface*, const Surface*>>>
188  makeBarrelStagger(int nPhi, int nZ, double shift = 0, double incl = M_PI / 9.,
189  double w = 2, double h = 1.5) {
190  double z0 = -(nZ - 1) * w;
191  SrfVec res;
192 
193  std::vector<std::pair<const Surface*, const Surface*>> pairs;
194 
195  for (int i = 0; i < nZ; i++) {
196  double z = i * w * 2 + z0;
197 
198  double phiStep = 2 * M_PI / nPhi;
199  for (int j = 0; j < nPhi; ++j) {
200  Transform3D trans;
201  trans.setIdentity();
202  trans.rotate(Eigen::AngleAxisd(j * phiStep + shift, Vector3D(0, 0, 1)));
203  trans.translate(Vector3D(10, 0, z));
204  trans.rotate(Eigen::AngleAxisd(incl, Vector3D(0, 0, 1)));
205  trans.rotate(Eigen::AngleAxisd(M_PI / 2., Vector3D(0, 1, 0)));
206 
207  auto bounds = std::make_shared<const RectangleBounds>(w, h);
208 
209  auto transAptr = std::make_shared<const Transform3D>(trans);
210 
211  std::shared_ptr<PlaneSurface> srfA =
212  Surface::makeShared<PlaneSurface>(transAptr, bounds);
213 
214  Vector3D nrm = srfA->normal(tgContext);
215  Transform3D transB = trans;
216  transB.pretranslate(nrm * 0.1);
217  auto transBptr = std::make_shared<const Transform3D>(transB);
218  std::shared_ptr<PlaneSurface> srfB =
219  Surface::makeShared<PlaneSurface>(transBptr, bounds);
220 
221  pairs.push_back(std::make_pair(srfA.get(), srfB.get()));
222 
223  res.push_back(srfA);
224  res.push_back(srfB);
225  m_surfaces.push_back(std::move(srfA));
226  m_surfaces.push_back(std::move(srfB));
227  }
228  }
229 
230  return std::make_pair(res, pairs);
231  }
232 };
233 
234 BOOST_AUTO_TEST_SUITE(Tools)
235 
236 BOOST_FIXTURE_TEST_CASE(LayerCreator_createCylinderLayer, LayerCreatorFixture) {
237  std::vector<std::shared_ptr<const Surface>> srf;
238 
239  srf = makeBarrel(30, 7, 2, 1.5);
240  draw_surfaces(srf, "LayerCreator_createCylinderLayer_BRL_1.obj");
241 
242  // CASE I
243  double envR = 0.1, envZ = 0.5;
244  ProtoLayer pl(tgContext, srf);
245  pl.envR = {envR, envR};
246  pl.envZ = {envZ, envZ};
247  std::shared_ptr<CylinderLayer> layer =
248  std::dynamic_pointer_cast<CylinderLayer>(
249  p_LC->cylinderLayer(tgContext, srf, equidistant, equidistant, pl));
250 
251  double rMax = 10.6071, rMin = 9.59111; // empirical
252  CHECK_CLOSE_REL(layer->thickness(), (rMax - rMin) + 2 * envR, 1e-3);
253 
254  const CylinderBounds* bounds = &layer->bounds();
255  CHECK_CLOSE_REL(bounds->get(CylinderBounds::eR), (rMax + rMin) / 2., 1e-3);
256  CHECK_CLOSE_REL(bounds->get(CylinderBounds::eHalfLengthZ), 14 + envZ, 1e-3);
257  BOOST_CHECK(checkBinning(tgContext, *layer->surfaceArray()));
258  auto axes = layer->surfaceArray()->getAxes();
259  BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), 30u);
260  BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), 7u);
261  CHECK_CLOSE_REL(axes.at(0)->getMin(), -M_PI, 1e-3);
262  CHECK_CLOSE_REL(axes.at(0)->getMax(), M_PI, 1e-3);
263  CHECK_CLOSE_REL(axes.at(1)->getMin(), -14, 1e-3);
264  CHECK_CLOSE_REL(axes.at(1)->getMax(), 14, 1e-3);
265 
266  // CASE II
267 
268  ProtoLayer pl2(tgContext, srf);
269  pl2.envR = {envR, envR};
270  pl2.envZ = {envZ, envZ};
271  layer = std::dynamic_pointer_cast<CylinderLayer>(
272  p_LC->cylinderLayer(tgContext, srf, 30, 7, pl2));
273  CHECK_CLOSE_REL(layer->thickness(), (rMax - rMin) + 2 * envR, 1e-3);
274  bounds = &layer->bounds();
275  CHECK_CLOSE_REL(bounds->get(CylinderBounds::eR), (rMax + rMin) / 2., 1e-3);
276  CHECK_CLOSE_REL(bounds->get(CylinderBounds::eHalfLengthZ), 14 + envZ, 1e-3);
277  BOOST_CHECK(checkBinning(tgContext, *layer->surfaceArray()));
278  axes = layer->surfaceArray()->getAxes();
279  BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), 30u);
280  BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), 7u);
281  CHECK_CLOSE_REL(axes.at(0)->getMin(), -M_PI, 1e-3);
282  CHECK_CLOSE_REL(axes.at(0)->getMax(), M_PI, 1e-3);
283  CHECK_CLOSE_REL(axes.at(1)->getMin(), -14, 1e-3);
284  CHECK_CLOSE_REL(axes.at(1)->getMax(), 14, 1e-3);
285 
286  layer = std::dynamic_pointer_cast<CylinderLayer>(
287  p_LC->cylinderLayer(tgContext, srf, 13, 3, pl2));
288  CHECK_CLOSE_REL(layer->thickness(), (rMax - rMin) + 2 * envR, 1e-3);
289  bounds = &layer->bounds();
290  CHECK_CLOSE_REL(bounds->get(CylinderBounds::eR), (rMax + rMin) / 2., 1e-3);
291  CHECK_CLOSE_REL(bounds->get(CylinderBounds::eHalfLengthZ), 14 + envZ, 1e-3);
292  // this succeeds despite sub-optimal binning
293  // since we now have multientry bins
294  BOOST_CHECK(checkBinning(tgContext, *layer->surfaceArray()));
295  axes = layer->surfaceArray()->getAxes();
296  BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), 13u);
297  BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), 3u);
298  CHECK_CLOSE_REL(axes.at(0)->getMin(), -M_PI, 1e-3);
299  CHECK_CLOSE_REL(axes.at(0)->getMax(), M_PI, 1e-3);
300  CHECK_CLOSE_REL(axes.at(1)->getMin(), -14, 1e-3);
301  CHECK_CLOSE_REL(axes.at(1)->getMax(), 14, 1e-3);
302 
303  // CASE III
304  ProtoLayer pl3;
305  pl3.minR = 1;
306  pl3.maxR = 20;
307  pl3.minZ = -25;
308  pl3.maxZ = 25;
309  layer = std::dynamic_pointer_cast<CylinderLayer>(
310  p_LC->cylinderLayer(tgContext, srf, equidistant, equidistant, pl3));
311  CHECK_CLOSE_REL(layer->thickness(), 19, 1e-3);
312  bounds = &layer->bounds();
313  CHECK_CLOSE_REL(bounds->get(CylinderBounds::eR), 10.5, 1e-3);
314  CHECK_CLOSE_REL(bounds->get(CylinderBounds::eHalfLengthZ), 25, 1e-3);
315 
316  // this should fail, b/c it's a completely inconvenient binning
317  // but it succeeds despite sub-optimal binning
318  // since we now have multientry bins
319  BOOST_CHECK(checkBinning(tgContext, *layer->surfaceArray()));
320 
321  axes = layer->surfaceArray()->getAxes();
322  BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), 30u);
323  BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), 7u);
324  CHECK_CLOSE_REL(axes.at(0)->getMin(), -M_PI, 1e-3);
325  CHECK_CLOSE_REL(axes.at(0)->getMax(), M_PI, 1e-3);
326  CHECK_CLOSE_REL(axes.at(1)->getMin(), -25, 1e-3);
327  CHECK_CLOSE_REL(axes.at(1)->getMax(), 25, 1e-3);
328 }
329 
330 BOOST_FIXTURE_TEST_CASE(LayerCreator_createDiscLayer, LayerCreatorFixture) {
331  std::vector<std::shared_ptr<const Surface>> surfaces;
332  auto ringa = fullPhiTestSurfacesEC(30, 0, 0, 10);
333  surfaces.insert(surfaces.end(), ringa.begin(), ringa.end());
334  auto ringb = fullPhiTestSurfacesEC(30, 0, 0, 15);
335  surfaces.insert(surfaces.end(), ringb.begin(), ringb.end());
336  auto ringc = fullPhiTestSurfacesEC(30, 0, 0, 20);
337  surfaces.insert(surfaces.end(), ringc.begin(), ringc.end());
338  draw_surfaces(surfaces, "LayerCreator_createDiscLayer_EC_1.obj");
339 
340  ProtoLayer pl(tgContext, surfaces);
341  pl.minZ = -10;
342  pl.maxZ = 10;
343  pl.minR = 5;
344  pl.maxR = 25;
345  std::shared_ptr<DiscLayer> layer = std::dynamic_pointer_cast<DiscLayer>(
346  p_LC->discLayer(tgContext, surfaces, equidistant, equidistant, pl));
347  CHECK_CLOSE_REL(layer->thickness(), 20, 1e-3);
348  const RadialBounds* bounds =
349  dynamic_cast<const RadialBounds*>(&layer->bounds());
350  CHECK_CLOSE_REL(bounds->rMin(), 5, 1e-3);
351  CHECK_CLOSE_REL(bounds->rMax(), 25, 1e-3);
352  BOOST_CHECK(checkBinning(tgContext, *layer->surfaceArray()));
353  auto axes = layer->surfaceArray()->getAxes();
354  BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), 3u);
355  BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), 30u);
356  CHECK_CLOSE_REL(axes.at(0)->getMin(), 5, 1e-3);
357  CHECK_CLOSE_REL(axes.at(0)->getMax(), 25, 1e-3);
358  CHECK_CLOSE_REL(axes.at(1)->getMin(), -M_PI, 1e-3);
359  CHECK_CLOSE_REL(axes.at(1)->getMax(), M_PI, 1e-3);
360  checkBinContentSize(layer->surfaceArray(), 1);
361 
362  // check that it's applying a rotation transform to improve phi binning
363  // BOOST_CHECK_NE(bu->transform(), nullptr);
364  // double actAngle = ((*bu->transform()) * Vector3D(1, 0, 0)).phi();
365  // double expAngle = -2 * M_PI / 30 / 2.;
366  // CHECK_CLOSE_REL(actAngle, expAngle, 1e-3);
367 
368  double envMinR = 1, envMaxR = 1, envZ = 5;
369  size_t nBinsR = 3, nBinsPhi = 30;
370  ProtoLayer pl2(tgContext, surfaces);
371  pl2.envR = {envMinR, envMaxR};
372  pl2.envZ = {envZ, envZ};
373  layer = std::dynamic_pointer_cast<DiscLayer>(
374  p_LC->discLayer(tgContext, surfaces, nBinsR, nBinsPhi, pl2));
375 
376  double rMin = 8, rMax = 22.0227;
377  CHECK_CLOSE_REL(layer->thickness(), 0.4 + 2 * envZ, 1e-3);
378  bounds = dynamic_cast<const RadialBounds*>(&layer->bounds());
379  CHECK_CLOSE_REL(bounds->rMin(), rMin - envMinR, 1e-3);
380  CHECK_CLOSE_REL(bounds->rMax(), rMax + envMaxR, 1e-3);
381  BOOST_CHECK(checkBinning(tgContext, *layer->surfaceArray()));
382  axes = layer->surfaceArray()->getAxes();
383  BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), nBinsR);
384  BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), nBinsPhi);
385  CHECK_CLOSE_REL(axes.at(0)->getMin(), rMin, 1e-3);
386  CHECK_CLOSE_REL(axes.at(0)->getMax(), rMax, 1e-3);
387  CHECK_CLOSE_REL(axes.at(1)->getMin(), -M_PI, 1e-3);
388  CHECK_CLOSE_REL(axes.at(1)->getMax(), M_PI, 1e-3);
389  checkBinContentSize(layer->surfaceArray(), 1);
390 
391  // check that it's applying a rotation transform to improve phi binning
392  // BOOST_CHECK_NE(bu->transform(), nullptr);
393  // actAngle = ((*bu->transform()) * Vector3D(1, 0, 0)).phi();
394  // expAngle = -2 * M_PI / 30 / 2.;
395  // CHECK_CLOSE_REL(actAngle, expAngle, 1e-3);
396 
397  layer = std::dynamic_pointer_cast<DiscLayer>(
398  p_LC->discLayer(tgContext, surfaces, equidistant, equidistant, pl2));
399  CHECK_CLOSE_REL(layer->thickness(), 0.4 + 2 * envZ, 1e-3);
400  bounds = dynamic_cast<const RadialBounds*>(&layer->bounds());
401  CHECK_CLOSE_REL(bounds->rMin(), rMin - envMinR, 1e-3);
402  CHECK_CLOSE_REL(bounds->rMax(), rMax + envMaxR, 1e-3);
403  BOOST_CHECK(checkBinning(tgContext, *layer->surfaceArray()));
404  axes = layer->surfaceArray()->getAxes();
405  BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), nBinsR);
406  BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), nBinsPhi);
407  CHECK_CLOSE_REL(axes.at(0)->getMin(), rMin, 1e-3);
408  CHECK_CLOSE_REL(axes.at(0)->getMax(), rMax, 1e-3);
409  CHECK_CLOSE_REL(axes.at(1)->getMin(), -M_PI, 1e-3);
410  CHECK_CLOSE_REL(axes.at(1)->getMax(), M_PI, 1e-3);
411  checkBinContentSize(layer->surfaceArray(), 1);
412 
413  // check that it's applying a rotation transform to improve phi binning
414  // BOOST_CHECK_NE(bu->transform(), nullptr);
415  // actAngle = ((*bu->transform()) * Vector3D(1, 0, 0)).phi();
416  // expAngle = -2 * M_PI / 30 / 2.;
417  // CHECK_CLOSE_REL(actAngle, expAngle, 1e-3);
418 }
419 
420 BOOST_FIXTURE_TEST_CASE(LayerCreator_barrelStagger, LayerCreatorFixture) {
421  auto barrel = makeBarrelStagger(30, 7, 0, M_PI / 9.);
422  auto brl = barrel.first;
423  draw_surfaces(brl, "LayerCreator_barrelStagger.obj");
424 
425  double envR = 0, envZ = 0;
426  ProtoLayer pl(tgContext, brl);
427  pl.envR = {envR, envR};
428  pl.envZ = {envZ, envZ};
429  std::shared_ptr<CylinderLayer> layer =
430  std::dynamic_pointer_cast<CylinderLayer>(
431  p_LC->cylinderLayer(tgContext, brl, equidistant, equidistant, pl));
432 
433  auto axes = layer->surfaceArray()->getAxes();
434  BOOST_CHECK_EQUAL(axes.at(0)->getNBins(), 30u);
435  BOOST_CHECK_EQUAL(axes.at(1)->getNBins(), 7u);
436 
437  // check if binning is good!
438  for (const auto& pr : barrel.second) {
439  auto A = pr.first;
440  auto B = pr.second;
441 
442  // std::cout << A->center().phi() << " ";
443  // std::cout << B->center().phi() << std::endl;
444  // std::cout << "dPHi = " << A->center().phi() - B->center().phi() <<
445  // std::endl;
446 
447  Vector3D ctr = A->binningPosition(tgContext, binR);
448  auto binContent = layer->surfaceArray()->at(ctr);
449  BOOST_CHECK_EQUAL(binContent.size(), 2u);
450  std::set<const Surface*> act;
451  act.insert(binContent[0]);
452  act.insert(binContent[1]);
453 
454  std::set<const Surface*> exp;
455  exp.insert(A);
456  exp.insert(B);
457  BOOST_CHECK(exp == act);
458  }
459 
460  // checkBinning should also report everything is fine
461  checkBinning(tgContext, *layer->surfaceArray());
462 }
463 
464 BOOST_AUTO_TEST_SUITE_END()
465 } // namespace Test
466 
467 } // namespace Acts