diff --git a/algo/CMakeLists.txt b/algo/CMakeLists.txt
index 4dd79e72e3ddc13d0b463252f11a3374ef8e3a72..4abc148e95b04a5ef0aab0ee5332601f23703194 100644
--- a/algo/CMakeLists.txt
+++ b/algo/CMakeLists.txt
@@ -162,6 +162,7 @@ target_include_directories(Algo
${CMAKE_CURRENT_SOURCE_DIR}/unpack
${CMAKE_CURRENT_SOURCE_DIR}/detectors
${CMAKE_CURRENT_SOURCE_DIR}/qa
+ ${CMAKE_CURRENT_SOURCE_DIR}/ca
${CMAKE_CURRENT_SOURCE_DIR}/ca/qa
${CMAKE_CURRENT_SOURCE_DIR}
${CMAKE_SOURCE_DIR}/core/data/global
diff --git a/algo/ca/TrackingChain.cxx b/algo/ca/TrackingChain.cxx
index ea1acc6fcf6f64d993f19b6662c76667890909f1..4136c65fcd0afad876b77d3136fefe7a90557fb2 100644
--- a/algo/ca/TrackingChain.cxx
+++ b/algo/ca/TrackingChain.cxx
@@ -66,10 +66,12 @@ void TrackingChain::Init()
fCaFramework.ReceiveParameters(std::move(parameters));
// ------ Initialize QA modules
- if (fOutputQa.IsSenderDefined()) {
+ if (fInputQa.IsSenderDefined()) {
fInputQa.RegisterParameters(&fCaFramework.GetParameters());
fInputQa.Init();
+ }
+ if (fOutputQa.IsSenderDefined()) {
fOutputQa.RegisterParameters(&fCaFramework.GetParameters());
fOutputQa.Init();
}
@@ -157,12 +159,14 @@ TrackingChain::Output_t TrackingChain::PrepareOutput()
output.monitorData = fCaMonitorData;
// QA
- if (fOutputQa.IsSenderDefined()) {
+ if (fInputQa.IsSenderDefined()) {
fInputQa.RegisterInputData(&fCaFramework.GetInputData());
fInputQa.RegisterTracks(&output.tracks);
fInputQa.RegisterRecoHitIndices(&fCaFramework.fRecoHits);
fInputQa.Exec();
+ }
+ if (fOutputQa.IsSenderDefined()) {
fOutputQa.RegisterInputData(&fCaFramework.GetInputData());
fOutputQa.RegisterTracks(&output.tracks);
fOutputQa.RegisterRecoHitIndices(&fCaFramework.fRecoHits);
diff --git a/algo/ca/qa/CaInputQa.cxx b/algo/ca/qa/CaInputQa.cxx
index b8e6352854f1b03e059f5bc4723480bc678ff998..c63cb0cea9a9c37c565a4223ee03fb8bcb6b51d1 100644
--- a/algo/ca/qa/CaInputQa.cxx
+++ b/algo/ca/qa/CaInputQa.cxx
@@ -9,10 +9,10 @@
#include "CaInputQa.h"
-#include "../TrackingDefs.h"
#include "CaConstants.h"
#include "CaInputData.h"
#include "CaParameters.h"
+#include "TrackingDefs.h"
#include <algorithm>
#include <limits>
@@ -50,132 +50,130 @@ void InputQa::Init()
// ------ Init histograms
// Occupancy distributions
- {
- fvphHitOccupXY.resize(nSt);
- fvphHitOccupZX.resize(nSt);
- fvphHitOccupZY.resize(nSt);
- fvphHitUsageXY.resize(nSt);
-
- // Station sizes in transverse plane
- std::vector<double> vMinX(nSt);
- std::vector<double> vMaxX(nSt);
- std::vector<double> vMinY(nSt);
- std::vector<double> vMaxY(nSt);
-
- int nBinsXY = 200;
- int nBinsZ = 400;
-
- for (int iSt = 0; iSt < nSt; ++iSt) {
- const auto& station = fpParameters->GetStation(iSt);
- vMaxX[iSt] = station.GetXmax<double>();
- vMinX[iSt] = station.GetXmin<double>();
- vMaxY[iSt] = station.GetYmax<double>();
- vMinY[iSt] = station.GetYmin<double>();
- double dy = (vMaxY[iSt] - vMinY[iSt]) * kXYZMargin;
- double dx = (vMaxX[iSt] - vMinX[iSt]) * kXYZMargin;
- vMaxX[iSt] += dx;
- vMinX[iSt] -= dx;
- vMaxY[iSt] += dy;
- vMinY[iSt] -= dy;
- }
- // Station max
- double xMinA = *std::min_element(vMinX.begin(), vMinX.end());
- double xMaxA = *std::max_element(vMaxX.begin(), vMaxX.end());
- double yMinA = *std::min_element(vMinY.begin(), vMinY.end());
- double yMaxA = *std::max_element(vMaxY.begin(), vMaxY.end());
- double zMinA = fpParameters->GetStation(0).GetZ<double>();
- double zMaxA = fpParameters->GetStation(nSt - 1).GetZ<double>();
- {
- double dz = (zMaxA - zMinA) * kXYZMargin;
- zMinA -= dz;
- zMaxA += dz;
- }
- for (int iSt = 0; iSt < nSt; ++iSt) {
- int iStGeo = fpParameters->GetActiveToGeoMap()[iSt];
- auto [detID, iStLoc] = fpParameters->GetGeoToLocalIdMap()[iStGeo];
- for (auto hitSet : kHitSets) {
- auto setNm = EHitSet::Input == hitSet ? "input" : "used";
- auto setTl = EHitSet::Input == hitSet ? "Input" : "Used";
- {
- auto name = format("hit_{}_occup_xy_sta_{}", setNm, iSt);
- auto titl = format("{} hit occupancy in XY plane for station {} ({}{});x [cm];y [cm]", setTl, iSt,
- kDetName[detID], iStLoc);
- fvphHitOccupXY[iSt][hitSet] =
- fQaData.MakeObj<H2D>(name, titl, nBinsXY, vMinX[iSt], vMaxX[iSt], nBinsXY, vMinY[iSt], vMaxY[iSt]);
- }
- {
- auto name = format("hit_{}_occup_zx_sta_{}", setNm, iSt);
- auto titl = format("{} hit occupancy in ZX plane for station {} ({}{});z [cm];x [cm]", setTl, iSt,
- kDetName[detID], iStLoc);
- fvphHitOccupZX[iSt][hitSet] = fQaData.MakeObj<H2D>(name, titl, nBinsZ, zMinA, zMaxA, nBinsXY, xMinA, xMaxA);
- }
- {
- auto name = format("hit_{}_occup_zy_sta_{}", setNm, iSt);
- auto titl = format("{} hit occupancy in ZY plane for station {} ({}{});z [cm];y [cm]", setTl, iSt,
- kDetName[detID], iStLoc);
- fvphHitOccupZY[iSt][hitSet] = fQaData.MakeObj<H2D>(name, titl, nBinsZ, zMinA, zMaxA, nBinsXY, yMinA, yMaxA);
- }
- }
- {
- auto name = format("hit_usage_xy_sta_{}", iSt);
- auto titl = format("Hit usage in XY plane for station {} ({}{});x [cm];y [cm]", iSt, kDetName[detID], iStLoc);
- fvphHitUsageXY[iSt] =
- fQaData.MakeObj<Prof2D>(name, titl, nBinsXY, vMinX[iSt], vMaxX[iSt], nBinsXY, vMinY[iSt], vMaxY[iSt], 0., 1.);
- }
- }
+ fvphHitOccupXY.resize(nSt);
+ fvphHitOccupZX.resize(nSt);
+ fvphHitOccupZY.resize(nSt);
+ fvphHitUsageXY.resize(nSt);
+
+ // Station sizes in transverse plane
+ std::vector<double> vMinX(nSt);
+ std::vector<double> vMaxX(nSt);
+ std::vector<double> vMinY(nSt);
+ std::vector<double> vMaxY(nSt);
+
+ int nBinsXY = 200;
+ int nBinsZ = 400;
+
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ const auto& station = fpParameters->GetStation(iSt);
+ vMaxX[iSt] = station.GetXmax<double>();
+ vMinX[iSt] = station.GetXmin<double>();
+ vMaxY[iSt] = station.GetYmax<double>();
+ vMinY[iSt] = station.GetYmin<double>();
+ double dy = (vMaxY[iSt] - vMinY[iSt]) * kXYZMargin;
+ double dx = (vMaxX[iSt] - vMinX[iSt]) * kXYZMargin;
+ vMaxX[iSt] += dx;
+ vMinX[iSt] -= dx;
+ vMaxY[iSt] += dy;
+ vMinY[iSt] -= dy;
}
-
- // ----- Init canvases
+ // Station max
+ double xMinA = *std::min_element(vMinX.begin(), vMinX.end());
+ double xMaxA = *std::max_element(vMaxX.begin(), vMaxX.end());
+ double yMinA = *std::min_element(vMinY.begin(), vMinY.end());
+ double yMaxA = *std::max_element(vMaxY.begin(), vMaxY.end());
+ double zMinA = fpParameters->GetStation(0).GetZ<double>();
+ double zMaxA = fpParameters->GetStation(nSt - 1).GetZ<double>();
{
- // Hit occupancies
+ double dz = (zMaxA - zMinA) * kXYZMargin;
+ zMinA -= dz;
+ zMaxA += dz;
+ }
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ int iStGeo = fpParameters->GetActiveToGeoMap()[iSt];
+ auto [detID, iStLoc] = fpParameters->GetGeoToLocalIdMap()[iStGeo];
for (auto hitSet : kHitSets) {
auto setNm = EHitSet::Input == hitSet ? "input" : "used";
auto setTl = EHitSet::Input == hitSet ? "Input" : "Used";
- { // XY
- auto name = format("ca_hit_{}_occupancy_xy", setNm);
- auto titl = format("{} hit occupancy in different stations in XY plane", setNm);
- auto canv = CanvasConfig(name, titl);
- for (int iSt = 0; iSt < nSt; ++iSt) {
- auto pad = PadConfig();
- pad.RegisterHistogram(fvphHitOccupXY[iSt][hitSet], "colz");
- canv.AddPadConfig(pad);
- }
- fQaData.AddCanvasConfig(canv);
+ {
+ auto name = format("hit_{}_occup_xy_sta_{}", setNm, iSt);
+ auto titl = format("{} hit occupancy in XY plane for station {} ({}{});x [cm];y [cm]", setTl, iSt,
+ kDetName[detID], iStLoc);
+ fvphHitOccupXY[iSt][hitSet] =
+ fQaData.MakeObj<H2D>(name, titl, nBinsXY, vMinX[iSt], vMaxX[iSt], nBinsXY, vMinY[iSt], vMaxY[iSt]);
}
- { // ZX and ZY
- auto name = format("ca_hit_{}_occupancy_zx_zy", setNm);
- auto titl = format("{} hit occupancy in different stations in ZX and ZY planes", setTl);
- auto canv = CanvasConfig(name, titl);
- { // ZX
- auto pad = PadConfig();
- for (int iSt = 0; iSt < nSt; ++iSt) {
- pad.RegisterHistogram(fvphHitOccupZX[iSt][hitSet], (iSt == 0 ? "colz" : "cols same"));
- }
- canv.AddPadConfig(pad);
- }
- { // ZY
- auto pad = PadConfig();
- for (int iSt = 0; iSt < nSt; ++iSt) {
- pad.RegisterHistogram(fvphHitOccupZY[iSt][hitSet], (iSt == 0 ? "colz" : "cols same"));
- }
- canv.AddPadConfig(pad);
- }
- fQaData.AddCanvasConfig(canv);
+ {
+ auto name = format("hit_{}_occup_zx_sta_{}", setNm, iSt);
+ auto titl = format("{} hit occupancy in ZX plane for station {} ({}{});z [cm];x [cm]", setTl, iSt,
+ kDetName[detID], iStLoc);
+ fvphHitOccupZX[iSt][hitSet] = fQaData.MakeObj<H2D>(name, titl, nBinsZ, zMinA, zMaxA, nBinsXY, xMinA, xMaxA);
+ }
+ {
+ auto name = format("hit_{}_occup_zy_sta_{}", setNm, iSt);
+ auto titl = format("{} hit occupancy in ZY plane for station {} ({}{});z [cm];y [cm]", setTl, iSt,
+ kDetName[detID], iStLoc);
+ fvphHitOccupZY[iSt][hitSet] = fQaData.MakeObj<H2D>(name, titl, nBinsZ, zMinA, zMaxA, nBinsXY, yMinA, yMaxA);
}
}
- // Hit usage profiles
{
- auto name = format("ca_hit_usage_xy");
- auto titl = format("Hit usage in different stations in XY plane");
+ auto name = format("hit_usage_xy_sta_{}", iSt);
+ auto titl = format("Hit usage in XY plane for station {} ({}{});x [cm];y [cm]", iSt, kDetName[detID], iStLoc);
+ fvphHitUsageXY[iSt] =
+ fQaData.MakeObj<Prof2D>(name, titl, nBinsXY, vMinX[iSt], vMaxX[iSt], nBinsXY, vMinY[iSt], vMaxY[iSt], 0., 1.);
+ }
+ }
+
+
+ // ----- Init canvases
+ // Hit occupancies
+ for (auto hitSet : kHitSets) {
+ auto setNm = EHitSet::Input == hitSet ? "input" : "used";
+ auto setTl = EHitSet::Input == hitSet ? "Input" : "Used";
+ { // XY
+ auto name = format("ca_hit_{}_occupancy_xy", setNm);
+ auto titl = format("{} hit occupancy in different stations in XY plane", setNm);
auto canv = CanvasConfig(name, titl);
for (int iSt = 0; iSt < nSt; ++iSt) {
auto pad = PadConfig();
- pad.RegisterHistogram(fvphHitUsageXY[iSt], "colz");
+ pad.RegisterHistogram(fvphHitOccupXY[iSt][hitSet], "colz");
canv.AddPadConfig(pad);
}
fQaData.AddCanvasConfig(canv);
}
+ { // ZX and ZY
+ auto name = format("ca_hit_{}_occupancy_zx_zy", setNm);
+ auto titl = format("{} hit occupancy in different stations in ZX and ZY planes", setTl);
+ auto canv = CanvasConfig(name, titl);
+ { // ZX
+ auto pad = PadConfig();
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pad.RegisterHistogram(fvphHitOccupZX[iSt][hitSet], (iSt == 0 ? "colz" : "cols same"));
+ }
+ canv.AddPadConfig(pad);
+ }
+ { // ZY
+ auto pad = PadConfig();
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ pad.RegisterHistogram(fvphHitOccupZY[iSt][hitSet], (iSt == 0 ? "colz" : "cols same"));
+ }
+ canv.AddPadConfig(pad);
+ }
+ fQaData.AddCanvasConfig(canv);
+ }
+ }
+ // Hit usage profiles
+ {
+ auto name = format("ca_hit_usage_xy");
+ auto titl = format("Hit usage in different stations in XY plane");
+ auto canv = CanvasConfig(name, titl);
+ for (int iSt = 0; iSt < nSt; ++iSt) {
+ auto pad = PadConfig();
+ pad.RegisterHistogram(fvphHitUsageXY[iSt], "colz");
+ canv.AddPadConfig(pad);
+ }
+ fQaData.AddCanvasConfig(canv);
}
+
fQaData.Init(fpSender);
}
diff --git a/algo/qa/DigiEventQa.h b/algo/qa/DigiEventQa.h
index 91961d46493f5f0fe9c540b809f32d01f9d0fd77..897cda8ebcb5cbf770e71b55f59f61aaa3275dc6 100644
--- a/algo/qa/DigiEventQa.h
+++ b/algo/qa/DigiEventQa.h
@@ -11,7 +11,7 @@
#include "evbuild/EventBuilderConfig.h"
#include <gsl/span>
-#include <map>
+#include <unordered_map>
#include <vector>
@@ -25,8 +25,8 @@ namespace cbm::algo::evbuild
**/
struct DigiEventQaData {
qa::HistogramContainer fHistContainer;
- std::map<ECbmModuleId, qa::H1D*> fDigiTimeHistos = {};
- size_t fNumEvents = 0;
+ std::unordered_map<ECbmModuleId, qa::H1D*> fDigiTimeHistos = {};
+ size_t fNumEvents = 0;
};
diff --git a/algo/qa/Histogram.h b/algo/qa/Histogram.h
index f2a29770e87c8ef58d90d9ac58c75958f0005d10..d136f46b5fb521f5795d62c35c59900a6a613312 100644
--- a/algo/qa/Histogram.h
+++ b/algo/qa/Histogram.h
@@ -36,10 +36,9 @@ namespace cbm::algo::qa
using Axes1D_t = std::tuple<RegularAxis_t>;
using Axes2D_t = std::tuple<RegularAxis_t, RegularAxis_t>;
- using HistStorage_t = bh::default_storage;
+ using HistStorage_t = bh::weight_storage;
using ProfStorage_t = bh::RootStyleProfileStorage<double>;
-
/// \enum EHistAxis
/// \brief Enumeration for histogram axes
enum class EAxis : unsigned
@@ -49,12 +48,122 @@ namespace cbm::algo::qa
Z
};
+ /// \class TotalSums1D
+ /// \brief Storage for total sums of weights, squared weights, weights over x, weights over squared x
+ class TotalSums1D {
+ public:
+ /// \brief Gets total sum of weights
+ double GetTotSumW() const { return fTotSumW; }
+
+ /// \brief Gets total sum of squared weights
+ double GetTotSumW2() const { return fTotSumW2; }
+
+ /// \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; }
+
+ private:
+ /// \brief Serialization rule
+ friend class boost::serialization::access;
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int /*version*/)
+ {
+ ar& fTotSumW;
+ ar& fTotSumW2;
+ ar& fTotSumWX;
+ ar& fTotSumWX2;
+ }
+
+ protected:
+ /// \brief Updates the sums
+ /// \param x X value
+ /// \param w weight
+ void UpdateTotalSums(double x, double w)
+ {
+ fTotSumW += w;
+ fTotSumW2 += w * w;
+ fTotSumWX += w * x;
+ fTotSumWX2 += w * x * x;
+ }
+
+ /// \brief Resets the sums
+ void Reset()
+ {
+ fTotSumW = 0;
+ fTotSumW2 = 0;
+ fTotSumWX = 0;
+ fTotSumWX2 = 0;
+ }
+
+ double fTotSumW = 0.; ///< Total sum (over all bins) of weights
+ double fTotSumW2 = 0.; ///< Total sum (over all bins) of squared weights
+ 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 TotalSums2D
+ /// \brief TotalSums1D including storage for total sums of w*x*y, w*y, w*y*y products
+ class TotalSums2D : public TotalSums1D {
+ public:
+ /// \brief Gets total sum of weight over squared y products
+ double GetTotSumWXY() const { return fTotSumWXY; }
+
+ /// \brief Gets total sum of weight over y products
+ double GetTotSumWY() const { return fTotSumWY; }
+
+ /// \brief Gets total sum of weight over squared y products
+ double GetTotSumWY2() const { return fTotSumWY2; }
+
+ 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<TotalSums1D>(*this);
+ ar& fTotSumWXY;
+ ar& fTotSumWY;
+ ar& fTotSumWY2;
+ }
+
+ protected:
+ /// \brief Resets the sums
+ void Reset()
+ {
+ TotalSums1D::Reset();
+ fTotSumWXY = 0;
+ fTotSumWY = 0;
+ fTotSumWY2 = 0;
+ }
+
+ /// \brief Updates the sums
+ /// \param x X value
+ /// \param y Y value
+ /// \param w weight
+ 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.; ///< Total sum (over all bins) of weight over y 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
+ };
+
+
/// \class Histogram
/// \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 {
+ /// \tparam Axes The axes type
+ /// \tparam Storage Storage type for the histogram
+ /// \tparam TotalSums Class for storing the total sums of weights for axes (essential to calculate mean and std.dev.)
+ 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>;
@@ -144,7 +253,11 @@ namespace cbm::algo::qa
const std::string& GetTitle() const { return fTitle; }
/// \brief Resets the histogram
- void Reset() { fHistogram.reset(); }
+ void Reset()
+ {
+ TotalSums::Reset();
+ fHistogram.reset();
+ }
/// \brief Sets name
/// \param name Histogram name
@@ -168,6 +281,7 @@ namespace cbm::algo::qa
/// \brief Serialization rule
void serialize(Archive& ar, const unsigned int /*version*/)
{
+ ar& boost::serialization::base_object<TotalSums>(*this);
ar& fHistogram;
ar& fName;
ar& fTitle;
@@ -178,13 +292,15 @@ namespace cbm::algo::qa
/// \brief Default constructor
Histogram() = default;
- explicit Histogram(const Histogram<Axes, Storage>& h)
- : fHistogram(h.fHistogram)
- , fName(h.fName)
- , fTitle(h.fTitle)
- , fEntries(h.fEntries)
- {
- }
+ 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)
+ //{
+ //}
/// \brief Constructor
/// \param bhist Boost histogram object
@@ -209,10 +325,10 @@ namespace cbm::algo::qa
int fEntries = 0; ///< Number of entries
};
- 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>;
+ 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 cbm::algo::qa::H1D
/// \brief 1D-histogram
@@ -225,8 +341,7 @@ namespace cbm::algo::qa
/// \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
+ : BaseH1D(bh::make_weighted_histogram(RegularAxis_t(nBins, xMin, xMax)), name, title)
{
}
@@ -251,22 +366,24 @@ namespace cbm::algo::qa
/// \return Bin number TODO: implement
int Fill(double x, double w = 1.)
{
- fHistogram(x, bh::weight(w));
+ auto cellIt = fHistogram(bh::weight(w), x);
++fEntries;
- return -1;
+ 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;
}
/// \brief Gets bin content
/// \param iBin Bin index
- double GetBinContent(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)); }
+ double GetBinContent(uint32_t iBin) const { return fHistogram.at(Histogram::GetBinBH(iBin)).value(); }
/// \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))));
- }
+ double GetBinError(uint32_t iBin) const { return std::sqrt(fHistogram.at(Histogram::GetBinBH(iBin)).variance()); }
/// \brief Gets underlying bin accumulator
/// \param iBin Bin index
@@ -297,7 +414,8 @@ namespace cbm::algo::qa
/// \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)
+ : BaseH2D(bh::make_weighted_histogram(RegularAxis_t(nBinsX, xMin, xMax), RegularAxis_t(nBinsY, yMin, yMax)), name,
+ title)
{
}
@@ -323,9 +441,21 @@ namespace cbm::algo::qa
/// \return Bin number TODO: implement
int Fill(double x, double y, double w = 1.)
{
- fHistogram(x, y, bh::weight(w));
+ auto cellIt = fHistogram(x, y, bh::weight(w));
++fEntries;
- return -1;
+
+ 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;
}
/// \brief Gets underlying bin accumulator
@@ -341,7 +471,7 @@ namespace cbm::algo::qa
/// \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));
+ return fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY)).value();
}
/// \brief Gets bin error
@@ -349,7 +479,7 @@ namespace cbm::algo::qa
/// \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))));
+ return std::sqrt(fHistogram.at(Histogram::GetBinBH(iBinX), Histogram::GetBinBH(iBinY)).variance());
}
private:
@@ -412,11 +542,15 @@ namespace cbm::algo::qa
return -1;
}
- fHistogram(x, bh::sample(y), bh::weight(w));
- fTotSumWX += w * x;
- fTotSumWX2 += w * x * x;
+ auto cellIt = fHistogram(x, bh::sample(y), bh::weight(w));
++fEntries;
- return -1;
+ 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;
}
/// \brief Gets underlying bin accumulator
@@ -441,20 +575,6 @@ namespace cbm::algo::qa
/// \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;
@@ -464,14 +584,10 @@ namespace cbm::algo::qa
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
+ double fYmin = 0.;
+ double fYmax = 0.;
};
class Prof2D : public BaseProf2D {
@@ -527,14 +643,21 @@ namespace cbm::algo::qa
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;
+ auto cellIt = fHistogram(x, y, bh::sample(z), bh::weight(w));
++fEntries;
- return -1;
+
+ 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;
}
/// \brief Gets underlying bin accumulator
@@ -575,32 +698,6 @@ namespace cbm::algo::qa
/// \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;
@@ -610,20 +707,10 @@ namespace cbm::algo::qa
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
+ double fZmin = 0.;
+ double fZmax = 0.;
};
} // namespace cbm::algo::qa
diff --git a/core/qa/CbmQaOnlineInterface.cxx b/core/qa/CbmQaOnlineInterface.cxx
index 0c6cebba0c086dc3307534da55c6a055a0a7cb39..0b5823b10a0326365cefa128ecd0d7cef69b5d0e 100644
--- a/core/qa/CbmQaOnlineInterface.cxx
+++ b/core/qa/CbmQaOnlineInterface.cxx
@@ -20,7 +20,7 @@ using cbm::algo::qa::H2D;
using cbm::algo::qa::Prof1D;
using cbm::algo::qa::Prof2D;
using cbm::qa::OnlineInterface;
-using cbm::qa::TProfileAccessor;
+using cbm::qa::RootHistogramAccessor;
// ---------------------------------------------------------------------------------------------------------------------
@@ -32,13 +32,9 @@ TH1D* OnlineInterface::ROOTHistogram(const cbm::algo::qa::H1D& hist)
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);
+ auto* pRes = new RootHistogramAccessor<TH1D, H1D>(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());
+ pRes->SetFromQaHistogram(hist);
return pRes;
}
@@ -54,15 +50,10 @@ TH2D* OnlineInterface::ROOTHistogram(const cbm::algo::qa::H2D& hist)
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);
+ auto* pRes = new RootHistogramAccessor<TH2D, H2D>(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());
+ pRes->SetFromQaHistogram(hist);
return pRes;
}
@@ -79,9 +70,9 @@ TProfile* OnlineInterface::ROOTHistogram(const Prof1D& prof)
double yMax = prof.GetMaxY();
bool add = TH1::AddDirectoryStatus();
TH1::AddDirectory(false);
- auto* pRes = new TProfileAccessor<TProfile, Prof1D>(name, titl, nBinsX, xMin, xMax, yMin, yMax);
+ auto* pRes = new RootHistogramAccessor<TProfile, Prof1D>(name, titl, nBinsX, xMin, xMax, yMin, yMax);
TH1::AddDirectory(add);
- pRes->SetFromQaProf(prof);
+ pRes->SetFromQaHistogram(prof);
return pRes;
}
@@ -101,8 +92,9 @@ TProfile2D* OnlineInterface::ROOTHistogram(const Prof2D& prof)
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);
+ auto* pRes =
+ new RootHistogramAccessor<TProfile2D, Prof2D>(name, titl, nBinsX, xMin, xMax, nBinsY, yMin, yMax, zMin, zMax);
TH1::AddDirectory(add);
- pRes->SetFromQaProf(prof);
+ pRes->SetFromQaHistogram(prof);
return pRes;
}
diff --git a/core/qa/CbmQaOnlineInterface.h b/core/qa/CbmQaOnlineInterface.h
index db19c16aa117febc815e7f4c72c5f00bf240ea22..462c90fb20d6e997394014e3424f14e625947ae1 100644
--- a/core/qa/CbmQaOnlineInterface.h
+++ b/core/qa/CbmQaOnlineInterface.h
@@ -13,40 +13,34 @@
#include "Histogram.h"
#include "Logger.h"
+#include "TH1D.h"
+#include "TH2D.h"
#include "TProfile.h"
#include "TProfile2D.h"
+#include <log.hpp>
#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
+ /// \class RootHistogramAccessor
+ /// \brief Helper class to access protected fields of TH1, TH2, TProfile and TProfile2D
///
/// 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 {
+ /// fields directly. Also the access to TH1/TH2 protected fields is needed in order to set the
+ /// total sums of w*x, w*x*x, w*y, w*x*y, and w*y*y.
+ template<class RootHistogram, class QaHistogram>
+ class RootHistogramAccessor : public RootHistogram {
public:
template<class... Args>
- TProfileAccessor(Args... args) : RootProfile(args...)
+ RootHistogramAccessor(Args... args) : RootHistogram(args...)
{
}
- /// \brief Sets fields from cbm::algo::qa::Prof1D
- void SetFromQaProf(const QaProfile& prof);
+ /// \brief Sets fields from qa histogram
+ void SetFromQaHistogram(const QaHistogram& histo);
};
/// \class OnlineInterface
@@ -80,44 +74,75 @@ namespace cbm::qa
{
// -------------------------------------------------------------------------------------------------------------------
//
- template<class RootProfile, class QaProfile>
- void TProfileAccessor<RootProfile, QaProfile>::SetFromQaProf(const QaProfile& prof)
+ template<class RootHistogram, class QaHistogram>
+ void RootHistogramAccessor<RootHistogram, QaHistogram>::SetFromQaHistogram(const QaHistogram& hist)
{
- 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);
+ // Set individual bin properties
+ if constexpr (std::is_same_v<RootHistogram, TProfile> || std::is_same_v<RootHistogram, TProfile2D>) {
+ 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();
+ };
+ if constexpr (std::is_same_v<RootHistogram, TProfile>) {
+ for (int iBinX = 0; iBinX <= this->GetNbinsX() + 1; ++iBinX) {
+ int iBin = this->GetBin(iBinX);
+ auto bin = hist.GetBinAccumulator(iBinX);
+ SetBin(iBin, bin);
+ this->fTsumwy += bin.GetSumWV();
+ this->fTsumwy2 += bin.GetSumWV2();
+ }
+ }
+ else if constexpr (std::is_same_v<RootHistogram, 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 = hist.GetBinAccumulator(iBinX, iBinY);
+ SetBin(iBin, bin);
+ this->fTsumwz += bin.GetSumWV();
+ this->fTsumwz2 += bin.GetSumWV2();
+ }
+ }
}
}
- 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);
+ else if constexpr (std::is_same_v<RootHistogram, TH1D> || std::is_same_v<RootHistogram, TH2D>) {
+ if (hist.GetTotSumW() != hist.GetTotSumW2()) { // some of weights were not equal to 1.
+ if (!this->fSumw2.fN) {
+ this->Sumw2();
+ }
+ }
+ auto SetBin = [&](int iBinRoot, const auto& binQa) {
+ this->fArray[iBinRoot] = binQa.value();
+ if (this->fSumw2.fN) {
+ this->fSumw2.fArray[iBinRoot] = binQa.variance();
+ }
+ };
+ if constexpr (std::is_same_v<RootHistogram, TH1D>) {
+ for (int iBinX = 0; iBinX <= this->GetNbinsX() + 1; ++iBinX) {
+ SetBin(this->GetBin(iBinX), hist.GetBinAccumulator(iBinX));
+ }
+ }
+ else if constexpr (std::is_same_v<RootHistogram, TH2D>) {
+ for (int iBinX = 0; iBinX <= this->GetNbinsX() + 1; ++iBinX) {
+ for (int iBinY = 0; iBinY <= this->GetNbinsY() + 1; ++iBinY) {
+ SetBin(this->GetBin(iBinX, iBinY), hist.GetBinAccumulator(iBinX, iBinY));
+ }
}
}
- this->fTsumwy = prof.GetTotSumWY();
- this->fTsumwy2 = prof.GetTotSumWY2();
- this->fTsumwxy = prof.GetTotSumWXY();
}
- this->fTsumwx = prof.GetTotSumWX();
- this->fTsumwx2 = prof.GetTotSumWX2();
- this->SetEntries(prof.GetEntries());
- }
+ // Set total sums
+ if constexpr (std::is_base_of_v<TH2D, RootHistogram>) {
+ this->fTsumwy = hist.GetTotSumWY();
+ this->fTsumwy2 = hist.GetTotSumWY2();
+ this->fTsumwxy = hist.GetTotSumWXY();
+ }
+ this->fTsumw = hist.GetTotSumW();
+ this->fTsumw2 = hist.GetTotSumW2();
+ this->fTsumwx = hist.GetTotSumWX();
+ this->fTsumwx2 = hist.GetTotSumWX2();
+ this->SetEntries(hist.GetEntries());
+ }
} // namespace cbm::qa