Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-05 08:09:28

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