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