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