File indexing completed on 2025-06-05 08:29:44
0001
0002
0003
0004
0005
0006
0007
0008
0009 #pragma once
0010
0011 #include "Acts/Utilities/AxisFwd.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 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 static constexpr AxisType type = AxisType::Equidistant;
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114 Axis(ActsScalar xmin, ActsScalar 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, ActsScalar xmin, ActsScalar 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(ActsScalar x) const {
0312 return wrapBin(
0313 static_cast<int>(std::floor((x - getMin()) / getBinWidth()) + 1));
0314 }
0315
0316
0317
0318
0319 ActsScalar getBinWidth(std::size_t = 0) const { return m_width; }
0320
0321
0322
0323
0324
0325
0326
0327
0328
0329
0330
0331 ActsScalar 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 ActsScalar getBinUpperBound(std::size_t bin) const {
0346 return getMin() + bin * getBinWidth();
0347 }
0348
0349
0350
0351
0352
0353
0354
0355
0356 ActsScalar getBinCenter(std::size_t bin) const {
0357 return getMin() + (bin - 0.5) * getBinWidth();
0358 }
0359
0360
0361
0362
0363 ActsScalar getMax() const override { return m_max; }
0364
0365
0366
0367
0368 ActsScalar 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 bool isInside(ActsScalar x) const { return (m_min <= x) && (x < m_max); }
0383
0384
0385
0386 std::vector<ActsScalar> getBinEdges() const override {
0387 std::vector<ActsScalar> binEdges;
0388 for (std::size_t i = 1; i <= m_bins; i++) {
0389 binEdges.push_back(getBinLowerBound(i));
0390 }
0391 binEdges.push_back(getBinUpperBound(m_bins));
0392 return binEdges;
0393 }
0394
0395 friend std::ostream& operator<<(std::ostream& os, const Axis& axis) {
0396 os << "Axis<Equidistant, " << bdt << ">(";
0397 os << axis.m_min << ", ";
0398 os << axis.m_max << ", ";
0399 os << axis.m_bins << ")";
0400 return os;
0401 }
0402
0403 protected:
0404 void toStream(std::ostream& os) const override { os << *this; }
0405
0406 private:
0407
0408 ActsScalar m_min{};
0409
0410 ActsScalar m_max{};
0411
0412 ActsScalar m_width{};
0413
0414 std::size_t m_bins{};
0415 };
0416
0417
0418
0419
0420
0421 template <AxisBoundaryType bdt>
0422 class Axis<AxisType::Variable, bdt> final : public IAxis {
0423 public:
0424 static constexpr AxisType type = AxisType::Variable;
0425
0426
0427
0428
0429
0430
0431
0432
0433 explicit Axis(std::vector<ActsScalar> binEdges)
0434 : m_binEdges(std::move(binEdges)) {}
0435
0436
0437
0438
0439
0440
0441
0442
0443
0444 Axis(AxisBoundaryTypeTag<bdt> typeTag, std::vector<ActsScalar> binEdges)
0445 : Axis(std::move(binEdges)) {
0446 static_cast<void>(typeTag);
0447 }
0448
0449
0450
0451
0452 bool isEquidistant() const override { return false; }
0453
0454
0455
0456
0457 bool isVariable() const override { return true; }
0458
0459
0460
0461 AxisType getType() const override { return type; }
0462
0463
0464
0465
0466 AxisBoundaryType getBoundaryType() const override { return bdt; }
0467
0468
0469
0470
0471
0472
0473
0474
0475 NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0476 std::size_t size = 1) const {
0477 return neighborHoodIndices(idx,
0478 std::make_pair(-static_cast<int>(size), size));
0479 }
0480
0481
0482
0483
0484
0485
0486
0487
0488
0489
0490
0491 NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0492 std::pair<int, int> sizes = {-1,
0493 1}) const
0494 requires(bdt == AxisBoundaryType::Open)
0495 {
0496 constexpr int min = 0;
0497 const int max = getNBins() + 1;
0498 const int itmin = std::max(min, static_cast<int>(idx) + sizes.first);
0499 const int itmax = std::min(max, static_cast<int>(idx) + sizes.second);
0500 return NeighborHoodIndices(itmin, itmax + 1);
0501 }
0502
0503
0504
0505
0506
0507
0508
0509
0510
0511
0512 NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0513 std::pair<int, int> sizes = {-1,
0514 1}) const
0515 requires(bdt == AxisBoundaryType::Bound)
0516 {
0517 if (idx <= 0 || idx >= (getNBins() + 1)) {
0518 return NeighborHoodIndices();
0519 }
0520 constexpr int min = 1;
0521 const int max = getNBins();
0522 const int itmin = std::max(min, static_cast<int>(idx) + sizes.first);
0523 const int itmax = std::min(max, static_cast<int>(idx) + sizes.second);
0524 return NeighborHoodIndices(itmin, itmax + 1);
0525 }
0526
0527
0528
0529
0530
0531
0532
0533
0534
0535
0536 NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0537 std::pair<int, int> sizes = {-1,
0538 1}) const
0539 requires(bdt == AxisBoundaryType::Closed)
0540 {
0541
0542 if (idx <= 0 || idx >= (getNBins() + 1)) {
0543 return NeighborHoodIndices();
0544 }
0545
0546
0547
0548
0549 const int max = getNBins();
0550 sizes.first = std::clamp(sizes.first, -max, max);
0551 sizes.second = std::clamp(sizes.second, -max, max);
0552 if (std::abs(sizes.first - sizes.second) >= max) {
0553 sizes.first = 1 - idx;
0554 sizes.second = max - idx;
0555 }
0556
0557
0558
0559
0560
0561
0562
0563
0564 const int itmin = idx + sizes.first;
0565 const int itmax = idx + sizes.second;
0566 const std::size_t itfirst = wrapBin(itmin);
0567 const std::size_t itlast = wrapBin(itmax);
0568 if (itfirst <= itlast) {
0569 return NeighborHoodIndices(itfirst, itlast + 1);
0570 } else {
0571 return NeighborHoodIndices(itfirst, max + 1, 1, itlast + 1);
0572 }
0573 }
0574
0575
0576
0577
0578
0579
0580
0581 std::size_t wrapBin(int bin) const
0582 requires(bdt == AxisBoundaryType::Open)
0583 {
0584 return std::max(std::min(bin, static_cast<int>(getNBins()) + 1), 0);
0585 }
0586
0587
0588
0589
0590
0591
0592
0593 std::size_t wrapBin(int bin) const
0594 requires(bdt == AxisBoundaryType::Bound)
0595 {
0596 return std::max(std::min(bin, static_cast<int>(getNBins())), 1);
0597 }
0598
0599
0600
0601
0602
0603
0604
0605 std::size_t wrapBin(int bin) const
0606 requires(bdt == AxisBoundaryType::Closed)
0607 {
0608 const int w = getNBins();
0609 return 1 + (w + ((bin - 1) % w)) % w;
0610
0611 }
0612
0613
0614
0615
0616
0617
0618
0619
0620
0621
0622
0623 std::size_t getBin(ActsScalar x) const {
0624 const auto it =
0625 std::upper_bound(std::begin(m_binEdges), std::end(m_binEdges), x);
0626 return wrapBin(std::distance(std::begin(m_binEdges), it));
0627 }
0628
0629
0630
0631
0632
0633
0634
0635
0636 ActsScalar getBinWidth(std::size_t bin) const {
0637 return m_binEdges.at(bin) - m_binEdges.at(bin - 1);
0638 }
0639
0640
0641
0642
0643
0644
0645
0646
0647
0648
0649
0650 ActsScalar getBinLowerBound(std::size_t bin) const {
0651 return m_binEdges.at(bin - 1);
0652 }
0653
0654
0655
0656
0657
0658
0659
0660
0661
0662
0663
0664 ActsScalar getBinUpperBound(std::size_t bin) const {
0665 return m_binEdges.at(bin);
0666 }
0667
0668
0669
0670
0671
0672
0673
0674
0675 ActsScalar getBinCenter(std::size_t bin) const {
0676 return 0.5 * (getBinLowerBound(bin) + getBinUpperBound(bin));
0677 }
0678
0679
0680
0681
0682 ActsScalar getMax() const override { return m_binEdges.back(); }
0683
0684
0685
0686
0687 ActsScalar 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 bool isInside(ActsScalar x) const {
0702 return (m_binEdges.front() <= x) && (x < m_binEdges.back());
0703 }
0704
0705
0706
0707 std::vector<ActsScalar> getBinEdges() const override { return m_binEdges; }
0708
0709 friend std::ostream& operator<<(std::ostream& os, const Axis& axis) {
0710 os << "Axis<Variable, " << bdt << ">(";
0711 os << axis.m_binEdges.front();
0712 for (std::size_t i = 1; i < axis.m_binEdges.size(); i++) {
0713 os << ", " << axis.m_binEdges[i];
0714 }
0715 os << ")";
0716 return os;
0717 }
0718
0719 protected:
0720 void toStream(std::ostream& os) const override { os << *this; }
0721
0722 private:
0723
0724 std::vector<ActsScalar> m_binEdges;
0725 };
0726 }