THistData.hxx

Go to the documentation of this file.

/// \file ROOT/THistData.h
/// \ingroup Hist ROOT7
/// \author Axel Naumann <axel@cern.ch>
/// \date 2015-06-14
/// \warning This is part of the ROOT 7 prototype! It will change without notice. It might trigger earthquakes. Feedback is welcome!

/*************************************************************************
 * Copyright (C) 1995-2015, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

#ifndef ROOT7_THistData_h
#define ROOT7_THistData_h

#include <cmath>
#include <vector>
#include "ROOT/RArrayView.hxx"
#include "ROOT/THistUtils.hxx"

namespace ROOT {
namespace Experimental {

template<int DIMENSIONS, class PRECISION,
  template <int D_, class P_, template <class P__> class STORAGE> class... STAT>
class THist;

/**
 \class THistStatContent
 Basic histogram statistics, keeping track of the bin content and the total
 number of calls to Fill().
 */
template<int DIMENSIONS, class PRECISION,
  template <class PRECISION_> class STORAGE>
class THistStatContent {
public:
  /// The type of a (possibly multi-dimensional) coordinate.
  using CoordArray_t = Hist::CoordArray_t<DIMENSIONS>;
  /// The type of the weight and the bin content.
  using Weight_t = PRECISION;
  /// Type of the bin content array.
  using Content_t = STORAGE<PRECISION>;

  /**
   \class TConstBinStat
   Const view on a THistStatContent for a given bin.
  */
  class TConstBinStat {
  public:
    TConstBinStat(const THistStatContent& stat, int index):
      fContent(stat.GetBinContent(index)) {}
    PRECISION GetContent() const { return fContent; }

  private:
    PRECISION fContent; ///< The content of this bin.
  };

  /**
   \class TBinStat
   Modifying view on a THistStatContent for a given bin.
  */
  class TBinStat {
  public:
    TBinStat(THistStatContent& stat, int index):
      fContent(stat.GetBinContent(index)) {}
    PRECISION& GetContent() const { return fContent; }

  private:
    PRECISION& fContent; ///< The content of this bin.
  };

  using ConstBinStat_t = TConstBinStat;
  using BinStat_t =  TBinStat;

private:
  /// Number of calls to Fill().
  int64_t fEntries = 0;

  /// Bin content.
  Content_t fBinContent;

public:
  THistStatContent() = default;
  THistStatContent(size_t in_size): fBinContent(in_size) {}

  /// Add weight to the bin content at binidx.
  void Fill(const CoordArray_t& /*x*/, int binidx, Weight_t weight = 1.) {
    fBinContent[binidx] += weight;
    ++fEntries;
  }

  /// Get the number of entries filled into the histogram - i.e. the number of
  /// calls to Fill().
  int64_t GetEntries() const { return fEntries; }

  /// Get the number of bins.
  size_t size() const noexcept { return fBinContent.size(); }

  /// Get the bin content for the given bin.
  Weight_t operator[](int idx) const { return fBinContent[idx]; }
  /// Get the bin content for the given bin (non-const).
  Weight_t& operator[](int idx) { return fBinContent[idx]; }

  /// Get the bin content for the given bin.
  Weight_t GetBinContent(int idx) const { return fBinContent[idx]; }
  /// Get the bin content for the given bin (non-const).
  Weight_t& GetBinContent(int idx) { return fBinContent[idx]; }

  /// Retrieve the content array.
  const Content_t& GetContentArray() const { return fBinContent; }
  /// Retrieve the content array (non-const).
  Content_t& GetContentArray() { return fBinContent; }
};

/**
 \class THistStatTotalSumOfWeights
 Keeps track of the histogram's total sum of weights.
 */
template<int DIMENSIONS, class PRECISION, template <class P_> class STORAGE>
class THistStatTotalSumOfWeights {
public:
  /// The type of a (possibly multi-dimensional) coordinate.
  using CoordArray_t = Hist::CoordArray_t<DIMENSIONS>;
  /// The type of the weight and the bin content.
  using Weight_t = PRECISION;

  /**
   \class TBinStat
   No-op; this class does not provide per-bin statistics.
  */
  class TBinStat {
  public:
    TBinStat(const THistStatTotalSumOfWeights&, int) {}
  };

  using ConstBinStat_t = TBinStat;
  using BinStat_t =  TBinStat;

private:
  /// Sum of weights.
  PRECISION fSumWeights = 0;

public:
  THistStatTotalSumOfWeights() = default;
  THistStatTotalSumOfWeights(size_t) {}

  /// Add weight to the bin content at binidx.
  void Fill(const CoordArray_t& /*x*/, int, Weight_t weight = 1.) {
    fSumWeights += weight;
  }

  /// Get the sum of weights.
  Weight_t GetSumOfWeights() const { return fSumWeights; }

};


/**
 \class THistStatTotalSumOfSquaredWeights
 Keeps track of the histogram's total sum of squared weights.
 */
template<int DIMENSIONS, class PRECISION, template <class P_> class STORAGE>
class THistStatTotalSumOfSquaredWeights {
public:
  /// The type of a (possibly multi-dimensional) coordinate.
  using CoordArray_t = Hist::CoordArray_t<DIMENSIONS>;
  /// The type of the weight and the bin content.
  using Weight_t = PRECISION;

  /**
   \class TBinStat
   No-op; this class does not provide per-bin statistics.
  */
  class TBinStat {
  public:
    TBinStat(const THistStatTotalSumOfSquaredWeights&, int) {}
  };

  using ConstBinStat_t = TBinStat;
  using BinStat_t =  TBinStat;

private:
  /// Sum of (weights^2).
  PRECISION fSumWeights2 = 0;

public:
  THistStatTotalSumOfSquaredWeights() = default;
  THistStatTotalSumOfSquaredWeights(size_t) {}

  /// Add weight to the bin content at binidx.
  void Fill(const CoordArray_t& /*x*/, int /*binidx*/, Weight_t weight = 1.) {
    fSumWeights2 += weight * weight;
  }

  /// Get the sum of weights.
  Weight_t GetSumOfSquaredWeights() const { return fSumWeights2; }

};


/**
 \class THistStatUncertainty
 Histogram statistics to keep track of the Poisson uncertainty per bin.
 */
template<int DIMENSIONS, class PRECISION, template <class P_> class STORAGE>
class THistStatUncertainty {

public:
  /// The type of a (possibly multi-dimensional) coordinate.
  using CoordArray_t = Hist::CoordArray_t<DIMENSIONS>;
  /// The type of the weight and the bin content.
  using Weight_t = PRECISION;
  /// Type of the bin content array.
  using Content_t = STORAGE<PRECISION>;

  /**
   \class TConstBinStat
   Const view on a THistStatUncertainty for a given bin.
  */
  class TConstBinStat {
  public:
    TConstBinStat(const THistStatUncertainty& stat, int index):
      fSumW2(stat.GetSumOfSquaredWeights(index)) {}
    PRECISION GetSumW2() const { return fSumW2; }

    double GetUncertaintyImpl() const { return std::sqrt(std::abs(fSumW2)); }

  private:
    PRECISION fSumW2; ///< The bin's sum of square of weights.
  };

  /**
   \class TBinStat
   Modifying view on a THistStatUncertainty for a given bin.
  */
  class TBinStat {
  public:
    TBinStat(THistStatUncertainty& stat, int index):
      fSumW2(stat.GetSumOfSquaredWeights(index)) {}
    PRECISION& GetSumW2() const { return fSumW2; }
    // Can never modify this. Set GetSumW2() instead.
    double GetUncertaintyImpl() const { return std::sqrt(std::abs(fSumW2)); }

  private:
    PRECISION& fSumW2; ///< The bin's sum of square of weights.
  };

  using ConstBinStat_t = TConstBinStat;
  using BinStat_t =  TBinStat;

private:
  /// Uncertainty of the content for each bin.
  Content_t fSumWeightsSquared; ///< Sum of squared weights

public:
  THistStatUncertainty() = default;
  THistStatUncertainty(size_t size): fSumWeightsSquared(size) {}

  /// Add weight to the bin at binidx; the coordinate was x.
  void Fill(const CoordArray_t& /*x*/, int binidx, Weight_t weight = 1.) {
    fSumWeightsSquared[binidx] += weight * weight;
  }

  /// Calculate a bin's (Poisson) uncertainty of the bin content as the
  /// square-root of the bin's sum of squared weights.
  double GetBinUncertaintyImpl(int binidx) const {
    return std::sqrt(fSumWeightsSquared[binidx]);
  }

  /// Get a bin's sum of squared weights.
  Weight_t GetSumOfSquaredWeights(int binidx) const {
    return fSumWeightsSquared[binidx];
  }

  /// Get a bin's sum of squared weights.
  Weight_t& GetSumOfSquaredWeights(int binidx) {
    return fSumWeightsSquared[binidx];
  }

  /// Get the structure holding the sum of squares of weights.
  const std::vector<double>& GetSumOfSquaredWeights() const { return fSumWeightsSquared; }
  /// Get the structure holding the sum of squares of weights (non-const).
  std::vector<double>& GetSumOfSquaredWeights() { return fSumWeightsSquared; }
};

/** \class THistDataMomentUncert
  For now do as TH1: calculate first (xw) and second (x^2w) moment.
*/
template<int DIMENSIONS, class PRECISION, template <class P_> class STORAGE>
class THistDataMomentUncert {
public:
  /// The type of a (possibly multi-dimensional) coordinate.
  using CoordArray_t = Hist::CoordArray_t<DIMENSIONS>;
  /// The type of the weight and the bin content.
  using Weight_t = PRECISION;
  /// Type of the bin content array.
  using Content_t = STORAGE<PRECISION>;

  /**
   \class TBinStat
   No-op; this class does not provide per-bin statistics.
  */
  class TBinStat {
  public:
    TBinStat(const THistDataMomentUncert&, int) {}
  };

  using ConstBinStat_t = TBinStat;
  using BinStat_t =  TBinStat;

private:
  std::array<Weight_t, DIMENSIONS> fMomentXW;
  std::array<Weight_t, DIMENSIONS> fMomentX2W;

public:
  THistDataMomentUncert() = default;
  THistDataMomentUncert(size_t) {}

  /// Add weight to the bin at binidx; the coordinate was x.
  void Fill(const CoordArray_t &x, int /*binidx*/, Weight_t weight = 1.) {
    for (int idim = 0; idim < DIMENSIONS; ++idim) {
      const PRECISION xw = x[idim] * weight;
      fMomentXW[idim] += xw ;
      fMomentX2W[idim] += x[idim] * xw;
    }
  }
};


/** \class THistStatRuntime
  Interface implementing a pure virtual functions DoFill(), DoFillN().
  */
template<int DIMENSIONS, class PRECISION, template <class P_> class STORAGE>
class THistStatRuntime {
public:
  /// The type of a (possibly multi-dimensional) coordinate.
  using CoordArray_t = Hist::CoordArray_t<DIMENSIONS>;
  /// The type of the weight and the bin content.
  using Weight_t = PRECISION;
  /// Type of the bin content array.
  using Content_t = STORAGE<PRECISION>;

  /**
   \class TBinStat
   No-op; this class does not provide per-bin statistics.
  */
  class TBinStat {
  public:
    TBinStat(const THistStatRuntime&, int) {}
  };

  using ConstBinStat_t = TBinStat;
  using BinStat_t =  TBinStat;

  THistStatRuntime() = default;
  THistStatRuntime(size_t) {}
  virtual ~THistStatRuntime() = default;


  virtual void DoFill(const CoordArray_t &x, int binidx, Weight_t weightN) = 0;
  void Fill(const CoordArray_t &x, int binidx, Weight_t weight = 1.) {
    DoFill(x, binidx, weight);
  }
};


namespace Detail {

/// std::vector has more template arguments; for the default storage we don't
/// care about them, so use-decl them away:
template<class PRECISION>
using THistDataDefaultStorage = std::vector<PRECISION>;


/** \class THistBinStat
  Const view on a bin's statistical data. Combines all STATs' BinStat_t views.
  */
template <class DATA, class...BASES>
class THistBinStat: public BASES... {
private:
  /// Check whether `double T::GetBinUncertaintyImpl(int)` can be called.
  template <class T>
  static auto HaveUncertainty(const T* This) -> decltype(This->GetUncertaintyImpl());
  /// Fall-back case for check whether `double T::GetBinUncertaintyImpl(int)` can be called.
  template <class T>
  static char HaveUncertainty(...);

public:
  THistBinStat(DATA& data, int index):
    BASES(data, index)... {}

  /// Whether this provides storage for uncertainties, or whether uncertainties
  /// are determined as poisson uncertainty of the content.
  static constexpr bool HasBinUncertainty() {
    struct AllYourBaseAreBelongToUs: public BASES... {};
    return sizeof(HaveUncertainty<AllYourBaseAreBelongToUs>(nullptr)) == sizeof(double);
  }
  /// Calculate the bin content's uncertainty for the given bin, using base class information,
  /// i.e. forwarding to a base's `GetUncertaintyImpl()`.
  template <bool B = true,
    class = typename std::enable_if<B && HasBinUncertainty()>::type>
  double GetUncertainty() const {
    return this->GetUncertaintyImpl();
  }
  /// Calculate the bin content's uncertainty for the given bin, using Poisson
  /// statistics on the absolute bin content. Only available if no base provides
  /// this functionality. Requires GetContent().
  template <bool B = true,
    class = typename std::enable_if<B && !HasBinUncertainty()>::type>
  double GetUncertainty(...) const {
    auto content = this->GetContent();
    return std::sqrt(std::fabs(content));
  }
};


/** \class THistData
  A THistImplBase's data, provides accessors to all its statistics.
  */
template<int DIMENSIONS, class PRECISION,
  template <class P_> class STORAGE,
  template <int D_, class P_, template <class P__> class S_> class... STAT>
class THistData: public STAT<DIMENSIONS, PRECISION, STORAGE>... {
private:
  /// Check whether `double T::GetBinUncertaintyImpl(int)` can be called.
  template <class T>
  static auto HaveUncertainty(const T* This) -> decltype(This->GetBinUncertaintyImpl(12));
  /// Fall-back case for check whether `double T::GetBinUncertaintyImpl(int)` can be called.
  template <class T>
  static char HaveUncertainty(...);

public:
  /// Matching THist
  using Hist_t = THist<DIMENSIONS, PRECISION, STAT...>;

  /// The type of the weight and the bin content.
  using Weight_t = PRECISION;

  /// The type of a (possibly multi-dimensional) coordinate.
  using CoordArray_t = Hist::CoordArray_t<DIMENSIONS>;

  /// The type of a non-modifying view on a bin.
  using ConstHistBinStat_t
    = THistBinStat<const THistData, typename STAT<DIMENSIONS, PRECISION, STORAGE>::ConstBinStat_t...>;

  /// The type of a modifying view on a bin.
  using HistBinStat_t
    = THistBinStat<THistData, typename STAT<DIMENSIONS, PRECISION, STORAGE>::BinStat_t...>;

  /// Number of dimensions of the coordinates
  static constexpr int GetNDim() noexcept { return DIMENSIONS; }

  THistData() = default;


  /// Constructor providing the number of bins (incl under, overflow) to the
  /// base classes.
  THistData(size_t size): STAT<DIMENSIONS, PRECISION, STORAGE>(size)... {}

  /// Fill weight at x to the bin content at binidx.
  void Fill(const CoordArray_t& x, int binidx, Weight_t weight = 1.) {
    // Call Fill() on all base classes.
    // This combines a couple of C++ spells:
    // - "STAT": is a template parameter pack of template template arguments. It
    //           has multiple (or one or no) elements; each is a template name
    //           that needs to be instantiated before it can be used.
    // - "...":  template parameter pack expansion; the expression is evaluated
    //           for each STAT. The expression is
    //           (STAT<DIMENSIONS, PRECISION, STORAGE>::Fill(x, binidx, weight), 0)
    // - "trigger_base_fill{}":
    //           initialization, provides a context in which template parameter
    //           pack expansion happens.
    // - ", 0":  because Fill() returns void it cannot be used as initializer
    //           expression. The trailing ", 0" gives it the type of the trailing
    //           comma-separated expression - int.
    using trigger_base_fill = int[];
    (void)trigger_base_fill{ (STAT<DIMENSIONS, PRECISION, STORAGE>::Fill(x, binidx, weight), 0)... };
  }

  /// Whether this provides storage for uncertainties, or whether uncertainties
  /// are determined as poisson uncertainty of the content.
  static constexpr bool HasBinUncertainty() {
    struct AllYourBaseAreBelongToUs: public STAT<DIMENSIONS, PRECISION, STORAGE>... {};
    return sizeof(HaveUncertainty<AllYourBaseAreBelongToUs>(nullptr)) == sizeof(double);
  }

  /// Calculate the bin content's uncertainty for the given bin, using base class information,
  /// i.e. forwarding to a base's `GetBinUncertaintyImpl(binidx)`.
  template <bool B = true, class = typename std::enable_if<B && HasBinUncertainty()>::type>
  double GetBinUncertainty(int binidx) const {
    return this->GetBinUncertaintyImpl(binidx);
  }
  /// Calculate the bin content's uncertainty for the given bin, using Poisson
  /// statistics on the absolute bin content. Only available if no base provides
  /// this functionality. Requires GetContent().
  template <bool B = true, class = typename std::enable_if<B && !HasBinUncertainty()>::type>
  double GetBinUncertainty(int binidx, ...) const {
    auto content = this->GetBinContent(binidx);
    return std::sqrt(std::fabs(content));
  }

  /// Get a view on the statistics values of a bin.
  ConstHistBinStat_t GetView(int idx) const {
    return ConstHistBinStat_t(*this, idx);
  }
  /// Get a (non-const) view on the statistics values of a bin.
  HistBinStat_t GetView(int idx) {
    return HistBinStat_t(*this, idx);
  }
};
} // namespace Detail

} // namespace Experimental
} // namespace ROOT
#endif