Acts/Utilities/BinUtility.hpp

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/BinningData.hpp"
0014 #include "Acts/Utilities/BinningType.hpp"
0015 #include "Acts/Utilities/Enumerate.hpp"
0016 #include "Acts/Utilities/ProtoAxis.hpp"
0017
0018 #include <algorithm>
0019 #include <array>
0020 #include <cstddef>
0021 #include <iostream>
0022 #include <stdexcept>
0023 #include <string>
0024 #include <vector>
0025
0026 namespace Acts {
0027
0028 /// @class BinUtility
0029 ///
0030 /// The BinUtility class that translated global and local position into a bins
0031 /// of a BinnedArray, most performant is equidistant binning without a
0032 /// transform,
0033 /// however, optionally a transform can be provided, e.g. for binning on shifted
0034 /// object, the transform is usually shared with the geometric object the Array
0035 /// is
0036 /// defined on, for performance reasons, also the inverse transform is stored.
0037 ///
0038 class BinUtility {
0039  public:
0040   /// Constructor for equidistant
0041   BinUtility()
0042       : m_binningData(),
0043         m_transform(Transform3::Identity()),
0044         m_itransform(Transform3::Identity()) {
0045     m_binningData.reserve(3);
0046   }
0047
0048   /// Constructor with only a Transform3
0049   ///
0050   /// @param tForm is the local to global transform
0051   explicit BinUtility(const Transform3& tForm)
0052       : m_binningData(), m_transform(tForm), m_itransform(tForm.inverse()) {
0053     m_binningData.reserve(3);
0054   }
0055
0056   /// Constructor from BinningData directly
0057   ///
0058   /// @param bData is the provided binning data
0059   /// @param tForm is the (optional) transform
0060   explicit BinUtility(const BinningData& bData,
0061                       const Transform3& tForm = Transform3::Identity())
0062       : m_binningData(), m_transform(tForm), m_itransform(tForm.inverse()) {
0063     m_binningData.reserve(3);
0064     m_binningData.emplace_back(bData);
0065   }
0066
0067   /// Constructor for equidistant
0068   ///
0069   /// @param bins is the number of bins
0070   /// @param min in the minimal value
0071   /// @param max is the maximal value
0072   /// @param opt is the binning option : open, closed
0073   /// @param value is the axis direction : AxisX, AxisY, AxisZ, etc.
0074   /// @param tForm is the (optional) transform
0075   BinUtility(std::size_t bins, float min, float max, BinningOption opt = open,
0076              AxisDirection value = AxisDirection::AxisX,
0077              const Transform3& tForm = Transform3::Identity())
0078       : m_binningData(), m_transform(tForm), m_itransform(tForm.inverse()) {
0079     m_binningData.reserve(3);
0080     m_binningData.emplace_back(opt, value, bins, min, max);
0081   }
0082
0083   /// Constructor for arbitrary
0084   ///
0085   /// @param bValues is the boundary values of the binning
0086   /// @param opt is the binning option : open, closed
0087   /// @param value is the axis direction : AxisX, AxisY, AxisZ, etc.
0088   /// @param tForm is the (optional) transform
0089   explicit BinUtility(std::vector<float>& bValues, BinningOption opt = open,
0090                       AxisDirection value = AxisDirection::AxisPhi,
0091                       const Transform3& tForm = Transform3::Identity())
0092       : m_binningData(), m_transform(tForm), m_itransform(tForm.inverse()) {
0093     m_binningData.reserve(3);
0094     m_binningData.emplace_back(opt, value, bValues);
0095   }
0096
0097   /// Create from a DirectedProtoAxis
0098   ///
0099   /// @param dpAxis the DirectedProtoAxis to be used
0100   explicit BinUtility(const DirectedProtoAxis& dpAxis)
0101       : m_binningData(),
0102         m_transform(Transform3::Identity()),
0103         m_itransform(Transform3::Identity()) {
0104     m_binningData.reserve(3);
0105     m_binningData.emplace_back(dpAxis);
0106   }
0107
0108   /// Create from several DirectedProtoAxis objects
0109   ///
0110   /// @param dpAxes the DirectedProtoAxis to be used with axis directions
0111   explicit BinUtility(const std::vector<DirectedProtoAxis>& dpAxes)
0112       : m_binningData(),
0113         m_transform(Transform3::Identity()),
0114         m_itransform(Transform3::Identity()) {
0115     m_binningData.reserve(3);
0116     for (const auto& dpAxis : dpAxes) {
0117       m_binningData.emplace_back(dpAxis);
0118     }
0119   }
0120
0121   /// Operator+= to make multidimensional BinUtility
0122   ///
0123   /// @param gbu is the additional BinUtility to be chosen
0124   /// @return Reference to this BinUtility after addition
0125   BinUtility& operator+=(const BinUtility& gbu) {
0126     const std::vector<BinningData>& bData = gbu.binningData();
0127
0128     m_transform = m_transform * gbu.transform();
0129     m_itransform = m_transform.inverse();
0130     if (m_binningData.size() + bData.size() > 3) {
0131       throw std::runtime_error{"BinUtility does not support dim > 3"};
0132     }
0133     m_binningData.insert(m_binningData.end(), bData.begin(), bData.end());
0134     return (*this);
0135   }
0136
0137   /// Equality operator
0138   /// @param other The other BinUtility to compare with
0139   /// @return True if the BinUtilities are equal, false otherwise
0140   bool operator==(const BinUtility& other) const {
0141     return (m_transform.isApprox(other.m_transform) &&
0142             m_binningData == other.binningData());
0143   }
0144
0145   /// Return the binning data vector
0146   /// @return Reference to the vector of binning data
0147   const std::vector<BinningData>& binningData() const { return m_binningData; }
0148
0149   /// Return the total number of bins
0150   /// @return Total number of bins across all dimensions
0151   std::size_t bins() const { return bins(0) * bins(1) * bins(2); }
0152
0153   /// Bin-triple fast access
0154   ///
0155   /// - calculate the bin triple with one transform
0156   ///
0157   /// @param position is the 3D position to be evaluated
0158   ///
0159   /// @return is the bin value in 3D
0160   std::array<std::size_t, 3> binTriple(const Vector3& position) const {
0161     /// transform or not
0162     const Vector3 bPosition = m_itransform * position;
0163     // get the dimension
0164     std::size_t mdim = m_binningData.size();
0165     /// now get the bins
0166     std::size_t bin0 = m_binningData[0].searchGlobal(bPosition);
0167     std::size_t bin1 = mdim > 1 ? m_binningData[1].searchGlobal(bPosition) : 0;
0168     std::size_t bin2 = mdim > 2 ? m_binningData[2].searchGlobal(bPosition) : 0;
0169     /// return the triple
0170     return {{bin0, bin1, bin2}};
0171   }
0172
0173   /// Bin from a 3D vector (already in binning frame)
0174   ///
0175   /// @param position is the 3D position to be evaluated
0176   /// @param ba is the bin dimension
0177   ///
0178   /// @return is the bin value
0179   std::size_t bin(const Vector3& position, std::size_t ba = 0) const {
0180     if (ba >= m_binningData.size()) {
0181       return 0;
0182     }
0183     std::size_t bEval = m_binningData[ba].searchGlobal(m_itransform * position);
0184     return bEval;
0185   }
0186
0187   /// Bin from a 2D vector (following local parameters definitions)
0188   /// - no optional transform applied
0189   /// - USE WITH CARE !!
0190   ///
0191   /// You need to make sure that the local position is actually in the binning
0192   /// frame of the BinUtility
0193   ///
0194   /// @param lposition is the local position to be set
0195   /// @param ba is the bin dimension
0196   ///
0197   /// @return bin calculated from local
0198   std::size_t bin(const Vector2& lposition, std::size_t ba = 0) const {
0199     if (ba >= m_binningData.size()) {
0200       return 0;
0201     }
0202     return m_binningData[ba].searchLocal(lposition);
0203   }
0204
0205   /// Bin from a scalar (following local parameters definitions)
0206   /// - no optional transform applied
0207   /// - USE WITH CARE !!
0208   ///
0209   /// You need to make sure that the local position is actually in the binning
0210   /// frame of the BinUtility
0211   ///
0212   /// @param value is the scalar value to be evaluated
0213   /// @param ba is the bin dimension
0214   ///
0215   /// @return bin calculated from local
0216   std::size_t bin(float value, std::size_t ba = 0) const {
0217     if (ba >= m_binningData.size()) {
0218       return 0;
0219     }
0220     return m_binningData[ba].search(value);
0221   }
0222
0223   /// Return the other direction for fast interlinking
0224   ///
0225   /// @param position is the global position for the next search
0226   /// @param direction is the global position for the next search
0227   /// @param ba is the bin accessor
0228   ///
0229   /// @todo the
0230   ///
0231   /// @return the next bin
0232   int nextDirection(const Vector3& position, const Vector3& direction,
0233                     std::size_t ba = 0) const {
0234     if (ba >= m_binningData.size()) {
0235       return 0;
0236     }
0237     return m_binningData[ba].nextDirection(position, direction);
0238   }
0239
0240   /// Check if bin is inside from Vector2 - optional transform applied
0241   ///
0242   /// @param position is the global position to be evaluated
0243   /// @return is a boolean check
0244   bool inside(const Vector3& position) const {
0245     /// transform or not
0246     const Vector3& bPosition = m_itransform * position;
0247     return std::ranges::all_of(m_binningData, [&](const auto& bData) {
0248       return bData.inside(bPosition);
0249     });
0250   }
0251
0252   /// First bin maximal value
0253   /// @return the dimension of the binning data
0254   std::size_t dimensions() const { return m_binningData.size(); }
0255
0256   /// First bin maximal value
0257   ///
0258   /// @param ba is the binaccessor
0259   ///
0260   /// @return std::size_t is the maximal bin of the accessor entry
0261   std::size_t max(std::size_t ba = 0) const {
0262     if (ba >= m_binningData.size()) {
0263       return 0;
0264     }
0265     return (m_binningData[ba].bins() - 1);
0266   }
0267
0268   /// Number of bins
0269   ///
0270   /// @param ba is the binaccessor
0271   ///
0272   /// @return std::size_t is the bins of the accessor entry
0273   std::size_t bins(std::size_t ba) const {
0274     if (ba >= m_binningData.size()) {
0275       return 1;
0276     }
0277     return (m_binningData[ba].bins());
0278   }
0279
0280   /// Transform applied to global positions before lookup
0281   ///
0282   /// @return Shared pointer to transform
0283   const Transform3& transform() const { return m_transform; }
0284
0285   /// The type/value of the binning
0286   ///
0287   /// @param ba is the binaccessor
0288   ///
0289   /// @return the binning value of the accessor entry
0290   AxisDirection binningValue(std::size_t ba = 0) const {
0291     if (ba >= m_binningData.size()) {
0292       throw std::runtime_error{"Dimension out of bounds"};
0293     }
0294     return (m_binningData[ba].binvalue);
0295   }
0296
0297   /// Serialize the bin triple
0298   /// - this creates a simple std::size_t from a triple object
0299   ///
0300   /// @param bin is the bin to be serialized
0301   /// @return Serialized bin index as a single std::size_t value
0302   std::size_t serialize(const std::array<std::size_t, 3>& bin) const {
0303     std::size_t serializedBin = bin[0];
0304     if (m_binningData.size() == 2) {
0305       serializedBin += bin[1] * m_binningData[0].bins();
0306     } else if (m_binningData.size() == 3) {
0307       serializedBin +=
0308           (bin[1] * m_binningData[0].bins() * bin[2] * m_binningData[1].bins());
0309     }
0310     return serializedBin;
0311   }
0312
0313   /// Output Method for std::ostream, to be overloaded by child classes
0314   ///
0315   /// @param sl is the ostream to be dumped into
0316   /// @param indent the current indentation
0317   ///
0318   /// @return the input stream
0319   std::ostream& toStream(std::ostream& sl,
0320                          const std::string& indent = "") const {
0321     sl << indent << "BinUtility for " << m_binningData.size()
0322        << "- dimensional array:" << std::endl;
0323     for (auto [ibd, bd] : enumerate(m_binningData)) {
0324       sl << indent << "dimension     : " << ibd << std::endl;
0325       sl << bd.toString(indent) << std::endl;
0326     }
0327     return sl;
0328   }
0329
0330   /// Output into a string
0331   ///
0332   /// @param indent the current indentation
0333   ///
0334   /// @return a string with the stream information
0335   std::string toString(const std::string& indent = "") const {
0336     std::stringstream ss;
0337     toStream(ss, indent);
0338     return ss.str();
0339   }
0340
0341   /// Overload of << operator for std::ostream for debug output
0342   friend std::ostream& operator<<(std::ostream& sl, const BinUtility& bgen) {
0343     return bgen.toStream(sl);
0344   }
0345
0346  private:
0347   std::vector<BinningData> m_binningData;  /// vector of BinningData
0348   Transform3 m_transform;                  /// shared transform
0349   Transform3 m_itransform;                 /// unique inverse transform
0350 };
0351
0352 }  // namespace Acts