Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LayerCreator.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file LayerCreator.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2016-2020 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 "Acts/Geometry/LayerCreator.hpp"
10 
11 #include "Acts/Definitions/Algebra.hpp"
12 #include "Acts/Geometry/CylinderLayer.hpp"
13 #include "Acts/Geometry/DiscLayer.hpp"
14 #include "Acts/Geometry/Extent.hpp"
15 #include "Acts/Geometry/Layer.hpp"
16 #include "Acts/Geometry/PlaneLayer.hpp"
17 #include "Acts/Geometry/ProtoLayer.hpp"
18 #include "Acts/Geometry/SurfaceArrayCreator.hpp"
19 #include "Acts/Surfaces/CylinderBounds.hpp"
20 #include "Acts/Surfaces/RadialBounds.hpp"
21 #include "Acts/Surfaces/RectangleBounds.hpp"
22 #include "Acts/Surfaces/Surface.hpp"
23 
24 #include <algorithm>
25 #include <array>
26 #include <iterator>
27 #include <ostream>
28 #include <set>
29 #include <utility>
30 
31 namespace Acts {
32 class PlanarBounds;
33 } // namespace Acts
34 
35 using Acts::VectorHelpers::perp;
36 using Acts::VectorHelpers::phi;
37 
38 Acts::LayerCreator::LayerCreator(const Acts::LayerCreator::Config& lcConfig,
39  std::unique_ptr<const Logger> logger)
40  : m_cfg(lcConfig), m_logger(std::move(logger)) {}
41 
42 void Acts::LayerCreator::setConfiguration(
43  const Acts::LayerCreator::Config& lcConfig) {
44  // @todo check consistency
45  // copy the configuration
46  m_cfg = lcConfig;
47 }
48 
49 void Acts::LayerCreator::setLogger(std::unique_ptr<const Logger> newLogger) {
50  m_logger = std::move(newLogger);
51 }
52 
53 Acts::MutableLayerPtr Acts::LayerCreator::cylinderLayer(
54  const GeometryContext& gctx,
55  std::vector<std::shared_ptr<const Surface>> surfaces, size_t binsPhi,
56  size_t binsZ, std::optional<ProtoLayer> _protoLayer,
57  const Transform3& transform, std::unique_ptr<ApproachDescriptor> ad) const {
58  ProtoLayer protoLayer =
59  _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
60 
61  // Remaining layer parameters - they include the envelopes
62  double layerR = protoLayer.medium(binR);
63  double layerZ = protoLayer.medium(binZ);
64  double layerHalfZ = 0.5 * protoLayer.range(binZ);
65  double layerThickness = protoLayer.range(binR);
66 
67  ACTS_VERBOSE("Creating a cylindrical Layer:");
68  ACTS_VERBOSE(" - with layer R = " << layerR);
69  ACTS_VERBOSE(" - from R min/max = " << protoLayer.min(binR, false) << " / "
70  << protoLayer.max(binR, false));
71  ACTS_VERBOSE(" - with R thickness = " << layerThickness);
72  ACTS_VERBOSE(" - incl envelope = " << protoLayer.envelope[binR][0u]
73  << " / "
74  << protoLayer.envelope[binR][1u]);
75 
76  ACTS_VERBOSE(" - with z min/max = "
77  << protoLayer.min(binZ, false) << " (-"
78  << protoLayer.envelope[binZ][0u] << ") / "
79  << protoLayer.max(binZ, false) << " (+"
80  << protoLayer.envelope[binZ][1u] << ")");
81 
82  ACTS_VERBOSE(" - z center = " << layerZ);
83  ACTS_VERBOSE(" - halflength z = " << layerHalfZ);
84 
85  // create the layer transforms if not given
86  // we need to transform in case layerZ != 0, so that the layer will be
87  // correctly defined using the halflength
88  Translation3 addTranslation(0., 0., 0.);
89  if (transform.isApprox(Transform3::Identity())) {
90  addTranslation = Translation3(0., 0., layerZ);
91  ACTS_VERBOSE(" - layer z shift = " << -layerZ);
92  }
93 
94  ACTS_VERBOSE(" - with phi min/max = " << protoLayer.min(binPhi, false)
95  << " / "
96  << protoLayer.max(binPhi, false));
97  ACTS_VERBOSE(" - # of modules = " << surfaces.size() << " ordered in ( "
98  << binsPhi << " x " << binsZ << ")");
99  std::unique_ptr<SurfaceArray> sArray;
100  if (!surfaces.empty()) {
101  sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnCylinder(
102  gctx, std::move(surfaces), binsPhi, binsZ, protoLayer);
103 
104  checkBinning(gctx, *sArray);
105  }
106 
107  // create the layer and push it back
108  std::shared_ptr<const CylinderBounds> cBounds(
109  new CylinderBounds(layerR, layerHalfZ));
110 
111  // create the layer
112  MutableLayerPtr cLayer = CylinderLayer::create(
113  addTranslation * transform, cBounds, std::move(sArray), layerThickness,
114  std::move(ad), active);
115 
116  if (!cLayer) {
117  ACTS_ERROR("Creation of cylinder layer did not succeed!");
118  }
119  associateSurfacesToLayer(*cLayer);
120 
121  // now return
122  return cLayer;
123 }
124 
125 Acts::MutableLayerPtr Acts::LayerCreator::cylinderLayer(
126  const GeometryContext& gctx,
127  std::vector<std::shared_ptr<const Surface>> surfaces, BinningType bTypePhi,
128  BinningType bTypeZ, std::optional<ProtoLayer> _protoLayer,
129  const Transform3& transform, std::unique_ptr<ApproachDescriptor> ad) const {
130  ProtoLayer protoLayer =
131  _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
132 
133  // remaining layer parameters
134  double layerR = protoLayer.medium(binR);
135  double layerZ = protoLayer.medium(binZ);
136  double layerHalfZ = 0.5 * protoLayer.range(binZ);
137  double layerThickness = protoLayer.range(binR);
138 
139  // adjust the layer radius
140  ACTS_VERBOSE("Creating a cylindrical Layer:");
141  ACTS_VERBOSE(" - with layer R = " << layerR);
142  ACTS_VERBOSE(" - from R min/max = " << protoLayer.min(binR, false) << " / "
143  << protoLayer.max(binR, false));
144  ACTS_VERBOSE(" - with R thickness = " << layerThickness);
145  ACTS_VERBOSE(" - incl envelope = " << protoLayer.envelope[binR][0u]
146  << " / "
147  << protoLayer.envelope[binR][1u]);
148  ACTS_VERBOSE(" - with z min/max = "
149  << protoLayer.min(binZ, false) << " (-"
150  << protoLayer.envelope[binZ][0u] << ") / "
151  << protoLayer.max(binZ, false) << " (+"
152  << protoLayer.envelope[binZ][1u] << ")");
153  ACTS_VERBOSE(" - z center = " << layerZ);
154  ACTS_VERBOSE(" - halflength z = " << layerHalfZ);
155 
156  // create the layer transforms if not given
157  // we need to transform in case layerZ != 0, so that the layer will be
158  // correctly defined using the halflength
159  // create the layer transforms if not given
160  Translation3 addTranslation(0., 0., 0.);
161  if (transform.isApprox(Transform3::Identity()) && bTypeZ == equidistant) {
162  addTranslation = Translation3(0., 0., layerZ);
163  ACTS_VERBOSE(" - layer z shift = " << -layerZ);
164  }
165 
166  ACTS_VERBOSE(" - with phi min/max = " << protoLayer.min(binPhi, false)
167  << " / "
168  << protoLayer.max(binPhi, false));
169  ACTS_VERBOSE(" - # of modules = " << surfaces.size() << "");
170 
171  // create the surface array
172  std::unique_ptr<SurfaceArray> sArray;
173  if (!surfaces.empty()) {
174  sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnCylinder(
175  gctx, std::move(surfaces), bTypePhi, bTypeZ, protoLayer);
176 
177  checkBinning(gctx, *sArray);
178  }
179 
180  // create the layer and push it back
181  std::shared_ptr<const CylinderBounds> cBounds(
182  new CylinderBounds(layerR, layerHalfZ));
183 
184  // create the layer
185  MutableLayerPtr cLayer = CylinderLayer::create(
186  addTranslation * transform, cBounds, std::move(sArray), layerThickness,
187  std::move(ad), active);
188 
189  if (!cLayer) {
190  ACTS_ERROR("Creation of cylinder layer did not succeed!");
191  }
192  associateSurfacesToLayer(*cLayer);
193 
194  // now return
195  return cLayer;
196 }
197 
198 Acts::MutableLayerPtr Acts::LayerCreator::discLayer(
199  const GeometryContext& gctx,
200  std::vector<std::shared_ptr<const Surface>> surfaces, size_t binsR,
201  size_t binsPhi, std::optional<ProtoLayer> _protoLayer,
202  const Transform3& transform, std::unique_ptr<ApproachDescriptor> ad) const {
203  ProtoLayer protoLayer =
204  _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
205 
206  double layerZ = protoLayer.medium(binZ);
207  double layerThickness = protoLayer.range(binZ);
208 
209  // adjust the layer radius
210  ACTS_VERBOSE("Creating a disk Layer:");
211  ACTS_VERBOSE(" - at Z position = " << layerZ);
212  ACTS_VERBOSE(" - from Z min/max = " << protoLayer.min(binZ, false) << " / "
213  << protoLayer.max(binZ, false));
214  ACTS_VERBOSE(" - with Z thickness = " << layerThickness);
215  ACTS_VERBOSE(" - incl envelope = " << protoLayer.envelope[binZ][0u]
216  << " / "
217  << protoLayer.envelope[binZ][1u]);
218  ACTS_VERBOSE(" - with R min/max = "
219  << protoLayer.min(binR, false) << " (-"
220  << protoLayer.envelope[binR][0u] << ") / "
221  << protoLayer.max(binR, false) << " (+"
222  << protoLayer.envelope[binR][1u] << ")");
223  ACTS_VERBOSE(" - with phi min/max = " << protoLayer.min(binPhi, false)
224  << " / "
225  << protoLayer.max(binPhi, false));
226  ACTS_VERBOSE(" - # of modules = " << surfaces.size() << " ordered in ( "
227  << binsR << " x " << binsPhi << ")");
228 
229  // create the layer transforms if not given
230  Translation3 addTranslation(0., 0., 0.);
231  if (transform.isApprox(Transform3::Identity())) {
232  addTranslation = Translation3(0., 0., layerZ);
233  }
234  // create the surface array
235  std::unique_ptr<SurfaceArray> sArray;
236  if (!surfaces.empty()) {
237  sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnDisc(
238  gctx, std::move(surfaces), binsR, binsPhi, protoLayer, transform);
239 
240  checkBinning(gctx, *sArray);
241  }
242 
243  // create the share disc bounds
244  auto dBounds = std::make_shared<const RadialBounds>(protoLayer.min(binR),
245  protoLayer.max(binR));
246 
247  // create the layers
248  // we use the same transform here as for the layer itself
249  // for disk this is fine since we don't bin in Z, so does not matter
250  MutableLayerPtr dLayer =
251  DiscLayer::create(addTranslation * transform, dBounds, std::move(sArray),
252  layerThickness, std::move(ad), active);
253 
254  if (!dLayer) {
255  ACTS_ERROR("Creation of disc layer did not succeed!");
256  }
257  associateSurfacesToLayer(*dLayer);
258  // return the layer
259  return dLayer;
260 }
261 
262 Acts::MutableLayerPtr Acts::LayerCreator::discLayer(
263  const GeometryContext& gctx,
264  std::vector<std::shared_ptr<const Surface>> surfaces, BinningType bTypeR,
265  BinningType bTypePhi, std::optional<ProtoLayer> _protoLayer,
266  const Transform3& transform, std::unique_ptr<ApproachDescriptor> ad) const {
267  ProtoLayer protoLayer =
268  _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
269 
270  double layerZ = protoLayer.medium(binZ);
271  double layerThickness = protoLayer.range(binZ);
272 
273  // adjust the layer radius
274  ACTS_VERBOSE("Creating a disk Layer:");
275  ACTS_VERBOSE(" - at Z position = " << layerZ);
276  ACTS_VERBOSE(" - from Z min/max = " << protoLayer.min(binZ, false) << " / "
277  << protoLayer.max(binZ, false));
278  ACTS_VERBOSE(" - with Z thickness = " << layerThickness);
279  ACTS_VERBOSE(" - incl envelope = " << protoLayer.envelope[binZ][0u]
280  << " / "
281  << protoLayer.envelope[binZ][1u]);
282  ACTS_VERBOSE(" - with R min/max = "
283  << protoLayer.min(binR, false) << " (-"
284  << protoLayer.envelope[binR][0u] << ") / "
285  << protoLayer.max(binR, false) << " (+"
286  << protoLayer.envelope[binR][1u] << ")");
287  ACTS_VERBOSE(" - with phi min/max = " << protoLayer.min(binPhi, false)
288  << " / "
289  << protoLayer.max(binPhi, false));
290  ACTS_VERBOSE(" - # of modules = " << surfaces.size());
291 
292  // create the layer transforms if not given
293  Translation3 addTranslation(0., 0., 0.);
294  if (transform.isApprox(Transform3::Identity())) {
295  addTranslation = Translation3(0., 0., layerZ);
296  }
297 
298  // create the surface array
299  std::unique_ptr<SurfaceArray> sArray;
300  if (!surfaces.empty()) {
301  sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnDisc(
302  gctx, std::move(surfaces), bTypeR, bTypePhi, protoLayer, transform);
303 
304  checkBinning(gctx, *sArray);
305  }
306 
307  // create the shared disc bounds
308  auto dBounds = std::make_shared<const RadialBounds>(protoLayer.min(binR),
309  protoLayer.max(binR));
310 
311  // create the layers
312  MutableLayerPtr dLayer =
313  DiscLayer::create(addTranslation * transform, dBounds, std::move(sArray),
314  layerThickness, std::move(ad), active);
315  if (!dLayer) {
316  ACTS_ERROR("Creation of disc layer did not succeed!");
317  }
318  associateSurfacesToLayer(*dLayer);
319  // return the layer
320  return dLayer;
321 }
322 
323 Acts::MutableLayerPtr Acts::LayerCreator::planeLayer(
324  const GeometryContext& gctx,
325  std::vector<std::shared_ptr<const Surface>> surfaces, size_t bins1,
326  size_t bins2, BinningValue bValue, std::optional<ProtoLayer> _protoLayer,
327  const Transform3& transform, std::unique_ptr<ApproachDescriptor> ad) const {
328  ProtoLayer protoLayer =
329  _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
330 
331  // remaining layer parameters
332  double layerHalf1 = 0, layerHalf2 = 0, layerThickness = 0;
333  switch (bValue) {
334  case BinningValue::binX: {
335  layerHalf1 = 0.5 * (protoLayer.max(binY) - protoLayer.min(binY));
336  layerHalf2 = 0.5 * (protoLayer.max(binZ) - protoLayer.min(binZ));
337  layerThickness = (protoLayer.max(binX) - protoLayer.min(binX));
338  break;
339  }
340  case BinningValue::binY: {
341  layerHalf1 = 0.5 * (protoLayer.max(binX) - protoLayer.min(binX));
342  layerHalf2 = 0.5 * (protoLayer.max(binZ) - protoLayer.min(binZ));
343  layerThickness = (protoLayer.max(binY) - protoLayer.min(binY));
344  break;
345  }
346  case BinningValue::binZ: {
347  layerHalf1 = 0.5 * (protoLayer.max(binX) - protoLayer.min(binX));
348  layerHalf2 = 0.5 * (protoLayer.max(binY) - protoLayer.min(binY));
349  layerThickness = (protoLayer.max(binZ) - protoLayer.min(binZ));
350  break;
351  }
352  default:
353  throw std::invalid_argument("Invalid binning value");
354  }
355 
356  double centerX = 0.5 * (protoLayer.max(binX) + protoLayer.min(binX));
357  double centerY = 0.5 * (protoLayer.max(binY) + protoLayer.min(binY));
358  double centerZ = 0.5 * (protoLayer.max(binZ) + protoLayer.min(binZ));
359 
360  ACTS_VERBOSE("Creating a plane Layer:");
361  ACTS_VERBOSE(" - with layer center = "
362  << "(" << centerX << ", " << centerY << ", " << centerZ << ")");
363  ACTS_VERBOSE(" - from X min/max = " << protoLayer.min(binX) << " / "
364  << protoLayer.max(binX));
365  ACTS_VERBOSE(" - from Y min/max = " << protoLayer.min(binY) << " / "
366  << protoLayer.max(binY));
367  ACTS_VERBOSE(" - with Z thickness = " << layerThickness);
368  ACTS_VERBOSE(" - incl envelope = " << protoLayer.envelope[bValue][0u]
369  << " / "
370  << protoLayer.envelope[bValue][1u]);
371 
372  // create the layer transforms if not given
373  // we need to transform in case centerX/centerY/centerZ != 0, so that the
374  // layer will be correctly defined
375  Translation3 addTranslation(0., 0., 0.);
376  if (transform.isApprox(Transform3::Identity())) {
377  addTranslation = Translation3(centerX, centerY, centerZ);
378  ACTS_VERBOSE(" - layer shift = "
379  << "(" << centerX << ", " << centerY << ", " << centerZ
380  << ")");
381  }
382 
383  std::unique_ptr<SurfaceArray> sArray;
384  if (!surfaces.empty()) {
385  sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnPlane(
386  gctx, std::move(surfaces), bins1, bins2, bValue, protoLayer, transform);
387 
388  checkBinning(gctx, *sArray);
389  }
390 
391  // create the layer and push it back
392  std::shared_ptr<const PlanarBounds> pBounds(
393  new RectangleBounds(layerHalf1, layerHalf2));
394 
395  // create the layer
396  MutableLayerPtr pLayer =
397  PlaneLayer::create(addTranslation * transform, pBounds, std::move(sArray),
398  layerThickness, std::move(ad), active);
399 
400  if (!pLayer) {
401  ACTS_ERROR("Creation of plane layer did not succeed!");
402  }
403  associateSurfacesToLayer(*pLayer);
404 
405  // now return
406  return pLayer;
407 }
408 
409 void Acts::LayerCreator::associateSurfacesToLayer(Layer& layer) const {
410  if (layer.surfaceArray() != nullptr) {
411  auto surfaces = layer.surfaceArray()->surfaces();
412 
413  for (auto& surface : surfaces) {
414  auto mutableSurface = const_cast<Surface*>(surface);
415  mutableSurface->associateLayer(layer);
416  }
417  }
418 }
419 
420 bool Acts::LayerCreator::checkBinning(const GeometryContext& gctx,
421  const SurfaceArray& sArray) const {
422  // do consistency check: can we access all sensitive surfaces
423  // through the binning? If not, surfaces get lost and the binning does not
424  // work
425 
426  ACTS_VERBOSE("Performing consistency check")
427 
428  std::vector<const Surface*> surfaces = sArray.surfaces();
429  std::set<const Surface*> sensitiveSurfaces(surfaces.begin(), surfaces.end());
430  std::set<const Surface*> accessibleSurfaces;
431  size_t nEmptyBins = 0;
432  size_t nBinsChecked = 0;
433 
434  // iterate over all bins
435  size_t size = sArray.size();
436  for (size_t b = 0; b < size; ++b) {
437  std::vector<const Surface*> binContent = sArray.at(b);
438  // we don't check under/overflow bins
439  if (!sArray.isValidBin(b)) {
440  continue;
441  }
442  for (const auto& srf : binContent) {
443  accessibleSurfaces.insert(srf);
444  }
445  if (binContent.empty()) {
446  nEmptyBins++;
447  }
448  nBinsChecked++;
449  }
450 
451  std::vector<const Acts::Surface*> diff;
452  std::set_difference(sensitiveSurfaces.begin(), sensitiveSurfaces.end(),
453  accessibleSurfaces.begin(), accessibleSurfaces.end(),
454  std::inserter(diff, diff.begin()));
455 
456  ACTS_VERBOSE(" - Checked " << nBinsChecked << " valid bins");
457 
458  if (nEmptyBins > 0) {
459  ACTS_ERROR(" -- Not all bins point to surface. " << nEmptyBins << " empty");
460  } else {
461  ACTS_VERBOSE(" -- All bins point to a surface");
462  }
463 
464  if (!diff.empty()) {
465  ACTS_ERROR(
466  " -- Not all sensitive surfaces are accessible through binning. "
467  "sensitive: "
468  << sensitiveSurfaces.size()
469  << " accessible: " << accessibleSurfaces.size());
470 
471  // print all inaccessibles
472  ACTS_ERROR(" -- Inaccessible surfaces: ");
473  for (const auto& srf : diff) {
474  // have to choose BinningValue here
475  Vector3 ctr = srf->binningPosition(gctx, binR);
476  ACTS_ERROR(" Surface(x=" << ctr.x() << ", y=" << ctr.y()
477  << ", z=" << ctr.z() << ", r=" << perp(ctr)
478  << ", phi=" << phi(ctr) << ")");
479  }
480 
481  } else {
482  ACTS_VERBOSE(" -- All sensitive surfaces are accessible through binning.");
483  }
484 
485  return nEmptyBins == 0 && diff.empty();
486 }