Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-10-21 08:02:37 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!