RHistEngine.hxx

Go to the documentation of this file.

/// \file
/// \warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes.
/// Feedback is welcome!
 
#ifndef ROOT_RHistEngine
#define ROOT_RHistEngine
 
#include "RAxes.hxx"
#include "RAxisVariant.hxx"
#include "RBinIndex.hxx"
#include "RBinIndexMultiDimRange.hxx"
#include "RHistUtils.hxx"
#include "RLinearizedIndex.hxx"
#include "RRegularAxis.hxx"
#include "RSliceBinIndexMapper.hxx"
#include "RSliceSpec.hxx"
#include "RWeight.hxx"
 
#include <array>
#include <atomic>
#include <cassert>
#include <cstddef>
#include <cstdint>
#include <cstring>
#include <initializer_list>
#include <stdexcept>
#include <tuple>
#include <type_traits>
#include <utility>
#include <vector>
 
class TBuffer;
 
namespace ROOT {
namespace Experimental {
 
// forward declaration for friend declaration
template <typename T>
class RHist;
 
/**
A histogram data structure to bin data along multiple dimensions.
 
Every call to \ref Fill(const A &... args) "Fill" bins the data according to the axis configuration and increments the
bin content:
\code
ROOT::Experimental::RHistEngine<int> hist(10, {5, 15});
hist.Fill(8.5);
// hist.GetBinContent(ROOT::Experimental::RBinIndex(3)) will return 1
\endcode
 
The class is templated on the bin content type. For counting, as in the example above, it may be an integral type such
as `int` or `long`. Narrower types such as `unsigned char` or `short` are supported, but may overflow due to their
limited range and must be used with care. For weighted filling, the bin content type must not be an integral type, but
a floating-point type such as `float` or `double`, or the special type RBinWithError. Note that `float` has a limited
significand precision of 24 bits.
 
An object can have arbitrary dimensionality determined at run-time. The axis configuration is passed as a vector of
RAxisVariant:
\code
std::vector<ROOT::Experimental::RAxisVariant> axes;
axes.push_back(ROOT::Experimental::RRegularAxis(10, {5, 15}));
axes.push_back(ROOT::Experimental::RVariableBinAxis({1, 10, 100, 1000}));
ROOT::Experimental::RHistEngine<int> hist(axes);
// hist.GetNDimensions() will return 2
\endcode
 
\warning This is part of the %ROOT 7 prototype! It will change without notice. It might trigger earthquakes.
Feedback is welcome!
*/
template <typename BinContentType>
class RHistEngine final {
   // For conversion, all other template instantiations must be a friend.
   template <typename U>
   friend class RHistEngine;
 
   // For slicing, RHist needs to call SliceImpl.
   friend class RHist<BinContentType>;
 
   friend class RProfile;
   friend class RProfileFillContext;
 
   /// The axis configuration for this histogram. Relevant methods are forwarded from the public interface.
   Internal::RAxes fAxes;
   /// The bin contents for this histogram
   std::vector<BinContentType> fBinContents;
 
   /// Flag to pause filling while a snapshot is ongoing
   mutable std::atomic<bool> fSnapshotInProgress{false}; //!
 
public:
   /// Construct a histogram engine.
   ///
   /// \param[in] axes the axis objects, must have size > 0
   explicit RHistEngine(std::vector<RAxisVariant> axes) : fAxes(std::move(axes))
   {
      fBinContents.resize(fAxes.ComputeTotalNBins());
   }
 
   /// Construct a histogram engine.
   ///
   /// Note that there is no perfect forwarding of the axis objects. If that is needed, use the
   /// \ref RHistEngine(std::vector<RAxisVariant> axes) "overload accepting a std::vector".
   ///
   /// \param[in] axes the axis objects, must have size > 0
   explicit RHistEngine(std::initializer_list<RAxisVariant> axes) : RHistEngine(std::vector(axes)) {}
 
   /// Construct a histogram engine.
   ///
   /// Note that there is no perfect forwarding of the axis objects. If that is needed, use the
   /// \ref RHistEngine(std::vector<RAxisVariant> axes) "overload accepting a std::vector".
   ///
   /// \param[in] axis1 the first axis object
   /// \param[in] axes the remaining axis objects
   template <typename... Axes>
   explicit RHistEngine(const RAxisVariant &axis1, const Axes &...axes)
      : RHistEngine(std::vector<RAxisVariant>{axis1, axes...})
   {
   }
 
   /// Construct a one-dimensional histogram engine with a regular axis.
   ///
   /// \param[in] nNormalBins the number of normal bins, must be > 0
   /// \param[in] interval the axis interval (lower end inclusive, upper end exclusive)
   /// \par See also
   /// the \ref RRegularAxis::RRegularAxis(std::uint64_t nNormalBins, std::pair<double, double> interval, bool
   /// enableFlowBins) "constructor of RRegularAxis"
   RHistEngine(std::uint64_t nNormalBins, std::pair<double, double> interval)
      : RHistEngine(std::vector<RAxisVariant>{RRegularAxis(nNormalBins, interval)})
   {
   }
 
   /// The copy constructor is deleted.
   ///
   /// Copying all bin contents can be an expensive operation, depending on the number of bins. If required, users can
   /// explicitly call Clone().
   RHistEngine(const RHistEngine &) = delete;
   /// Efficiently move construct a histogram engine.
   ///
   /// After this operation, the moved-from object is invalid.
   RHistEngine(RHistEngine &&rhs) noexcept : fAxes(std::move(rhs.fAxes)), fBinContents(std::move(rhs.fBinContents)) {}
 
   /// The copy assignment operator is deleted.
   ///
   /// Copying all bin contents can be an expensive operation, depending on the number of bins. If required, users can
   /// explicitly call Clone().
   RHistEngine &operator=(const RHistEngine &) = delete;
   /// Efficiently move a histogram engine.
   ///
   /// After this operation, the moved-from object is invalid.
   RHistEngine &operator=(RHistEngine &&rhs) noexcept
   {
      std::swap(fAxes, rhs.fAxes);
      std::swap(fBinContents, rhs.fBinContents);
      return *this;
   }
 
   ~RHistEngine() = default;
 
   /// \name Accessors
   /// \{
 
   const std::vector<RAxisVariant> &GetAxes() const { return fAxes.Get(); }
   std::size_t GetNDimensions() const { return fAxes.GetNDimensions(); }
   std::uint64_t GetTotalNBins() const { return fBinContents.size(); }
 
   /// Get the content of a single bin.
   ///
   /// \code
   /// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
   /// std::array<ROOT::Experimental::RBinIndex, 2> indices = {3, 5};
   /// int content = hist.GetBinContent(indices);
   /// \endcode
   ///
   /// \note Compared to TH1 conventions, the first normal bin has index 0 and underflow and overflow bins are special
   /// values. See also the class documentation of RBinIndex.
   ///
   /// Throws an exception if the number of indices does not match the axis configuration or the bin is not found.
   ///
   /// \param[in] indices the array of indices for each axis
   /// \return the bin content
   /// \par See also
   /// the \ref GetBinContent(const A &... args) const "variadic function template overload" accepting arguments
   /// directly
   template <std::size_t N>
   const BinContentType &GetBinContent(const std::array<RBinIndex, N> &indices) const
   {
      // We could rely on RAxes::ComputeGlobalIndex to check the number of arguments, but its exception message might
      // be confusing for users.
      if (N != GetNDimensions()) {
         throw std::invalid_argument("invalid number of indices passed to GetBinContent");
      }
      RLinearizedIndex index = fAxes.ComputeGlobalIndex(indices);
      if (!index.fValid) {
         throw std::invalid_argument("bin not found in GetBinContent");
      }
      assert(index.fIndex < fBinContents.size());
      return fBinContents[index.fIndex];
   }
 
   /// Get the content of a single bin.
   ///
   /// \code
   /// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
   /// std::vector<ROOT::Experimental::RBinIndex> indices = {3, 5};
   /// int content = hist.GetBinContent(indices);
   /// \endcode
   ///
   /// \note Compared to TH1 conventions, the first normal bin has index 0 and underflow and overflow bins are special
   /// values. See also the class documentation of RBinIndex.
   ///
   /// Throws an exception if the number of indices does not match the axis configuration or the bin is not found.
   ///
   /// \param[in] indices the vector of indices for each axis
   /// \return the bin content
   /// \par See also
   /// the \ref GetBinContent(const A &... args) const "variadic function template overload" accepting arguments
   /// directly
   const BinContentType &GetBinContent(const std::vector<RBinIndex> &indices) const
   {
      // We could rely on RAxes::ComputeGlobalIndex to check the number of arguments, but its exception message might
      // be confusing for users.
      if (indices.size() != GetNDimensions()) {
         throw std::invalid_argument("invalid number of indices passed to GetBinContent");
      }
      RLinearizedIndex index = fAxes.ComputeGlobalIndex(indices);
      if (!index.fValid) {
         throw std::invalid_argument("bin not found in GetBinContent");
      }
      assert(index.fIndex < fBinContents.size());
      return fBinContents[index.fIndex];
   }
 
   /// Get the content of a single bin.
   ///
   /// \code
   /// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
   /// int content = hist.GetBinContent(ROOT::Experimental::RBinIndex(3), ROOT::Experimental::RBinIndex(5));
   /// // ... or construct the RBinIndex arguments implicitly from integers:
   /// content = hist.GetBinContent(3, 5);
   /// \endcode
   ///
   /// \note Compared to TH1 conventions, the first normal bin has index 0 and underflow and overflow bins are special
   /// values. See also the class documentation of RBinIndex.
   ///
   /// Throws an exception if the number of arguments does not match the axis configuration or the bin is not found.
   ///
   /// \param[in] args the arguments for each axis
   /// \return the bin content
   /// \par See also
   /// the function overloads accepting \ref GetBinContent(const std::array<RBinIndex, N> &indices) const "`std::array`"
   /// or \ref GetBinContent(const std::vector<RBinIndex> &indices) const "`std::vector`"
   template <typename... A>
   const BinContentType &GetBinContent(const A &...args) const
   {
      std::array<RBinIndex, sizeof...(A)> indices{args...};
      return GetBinContent(indices);
   }
 
   /// Get the multidimensional range of all bins.
   ///
   /// \return the multidimensional range
   RBinIndexMultiDimRange GetFullMultiDimRange() const { return fAxes.GetFullMultiDimRange(); }
 
   /// Set the content of a single bin.
   ///
   /// \code
   /// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
   /// std::array<ROOT::Experimental::RBinIndex, 2> indices = {3, 5};
   /// int value = /* ... */;
   /// hist.SetBinContent(indices, value);
   /// \endcode
   ///
   /// \note Compared to TH1 conventions, the first normal bin has index 0 and underflow and overflow bins are special
   /// values. See also the class documentation of RBinIndex.
   ///
   /// Throws an exception if the number of indices does not match the axis configuration or the bin is not found.
   ///
   /// \param[in] indices the array of indices for each axis
   /// \param[in] value the new value of the bin content
   /// \par See also
   /// the \ref SetBinContent(const A &... args) "variadic function template overload" accepting arguments directly
   template <std::size_t N, typename V>
   void SetBinContent(const std::array<RBinIndex, N> &indices, const V &value)
   {
      // We could rely on RAxes::ComputeGlobalIndex to check the number of arguments, but its exception message might
      // be confusing for users.
      if (N != GetNDimensions()) {
         throw std::invalid_argument("invalid number of indices passed to SetBinContent");
      }
      RLinearizedIndex index = fAxes.ComputeGlobalIndex(indices);
      if (!index.fValid) {
         throw std::invalid_argument("bin not found in SetBinContent");
      }
      assert(index.fIndex < fBinContents.size());
      // To allow conversion, we have to accept value with a template type V to capture any argument. Otherwise it would
      // select the variadic function template...
      fBinContents[index.fIndex] = value;
   }
 
   /// \}
   // End the group to ensure that all contained member functions are public.
 
private:
   template <typename... A, std::size_t... I>
   void SetBinContentImpl(const std::tuple<A...> &args, std::index_sequence<I...>)
   {
      std::array<RBinIndex, sizeof...(A) - 1> indices{std::get<I>(args)...};
      SetBinContent(indices, std::get<sizeof...(A) - 1>(args));
   }
 
public:
   /// \name Accessors
   /// \{
 
   /// Set the content of a single bin.
   ///
   /// \code
   /// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
   /// int value = /* ... */;
   /// hist.SetBinContent(ROOT::Experimental::RBinIndex(3), ROOT::Experimental::RBinIndex(5), value);
   /// // ... or construct the RBinIndex arguments implicitly from integers:
   /// hist.SetBinContent(3, 5, value);
   /// \endcode
   ///
   /// \note Compared to TH1 conventions, the first normal bin has index 0 and underflow and overflow bins are special
   /// values. See also the class documentation of RBinIndex.
   ///
   /// Throws an exception if the number of arguments does not match the axis configuration or the bin is not found.
   ///
   /// \param[in] args the arguments for each axis and the new value of the bin content
   /// \par See also
   /// the \ref SetBinContent(const std::array<RBinIndex, N> &indices, const V &value) "function overload" accepting
   /// `std::array`
   template <typename... A>
   void SetBinContent(const A &...args)
   {
      auto t = std::forward_as_tuple(args...);
      SetBinContentImpl(t, std::make_index_sequence<sizeof...(A) - 1>());
   }
 
   /// \}
 
   /// Whether this histogram engine type supports weighted filling.
   static constexpr bool SupportsWeightedFilling = !std::is_integral_v<BinContentType>;
 
   // SupportsWeightedFilling is not included because it is static, which would mess up the subgrouping below "Public
   // Member Functions".
   /// \name Filling
   /// \{
 
   /// Fill an entry into the histogram.
   ///
   /// \code
   /// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
   /// auto args = std::make_tuple(8.5, 10.5);
   /// hist.Fill(args);
   /// \endcode
   ///
   /// If one of the arguments is outside the corresponding axis and flow bins are disabled, the entry will be silently
   /// discarded.
   ///
   /// Throws an exception if the number of arguments does not match the axis configuration, or if an argument cannot be
   /// converted for the axis type at run-time.
   ///
   /// \param[in] args the arguments for each axis
   /// \par See also
   /// the \ref Fill(const A &... args) "variadic function template overload" accepting arguments directly and the
   /// \ref Fill(const std::tuple<A...> &args, RWeight weight) "overload for weighted filling"
   template <typename... A>
   void Fill(const std::tuple<A...> &args)
   {
      // We could rely on RAxes::ComputeGlobalIndex to check the number of arguments, but its exception message might
      // be confusing for users.
      if (sizeof...(A) != GetNDimensions()) {
         throw std::invalid_argument("invalid number of arguments to Fill");
      }
      RLinearizedIndex index = fAxes.ComputeGlobalIndexImpl<sizeof...(A)>(args);
      if (index.fValid) {
         assert(index.fIndex < fBinContents.size());
         fBinContents[index.fIndex]++;
      }
   }
 
   /// Fill an entry into the histogram with a weight.
   ///
   /// This overload is not available for integral bin content types (see \ref SupportsWeightedFilling).
   ///
   /// \code
   /// ROOT::Experimental::RHistEngine<float> hist({/* two dimensions */});
   /// auto args = std::make_tuple(8.5, 10.5);
   /// hist.Fill(args, ROOT::Experimental::RWeight(0.8));
   /// \endcode
   ///
   /// If one of the arguments is outside the corresponding axis and flow bins are disabled, the entry will be silently
   /// discarded.
   ///
   /// Throws an exception if the number of arguments does not match the axis configuration, or if an argument cannot be
   /// converted for the axis type at run-time.
   ///
   /// \param[in] args the arguments for each axis
   /// \param[in] weight the weight for this entry
   /// \par See also
   /// the \ref Fill(const A &... args) "variadic function template overload" accepting arguments directly and the
   /// \ref Fill(const std::tuple<A...> &args) "overload for unweighted filling"
   template <typename... A>
   void Fill(const std::tuple<A...> &args, RWeight weight)
   {
      static_assert(SupportsWeightedFilling, "weighted filling is not supported for integral bin content types");
 
      // We could rely on RAxes::ComputeGlobalIndex to check the number of arguments, but its exception message might
      // be confusing for users.
      if (sizeof...(A) != GetNDimensions()) {
         throw std::invalid_argument("invalid number of arguments to Fill");
      }
      RLinearizedIndex index = fAxes.ComputeGlobalIndexImpl<sizeof...(A)>(args);
      if (index.fValid) {
         assert(index.fIndex < fBinContents.size());
         fBinContents[index.fIndex] += weight.fValue;
      }
   }
 
   /// \}
   // End the group to ensure that all contained member functions are public.
 
private:
   // Also used by RProfile::Fill(const A &...args) - similar to the variadic RHistEngine::Fill(const A &...args) below,
   // is has all arguments in the forwarded std::tuple and needs to explicitly specify how many of them should be used
   // by RAxes::ComputeGlobalIndexImpl<N>(args).
   template <std::size_t N, typename... A, typename W>
   void FillImpl(const std::tuple<A...> &args, const W &weight)
   {
      RLinearizedIndex index = fAxes.ComputeGlobalIndexImpl<N>(args);
      if (index.fValid) {
         assert(index.fIndex < fBinContents.size());
         fBinContents[index.fIndex] += weight;
      }
   }
 
public:
   /// \name Filling
   /// \{
 
   /// Fill an entry into the histogram with a user-defined weight.
   ///
   /// This overload is only available for user-defined bin content types.
   ///
   /// If one of the arguments is outside the corresponding axis and flow bins are disabled, the entry will be silently
   /// discarded.
   ///
   /// Throws an exception if the number of arguments does not match the axis configuration, or if an argument cannot be
   /// converted for the axis type at run-time.
   ///
   /// \param[in] args the arguments for each axis
   /// \param[in] weight the weight for this entry
   template <typename... A, typename W>
   void Fill(const std::tuple<A...> &args, const W &weight)
   {
      static_assert(std::is_class_v<BinContentType>,
                    "user-defined weight types are only supported for user-defined bin content types");
 
      // We could rely on RAxes::ComputeGlobalIndex to check the number of arguments, but its exception message might
      // be confusing for users.
      if (sizeof...(A) != GetNDimensions()) {
         throw std::invalid_argument("invalid number of arguments to Fill");
      }
      FillImpl<sizeof...(A)>(args, weight);
   }
 
   /// Fill an entry into the histogram.
   ///
   /// \code
   /// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
   /// hist.Fill(8.5, 10.5);
   /// \endcode
   ///
   /// For weighted filling, pass an RWeight as the last argument:
   /// \code
   /// ROOT::Experimental::RHistEngine<float> hist({/* two dimensions */});
   /// hist.Fill(8.5, 10.5, ROOT::Experimental::RWeight(0.8));
   /// \endcode
   /// This is not available for integral bin content types (see \ref SupportsWeightedFilling).
   ///
   /// If one of the arguments is outside the corresponding axis and flow bins are disabled, the entry will be silently
   /// discarded.
   ///
   /// Throws an exception if the number of arguments does not match the axis configuration, or if an argument cannot be
   /// converted for the axis type at run-time.
   ///
   /// \param[in] args the arguments for each axis
   /// \par See also
   /// the function overloads accepting `std::tuple` \ref Fill(const std::tuple<A...> &args) "for unweighted filling"
   /// and \ref Fill(const std::tuple<A...> &args, RWeight) "for weighted filling"
   template <typename... A>
   void Fill(const A &...args)
   {
      static_assert(sizeof...(A) >= 1, "need at least one argument to Fill");
      if constexpr (sizeof...(A) >= 1) {
         auto t = std::forward_as_tuple(args...);
         if constexpr (std::is_same_v<typename Internal::LastType<A...>::type, RWeight>) {
            static_assert(SupportsWeightedFilling, "weighted filling is not supported for integral bin content types");
            static constexpr std::size_t N = sizeof...(A) - 1;
            if (N != GetNDimensions()) {
               throw std::invalid_argument("invalid number of arguments to Fill");
            }
            RWeight weight = std::get<N>(t);
            RLinearizedIndex index = fAxes.ComputeGlobalIndexImpl<N>(t);
            if (index.fValid) {
               assert(index.fIndex < fBinContents.size());
               fBinContents[index.fIndex] += weight.fValue;
            }
         } else {
            Fill(t);
         }
      }
   }
 
   /// Fill an entry into the histogram using atomic instructions.
   ///
   /// \param[in] args the arguments for each axis
   /// \see Fill(const std::tuple<A...> &args)
   template <typename... A>
   void FillAtomic(const std::tuple<A...> &args)
   {
      while (fSnapshotInProgress.load(std::memory_order_relaxed)) {
         // Spin while a snapshot is running
      }
 
      // We could rely on RAxes::ComputeGlobalIndex to check the number of arguments, but its exception message might
      // be confusing for users.
      if (sizeof...(A) != GetNDimensions()) {
         throw std::invalid_argument("invalid number of arguments to Fill");
      }
      RLinearizedIndex index = fAxes.ComputeGlobalIndexImpl<sizeof...(A)>(args);
      if (index.fValid) {
         assert(index.fIndex < fBinContents.size());
         Internal::AtomicIncRelease(&fBinContents[index.fIndex]);
      }
   }
 
   /// Fill an entry into the histogram with a weight using atomic instructions.
   ///
   /// This overload is not available for integral bin content types (see \ref SupportsWeightedFilling).
   ///
   /// \param[in] args the arguments for each axis
   /// \param[in] weight the weight for this entry
   /// \see Fill(const std::tuple<A...> &args, RWeight weight)
   template <typename... A>
   void FillAtomic(const std::tuple<A...> &args, RWeight weight)
   {
      static_assert(SupportsWeightedFilling, "weighted filling is not supported for integral bin content types");
 
      while (fSnapshotInProgress.load(std::memory_order_relaxed)) {
         // Spin while a snapshot is running
      }
 
      // We could rely on RAxes::ComputeGlobalIndex to check the number of arguments, but its exception message might
      // be confusing for users.
      if (sizeof...(A) != GetNDimensions()) {
         throw std::invalid_argument("invalid number of arguments to Fill");
      }
      RLinearizedIndex index = fAxes.ComputeGlobalIndexImpl<sizeof...(A)>(args);
      if (index.fValid) {
         assert(index.fIndex < fBinContents.size());
         Internal::AtomicAddRelease(&fBinContents[index.fIndex], weight.fValue);
      }
   }
 
   /// \}
   // End the group to ensure that all contained member functions are public.
 
private:
   // Also used by RProfileFillContext::Fill(const A &...args) - similar to the variadic
   // RHistEngine::FillAtomic(const A &...args) below, is has all arguments in the forwarded std::tuple and needs to
   // explicitly specify how many of them should be used by RAxes::ComputeGlobalIndexImpl<N>(args).
   template <std::size_t N, typename... A, typename W>
   void FillAtomicImpl(const std::tuple<A...> &args, const W &weight)
   {
      RLinearizedIndex index = fAxes.ComputeGlobalIndexImpl<N>(args);
      if (index.fValid) {
         assert(index.fIndex < fBinContents.size());
         Internal::AtomicAddRelease(&fBinContents[index.fIndex], weight);
      }
   }
 
public:
   /// \name Filling
   /// \{
 
   /// Fill an entry into the histogram with a user-defined weight using atomic instructions.
   ///
   /// This overload is only available for user-defined bin content types.
   ///
   /// \param[in] args the arguments for each axis
   /// \param[in] weight the weight for this entry
   /// \see Fill(const std::tuple<A...> &args, const W &weight)
   template <typename... A, typename W>
   void FillAtomic(const std::tuple<A...> &args, const W &weight)
   {
      static_assert(std::is_class_v<BinContentType>,
                    "user-defined weight types are only supported for user-defined bin content types");
 
      while (fSnapshotInProgress.load(std::memory_order_relaxed)) {
         // Spin while a snapshot is running
      }
 
      // We could rely on RAxes::ComputeGlobalIndex to check the number of arguments, but its exception message might
      // be confusing for users.
      if (sizeof...(A) != GetNDimensions()) {
         throw std::invalid_argument("invalid number of arguments to Fill");
      }
      FillAtomicImpl<sizeof...(A)>(args, weight);
   }
 
   /// Fill an entry into the histogram using atomic instructions.
   ///
   /// \param[in] args the arguments for each axis
   /// \see Fill(const A &...args)
   template <typename... A>
   void FillAtomic(const A &...args)
   {
      static_assert(sizeof...(A) >= 1, "need at least one argument to Fill");
      if constexpr (sizeof...(A) >= 1) {
         while (fSnapshotInProgress.load(std::memory_order_relaxed)) {
            // Spin while a snapshot is running
         }
 
         auto t = std::forward_as_tuple(args...);
         if constexpr (std::is_same_v<typename Internal::LastType<A...>::type, RWeight>) {
            static_assert(SupportsWeightedFilling, "weighted filling is not supported for integral bin content types");
            static constexpr std::size_t N = sizeof...(A) - 1;
            if (N != GetNDimensions()) {
               throw std::invalid_argument("invalid number of arguments to Fill");
            }
            RWeight weight = std::get<N>(t);
            RLinearizedIndex index = fAxes.ComputeGlobalIndexImpl<N>(t);
            if (index.fValid) {
               assert(index.fIndex < fBinContents.size());
               Internal::AtomicAddRelease(&fBinContents[index.fIndex], weight.fValue);
            }
         } else {
            FillAtomic(t);
         }
      }
   }
 
   /// \}
   /// \name Operations
   /// \{
 
   /// Add all bin contents of another histogram.
   ///
   /// Throws an exception if the axes configurations are not identical.
   ///
   /// \param[in] other another histogram
   void Add(const RHistEngine &other)
   {
      if (fAxes != other.fAxes) {
         throw std::invalid_argument("axes configurations not identical in Add");
      }
      for (std::size_t i = 0; i < fBinContents.size(); i++) {
         fBinContents[i] += other.fBinContents[i];
      }
   }
 
   /// Add all bin contents of another histogram using atomic instructions.
   ///
   /// Throws an exception if the axes configurations are not identical.
   ///
   /// \param[in] other another histogram that must not be modified during the operation
   void AddAtomic(const RHistEngine &other)
   {
      if (fAxes != other.fAxes) {
         throw std::invalid_argument("axes configurations not identical in AddAtomic");
      }
      for (std::size_t i = 0; i < fBinContents.size(); i++) {
         Internal::AtomicAdd(&fBinContents[i], other.fBinContents[i]);
      }
   }
 
   /// Clear all bin contents.
   void Clear()
   {
      for (std::size_t i = 0; i < fBinContents.size(); i++) {
         fBinContents[i] = {};
      }
   }
 
   /// Clone this histogram engine.
   ///
   /// Copying all bin contents can be an expensive operation, depending on the number of bins.
   ///
   /// \return the cloned object
   RHistEngine Clone() const
   {
      RHistEngine h(fAxes.Get());
      for (std::size_t i = 0; i < fBinContents.size(); i++) {
         h.fBinContents[i] = fBinContents[i];
      }
      return h;
   }
 
   /// Convert this histogram engine to a different bin content type.
   ///
   /// There is no bounds checking to make sure that the converted values can be represented. Note that it is not
   /// possible to convert to RBinWithError since the information about individual weights has been lost since filling.
   ///
   /// Converting all bin contents can be an expensive operation, depending on the number of bins.
   ///
   /// \return the converted object
   template <typename U>
   RHistEngine<U> Convert() const
   {
      RHistEngine<U> h(fAxes.Get());
      for (std::size_t i = 0; i < fBinContents.size(); i++) {
         h.fBinContents[i] = static_cast<U>(fBinContents[i]);
      }
      return h;
   }
 
   /// Scale all histogram bin contents.
   ///
   /// This method is not available for integral bin content types.
   ///
   /// \param[in] factor the scale factor
   void Scale(double factor)
   {
      static_assert(!std::is_integral_v<BinContentType>, "scaling is not supported for integral bin content types");
      for (std::size_t i = 0; i < fBinContents.size(); i++) {
         fBinContents[i] *= factor;
      }
   }
 
   /// \}
   // End the group to ensure that all contained member functions are public.
 
private:
   RHistEngine SliceImpl(const std::vector<RSliceSpec> &sliceSpecs, bool &dropped) const
   {
      if (sliceSpecs.size() != GetNDimensions()) {
         throw std::invalid_argument("invalid number of specifications passed to Slice");
      }
 
      // Slice the axes.
      std::vector<RAxisVariant> axes;
      for (std::size_t i = 0; i < sliceSpecs.size(); i++) {
         // A sum operation makes the dimension disappear.
         if (sliceSpecs[i].GetOperationSum() == nullptr) {
            axes.push_back(fAxes.Get()[i].Slice(sliceSpecs[i]));
         }
      }
      if (axes.empty()) {
         throw std::invalid_argument("summing across all dimensions is not supported");
      }
 
      RHistEngine sliced(std::move(axes));
 
      // Create the helper objects to map the bin contents to the sliced histogram.
      Internal::RSliceBinIndexMapper mapper(sliceSpecs);
      assert(mapper.GetMappedDimensionality() == sliced.GetNDimensions());
      std::vector<RBinIndex> mappedIndices(mapper.GetMappedDimensionality());
 
      auto origRange = fAxes.GetFullMultiDimRange();
      auto origRangeIt = origRange.begin();
 
      for (std::size_t i = 0; i < fBinContents.size(); i++) {
         const auto &origIndices = *origRangeIt;
#ifndef NDEBUG
         // Verify that the original indices correspond to the iteration variable.
         RLinearizedIndex origIndex = fAxes.ComputeGlobalIndex(origIndices);
         assert(origIndex.fValid);
         assert(origIndex.fIndex == i);
#endif
 
         bool success = mapper.Map(origIndices, mappedIndices);
         if (success) {
            RLinearizedIndex mappedIndex = sliced.fAxes.ComputeGlobalIndex(mappedIndices);
            assert(mappedIndex.fValid);
            sliced.fBinContents[mappedIndex.fIndex] += fBinContents[i];
         } else {
            dropped = true;
         }
         ++origRangeIt;
      }
 
      return sliced;
   }
 
public:
   /// \name Operations
   /// \{
 
   /// Slice this histogram with an RSliceSpec per dimension.
   ///
   /// With a range, only the specified bins are retained. All other bin contents are transferred to the underflow and
   /// overflow bins:
   /// \code
   /// ROOT::Experimental::RHistEngine<int> hist(/* one dimension */);
   /// // Fill the histogram with a number of entries...
   /// auto sliced = hist.Slice({hist.GetAxes()[0].GetNormalRange(1, 5)});
   /// // The returned histogram will have 4 normal bins, an underflow and an overflow bin.
   /// \endcode
   ///
   /// Slicing can also perform operations per dimension, see RSliceSpec. RSliceSpec::ROperationRebin allows to rebin
   /// the histogram axis, grouping a number of normal bins into a new one:
   /// \code
   /// ROOT::Experimental::RHistEngine<int> hist(/* one dimension */);
   /// // Fill the histogram with a number of entries...
   /// auto rebinned = hist.Slice(ROOT::Experimental::RSliceSpec::ROperationRebin(2));
   /// // The returned histogram has groups of two normal bins merged.
   /// \endcode
   ///
   /// RSliceSpec::ROperationSum sums the bin contents along that axis, which allows to project to a lower-dimensional
   /// histogram:
   /// \code
   /// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
   /// // Fill the histogram with a number of entries...
   /// auto projected = hist.Slice(ROOT::Experimental::RSliceSpec{}, ROOT::Experimental::RSliceSpec::ROperationSum{});
   /// // The returned histogram has one dimension, with bin contents summed along the second axis.
   /// \endcode
   /// Note that it is not allowed to sum along all histogram axes because the return value would be a scalar.
   ///
   /// Ranges and operations can be combined. In that case, the range is applied before the operation.
   ///
   /// \param[in] sliceSpecs the slice specifications for each axis
   /// \return the sliced histogram
   /// \par See also
   /// the \ref Slice(const A &... args) const "variadic function template overload" accepting arguments directly
   RHistEngine Slice(const std::vector<RSliceSpec> &sliceSpecs) const
   {
      bool dropped = false;
      return SliceImpl(sliceSpecs, dropped);
   }
 
   /// Slice this histogram with an RSliceSpec per dimension.
   ///
   /// With a range, only the specified bins are retained. All other bin contents are transferred to the underflow and
   /// overflow bins:
   /// \code
   /// ROOT::Experimental::RHistEngine<int> hist(/* one dimension */);
   /// // Fill the histogram with a number of entries...
   /// auto sliced = hist.Slice(hist.GetAxes()[0].GetNormalRange(1, 5));
   /// // The returned histogram will have 4 normal bins, an underflow and an overflow bin.
   /// \endcode
   ///
   /// Slicing can also perform operations per dimension, see RSliceSpec. RSliceSpec::ROperationRebin allows to rebin
   /// the histogram axis, grouping a number of normal bins into a new one:
   /// \code
   /// ROOT::Experimental::RHistEngine<int> hist(/* one dimension */);
   /// // Fill the histogram with a number of entries...
   /// auto rebinned = hist.Slice(ROOT::Experimental::RSliceSpec::ROperationRebin(2));
   /// // The returned histogram has groups of two normal bins merged.
   /// \endcode
   ///
   /// RSliceSpec::ROperationSum sums the bin contents along that axis, which allows to project to a lower-dimensional
   /// histogram:
   /// \code
   /// ROOT::Experimental::RHistEngine<int> hist({/* two dimensions */});
   /// // Fill the histogram with a number of entries...
   /// auto projected = hist.Slice(ROOT::Experimental::RSliceSpec{}, ROOT::Experimental::RSliceSpec::ROperationSum{});
   /// // The returned histogram has one dimension, with bin contents summed along the second axis.
   /// \endcode
   /// Note that it is not allowed to sum along all histogram axes because the return value would be a scalar.
   ///
   /// Ranges and operations can be combined. In that case, the range is applied before the operation.
   ///
   /// \param[in] args the arguments for each axis
   /// \return the sliced histogram
   /// \par See also
   /// the \ref Slice(const std::vector<RSliceSpec> &sliceSpecs) const "function overload" accepting `std::vector`
   template <typename... A>
   RHistEngine Slice(const A &...args) const
   {
      std::vector<RSliceSpec> sliceSpecs{args...};
      return Slice(sliceSpecs);
   }
 
   /// Create an atomic snapshot of this histogram engine.
   ///
   /// A snapshot is a consistent copy of the histogram, during concurrent filling. It is guaranteed that the returned
   /// copy represents a state between the begin and end of the snapshot operation.
   ///
   /// Snapshotting a histogram engine with many bins can be an expensive operation.
   ///
   /// \return the atomic snapshot
   RHistEngine SnapshotAtomic() const
   {
      static_assert(std::is_trivially_copyable_v<BinContentType>,
                    "snapshotting requires a trivially copyable bin content type");
 
      do {
         while (fSnapshotInProgress.load(std::memory_order_relaxed)) {
            // Spin while another snapshot is running
         }
      } while (fSnapshotInProgress.exchange(true, std::memory_order_relaxed));
 
      RHistEngine snapshot(fAxes.Get());
      // Do a first collect.
      for (std::size_t i = 0; i < fBinContents.size(); i++) {
         Internal::AtomicLoad(&fBinContents[i], &snapshot.fBinContents[i]);
      }
 
      // Now do another collect. If no change is detected, the snapshot is consistent. Otherwise update the bin contents
      // and try again.
      BinContentType tmp;
      bool changed;
      do {
         // To guarantee correctness, we let the release operation(s) in FillAtomic synchronize with this acquire fence.
         // This ensures that all previous writes become visible side-effects and the atomic loads will see them.
         std::atomic_thread_fence(std::memory_order_acquire);
 
         changed = false;
         for (std::size_t i = 0; i < fBinContents.size(); i++) {
            Internal::AtomicLoad(&fBinContents[i], &tmp);
            if (std::memcmp(&tmp, &snapshot.fBinContents[i], sizeof(BinContentType))) {
               std::memcpy(&snapshot.fBinContents[i], &tmp, sizeof(BinContentType));
               changed = true;
            }
         }
      } while (changed);
 
      fSnapshotInProgress.store(false, std::memory_order_relaxed);
 
      return snapshot;
   }
 
   /// \}
 
   /// %ROOT Streamer function to throw when trying to store an object of this class.
   void Streamer(TBuffer &) { throw std::runtime_error("unable to store RHistEngine"); }
};
 
} // namespace Experimental
} // namespace ROOT
 
#endif