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