src/Geometry/SurfaceArrayCreator.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2016-2020 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Geometry/SurfaceArrayCreator.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Surfaces/PlanarBounds.hpp"
0013 #include "Acts/Surfaces/Surface.hpp"
0014 #include "Acts/Surfaces/SurfaceArray.hpp"
0015 #include "Acts/Surfaces/SurfaceBounds.hpp"
0016 #include "Acts/Utilities/BinningType.hpp"
0017 #include "Acts/Utilities/Grid.hpp"
0018 #include "Acts/Utilities/Helpers.hpp"
0019 #include "Acts/Utilities/IAxis.hpp"
0020 #include "Acts/Utilities/detail/grid_helper.hpp"
0021
0022 #include <algorithm>
0023 #include <cmath>
0024 #include <stdexcept>
0025
0026 using Acts::VectorHelpers::perp;
0027 using Acts::VectorHelpers::phi;
0028
0029 std::unique_ptr<Acts::SurfaceArray>
0030 Acts::SurfaceArrayCreator::surfaceArrayOnCylinder(
0031     const GeometryContext& gctx,
0032     std::vector<std::shared_ptr<const Surface>> surfaces, std::size_t binsPhi,
0033     std::size_t binsZ, std::optional<ProtoLayer> protoLayerOpt,
0034     const Transform3& transform) const {
0035   std::vector<const Surface*> surfacesRaw = unpack_shared_vector(surfaces);
0036   // Check if we have proto layer, else build it
0037   ProtoLayer protoLayer =
0038       protoLayerOpt ? *protoLayerOpt : ProtoLayer(gctx, surfacesRaw);
0039
0040   ACTS_VERBOSE("Creating a SurfaceArray on a cylinder");
0041   ACTS_VERBOSE(" -- with " << surfaces.size() << " surfaces.")
0042   ACTS_VERBOSE(" -- with phi x z  = " << binsPhi << " x " << binsZ << " = "
0043                                       << binsPhi * binsZ << " bins.");
0044
0045   Transform3 ftransform = transform;
0046   ProtoAxis pAxisPhi = createEquidistantAxis(gctx, surfacesRaw, binPhi,
0047                                              protoLayer, ftransform, binsPhi);
0048   ProtoAxis pAxisZ = createEquidistantAxis(gctx, surfacesRaw, binZ, protoLayer,
0049                                            ftransform, binsZ);
0050
0051   double R = protoLayer.medium(binR, true);
0052
0053   Transform3 itransform = ftransform.inverse();
0054   // Transform lambda captures the transform matrix
0055   auto globalToLocal = [ftransform](const Vector3& pos) {
0056     Vector3 loc = ftransform * pos;
0057     return Vector2(phi(loc), loc.z());
0058   };
0059   auto localToGlobal = [itransform, R](const Vector2& loc) {
0060     return itransform *
0061            Vector3(R * std::cos(loc[0]), R * std::sin(loc[0]), loc[1]);
0062   };
0063
0064   std::unique_ptr<SurfaceArray::ISurfaceGridLookup> sl =
0065       makeSurfaceGridLookup2D<detail::AxisBoundaryType::Closed,
0066                               detail::AxisBoundaryType::Bound>(
0067           globalToLocal, localToGlobal, pAxisPhi, pAxisZ);
0068
0069   sl->fill(gctx, surfacesRaw);
0070   completeBinning(gctx, *sl, surfacesRaw);
0071
0072   return std::make_unique<SurfaceArray>(std::move(sl), std::move(surfaces),
0073                                         ftransform);
0074 }
0075
0076 std::unique_ptr<Acts::SurfaceArray>
0077 Acts::SurfaceArrayCreator::surfaceArrayOnCylinder(
0078     const GeometryContext& gctx,
0079     std::vector<std::shared_ptr<const Surface>> surfaces, BinningType bTypePhi,
0080     BinningType bTypeZ, std::optional<ProtoLayer> protoLayerOpt,
0081     const Transform3& transform) const {
0082   std::vector<const Surface*> surfacesRaw = unpack_shared_vector(surfaces);
0083   // check if we have proto layer, else build it
0084   ProtoLayer protoLayer =
0085       protoLayerOpt ? *protoLayerOpt : ProtoLayer(gctx, surfacesRaw);
0086
0087   double R = protoLayer.medium(binR, true);
0088
0089   ProtoAxis pAxisPhi;
0090   ProtoAxis pAxisZ;
0091
0092   Transform3 ftransform = transform;
0093
0094   if (bTypePhi == equidistant) {
0095     pAxisPhi = createEquidistantAxis(gctx, surfacesRaw, binPhi, protoLayer,
0096                                      ftransform, 0);
0097   } else {
0098     pAxisPhi =
0099         createVariableAxis(gctx, surfacesRaw, binPhi, protoLayer, ftransform);
0100   }
0101
0102   if (bTypeZ == equidistant) {
0103     pAxisZ =
0104         createEquidistantAxis(gctx, surfacesRaw, binZ, protoLayer, ftransform);
0105   } else {
0106     pAxisZ =
0107         createVariableAxis(gctx, surfacesRaw, binZ, protoLayer, ftransform);
0108   }
0109
0110   Transform3 itransform = ftransform.inverse();
0111   auto globalToLocal = [ftransform](const Vector3& pos) {
0112     Vector3 loc = ftransform * pos;
0113     return Vector2(phi(loc), loc.z());
0114   };
0115   auto localToGlobal = [itransform, R](const Vector2& loc) {
0116     return itransform *
0117            Vector3(R * std::cos(loc[0]), R * std::sin(loc[0]), loc[1]);
0118   };
0119
0120   std::unique_ptr<SurfaceArray::ISurfaceGridLookup> sl =
0121       makeSurfaceGridLookup2D<detail::AxisBoundaryType::Closed,
0122                               detail::AxisBoundaryType::Bound>(
0123           globalToLocal, localToGlobal, pAxisPhi, pAxisZ);
0124
0125   sl->fill(gctx, surfacesRaw);
0126   completeBinning(gctx, *sl, surfacesRaw);
0127
0128   // get the number of bins
0129   auto axes = sl->getAxes();
0130   std::size_t bins0 = axes.at(0)->getNBins();
0131   std::size_t bins1 = axes.at(1)->getNBins();
0132
0133   ACTS_VERBOSE("Creating a SurfaceArray on a cylinder");
0134   ACTS_VERBOSE(" -- with " << surfaces.size() << " surfaces.")
0135   ACTS_VERBOSE(" -- with phi x z  = " << bins0 << " x " << bins1 << " = "
0136                                       << bins0 * bins1 << " bins.");
0137
0138   return std::make_unique<SurfaceArray>(std::move(sl), std::move(surfaces),
0139                                         ftransform);
0140 }
0141
0142 std::unique_ptr<Acts::SurfaceArray>
0143 Acts::SurfaceArrayCreator::surfaceArrayOnDisc(
0144     const GeometryContext& gctx,
0145     std::vector<std::shared_ptr<const Surface>> surfaces, std::size_t binsR,
0146     std::size_t binsPhi, std::optional<ProtoLayer> protoLayerOpt,
0147     const Transform3& transform) const {
0148   std::vector<const Surface*> surfacesRaw = unpack_shared_vector(surfaces);
0149   // check if we have proto layer, else build it
0150   ProtoLayer protoLayer =
0151       protoLayerOpt ? *protoLayerOpt : ProtoLayer(gctx, surfacesRaw);
0152
0153   ACTS_VERBOSE("Creating a SurfaceArray on a disc");
0154
0155   Transform3 ftransform = transform;
0156   ProtoAxis pAxisR = createEquidistantAxis(gctx, surfacesRaw, binR, protoLayer,
0157                                            ftransform, binsR);
0158   ProtoAxis pAxisPhi = createEquidistantAxis(gctx, surfacesRaw, binPhi,
0159                                              protoLayer, ftransform, binsPhi);
0160
0161   double Z = protoLayer.medium(binZ, true);
0162   ACTS_VERBOSE("- z-position of disk estimated as " << Z);
0163
0164   Transform3 itransform = transform.inverse();
0165   // transform lambda captures the transform matrix
0166   auto globalToLocal = [ftransform](const Vector3& pos) {
0167     Vector3 loc = ftransform * pos;
0168     return Vector2(perp(loc), phi(loc));
0169   };
0170   auto localToGlobal = [itransform, Z](const Vector2& loc) {
0171     return itransform *
0172            Vector3(loc[0] * std::cos(loc[1]), loc[0] * std::sin(loc[1]), Z);
0173   };
0174
0175   std::unique_ptr<SurfaceArray::ISurfaceGridLookup> sl =
0176       makeSurfaceGridLookup2D<detail::AxisBoundaryType::Bound,
0177                               detail::AxisBoundaryType::Closed>(
0178           globalToLocal, localToGlobal, pAxisR, pAxisPhi);
0179
0180   // get the number of bins
0181   auto axes = sl->getAxes();
0182   std::size_t bins0 = axes.at(0)->getNBins();
0183   std::size_t bins1 = axes.at(1)->getNBins();
0184
0185   ACTS_VERBOSE(" -- with " << surfaces.size() << " surfaces.")
0186   ACTS_VERBOSE(" -- with r x phi  = " << bins0 << " x " << bins1 << " = "
0187                                       << bins0 * bins1 << " bins.");
0188   sl->fill(gctx, surfacesRaw);
0189   completeBinning(gctx, *sl, surfacesRaw);
0190
0191   return std::make_unique<SurfaceArray>(std::move(sl), std::move(surfaces),
0192                                         ftransform);
0193 }
0194
0195 std::unique_ptr<Acts::SurfaceArray>
0196 Acts::SurfaceArrayCreator::surfaceArrayOnDisc(
0197     const GeometryContext& gctx,
0198     std::vector<std::shared_ptr<const Surface>> surfaces, BinningType bTypeR,
0199     BinningType bTypePhi, std::optional<ProtoLayer> protoLayerOpt,
0200     const Transform3& transform) const {
0201   std::vector<const Surface*> surfacesRaw = unpack_shared_vector(surfaces);
0202   // check if we have proto layer, else build it
0203   ProtoLayer protoLayer =
0204       protoLayerOpt ? *protoLayerOpt : ProtoLayer(gctx, surfacesRaw);
0205
0206   ACTS_VERBOSE("Creating a SurfaceArray on a disc");
0207
0208   ProtoAxis pAxisPhi;
0209   ProtoAxis pAxisR;
0210
0211   Transform3 ftransform = transform;
0212
0213   if (bTypeR == equidistant) {
0214     pAxisR =
0215         createEquidistantAxis(gctx, surfacesRaw, binR, protoLayer, ftransform);
0216   } else {
0217     pAxisR =
0218         createVariableAxis(gctx, surfacesRaw, binR, protoLayer, ftransform);
0219   }
0220
0221   // if we have more than one R ring, we need to figure out
0222   // the number of phi bins.
0223   if (pAxisR.nBins > 1) {
0224     // more than one R-Ring, we need to adjust
0225     // this FORCES equidistant binning
0226     std::vector<std::vector<const Surface*>> phiModules(pAxisR.nBins);
0227     for (const auto& srf : surfacesRaw) {
0228       Vector3 bpos = srf->binningPosition(gctx, binR);
0229       std::size_t bin = pAxisR.getBin(perp(bpos));
0230       phiModules.at(bin).push_back(srf);
0231     }
0232
0233     std::vector<std::size_t> nPhiModules;
0234     auto matcher = m_cfg.surfaceMatcher;
0235     auto equal = [&gctx, &matcher](const Surface* a, const Surface* b) {
0236       return matcher(gctx, binPhi, a, b);
0237     };
0238
0239     std::transform(
0240         phiModules.begin(), phiModules.end(), std::back_inserter(nPhiModules),
0241         [&equal,
0242          this](const std::vector<const Surface*>& surfaces_) -> std::size_t {
0243           return this->findKeySurfaces(surfaces_, equal).size();
0244         });
0245
0246     // @FIXME: Problem: phi binning runs rotation to optimize
0247     // for bin edges. This FAILS after this modification, since
0248     // the bin count is the one from the lowest module-count bin,
0249     // but the rotation is done considering all bins.
0250     // This might be resolved through bin completion, but not sure.
0251     // @TODO: check in extrapolation
0252     std::size_t nBinsPhi =
0253         (*std::min_element(nPhiModules.begin(), nPhiModules.end()));
0254     pAxisPhi = createEquidistantAxis(gctx, surfacesRaw, binPhi, protoLayer,
0255                                      ftransform, nBinsPhi);
0256
0257   } else {
0258     // use regular determination
0259     if (bTypePhi == equidistant) {
0260       pAxisPhi = createEquidistantAxis(gctx, surfacesRaw, binPhi, protoLayer,
0261                                        ftransform, 0);
0262     } else {
0263       pAxisPhi =
0264           createVariableAxis(gctx, surfacesRaw, binPhi, protoLayer, ftransform);
0265     }
0266   }
0267
0268   double Z = protoLayer.medium(binZ, true);
0269   ACTS_VERBOSE("- z-position of disk estimated as " << Z);
0270
0271   Transform3 itransform = ftransform.inverse();
0272   // transform lambda captures the transform matrix
0273   auto globalToLocal = [ftransform](const Vector3& pos) {
0274     Vector3 loc = ftransform * pos;
0275     return Vector2(perp(loc), phi(loc));
0276   };
0277   auto localToGlobal = [itransform, Z](const Vector2& loc) {
0278     return itransform *
0279            Vector3(loc[0] * std::cos(loc[1]), loc[0] * std::sin(loc[1]), Z);
0280   };
0281
0282   std::unique_ptr<SurfaceArray::ISurfaceGridLookup> sl =
0283       makeSurfaceGridLookup2D<detail::AxisBoundaryType::Bound,
0284                               detail::AxisBoundaryType::Closed>(
0285           globalToLocal, localToGlobal, pAxisR, pAxisPhi);
0286
0287   // get the number of bins
0288   auto axes = sl->getAxes();
0289   std::size_t bins0 = axes.at(0)->getNBins();
0290   std::size_t bins1 = axes.at(1)->getNBins();
0291
0292   ACTS_VERBOSE(" -- with " << surfaces.size() << " surfaces.")
0293   ACTS_VERBOSE(" -- with r x phi  = " << bins0 << " x " << bins1 << " = "
0294                                       << bins0 * bins1 << " bins.");
0295
0296   sl->fill(gctx, surfacesRaw);
0297   completeBinning(gctx, *sl, surfacesRaw);
0298
0299   return std::make_unique<SurfaceArray>(std::move(sl), std::move(surfaces),
0300                                         ftransform);
0301 }
0302
0303 /// SurfaceArrayCreator interface method - create an array on a plane
0304 std::unique_ptr<Acts::SurfaceArray>
0305 Acts::SurfaceArrayCreator::surfaceArrayOnPlane(
0306     const GeometryContext& gctx,
0307     std::vector<std::shared_ptr<const Surface>> surfaces, std::size_t bins1,
0308     std::size_t bins2, BinningValue bValue,
0309     std::optional<ProtoLayer> protoLayerOpt,
0310     const Transform3& transform) const {
0311   std::vector<const Surface*> surfacesRaw = unpack_shared_vector(surfaces);
0312   // check if we have proto layer, else build it
0313   ProtoLayer protoLayer =
0314       protoLayerOpt ? *protoLayerOpt : ProtoLayer(gctx, surfacesRaw);
0315
0316   ACTS_VERBOSE("Creating a SurfaceArray on a plance");
0317   ACTS_VERBOSE(" -- with " << surfaces.size() << " surfaces.")
0318   ACTS_VERBOSE(" -- with " << bins1 << " x " << bins2 << " = " << bins1 * bins2
0319                            << " bins.");
0320   Transform3 ftransform = transform;
0321   Transform3 itransform = transform.inverse();
0322   // transform lambda captures the transform matrix
0323   auto globalToLocal = [ftransform](const Vector3& pos) {
0324     Vector3 loc = ftransform * pos;
0325     return Vector2(loc.x(), loc.y());
0326   };
0327   auto localToGlobal = [itransform](const Vector2& loc) {
0328     return itransform * Vector3(loc.x(), loc.y(), 0.);
0329   };
0330   // Build the grid
0331   std::unique_ptr<SurfaceArray::ISurfaceGridLookup> sl;
0332
0333   // Axis along the binning
0334   switch (bValue) {
0335     case BinningValue::binX: {
0336       ProtoAxis pAxis1 = createEquidistantAxis(gctx, surfacesRaw, binY,
0337                                                protoLayer, ftransform, bins1);
0338       ProtoAxis pAxis2 = createEquidistantAxis(gctx, surfacesRaw, binZ,
0339                                                protoLayer, ftransform, bins2);
0340       sl = makeSurfaceGridLookup2D<detail::AxisBoundaryType::Bound,
0341                                    detail::AxisBoundaryType::Bound>(
0342           globalToLocal, localToGlobal, pAxis1, pAxis2);
0343       break;
0344     }
0345     case BinningValue::binY: {
0346       ProtoAxis pAxis1 = createEquidistantAxis(gctx, surfacesRaw, binX,
0347                                                protoLayer, ftransform, bins1);
0348       ProtoAxis pAxis2 = createEquidistantAxis(gctx, surfacesRaw, binZ,
0349                                                protoLayer, ftransform, bins2);
0350       sl = makeSurfaceGridLookup2D<detail::AxisBoundaryType::Bound,
0351                                    detail::AxisBoundaryType::Bound>(
0352           globalToLocal, localToGlobal, pAxis1, pAxis2);
0353       break;
0354     }
0355     case BinningValue::binZ: {
0356       ProtoAxis pAxis1 = createEquidistantAxis(gctx, surfacesRaw, binX,
0357                                                protoLayer, ftransform, bins1);
0358       ProtoAxis pAxis2 = createEquidistantAxis(gctx, surfacesRaw, binY,
0359                                                protoLayer, ftransform, bins2);
0360       sl = makeSurfaceGridLookup2D<detail::AxisBoundaryType::Bound,
0361                                    detail::AxisBoundaryType::Bound>(
0362           globalToLocal, localToGlobal, pAxis1, pAxis2);
0363       break;
0364     }
0365     default: {
0366       throw std::invalid_argument(
0367           "Acts::SurfaceArrayCreator::"
0368           "surfaceArrayOnPlane: Invalid binning "
0369           "direction");
0370     }
0371   }
0372
0373   sl->fill(gctx, surfacesRaw);
0374   completeBinning(gctx, *sl, surfacesRaw);
0375
0376   return std::make_unique<SurfaceArray>(std::move(sl), std::move(surfaces),
0377                                         ftransform);
0378   //!< @todo implement - take from ATLAS complex TRT builder
0379 }
0380
0381 std::vector<const Acts::Surface*> Acts::SurfaceArrayCreator::findKeySurfaces(
0382     const std::vector<const Surface*>& surfaces,
0383     const std::function<bool(const Surface*, const Surface*)>& equal) const {
0384   std::vector<const Surface*> keys;
0385   for (const auto& srfA : surfaces) {
0386     bool exists = false;
0387     for (const auto& srfB : keys) {
0388       if (equal(srfA, srfB)) {
0389         exists = true;
0390         break;
0391       }
0392     }
0393     if (!exists) {
0394       keys.push_back(srfA);
0395     }
0396   }
0397
0398   return keys;
0399 }
0400
0401 std::size_t Acts::SurfaceArrayCreator::determineBinCount(
0402     const GeometryContext& gctx, const std::vector<const Surface*>& surfaces,
0403     BinningValue bValue) const {
0404   auto matcher = m_cfg.surfaceMatcher;
0405   auto equal = [&gctx, &bValue, &matcher](const Surface* a, const Surface* b) {
0406     return matcher(gctx, bValue, a, b);
0407   };
0408   std::vector<const Surface*> keys = findKeySurfaces(surfaces, equal);
0409
0410   return keys.size();
0411 }
0412
0413 Acts::SurfaceArrayCreator::ProtoAxis
0414 Acts::SurfaceArrayCreator::createVariableAxis(
0415     const GeometryContext& gctx, const std::vector<const Surface*>& surfaces,
0416     BinningValue bValue, const ProtoLayer& protoLayer,
0417     Transform3& transform) const {
0418   if (surfaces.empty()) {
0419     throw std::logic_error(
0420         "No surfaces handed over for creating arbitrary bin utility!");
0421   }
0422   // BinningOption is open for z and r, in case of phi binning reset later
0423   // the vector with the binning Values (boundaries for each bin)
0424
0425   // bind matcher with binning type
0426   auto matcher = m_cfg.surfaceMatcher;
0427   // find the key surfaces
0428   auto equal = [&gctx, &bValue, &matcher](const Surface* a, const Surface* b) {
0429     return matcher(gctx, bValue, a, b);
0430   };
0431   std::vector<const Acts::Surface*> keys = findKeySurfaces(surfaces, equal);
0432
0433   std::vector<AxisScalar> bValues;
0434   if (bValue == Acts::binPhi) {
0435     std::stable_sort(keys.begin(), keys.end(),
0436                      [&gctx](const Acts::Surface* a, const Acts::Surface* b) {
0437                        return (phi(a->binningPosition(gctx, binPhi)) <
0438                                phi(b->binningPosition(gctx, binPhi)));
0439                      });
0440
0441     AxisScalar maxPhi = 0.5 * (phi(keys.at(0)->binningPosition(gctx, binPhi)) +
0442                                phi(keys.at(1)->binningPosition(gctx, binPhi)));
0443
0444     // create rotation, so that maxPhi is +pi
0445     AxisScalar angle = -(M_PI + maxPhi);
0446     transform = (transform)*AngleAxis3(angle, Vector3::UnitZ());
0447
0448     // iterate over all key surfaces, and use their mean position as bValues,
0449     // but
0450     // rotate using transform from before
0451     AxisScalar previous = phi(keys.at(0)->binningPosition(gctx, binPhi));
0452     // go through key surfaces
0453     for (std::size_t i = 1; i < keys.size(); i++) {
0454       const Surface* surface = keys.at(i);
0455       // create central binning values which is the mean of the center
0456       // positions in the binning direction of the current and previous
0457       // surface
0458       AxisScalar edge =
0459           0.5 * (previous + phi(surface->binningPosition(gctx, binPhi))) +
0460           angle;
0461       bValues.push_back(edge);
0462       previous = phi(surface->binningPosition(gctx, binPhi));
0463     }
0464
0465     // segments
0466     unsigned int segments = 72;
0467
0468     // get the bounds of the last surfaces
0469     const Acts::Surface* backSurface = keys.back();
0470     const Acts::PlanarBounds* backBounds =
0471         dynamic_cast<const Acts::PlanarBounds*>(&(backSurface->bounds()));
0472     if (backBounds == nullptr) {
0473       ACTS_ERROR(
0474           "Given SurfaceBounds are not planar - not implemented for "
0475           "other bounds yet! ");
0476     }
0477     // get the global vertices
0478     std::vector<Acts::Vector3> backVertices =
0479         makeGlobalVertices(gctx, *backSurface, backBounds->vertices(segments));
0480     AxisScalar maxBValue = phi(
0481         *std::max_element(backVertices.begin(), backVertices.end(),
0482                           [](const Acts::Vector3& a, const Acts::Vector3& b) {
0483                             return phi(a) < phi(b);
0484                           }));
0485
0486     bValues.push_back(maxBValue);
0487
0488     bValues.push_back(M_PI);
0489
0490   } else if (bValue == Acts::binZ) {
0491     std::stable_sort(keys.begin(), keys.end(),
0492                      [&gctx](const Acts::Surface* a, const Acts::Surface* b) {
0493                        return (a->binningPosition(gctx, binZ).z() <
0494                                b->binningPosition(gctx, binZ).z());
0495                      });
0496
0497     bValues.push_back(protoLayer.min(binZ));
0498     bValues.push_back(protoLayer.max(binZ));
0499
0500     // the z-center position of the previous surface
0501     AxisScalar previous = keys.front()->binningPosition(gctx, binZ).z();
0502     // go through key surfaces
0503     for (auto surface = keys.begin() + 1; surface != keys.end(); surface++) {
0504       // create central binning values which is the mean of the center
0505       // positions in the binning direction of the current and previous
0506       // surface
0507       bValues.push_back(
0508           0.5 * (previous + (*surface)->binningPosition(gctx, binZ).z()));
0509       previous = (*surface)->binningPosition(gctx, binZ).z();
0510     }
0511   } else {  // binR
0512     std::stable_sort(keys.begin(), keys.end(),
0513                      [&gctx](const Acts::Surface* a, const Acts::Surface* b) {
0514                        return (perp(a->binningPosition(gctx, binR)) <
0515                                perp(b->binningPosition(gctx, binR)));
0516                      });
0517
0518     bValues.push_back(protoLayer.min(binR));
0519     bValues.push_back(protoLayer.max(binR));
0520
0521     // the r-center position of the previous surface
0522     AxisScalar previous = perp(keys.front()->binningPosition(gctx, binR));
0523
0524     // go through key surfaces
0525     for (auto surface = keys.begin() + 1; surface != keys.end(); surface++) {
0526       // create central binning values which is the mean of the center
0527       // positions in the binning direction of the current and previous
0528       // surface
0529       bValues.push_back(
0530           0.5 * (previous + perp((*surface)->binningPosition(gctx, binR))));
0531       previous = perp((*surface)->binningPosition(gctx, binR));
0532     }
0533   }
0534   std::sort(bValues.begin(), bValues.end());
0535   ACTS_VERBOSE("Create variable binning Axis for binned SurfaceArray");
0536   ACTS_VERBOSE("    BinningValue: " << bValue);
0537   ACTS_VERBOSE(
0538       " (binX = 0, binY = 1, binZ = 2, binR = 3, binPhi = 4, "
0539       "binRPhi = 5, binH = 6, binEta = 7)");
0540   ACTS_VERBOSE("    Number of bins: " << (bValues.size() - 1));
0541   ACTS_VERBOSE("    (Min/Max) = (" << bValues.front() << "/"
0542                                        << bValues.back() << ")");
0543
0544   ProtoAxis pAxis;
0545   pAxis.bType = arbitrary;
0546   pAxis.bValue = bValue;
0547   pAxis.binEdges = bValues;
0548   pAxis.nBins = bValues.size() - 1;
0549
0550   return pAxis;
0551 }
0552
0553 Acts::SurfaceArrayCreator::ProtoAxis
0554 Acts::SurfaceArrayCreator::createEquidistantAxis(
0555     const GeometryContext& gctx, const std::vector<const Surface*>& surfaces,
0556     BinningValue bValue, const ProtoLayer& protoLayer, Transform3& transform,
0557     std::size_t nBins) const {
0558   if (surfaces.empty()) {
0559     throw std::logic_error(
0560         "No surfaces handed over for creating equidistant axis!");
0561   }
0562   // check the binning type first
0563
0564   double minimum = 0.;
0565   double maximum = 0.;
0566
0567   // binning option is open for z and r, in case of phi binning reset later
0568   // Acts::BinningOption bOption = Acts::open;
0569
0570   // the key surfaces - placed in different bins in the given binning
0571   // direction
0572   std::vector<const Acts::Surface*> keys;
0573
0574   std::size_t binNumber = 0;
0575   if (nBins == 0) {
0576     // determine bin count
0577     binNumber = determineBinCount(gctx, surfaces, bValue);
0578   } else {
0579     // use bin count
0580     binNumber = nBins;
0581   }
0582
0583   // bind matcher & context with binning type
0584   auto matcher = m_cfg.surfaceMatcher;
0585
0586   // now check the binning value
0587   if (bValue == binPhi) {
0588     if (m_cfg.doPhiBinningOptimization) {
0589       // Phi binning
0590       // set the binning option for phi
0591       // sort first in phi
0592       const Acts::Surface* maxElem = *std::max_element(
0593           surfaces.begin(), surfaces.end(),
0594           [&gctx](const Acts::Surface* a, const Acts::Surface* b) {
0595             return phi(a->binningPosition(gctx, binR)) <
0596                    phi(b->binningPosition(gctx, binR));
0597           });
0598
0599       // get the key surfaces at the different phi positions
0600       auto equal = [&gctx, &bValue, &matcher](const Surface* a,
0601                                               const Surface* b) {
0602         return matcher(gctx, bValue, a, b);
0603       };
0604       keys = findKeySurfaces(surfaces, equal);
0605
0606       // multiple surfaces, we bin from -pi to pi closed
0607       if (keys.size() > 1) {
0608         // bOption = Acts::closed;
0609
0610         minimum = -M_PI;
0611         maximum = M_PI;
0612
0613         // double step = 2 * M_PI / keys.size();
0614         double step = 2 * M_PI / binNumber;
0615         // rotate to max phi module plus one half step
0616         // this should make sure that phi wrapping at +- pi
0617         // never falls on a module center
0618         double max = phi(maxElem->binningPosition(gctx, binR));
0619         double angle = M_PI - (max + 0.5 * step);
0620
0621         // replace given transform ref
0622         transform = (transform)*AngleAxis3(angle, Vector3::UnitZ());
0623
0624       } else {
0625         minimum = protoLayer.min(binPhi, true);
0626         maximum = protoLayer.max(binPhi, true);
0627         // we do not need a transform in this case
0628       }
0629     } else {
0630       minimum = -M_PI;
0631       maximum = M_PI;
0632     }
0633   } else {
0634     maximum = protoLayer.max(bValue, false);
0635     minimum = protoLayer.min(bValue, false);
0636   }
0637
0638   // assign the bin size
0639   ACTS_VERBOSE("Create equidistant binning Axis for binned SurfaceArray");
0640   ACTS_VERBOSE("    BinningValue: " << bValue);
0641   ACTS_VERBOSE(
0642       " (binX = 0, binY = 1, binZ = 2, binR = 3, binPhi = 4, "
0643       "binRPhi = 5, binH = 6, binEta = 7)");
0644   ACTS_VERBOSE("    Number of bins: " << binNumber);
0645   ACTS_VERBOSE("    (Min/Max) = (" << minimum << "/" << maximum << ")");
0646
0647   ProtoAxis pAxis;
0648   pAxis.max = maximum;
0649   pAxis.min = minimum;
0650   pAxis.bType = equidistant;
0651   pAxis.bValue = bValue;
0652   pAxis.nBins = binNumber;
0653
0654   return pAxis;
0655 }
0656
0657 std::vector<Acts::Vector3> Acts::SurfaceArrayCreator::makeGlobalVertices(
0658     const GeometryContext& gctx, const Acts::Surface& surface,
0659     const std::vector<Acts::Vector2>& locVertices) const {
0660   std::vector<Acts::Vector3> globVertices;
0661   for (auto& vertex : locVertices) {
0662     Acts::Vector3 globVertex =
0663         surface.localToGlobal(gctx, vertex, Acts::Vector3());
0664     globVertices.push_back(globVertex);
0665   }
0666   return globVertices;
0667 }