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