Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-06-05 08:29:44

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2017-2018 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
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 // This object can be iterated to produce up to two sequences of integer
0021 // indices, corresponding to the half-open integer ranges [begin1, end1[ and
0022 // [begin2, end2[.
0023 //
0024 // The goal is to emulate the effect of enumerating a range of neighbor
0025 // indices on an axis (which may go out of bounds and wrap around since we
0026 // have AxisBoundaryType::Closed), inserting them into an std::vector, and
0027 // discarding duplicates, without paying the price of duplicate removal
0028 // and dynamic memory allocation in hot magnetic field interpolation code.
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     // Specialized constructor for end() iterator
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   // Number of indices that will be produced if this sequence is iterated
0081   std::size_t size() const { return (m_end1 - m_begin1) + (m_end2 - m_begin2); }
0082 
0083   // Collect the sequence of indices into an std::vector
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 /// @brief calculate bin indices for an equidistant binning
0098 ///
0099 /// This class provides some basic functionality for calculating bin indices
0100 /// for a given equidistant binning.
0101 template <AxisBoundaryType bdt>
0102 class Axis<AxisType::Equidistant, bdt> final : public IAxis {
0103  public:
0104   static constexpr AxisType type = AxisType::Equidistant;
0105 
0106   /// @brief default constructor
0107   ///
0108   /// @param [in] xmin lower boundary of axis range
0109   /// @param [in] xmax upper boundary of axis range
0110   /// @param [in] nBins number of bins to divide the axis range into
0111   ///
0112   /// Divide the range \f$[\text{xmin},\text{xmax})\f$ into \f$\text{nBins}\f$
0113   /// equidistant bins.
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   /// Constructor with a tag for the boundary type
0121   ///
0122   /// @param [in] typeTag boundary type tag
0123   /// @param [in] xmin lower boundary of axis range
0124   /// @param [in] xmax upper boundary of axis range
0125   /// @param [in] nBins number of bins to divide the axis range into
0126   ///
0127   /// Divide the range \f$[\text{xmin},\text{xmax})\f$ into \f$\text{nBins}\f$
0128   /// equidistant bins.
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   /// @brief returns whether the axis is equidistant
0136   ///
0137   /// @return bool is equidistant
0138   bool isEquidistant() const override { return true; }
0139 
0140   /// @brief returns whether the axis is variable
0141   ///
0142   /// @return bool is variable
0143   bool isVariable() const override { return false; }
0144 
0145   /// @brief returns the type of the axis
0146   /// @return @c AxisType of this axis
0147   AxisType getType() const override { return type; }
0148 
0149   /// @brief returns the boundary type set in the template param
0150   ///
0151   /// @return @c AxisBoundaryType of this axis
0152   AxisBoundaryType getBoundaryType() const override { return bdt; }
0153 
0154   /// @brief Get #size bins which neighbor the one given
0155   ///
0156   /// Generic overload with symmetric size
0157   ///
0158   /// @param [in] idx requested bin index
0159   /// @param [in] size how many neighboring bins (up/down)
0160   /// @return Set of neighboring bin indices (global)
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   /// @brief Get #size bins which neighbor the one given
0168   ///
0169   /// This is the version for Open
0170   ///
0171   /// @param [in] idx requested bin index
0172   /// @param [in] sizes how many neighboring bins (up/down)
0173   /// @return Set of neighboring bin indices (global)
0174   /// @note Open varies given bin and allows 0 and NBins+1 (underflow,
0175   /// overflow)
0176   ///       as neighbors
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   /// @brief Get #size bins which neighbor the one given
0191   ///
0192   /// This is the version for Bound
0193   ///
0194   /// @param [in] idx requested bin index
0195   /// @param [in] sizes how many neighboring bins (up/down)
0196   /// @return Set of neighboring bin indices (global)
0197   /// @note Bound varies given bin and allows 1 and NBins (regular bins)
0198   ///       as neighbors
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   /// @brief Get #size bins which neighbor the one given
0216   ///
0217   /// This is the version for Closed (i.e. Wrapping)
0218   ///
0219   /// @param [in] idx requested bin index
0220   /// @param [in] sizes how many neighboring bins (up/down)
0221   /// @return Set of neighboring bin indices (global)
0222   /// @note Closed varies given bin and allows bins on the opposite
0223   ///       side of the axis as neighbors. (excludes underflow / overflow)
0224   NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0225                                           std::pair<int, int> sizes = {-1,
0226                                                                        1}) const
0227     requires(bdt == AxisBoundaryType::Closed)
0228   {
0229     // Handle invalid indices
0230     if (idx <= 0 || idx >= (getNBins() + 1)) {
0231       return NeighborHoodIndices();
0232     }
0233 
0234     // Handle corner case where user requests more neighbours than the number
0235     // of bins on the axis. All bins are returned in this case.
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     // If the entire index range is not covered, we must wrap the range of
0246     // targeted neighbor indices into the range of valid bin indices. This may
0247     // split the range of neighbor indices in two parts:
0248     //
0249     // Before wraparound - [        XXXXX]XXX
0250     // After wraparound  - [ XXXX   XXXX ]
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   /// @brief Converts bin index into a valid one for this axis.
0264   ///
0265   /// @note Open: bin index is clamped to [0, nBins+1]
0266   ///
0267   /// @param [in] bin The bin to wrap
0268   /// @return valid bin index
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   /// @brief Converts bin index into a valid one for this axis.
0276   ///
0277   /// @note Bound: bin index is clamped to [1, nBins]
0278   ///
0279   /// @param [in] bin The bin to wrap
0280   /// @return valid bin index
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   /// @brief Converts bin index into a valid one for this axis.
0288   ///
0289   /// @note Closed: bin index wraps around to other side
0290   ///
0291   /// @param [in] bin The bin to wrap
0292   /// @return valid bin index
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     // return int(bin<1)*w - int(bin>w)*w + bin;
0299   }
0300 
0301   /// @brief get corresponding bin index for given coordinate
0302   ///
0303   /// @param  [in] x input coordinate
0304   /// @return index of bin containing the given value
0305   ///
0306   /// @note Bin intervals are defined with closed lower bounds and open upper
0307   ///       bounds, that is \f$l <= x < u\f$ if the value @c x lies within a
0308   ///       bin with lower bound @c l and upper bound @c u.
0309   /// @note Bin indices start at @c 1. The underflow bin has the index @c 0
0310   ///       while the index <tt>nBins + 1</tt> indicates the overflow bin .
0311   std::size_t getBin(ActsScalar x) const {
0312     return wrapBin(
0313         static_cast<int>(std::floor((x - getMin()) / getBinWidth()) + 1));
0314   }
0315 
0316   /// @brief get bin width
0317   ///
0318   /// @return constant width for all bins
0319   ActsScalar getBinWidth(std::size_t /*bin*/ = 0) const { return m_width; }
0320 
0321   /// @brief get lower bound of bin
0322   ///
0323   /// @param  [in] bin index of bin
0324   /// @return lower bin boundary
0325   ///
0326   /// @pre @c bin must be a valid bin index (excluding the underflow bin),
0327   ///      i.e. \f$1 \le \text{bin} \le \text{nBins} + 1\f$
0328   ///
0329   /// @note Bin intervals have a closed lower bound, i.e. the lower boundary
0330   ///       belongs to the bin with the given bin index.
0331   ActsScalar getBinLowerBound(std::size_t bin) const {
0332     return getMin() + (bin - 1) * getBinWidth();
0333   }
0334 
0335   /// @brief get upper bound of bin
0336   ///
0337   /// @param  [in] bin index of bin
0338   /// @return upper bin boundary
0339   ///
0340   /// @pre @c bin must be a valid bin index (excluding the overflow bin),
0341   ///      i.e. \f$0 \le \text{bin} \le \text{nBins}\f$
0342   ///
0343   /// @note Bin intervals have an open upper bound, i.e. the upper boundary
0344   ///       does @b not belong to the bin with the given bin index.
0345   ActsScalar getBinUpperBound(std::size_t bin) const {
0346     return getMin() + bin * getBinWidth();
0347   }
0348 
0349   /// @brief get bin center
0350   ///
0351   /// @param  [in] bin index of bin
0352   /// @return bin center position
0353   ///
0354   /// @pre @c bin must be a valid bin index (excluding under-/overflow bins),
0355   ///      i.e. \f$1 \le \text{bin} \le \text{nBins}\f$
0356   ActsScalar getBinCenter(std::size_t bin) const {
0357     return getMin() + (bin - 0.5) * getBinWidth();
0358   }
0359 
0360   /// @brief get maximum of binning range
0361   ///
0362   /// @return maximum of binning range
0363   ActsScalar getMax() const override { return m_max; }
0364 
0365   /// @brief get minimum of binning range
0366   ///
0367   /// @return minimum of binning range
0368   ActsScalar getMin() const override { return m_min; }
0369 
0370   /// @brief get total number of bins
0371   ///
0372   /// @return total number of bins (excluding under-/overflow bins)
0373   std::size_t getNBins() const override { return m_bins; }
0374 
0375   /// @brief check whether value is inside axis limits
0376   ///
0377   /// @return @c true if \f$\text{xmin} \le x < \text{xmax}\f$, otherwise
0378   ///         @c false
0379   ///
0380   /// @post If @c true is returned, the bin containing the given value is a
0381   ///       valid bin, i.e. it is neither the underflow nor the overflow bin.
0382   bool isInside(ActsScalar x) const { return (m_min <= x) && (x < m_max); }
0383 
0384   /// @brief Return a vector of bin edges
0385   /// @return Vector which contains the bin edges
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   /// minimum of binning range
0408   ActsScalar m_min{};
0409   /// maximum of binning range
0410   ActsScalar m_max{};
0411   /// constant bin width
0412   ActsScalar m_width{};
0413   /// number of bins (excluding under-/overflow bins)
0414   std::size_t m_bins{};
0415 };
0416 
0417 /// @brief calculate bin indices for a variable binning
0418 ///
0419 /// This class provides some basic functionality for calculating bin indices
0420 /// for a given binning with variable bin sizes.
0421 template <AxisBoundaryType bdt>
0422 class Axis<AxisType::Variable, bdt> final : public IAxis {
0423  public:
0424   static constexpr AxisType type = AxisType::Variable;
0425 
0426   /// @param [in] binEdges vector of bin edges
0427   /// @pre @c binEdges must be strictly sorted in ascending order.
0428   /// @pre @c binEdges must contain at least two entries.
0429   ///
0430   /// Create a binning structure with @c nBins variable-sized bins from the
0431   /// given bin boundaries. @c nBins is given by the number of bin edges
0432   /// reduced by one.
0433   explicit Axis(std::vector<ActsScalar> binEdges)
0434       : m_binEdges(std::move(binEdges)) {}
0435 
0436   /// @param [in] typeTag boundary type tag
0437   /// @param [in] binEdges vector of bin edges
0438   /// @pre @c binEdges must be strictly sorted in ascending order.
0439   /// @pre @c binEdges must contain at least two entries.
0440   ///
0441   /// Create a binning structure with @c nBins variable-sized bins from the
0442   /// given bin boundaries. @c nBins is given by the number of bin edges
0443   /// reduced by one.
0444   Axis(AxisBoundaryTypeTag<bdt> typeTag, std::vector<ActsScalar> binEdges)
0445       : Axis(std::move(binEdges)) {
0446     static_cast<void>(typeTag);
0447   }
0448 
0449   /// @brief returns whether the axis is equidistante
0450   ///
0451   /// @return bool is equidistant
0452   bool isEquidistant() const override { return false; }
0453 
0454   /// @brief returns whether the axis is variable
0455   ///
0456   /// @return bool is variable
0457   bool isVariable() const override { return true; }
0458 
0459   /// @brief returns the type of the axis
0460   /// @return @c AxisType of this axis
0461   AxisType getType() const override { return type; }
0462 
0463   /// @brief returns the boundary type set in the template param
0464   ///
0465   /// @return @c AxisBoundaryType of this axis
0466   AxisBoundaryType getBoundaryType() const override { return bdt; }
0467 
0468   /// @brief Get #size bins which neighbor the one given
0469   ///
0470   /// Generic overload with symmetric size
0471   ///
0472   /// @param [in] idx requested bin index
0473   /// @param [in] size how many neighboring bins
0474   /// @return Set of neighboring bin indices (global)
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   /// @brief Get #size bins which neighbor the one given
0482   ///
0483   /// This is the version for Open
0484   ///
0485   /// @param [in] idx requested bin index
0486   /// @param [in] sizes how many neighboring bins (up/down)
0487   /// @return Set of neighboring bin indices (global)
0488   /// @note Open varies given bin and allows 0 and NBins+1 (underflow,
0489   /// overflow)
0490   ///       as neighbors
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   /// @brief Get #size bins which neighbor the one given
0504   ///
0505   /// This is the version for Bound
0506   ///
0507   /// @param [in] idx requested bin index
0508   /// @param [in] sizes how many neighboring bins (up/down)
0509   /// @return Set of neighboring bin indices (global)
0510   /// @note Bound varies given bin and allows 1 and NBins (regular bins)
0511   ///       as neighbors
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   /// @brief Get #size bins which neighbor the one given
0528   ///
0529   /// This is the version for Closed
0530   ///
0531   /// @param [in] idx requested bin index
0532   /// @param [in] sizes how many neighboring bins (up/down)
0533   /// @return Set of neighboring bin indices (global)
0534   /// @note Closed varies given bin and allows bins on the opposite
0535   ///       side of the axis as neighbors. (excludes underflow / overflow)
0536   NeighborHoodIndices neighborHoodIndices(std::size_t idx,
0537                                           std::pair<int, int> sizes = {-1,
0538                                                                        1}) const
0539     requires(bdt == AxisBoundaryType::Closed)
0540   {
0541     // Handle invalid indices
0542     if (idx <= 0 || idx >= (getNBins() + 1)) {
0543       return NeighborHoodIndices();
0544     }
0545 
0546     // Handle corner case where user requests more neighbours than the number
0547     // of bins on the axis. All bins are returned in this case
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     // If the entire index range is not covered, we must wrap the range of
0558     // targeted neighbor indices into the range of valid bin indices. This may
0559     // split the range of neighbor indices in two parts:
0560     //
0561     // Before wraparound - [        XXXXX]XXX
0562     // After wraparound  - [ XXXX   XXXX ]
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   /// @brief Converts bin index into a valid one for this axis.
0576   ///
0577   /// @note Open: bin index is clamped to [0, nBins+1]
0578   ///
0579   /// @param [in] bin The bin to wrap
0580   /// @return valid bin index
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   /// @brief Converts bin index into a valid one for this axis.
0588   ///
0589   /// @note Bound: bin index is clamped to [1, nBins]
0590   ///
0591   /// @param [in] bin The bin to wrap
0592   /// @return valid bin index
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   /// @brief Converts bin index into a valid one for this axis.
0600   ///
0601   /// @note Closed: bin index wraps around to other side
0602   ///
0603   /// @param [in] bin The bin to wrap
0604   /// @return valid bin index
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     // return int(bin<1)*w - int(bin>w)*w + bin;
0611   }
0612 
0613   /// @brief get corresponding bin index for given coordinate
0614   ///
0615   /// @param  [in] x input coordinate
0616   /// @return index of bin containing the given value
0617   ///
0618   /// @note Bin intervals are defined with closed lower bounds and open upper
0619   ///       bounds, that is \f$l <= x < u\f$ if the value @c x lies within a
0620   ///       bin with lower bound @c l and upper bound @c u.
0621   /// @note Bin indices start at @c 1. The underflow bin has the index @c 0
0622   ///       while the index <tt>nBins + 1</tt> indicates the overflow bin .
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   /// @brief get bin width
0630   ///
0631   /// @param  [in] bin index of bin
0632   /// @return width of given bin
0633   ///
0634   /// @pre @c bin must be a valid bin index (excluding under-/overflow bins),
0635   ///      i.e. \f$1 \le \text{bin} \le \text{nBins}\f$
0636   ActsScalar getBinWidth(std::size_t bin) const {
0637     return m_binEdges.at(bin) - m_binEdges.at(bin - 1);
0638   }
0639 
0640   /// @brief get lower bound of bin
0641   ///
0642   /// @param  [in] bin index of bin
0643   /// @return lower bin boundary
0644   ///
0645   /// @pre @c bin must be a valid bin index (excluding the underflow bin),
0646   ///      i.e. \f$1 \le \text{bin} \le \text{nBins} + 1\f$
0647   ///
0648   /// @note Bin intervals have a closed lower bound, i.e. the lower boundary
0649   ///       belongs to the bin with the given bin index.
0650   ActsScalar getBinLowerBound(std::size_t bin) const {
0651     return m_binEdges.at(bin - 1);
0652   }
0653 
0654   /// @brief get upper bound of bin
0655   ///
0656   /// @param  [in] bin index of bin
0657   /// @return upper bin boundary
0658   ///
0659   /// @pre @c bin must be a valid bin index (excluding the overflow bin),
0660   ///      i.e. \f$0 \le \text{bin} \le \text{nBins}\f$
0661   ///
0662   /// @note Bin intervals have an open upper bound, i.e. the upper boundary
0663   ///       does @b not belong to the bin with the given bin index.
0664   ActsScalar getBinUpperBound(std::size_t bin) const {
0665     return m_binEdges.at(bin);
0666   }
0667 
0668   /// @brief get bin center
0669   ///
0670   /// @param  [in] bin index of bin
0671   /// @return bin center position
0672   ///
0673   /// @pre @c bin must be a valid bin index (excluding under-/overflow bins),
0674   ///      i.e. \f$1 \le \text{bin} \le \text{nBins}\f$
0675   ActsScalar getBinCenter(std::size_t bin) const {
0676     return 0.5 * (getBinLowerBound(bin) + getBinUpperBound(bin));
0677   }
0678 
0679   /// @brief get maximum of binning range
0680   ///
0681   /// @return maximum of binning range
0682   ActsScalar getMax() const override { return m_binEdges.back(); }
0683 
0684   /// @brief get minimum of binning range
0685   ///
0686   /// @return minimum of binning range
0687   ActsScalar getMin() const override { return m_binEdges.front(); }
0688 
0689   /// @brief get total number of bins
0690   ///
0691   /// @return total number of bins (excluding under-/overflow bins)
0692   std::size_t getNBins() const override { return m_binEdges.size() - 1; }
0693 
0694   /// @brief check whether value is inside axis limits
0695   ///
0696   /// @return @c true if \f$\text{xmin} \le x < \text{xmax}\f$, otherwise
0697   ///         @c false
0698   ///
0699   /// @post If @c true is returned, the bin containing the given value is a
0700   ///       valid bin, i.e. it is neither the underflow nor the overflow bin.
0701   bool isInside(ActsScalar x) const {
0702     return (m_binEdges.front() <= x) && (x < m_binEdges.back());
0703   }
0704 
0705   /// @brief Return a vector of bin edges
0706   /// @return Vector which contains the bin edges
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   /// vector of bin edges (sorted in ascending order)
0724   std::vector<ActsScalar> m_binEdges;
0725 };
0726 }  // namespace Acts