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 <array>
0019 #include <cstddef>
0020 #include <iostream>
0021 #include <iterator>
0022 #include <memory>
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   /// Copy constructor
0099   ///
0100   /// @param sbu is the source bin utility
0101   BinUtility(const BinUtility& sbu) = default;
0102
0103   /// Move constructor
0104   /// @param sbu is the source bin utility
0105   BinUtility(BinUtility&& sbu) = default;
0106
0107   /// Create from a DirectedProtoAxis
0108   ///
0109   /// @param dpAxis the DirectedProtoAxis to be used
0110   explicit BinUtility(const DirectedProtoAxis& dpAxis)
0111       : m_binningData(),
0112         m_transform(Transform3::Identity()),
0113         m_itransform(Transform3::Identity()) {
0114     m_binningData.reserve(3);
0115     m_binningData.emplace_back(dpAxis);
0116   }
0117
0118   /// Create from several DirectedProtoAxis objects
0119   ///
0120   /// @param dpAxes the DirectedProtoAxis to be used with axis directions
0121   explicit BinUtility(const std::vector<DirectedProtoAxis>& dpAxes)
0122       : m_binningData(),
0123         m_transform(Transform3::Identity()),
0124         m_itransform(Transform3::Identity()) {
0125     m_binningData.reserve(3);
0126     for (const auto& dpAxis : dpAxes) {
0127       m_binningData.emplace_back(dpAxis);
0128     }
0129   }
0130
0131   /// Assignment operator
0132   ///
0133   /// @param sbu is the source bin utility
0134   /// @return Reference to this BinUtility after assignment
0135   BinUtility& operator=(const BinUtility& sbu) {
0136     if (this != &sbu) {
0137       m_binningData = sbu.m_binningData;
0138       m_transform = sbu.m_transform;
0139       m_itransform = sbu.m_itransform;
0140     }
0141     return (*this);
0142   }
0143
0144   /// Move assignment operator
0145   /// @return Reference to this BinUtility after move assignment
0146   BinUtility& operator=(BinUtility&&) = default;
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   /// Virtual Destructor
0165   ~BinUtility() = default;
0166
0167   /// Equality operator
0168   /// @param other The other BinUtility to compare with
0169   /// @return True if the BinUtilities are equal, false otherwise
0170   bool operator==(const BinUtility& other) const {
0171     return (m_transform.isApprox(other.m_transform) &&
0172             m_binningData == other.binningData());
0173   }
0174
0175   /// Return the binning data vector
0176   /// @return Reference to the vector of binning data
0177   const std::vector<BinningData>& binningData() const { return m_binningData; }
0178
0179   /// Return the total number of bins
0180   /// @return Total number of bins across all dimensions
0181   std::size_t bins() const { return bins(0) * bins(1) * bins(2); }
0182
0183   /// Bin-triple fast access
0184   ///
0185   /// - calculate the bin triple with one transform
0186   ///
0187   /// @param position is the 3D position to be evaluated
0188   ///
0189   /// @return is the bin value in 3D
0190   std::array<std::size_t, 3> binTriple(const Vector3& position) const {
0191     /// transform or not
0192     const Vector3 bPosition = m_itransform * position;
0193     // get the dimension
0194     std::size_t mdim = m_binningData.size();
0195     /// now get the bins
0196     std::size_t bin0 = m_binningData[0].searchGlobal(bPosition);
0197     std::size_t bin1 = mdim > 1 ? m_binningData[1].searchGlobal(bPosition) : 0;
0198     std::size_t bin2 = mdim > 2 ? m_binningData[2].searchGlobal(bPosition) : 0;
0199     /// return the triple
0200     return {{bin0, bin1, bin2}};
0201   }
0202
0203   /// Bin from a 3D vector (already in binning frame)
0204   ///
0205   /// @param position is the 3D position to be evaluated
0206   /// @param ba is the bin dimension
0207   ///
0208   /// @return is the bin value
0209   std::size_t bin(const Vector3& position, std::size_t ba = 0) const {
0210     if (ba >= m_binningData.size()) {
0211       return 0;
0212     }
0213     std::size_t bEval = m_binningData[ba].searchGlobal(m_itransform * position);
0214     return bEval;
0215   }
0216
0217   /// Return the other direction for fast interlinking
0218   ///
0219   /// @param position is the global position for the next search
0220   /// @param direction is the global position for the next search
0221   /// @param ba is the bin accessor
0222   ///
0223   /// @todo the
0224   ///
0225   /// @return the next bin
0226   int nextDirection(const Vector3& position, const Vector3& direction,
0227                     std::size_t ba = 0) const {
0228     if (ba >= m_binningData.size()) {
0229       return 0;
0230     }
0231     return m_binningData[ba].nextDirection(position, direction);
0232   }
0233
0234   /// Bin from a 2D vector (following local parameters defintitions)
0235   /// - no optional transform applied
0236   /// - USE WITH CARE !!
0237   ///
0238   /// You need to make sure that the local position is actually in the binning
0239   /// frame of the BinUtility
0240   ///
0241   /// @param lposition is the local position to be set
0242   /// @param ba is the bin dimension
0243   ///
0244   /// @return bin calculated from local
0245   std::size_t bin(const Vector2& lposition, std::size_t ba = 0) const {
0246     if (ba >= m_binningData.size()) {
0247       return 0;
0248     }
0249     return m_binningData[ba].searchLocal(lposition);
0250   }
0251   /// Check if bin is inside from Vector2 - optional transform applied
0252   ///
0253   /// @param position is the global position to be evaluated
0254   /// @return is a boolean check
0255   bool inside(const Vector3& position) const {
0256     /// transform or not
0257     const Vector3& bPosition = m_itransform * position;
0258     // loop and break
0259     for (auto& bData : m_binningData) {
0260       if (!(bData.inside(bPosition))) {
0261         return false;
0262       }
0263     }
0264     // survived all the checks
0265     return true;
0266   }
0267
0268   /// First bin maximal value
0269   /// @return the dimension of the binning data
0270   std::size_t dimensions() const { return m_binningData.size(); }
0271
0272   /// First bin maximal value
0273   ///
0274   /// @param ba is the binaccessor
0275   ///
0276   /// @return std::size_t is the maximal bin of the accessor entry
0277   std::size_t max(std::size_t ba = 0) const {
0278     if (ba >= m_binningData.size()) {
0279       return 0;
0280     }
0281     return (m_binningData[ba].bins() - 1);
0282   }
0283
0284   /// Number of bins
0285   ///
0286   /// @param ba is the binaccessor
0287   ///
0288   /// @return std::size_t is the bins of the accessor entry
0289   std::size_t bins(std::size_t ba) const {
0290     if (ba >= m_binningData.size()) {
0291       return 1;
0292     }
0293     return (m_binningData[ba].bins());
0294   }
0295
0296   /// Transform applied to global positions before lookup
0297   ///
0298   /// @return Shared pointer to transform
0299   const Transform3& transform() const { return m_transform; }
0300
0301   /// The type/value of the binning
0302   ///
0303   /// @param ba is the binaccessor
0304   ///
0305   /// @return the binning value of the accessor entry
0306   AxisDirection binningValue(std::size_t ba = 0) const {
0307     if (ba >= m_binningData.size()) {
0308       throw std::runtime_error{"Dimension out of bounds"};
0309     }
0310     return (m_binningData[ba].binvalue);
0311   }
0312
0313   /// Serialize the bin triple
0314   /// - this creates a simple std::size_t from a triple object
0315   ///
0316   /// @param bin is the bin to be serialized
0317   /// @return Serialized bin index as a single std::size_t value
0318   std::size_t serialize(const std::array<std::size_t, 3>& bin) const {
0319     std::size_t serializedBin = bin[0];
0320     if (m_binningData.size() == 2) {
0321       serializedBin += bin[1] * m_binningData[0].bins();
0322     } else if (m_binningData.size() == 3) {
0323       serializedBin +=
0324           (bin[1] * m_binningData[0].bins() * bin[2] * m_binningData[1].bins());
0325     }
0326     return serializedBin;
0327   }
0328
0329   /// Output Method for std::ostream, to be overloaded by child classes
0330   ///
0331   /// @param sl is the ostream to be dumped into
0332   /// @param indent the current indentation
0333   ///
0334   /// @return the input stream
0335   std::ostream& toStream(std::ostream& sl,
0336                          const std::string& indent = "") const {
0337     sl << indent << "BinUtility for " << m_binningData.size()
0338        << "- dimensional array:" << std::endl;
0339     for (auto [ibd, bd] : enumerate(m_binningData)) {
0340       sl << indent << "dimension     : " << ibd << std::endl;
0341       sl << bd.toString(indent) << std::endl;
0342     }
0343     return sl;
0344   }
0345
0346   /// Output into a string
0347   ///
0348   /// @param indent the current indentation
0349   ///
0350   /// @return a string with the stream information
0351   std::string toString(const std::string& indent = "") const {
0352     std::stringstream ss;
0353     toStream(ss, indent);
0354     return ss.str();
0355   }
0356
0357   /// Overload of << operator for std::ostream for debug output
0358   friend std::ostream& operator<<(std::ostream& sl, const BinUtility& bgen) {
0359     return bgen.toStream(sl);
0360   }
0361
0362  private:
0363   std::vector<BinningData> m_binningData;  /// vector of BinningData
0364   Transform3 m_transform;                  /// shared transform
0365   Transform3 m_itransform;                 /// unique inverse transform
0366 };
0367
0368 }  // namespace Acts