File indexing completed on 2025-01-18 09:27:55
0001
0002
0003
0004
0005
0006
0007
0008
0009 #pragma once
0010
0011 #include "Acts/Utilities/IAxis.hpp"
0012 #include "Acts/Utilities/detail/AxisFwd.hpp"
0013
0014 #include <algorithm>
0015 #include <cmath>
0016 #include <vector>
0017
0018 namespace Acts::detail {
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030 class NeighborHoodIndices {
0031 public:
0032 NeighborHoodIndices() = default;
0033
0034 NeighborHoodIndices(std::size_t begin, std::size_t end)
0035 : m_begin1(begin), m_end1(end), m_begin2(end), m_end2(end) {}
0036
0037 NeighborHoodIndices(std::size_t begin1, std::size_t end1, std::size_t begin2,
0038 std::size_t end2)
0039 : m_begin1(begin1), m_end1(end1), m_begin2(begin2), m_end2(end2) {}
0040
0041 class iterator {
0042 public:
0043 iterator() = default;
0044
0045
0046 iterator(std::size_t current) : m_current(current), m_wrapped(true) {}
0047
0048 iterator(std::size_t begin1, std::size_t end1, std::size_t begin2)
0049 : m_current(begin1),
0050 m_end1(end1),
0051 m_begin2(begin2),
0052 m_wrapped(begin1 == begin2) {}
0053
0054 std::size_t operator*() const { return m_current; }
0055
0056 iterator& operator++() {
0057 ++m_current;
0058 if (m_current == m_end1) {
0059 m_current = m_begin2;
0060 m_wrapped = true;
0061 }
0062 return *this;
0063 }
0064
0065 bool operator==(const iterator& it) const {
0066 return (m_current == it.m_current) && (m_wrapped == it.m_wrapped);
0067 }
0068
0069 bool operator!=(const iterator& it) const { return !(*this == it); }
0070
0071 private:
0072 std::size_t m_current = 0, m_end1 = 0, m_begin2 = 0;
0073 bool m_wrapped = false;
0074 };
0075
0076 iterator begin() const { return iterator(m_begin1, m_end1, m_begin2); }
0077
0078 iterator end() const { return iterator(m_end2); }
0079
0080
0081 std::size_t size() const { return (m_end1 - m_begin1) + (m_end2 - m_begin2); }
0082
0083
0084 std::vector<std::size_t> collect() const {
0085 std::vector<std::size_t> result;
0086 result.reserve(this->size());
0087 for (std::size_t idx : *this) {
0088 result.push_back(idx);
0089 }
0090 return result;
0091 }
0092
0093 private:
0094 std::size_t m_begin1 = 0, m_end1 = 0, m_begin2 = 0, m_end2 = 0;
0095 };
0096
0097
0098
0099
0100
0101 template <AxisBoundaryType bdt>
0102 class Axis<AxisType::Equidistant, bdt> final : public IAxis {
0103 public:
0104
0105
0106
0107
0108
0109
0110
0111
0112 Axis(ActsScalar xmin, ActsScalar xmax, std::size_t nBins)
0113 : m_min(xmin),
0114 m_max(xmax),
0115 m_width((xmax - xmin) / nBins),
0116 m_bins(nBins) {}
0117
0118
0119
0120
0121 bool isEquidistant() const override { return true; }
0122
0123
0124
0125
0126 bool isVariable() const override { return false; }
0127
0128
0129
0130
0131 AxisBoundaryType getBoundaryType() const override { return bdt; }
0132
0133
0134
0135
0136
0137
0138
0139
0140 NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0141 std::size_t size = 1) const {
0142 return neighborHoodIndices(idx, std::make_pair(-size, size));
0143 }
0144
0145
0146
0147
0148
0149
0150
0151
0152
0153
0154
0155 template <AxisBoundaryType T = bdt,
0156 std::enable_if_t<T == AxisBoundaryType::Open, int> = 0>
0157 NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0158 std::pair<int, int> sizes = {
0159 -1, 1}) const {
0160 constexpr int min = 0;
0161 const int max = getNBins() + 1;
0162 const int itmin = std::clamp(static_cast<int>(idx + sizes.first), min, max);
0163 const int itmax =
0164 std::clamp(static_cast<int>(idx + sizes.second), min, max);
0165 return NeighborHoodIndices(itmin, itmax + 1);
0166 }
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177 template <AxisBoundaryType T = bdt,
0178 std::enable_if_t<T == AxisBoundaryType::Bound, int> = 0>
0179 NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0180 std::pair<int, int> sizes = {
0181 -1, 1}) const {
0182 if (idx <= 0 || idx >= (getNBins() + 1)) {
0183 return NeighborHoodIndices();
0184 }
0185 constexpr int min = 1;
0186 const int max = getNBins();
0187 const int itmin = std::clamp(static_cast<int>(idx) + sizes.first, min, max);
0188 const int itmax =
0189 std::clamp(static_cast<int>(idx) + sizes.second, min, max);
0190 return NeighborHoodIndices(itmin, itmax + 1);
0191 }
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202 template <AxisBoundaryType T = bdt,
0203 std::enable_if_t<T == AxisBoundaryType::Closed, int> = 0>
0204 NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0205 std::pair<int, int> sizes = {
0206 -1, 1}) const {
0207
0208 if (idx <= 0 || idx >= (getNBins() + 1)) {
0209 return NeighborHoodIndices();
0210 }
0211
0212
0213
0214
0215 const int max = getNBins();
0216 sizes.first = std::clamp(sizes.first, -max, max);
0217 sizes.second = std::clamp(sizes.second, -max, max);
0218 if (std::abs(sizes.first - sizes.second) >= max) {
0219 sizes.first = 1 - idx;
0220 sizes.second = max - idx;
0221 }
0222
0223
0224
0225
0226
0227
0228
0229
0230 const int itmin = idx + sizes.first;
0231 const int itmax = idx + sizes.second;
0232 const std::size_t itfirst = wrapBin(itmin);
0233 const std::size_t itlast = wrapBin(itmax);
0234 if (itfirst <= itlast) {
0235 return NeighborHoodIndices(itfirst, itlast + 1);
0236 } else {
0237 return NeighborHoodIndices(itfirst, max + 1, 1, itlast + 1);
0238 }
0239 }
0240
0241
0242
0243
0244
0245
0246
0247 template <AxisBoundaryType T = bdt,
0248 std::enable_if_t<T == AxisBoundaryType::Open, int> = 0>
0249 std::size_t wrapBin(int bin) const {
0250 return std::max(std::min(bin, static_cast<int>(getNBins()) + 1), 0);
0251 }
0252
0253
0254
0255
0256
0257
0258
0259 template <AxisBoundaryType T = bdt,
0260 std::enable_if_t<T == AxisBoundaryType::Bound, int> = 0>
0261 std::size_t wrapBin(int bin) const {
0262 return std::max(std::min(bin, static_cast<int>(getNBins())), 1);
0263 }
0264
0265
0266
0267
0268
0269
0270
0271 template <AxisBoundaryType T = bdt,
0272 std::enable_if_t<T == AxisBoundaryType::Closed, int> = 0>
0273 std::size_t wrapBin(int bin) const {
0274 const int w = getNBins();
0275 return 1 + (w + ((bin - 1) % w)) % w;
0276
0277 }
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
0289 std::size_t getBin(ActsScalar x) const {
0290 return wrapBin(
0291 static_cast<int>(std::floor((x - getMin()) / getBinWidth()) + 1));
0292 }
0293
0294
0295
0296
0297 ActsScalar getBinWidth(std::size_t = 0) const { return m_width; }
0298
0299
0300
0301
0302
0303
0304
0305
0306
0307
0308
0309 ActsScalar getBinLowerBound(std::size_t bin) const {
0310 return getMin() + (bin - 1) * getBinWidth();
0311 }
0312
0313
0314
0315
0316
0317
0318
0319
0320
0321
0322
0323 ActsScalar getBinUpperBound(std::size_t bin) const {
0324 return getMin() + bin * getBinWidth();
0325 }
0326
0327
0328
0329
0330
0331
0332
0333
0334 ActsScalar getBinCenter(std::size_t bin) const {
0335 return getMin() + (bin - 0.5) * getBinWidth();
0336 }
0337
0338
0339
0340
0341 ActsScalar getMax() const override { return m_max; }
0342
0343
0344
0345
0346 ActsScalar getMin() const override { return m_min; }
0347
0348
0349
0350
0351 std::size_t getNBins() const override { return m_bins; }
0352
0353
0354
0355
0356
0357
0358
0359
0360 bool isInside(ActsScalar x) const { return (m_min <= x) && (x < m_max); }
0361
0362
0363
0364 std::vector<ActsScalar> getBinEdges() const override {
0365 std::vector<ActsScalar> binEdges;
0366 for (std::size_t i = 1; i <= m_bins; i++) {
0367 binEdges.push_back(getBinLowerBound(i));
0368 }
0369 binEdges.push_back(getBinUpperBound(m_bins));
0370 return binEdges;
0371 }
0372
0373 private:
0374
0375 ActsScalar m_min;
0376
0377 ActsScalar m_max;
0378
0379 ActsScalar m_width;
0380
0381 std::size_t m_bins;
0382 };
0383
0384
0385
0386
0387
0388 template <AxisBoundaryType bdt>
0389 class Axis<AxisType::Variable, bdt> final : public IAxis {
0390 public:
0391
0392
0393
0394
0395
0396
0397
0398
0399
0400 Axis(std::vector<ActsScalar> binEdges) : m_binEdges(std::move(binEdges)) {}
0401
0402
0403
0404
0405 bool isEquidistant() const override { return false; }
0406
0407
0408
0409
0410 bool isVariable() const override { return true; }
0411
0412
0413
0414
0415 AxisBoundaryType getBoundaryType() const override { return bdt; }
0416
0417
0418
0419
0420
0421
0422
0423
0424 NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0425 std::size_t size = 1) const {
0426 return neighborHoodIndices(idx, std::make_pair(-size, size));
0427 }
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439 template <AxisBoundaryType T = bdt,
0440 std::enable_if_t<T == AxisBoundaryType::Open, int> = 0>
0441 NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0442 std::pair<int, int> sizes = {
0443 -1, 1}) const {
0444 constexpr int min = 0;
0445 const int max = getNBins() + 1;
0446 const int itmin = std::max(min, static_cast<int>(idx) + sizes.first);
0447 const int itmax = std::min(max, static_cast<int>(idx) + sizes.second);
0448 return NeighborHoodIndices(itmin, itmax + 1);
0449 }
0450
0451
0452
0453
0454
0455
0456
0457
0458
0459
0460 template <AxisBoundaryType T = bdt,
0461 std::enable_if_t<T == AxisBoundaryType::Bound, int> = 0>
0462 NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0463 std::pair<int, int> sizes = {
0464 -1, 1}) const {
0465 if (idx <= 0 || idx >= (getNBins() + 1)) {
0466 return NeighborHoodIndices();
0467 }
0468 constexpr int min = 1;
0469 const int max = getNBins();
0470 const int itmin = std::max(min, static_cast<int>(idx) + sizes.first);
0471 const int itmax = std::min(max, static_cast<int>(idx) + sizes.second);
0472 return NeighborHoodIndices(itmin, itmax + 1);
0473 }
0474
0475
0476
0477
0478
0479
0480
0481
0482
0483
0484 template <AxisBoundaryType T = bdt,
0485 std::enable_if_t<T == AxisBoundaryType::Closed, int> = 0>
0486 NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0487 std::pair<int, int> sizes = {
0488 -1, 1}) const {
0489
0490 if (idx <= 0 || idx >= (getNBins() + 1)) {
0491 return NeighborHoodIndices();
0492 }
0493
0494
0495
0496
0497 const int max = getNBins();
0498 sizes.first = std::clamp(sizes.first, -max, max);
0499 sizes.second = std::clamp(sizes.second, -max, max);
0500 if (std::abs(sizes.first - sizes.second) >= max) {
0501 sizes.first = 1 - idx;
0502 sizes.second = max - idx;
0503 }
0504
0505
0506
0507
0508
0509
0510
0511
0512 const int itmin = idx + sizes.first;
0513 const int itmax = idx + sizes.second;
0514 const std::size_t itfirst = wrapBin(itmin);
0515 const std::size_t itlast = wrapBin(itmax);
0516 if (itfirst <= itlast) {
0517 return NeighborHoodIndices(itfirst, itlast + 1);
0518 } else {
0519 return NeighborHoodIndices(itfirst, max + 1, 1, itlast + 1);
0520 }
0521 }
0522
0523
0524
0525
0526
0527
0528
0529 template <AxisBoundaryType T = bdt,
0530 std::enable_if_t<T == AxisBoundaryType::Open, int> = 0>
0531 std::size_t wrapBin(int bin) const {
0532 return std::max(std::min(bin, static_cast<int>(getNBins()) + 1), 0);
0533 }
0534
0535
0536
0537
0538
0539
0540
0541 template <AxisBoundaryType T = bdt,
0542 std::enable_if_t<T == AxisBoundaryType::Bound, int> = 0>
0543 std::size_t wrapBin(int bin) const {
0544 return std::max(std::min(bin, static_cast<int>(getNBins())), 1);
0545 }
0546
0547
0548
0549
0550
0551
0552
0553 template <AxisBoundaryType T = bdt,
0554 std::enable_if_t<T == AxisBoundaryType::Closed, int> = 0>
0555 std::size_t wrapBin(int bin) const {
0556 const int w = getNBins();
0557 return 1 + (w + ((bin - 1) % w)) % w;
0558
0559 }
0560
0561
0562
0563
0564
0565
0566
0567
0568
0569
0570
0571 std::size_t getBin(ActsScalar x) const {
0572 const auto it =
0573 std::upper_bound(std::begin(m_binEdges), std::end(m_binEdges), x);
0574 return wrapBin(std::distance(std::begin(m_binEdges), it));
0575 }
0576
0577
0578
0579
0580
0581
0582
0583
0584 ActsScalar getBinWidth(std::size_t bin) const {
0585 return m_binEdges.at(bin) - m_binEdges.at(bin - 1);
0586 }
0587
0588
0589
0590
0591
0592
0593
0594
0595
0596
0597
0598 ActsScalar getBinLowerBound(std::size_t bin) const {
0599 return m_binEdges.at(bin - 1);
0600 }
0601
0602
0603
0604
0605
0606
0607
0608
0609
0610
0611
0612 ActsScalar getBinUpperBound(std::size_t bin) const {
0613 return m_binEdges.at(bin);
0614 }
0615
0616
0617
0618
0619
0620
0621
0622
0623 ActsScalar getBinCenter(std::size_t bin) const {
0624 return 0.5 * (getBinLowerBound(bin) + getBinUpperBound(bin));
0625 }
0626
0627
0628
0629
0630 ActsScalar getMax() const override { return m_binEdges.back(); }
0631
0632
0633
0634
0635 ActsScalar getMin() const override { return m_binEdges.front(); }
0636
0637
0638
0639
0640 std::size_t getNBins() const override { return m_binEdges.size() - 1; }
0641
0642
0643
0644
0645
0646
0647
0648
0649 bool isInside(ActsScalar x) const {
0650 return (m_binEdges.front() <= x) && (x < m_binEdges.back());
0651 }
0652
0653
0654
0655 std::vector<ActsScalar> getBinEdges() const override { return m_binEdges; }
0656
0657 private:
0658
0659 std::vector<ActsScalar> m_binEdges;
0660 };
0661 }