Back to home page

EIC code displayed by LXR

 
 

    


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

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