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