Histogram.h

Go to the documentation of this file.

/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Sergei Zharko [committer] */
 
 
#pragma once
 
#include "Accumulators.h"
 
#include <boost/histogram.hpp>
//#include <boost/histogram/algorithm/sum.hpp>
#include <boost/histogram/serialization.hpp>
#include <boost/serialization/access.hpp>
#include <boost/serialization/base_object.hpp>
#include <boost/serialization/string.hpp>
#include <boost/serialization/vector.hpp>
 
#include <algorithm>
#include <cstdint>
#include <memory>
#include <numeric>
#include <string>
#include <tuple>
#include <type_traits>
 
#include <fmt/format.h>
 
namespace cbm::algo::qa
{
  namespace bh = boost::histogram;
 
  using RegularAxis_t = bh::axis::regular<>;
 
  using Axes1D_t = std::tuple<RegularAxis_t>;
  using Axes2D_t = std::tuple<RegularAxis_t, RegularAxis_t>;
 
  using HistStorage_t = bh::weight_storage;
  using ProfStorage_t = bh::RootStyleProfileStorage<double>;
 
  enum class EAxis : unsigned
  {
    X = 0,
    Y,
    Z
  };
 
  enum class EHistFlag : uint8_t
  {
    StoreVsTsId    = 0b00000001,  
    OmitIntegrated = 0b00000010,  
    SetMinimum     = 0b00000100   
  };
 
  class HistogramMetadata {
   public:
    using Flags_t = std::underlying_type_t<EHistFlag>;
 
    static constexpr std::string_view ksTsIdSuffix = "_ts_id";  
 
    HistogramMetadata() = default;
 
    explicit HistogramMetadata(const std::string& msg)
    {
      if (!msg.empty()) {
        fFlags = std::stoi(msg, nullptr, 16);
      }
    }
 
    ~HistogramMetadata() = default;
 
    bool CheckFlags() const
    {
      // The histogram must be plotted either vs TS, or over all TS
      return (GetFlag(EHistFlag::StoreVsTsId) || !GetFlag(EHistFlag::OmitIntegrated));
    }
 
    bool GetFlag(EHistFlag key) const { return static_cast<bool>(fFlags & static_cast<Flags_t>(key)); }
 
    void SetFlag(EHistFlag key, bool flag = true)
    {
      flag ? (fFlags |= static_cast<Flags_t>(key)) : (fFlags &= ~static_cast<Flags_t>(key));
    }
 
    std::string ToString() const { return fmt::format("{0:02x}", fFlags); }
 
    static std::pair<std::string, std::string> SeparateNameAndMetadata(const std::string& msg)
    {
      size_t pos = msg.find_last_of('!');
      if (pos != msg.npos) {
        return std::make_pair(msg.substr(0, pos), msg.substr(pos + 1));
      }
      else {
        return std::make_pair(msg, "");
      }
    }
 
   private:
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive& ar, const unsigned int /*version*/)
    {
      ar& fFlags;
    }
 
    Flags_t fFlags{};  
  };
 
  class TotalSums1D {
   public:
    double GetTotSumW() const { return fTotSumW; }
 
    double GetTotSumW2() const { return fTotSumW2; }
 
    double GetTotSumWX() const { return fTotSumWX; }
 
    double GetTotSumWX2() const { return fTotSumWX2; }
 
   protected:
    void UpdateTotalSums(double x, double w)
    {
      fTotSumW += w;
      fTotSumW2 += w * w;
      fTotSumWX += w * x;
      fTotSumWX2 += w * x * x;
    }
 
    void Reset()
    {
      fTotSumW   = 0;
      fTotSumW2  = 0;
      fTotSumWX  = 0;
      fTotSumWX2 = 0;
    }
 
    double fTotSumW   = 0.;  
    double fTotSumW2  = 0.;  
    double fTotSumWX  = 0.;  
    double fTotSumWX2 = 0.;  
 
   private:
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive& ar, const unsigned int /*version*/)
    {
      ar& fTotSumW;
      ar& fTotSumW2;
      ar& fTotSumWX;
      ar& fTotSumWX2;
    }
  };
 
  class TotalSums2D : public TotalSums1D {
   public:
    double GetTotSumWXY() const { return fTotSumWXY; }
 
    double GetTotSumWY() const { return fTotSumWY; }
 
    double GetTotSumWY2() const { return fTotSumWY2; }
 
   protected:
    void Reset()
    {
      TotalSums1D::Reset();
      fTotSumWXY = 0;
      fTotSumWY  = 0;
      fTotSumWY2 = 0;
    }
 
    void UpdateTotalSums(double x, double y, double w)
    {
      TotalSums1D::UpdateTotalSums(x, w);
      fTotSumWXY += w * x * y;
      fTotSumWY += w * y;
      fTotSumWY2 += w * y * y;
    }
 
    double fTotSumWY  = 0.;  
    double fTotSumWXY = 0.;  
    double fTotSumWY2 = 0.;  
 
   private:
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive& ar, const unsigned int /*version*/)
    {
      ar& boost::serialization::base_object<TotalSums1D>(*this);
      ar& fTotSumWXY;
      ar& fTotSumWY;
      ar& fTotSumWY2;
    }
  };
 
 
  template<class Axes, class Storage, class TotalSums>
  class Histogram : public TotalSums {
   protected:
    using Hist_t                   = bh::histogram<Axes, Storage>;  // TODO: Test in future
    static constexpr unsigned Rank = std::tuple_size_v<Axes>;
 
   public:
    ~Histogram() = default;
 
    template<EAxis A, unsigned IA = static_cast<unsigned>(A), std::enable_if_t<(IA < Rank), bool> = true>
    double GetAxisMin() const
    {
      return fHistogram.template axis<IA>().value(0.);
    }
 
    template<EAxis A, unsigned IA = static_cast<unsigned>(A), std::enable_if_t<(IA < Rank), bool> = true>
    double GetAxisMax() const
    {
      return fHistogram.template axis<IA>().value(static_cast<double>(GetAxisNbins<A>()));
    }
 
    template<EAxis A, unsigned IA = static_cast<unsigned>(A), std::enable_if_t<(IA < Rank), bool> = true>
    uint32_t GetAxisNbins() const
    {
      return fHistogram.template axis<IA>().size();
    }
 
    // TODO: Gives different results, if weights are not 1 (investigate!)
    double GetEntries() const
    {
      // return bh::algorithm::sum(fHistogram); // -> effective entries (but different to the ROOT ones: if w != 1)
      return fEntries;
    }
 
    bool GetFlag(EHistFlag key) const { return fMetadata.GetFlag(key); }
 
    const std::string& GetName() const { return fName; }
 
    const HistogramMetadata& GetMetadata() const { return fMetadata; }
 
    std::string GetMetadataString() const { return fMetadata.ToString(); }
 
    uint32_t GetNbinsX() const { return GetAxisNbins<EAxis::X>(); }
 
    double GetMinX() const { return GetAxisMin<EAxis::X>(); }
 
    double GetMaxX() const { return GetAxisMax<EAxis::X>(); }
 
    template<unsigned RankCheck = Rank, std::enable_if_t<(RankCheck > 1), bool> = true>
    uint32_t GetNbinsY() const
    {
      return GetAxisNbins<EAxis::Y>();
    }
 
    template<unsigned RankCheck = Rank, std::enable_if_t<(RankCheck > 1), bool> = true>
    double GetMinY() const
    {
      return GetAxisMin<EAxis::Y>();
    }
 
    template<unsigned RankCheck = Rank, std::enable_if_t<(RankCheck > 1), bool> = true>
    double GetMaxY() const
    {
      return GetAxisMax<EAxis::Y>();
    }
 
    double GetMaximum() const
    {
      return *std::max_element(bh::indexed(fHistogram).begin(), bh::indexed(fHistogram).end());
    }
 
    double GetMinimum() const
    {
      return *std::min_element(bh::indexed(fHistogram).begin(), bh::indexed(fHistogram).end());
    }
 
    const std::string& GetTitle() const { return fTitle; }
 
    void Reset()
    {
      TotalSums::Reset();
      fHistogram.reset();
      fEntries = 0;
    }
 
    void SetFlag(EHistFlag key, bool flag = true) { fMetadata.SetFlag(key, flag); }
 
    void SetName(const std::string& name) { fName = name; }
 
    void SetTitle(const std::string& title) { fTitle = title; }
 
    //std::string ToString() const
    //{
    //  std::stringstream msg;
    //  msg <<
    //}
 
   protected:
    Histogram() = default;
 
    explicit Histogram(const Histogram<Axes, Storage, TotalSums>& h) = default;
 
    //explicit Histogram(const Histogram<Axes, Storage>& h)
    //  : fHistogram(h.fHistogram)
    //  , fName(h.fName)
    //  , fTitle(h.fTitle)
    //  , fEntries(h.fEntries)
    //{
    //}
 
    Histogram(const Hist_t& bhist, const std::string& name, const std::string title)
      : fHistogram(bhist)
      , fName(name)
      , fTitle(title)
      , fEntries(0)
      , fMetadata(HistogramMetadata())
    {
    }
 
    inline static int GetBinBH(uint32_t iBin) { return (iBin > 0) ? iBin - 1 : -1; }
 
    Hist_t fHistogram;            
    std::string fName  = "";      
    std::string fTitle = "";      
    int fEntries       = 0;       
    HistogramMetadata fMetadata;  
 
 
   private:
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive& ar, const unsigned int /*version*/)
    {
      ar& boost::serialization::base_object<TotalSums>(*this);
      ar& fHistogram;
      ar& fName;
      ar& fTitle;
      ar& fEntries;
      ar& fMetadata;
    }
  };
 
  using BaseH1D    = Histogram<Axes1D_t, HistStorage_t, TotalSums1D>;
  using BaseH2D    = Histogram<Axes2D_t, HistStorage_t, TotalSums2D>;
  using BaseProf1D = Histogram<Axes1D_t, ProfStorage_t, TotalSums1D>;
  using BaseProf2D = Histogram<Axes2D_t, ProfStorage_t, TotalSums2D>;
 
  class H1D : public BaseH1D {
   public:
    H1D(const std::string& name, const std::string& title, uint32_t nBins, double xMin, double xMax)
      : BaseH1D(bh::make_weighted_histogram(RegularAxis_t(nBins, xMin, xMax)), name, title)
    {
    }
 
    H1D() = default;
 
    H1D(const H1D&) = default;
 
    H1D(H1D&&) = default;
 
    H1D& operator=(const H1D&) = default;
 
    H1D& operator=(H1D&&) = default;
 
    int Fill(double x, double w = 1.)
    {
      auto cellIt = fHistogram(bh::weight(w), x);
      ++fEntries;
      int iBin = cellIt - fHistogram.begin();
      if (iBin == 0 || iBin > static_cast<int>(GetNbinsX())) {
        return -1;  // ROOT TH1::Fill behaviour
      }
      // NOTE: In ROOT TH1 the total sums are updated only for non-overflow bins
      BaseH1D::UpdateTotalSums(x, w);
      return iBin;
    }
 
    auto GetBinAccumulator(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)); }
 
    double GetBinContent(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)).value(); }
 
    double GetBinError(uint32_t iBin) const { return std::sqrt(fHistogram.at(Histogram::GetBinBH(iBin)).variance()); }
 
   private:
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive& ar, const unsigned int /*version*/)
    {
      ar& boost::serialization::base_object<BaseH1D>(*this);
    }
  };
 
  class H2D : public BaseH2D {
   public:
    H2D(const std::string& name, const std::string& title, uint32_t nBinsX, double xMin, double xMax, uint32_t nBinsY,
        double yMin, double yMax)
      : BaseH2D(bh::make_weighted_histogram(RegularAxis_t(nBinsX, xMin, xMax), RegularAxis_t(nBinsY, yMin, yMax)), name,
                title)
    {
    }
 
    H2D() = default;
 
    H2D(const H2D&) = default;
 
    H2D(H2D&&) = default;
 
    H2D& operator=(const H2D&) = default;
 
    H2D& operator=(H2D&&) = default;
 
    int Fill(double x, double y, double w = 1.)
    {
      auto cellIt = fHistogram(x, y, bh::weight(w));
      ++fEntries;
 
      int iBin       = cellIt - fHistogram.begin();
      uint32_t iBinX = iBin % (GetNbinsX() + 2);
      if (iBinX == 0 || iBinX > GetNbinsX()) {
        return -1;
      }
      uint32_t iBinY = iBin / (GetNbinsX() + 2);
      if (iBinY == 0 || iBinY > GetNbinsY()) {
        return -1;
      }
      // NOTE: In ROOT TH2 the total sums are updated only for non-overflow bins
      BaseH2D::UpdateTotalSums(x, y, w);
      return iBin;
    }
 
    auto GetBinAccumulator(uint32_t iBinX, uint32_t iBinY) const
    {
      return fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY));
    }
 
    double GetBinContent(uint32_t iBinX, uint32_t iBinY) const
    {
      return fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY)).value();
    }
 
    double GetBinError(uint32_t iBinX, uint32_t iBinY) const
    {
      return std::sqrt(fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY)).variance());
    }
 
   private:
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive& ar, const unsigned int /*version*/)
    {
      ar& boost::serialization::base_object<BaseH2D>(*this);
    }
  };
 
  // TODO: Boost::histogram by default uses a different error calculation approach as ROOT TProfile. TODO: investigate
 
  class Prof1D : public BaseProf1D {
   public:
    Prof1D(const std::string& name, const std::string& title, uint32_t nBinsX, double xMin, double xMax,
           double yMin = 0., double yMax = 0.)
      : BaseProf1D(bh::MakeRootStyleProfile(RegularAxis_t(nBinsX, xMin, xMax)), name, title)
      , fYmin(yMin)
      , fYmax(yMax)
    {
    }
 
    Prof1D() = default;
 
    Prof1D(const Prof1D&) = default;
 
    Prof1D(Prof1D&&) = default;
 
    Prof1D& operator=(const Prof1D&) = default;
 
    Prof1D& operator=(Prof1D&&) = default;
 
    int Fill(double x, double y, double w = 1.)
    {
      if ((fYmin != fYmax) && (y < fYmin || y > fYmax || std::isnan(y))) {
        return -1;
      }
 
      auto cellIt = fHistogram(x, bh::sample(y), bh::weight(w));
      ++fEntries;
      int iBin = cellIt - fHistogram.begin();
      if (iBin == 0 || iBin > static_cast<int>(GetNbinsX())) {
        return -1;  // ROOT TH1::Fill behaviour
      }
      // NOTE: In ROOT TProfile the total sums are updated only for non-overflow bins
      BaseProf1D::UpdateTotalSums(x, w);
      return iBin;
    }
 
    auto GetBinAccumulator(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)); }
 
    double GetBinContent(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)).GetMean(); }
 
    double GetBinCount(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)).GetEffCount(); }
 
    double GetBinError(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)).GetSEM(); }
 
    double GetMinY() const { return fYmin; }
 
    double GetMaxY() const { return fYmax; }
 
   private:
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive& ar, const unsigned int /*version*/)
    {
      ar& boost::serialization::base_object<BaseProf1D>(*this);
      ar& fYmin;
      ar& fYmax;
    }
 
    double fYmin = 0.;  
    double fYmax = 0.;  
  };
 
  class Prof2D : public BaseProf2D {
   public:
    Prof2D(const std::string& name, const std::string& title, uint32_t nBinsX, double xMin, double xMax,
           uint32_t nBinsY, double yMin, double yMax, double zMin = 0., double zMax = 0.)
      : BaseProf2D(bh::MakeRootStyleProfile(RegularAxis_t(nBinsX, xMin, xMax), RegularAxis_t(nBinsY, yMin, yMax)), name,
                   title)
      , fZmin(zMin)
      , fZmax(zMax)
    {
    }
 
    Prof2D() = default;
 
    Prof2D(const Prof2D&) = default;
 
    Prof2D(Prof2D&&) = default;
 
    Prof2D& operator=(const Prof2D&) = default;
 
    Prof2D& operator=(Prof2D&&) = default;
 
    int Fill(double x, double y, double z, double w = 1.)
    {
      if ((fZmin != fZmax) && (z < fZmin || z > fZmax || std::isnan(z))) {
        return -1;
      }
 
      auto cellIt = fHistogram(x, y, bh::sample(z), bh::weight(w));
      ++fEntries;
 
      int iBin       = cellIt - fHistogram.begin();
      uint32_t iBinX = iBin % (GetNbinsX() + 2);
      if (iBinX == 0 || iBinX > GetNbinsX()) {
        return -1;
      }
      uint32_t iBinY = iBin / (GetNbinsX() + 2);
      if (iBinY == 0 || iBinY > GetNbinsY()) {
        return -1;
      }
      // NOTE: In ROOT TProfile2D the total sums are updated only for non-overflow bins
      BaseProf2D::UpdateTotalSums(x, y, w);
      return iBin;
    }
 
    auto GetBinAccumulator(uint32_t iBinX, uint32_t iBinY) const
    {
      return fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY));
    }
 
    double GetBinContent(uint32_t iBinX, uint32_t iBinY) const
    {
      return fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY)).GetMean();
    }
 
    double GetBinCount(uint32_t iBinX, uint32_t iBinY) const
    {
      return fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY)).GetEffCount();
    }
 
    double GetBinError(uint32_t iBinX, uint32_t iBinY) const
    {
      return fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY)).GetSEM();
    }
 
    double GetMinZ() const { return fZmin; }
 
    double GetMaxZ() const { return fZmax; }
 
   private:
    friend class boost::serialization::access;
    template<class Archive>
    void serialize(Archive& ar, const unsigned int /*version*/)
    {
      ar& boost::serialization::base_object<BaseProf2D>(*this);
      ar& fZmin;
      ar& fZmax;
    }
 
    double fZmin = 0.;  
    double fZmax = 0.;  
  };
 
}  // namespace cbm::algo::qa