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