From 86983827b6fec1472ef131a00d5f7ab9f6baa95a Mon Sep 17 00:00:00 2001 From: Mario Kruger Date: Mon, 12 Oct 2020 20:11:18 +0200 Subject: [PATCH 1/2] HistogramRegistry: support all root histogram types - renamed readableName to title (to be consistent with root naming conventions) - made nBins in AxisSpec std::optional to seamlessly distinguish between fixed and variable binning per axis -> removed binsEqual in AxisSpec and HistogramConfigSpec - added callSumw2 property to HistogramSpec - rename axis lable -> title (to be consistent with root naming conventions) - added optional name property to axis (this can be helpful in n dimensional histograms so they are not only called axis0, axis1, axis2...) - added possibility to insert axes also after ctor was called - store histograms in registry in form of a std::variant (HistPtr) that can take any root hist pointer type - added static HistFactory helper class that handles the creation of the actual root histograms - HistFactory provides createHist<>() function to instanciate arbitrary histogram based on the blueprint defined in HistogramSpec - HistFactory provides createHistVariant<>() function that creates the histogram and returns a properly casted HistPtr which can be stored in the registry - HistFactory provides a runtime version of the above, which creates the histogram based on the type information stored in HistogramConfigSpec - it was necessary to move from unique_ptrs to shared_ptrs because we must be able to 'downcast' the actual pointer to the few interface types stored in HistPtr (this can not be done with unique ptrs). As a beneficial side effect, this way the pointers created by the HistFactory can also be stored in OutputObj<> (which uses shared_ptrs) - extend HistogramCreationCallback to cover all root histogram types - moved HistogramCallbacks map to .cxx as it is now const static member variable in HistFactory that should be defined only once (otherwise linker will complain) - get() now also requires to specify the interface type (TH1, TH2, TH3, THn, THnSparse, TProfile, TProfile2D or TProfile3D) since type cannot be infered from the variant at compile time - allowed adding histograms to an existing registry - added HistFiller static helper class that contains functions to fill histogams of any type both simply via values or via filtered tables - extended existing fill-from-table functionality to arbitrary number of dimensions - added fill and fillWeight functions to registry that can fill the stored histograms (again providing a version to fill values and one to fill whole filtered table) - renamed HistogramType -> HistType to reduce code verbosity --- Analysis/Tutorials/src/histogramRegistry.cxx | 137 +++- Framework/Core/CMakeLists.txt | 1 + .../include/Framework/HistogramRegistry.h | 600 +++++++++++------- Framework/Core/src/HistogramRegistry.cxx | 37 ++ .../Core/test/benchmark_HistogramRegistry.cxx | 8 +- .../Core/test/test_HistogramRegistry.cxx | 36 +- 6 files changed, 572 insertions(+), 247 deletions(-) create mode 100644 Framework/Core/src/HistogramRegistry.cxx diff --git a/Analysis/Tutorials/src/histogramRegistry.cxx b/Analysis/Tutorials/src/histogramRegistry.cxx index 6a14f665de918..f7f2d92d37657 100644 --- a/Analysis/Tutorials/src/histogramRegistry.cxx +++ b/Analysis/Tutorials/src/histogramRegistry.cxx @@ -27,16 +27,16 @@ struct ATask { "registry", true, { - {"eta", "#eta", {HistogramType::kTH1F, {{102, -2.01, 2.01}}}}, // - {"phi", "#varphi", {HistogramType::kTH1F, {{100, 0., 2. * M_PI}}}} // - } // + {"eta", "#eta", {HistType::kTH1F, {{102, -2.01, 2.01}}}}, // + {"phi", "#varphi", {HistType::kTH1F, {{100, 0., 2. * M_PI}}}} // + } // }; void process(aod::Tracks const& tracks) { for (auto& track : tracks) { - registry.get("eta")->Fill(track.eta()); - registry.get("phi")->Fill(track.phi()); + registry.get("eta")->Fill(track.eta()); + registry.get("phi")->Fill(track.phi()); } } }; @@ -47,9 +47,9 @@ struct BTask { "registry", true, { - {"eta", "#eta", {HistogramType::kTH1F, {{102, -2.01, 2.01}}}}, // - {"ptToPt", "#ptToPt", {HistogramType::kTH2F, {{100, -0.01, 10.01}, {100, -0.01, 10.01}}}} // - } // + {"eta", "#eta", {HistType::kTH1F, {{102, -2.01, 2.01}}}}, // + {"ptToPt", "#ptToPt", {HistType::kTH2F, {{100, -0.01, 10.01}, {100, -0.01, 10.01}}}} // + } // }; void process(aod::Tracks const& tracks) @@ -59,9 +59,128 @@ struct BTask { } }; +struct CTask { + + HistogramRegistry registry{ + "registry", + true, + { + {"1d", "test 1d", {HistType::kTH1I, {{100, -10.0f, 10.0f}}}}, // + {"2d", "test 2d", {HistType::kTH2F, {{100, -10.0f, 10.01f}, {100, -10.0f, 10.01f}}}}, // + {"3d", "test 3d", {HistType::kTH3D, {{100, -10.0f, 10.01f}, {100, -10.0f, 10.01f}, {100, -10.0f, 10.01f}}}}, // + {"4d", "test 4d", {HistType::kTHnC, {{100, -10.0f, 10.01f}, {100, -10.0f, 10.01f}, {100, -10.0f, 10.01f}, {100, -10.0f, 10.01f}}}}, // + {"5d", "test 5d", {HistType::kTHnSparseL, {{10, -10.0f, 10.01f}, {10, -10.0f, 10.01f}, {10, -10.0f, 10.01f}, {10, -10.0f, 10.01f}, {10, -10.0f, 10.01f}}}}, // + } // + }; + + void init(o2::framework::InitContext&) + { + registry.add({"7d", "test 7d", {HistType::kTHnC, {{3, -10.0f, 10.01f}, {3, -10.0f, 10.01f}, {3, -10.0f, 10.01f}, {3, -10.0f, 10.01f}, {3, -10.0f, 10.01f}, {3, -10.0f, 10.01f}, {3, -10.0f, 10.01f}}}}); + + registry.add({"6d", "test 6d", {HistType::kTHnC, {{3, -10.0f, 10.01f}, {3, -10.0f, 10.01f}, {3, -10.0f, 10.01f}, {3, -10.0f, 10.01f}, {3, -10.0f, 10.01f}, {3, -10.0f, 10.01f}}}}); + + registry.add({"1d-profile", "test 1d profile", {HistType::kTProfile, {{20, 0.0f, 10.01f}}}}); + registry.add({"2d-profile", "test 2d profile", {HistType::kTProfile2D, {{20, 0.0f, 10.01f}, {20, 0.0f, 10.01f}}}}); + registry.add({"3d-profile", "test 3d profile", {HistType::kTProfile3D, {{20, 0.0f, 10.01f}, {20, 0.0f, 10.01f}, {20, 0.0f, 10.01f}}}}); + + registry.add({"2d-weight", "test 2d weight", {HistType::kTH2C, {{2, -10.0f, 10.01f}, {2, -10.0f, 10.01f}}}, true}); + + registry.add({"3d-weight", "test 3d weight", {HistType::kTH3C, {{2, -10.0f, 10.01f}, {2, -10.0f, 10.01f}, {2, -10.0f, 10.01f}}}, true}); + + registry.add({"4d-weight", "test 4d weight", {HistType::kTHnC, {{2, -10.0f, 10.01f}, {2, -10.0f, 10.01f}, {2, -10.0f, 10.01f}, {100, -10.0f, 10.01f}}}, true}); + + registry.add({"1d-profile-weight", "test 1d profile weight", {HistType::kTProfile, {{2, -10.0f, 10.01f}}}, true}); + registry.add({"2d-profile-weight", "test 2d profile weight", {HistType::kTProfile2D, {{2, -10.0f, 10.01f}, {2, -10.0f, 10.01f}}}, true}); + } + + void process(aod::Tracks const& tracks) + { + using namespace aod::track; + // does not work with dynamic columns (e.g. Charge, NormalizedPhi) + registry.fill("1d", tracks, eta > -0.7f); + registry.fill("3d", tracks, eta > 0.f); + registry.fill("5d", tracks, pt > 0.15f); + registry.fill("7d", tracks, pt > 0.15f); + registry.fill("2d-profile", tracks, eta > -0.5f); + + // fill 4d histogram with weight (column X) + registry.fillWeight("4d-weight", tracks, eta > 0.f); + + registry.fillWeight("2d-weight", tracks, eta > 0.f); + + registry.fillWeight("1d-profile-weight", tracks, eta > 0.f); + + for (auto& track : tracks) { + registry.fill("2d", track.eta(), track.pt()); + registry.fill("4d", track.pt(), track.eta(), track.phi(), track.signed1Pt()); + registry.fill("6d", track.pt(), track.eta(), track.phi(), track.snp(), track.tgl(), track.alpha()); + registry.fill("1d-profile", track.pt(), track.eta()); + registry.fill("3d-profile", track.pt(), track.eta(), track.phi(), track.snp()); + + // fill 3d histogram with weight (2.) + registry.fillWeight("3d-weight", track.pt(), track.eta(), track.phi(), 2.); + + registry.fillWeight("2d-profile-weight", track.pt(), track.eta(), track.phi(), 5.); + } + } +}; + +struct DTask { + HistogramRegistry spectra{"spectra", true, {}}; + HistogramRegistry etaStudy{"etaStudy", true, {}}; + + void init(o2::framework::InitContext&) + { + std::vector ptBinning = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, + 1.1, 1.2, 1.3, 1.4, 1.5, 2.0, 5.0, 10.0, 20.0, 50.0}; + std::vector centBinning = {0., 30., 60., 90.}; + + AxisSpec ptAxis = {ptBinning, "#it{p}_{T} (GeV/c)"}; + AxisSpec centAxis = {centBinning, "centrality"}; + AxisSpec etaAxis = {5, -0.8, 0.8, "#eta"}; + AxisSpec phiAxis = {4, 0., 2. * M_PI, "#phi"}; + const int nCuts = 5; + AxisSpec cutAxis = {nCuts, -0.5, nCuts - 0.5, "cut setting"}; + + HistogramConfigSpec defaultParticleHist({HistType::kTHnF, {ptAxis, etaAxis, centAxis, cutAxis}}); + + spectra.add({"myControlHist", "a", {HistType::kTH2F, {ptAxis, etaAxis}}}); + spectra.get("myControlHist")->GetYaxis()->SetTitle("my-y-axis"); + spectra.get("myControlHist")->SetTitle("something meaningful"); + + spectra.add({"pions", "Pions", defaultParticleHist}); + spectra.add({"kaons", "Kaons", defaultParticleHist}); + spectra.add({"sigmas", "Sigmas", defaultParticleHist}); + spectra.add({"lambdas", "Lambd", defaultParticleHist}); + + spectra.get("lambdas")->SetTitle("Lambdas"); + + etaStudy.add({"positive", "A side spectra", {HistType::kTH1I, {ptAxis}}}); + etaStudy.add({"negative", "C side spectra", {HistType::kTH1I, {ptAxis}}}); + } + + void process(aod::Tracks const& tracks) + { + using namespace aod::track; + + etaStudy.fill("positive", tracks, eta > 0.f); + etaStudy.fill("negative", tracks, eta < 0.f); + + for (auto& track : tracks) { + spectra.fill("myControlHist", track.pt(), track.eta()); + spectra.fill("pions", track.pt(), track.eta(), 50., 0.); + spectra.fill("kaons", track.pt(), track.eta(), 50., 0.); + spectra.fill("sigmas", track.pt(), track.eta(), 50., 0.); + spectra.fill("lambdas", track.pt(), track.eta(), 50., 0.); + } + } +}; + WorkflowSpec defineDataProcessing(ConfigContext const&) { return WorkflowSpec{ adaptAnalysisTask("eta-and-phi-histograms"), - adaptAnalysisTask("filtered-histograms")}; + adaptAnalysisTask("filtered-histograms"), + adaptAnalysisTask("dimension-test"), + adaptAnalysisTask("realistic-example")}; } diff --git a/Framework/Core/CMakeLists.txt b/Framework/Core/CMakeLists.txt index d08fae78bf48d..badffa9815841 100644 --- a/Framework/Core/CMakeLists.txt +++ b/Framework/Core/CMakeLists.txt @@ -103,6 +103,7 @@ o2_add_library(Framework src/WorkflowSpec.cxx src/runDataProcessing.cxx src/ExternalFairMQDeviceProxy.cxx + src/HistogramRegistry.cxx test/TestClasses.cxx PRIVATE_INCLUDE_DIRECTORIES ${CMAKE_CURRENT_LIST_DIR}/src PUBLIC_LINK_LIBRARIES ${DEBUGGUI_TARGET} diff --git a/Framework/Core/include/Framework/HistogramRegistry.h b/Framework/Core/include/Framework/HistogramRegistry.h index 93f21ec50479d..b40577ccdf04e 100644 --- a/Framework/Core/include/Framework/HistogramRegistry.h +++ b/Framework/Core/include/Framework/HistogramRegistry.h @@ -23,96 +23,155 @@ #include "Framework/RuntimeError.h" #include "TClass.h" -#include "TH1.h" -#include "TH2.h" -#include "TH3.h" -#include "THn.h" -#include "THnSparse.h" -#include "TFolder.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +//#include +//#include #include #include -namespace o2 -{ -namespace framework +namespace o2::framework { -// Most common histogram types -enum HistogramType : unsigned int { - kTH1D = 0, - kTH1F = 1, - kTH1I = 2, - kTH2D = 3, - kTH2F = 4, - kTH2I = 5, - kTH3D = 6, - kTH3F = 7, - kTH3I = 8 +// Available root histogram types +enum HistType : unsigned int { + kUndefinedHist = 0, + kTH1D, + kTH1F, + kTH1I, + kTH1C, + kTH1S, + kTH2D, + kTH2F, + kTH2I, + kTH2C, + kTH2S, + kTH3D, + kTH3F, + kTH3I, + kTH3C, + kTH3S, + kTHnD, + kTHnF, + kTHnI, + kTHnC, + kTHnS, + kTHnL, + kTHnSparseD, + kTHnSparseF, + kTHnSparseI, + kTHnSparseC, + kTHnSparseS, + kTHnSparseL, + kTProfile, + kTProfile2D, + kTProfile3D, + kStepTHnF, // FIXME: for these two to work we need to align StepTHn ctors with the root THn ones + kStepTHnD }; -/// Description of a single histogram axis +// variant of all possible root pointers; here we use only the interface types since the underlying data representation (int,float,double,long,char) is irrelevant +using HistPtr = std::variant, std::shared_ptr, std::shared_ptr, std::shared_ptr, std::shared_ptr, std::shared_ptr, std::shared_ptr, std::shared_ptr>; + +//************************************************************************************************** +/** + * Specification of an Axis. + */ +//************************************************************************************************** struct AxisSpec { - AxisSpec(int nBins_, std::vector bins_, std::string label_ = "") - : nBins(nBins_), - bins(bins_), - binsEqual(false), - label(label_) + AxisSpec(std::vector binEdges_, std::string title_ = "", std::optional name_ = std::nullopt) + : nBins(std::nullopt), + binEdges(binEdges_), + title(title_), + name(name_) { } - AxisSpec(int nBins_, double binMin_, double binMax_, std::string label_ = "") + AxisSpec(int nBins_, double binMin_, double binMax_, std::string title_ = "", std::optional name_ = std::nullopt) : nBins(nBins_), - bins({binMin_, binMax_}), - binsEqual(true), - label(label_) + binEdges({binMin_, binMax_}), + title(title_), + name(name_) { } - AxisSpec() : nBins(1), binsEqual(false), bins(), label("") {} - - int nBins; - std::vector bins; - bool binsEqual; // if true, then bins specify min and max for equidistant binning - std::string label; + std::optional nBins{}; + std::vector binEdges{}; + std::string title{}; + std::optional name{}; // optional axis name for ndim histograms }; -/// Data sctructure that will allow to construct a fully qualified TH* histogram +//************************************************************************************************** +/** + * Specification of a histogram configuration. + */ +//************************************************************************************************** struct HistogramConfigSpec { - HistogramConfigSpec(HistogramType type_, std::vector axes_) + HistogramConfigSpec(HistType type_, std::vector axes_) : type(type_), - axes(axes_), - binsEqual(axes.size() > 0 ? axes[0].binsEqual : false) + axes(axes_) { } + HistogramConfigSpec() = default; + HistogramConfigSpec(HistogramConfigSpec const& other) = default; + HistogramConfigSpec(HistogramConfigSpec&& other) = default; - HistogramConfigSpec() - : type(HistogramType::kTH1F), - axes(), - binsEqual(false) + void addAxis(const AxisSpec& axis) { + axes.push_back(axis); + } + + void addAxis(int nBins_, const double binMin_, const double binMax_, const std::string& title_ = "", const std::optional& name_ = std::nullopt) + { + axes.push_back({nBins_, binMin_, binMax_, title_, name_}); + } + + void addAxis(const std::vector& binEdges_, const std::string& title_, const std::string& name_) + { + axes.push_back({binEdges_, title_, name_}); + } + + void addAxes(const std::vector& axes_) + { + axes.insert(axes.end(), axes_.begin(), axes_.end()); + } + + // add axes defined in other HistogramConfigSpec object + void addAxes(const HistogramConfigSpec& other) + { + axes.insert(axes.end(), other.axes.begin(), other.axes.end()); } - HistogramConfigSpec(HistogramConfigSpec const& other) = default; - HistogramConfigSpec(HistogramConfigSpec&& other) = default; - HistogramType type; - std::vector axes; - bool binsEqual; + HistType type{HistType::kUndefinedHist}; + std::vector axes{}; }; -/// Data structure containing histogram specification for the HistogramRegistry -/// Contains hashed name as id for fast lookup +//************************************************************************************************** +/** + * Specification of a histogram. + */ +//************************************************************************************************** struct HistogramSpec { - HistogramSpec(char const* const name_, char const* const readableName_, HistogramConfigSpec config_) + HistogramSpec(char const* const name_, char const* const title_, HistogramConfigSpec config_, bool callSumw2_ = false) : name(name_), - readableName(readableName_), id(compile_time_hash(name_)), - config(config_) + title(title_), + config(config_), + callSumw2(callSumw2_) { } HistogramSpec() : name(""), - readableName(""), id(0), config() { @@ -120,191 +179,265 @@ struct HistogramSpec { HistogramSpec(HistogramSpec const& other) = default; HistogramSpec(HistogramSpec&& other) = default; - std::string name; - std::string readableName; - uint32_t id; - HistogramConfigSpec config; + std::string name{}; + uint32_t id{}; + std::string title{}; + HistogramConfigSpec config{}; + bool callSumw2{}; // wether or not hist needs heavy error structure produced by Sumw2() }; -template -std::unique_ptr createTH1FromSpec(HistogramSpec const& spec) -{ - if (spec.config.axes.size() == 0) { - throw runtime_error("No arguments available in spec to create a histogram"); - } - if (spec.config.binsEqual) { - return std::make_unique(spec.name.data(), spec.readableName.data(), spec.config.axes[0].nBins, spec.config.axes[0].bins[0], spec.config.axes[0].bins[1]); - } - return std::make_unique(spec.name.data(), spec.readableName.data(), spec.config.axes[0].nBins, spec.config.axes[0].bins.data()); -} +//************************************************************************************************** +/** + * Static helper class to generate histograms from the specifications. + * Also provides functions to obtain pointer to the created histogram casted to the correct alternative of the std::variant HistPtr that is used in HistogramRegistry. + */ +//************************************************************************************************** +struct HistFactory { + + // create histogram of type T with the axes defined in HistogramConfigSpec + template + static std::shared_ptr createHist(const HistogramSpec& histSpec) + { + constexpr std::size_t MAX_DIM{10}; + const std::size_t nAxes{histSpec.config.axes.size()}; + if (nAxes == 0 || nAxes > MAX_DIM) { + LOGF(FATAL, "The histogram specification contains no (or too many) axes."); + return nullptr; + } -template -std::unique_ptr createTH2FromSpec(HistogramSpec const& spec) -{ - if (spec.config.axes.size() == 0) { - throw runtime_error("No arguments available in spec to create a histogram"); - } - if (spec.config.binsEqual) { - return std::make_unique(spec.name.data(), spec.readableName.data(), spec.config.axes[0].nBins, spec.config.axes[0].bins[0], spec.config.axes[0].bins[1], spec.config.axes[1].nBins, spec.config.axes[1].bins[0], spec.config.axes[1].bins[1]); - } - return std::make_unique(spec.name.data(), spec.readableName.data(), spec.config.axes[0].nBins, spec.config.axes[0].bins.data(), spec.config.axes[1].nBins, spec.config.axes[1].bins.data()); -} + int nBins[MAX_DIM]{0}; + double lowerBounds[MAX_DIM]{0.}; + double upperBounds[MAX_DIM]{0.}; -template -std::unique_ptr createTH3FromSpec(HistogramSpec const& spec) -{ - if (spec.config.axes.size() == 0) { - throw runtime_error("No arguments available in spec to create a histogram"); - } - if (spec.config.binsEqual) { - return std::make_unique(spec.name.data(), spec.readableName.data(), spec.config.axes[0].nBins, spec.config.axes[0].bins[0], spec.config.axes[0].bins[1], spec.config.axes[1].nBins, spec.config.axes[1].bins[0], spec.config.axes[1].bins[1], spec.config.axes[2].nBins, spec.config.axes[2].bins[0], spec.config.axes[2].bins[1]); - } - return std::make_unique(spec.name.data(), spec.readableName.data(), spec.config.axes[0].nBins, spec.config.axes[0].bins.data(), spec.config.axes[1].nBins, spec.config.axes[1].bins.data(), spec.config.axes[2].nBins, spec.config.axes[2].bins.data()); -} + // first figure out number of bins and dimensions + for (std::size_t i = 0; i < nAxes; i++) { + nBins[i] = (histSpec.config.axes[i].nBins) ? *histSpec.config.axes[i].nBins : histSpec.config.axes[i].binEdges.size() - 1; + lowerBounds[i] = histSpec.config.axes[i].binEdges.front(); + upperBounds[i] = histSpec.config.axes[i].binEdges.back(); + } -/// Helper functions to fill histograms with expressions -template -void fill(TH1* hist, const T& table, const o2::framework::expressions::Filter& filter) -{ - auto filtered = o2::soa::Filtered{{table.asArrowTable()}, o2::framework::expressions::createSelection(table.asArrowTable(), filter)}; - for (auto& t : filtered) { - hist->Fill(*(static_cast(t).getIterator()), *(static_cast(t).getIterator()), *(static_cast(t).getIterator())); - } -} + // create histogram + std::shared_ptr hist{generateHist(histSpec.name, histSpec.title, nAxes, nBins, lowerBounds, upperBounds)}; + if (!hist) { + LOGF(FATAL, "The number of specified dimensions does not match the type."); + return nullptr; + } -template -void fill(TH1* hist, const T& table, const o2::framework::expressions::Filter& filter) -{ - auto filtered = o2::soa::Filtered{{table.asArrowTable()}, o2::framework::expressions::createSelection(table.asArrowTable(), filter)}; - for (auto& t : filtered) { - hist->Fill(*(static_cast(t).getIterator()), *(static_cast(t).getIterator())); - } -} + // set axis properties + for (std::size_t i = 0; i < nAxes; i++) { + TAxis* axis{getAxis(i, hist)}; + if (axis) { + axis->SetTitle(histSpec.config.axes[i].title.data()); + + // this helps to have axes not only called 0,1,2... in ndim histos + if constexpr (std::is_base_of_v) { + if (histSpec.config.axes[i].name) + axis->SetName((std::string(axis->GetName()) + "-" + *histSpec.config.axes[i].name).data()); + } + + // move the bin edges in case a variable binning was requested + if (!histSpec.config.axes[i].nBins) { + if (!std::is_sorted(std::begin(histSpec.config.axes[i].binEdges), std::end(histSpec.config.axes[i].binEdges))) { + LOGF(FATAL, "The bin edges specified for axis %s in histogram %s are not in increasing order!", (histSpec.config.axes[i].name) ? *histSpec.config.axes[i].name : histSpec.config.axes[i].title, histSpec.name); + return nullptr; + } + axis->Set(nBins[i], histSpec.config.axes[i].binEdges.data()); + } + } + } + if (histSpec.callSumw2) + hist->Sumw2(); -template -void fill(TH1* hist, const T& table, const o2::framework::expressions::Filter& filter) -{ - auto filtered = o2::soa::Filtered{{table.asArrowTable()}, o2::framework::expressions::createSelection(table.asArrowTable(), filter)}; - for (auto& t : filtered) { - hist->Fill(*(static_cast(t).getIterator())); + return hist; } -} -using HistogramCreationCallback = std::function(HistogramSpec const& spec)>; - -// Wrapper to avoid multiple function definitinions error -struct HistogramCallbacks { - static HistogramCreationCallback createTH1D() + // create histogram and return it casted to the correct alternative held in HistPtr variant + template + static HistPtr createHistVariant(HistogramSpec const& histSpec) { - return [](HistogramSpec const& spec) -> std::unique_ptr { - return createTH1FromSpec(spec); - }; + if (auto hist = castToVariant(createHist(histSpec))) + return *hist; + else + throw runtime_error("Histogram was not created properly."); } - static HistogramCreationCallback createTH1F() + // runtime version of the above + static HistPtr createHistVariant(HistogramSpec const& histSpec) { - return [](HistogramSpec const& spec) -> std::unique_ptr { - return createTH1FromSpec(spec); - }; + if (histSpec.config.type == HistType::kUndefinedHist) + throw runtime_error("Histogram type was not specified."); + else + return HistogramCreationCallbacks.at(histSpec.config.type)(histSpec); } - static HistogramCreationCallback createTH1I() + private: + static const std::map> HistogramCreationCallbacks; + + // helper function to generate the actual histograms + template + static T* generateHist(const std::string& name, const std::string& title, const std::size_t nDim, + const int nBins[], const double lowerBounds[], const double upperBounds[]) { - return [](HistogramSpec const& spec) -> std::unique_ptr { - return createTH1FromSpec(spec); - }; + if constexpr (std::is_base_of_v) { + return new T(name.data(), title.data(), nDim, nBins, lowerBounds, upperBounds); + } else if constexpr (std::is_base_of_v) { + return (nDim != 3) ? nullptr + : new T(name.data(), title.data(), nBins[0], lowerBounds[0], + upperBounds[0], nBins[1], lowerBounds[1], upperBounds[1], + nBins[2], lowerBounds[2], upperBounds[2]); + } else if constexpr (std::is_base_of_v) { + return (nDim != 2) ? nullptr + : new T(name.data(), title.data(), nBins[0], lowerBounds[0], + upperBounds[0], nBins[1], lowerBounds[1], upperBounds[1]); + } else if constexpr (std::is_base_of_v) { + return (nDim != 1) + ? nullptr + : new T(name.data(), title.data(), nBins[0], lowerBounds[0], upperBounds[0]); + } + return nullptr; } - static HistogramCreationCallback createTH2D() + // helper function to get the axis via index for any type of root histogram + template + static TAxis* getAxis(const int i, std::shared_ptr& hist) { - return [](HistogramSpec const& spec) -> std::unique_ptr { - return createTH2FromSpec(spec); - }; + if constexpr (std::is_base_of_v) { + return hist->GetAxis(i); + } else { + return (i == 0) ? hist->GetXaxis() + : (i == 1) ? hist->GetYaxis() : (i == 2) ? hist->GetZaxis() : nullptr; + } } - static HistogramCreationCallback createTH2F() + // helper function to cast the actual histogram type (e.g. TH2F) to the correct interface type (e.g. TH2) that is stored in the HistPtr variant + template + static std::optional castToVariant(std::shared_ptr obj) { - return [](HistogramSpec const& spec) -> std::unique_ptr { - return createTH2FromSpec(spec); - }; + if (obj->InheritsFrom(T::Class())) { + return std::static_pointer_cast(obj); + } + return std::nullopt; } - - static HistogramCreationCallback createTH2I() + template + static std::optional castToVariant(std::shared_ptr obj) { - return [](HistogramSpec const& spec) -> std::unique_ptr { - return createTH2FromSpec(spec); - }; + if (auto hist = castToVariant(obj)) { + return hist; + } + return castToVariant(obj); } - - static HistogramCreationCallback createTH3D() + static std::optional castToVariant(std::shared_ptr obj) { - return [](HistogramSpec const& spec) -> std::unique_ptr { - return createTH3FromSpec(spec); - }; + if (obj) { + // TProfile3D is TH3, TProfile2D is TH2, TH3 is TH1, TH2 is TH1, TProfile is TH1 + return castToVariant(obj); + } + return std::nullopt; } +}; - static HistogramCreationCallback createTH3F() +//************************************************************************************************** +/** + * Static helper class to fill existing root histograms of any type. + * Contains functionality to fill once per call or a whole (filtered) table at once. + */ +//************************************************************************************************** +struct HistFiller { + // fill any type of histogram (if weight was requested it must be the last argument) + template + static void fillHistAny(std::shared_ptr& hist, const Ts&... positionAndWeight) { - return [](HistogramSpec const& spec) -> std::unique_ptr { - return createTH3FromSpec(spec); - }; + constexpr int nDim = sizeof...(Ts) - useWeight; + + constexpr bool validTH3 = (std::is_same_v && nDim == 3); + constexpr bool validTH2 = (std::is_same_v && nDim == 2); + constexpr bool validTH1 = (std::is_same_v && nDim == 1); + constexpr bool validTProfile3D = (std::is_same_v && nDim == 4); + constexpr bool validTProfile2D = (std::is_same_v && nDim == 3); + constexpr bool validTProfile = (std::is_same_v && nDim == 2); + + constexpr bool validSimpleFill = validTH1 || validTH2 || validTH3 || validTProfile || validTProfile2D || validTProfile3D; + // unfortunately we dont know at compile the dimension of THn(Sparse) + constexpr bool validComplexFill = std::is_base_of_v; + + if constexpr (validSimpleFill) { + hist->Fill(static_cast(positionAndWeight)...); + } else if constexpr (validComplexFill) { + // savety check for n dimensional histograms (runtime overhead) + if (hist->GetNdimensions() != nDim) + throw runtime_error("The number of position (and weight) arguments does not match the histogram dimensions!"); + + double tempArray[sizeof...(Ts)] = {static_cast(positionAndWeight)...}; + if constexpr (useWeight) + hist->Fill(tempArray, tempArray[sizeof...(Ts) - 1]); + else + hist->Fill(tempArray); + } else { + throw runtime_error("The number of position (and weight) arguments does not match the histogram dimensions!"); + } } - static HistogramCreationCallback createTH3I() + // fill any type of histogram with columns (Cs) of a filtered table (if weight is requested it must reside the last specified column) + template + static void fillHistAnyTable(std::shared_ptr& hist, const T& table, const o2::framework::expressions::Filter& filter) { - return [](HistogramSpec const& spec) -> std::unique_ptr { - return createTH3FromSpec(spec); - }; + auto filtered = o2::soa::Filtered{{table.asArrowTable()}, o2::framework::expressions::createSelection(table.asArrowTable(), filter)}; + for (auto& t : filtered) { + fillHistAny(hist, (*(static_cast(t).getIterator()))...); + } } }; -/// Histogram registry for an analysis task that allows to define needed histograms -/// and serves as the container/wrapper to fill them +//************************************************************************************************** +/** + * Histogram registry that can be used to store and fill histograms of any type. + */ +//************************************************************************************************** class HistogramRegistry { public: - HistogramRegistry(char const* const name_, bool enable, std::vector specs, OutputObjHandlingPolicy policy_ = OutputObjHandlingPolicy::AnalysisObject) - : name(name_), - policy(policy_), - enabled(enable), - mRegistryKey(), - mRegistryValue(), - mHistogramCreationCallbacks({HistogramCallbacks::createTH1D(), - HistogramCallbacks::createTH1F(), - HistogramCallbacks::createTH1I(), - HistogramCallbacks::createTH2D(), - HistogramCallbacks::createTH2F(), - HistogramCallbacks::createTH2I(), - HistogramCallbacks::createTH3D(), - HistogramCallbacks::createTH3F(), - HistogramCallbacks::createTH3I()}) + HistogramRegistry(char const* const name_, bool enable, std::vector histSpecs_, OutputObjHandlingPolicy policy_ = OutputObjHandlingPolicy::AnalysisObject) : name(name_), + policy(policy_), + enabled(enable), + mRegistryKey(), + mRegistryValue() { mRegistryKey.fill(0u); - for (auto& spec : specs) { - insert(spec); + for (auto& histSpec : histSpecs_) { + insert(histSpec); } } + void add(HistogramSpec&& histSpec_) + { + insert(histSpec_); + } + + // gets the underlying histogram pointer + // we cannot automatically infer type here so it has to be explicitly specified + // -> get(), get(), get(), get(), get(), get(), get(), get() /// @return the histogram registered with name @a name - auto& get(char const* const name) const + template + auto& get(char const* const name) { const uint32_t id = compile_time_hash(name); const uint32_t i = imask(id); if (O2_BUILTIN_LIKELY(id == mRegistryKey[i])) { - return mRegistryValue[i]; + return *std::get_if>(&mRegistryValue[i]); } for (auto j = 1u; j < MAX_REGISTRY_SIZE; ++j) { if (id == mRegistryKey[imask(j + i)]) { - return mRegistryValue[imask(j + i)]; + return *std::get_if>(&mRegistryValue[imask(j + i)]); } } - throw runtime_error("No match found!"); + throw runtime_error("No matching histogram found in HistogramRegistry!"); } /// @return the histogram registered with name @a name - auto& operator()(char const* const name) const + template + auto& operator()(char const* const name) { - return get(name); + return get(name); } /// @return the associated OutputSpec @@ -324,79 +457,114 @@ class HistogramRegistry taskHash = hash; } + // TODO: maybe support also TDirectory,TList,TObjArray?, find a way to write to file in 'single key' mode (arg in WriteObjAny) TFolder* operator*() { TFolder* folder = new TFolder(this->name.c_str(), this->name.c_str()); for (auto j = 0u; j < MAX_REGISTRY_SIZE; ++j) { - if (mRegistryValue[j].get() != nullptr) { - auto hist = mRegistryValue[j].get(); - folder->Add(hist); + TObject* rawPtr = nullptr; + std::visit([&](const auto& sharedPtr) { rawPtr = sharedPtr.get(); }, mRegistryValue[j]); + if (rawPtr) { + folder->Add(rawPtr); } } - folder->SetOwner(); + folder->SetOwner(false); // object deletion will be handled by shared_ptrs return folder; } - /// fill the histogram with an expression - template - void fill(char const* const name, const T& table, const o2::framework::expressions::Filter& filter) + // fill hist with values + template + void fill(char const* const name, Ts&&... positionAndWeight) { - TH1* hist = get(name).get(); - framework::fill(hist, table, filter); + const uint32_t id = compile_time_hash(name); + const uint32_t i = imask(id); + if (O2_BUILTIN_LIKELY(id == mRegistryKey[i])) { + std::visit([this, &positionAndWeight...](auto&& hist) { HistFiller::fillHistAny(hist, std::forward(positionAndWeight)...); }, mRegistryValue[i]); + return; + } + for (auto j = 1u; j < MAX_REGISTRY_SIZE; ++j) { + if (id == mRegistryKey[imask(j + i)]) { + std::visit([this, &positionAndWeight...](auto&& hist) { HistFiller::fillHistAny(hist, std::forward(positionAndWeight)...); }, mRegistryValue[imask(j + i)]); + return; + } + } + throw runtime_error("No matching histogram found in HistogramRegistry!"); } - template - void fill(char const* const name, const T& table, const o2::framework::expressions::Filter& filter) + template + void fillWeight(char const* const name, Ts&&... positionAndWeight) { - TH1* hist = get(name).get(); - framework::fill(hist, table, filter); + fill(name, std::forward(positionAndWeight)...); } - template + // fill hist with content of (filtered) table columns + template void fill(char const* const name, const T& table, const o2::framework::expressions::Filter& filter) { - TH1* hist = get(name).get(); - framework::fill(hist, table, filter); + fillTable(name, table, filter); + } + template + void fillWeight(char const* const name, const T& table, const o2::framework::expressions::Filter& filter) + { + fillTable(name, table, filter); + } + + // this is for internal use only because we dont want user to have to specify 'useWeight' argument + template + void fillTable(char const* const name, const T& table, const o2::framework::expressions::Filter& filter) + { + const uint32_t id = compile_time_hash(name); + const uint32_t i = imask(id); + if (O2_BUILTIN_LIKELY(id == mRegistryKey[i])) { + std::visit([this, &table, &filter](auto&& hist) { HistFiller::fillHistAnyTable(hist, table, filter); }, mRegistryValue[i]); + return; + } + for (auto j = 1u; j < MAX_REGISTRY_SIZE; ++j) { + if (id == mRegistryKey[imask(j + i)]) { + std::visit([this, &table, &filter](auto&& hist) { HistFiller::fillHistAnyTable(hist, table, filter); }, mRegistryValue[imask(j + i)]); + return; + } + } + throw runtime_error("No matching histogram found in HistogramRegistry!"); } /// lookup distance counter for benchmarking mutable uint32_t lookup = 0; private: - void insert(HistogramSpec& spec) + void insert(HistogramSpec& histSpec) { - uint32_t i = imask(spec.id); + uint32_t i = imask(histSpec.id); for (auto j = 0u; j < MAX_REGISTRY_SIZE; ++j) { - if (mRegistryValue[imask(j + i)].get() == nullptr) { - mRegistryKey[imask(j + i)] = spec.id; - mRegistryValue[imask(j + i)] = mHistogramCreationCallbacks[spec.config.type](spec); + TObject* rawPtr = nullptr; + std::visit([&](const auto& sharedPtr) { rawPtr = sharedPtr.get(); }, mRegistryValue[imask(j + i)]); + if (!rawPtr) { + mRegistryKey[imask(j + i)] = histSpec.id; + mRegistryValue[imask(j + i)] = HistFactory::createHistVariant(histSpec); lookup += j; return; } } - throw runtime_error("Internal array is full."); + throw runtime_error("Internal array of HistogramRegistry is full."); } inline constexpr uint32_t imask(uint32_t i) const { return i & mask; } - std::string name; - bool enabled; - OutputObjHandlingPolicy policy; - uint32_t taskHash; + + std::string name{}; + bool enabled{}; + OutputObjHandlingPolicy policy{}; + uint32_t taskHash{}; /// The maximum number of histograms in buffer is currently set to 512 /// which seems to be both reasonably large and allowing for very fast lookup static constexpr uint32_t mask = 0x1FF; static constexpr uint32_t MAX_REGISTRY_SIZE = mask + 1; - std::array mRegistryKey; - std::array, MAX_REGISTRY_SIZE> mRegistryValue; - std::vector mHistogramCreationCallbacks; + std::array mRegistryKey{}; + std::array mRegistryValue{}; }; -} // namespace framework - -} // namespace o2 - +} // namespace o2::framework #endif // FRAMEWORK_HISTOGRAMREGISTRY_H_ diff --git a/Framework/Core/src/HistogramRegistry.cxx b/Framework/Core/src/HistogramRegistry.cxx new file mode 100644 index 0000000000000..ef336c52ce208 --- /dev/null +++ b/Framework/Core/src/HistogramRegistry.cxx @@ -0,0 +1,37 @@ +// Copyright CERN and copyright holders of ALICE O2. This software is +// distributed under the terms of the GNU General Public License v3 (GPL +// Version 3), copied verbatim in the file "COPYING". +// +// See http://alice-o2.web.cern.ch/license for full licensing information. +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +#include "Framework/HistogramRegistry.h" + +namespace o2::framework +{ + +// define histogram callbacks for runtime histogram creation +#define CALLB(HType) \ + { \ + k##HType, \ + [](HistogramSpec const& histSpec) { \ + return HistFactory::createHistVariant(histSpec); \ + } \ + } + +const std::map> HistFactory::HistogramCreationCallbacks{ + CALLB(TH1D), CALLB(TH1F), CALLB(TH1I), CALLB(TH1C), CALLB(TH1S), + CALLB(TH2D), CALLB(TH2F), CALLB(TH2I), CALLB(TH2C), CALLB(TH2S), + CALLB(TH3D), CALLB(TH3F), CALLB(TH3I), CALLB(TH3C), CALLB(TH3S), + CALLB(THnD), CALLB(THnF), CALLB(THnI), CALLB(THnC), CALLB(THnS), CALLB(THnL), + CALLB(THnSparseD), CALLB(THnSparseF), CALLB(THnSparseI), CALLB(THnSparseC), CALLB(THnSparseS), CALLB(THnSparseL), + CALLB(TProfile), CALLB(TProfile2D), CALLB(TProfile3D), + //CALLB(StepTHnF), CALLB(StepTHnD) +}; + +#undef CALLB + +} // namespace o2::framework diff --git a/Framework/Core/test/benchmark_HistogramRegistry.cxx b/Framework/Core/test/benchmark_HistogramRegistry.cxx index f91436dcc4cdd..05c4d17047ec8 100644 --- a/Framework/Core/test/benchmark_HistogramRegistry.cxx +++ b/Framework/Core/test/benchmark_HistogramRegistry.cxx @@ -28,15 +28,15 @@ static void BM_HashedNameLookup(benchmark::State& state) { for (auto _ : state) { state.PauseTiming(); - std::vector specs; + std::vector histSpecs; for (auto i = 0; i < state.range(0); ++i) { - specs.push_back({fmt::format("histo{}", i + 1).c_str(), fmt::format("Histo {}", i + 1).c_str(), {HistogramType::kTH1F, {{100, 0, 1}}}}); + histSpecs.push_back({fmt::format("histo{}", i + 1).c_str(), fmt::format("Histo {}", i + 1).c_str(), {HistType::kTH1F, {{100, 0, 1}}}}); } - HistogramRegistry registry{"registry", true, specs}; + HistogramRegistry registry{"registry", true, histSpecs}; state.ResumeTiming(); for (auto i = 0; i < nLookups; ++i) { - auto& x = registry.get("histo4"); + auto& x = registry.get("histo4"); benchmark::DoNotOptimize(x); } state.counters["Average lookup distance"] = ((double)registry.lookup / (double)(state.range(0))); diff --git a/Framework/Core/test/test_HistogramRegistry.cxx b/Framework/Core/test/test_HistogramRegistry.cxx index 0ecb7c367f16f..e4bb459a81bdb 100644 --- a/Framework/Core/test/test_HistogramRegistry.cxx +++ b/Framework/Core/test/test_HistogramRegistry.cxx @@ -25,7 +25,7 @@ DECLARE_SOA_COLUMN_FULL(Y, y, float, "y"); HistogramRegistry foo() { - return {"r", true, {{"histo", "histo", {HistogramType::kTH1F, {{100, 0, 1}}}}}}; + return {"r", true, {{"histo", "histo", {HistType::kTH1F, {{100, 0, 1}}}}}}; } BOOST_AUTO_TEST_CASE(HistogramRegistryLookup) @@ -33,27 +33,27 @@ BOOST_AUTO_TEST_CASE(HistogramRegistryLookup) /// Construct a registry object with direct declaration HistogramRegistry registry{ "registry", true, { - {"eta", "#Eta", {HistogramType::kTH1F, {{100, -2.0, 2.0}}}}, // - {"phi", "#Phi", {HistogramType::kTH1D, {{102, 0, 2 * M_PI}}}}, // - {"pt", "p_{T}", {HistogramType::kTH1D, {{1002, -0.01, 50.1}}}}, // - {"ptToPt", "#ptToPt", {HistogramType::kTH2F, {{100, -0.01, 10.01}, {100, -0.01, 10.01}}}} // - } // + {"eta", "#Eta", {HistType::kTH1F, {{100, -2.0, 2.0}}}}, // + {"phi", "#Phi", {HistType::kTH1D, {{102, 0, 2 * M_PI}}}}, // + {"pt", "p_{T}", {HistType::kTH1D, {{1002, -0.01, 50.1}}}}, // + {"ptToPt", "#ptToPt", {HistType::kTH2F, {{100, -0.01, 10.01}, {100, -0.01, 10.01}}}} // + } // }; /// Get histograms by name - BOOST_REQUIRE_EQUAL(registry.get("eta")->GetNbinsX(), 100); - BOOST_REQUIRE_EQUAL(registry.get("phi")->GetNbinsX(), 102); - BOOST_REQUIRE_EQUAL(registry.get("pt")->GetNbinsX(), 1002); - BOOST_REQUIRE_EQUAL(registry.get("ptToPt")->GetNbinsX(), 100); - BOOST_REQUIRE_EQUAL(registry.get("ptToPt")->GetNbinsY(), 100); + BOOST_REQUIRE_EQUAL(registry.get("eta")->GetNbinsX(), 100); + BOOST_REQUIRE_EQUAL(registry.get("phi")->GetNbinsX(), 102); + BOOST_REQUIRE_EQUAL(registry.get("pt")->GetNbinsX(), 1002); + BOOST_REQUIRE_EQUAL(registry.get("ptToPt")->GetNbinsX(), 100); + BOOST_REQUIRE_EQUAL(registry.get("ptToPt")->GetNbinsY(), 100); /// Get a pointer to the histogram - auto histo = registry.get("pt").get(); + auto histo = registry.get("pt").get(); BOOST_REQUIRE_EQUAL(histo->GetNbinsX(), 1002); /// Get registry object from a function auto r = foo(); - auto histo2 = r.get("histo").get(); + auto histo2 = r.get("histo").get(); BOOST_REQUIRE_EQUAL(histo2->GetNbinsX(), 100); } @@ -78,16 +78,16 @@ BOOST_AUTO_TEST_CASE(HistogramRegistryExpressionFill) /// Construct a registry object with direct declaration HistogramRegistry registry{ "registry", true, { - {"x", "test x", {HistogramType::kTH1F, {{100, 0.0f, 10.0f}}}}, // - {"xy", "test xy", {HistogramType::kTH2F, {{100, -10.0f, 10.01f}, {100, -10.0f, 10.01f}}}} // - } // + {"x", "test x", {HistType::kTH1F, {{100, 0.0f, 10.0f}}}}, // + {"xy", "test xy", {HistType::kTH2F, {{100, -10.0f, 10.01f}, {100, -10.0f, 10.01f}}}} // + } // }; /// Fill histogram with expression and table registry.fill("x", tests, test::x > 3.0f); - BOOST_CHECK_EQUAL(registry.get("x")->GetEntries(), 4); + BOOST_CHECK_EQUAL(registry.get("x")->GetEntries(), 4); /// Fill histogram with expression and table registry.fill("xy", tests, test::x > 3.0f && test::y > -5.0f); - BOOST_CHECK_EQUAL(registry.get("xy")->GetEntries(), 2); + BOOST_CHECK_EQUAL(registry.get("xy")->GetEntries(), 2); } From fb22bfcb2b078a72ad9090bb5499e8fb706fd776 Mon Sep 17 00:00:00 2001 From: Mario Kruger Date: Sat, 17 Oct 2020 09:20:29 +0200 Subject: [PATCH 2/2] fix typos in auto test case --- Framework/Core/test/test_HistogramRegistry.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Framework/Core/test/test_HistogramRegistry.cxx b/Framework/Core/test/test_HistogramRegistry.cxx index e4bb459a81bdb..59227a8ba0388 100644 --- a/Framework/Core/test/test_HistogramRegistry.cxx +++ b/Framework/Core/test/test_HistogramRegistry.cxx @@ -44,8 +44,8 @@ BOOST_AUTO_TEST_CASE(HistogramRegistryLookup) BOOST_REQUIRE_EQUAL(registry.get("eta")->GetNbinsX(), 100); BOOST_REQUIRE_EQUAL(registry.get("phi")->GetNbinsX(), 102); BOOST_REQUIRE_EQUAL(registry.get("pt")->GetNbinsX(), 1002); - BOOST_REQUIRE_EQUAL(registry.get("ptToPt")->GetNbinsX(), 100); - BOOST_REQUIRE_EQUAL(registry.get("ptToPt")->GetNbinsY(), 100); + BOOST_REQUIRE_EQUAL(registry.get("ptToPt")->GetNbinsX(), 100); + BOOST_REQUIRE_EQUAL(registry.get("ptToPt")->GetNbinsY(), 100); /// Get a pointer to the histogram auto histo = registry.get("pt").get();