Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-10-19 07:57:28

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/Definitions/Algebra.hpp"
0012 #include "Acts/Utilities/AxisDefinitions.hpp"
0013 #include "Acts/Utilities/BinningType.hpp"
0014 #include "Acts/Utilities/Helpers.hpp"
0015 #include "Acts/Utilities/ProtoAxis.hpp"
0016 #include "Acts/Utilities/ThrowAssert.hpp"
0017 #include "Acts/Utilities/VectorHelpers.hpp"
0018 
0019 #include <algorithm>
0020 #include <cmath>
0021 #include <memory>
0022 #include <sstream>
0023 #include <string>
0024 #include <utility>
0025 #include <vector>
0026 
0027 namespace Acts {
0028 
0029 /// @class BinningData
0030 ///
0031 ///   This class holds all the data necessary for the bin calculation
0032 ///
0033 ///   phi has a very particular behaviour:
0034 ///   - there's the change around +/- PI
0035 ///
0036 ///   - it can be multiplicative or additive
0037 ///   multiplicative : each major bin has the same sub structure
0038 ///                    i.e. first binnning
0039 ///
0040 /// structure is equidistant
0041 ///   additive : sub structure replaces one bin (and one bin only)
0042 ///
0043 ///
0044 class BinningData {
0045  public:
0046   BinningType type{};        ///< binning type: equidistant, arbitrary
0047   BinningOption option{};    ///< binning option: open, closed
0048   AxisDirection binvalue{};  ///< axis direction: AxisX, AxisY, AxisZ, ...
0049   float min{};               ///< minimum value
0050   float max{};               ///< maximum value
0051   float step{};              ///< binning step
0052   bool zdim{};               ///< zero dimensional binning : direct access
0053 
0054   /// sub structure: describe some sub binning
0055   std::unique_ptr<const BinningData> subBinningData;
0056   /// sub structure: additive or multipicative
0057   bool subBinningAdditive{};
0058 
0059   /// Constructor for 0D binning
0060   ///
0061   /// @param bValue is the axis direction AxisX, AxisY, etc.
0062   /// @param bMin is the minimum value
0063   /// @param bMax is the maximum value
0064   BinningData(AxisDirection bValue, float bMin, float bMax)
0065       : type(equidistant),
0066         option(open),
0067         binvalue(bValue),
0068         min(bMin),
0069         max(bMax),
0070         step((bMax - bMin)),
0071         zdim(true),
0072         subBinningData(nullptr),
0073         m_bins(1),
0074         m_boundaries({{min, max}}),
0075         m_totalBins(1),
0076         m_totalBoundaries(std::vector<float>()),
0077         m_functionPtr(&searchEquidistantWithBoundary) {}
0078 
0079   /// Constructor for equidistant binning
0080   /// and optional sub structure can be
0081   /// multiplicative or additive
0082   ///
0083   /// @param bOption is the binning option : open, closed
0084   /// @param bValue is the axis direction: Axis, AxisY, etc.
0085   /// @param bBins is number of equidistant bins
0086   /// @param bMin is the minimum value
0087   /// @param bMax is the maximum value
0088   /// @param sBinData is (optional) sub structure
0089   /// @param sBinAdditive is the prescription for the sub structure
0090   BinningData(BinningOption bOption, AxisDirection bValue, std::size_t bBins,
0091               float bMin, float bMax,
0092               std::unique_ptr<const BinningData> sBinData = nullptr,
0093               bool sBinAdditive = false)
0094       : type(equidistant),
0095         option(bOption),
0096         binvalue(bValue),
0097         min(bMin),
0098         max(bMax),
0099         step((bMax - bMin) / bBins),
0100         zdim(bBins == 1 ? true : false),
0101         subBinningData(std::move(sBinData)),
0102         subBinningAdditive(sBinAdditive),
0103         m_bins(bBins),
0104         m_boundaries(std::vector<float>()),
0105         m_totalBins(bBins),
0106         m_totalBoundaries(std::vector<float>()) {
0107     // set to equidistant search
0108     m_functionPtr = &searchEquidistantWithBoundary;
0109     // fill the boundary vector for fast access to center & boundaries
0110     m_boundaries.reserve(m_bins + 1);
0111     for (std::size_t ib = 0; ib < m_bins + 1; ++ib) {
0112       m_boundaries.push_back(min + ib * step);
0113     }
0114     // the binning data has sub structure - multiplicative or additive
0115     checkSubStructure();
0116   }
0117 
0118   /// Constructor for non-equidistant binning
0119   ///
0120   /// @param bOption is the binning option : open / closed
0121   /// @param bValue is the axis direction : AxisX, AxisY, etc.
0122   /// @param bBoundaries are the bin boundaries
0123   /// @param sBinData is (optional) sub structure
0124   BinningData(BinningOption bOption, AxisDirection bValue,
0125               const std::vector<float>& bBoundaries,
0126               std::unique_ptr<const BinningData> sBinData = nullptr)
0127       : type(arbitrary),
0128         option(bOption),
0129         binvalue(bValue),
0130         zdim(bBoundaries.size() == 2 ? true : false),
0131         subBinningData(std::move(sBinData)),
0132         subBinningAdditive(true),
0133         m_bins(bBoundaries.size() - 1),
0134         m_boundaries(bBoundaries),
0135         m_totalBins(bBoundaries.size() - 1),
0136         m_totalBoundaries(bBoundaries) {
0137     // assert a no-size case
0138     throw_assert(m_boundaries.size() > 1, "Must have more than one boundary");
0139     min = m_boundaries[0];
0140     max = m_boundaries[m_boundaries.size() - 1];
0141     // set to equidistant search
0142     m_functionPtr = &searchInVectorWithBoundary;
0143     // the binning data has sub structure - multiplicative
0144     checkSubStructure();
0145   }
0146 
0147   /// Copy constructor
0148   ///
0149   /// @param bdata is the source object
0150   BinningData(const BinningData& bdata)
0151       : type(bdata.type),
0152         option(bdata.option),
0153         binvalue(bdata.binvalue),
0154         min(bdata.min),
0155         max(bdata.max),
0156         step(bdata.step),
0157         zdim(bdata.zdim),
0158         subBinningData(nullptr),
0159         subBinningAdditive(bdata.subBinningAdditive),
0160         m_bins(bdata.m_bins),
0161         m_boundaries(bdata.m_boundaries),
0162         m_totalBins(bdata.m_totalBins),
0163         m_totalBoundaries(bdata.m_totalBoundaries) {
0164     // get the binning data
0165     subBinningData =
0166         bdata.subBinningData
0167             ? std::make_unique<const BinningData>(*bdata.subBinningData)
0168             : nullptr;
0169     // set the pointer depending on the type
0170     // set the correct function pointer
0171     if (type == equidistant) {
0172       m_functionPtr = &searchEquidistantWithBoundary;
0173     } else {
0174       m_functionPtr = &searchInVectorWithBoundary;
0175     }
0176   }
0177 
0178   /// Constructor from DirectedProtoAxis
0179   ///
0180   /// @param dpAxis is the ProtoAxis object
0181   ///
0182   explicit BinningData(const DirectedProtoAxis& dpAxis)
0183       : binvalue(dpAxis.getAxisDirection()), subBinningData(nullptr) {
0184     const auto& axis = dpAxis.getAxis();
0185     type = axis.getType() == AxisType::Equidistant ? equidistant : arbitrary;
0186     option = axis.getBoundaryType() == AxisBoundaryType::Closed ? closed : open;
0187     min = static_cast<float>(axis.getMin());
0188     max = static_cast<float>(axis.getMax());
0189     m_bins = axis.getNBins();
0190     step = (max - min) / static_cast<float>(m_bins);
0191     zdim = (m_bins == 1);
0192     m_boundaries.reserve(axis.getBinEdges().size());
0193     for (const auto& edge : axis.getBinEdges()) {
0194       m_boundaries.push_back(static_cast<float>(edge));
0195     }
0196     m_totalBins = m_bins;
0197     m_totalBoundaries = m_boundaries;
0198     // Set the search function pointer based on axis type
0199     m_functionPtr = (type == equidistant) ? &searchEquidistantWithBoundary
0200                                           : &searchInVectorWithBoundary;
0201   }
0202 
0203   /// Assignment operator
0204   ///
0205   /// @param bdata is the source object
0206   /// @return Reference to this BinningData after assignment
0207   BinningData& operator=(const BinningData& bdata) {
0208     if (this != &bdata) {
0209       type = bdata.type;
0210       option = bdata.option;
0211       binvalue = bdata.binvalue;
0212       min = bdata.min;
0213       max = bdata.max;
0214       step = bdata.step;
0215       zdim = bdata.zdim;
0216       subBinningAdditive = bdata.subBinningAdditive;
0217       subBinningData =
0218           bdata.subBinningData
0219               ? std::make_unique<const BinningData>(*bdata.subBinningData)
0220               : nullptr;
0221       m_bins = bdata.m_bins;
0222       m_boundaries = bdata.m_boundaries;
0223       m_totalBins = bdata.m_totalBins;
0224       m_totalBoundaries = bdata.m_totalBoundaries;
0225       // set the correct function pointer
0226       if (type == equidistant) {
0227         m_functionPtr = &searchEquidistantWithBoundary;
0228       } else {
0229         m_functionPtr = &searchInVectorWithBoundary;
0230       }
0231     }
0232     return (*this);
0233   }
0234 
0235   BinningData() = default;
0236   ~BinningData() = default;
0237 
0238   /// Equality operator
0239   ///
0240   /// @param bData is the binning data to be checked against
0241   ///
0242   /// @return a boolean indicating if they are the same
0243   bool operator==(const BinningData& bData) const {
0244     return (type == bData.type && option == bData.option &&
0245             binvalue == bData.binvalue && min == bData.min &&
0246             max == bData.max && step == bData.step && zdim == bData.zdim &&
0247             ((subBinningData == nullptr && bData.subBinningData == nullptr) ||
0248              (subBinningData != nullptr && bData.subBinningData != nullptr &&
0249               (*subBinningData == *bData.subBinningData))) &&
0250             subBinningAdditive == bData.subBinningAdditive);
0251   }
0252 
0253   /// Return the number of bins - including sub bins
0254   /// @return Total number of bins including sub-bins
0255   std::size_t bins() const { return m_totalBins; }
0256 
0257   /// Return the boundaries  - including sub boundaries
0258   /// @return vector of floats indicating the boundary values
0259   const std::vector<float>& boundaries() const {
0260     if (subBinningData) {
0261       return m_totalBoundaries;
0262     }
0263     return m_boundaries;
0264   }
0265 
0266   /// Take the right float value
0267   ///
0268   /// @param lposition assumes the correct local position expression
0269   ///
0270   /// @return float value according to the binning setup
0271   float value(const Vector2& lposition) const {
0272     // ordered after occurrence
0273     if (binvalue == AxisDirection::AxisR ||
0274         binvalue == AxisDirection::AxisRPhi ||
0275         binvalue == AxisDirection::AxisX ||
0276         binvalue == AxisDirection::AxisTheta) {
0277       return lposition[0];
0278     }
0279 
0280     return lposition[1];
0281   }
0282 
0283   /// Take the right float value
0284   ///
0285   /// @param position is the global position
0286   ///
0287   /// @return float value according to the binning setup
0288   float value(const Vector3& position) const {
0289     using VectorHelpers::eta;
0290     using VectorHelpers::perp;
0291     using VectorHelpers::phi;
0292     // ordered after occurrence
0293     if (binvalue == AxisDirection::AxisR ||
0294         binvalue == AxisDirection::AxisTheta) {
0295       return (perp(position));
0296     }
0297     if (binvalue == AxisDirection::AxisRPhi) {
0298       return (perp(position) * phi(position));
0299     }
0300     if (binvalue == AxisDirection::AxisEta) {
0301       return (eta(position));
0302     }
0303     if (toUnderlying(binvalue) < 3) {
0304       return static_cast<float>(position[toUnderlying(binvalue)]);
0305     }
0306     // phi gauging
0307     return phi(position);
0308   }
0309 
0310   /// Get the center value of a bin
0311   ///
0312   /// @param bin is the bin for which the center value is requested
0313   ///
0314   /// @return float value according to the bin center
0315   float center(std::size_t bin) const {
0316     const std::vector<float>& bvals = boundaries();
0317     // take the center between bin boundaries
0318     float value =
0319         bin < (bvals.size() - 1) ? 0.5 * (bvals[bin] + bvals[bin + 1]) : 0.;
0320     return value;
0321   }
0322 
0323   /// Get the width of a bin
0324   ///
0325   /// @param bin is the bin for which the width is requested
0326   ///
0327   /// @return float value of width
0328   float width(std::size_t bin) const {
0329     const std::vector<float>& bvals = boundaries();
0330     // take the center between bin boundaries
0331     float value = bin < (bvals.size() - 1) ? bvals[bin + 1] - bvals[bin] : 0.;
0332     return value;
0333   }
0334 
0335   /// Check if bin is inside from Vector3
0336   ///
0337   /// @param position is the search position in global coordinated
0338   ///
0339   /// @return boolean if this is inside() method is true
0340   bool inside(const Vector3& position) const {
0341     // closed one is always inside
0342     if (option == closed) {
0343       return true;
0344     }
0345     // all other options
0346     // @todo remove hard-coded tolerance parameters
0347     float val = value(position);
0348     return (val > min - 0.001 && val < max + 0.001);
0349   }
0350 
0351   /// Check if bin is inside from Vector2
0352   ///
0353   /// @param lposition is the search position in global coordinated
0354   ///
0355   /// @return boolean if this is inside() method is true
0356   bool inside(const Vector2& lposition) const {
0357     // closed one is always inside
0358     if (option == closed) {
0359       return true;
0360     }
0361     // all other options
0362     // @todo remove hard-coded tolerance parameters
0363     float val = value(lposition);
0364     return (val > min - 0.001 && val < max + 0.001);
0365   }
0366 
0367   /// Generic search from a 2D position
0368   /// -- corresponds to local coordinate schema
0369   /// @param lposition is the search position in local coordinated
0370   ///
0371   /// @return bin according tot this
0372   std::size_t searchLocal(const Vector2& lposition) const {
0373     if (zdim) {
0374       return 0;
0375     }
0376     return search(value(lposition));
0377   }
0378 
0379   /// Generic search from a 3D position
0380   /// -- corresponds to global coordinate schema
0381   /// @param position is the search position in global coordinated
0382   ///
0383   /// @return bin according tot this
0384   std::size_t searchGlobal(const Vector3& position) const {
0385     if (zdim) {
0386       return 0;
0387     }
0388     return search(value(position));
0389   }
0390 
0391   /// Generic search - forwards to correct function pointer
0392   ///
0393   /// @param value is the searchvalue as float
0394   ///
0395   /// @return bin according tot this
0396   std::size_t search(float value) const {
0397     if (zdim) {
0398       return 0;
0399     }
0400     assert(m_functionPtr != nullptr);
0401     return (!subBinningData) ? (*m_functionPtr)(value, *this)
0402                              : searchWithSubStructure(value);
0403   }
0404 
0405   ///  Generic search with sub structure
0406   /// - forwards to correct function pointer
0407   ///
0408   /// @param value is the searchvalue as float
0409   ///
0410   /// @return bin according tot this
0411   std::size_t searchWithSubStructure(float value) const {
0412     // find the masterbin with the correct function pointer
0413     std::size_t masterbin = (*m_functionPtr)(value, *this);
0414     // additive sub binning -
0415     if (subBinningAdditive) {
0416       // no gauging done, for additive sub structure
0417       return masterbin + subBinningData->search(value);
0418     }
0419     // gauge the value to the subBinData
0420     float gvalue =
0421         value - masterbin * (subBinningData->max - subBinningData->min);
0422     // now go / additive or multiplicative
0423     std::size_t subbin = subBinningData->search(gvalue);
0424     // now return
0425     return masterbin * subBinningData->bins() + subbin;
0426   }
0427 
0428   /// Layer next direction is needed
0429   ///
0430   /// @param position is the start search position
0431   /// @param dir is the direction
0432   /// @todo check if this can be changed
0433   ///
0434   /// @return integer that indicates which direction to move
0435   int nextDirection(const Vector3& position, const Vector3& dir) const {
0436     if (zdim) {
0437       return 0;
0438     }
0439     float val = value(position);
0440     Vector3 probe = position + dir.normalized();
0441     float nextval = value(probe);
0442     return (nextval > val) ? 1 : -1;
0443   }
0444 
0445   /// access to the center value
0446   /// this uses the bin boundary vector, it also works with sub structure
0447   ///
0448   /// @param bin is the bin for which the value is requested, if bin > nbins
0449   /// it is set to max
0450   ///
0451   /// @return the center value of the bin is given
0452   float centerValue(std::size_t bin) const {
0453     if (zdim) {
0454       return 0.5 * (min + max);
0455     }
0456     float bmin = m_boundaries[bin];
0457     float bmax = bin < m_boundaries.size() ? m_boundaries[bin + 1] : max;
0458     return 0.5 * (bmin + bmax);
0459   }
0460 
0461  private:
0462   std::size_t m_bins{};             ///< number of bins
0463   std::vector<float> m_boundaries;  ///< vector of holding the bin boundaries
0464   std::size_t m_totalBins{};        ///< including potential substructure
0465   std::vector<float> m_totalBoundaries;  ///< including potential substructure
0466 
0467   std::size_t (*m_functionPtr)(float,
0468                                const BinningData&){};  /// function pointer
0469 
0470   /// helper method to set the sub structure
0471   void checkSubStructure() {
0472     // sub structure is only checked when sBinData is defined
0473     if (subBinningData) {
0474       m_totalBoundaries.clear();
0475       // (A) additive sub structure
0476       if (subBinningAdditive) {
0477         // one bin is replaced by the sub bins
0478         m_totalBins = m_bins + subBinningData->bins() - 1;
0479         // the tricky one - exchange one bin by many others
0480         m_totalBoundaries.reserve(m_totalBins + 1);
0481         // get the sub bin boundaries
0482         const std::vector<float>& subBinBoundaries =
0483             subBinningData->boundaries();
0484         float sBinMin = subBinBoundaries[0];
0485         // get the min value of the sub bin boundaries
0486         std::vector<float>::const_iterator mbvalue = m_boundaries.begin();
0487         for (; mbvalue != m_boundaries.end(); ++mbvalue) {
0488           // should define numerically stable
0489           if (std::abs((*mbvalue) - sBinMin) < 10e-10) {
0490             // copy the sub bin boundaries into the vector
0491             m_totalBoundaries.insert(m_totalBoundaries.begin(),
0492                                      subBinBoundaries.begin(),
0493                                      subBinBoundaries.end());
0494             ++mbvalue;
0495           } else {
0496             m_totalBoundaries.push_back(*mbvalue);
0497           }
0498         }
0499       } else {  // (B) multiplicative sub structure
0500         // every bin is just replaced by the sub binning structure
0501         m_totalBins = m_bins * subBinningData->bins();
0502         m_totalBoundaries.reserve(m_totalBins + 1);
0503         // get the sub bin boundaries if there are any
0504         const std::vector<float>& subBinBoundaries =
0505             subBinningData->boundaries();
0506         // create the boundary vector
0507         m_totalBoundaries.push_back(min);
0508         for (std::size_t ib = 0; ib < m_bins; ++ib) {
0509           float offset = ib * step;
0510           for (std::size_t isb = 1; isb < subBinBoundaries.size(); ++isb) {
0511             m_totalBoundaries.push_back(offset + subBinBoundaries[isb]);
0512           }
0513         }
0514       }
0515       // sort the total boundary vector
0516       std::ranges::sort(m_totalBoundaries);
0517     }
0518   }
0519 
0520   // Equidistant search
0521   // - fastest method
0522   static std::size_t searchEquidistantWithBoundary(float value,
0523                                                    const BinningData& bData) {
0524     // vanilla
0525 
0526     int bin = static_cast<int>((value - bData.min) / bData.step);
0527     // special treatment of the 0 bin for closed
0528     if (bData.option == closed) {
0529       if (value < bData.min) {
0530         return (bData.m_bins - 1);
0531       }
0532       if (value > bData.max) {
0533         return 0;
0534       }
0535     }
0536     // if outside boundary : return boundary for open, opposite bin for closed
0537     bin = bin < 0 ? ((bData.option == open) ? 0 : (bData.m_bins - 1)) : bin;
0538     return static_cast<std::size_t>(
0539         (bin <= static_cast<int>(bData.m_bins - 1))
0540             ? bin
0541             : ((bData.option == open) ? (bData.m_bins - 1) : 0));
0542   }
0543 
0544   // Search in arbitrary boundary
0545   static std::size_t searchInVectorWithBoundary(float value,
0546                                                 const BinningData& bData) {
0547     // lower boundary
0548     if (value <= bData.m_boundaries[0]) {
0549       return (bData.option == closed) ? (bData.m_bins - 1) : 0;
0550     }
0551     // higher boundary
0552     if (value >= bData.max) {
0553       return (bData.option == closed) ? 0 : (bData.m_bins - 1);
0554     }
0555 
0556     auto lb = std::lower_bound(bData.m_boundaries.begin(),
0557                                bData.m_boundaries.end(), value);
0558     return static_cast<std::size_t>(
0559         std::distance(bData.m_boundaries.begin(), lb) - 1);
0560   }
0561 
0562  public:
0563   /// String screen output method
0564   /// @param indent the current indentation
0565   /// @return a string containing the screen information
0566   std::string toString(const std::string& indent = "") const {
0567     std::stringstream sl;
0568     sl << indent << "BinningData object:" << '\n';
0569     sl << indent << "  - type       : " << static_cast<std::size_t>(type)
0570        << '\n';
0571     sl << indent << "  - option     : " << static_cast<std::size_t>(option)
0572        << '\n';
0573     sl << indent << "  - value      : " << static_cast<std::size_t>(binvalue)
0574        << '\n';
0575     sl << indent << "  - bins       : " << bins() << '\n';
0576     sl << indent << "  - min/max    : " << min << " / " << max << '\n';
0577     if (type == equidistant) {
0578       sl << indent << "  - step       : " << step << '\n';
0579     }
0580     sl << indent << "  - boundaries : | ";
0581     for (const auto& b : boundaries()) {
0582       sl << b << " | ";
0583     }
0584     sl << '\n';
0585     return sl.str();
0586   }
0587 };
0588 }  // namespace Acts