File indexing completed on 2025-10-28 07:53:51
0001
0002
0003
0004
0005
0006
0007
0008
0009 #pragma once
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Geometry/GeometryContext.hpp"
0013 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0014 #include "Acts/Surfaces/CylinderBounds.hpp"
0015 #include "Acts/Surfaces/RegularSurface.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017 #include "Acts/Utilities/AnyGridView.hpp"
0018 #include "Acts/Utilities/AxisDefinitions.hpp"
0019 #include "Acts/Utilities/Grid.hpp"
0020 #include "Acts/Utilities/IAxis.hpp"
0021 #include "Acts/Utilities/Intersection.hpp"
0022
0023 #include <iostream>
0024 #include <limits>
0025 #include <vector>
0026
0027 namespace Acts {
0028
0029 using SurfaceVector = std::vector<const Surface*>;
0030
0031
0032
0033
0034
0035
0036
0037 class SurfaceArray {
0038 public:
0039
0040 struct ISurfaceGridLookup {
0041
0042
0043
0044 virtual void fill(const GeometryContext& gctx,
0045 const SurfaceVector& surfaces) = 0;
0046
0047
0048
0049
0050
0051
0052 virtual const SurfaceVector& lookup(const Vector3& position,
0053 const Vector3& direction) const = 0;
0054
0055
0056
0057
0058
0059 virtual SurfaceVector& lookup(std::size_t bin) = 0;
0060
0061
0062
0063
0064
0065 virtual const SurfaceVector& lookup(std::size_t bin) const = 0;
0066
0067
0068
0069
0070
0071
0072 virtual const SurfaceVector& neighbors(const Vector3& position,
0073 const Vector3& direction) const = 0;
0074
0075
0076
0077
0078 virtual std::size_t size() const = 0;
0079
0080
0081
0082
0083 virtual Vector3 getBinCenter(std::size_t bin) const = 0;
0084
0085
0086
0087
0088 virtual std::vector<const IAxis*> getAxes() const = 0;
0089
0090
0091
0092 virtual std::optional<AnyGridConstView<SurfaceVector>> getGridView()
0093 const = 0;
0094
0095
0096
0097 virtual const Surface* surfaceRepresentation() const = 0;
0098
0099
0100
0101
0102
0103
0104 virtual bool isValidBin(std::size_t bin) const = 0;
0105
0106
0107
0108
0109 virtual std::vector<AxisDirection> binningValues() const { return {}; };
0110
0111
0112 virtual ~ISurfaceGridLookup() = 0;
0113 };
0114
0115
0116
0117 template <class Axis1, class Axis2>
0118 struct SurfaceGridLookup : ISurfaceGridLookup {
0119 using Grid_t = Grid<SurfaceVector, Axis1, Axis2>;
0120
0121
0122
0123
0124
0125
0126 SurfaceGridLookup(std::shared_ptr<RegularSurface> representative,
0127 double tolerance, std::tuple<Axis1, Axis2> axes,
0128 std::vector<AxisDirection> bValues = {})
0129 : m_representative(std::move(representative)),
0130 m_tolerance(tolerance),
0131 m_grid(std::move(axes)),
0132 m_binValues(std::move(bValues)) {
0133 m_neighborMap.resize(m_grid.size());
0134 }
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145 void fill(const GeometryContext& gctx,
0146 const SurfaceVector& surfaces) override {
0147 for (const Surface* surface : surfaces) {
0148 const std::size_t globalBin = fillSurfaceToBinMapping(gctx, *surface);
0149 if (globalBin == 0) {
0150 continue;
0151 }
0152
0153 fillBinToSurfaceMapping(gctx, *surface, globalBin);
0154 }
0155
0156 populateNeighborCache();
0157 }
0158
0159 const SurfaceVector& lookup(const Vector3& position,
0160 const Vector3& direction) const override {
0161 return m_grid.at(findGlobalBin(position, direction,
0162 std::numeric_limits<double>::infinity()));
0163 }
0164
0165
0166
0167
0168
0169 SurfaceVector& lookup(std::size_t bin) override { return m_grid.at(bin); }
0170
0171
0172
0173
0174
0175 const SurfaceVector& lookup(std::size_t bin) const override {
0176 return m_grid.at(bin);
0177 }
0178
0179
0180
0181
0182
0183
0184 const SurfaceVector& neighbors(const Vector3& position,
0185 const Vector3& direction) const override {
0186 return m_neighborMap.at(findGlobalBin(
0187 position, direction, std::numeric_limits<double>::infinity()));
0188 }
0189
0190
0191
0192
0193 std::size_t size() const override { return m_grid.size(); }
0194
0195
0196
0197
0198 std::vector<AxisDirection> binningValues() const override {
0199 return m_binValues;
0200 }
0201
0202
0203
0204
0205 Vector3 getBinCenter(std::size_t bin) const override {
0206 GeometryContext gctx;
0207 return getBinCenterImpl(gctx, bin);
0208 }
0209
0210
0211
0212
0213 std::vector<const IAxis*> getAxes() const override {
0214 auto arr = m_grid.axes();
0215 return std::vector<const IAxis*>(arr.begin(), arr.end());
0216 }
0217
0218 std::optional<AnyGridConstView<SurfaceVector>> getGridView()
0219 const override {
0220 return AnyGridConstView<SurfaceVector>{m_grid};
0221 }
0222
0223 const Surface* surfaceRepresentation() const override {
0224 return m_representative.get();
0225 }
0226
0227
0228
0229
0230
0231
0232 bool isValidBin(std::size_t bin) const override {
0233 std::array<std::size_t, 2> indices = m_grid.localBinsFromGlobalBin(bin);
0234 std::array<std::size_t, 2> nBins = m_grid.numLocalBins();
0235 for (std::size_t i = 0; i < indices.size(); ++i) {
0236 std::size_t idx = indices.at(i);
0237 if (idx <= 0 || idx >= nBins.at(i) + 1) {
0238 return false;
0239 }
0240 }
0241 return true;
0242 }
0243
0244 private:
0245
0246 std::size_t fillSurfaceToBinMapping(const GeometryContext& gctx,
0247 const Surface& surface) {
0248 const Vector3 pos = surface.referencePosition(gctx, AxisDirection::AxisR);
0249 const Vector3 normal = m_representative->normal(gctx, pos);
0250 const std::size_t globalBin = findGlobalBin(pos, normal, m_tolerance);
0251 if (globalBin != 0) {
0252 m_grid.at(globalBin).push_back(&surface);
0253 }
0254 return globalBin;
0255 };
0256
0257
0258 void fillBinToSurfaceMapping(const GeometryContext& gctx,
0259 const Surface& surface, std::size_t startBin) {
0260 const std::array<std::size_t, 2> startIndices =
0261 m_grid.localBinsFromGlobalBin(startBin);
0262 const auto startNeighborIndices =
0263 m_grid.neighborHoodIndices(startIndices, 1u);
0264
0265 std::set<std::size_t> visited({startBin});
0266 std::vector<std::size_t> queue(startNeighborIndices.begin(),
0267 startNeighborIndices.end());
0268
0269 while (!queue.empty()) {
0270 const std::size_t current = queue.back();
0271 queue.pop_back();
0272 if (visited.contains(current)) {
0273 continue;
0274 }
0275
0276 const std::array<std::size_t, 2> currentIndices =
0277 m_grid.localBinsFromGlobalBin(current);
0278 visited.insert(current);
0279
0280 const std::array<double, 2> gridLocal =
0281 m_grid.binCenter(currentIndices);
0282 const Vector2 surfaceLocal = gridToSurfaceLocal(gridLocal);
0283 const Vector3 normal = m_representative->normal(gctx, surfaceLocal);
0284 const Vector3 global =
0285 m_representative->localToGlobal(gctx, surfaceLocal, normal);
0286
0287 const Intersection3D intersection =
0288 surface.intersect(gctx, global, normal, BoundaryTolerance::None())
0289 .closest();
0290 if (!intersection.isValid() ||
0291 std::abs(intersection.pathLength()) > m_tolerance) {
0292 continue;
0293 }
0294 m_grid.at(current).push_back(&surface);
0295
0296 const auto neighborIndices =
0297 m_grid.neighborHoodIndices(currentIndices, 1u);
0298 queue.insert(queue.end(), neighborIndices.begin(),
0299 neighborIndices.end());
0300 }
0301 };
0302
0303 void populateNeighborCache() {
0304
0305 for (std::size_t i = 0; i < m_grid.size(); i++) {
0306 if (!isValidBin(i)) {
0307 continue;
0308 }
0309 const std::array<std::size_t, 2> indices =
0310 m_grid.localBinsFromGlobalBin(i);
0311 std::vector<const Surface*>& neighbors = m_neighborMap.at(i);
0312 neighbors.clear();
0313
0314 for (std::size_t idx : m_grid.neighborHoodIndices(indices, 1u)) {
0315 const std::vector<const Surface*>& binContent = m_grid.at(idx);
0316 std::copy(binContent.begin(), binContent.end(),
0317 std::back_inserter(neighbors));
0318 }
0319
0320 std::ranges::sort(neighbors);
0321 auto last = std::ranges::unique(neighbors);
0322 neighbors.erase(last.begin(), last.end());
0323 neighbors.shrink_to_fit();
0324 }
0325 }
0326
0327 Vector3 getBinCenterImpl(const GeometryContext& gctx,
0328 std::size_t bin) const {
0329 const std::array<double, 2> gridLocal =
0330 m_grid.binCenter(m_grid.localBinsFromGlobalBin(bin));
0331 const Vector2 surfaceLocal = gridToSurfaceLocal(gridLocal);
0332 return m_representative->localToGlobal(gctx, surfaceLocal);
0333 }
0334
0335 const CylinderBounds* getCylinderBounds() const {
0336 return dynamic_cast<const CylinderBounds*>(&m_representative->bounds());
0337 }
0338
0339 Vector2 gridToSurfaceLocal(std::array<double, 2> gridLocal) const {
0340 Vector2 surfaceLocal = Eigen::Map<Vector2>(gridLocal.data());
0341 if (const CylinderBounds* bounds = getCylinderBounds();
0342 bounds != nullptr) {
0343 surfaceLocal[0] *= bounds->get(CylinderBounds::eR);
0344 }
0345 return surfaceLocal;
0346 }
0347 std::array<double, 2> surfaceToGridLocal(Vector2 local) const {
0348 std::array<double, 2> gridLocal = {local[0], local[1]};
0349 if (const CylinderBounds* bounds = getCylinderBounds();
0350 bounds != nullptr) {
0351 gridLocal[0] /= bounds->get(CylinderBounds::eR);
0352 }
0353 return gridLocal;
0354 }
0355
0356 std::size_t findGlobalBin(const Vector3& position, const Vector3& direction,
0357 double tolerance) const {
0358 GeometryContext gctx;
0359
0360 const Intersection3D intersection =
0361 m_representative
0362 ->intersect(gctx, position, direction,
0363 BoundaryTolerance::Infinite())
0364 .closest();
0365 if (!intersection.isValid() ||
0366 std::abs(intersection.pathLength()) > tolerance) {
0367 return 0;
0368 }
0369 const Vector2 surfaceLocal =
0370 m_representative
0371 ->globalToLocal(gctx, intersection.position(), direction)
0372 .value();
0373 const std::array<double, 2> gridLocal = surfaceToGridLocal(surfaceLocal);
0374 return m_grid.globalBinFromPosition(gridLocal);
0375 }
0376
0377 std::shared_ptr<RegularSurface> m_representative;
0378 double m_tolerance{};
0379 Grid_t m_grid;
0380 std::vector<AxisDirection> m_binValues;
0381 std::vector<SurfaceVector> m_neighborMap;
0382 };
0383
0384
0385
0386 struct SingleElementLookup : ISurfaceGridLookup {
0387
0388
0389 explicit SingleElementLookup(SurfaceVector::value_type element)
0390 : m_element({element}) {}
0391
0392
0393
0394 explicit SingleElementLookup(const SurfaceVector& elements)
0395 : m_element(elements) {}
0396
0397
0398
0399 const SurfaceVector& lookup(const Vector3& ,
0400 const Vector3& ) const override {
0401 return m_element;
0402 }
0403
0404
0405
0406 SurfaceVector& lookup(std::size_t ) override { return m_element; }
0407
0408
0409
0410 const SurfaceVector& lookup(std::size_t ) const override {
0411 return m_element;
0412 }
0413
0414
0415
0416 const SurfaceVector& neighbors(
0417 const Vector3& ,
0418 const Vector3& ) const override {
0419 return m_element;
0420 }
0421
0422
0423
0424 std::size_t size() const override { return 1; }
0425
0426
0427
0428 Vector3 getBinCenter(std::size_t ) const override {
0429 return Vector3(0, 0, 0);
0430 }
0431
0432
0433
0434 std::vector<const IAxis*> getAxes() const override { return {}; }
0435
0436 std::optional<AnyGridConstView<SurfaceVector>> getGridView()
0437 const override {
0438 return std::nullopt;
0439 }
0440
0441 const Surface* surfaceRepresentation() const override { return nullptr; }
0442
0443
0444
0445 void fill(const GeometryContext& ,
0446 const SurfaceVector& ) override {}
0447
0448
0449
0450 bool isValidBin(std::size_t ) const override { return true; }
0451
0452 private:
0453 SurfaceVector m_element;
0454 };
0455
0456
0457
0458
0459
0460
0461
0462
0463 explicit SurfaceArray(std::unique_ptr<ISurfaceGridLookup> gridLookup,
0464 std::vector<std::shared_ptr<const Surface>> surfaces,
0465 const Transform3& transform = Transform3::Identity());
0466
0467
0468
0469 explicit SurfaceArray(std::shared_ptr<const Surface> srf);
0470
0471
0472
0473
0474
0475
0476 const SurfaceVector& at(const Vector3& position,
0477 const Vector3& direction) const {
0478 return p_gridLookup->lookup(position, direction);
0479 }
0480
0481
0482
0483
0484 SurfaceVector& at(std::size_t bin) { return p_gridLookup->lookup(bin); }
0485
0486
0487
0488
0489 const SurfaceVector& at(std::size_t bin) const {
0490 return p_gridLookup->lookup(bin);
0491 }
0492
0493
0494
0495
0496
0497
0498
0499
0500 const SurfaceVector& neighbors(const Vector3& position,
0501 const Vector3& direction) const {
0502 return p_gridLookup->neighbors(position, direction);
0503 }
0504
0505
0506
0507
0508 std::size_t size() const { return p_gridLookup->size(); }
0509
0510
0511
0512
0513 Vector3 getBinCenter(std::size_t bin) const {
0514 return p_gridLookup->getBinCenter(bin);
0515 }
0516
0517
0518
0519
0520
0521
0522 const SurfaceVector& surfaces() const { return m_surfacesRawPointers; }
0523
0524
0525
0526
0527
0528 std::vector<const IAxis*> getAxes() const { return p_gridLookup->getAxes(); }
0529
0530
0531
0532
0533
0534
0535 bool isValidBin(std::size_t bin) const {
0536 return p_gridLookup->isValidBin(bin);
0537 }
0538
0539
0540
0541 const Transform3& transform() const { return m_transform; }
0542
0543
0544
0545
0546 std::vector<AxisDirection> binningValues() const {
0547 return p_gridLookup->binningValues();
0548 };
0549
0550
0551
0552
0553
0554 std::ostream& toStream(const GeometryContext& gctx, std::ostream& sl) const;
0555
0556
0557
0558 const ISurfaceGridLookup& gridLookup() const { return *p_gridLookup; }
0559
0560 private:
0561 std::unique_ptr<ISurfaceGridLookup> p_gridLookup;
0562
0563 std::vector<std::shared_ptr<const Surface>> m_surfaces;
0564
0565
0566 SurfaceVector m_surfacesRawPointers;
0567
0568
0569 Transform3 m_transform;
0570 };
0571
0572 }