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