diff --git a/algo/CMakeLists.txt b/algo/CMakeLists.txt index 0d419051e43f3529ede26dd5026da8dd62674455..26ef7d4404885a88806f80d7d74286aa56979c1b 100644 --- a/algo/CMakeLists.txt +++ b/algo/CMakeLists.txt @@ -242,6 +242,7 @@ install( global/RecoResultsInputArchive.h global/RecoResultsOutputArchive.h global/StorableRecoResults.h + qa/Histogram.h ca/TrackingChain.h ca/TrackingChainConfig.h # NOTE: SZh 20.11.2023: diff --git a/algo/qa/Histogram.h b/algo/qa/Histogram.h new file mode 100644 index 0000000000000000000000000000000000000000..273b6578c097e35ef51307e3f4b7052994079c0a --- /dev/null +++ b/algo/qa/Histogram.h @@ -0,0 +1,315 @@ +/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file Histogram.h +/// \date 27.02.2024 +/// \brief Interface to 1D histogram +/// \author Sergei Zharko <s.zharko@gsi.de> + +#pragma once + +#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 <cstdint> +#include <memory> +#include <numeric> +#include <string> +#include <tuple> +#include <type_traits> + +namespace cbm::algo::qa +{ + namespace bh = boost::histogram; + + using RegularAxis_t = bh::axis::regular<>; + + /// \enum EHistAxis + /// \brief Enumeration for histogram axes + enum class EAxis : unsigned + { + X = 0, + Y, + Z + }; + + /// \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> + class Histogram { + //using Hist_t = bh::histogram<Axes, bh::weight_storage>; // TODO: Test in future + using Hist_t = bh::histogram<Axes, bh::default_storage>; + + public: + /// \brief Destructor + ~Histogram() = default; + + /// \brief Gets number of entries + // 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; + } + + /// \brief Gets name + const std::string& GetName() const { return fName; } + + /// \brief Gets number of bins in axis + /// \tparam AxisT An axis + /// \tparam D + 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(); + } + + /// \brief Gets range lower bound of the selected axis + /// \tparam AxisType An axis + 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.); + } + + /// \brief Gets range upper bound of the selected axis + /// \tparam AxisType An axis + 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>())); + } + + /// \brief Gets number of bins for x axis + uint32_t GetNbinsX() const { return GetAxisNbins<EAxis::X>(); } + + /// \brief Gets x-axis lower bound + double GetMinX() const { return GetAxisMin<EAxis::X>(); } + + /// \brief Gets x-axis lower bound + double GetMaxX() const { return GetAxisMax<EAxis::X>(); } + + /// \brief Gets number of bins for y axis + template<unsigned RankCheck = Rank, std::enable_if_t<(RankCheck > 1), bool> = true> + uint32_t GetNbinsY() const + { + return GetAxisNbins<EAxis::Y>(); + } + + /// \brief Gets y-axis lower bound + template<unsigned RankCheck = Rank, std::enable_if_t<(RankCheck > 1), bool> = true> + double GetMinY() const + { + return GetAxisMin<EAxis::Y>(); + } + + /// \brief Gets y-axis lower bound + template<unsigned RankCheck = Rank, std::enable_if_t<(RankCheck > 1), bool> = true> + double GetMaxY() const + { + return GetAxisMax<EAxis::Y>(); + } + + /// \brief Gets title + const std::string& GetTitle() const { return fTitle; } + + /// \brief Resets the histogram + void Reset() { fHistogram.reset(); } + + /// \brief Sets name + /// \param name Histogram name + void SetName(const std::string& name) { fName = name; } + + /// \brief Sets title + /// \param title Histogram title + /// \note ROOT axes convention is applicable: "<main title>;<x-axis title>;<y-axis title>" + void SetTitle(const std::string& title) { fTitle = title; } + + private: + friend class boost::serialization::access; + template<class Archive> + void serialize(Archive& ar, const unsigned int /*version*/) + { + ar& fHistogram; + ar& fName; + ar& fTitle; + ar& fEntries; + } + + protected: + /// \brief Default constructor + Histogram() = default; + + explicit Histogram(const Histogram<Rank, Axes>& h) + : fHistogram(h.fHistogram) + , fName(h.fName) + , fTitle(h.fTitle) + , fEntries(h.fEntries) + { + } + + /// \brief Constructor + /// \param bhist Boost histogram object + /// \param name Name of the histogram + /// \param title Title of the histogram + Histogram(const Hist_t& bhist, const std::string& name, const std::string title) + : fHistogram(bhist) + , fName(name) + , fTitle(title) + , fEntries(0) + { + } + + /// \brief Gets bin index in underlying boost histogram + /// \param iBin Bin in histogram + /// \return Bin in boost histogram + 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; ///< Number of entries + }; + + using BaseH1D = Histogram<1, std::tuple<RegularAxis_t>>; + using BaseH2D = Histogram<2, std::tuple<RegularAxis_t, RegularAxis_t>>; + + /// \class cbm::algo::qa::H1D + /// \brief 1D-histogram + class H1D : public BaseH1D { + public: + /// \brief Constructor for 1D-histogram + /// \param name Name of histogram + /// \param title Title of histogram + /// \param nBins Number of bins + /// \param xMin Lower x bound + /// \param xMax Upper x bound + H1D(const std::string& name, const std::string& title, uint32_t nBins, double xMin, double xMax) + : BaseH1D(bh::make_histogram(RegularAxis_t(nBins, xMin, xMax)), name, + title) // TODO: test bh::make_weighted_histogram + { + } + + /// \brief Default constructor + H1D() = default; + + /// \brief Copy constructor + explicit H1D(const H1D&) = default; + + /// \brief Move constructor + explicit H1D(H1D&&) = default; + + /// \brief Fills histogram + /// \param value Value + /// \param weight Weight + /// \return Bin number TODO: implement + int Fill(double value, double weight = 1.) + { + fHistogram(value, bh::weight(weight)); + ++fEntries; + return -1; + } + + /// \brief Gets bin content + /// \param iBin Bin index + double GetBinContent(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)); } + + /// \brief Gets bin error + /// \param iBin Bin index + double GetBinError(uint32_t iBin) const + { + // TODO: implement Sumw2! + return std::sqrt(std::fabs(fHistogram.at(Histogram::GetBinBH(iBin)))); + } + + 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 cbm::algo::qa::H2D + /// \brief 2D-histogram + class H2D : public BaseH2D { + 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 + 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_histogram(RegularAxis_t(nBinsX, xMin, xMax), RegularAxis_t(nBinsY, yMin, yMax)), name, title) + { + } + + /// \brief Default constructor + H2D() = default; + + /// \brief Copy constructor + explicit H2D(const H2D&) = default; + + /// \brief Move constructor + explicit H2D(H2D&&) = default; + + /// \brief Fills histogram + /// \param valueX Value along x-axis + /// \param valueY Value along y-axis + /// \param weight Weight + /// \return Bin number TODO: implement + int Fill(double valueX, double valueY, double weight = 1.) + { + fHistogram(valueX, valueY, bh::weight(weight)); + ++fEntries; + return -1; + } + + /// \brief Gets bin content + /// \param iBinX Bin index in x-axis + /// \param iBinY Bin index in y-axis + double GetBinContent(uint32_t iBinX, uint32_t iBinY) const + { + return fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY)); + } + + /// \brief Gets bin error + /// \param iBinX Bin index in x-axis + /// \param iBinY Bin index in y-axis + double GetBinError(uint32_t iBinX, uint32_t iBinY) const + { + return std::sqrt(std::fabs(fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY)))); + } + + private: + friend class boost::serialization::access; + template<class Archive> + void serialize(Archive& ar, const unsigned int /*version*/) + { + ar& boost::serialization::base_object<BaseH2D>(*this); + } + }; + + //class Prof1D: public Histogram<1, std::tuple<RegularAxis_t>> { + + //}; + + //class Prof2D: public Histogram<2, std::tuple<RegularAxis_t, RegularAxis_t>> { + + //}; + +} // namespace cbm::algo::qa diff --git a/core/qa/CMakeLists.txt b/core/qa/CMakeLists.txt index 69c773be0852e1f8261675e9cd0aa5cf92e20d2a..456a3b7adf1dd2b256972f835ee8655cf342ca97 100644 --- a/core/qa/CMakeLists.txt +++ b/core/qa/CMakeLists.txt @@ -4,6 +4,7 @@ set(INCLUDE_DIRECTORIES ) set(SRCS + CbmQaOnlineInterface.cxx CbmQaCanvas.cxx CbmQaManager.cxx CbmQaEff.cxx @@ -30,6 +31,7 @@ set(HEADERS set(LIBRARY_NAME CbmQaBase) set(LINKDEF ${LIBRARY_NAME}LinkDef.h) set(PUBLIC_DEPENDENCIES + Algo FairRoot::Base FairLogger::FairLogger external::yaml-cpp @@ -50,6 +52,7 @@ set(INTERFACE_DEPENDENCIES generate_cbm_library() Install(FILES + CbmQaOnlineInterface.h CbmQaCompare.h CbmQaManager.h CbmQaTask.h diff --git a/core/qa/CbmQaBaseLinkDef.h b/core/qa/CbmQaBaseLinkDef.h index 586d32b8265ae4ee651e71ee6cf41da42e098095..9e903ee0005238f3c6c45944dd04952c300b14d4 100644 --- a/core/qa/CbmQaBaseLinkDef.h +++ b/core/qa/CbmQaBaseLinkDef.h @@ -16,6 +16,9 @@ #pragma link C++ class CbmQaPie - ; // create streamers automatically +#pragma link C++ class cbm::algo::qa::H1D - ; +#pragma link C++ class cbm::algo::qa::H2D - ; +#pragma link C++ class cbm::qa::OnlineInterface + ; #pragma link C++ class CbmQaEff + ; #pragma link C++ class CbmQaPieSlice + ; #pragma link C++ class CbmQaHist < TH1F> + ; diff --git a/core/qa/CbmQaOnlineInterface.cxx b/core/qa/CbmQaOnlineInterface.cxx new file mode 100644 index 0000000000000000000000000000000000000000..e479bb282c56f20c7c7f745f48af7703452d1f08 --- /dev/null +++ b/core/qa/CbmQaOnlineInterface.cxx @@ -0,0 +1,59 @@ +/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CbmQaOnlineInterface.cxx +/// \date 28.02.2024 +/// \brief Set of tools for online->ROOT QA-objects conversions (source) +/// \author Sergei Zharko <s.zharko@gsi.de> + +#include "CbmQaOnlineInterface.h" +//#include "Histogram.h" + +#include "TH1D.h" +#include "TH2D.h" + +using cbm::qa::OnlineInterface; + +// --------------------------------------------------------------------------------------------------------------------- +// +TH1D* OnlineInterface::ROOTHistogram(const cbm::algo::qa::H1D& hist) +{ + int nBins = hist.GetNbinsX(); + double xMin = hist.GetMinX(); + double xMax = hist.GetMaxX(); + bool add = TH1::AddDirectoryStatus(); + TH1::AddDirectory(false); + auto* pRes = new TH1D(hist.GetName().c_str(), hist.GetTitle().c_str(), nBins, xMin, xMax); + TH1::AddDirectory(add); + for (int iBinX = 0; iBinX <= pRes->GetNbinsX() + 1; ++iBinX) { + pRes->SetBinContent(iBinX, hist.GetBinContent(iBinX)); + //pRes->SetBinError(iBinX, hist.GetBinError(iBinX)); + } + pRes->SetEntries(hist.GetEntries()); + return pRes; +} + +// --------------------------------------------------------------------------------------------------------------------- +// +TH2D* OnlineInterface::ROOTHistogram(const cbm::algo::qa::H2D& hist) +{ + int nBinsX = hist.GetNbinsX(); + double xMin = hist.GetMinX(); + double xMax = hist.GetMaxX(); + int nBinsY = hist.GetNbinsY(); + double yMin = hist.GetMinY(); + double yMax = hist.GetMaxY(); + bool add = TH1::AddDirectoryStatus(); + TH1::AddDirectory(false); + auto* pRes = new TH2D(hist.GetName().c_str(), hist.GetTitle().c_str(), nBinsX, xMin, xMax, nBinsY, yMin, yMax); + TH1::AddDirectory(add); + for (int iBinX = 0; iBinX <= pRes->GetNbinsX() + 1; ++iBinX) { + for (int iBinY = 0; iBinY <= pRes->GetNbinsY() + 1; ++iBinY) { + pRes->SetBinContent(iBinX, iBinY, hist.GetBinContent(iBinX, iBinY)); + //pRes->SetBinError(iBinX, iBinY, hist.GetBinError(iBinX, iBinY)); + } + } + pRes->SetEntries(hist.GetEntries()); + return pRes; +} diff --git a/core/qa/CbmQaOnlineInterface.h b/core/qa/CbmQaOnlineInterface.h new file mode 100644 index 0000000000000000000000000000000000000000..2b995df8594061278fed69e3eef7cbab7efe8db3 --- /dev/null +++ b/core/qa/CbmQaOnlineInterface.h @@ -0,0 +1,40 @@ +/* Copyright (C) 2024 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt + SPDX-License-Identifier: GPL-3.0-only + Authors: Sergei Zharko [committer] */ + +/// \file CbmQaOnlineInterface.h +/// \date 28.02.2024 +/// \brief Set of tools for online->ROOT QA-objects conversions (header) +/// \author Sergei Zharko <s.zharko@gsi.de> + +// TODO: move to more suitable place (histoserv, for example) + +#pragma once + +#include "Histogram.h" + +class TH1D; +class TH2D; +namespace cbm::algo::qa +{ + class H1D; + class H2D; +} // namespace cbm::algo::qa + +namespace cbm::qa +{ + /// \class OnlineInterface + /// \brief A collection of tools for online QA object conversions + class OnlineInterface { + public: + /// \brief Converts histogram H1 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 + /// \param hist 1D-histogram + /// \note Allocates memory on heap for the ROOT histogram + static TH2D* ROOTHistogram(const cbm::algo::qa::H2D& hist); + }; +} // namespace cbm::qa