diff --git a/algo/qa/Accumulators.h b/algo/qa/Accumulators.h
new file mode 100644
index 0000000000000000000000000000000000000000..67867da4ca2bff9f1d182482f10b4b7a64743711
--- /dev/null
+++ b/algo/qa/Accumulators.h
@@ -0,0 +1,172 @@
+/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergei Zharko [committer] */
+
+/// \file Accumulators.h
+/// \date 03.03.2024
+/// \brief Custom accumulators for boost::histogram (header)
+/// \author Sergei Zharko <s.zharko@gsi.de>
+
+#pragma once
+
+#include <boost/histogram/fwd.hpp> // for weighted_mean<>
+#include <boost/histogram/make_histogram.hpp>
+#include <boost/serialization/access.hpp>
+
+#include <cmath>
+#include <type_traits>
+
+namespace boost::histogram
+{
+ /// \class RootStyleProfileAccumulator
+ /// \brief A ROOT-style accumulator for the Boost-histogram profiles
+ /// \tparam ValueType Type of the values
+ template<class ValueType>
+ class RootStyleProfileAccumulator {
+ public:
+ using Value_t = ValueType;
+ using ConstRef_t = const Value_t&;
+ using Weight_t = weight_type<Value_t>;
+
+ /// \brief Default constructor
+ RootStyleProfileAccumulator() = default;
+
+ /// \brief Copy constructor from the instances with other types
+ template<class OtherValueType>
+ RootStyleProfileAccumulator(const RootStyleProfileAccumulator<OtherValueType>& other)
+ : fSumWV(other.fSumWV)
+ , fSumWV2(other.fSumWV2)
+ , fSumW(other.fSumW)
+ , fSumW2(other.fSumW2)
+ {
+ }
+
+ /// \brief Constructor from the external counters
+ /// \param sumWV Sum of weight * value
+ /// \param sumWV Sum of weight * value * value
+ /// \param sumW Sum of weights
+ /// \param sumW2 Sum of weights
+ RootStyleProfileAccumulator(ConstRef_t sumWV, ConstRef_t sumWV2, ConstRef_t sumW, ConstRef_t sumW2)
+ : fSumWV(sumWV)
+ , fSumWV2(sumWV2)
+ , fSumW(sumW)
+ , fSumW2(sumW2)
+ {
+ }
+
+ /// \brief Inserts sample with 1. weight
+ /// \param v Sample
+ void operator()(ConstRef_t v) { operator()(weight(1.), v); }
+
+ /// \brief Inserts sample with weight
+ /// \param v Sample
+ /// \param w Sample weight
+ void operator()(const Weight_t& w, ConstRef_t v)
+ {
+ fSumWV += w.value * v;
+ fSumWV2 += w.value * v * v;
+ fSumW += w.value;
+ fSumW2 += w.value * w.value;
+ }
+
+ /// \brief Adds another accumulator
+ RootStyleProfileAccumulator& operator+=(const RootStyleProfileAccumulator& rhs)
+ {
+ if (rhs.fSumW == 0) {
+ return *this;
+ }
+ fSumWV += rhs.fSumWV;
+ fSumWV2 += rhs.fSumWV2;
+ fSumW += rhs.fSumW;
+ fSumW2 += rhs.fSumW2;
+ return *this;
+ }
+
+ /// \brief Scales all entries by value
+ RootStyleProfileAccumulator& operator*=(ConstRef_t s) noexcept
+ {
+ fSumWV *= s;
+ fSumWV2 *= s * s;
+ return *this;
+ }
+
+ /// \brief Equality operator
+ bool operator==(const RootStyleProfileAccumulator& rhs) const noexcept
+ {
+ return fSumWV == rhs.fSumWV && fSumWV2 == rhs.fSumWV2 && fSumW == rhs.fSumW && fSumW2 == rhs.fSumW2;
+ }
+
+ /// \brief Non-equality operator
+ bool operator!=(const RootStyleProfileAccumulator& rhs) const noexcept { return !operator==(rhs); }
+
+ /// \brief Returns sum of weights
+ ConstRef_t GetSumW() const noexcept { return fSumW; }
+
+ /// \brief Returns sum of weight squares
+ ConstRef_t GetSumW2() const noexcept { return fSumW2; }
+
+ /// \brief Returns sum of products of weight multiplied by value
+ ConstRef_t GetSumWV() const noexcept { return fSumWV; }
+
+ /// \brief Returns sum of products of weight multiplied by squared value
+ ConstRef_t GetSumWV2() const noexcept { return fSumWV2; }
+
+ /// \brief Returns the mean value
+ Value_t GetMean() const noexcept { return fSumW ? (fSumWV / fSumW) : 0; }
+
+ /// \brief Returns effective sum of values
+ Value_t GetEffCount() const noexcept { return fSumW2 ? (fSumW * fSumW / fSumW2) : 0; }
+
+ /// \brief Returns variance
+ Value_t GetVariance() const noexcept { return fSumW ? (fSumWV2 / fSumW - GetMean() * GetMean()) : 0; }
+
+ /// \brief Returns standard deviation of the sample
+ Value_t GetStdDev() const noexcept { return std::sqrt(GetVariance()); }
+
+ /// \brief Returns standard error of the mean
+ Value_t GetSEM() const noexcept { return fSumW ? GetStdDev() / std::sqrt(GetEffCount()) : 0; }
+
+ private:
+ friend class boost::serialization::access;
+ /// \brief Serialization rule
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int /*version*/)
+ {
+ ar& make_nvp("fSumWV", fSumWV);
+ ar& make_nvp("fSumWV2", fSumWV2);
+ ar& make_nvp("fSumW", fSumW);
+ ar& make_nvp("fSumW2", fSumW2);
+ }
+
+ Value_t fSumWV{}; ///< Sum of weight * value
+ Value_t fSumWV2{}; ///< Sum of weight * value * value
+ Value_t fSumW{}; ///< Sum of weights
+ Value_t fSumW2{}; ///< Sum of squared weights
+ };
+} // namespace boost::histogram
+
+#ifndef BOOST_HISTOGRAM_DOXYGEN_INVOKED
+namespace std
+{
+ template<class T, class U>
+ /// Specialization for boost::histogram::accumulators::weighted_mean.
+ struct common_type<boost::histogram::RootStyleProfileAccumulator<T>,
+ boost::histogram::RootStyleProfileAccumulator<U>> {
+ using type = boost::histogram::RootStyleProfileAccumulator<common_type_t<T, U>>;
+ };
+} // namespace std
+#endif
+
+namespace boost::histogram
+{
+ template<typename T>
+ using RootStyleProfileStorage = dense_storage<RootStyleProfileAccumulator<T>>;
+
+ /// \brief Maker for RootStyleProfileAccumulator
+ template<class Axis, class... Axes, class = detail::requires_axis<Axis>>
+ auto MakeRootStyleProfile(Axis&& axis, Axes&&... axes)
+ {
+ return make_histogram_with(RootStyleProfileStorage<double>(), std::forward<Axis>(axis),
+ std::forward<Axes>(axes)...);
+ }
+} // namespace boost::histogram
diff --git a/algo/qa/Histogram.h b/algo/qa/Histogram.h
index adb2f5b80af74585359972a0691031bf3c1f2ab3..f2a29770e87c8ef58d90d9ac58c75958f0005d10 100644
--- a/algo/qa/Histogram.h
+++ b/algo/qa/Histogram.h
@@ -9,6 +9,8 @@
#pragma once
+#include "Accumulators.h"
+
#include <boost/histogram.hpp>
//#include <boost/histogram/algorithm/sum.hpp>
#include <boost/histogram/serialization.hpp>
@@ -17,6 +19,7 @@
#include <boost/serialization/string.hpp>
#include <boost/serialization/vector.hpp>
+#include <algorithm>
#include <cstdint>
#include <memory>
#include <numeric>
@@ -30,6 +33,13 @@ namespace cbm::algo::qa
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::default_storage;
+ using ProfStorage_t = bh::RootStyleProfileStorage<double>;
+
+
/// \enum EHistAxis
/// \brief Enumeration for histogram axes
enum class EAxis : unsigned
@@ -40,13 +50,13 @@ namespace cbm::algo::qa
};
/// \class Histogram
- /// \brief Interface to one-dimentional histogram
- /// \tparam Rank Rank of the histogram (1 - 1D, 2 - 2D, ...)
- /// \tparam Axes The axes type
- template<unsigned Rank, class Axes, std::enable_if_t<(Rank > 0) && (Rank < 3), bool> = true>
+ /// \brief Interface to a histogram/profile
+ /// \tparam Axes The axes type
+ /// \tparam Storage Storage type for the histogram
+ template<class Axes, class Storage>
class Histogram {
- //using Hist_t = bh::histogram<Axes, bh::weight_storage>; // TODO: Test in future
- using Hist_t = bh::histogram<Axes, bh::default_storage>;
+ using Hist_t = bh::histogram<Axes, Storage>; // TODO: Test in future
+ static constexpr unsigned Rank = std::tuple_size_v<Axes>;
public:
/// \brief Destructor
@@ -118,6 +128,18 @@ namespace cbm::algo::qa
return GetAxisMax<EAxis::Y>();
}
+ /// \brief Gets maximum value
+ double GetMaximum() const
+ {
+ return *std::max_element(bh::indexed(fHistogram).begin(), bh::indexed(fHistogram).end());
+ }
+
+ /// \brief Gets minimum value
+ double GetMinimum() const
+ {
+ return *std::min_element(bh::indexed(fHistogram).begin(), bh::indexed(fHistogram).end());
+ }
+
/// \brief Gets title
const std::string& GetTitle() const { return fTitle; }
@@ -133,9 +155,17 @@ namespace cbm::algo::qa
/// \note ROOT axes convention is applicable: "<main title>;<x-axis title>;<y-axis title>"
void SetTitle(const std::string& title) { fTitle = title; }
+ /// \brief String representation of the histogram/profile
+ //std::string ToString() const
+ //{
+ // std::stringstream msg;
+ // msg <<
+ //}
+
private:
friend class boost::serialization::access;
template<class Archive>
+ /// \brief Serialization rule
void serialize(Archive& ar, const unsigned int /*version*/)
{
ar& fHistogram;
@@ -148,7 +178,7 @@ namespace cbm::algo::qa
/// \brief Default constructor
Histogram() = default;
- explicit Histogram(const Histogram<Rank, Axes>& h)
+ explicit Histogram(const Histogram<Axes, Storage>& h)
: fHistogram(h.fHistogram)
, fName(h.fName)
, fTitle(h.fTitle)
@@ -179,8 +209,10 @@ namespace cbm::algo::qa
int fEntries = 0; ///< Number of entries
};
- using BaseH1D = Histogram<1, std::tuple<RegularAxis_t>>;
- using BaseH2D = Histogram<2, std::tuple<RegularAxis_t, RegularAxis_t>>;
+ using BaseH1D = Histogram<Axes1D_t, HistStorage_t>;
+ using BaseH2D = Histogram<Axes2D_t, HistStorage_t>;
+ using BaseProf1D = Histogram<Axes1D_t, ProfStorage_t>;
+ using BaseProf2D = Histogram<Axes2D_t, ProfStorage_t>;
/// \class cbm::algo::qa::H1D
/// \brief 1D-histogram
@@ -214,12 +246,12 @@ namespace cbm::algo::qa
H1D& operator=(H1D&&) = default;
/// \brief Fills histogram
- /// \param value Value
- /// \param weight Weight
+ /// \param x Value
+ /// \param w Weight
/// \return Bin number TODO: implement
- int Fill(double value, double weight = 1.)
+ int Fill(double x, double w = 1.)
{
- fHistogram(value, bh::weight(weight));
+ fHistogram(x, bh::weight(w));
++fEntries;
return -1;
}
@@ -236,7 +268,12 @@ namespace cbm::algo::qa
return std::sqrt(std::fabs(fHistogram.at(Histogram::GetBinBH(iBin))));
}
+ /// \brief Gets underlying bin accumulator
+ /// \param iBin Bin index
+ auto GetBinAccumulator(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)); }
+
private:
+ /// \brief Serialization function
friend class boost::serialization::access;
template<class Archive>
void serialize(Archive& ar, const unsigned int /*version*/)
@@ -280,17 +317,25 @@ namespace cbm::algo::qa
H2D& operator=(H2D&&) = default;
/// \brief Fills histogram
- /// \param valueX Value along x-axis
- /// \param valueY Value along y-axis
- /// \param weight Weight
+ /// \param x Value along x-axis
+ /// \param y Value along y-axis
+ /// \param w Weight
/// \return Bin number TODO: implement
- int Fill(double valueX, double valueY, double weight = 1.)
+ int Fill(double x, double y, double w = 1.)
{
- fHistogram(valueX, valueY, bh::weight(weight));
+ fHistogram(x, y, bh::weight(w));
++fEntries;
return -1;
}
+ /// \brief Gets underlying bin accumulator
+ /// \param iBinX Bin index along x-axis
+ /// \param iBinY Bin index along y-axis
+ auto GetBinAccumulator(uint32_t iBinX, uint32_t iBinY) const
+ {
+ return fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY));
+ }
+
/// \brief Gets bin content
/// \param iBinX Bin index in x-axis
/// \param iBinY Bin index in y-axis
@@ -308,6 +353,7 @@ namespace cbm::algo::qa
}
private:
+ /// \brief Serialization rule
friend class boost::serialization::access;
template<class Archive>
void serialize(Archive& ar, const unsigned int /*version*/)
@@ -316,12 +362,268 @@ namespace cbm::algo::qa
}
};
- //class Prof1D: public Histogram<1, std::tuple<RegularAxis_t>> {
+ // NOTE: Boost::histogram by default uses a different error calculation approach as ROOT TProfile. TODO: investigate
+
+ class Prof1D : public BaseProf1D {
+ public:
+ /// \brief Constructor for 2D-histogram
+ /// \param name Name of histogram
+ /// \param title Title of histogram
+ /// \param nBinsX Number of x bins
+ /// \param xMin Lower x bound
+ /// \param xMax Upper x bound
+ /// \param yMin Lower y bound (optional)
+ /// \param yMax Upper y bound (optional)
+ ///
+ /// \note If the yMin and yMax are defined (yMin != yMax), all the passed values outside the [yMin, yMax] range
+ /// will be ignored.
+ 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)
+ {
+ }
+
+ /// \brief Default constructor
+ Prof1D() = default;
+
+ /// \brief Copy constructor
+ Prof1D(const Prof1D&) = default;
+
+ /// \brief Move constructor
+ Prof1D(Prof1D&&) = default;
+
+ /// \brief Copy assignment operator
+ Prof1D& operator=(const Prof1D&) = default;
+
+ /// \brief Move assignment operator
+ Prof1D& operator=(Prof1D&&) = default;
+
+ /// \brief Fills histogram
+ /// \param x Value along x-axis
+ /// \param y Value along y-axis
+ /// \param w Weight
+ /// \return Bin number TODO: implement
+ int Fill(double x, double y, double w = 1.)
+ {
+ /// Skip measurement, if it goes beyond the defined [fYmin, fYmax] range
+ if ((fYmin != fYmax) && (y < fYmin || y > fYmax || std::isnan(y))) {
+ return -1;
+ }
+
+ fHistogram(x, bh::sample(y), bh::weight(w));
+ fTotSumWX += w * x;
+ fTotSumWX2 += w * x * x;
+ ++fEntries;
+ return -1;
+ }
+
+ /// \brief Gets underlying bin accumulator
+ /// \param iBin Bin index along x-axis
+ auto GetBinAccumulator(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)); }
+
+ /// \brief Gets bin content
+ /// \param iBin Bin index
+ double GetBinContent(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)).GetMean(); }
+
+ /// \brief Gets bin entries
+ /// \param iBin Bin index
+ double GetBinCount(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)).GetEffCount(); }
+
+ /// \brief Gets bin error
+ /// \param iBin Bin index
+ double GetBinError(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)).GetSEM(); }
+
+ /// \brief Gets y-axis lower bound
+ double GetMinY() const { return fYmin; }
+
+ /// \brief Gets y-axis lower bound
+ double GetMaxY() const { return fYmax; }
+
+ /// \brief Gets total sum of weight over x products
+ double GetTotSumWX() const { return fTotSumWX; }
+
+ /// \brief Gets total sum of weight over squared x products
+ double GetTotSumWX2() const { return fTotSumWX2; }
+
+ /// \brief Resets the profile
+ void Reset()
+ {
+ BaseProf1D::Reset();
+ fTotSumWX = 0;
+ fTotSumWX2 = 0;
+ }
+
+ private:
+ /// \brief Serialization rule
+ 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;
+ ar& fTotSumWX;
+ ar& fTotSumWX2;
+ }
+
+ double fYmin = 0.;
+ double fYmax = 0.;
+ double fTotSumWX = 0.; ///< Total sum (over all bins) of weight over x products
+ double fTotSumWX2 = 0.; ///< Total sum (over all bins) of weight over square x products
+ };
+
+ class Prof2D : public BaseProf2D {
+ public:
+ /// \brief Constructor for 2D-histogram
+ /// \param name Name of histogram
+ /// \param title Title of histogram
+ /// \param nBinsX Number of x bins
+ /// \param xMin Lower x bound
+ /// \param xMax Upper x bound
+ /// \param nBinsY Number of y bins
+ /// \param yMin Lower y bound
+ /// \param yMax Upper y bound
+ /// \param zMin Lower z bound (optional)
+ /// \param zMax Upper z bound (optional)
+ ///
+ /// \note If the zMin and zMax are defined (zMin != zMax), all the passed values outside the [zMin, zMax] range
+ /// will be ignored.
+ 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)
+ {
+ }
+
+ /// \brief Default constructor
+ Prof2D() = default;
- //};
+ /// \brief Copy constructor
+ Prof2D(const Prof2D&) = default;
- //class Prof2D: public Histogram<2, std::tuple<RegularAxis_t, RegularAxis_t>> {
+ /// \brief Move constructor
+ Prof2D(Prof2D&&) = default;
- //};
+ /// \brief Copy assignment operator
+ Prof2D& operator=(const Prof2D&) = default;
+
+ /// \brief Move assignment operator
+ Prof2D& operator=(Prof2D&&) = default;
+
+ /// \brief Fills histogram
+ /// \param x Value along x-axis
+ /// \param y Value along y-axis
+ /// \param z Value along z-axis
+ /// \param w Weight
+ /// \return Bin number TODO: implement
+ int Fill(double x, double y, double z, double w = 1.)
+ {
+ /// Skip measurement, if it goes beyond the defined [fYmin, fYmax] range
+ if ((fZmin != fZmax) && (z < fZmin || z > fZmax || std::isnan(z))) {
+ return -1;
+ }
+
+ fHistogram(x, y, bh::sample(z), bh::weight(w));
+ fTotSumWX += w * x;
+ fTotSumWY += w * y;
+ fTotSumWX2 += w * x * x;
+ fTotSumWXY += w * x * y;
+ fTotSumWY2 += w * y * y;
+ ++fEntries;
+ return -1;
+ }
+
+ /// \brief Gets underlying bin accumulator
+ /// \param iBinX Bin index along x-axis
+ /// \param iBinY Bin index along y-axis
+ auto GetBinAccumulator(uint32_t iBinX, uint32_t iBinY) const
+ {
+ return fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY));
+ }
+
+ /// \brief Gets bin content
+ /// \param iBinX Bin index along x-axis
+ /// \param iBinY Bin index along y-axis
+ double GetBinContent(uint32_t iBinX, uint32_t iBinY) const
+ {
+ return fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY)).GetMean();
+ }
+
+ /// \brief Gets bin entries
+ /// \param iBinX Bin index along x-axis
+ /// \param iBinY Bin index along y-axis
+ double GetBinCount(uint32_t iBinX, uint32_t iBinY) const
+ {
+ return fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY)).GetEffCount();
+ }
+
+ /// \brief Gets bin error
+ /// \param iBinX Bin index along x-axis
+ /// \param iBinY Bin index along y-axis
+ double GetBinError(uint32_t iBinX, uint32_t iBinY) const
+ {
+ return fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY)).GetSEM();
+ }
+
+ /// \brief Gets z-axis lower bound
+ double GetMinZ() const { return fZmin; }
+
+ /// \brief Gets z-axis lower bound
+ double GetMaxZ() const { return fZmax; }
+
+ /// \brief Gets total sum of weight over x products
+ double GetTotSumWX() const { return fTotSumWX; }
+
+ /// \brief Gets total sum of weight over y products
+ double GetTotSumWY() const { return fTotSumWY; }
+
+ /// \brief Gets total sum of weight over squared x products
+ double GetTotSumWX2() const { return fTotSumWX2; }
+
+ /// \brief Gets total sum of weight over x over y products
+ double GetTotSumWXY() const { return fTotSumWXY; }
+
+ /// \brief Gets total sum of weight over squared y products
+ double GetTotSumWY2() const { return fTotSumWY2; }
+
+ /// \brief Resets the profile
+ void Reset()
+ {
+ BaseProf2D::Reset();
+ fTotSumWX = 0;
+ fTotSumWY = 0;
+ fTotSumWX2 = 0;
+ fTotSumWXY = 0;
+ fTotSumWY2 = 0;
+ }
+
+ private:
+ /// \brief Serialization rule
+ 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;
+ ar& fTotSumWX;
+ ar& fTotSumWY;
+ ar& fTotSumWX2;
+ ar& fTotSumWXY;
+ ar& fTotSumWY2;
+ }
+
+ double fZmin = 0.;
+ double fZmax = 0.;
+ double fTotSumWX = 0.; ///< Total sum (over all bins) of weight over x products
+ double fTotSumWY = 0.; ///< Total sum (over all bins) of weight over y products
+ double fTotSumWX2 = 0.; ///< Total sum (over all bins) of weight over square x products
+ double fTotSumWXY = 0.; ///< Total sum (over all bins) of weight over square x products
+ double fTotSumWY2 = 0.; ///< Total sum (over all bins) of weight over x over y products
+ };
} // namespace cbm::algo::qa
diff --git a/algo/qa/HistogramContainer.cxx b/algo/qa/HistogramContainer.cxx
index a0fe588738daa6f7afd21b3648cf98509c0930b0..08ef1af72af188515c797512ca27f8d9d2a58835 100644
--- a/algo/qa/HistogramContainer.cxx
+++ b/algo/qa/HistogramContainer.cxx
@@ -21,4 +21,10 @@ void HistogramContainer::Reset()
for (auto& h : fvH2) {
h.Reset();
}
+ //for (auto& h : fvP1) {
+ // h.Reset();
+ //}
+ //for (auto& h : fvP2) {
+ // h.Reset();
+ //}
}
diff --git a/algo/qa/HistogramContainer.h b/algo/qa/HistogramContainer.h
index 11dc4316cd80341741a29df198737ec92151c810..0420d87e8dff49f08f207be8065515c4b70062f5 100644
--- a/algo/qa/HistogramContainer.h
+++ b/algo/qa/HistogramContainer.h
@@ -22,6 +22,8 @@ namespace cbm::algo::qa
struct HistogramContainer {
std::forward_list<qa::H1D> fvH1 = {}; ///< List of 1D-histograms
std::forward_list<qa::H2D> fvH2 = {}; ///< List of 2D-histograms
+ //std::forward_list<qa::Prof1D> fvP1 = {}; ///< List of 1D-profiles
+ //std::forward_list<qa::Prof2D> fvP2 = {}; ///< List of 2D-profiles
/// \brief Resets the histograms
void Reset();
@@ -33,6 +35,8 @@ namespace cbm::algo::qa
{
ar& fvH1;
ar& fvH2;
+ //ar& fvP1;
+ //ar& fvP2;
}
};
} // namespace cbm::algo::qa
diff --git a/algo/qa/PadConfig.h b/algo/qa/PadConfig.h
index 9c2e7bcf1a5b52bf5de08524b5cb087782281e9b..1c4d3cb9926e28581897d1ad96783b18304152c2 100644
--- a/algo/qa/PadConfig.h
+++ b/algo/qa/PadConfig.h
@@ -63,11 +63,15 @@ namespace cbm::algo::qa
}
/// \brief Registers a histogram in the pad
+ /// \tparam Hist Histogram class
/// \param hist Histogram object
/// \param opt Draw options for the histogram
// TODO: implement a single function
- void RegisterHistogram(const H1D* hist, std::string_view opt) { RegisterObject(hist->GetName(), opt); }
- void RegisterHistogram(const H2D* hist, std::string_view opt) { RegisterObject(hist->GetName(), opt); }
+ template<class Hist>
+ void RegisterHistogram(const Hist* hist, std::string_view opt)
+ {
+ RegisterObject(hist->GetName(), opt);
+ }
/// \brief Returns message config
std::string ToString() const;
diff --git a/core/qa/CbmQaBaseLinkDef.h b/core/qa/CbmQaBaseLinkDef.h
index 9e903ee0005238f3c6c45944dd04952c300b14d4..2d33214ac1ed04ddfd13d9b57853a683aacdbfd3 100644
--- a/core/qa/CbmQaBaseLinkDef.h
+++ b/core/qa/CbmQaBaseLinkDef.h
@@ -18,6 +18,8 @@
// create streamers automatically
#pragma link C++ class cbm::algo::qa::H1D - ;
#pragma link C++ class cbm::algo::qa::H2D - ;
+#pragma link C++ class cbm::algo::qa::Prof1D - ;
+#pragma link C++ class cbm::algo::qa::Prof2D - ;
#pragma link C++ class cbm::qa::OnlineInterface + ;
#pragma link C++ class CbmQaEff + ;
#pragma link C++ class CbmQaPieSlice + ;
diff --git a/core/qa/CbmQaOnlineInterface.cxx b/core/qa/CbmQaOnlineInterface.cxx
index e479bb282c56f20c7c7f745f48af7703452d1f08..0c6cebba0c086dc3307534da55c6a055a0a7cb39 100644
--- a/core/qa/CbmQaOnlineInterface.cxx
+++ b/core/qa/CbmQaOnlineInterface.cxx
@@ -12,8 +12,16 @@
#include "TH1D.h"
#include "TH2D.h"
+#include "TProfile.h"
+#include "TProfile2D.h"
+using cbm::algo::qa::H1D;
+using cbm::algo::qa::H2D;
+using cbm::algo::qa::Prof1D;
+using cbm::algo::qa::Prof2D;
using cbm::qa::OnlineInterface;
+using cbm::qa::TProfileAccessor;
+
// ---------------------------------------------------------------------------------------------------------------------
//
@@ -57,3 +65,44 @@ TH2D* OnlineInterface::ROOTHistogram(const cbm::algo::qa::H2D& hist)
pRes->SetEntries(hist.GetEntries());
return pRes;
}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TProfile* OnlineInterface::ROOTHistogram(const Prof1D& prof)
+{
+ const char* name = prof.GetName().c_str();
+ const char* titl = prof.GetTitle().c_str();
+ int nBinsX = prof.GetNbinsX();
+ double xMin = prof.GetMinX();
+ double xMax = prof.GetMaxX();
+ double yMin = prof.GetMinY();
+ double yMax = prof.GetMaxY();
+ bool add = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(false);
+ auto* pRes = new TProfileAccessor<TProfile, Prof1D>(name, titl, nBinsX, xMin, xMax, yMin, yMax);
+ TH1::AddDirectory(add);
+ pRes->SetFromQaProf(prof);
+ return pRes;
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+TProfile2D* OnlineInterface::ROOTHistogram(const Prof2D& prof)
+{
+ const char* name = prof.GetName().c_str();
+ const char* titl = prof.GetTitle().c_str();
+ int nBinsX = prof.GetNbinsX();
+ double xMin = prof.GetMinX();
+ double xMax = prof.GetMaxX();
+ int nBinsY = prof.GetNbinsY();
+ double yMin = prof.GetMinY();
+ double yMax = prof.GetMaxY();
+ double zMin = prof.GetMinZ();
+ double zMax = prof.GetMaxZ();
+ bool add = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(false);
+ auto* pRes = new TProfileAccessor<TProfile2D, Prof2D>(name, titl, nBinsX, xMin, xMax, nBinsY, yMin, yMax, zMin, zMax);
+ TH1::AddDirectory(add);
+ pRes->SetFromQaProf(prof);
+ return pRes;
+}
diff --git a/core/qa/CbmQaOnlineInterface.h b/core/qa/CbmQaOnlineInterface.h
index 2b995df8594061278fed69e3eef7cbab7efe8db3..db19c16aa117febc815e7f4c72c5f00bf240ea22 100644
--- a/core/qa/CbmQaOnlineInterface.h
+++ b/core/qa/CbmQaOnlineInterface.h
@@ -12,29 +12,112 @@
#pragma once
#include "Histogram.h"
+#include "Logger.h"
+#include "TProfile.h"
+#include "TProfile2D.h"
+
+#include <type_traits>
class TH1D;
class TH2D;
+class TProfile2D;
namespace cbm::algo::qa
{
class H1D;
class H2D;
+ class Prof1D;
+ class Prof2D;
} // namespace cbm::algo::qa
namespace cbm::qa
{
+ /// \class TProfileAccessor
+ /// \brief Helper class to access protected fields of TProfile
+ ///
+ /// The ROOT TProfile classes do not alow to set a profile from outside, so we need to set all the
+ /// fields directly.
+ template<class RootProfile, class QaProfile>
+ class TProfileAccessor : public RootProfile {
+
+ public:
+ template<class... Args>
+ TProfileAccessor(Args... args) : RootProfile(args...)
+ {
+ }
+
+ /// \brief Sets fields from cbm::algo::qa::Prof1D
+ void SetFromQaProf(const QaProfile& prof);
+ };
+
/// \class OnlineInterface
/// \brief A collection of tools for online QA object conversions
class OnlineInterface {
public:
- /// \brief Converts histogram H1 to ROOT histogram TH1D
+ /// \brief Converts histogram H1D to ROOT histogram TH1D
/// \param hist 1D-histogram
/// \note Allocates memory on heap for the ROOT histogram
static TH1D* ROOTHistogram(const cbm::algo::qa::H1D& hist);
- /// \brief Converts histogram H2 to ROOT histogram TH2D
+ /// \brief Converts histogram H2D to ROOT histogram TH2D
/// \param hist 1D-histogram
/// \note Allocates memory on heap for the ROOT histogram
static TH2D* ROOTHistogram(const cbm::algo::qa::H2D& hist);
+
+ /// \brief Converts profile Prof1D to ROOT TProfile
+ /// \param prof 1D-profile
+ /// \note Allocates memory on heap for the ROOT histogram
+ static TProfile* ROOTHistogram(const cbm::algo::qa::Prof1D& prof);
+
+ /// \brief Converts profile Prof2D to ROOT TProfile2D
+ /// \param prof 2D-profile
+ /// \note Allocates memory on heap for the ROOT histogram
+ static TProfile2D* ROOTHistogram(const cbm::algo::qa::Prof2D& prof);
};
} // namespace cbm::qa
+
+
+namespace cbm::qa
+{
+ // -------------------------------------------------------------------------------------------------------------------
+ //
+ template<class RootProfile, class QaProfile>
+ void TProfileAccessor<RootProfile, QaProfile>::SetFromQaProf(const QaProfile& prof)
+ {
+ this->Sumw2();
+ auto SetBin = [&](int iBinRoot, const auto& binQa) {
+ this->fArray[iBinRoot] = binQa.GetSumWV();
+ this->fBinEntries.fArray[iBinRoot] = binQa.GetSumW();
+ this->fSumw2.fArray[iBinRoot] = binQa.GetSumWV2();
+ this->fBinSumw2.fArray[iBinRoot] = binQa.GetSumW2();
+
+ this->fTsumw += binQa.GetSumW();
+ this->fTsumw2 += binQa.GetSumW2();
+ this->fTsumwy += binQa.GetSumWV();
+ this->fTsumwy2 += binQa.GetSumWV2();
+ };
+
+ if constexpr (std::is_same_v<RootProfile, TProfile>) {
+ for (int iBinX = 0; iBinX <= this->GetNbinsX() + 1; ++iBinX) {
+ int iBin = this->GetBin(iBinX);
+ auto bin = prof.GetBinAccumulator(iBinX);
+ SetBin(iBin, bin);
+ }
+ }
+ else if constexpr (std::is_same_v<RootProfile, TProfile2D>) {
+ for (int iBinX = 0; iBinX <= this->GetNbinsX() + 1; ++iBinX) {
+ for (int iBinY = 0; iBinY <= this->GetNbinsY() + 1; ++iBinY) {
+ int iBin = this->GetBin(iBinX, iBinY);
+ auto bin = prof.GetBinAccumulator(iBinX, iBinY);
+ SetBin(iBin, bin);
+ }
+ }
+ this->fTsumwy = prof.GetTotSumWY();
+ this->fTsumwy2 = prof.GetTotSumWY2();
+ this->fTsumwxy = prof.GetTotSumWXY();
+ }
+ this->fTsumwx = prof.GetTotSumWX();
+ this->fTsumwx2 = prof.GetTotSumWX2();
+ this->SetEntries(prof.GetEntries());
+ }
+
+} // namespace cbm::qa