Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:11:08

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