src/Geometry/SurfaceArrayCreator.cpp

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