Back to home page

EIC code displayed by LXR

 
 

    


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

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