Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-05-28 07:11:21 

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