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