Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-07-05 08:11:17

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