File indexing completed on 2025-10-31 08:16:24
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 }