Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-12 07:51:57

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008 
0009 #include "Acts/Geometry/LayerCreator.hpp"
0010 
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/CylinderLayer.hpp"
0013 #include "Acts/Geometry/DiscLayer.hpp"
0014 #include "Acts/Geometry/Extent.hpp"
0015 #include "Acts/Geometry/Layer.hpp"
0016 #include "Acts/Geometry/PlaneLayer.hpp"
0017 #include "Acts/Geometry/ProtoLayer.hpp"
0018 #include "Acts/Geometry/SurfaceArrayCreator.hpp"
0019 #include "Acts/Surfaces/CylinderBounds.hpp"
0020 #include "Acts/Surfaces/RadialBounds.hpp"
0021 #include "Acts/Surfaces/RectangleBounds.hpp"
0022 #include "Acts/Surfaces/Surface.hpp"
0023 
0024 #include <algorithm>
0025 #include <iterator>
0026 #include <ostream>
0027 #include <set>
0028 #include <utility>
0029 
0030 namespace Acts {
0031 
0032 using VectorHelpers::perp;
0033 using VectorHelpers::phi;
0034 
0035 LayerCreator::LayerCreator(const LayerCreator::Config& lcConfig,
0036                            std::unique_ptr<const Logger> logger)
0037     : m_cfg(lcConfig), m_logger(std::move(logger)) {}
0038 
0039 void LayerCreator::setConfiguration(const LayerCreator::Config& lcConfig) {
0040   // @todo check consistency
0041   // copy the configuration
0042   m_cfg = lcConfig;
0043 }
0044 
0045 void LayerCreator::setLogger(std::unique_ptr<const Logger> newLogger) {
0046   m_logger = std::move(newLogger);
0047 }
0048 
0049 MutableLayerPtr LayerCreator::cylinderLayer(
0050     const GeometryContext& gctx,
0051     std::vector<std::shared_ptr<const Surface>> surfaces, std::size_t binsPhi,
0052     std::size_t binsZ, std::optional<ProtoLayer> _protoLayer,
0053     const Transform3& transform, std::unique_ptr<ApproachDescriptor> ad) const {
0054   ProtoLayer protoLayer =
0055       _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
0056   if (!_protoLayer) {
0057     protoLayer.envelope[AxisDirection::AxisR] = m_cfg.defaultEnvelopeR;
0058     protoLayer.envelope[AxisDirection::AxisZ] = m_cfg.defaultEnvelopeZ;
0059   }
0060 
0061   // Remaining layer parameters - they include the envelopes
0062   double layerR = protoLayer.medium(AxisDirection::AxisR);
0063   double layerZ = protoLayer.medium(AxisDirection::AxisZ);
0064   double layerHalfZ = 0.5 * protoLayer.range(AxisDirection::AxisZ);
0065   double layerThickness = protoLayer.range(AxisDirection::AxisR);
0066 
0067   ACTS_VERBOSE("Creating a cylindrical Layer:");
0068   ACTS_VERBOSE(" - with layer R     = " << layerR);
0069   ACTS_VERBOSE(" - from R min/max   = "
0070                << protoLayer.min(AxisDirection::AxisR, false) << " / "
0071                << protoLayer.max(AxisDirection::AxisR, false));
0072   ACTS_VERBOSE(" - with R thickness = " << layerThickness);
0073   ACTS_VERBOSE("   - incl envelope  = "
0074                << protoLayer.envelope[AxisDirection::AxisR][0u] << " / "
0075                << protoLayer.envelope[AxisDirection::AxisR][1u]);
0076 
0077   ACTS_VERBOSE(" - with z min/max   = "
0078                << protoLayer.min(AxisDirection::AxisZ, false) << " (-"
0079                << protoLayer.envelope[AxisDirection::AxisZ][0u] << ") / "
0080                << protoLayer.max(AxisDirection::AxisZ, false) << " (+"
0081                << protoLayer.envelope[AxisDirection::AxisZ][1u] << ")");
0082 
0083   ACTS_VERBOSE(" - z center         = " << layerZ);
0084   ACTS_VERBOSE(" - halflength z     = " << layerHalfZ);
0085 
0086   // create the layer transforms if not given
0087   // we need to transform in case layerZ != 0, so that the layer will be
0088   // correctly defined using the halflength
0089   Translation3 addTranslation(0., 0., 0.);
0090   if (transform.isApprox(Transform3::Identity())) {
0091     addTranslation = Translation3(0., 0., layerZ);
0092     ACTS_VERBOSE(" - layer z shift  = " << -layerZ);
0093   }
0094 
0095   ACTS_VERBOSE(" - with phi min/max = "
0096                << protoLayer.min(AxisDirection::AxisPhi, false) << " / "
0097                << protoLayer.max(AxisDirection::AxisPhi, false));
0098   ACTS_VERBOSE(" - # of modules     = " << surfaces.size() << " ordered in ( "
0099                                         << binsPhi << " x " << binsZ << ")");
0100   std::unique_ptr<SurfaceArray> sArray;
0101   if (!surfaces.empty()) {
0102     sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnCylinder(
0103         gctx, std::move(surfaces), binsPhi, binsZ, protoLayer);
0104 
0105     checkBinning(gctx, *sArray);
0106   }
0107 
0108   // create the layer and push it back
0109   std::shared_ptr<const CylinderBounds> cBounds(
0110       new CylinderBounds(layerR, layerHalfZ));
0111 
0112   // create the layer
0113   MutableLayerPtr cLayer = CylinderLayer::create(
0114       addTranslation * transform, cBounds, std::move(sArray), layerThickness,
0115       std::move(ad), active);
0116 
0117   if (!cLayer) {
0118     ACTS_ERROR("Creation of cylinder layer did not succeed!");
0119   }
0120   associateSurfacesToLayer(*cLayer);
0121 
0122   // now return
0123   return cLayer;
0124 }
0125 
0126 MutableLayerPtr LayerCreator::cylinderLayer(
0127     const GeometryContext& gctx,
0128     std::vector<std::shared_ptr<const Surface>> surfaces, BinningType bTypePhi,
0129     BinningType bTypeZ, std::optional<ProtoLayer> _protoLayer,
0130     const Transform3& transform, std::unique_ptr<ApproachDescriptor> ad) const {
0131   ProtoLayer protoLayer =
0132       _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
0133   if (!_protoLayer) {
0134     protoLayer.envelope[AxisDirection::AxisR] = m_cfg.defaultEnvelopeR;
0135     protoLayer.envelope[AxisDirection::AxisZ] = m_cfg.defaultEnvelopeZ;
0136   }
0137 
0138   // remaining layer parameters
0139   double layerR = protoLayer.medium(AxisDirection::AxisR);
0140   double layerZ = protoLayer.medium(AxisDirection::AxisZ);
0141   double layerHalfZ = 0.5 * protoLayer.range(AxisDirection::AxisZ);
0142   double layerThickness = protoLayer.range(AxisDirection::AxisR);
0143 
0144   // adjust the layer radius
0145   ACTS_VERBOSE("Creating a cylindrical Layer:");
0146   ACTS_VERBOSE(" - with layer R     = " << layerR);
0147   ACTS_VERBOSE(" - from R min/max   = "
0148                << protoLayer.min(AxisDirection::AxisR, false) << " / "
0149                << protoLayer.max(AxisDirection::AxisR, false));
0150   ACTS_VERBOSE(" - with R thickness = " << layerThickness);
0151   ACTS_VERBOSE("   - incl envelope  = "
0152                << protoLayer.envelope[AxisDirection::AxisR][0u] << " / "
0153                << protoLayer.envelope[AxisDirection::AxisR][1u]);
0154   ACTS_VERBOSE(" - with z min/max   = "
0155                << protoLayer.min(AxisDirection::AxisZ, false) << " (-"
0156                << protoLayer.envelope[AxisDirection::AxisZ][0u] << ") / "
0157                << protoLayer.max(AxisDirection::AxisZ, false) << " (+"
0158                << protoLayer.envelope[AxisDirection::AxisZ][1u] << ")");
0159   ACTS_VERBOSE(" - z center         = " << layerZ);
0160   ACTS_VERBOSE(" - halflength z     = " << layerHalfZ);
0161 
0162   // create the layer transforms if not given
0163   // we need to transform in case layerZ != 0, so that the layer will be
0164   // correctly defined using the halflength
0165   // create the layer transforms if not given
0166   Translation3 addTranslation(0., 0., 0.);
0167   if (transform.isApprox(Transform3::Identity()) && bTypeZ == equidistant) {
0168     addTranslation = Translation3(0., 0., layerZ);
0169     ACTS_VERBOSE(" - layer z shift    = " << -layerZ);
0170   }
0171 
0172   ACTS_VERBOSE(" - with phi min/max = "
0173                << protoLayer.min(AxisDirection::AxisPhi, false) << " / "
0174                << protoLayer.max(AxisDirection::AxisPhi, false));
0175   ACTS_VERBOSE(" - # of modules     = " << surfaces.size() << "");
0176 
0177   // create the surface array
0178   std::unique_ptr<SurfaceArray> sArray;
0179   if (!surfaces.empty()) {
0180     sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnCylinder(
0181         gctx, std::move(surfaces), bTypePhi, bTypeZ, protoLayer);
0182 
0183     checkBinning(gctx, *sArray);
0184   }
0185 
0186   // create the layer and push it back
0187   std::shared_ptr<const CylinderBounds> cBounds(
0188       new CylinderBounds(layerR, layerHalfZ));
0189 
0190   // create the layer
0191   MutableLayerPtr cLayer = CylinderLayer::create(
0192       addTranslation * transform, cBounds, std::move(sArray), layerThickness,
0193       std::move(ad), active);
0194 
0195   if (!cLayer) {
0196     ACTS_ERROR("Creation of cylinder layer did not succeed!");
0197   }
0198   associateSurfacesToLayer(*cLayer);
0199 
0200   // now return
0201   return cLayer;
0202 }
0203 
0204 MutableLayerPtr LayerCreator::discLayer(
0205     const GeometryContext& gctx,
0206     std::vector<std::shared_ptr<const Surface>> surfaces, std::size_t binsR,
0207     std::size_t binsPhi, std::optional<ProtoLayer> _protoLayer,
0208     const Transform3& transform, std::unique_ptr<ApproachDescriptor> ad) const {
0209   ProtoLayer protoLayer =
0210       _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
0211   if (!_protoLayer) {
0212     protoLayer.envelope[AxisDirection::AxisR] = m_cfg.defaultEnvelopeR;
0213     protoLayer.envelope[AxisDirection::AxisZ] = m_cfg.defaultEnvelopeZ;
0214   }
0215 
0216   double layerZ = protoLayer.medium(AxisDirection::AxisZ);
0217   double layerThickness = protoLayer.range(AxisDirection::AxisZ);
0218 
0219   // adjust the layer radius
0220   ACTS_VERBOSE("Creating a disk Layer:");
0221   ACTS_VERBOSE(" - at Z position    = " << layerZ);
0222   ACTS_VERBOSE(" - from Z min/max   = "
0223                << protoLayer.min(AxisDirection::AxisZ, false) << " / "
0224                << protoLayer.max(AxisDirection::AxisZ, false));
0225   ACTS_VERBOSE(" - with Z thickness = " << layerThickness);
0226   ACTS_VERBOSE("   - incl envelope  = "
0227                << protoLayer.envelope[AxisDirection::AxisZ][0u] << " / "
0228                << protoLayer.envelope[AxisDirection::AxisZ][1u]);
0229   ACTS_VERBOSE(" - with R min/max   = "
0230                << protoLayer.min(AxisDirection::AxisR, false) << " (-"
0231                << protoLayer.envelope[AxisDirection::AxisR][0u] << ") / "
0232                << protoLayer.max(AxisDirection::AxisR, false) << " (+"
0233                << protoLayer.envelope[AxisDirection::AxisR][1u] << ")");
0234   ACTS_VERBOSE(" - with phi min/max = "
0235                << protoLayer.min(AxisDirection::AxisPhi, false) << " / "
0236                << protoLayer.max(AxisDirection::AxisPhi, false));
0237   ACTS_VERBOSE(" - # of modules    = " << surfaces.size() << " ordered in ( "
0238                                        << binsR << " x " << binsPhi << ")");
0239 
0240   // create the layer transforms if not given
0241   Translation3 addTranslation(0., 0., 0.);
0242   if (transform.isApprox(Transform3::Identity())) {
0243     addTranslation = Translation3(0., 0., layerZ);
0244   }
0245   // create the surface array
0246   std::unique_ptr<SurfaceArray> sArray;
0247   if (!surfaces.empty()) {
0248     sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnDisc(
0249         gctx, std::move(surfaces), binsR, binsPhi, protoLayer, transform);
0250 
0251     checkBinning(gctx, *sArray);
0252   }
0253 
0254   // create the share disc bounds
0255   auto dBounds = std::make_shared<const RadialBounds>(
0256       protoLayer.min(AxisDirection::AxisR),
0257       protoLayer.max(AxisDirection::AxisR));
0258 
0259   // create the layers
0260   // we use the same transform here as for the layer itself
0261   // for disk this is fine since we don't bin in Z, so does not matter
0262   MutableLayerPtr dLayer =
0263       DiscLayer::create(addTranslation * transform, dBounds, std::move(sArray),
0264                         layerThickness, std::move(ad), active);
0265 
0266   if (!dLayer) {
0267     ACTS_ERROR("Creation of disc layer did not succeed!");
0268   }
0269   associateSurfacesToLayer(*dLayer);
0270   // return the layer
0271   return dLayer;
0272 }
0273 
0274 MutableLayerPtr LayerCreator::discLayer(
0275     const GeometryContext& gctx,
0276     std::vector<std::shared_ptr<const Surface>> surfaces, BinningType bTypeR,
0277     BinningType bTypePhi, std::optional<ProtoLayer> _protoLayer,
0278     const Transform3& transform, std::unique_ptr<ApproachDescriptor> ad) const {
0279   ProtoLayer protoLayer =
0280       _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
0281   if (!_protoLayer) {
0282     protoLayer.envelope[AxisDirection::AxisR] = m_cfg.defaultEnvelopeR;
0283     protoLayer.envelope[AxisDirection::AxisZ] = m_cfg.defaultEnvelopeZ;
0284   }
0285 
0286   double layerZ = protoLayer.medium(AxisDirection::AxisZ);
0287   double layerThickness = protoLayer.range(AxisDirection::AxisZ);
0288 
0289   // adjust the layer radius
0290   ACTS_VERBOSE("Creating a disk Layer:");
0291   ACTS_VERBOSE(" - at Z position    = " << layerZ);
0292   ACTS_VERBOSE(" - from Z min/max   = "
0293                << protoLayer.min(AxisDirection::AxisZ, false) << " / "
0294                << protoLayer.max(AxisDirection::AxisZ, false));
0295   ACTS_VERBOSE(" - with Z thickness = " << layerThickness);
0296   ACTS_VERBOSE("   - incl envelope  = "
0297                << protoLayer.envelope[AxisDirection::AxisZ][0u] << " / "
0298                << protoLayer.envelope[AxisDirection::AxisZ][1u]);
0299   ACTS_VERBOSE(" - with R min/max   = "
0300                << protoLayer.min(AxisDirection::AxisR, false) << " (-"
0301                << protoLayer.envelope[AxisDirection::AxisR][0u] << ") / "
0302                << protoLayer.max(AxisDirection::AxisR, false) << " (+"
0303                << protoLayer.envelope[AxisDirection::AxisR][1u] << ")");
0304   ACTS_VERBOSE(" - with phi min/max = "
0305                << protoLayer.min(AxisDirection::AxisPhi, false) << " / "
0306                << protoLayer.max(AxisDirection::AxisPhi, false));
0307   ACTS_VERBOSE(" - # of modules     = " << surfaces.size());
0308 
0309   // create the layer transforms if not given
0310   Translation3 addTranslation(0., 0., 0.);
0311   if (transform.isApprox(Transform3::Identity())) {
0312     addTranslation = Translation3(0., 0., layerZ);
0313   }
0314 
0315   // create the surface array
0316   std::unique_ptr<SurfaceArray> sArray;
0317   if (!surfaces.empty()) {
0318     sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnDisc(
0319         gctx, std::move(surfaces), bTypeR, bTypePhi, protoLayer, transform);
0320 
0321     checkBinning(gctx, *sArray);
0322   }
0323 
0324   // create the shared disc bounds
0325   auto dBounds = std::make_shared<const RadialBounds>(
0326       protoLayer.min(AxisDirection::AxisR),
0327       protoLayer.max(AxisDirection::AxisR));
0328 
0329   // create the layers
0330   MutableLayerPtr dLayer =
0331       DiscLayer::create(addTranslation * transform, dBounds, std::move(sArray),
0332                         layerThickness, std::move(ad), active);
0333   if (!dLayer) {
0334     ACTS_ERROR("Creation of disc layer did not succeed!");
0335   }
0336   associateSurfacesToLayer(*dLayer);
0337   // return the layer
0338   return dLayer;
0339 }
0340 
0341 MutableLayerPtr LayerCreator::planeLayer(
0342     const GeometryContext& gctx,
0343     std::vector<std::shared_ptr<const Surface>> surfaces, std::size_t bins1,
0344     std::size_t bins2, AxisDirection aDir,
0345     std::optional<ProtoLayer> _protoLayer, const Transform3& transform,
0346     std::unique_ptr<ApproachDescriptor> ad) const {
0347   ProtoLayer protoLayer =
0348       _protoLayer ? *_protoLayer : ProtoLayer(gctx, surfaces);
0349   if (!_protoLayer) {
0350     protoLayer.envelope[AxisDirection::AxisR] = m_cfg.defaultEnvelopeR;
0351     protoLayer.envelope[AxisDirection::AxisZ] = m_cfg.defaultEnvelopeZ;
0352   }
0353 
0354   // remaining layer parameters
0355   double layerHalf1 = 0, layerHalf2 = 0, layerThickness = 0;
0356   switch (aDir) {
0357     case AxisDirection::AxisX: {
0358       layerHalf1 = 0.5 * (protoLayer.max(AxisDirection::AxisY) -
0359                           protoLayer.min(AxisDirection::AxisY));
0360       layerHalf2 = 0.5 * (protoLayer.max(AxisDirection::AxisZ) -
0361                           protoLayer.min(AxisDirection::AxisZ));
0362       layerThickness = (protoLayer.max(AxisDirection::AxisX) -
0363                         protoLayer.min(AxisDirection::AxisX));
0364       break;
0365     }
0366     case AxisDirection::AxisY: {
0367       layerHalf1 = 0.5 * (protoLayer.max(AxisDirection::AxisX) -
0368                           protoLayer.min(AxisDirection::AxisX));
0369       layerHalf2 = 0.5 * (protoLayer.max(AxisDirection::AxisZ) -
0370                           protoLayer.min(AxisDirection::AxisZ));
0371       layerThickness = (protoLayer.max(AxisDirection::AxisY) -
0372                         protoLayer.min(AxisDirection::AxisY));
0373       break;
0374     }
0375     case AxisDirection::AxisZ: {
0376       layerHalf1 = 0.5 * (protoLayer.max(AxisDirection::AxisX) -
0377                           protoLayer.min(AxisDirection::AxisX));
0378       layerHalf2 = 0.5 * (protoLayer.max(AxisDirection::AxisY) -
0379                           protoLayer.min(AxisDirection::AxisY));
0380       layerThickness = (protoLayer.max(AxisDirection::AxisZ) -
0381                         protoLayer.min(AxisDirection::AxisZ));
0382       break;
0383     }
0384     default:
0385       throw std::invalid_argument("Invalid binning value");
0386   }
0387 
0388   double centerX = 0.5 * (protoLayer.max(AxisDirection::AxisX) +
0389                           protoLayer.min(AxisDirection::AxisX));
0390   double centerY = 0.5 * (protoLayer.max(AxisDirection::AxisY) +
0391                           protoLayer.min(AxisDirection::AxisY));
0392   double centerZ = 0.5 * (protoLayer.max(AxisDirection::AxisZ) +
0393                           protoLayer.min(AxisDirection::AxisZ));
0394 
0395   ACTS_VERBOSE("Creating a plane Layer:");
0396   ACTS_VERBOSE(" - with layer center     = "
0397                << "(" << centerX << ", " << centerY << ", " << centerZ << ")");
0398   ACTS_VERBOSE(" - from X min/max   = "
0399                << protoLayer.min(AxisDirection::AxisX) << " / "
0400                << protoLayer.max(AxisDirection::AxisX));
0401   ACTS_VERBOSE(" - from Y min/max   = "
0402                << protoLayer.min(AxisDirection::AxisY) << " / "
0403                << protoLayer.max(AxisDirection::AxisY));
0404   ACTS_VERBOSE(" - with Z thickness = " << layerThickness);
0405   ACTS_VERBOSE("   - incl envelope  = " << protoLayer.envelope[aDir][0u]
0406                                         << " / "
0407                                         << protoLayer.envelope[aDir][1u]);
0408 
0409   // create the layer transforms if not given
0410   // we need to transform in case centerX/centerY/centerZ != 0, so that the
0411   // layer will be correctly defined
0412   Translation3 addTranslation(0., 0., 0.);
0413   if (transform.isApprox(Transform3::Identity())) {
0414     addTranslation = Translation3(centerX, centerY, centerZ);
0415     ACTS_VERBOSE(" - layer shift  = " << "(" << centerX << ", " << centerY
0416                                       << ", " << centerZ << ")");
0417   }
0418 
0419   std::unique_ptr<SurfaceArray> sArray;
0420   if (!surfaces.empty()) {
0421     sArray = m_cfg.surfaceArrayCreator->surfaceArrayOnPlane(
0422         gctx, std::move(surfaces), bins1, bins2, aDir, protoLayer, transform);
0423 
0424     checkBinning(gctx, *sArray);
0425   }
0426 
0427   // create the layer and push it back
0428   auto pBounds = std::make_shared<RectangleBounds>(layerHalf1, layerHalf2);
0429 
0430   // create the layer
0431   MutableLayerPtr pLayer =
0432       PlaneLayer::create(addTranslation * transform, pBounds, std::move(sArray),
0433                          layerThickness, std::move(ad), active);
0434 
0435   if (!pLayer) {
0436     ACTS_ERROR("Creation of plane layer did not succeed!");
0437   }
0438   associateSurfacesToLayer(*pLayer);
0439 
0440   // now return
0441   return pLayer;
0442 }
0443 
0444 void LayerCreator::associateSurfacesToLayer(Layer& layer) const {
0445   if (layer.surfaceArray() != nullptr) {
0446     auto surfaces = layer.surfaceArray()->surfaces();
0447 
0448     for (auto& surface : surfaces) {
0449       auto mutableSurface = const_cast<Surface*>(surface);
0450       mutableSurface->associateLayer(layer);
0451     }
0452   }
0453 }
0454 
0455 bool LayerCreator::checkBinning(const GeometryContext& gctx,
0456                                 const SurfaceArray& sArray) const {
0457   // do consistency check: can we access all sensitive surfaces
0458   // through the binning? If not, surfaces get lost and the binning does not
0459   // work
0460 
0461   ACTS_VERBOSE("Performing consistency check");
0462 
0463   std::vector<const Surface*> surfaces = sArray.surfaces();
0464   std::set<const Surface*> sensitiveSurfaces(surfaces.begin(), surfaces.end());
0465   std::set<const Surface*> accessibleSurfaces;
0466   std::size_t nEmptyBins = 0;
0467   std::size_t nBinsChecked = 0;
0468 
0469   // iterate over all bins
0470   std::size_t size = sArray.size();
0471   for (std::size_t b = 0; b < size; ++b) {
0472     std::vector<const Surface*> binContent = sArray.at(b);
0473     // we don't check under/overflow bins
0474     if (!sArray.isValidBin(b)) {
0475       continue;
0476     }
0477     for (const auto& srf : binContent) {
0478       accessibleSurfaces.insert(srf);
0479     }
0480     if (binContent.empty()) {
0481       nEmptyBins++;
0482     }
0483     nBinsChecked++;
0484   }
0485 
0486   std::vector<const Surface*> diff;
0487   std::set_difference(sensitiveSurfaces.begin(), sensitiveSurfaces.end(),
0488                       accessibleSurfaces.begin(), accessibleSurfaces.end(),
0489                       std::inserter(diff, diff.begin()));
0490 
0491   ACTS_VERBOSE(" - Checked " << nBinsChecked << " valid bins");
0492 
0493   if (nEmptyBins > 0) {
0494     ACTS_ERROR(" -- Not all bins point to surface. " << nEmptyBins << " empty");
0495   } else {
0496     ACTS_VERBOSE(" -- All bins point to a surface");
0497   }
0498 
0499   if (!diff.empty()) {
0500     ACTS_ERROR(
0501         " -- Not all sensitive surfaces are accessible through binning. "
0502         "sensitive: "
0503         << sensitiveSurfaces.size()
0504         << "    accessible: " << accessibleSurfaces.size());
0505 
0506     // print all inaccessibles
0507     ACTS_ERROR(" -- Inaccessible surfaces: ");
0508     for (const auto& srf : diff) {
0509       // have to choose AxisDirection here
0510       Vector3 ctr = srf->referencePosition(gctx, AxisDirection::AxisR);
0511       ACTS_ERROR(" Surface(x=" << ctr.x() << ", y=" << ctr.y()
0512                                << ", z=" << ctr.z() << ", r=" << perp(ctr)
0513                                << ", phi=" << phi(ctr) << ")");
0514     }
0515 
0516   } else {
0517     ACTS_VERBOSE(" -- All sensitive surfaces are accessible through binning.");
0518   }
0519 
0520   return nEmptyBins == 0 && diff.empty();
0521 }
0522 
0523 }  // namespace Acts