diff --git a/Analysis/Tutorials/src/histogramRegistry.cxx b/Analysis/Tutorials/src/histogramRegistry.cxx index fa79c947733d6..2c74ce051d91d 100644 --- a/Analysis/Tutorials/src/histogramRegistry.cxx +++ b/Analysis/Tutorials/src/histogramRegistry.cxx @@ -91,6 +91,10 @@ struct CTask { 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); + + registry.add("2d-step", "test 2d step", {HistType::kStepTHnD, {{2, -10.0f, 10.01f}, {2, -10.0f, 10.01f}}, 3}); + + registry.add("2d-step-weight", "test 2d step weight", {HistType::kStepTHnF, {{2, -10.0f, 10.01f}, {2, -10.0f, 10.01f}}, 4}, true); } void process(aod::Tracks const& tracks) @@ -121,6 +125,9 @@ struct CTask { registry.fill(HIST("3d-weight"), track.pt(), track.eta(), track.phi(), 2.); registry.fill(HIST("2d-profile-weight"), track.pt(), track.eta(), track.phi(), 5.); + + registry.fill(HIST("2d-step"), 1, track.pt(), track.eta()); + registry.fill(HIST("2d-step-weight"), 2, track.pt(), track.eta(), track.phi()); } } }; diff --git a/Framework/Core/include/Framework/HistogramRegistry.h b/Framework/Core/include/Framework/HistogramRegistry.h index d201c5d6f6d5f..c40e506d44c7b 100644 --- a/Framework/Core/include/Framework/HistogramRegistry.h +++ b/Framework/Core/include/Framework/HistogramRegistry.h @@ -76,12 +76,12 @@ enum HistType : unsigned int { kTProfile, kTProfile2D, kTProfile3D, - kStepTHnF, // FIXME: for these two to work we need to align StepTHn ctors with the root THn ones + kStepTHnF, kStepTHnD }; // 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>; +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, std::shared_ptr>; //************************************************************************************************** /** @@ -117,9 +117,10 @@ struct AxisSpec { */ //************************************************************************************************** struct HistogramConfigSpec { - HistogramConfigSpec(HistType type_, std::vector axes_) + HistogramConfigSpec(HistType type_, std::vector axes_, uint8_t nSteps_ = 1) : type(type_), - axes(axes_) + axes(axes_), + nSteps(nSteps_) { } HistogramConfigSpec() = default; @@ -154,6 +155,7 @@ struct HistogramConfigSpec { HistType type{HistType::kUndefinedHist}; std::vector axes{}; + uint32_t nSteps{1}; // variable used only in StepTHn }; //************************************************************************************************** @@ -218,7 +220,7 @@ struct HistFactory { } // create histogram - std::shared_ptr hist{generateHist(histSpec.name, histSpec.title, nAxes, nBins, lowerBounds, upperBounds)}; + std::shared_ptr hist{generateHist(histSpec.name, histSpec.title, nAxes, nBins, lowerBounds, upperBounds, histSpec.config.nSteps)}; if (!hist) { LOGF(FATAL, "The number of dimensions specified for histogram %s does not match the type.", histSpec.name); return nullptr; @@ -280,7 +282,7 @@ struct HistFactory { template static TAxis* getAxis(const int i, std::shared_ptr& hist) { - if constexpr (std::is_base_of_v) { + if constexpr (std::is_base_of_v || std::is_base_of_v) { return hist->GetAxis(i); } else { return (i == 0) ? hist->GetXaxis() @@ -296,9 +298,11 @@ struct HistFactory { // 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[]) + const int nBins[], const double lowerBounds[], const double upperBounds[], const int nSteps = 1) { - if constexpr (std::is_base_of_v) { + if constexpr (std::is_base_of_v) { + return new T(name.data(), title.data(), nSteps, nDim, nBins, lowerBounds, upperBounds); + } else 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 @@ -340,7 +344,7 @@ struct HistFactory { { if (obj) { // TProfile3D is TH3, TProfile2D is TH2, TH3 is TH1, TH2 is TH1, TProfile is TH1 - return castToVariant(obj); + return castToVariant(obj); } return std::nullopt; } @@ -369,9 +373,12 @@ struct HistFiller { 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; + constexpr bool validComplexFillStep = std::is_base_of_v; if constexpr (validSimpleFill) { hist->Fill(static_cast(positionAndWeight)...); + } else if constexpr (validComplexFillStep) { + hist->Fill(positionAndWeight...); // first argument in pack is iStep, dimension check is done in StepTHn itself } else if constexpr (validComplexFill) { double tempArray[] = {static_cast(positionAndWeight)...}; double weight{1.}; @@ -391,6 +398,10 @@ struct HistFiller { template static void fillHistAny(std::shared_ptr& hist, const T& table, const o2::framework::expressions::Filter& filter) { + if constexpr (std::is_base_of_v) { + LOGF(FATAL, "Table filling is not (yet?) supported for StepTHn."); + return; + } auto filtered = o2::soa::Filtered{{table.asArrowTable()}, o2::framework::expressions::createSelection(table.asArrowTable(), filter)}; for (auto& t : filtered) { fillHistAny(hist, (*(static_cast(t).getIterator()))...); @@ -615,7 +626,7 @@ class HistogramRegistry template void insertClone(const HistName& histName, const std::shared_ptr& originalHist) { - validateHash(histName.hash, histName.str); + validateHistName(histName.str, histName.hash); for (auto i = 0u; i < MAX_REGISTRY_SIZE; ++i) { TObject* rawPtr = nullptr; std::visit([&](const auto& sharedPtr) { rawPtr = sharedPtr.get(); }, mRegistryValue[imask(histName.idx + i)]); @@ -627,10 +638,11 @@ class HistogramRegistry return; } } - LOGF(FATAL, "Internal array of HistogramRegistry %s is full.", mName); + LOGF(FATAL, R"(Internal array of HistogramRegistry "%s" is full.)", mName); } - void validateHash(const uint32_t hash, const char* name); + // helper function that checks if histogram name can be used in registry + void validateHistName(const char* name, const uint32_t hash); // helper function to find the histogram position in the registry template diff --git a/Framework/Core/include/Framework/StepTHn.h b/Framework/Core/include/Framework/StepTHn.h index 33a06d260ef77..71f1288b06533 100644 --- a/Framework/Core/include/Framework/StepTHn.h +++ b/Framework/Core/include/Framework/StepTHn.h @@ -30,7 +30,7 @@ class StepTHn : public TNamed { public: StepTHn(); - StepTHn(const Char_t* name, const Char_t* title, const Int_t nSteps, const Int_t nAxis, Int_t* nBins, std::vector binLimits[], const char** axisTitles); + StepTHn(const Char_t* name, const Char_t* title, const Int_t nSteps, const Int_t nAxes); ~StepTHn() override; template @@ -55,6 +55,9 @@ class StepTHn : public TNamed virtual Long64_t Merge(TCollection* list) = 0; + TAxis* GetAxis(int i) { return mPrototype->GetAxis(i); } + void Sumw2(){}; // TODO: added for compatibiltiy with registry, but maybe it would be useful also in StepTHn as toggle for error weights + protected: void init(); virtual TArray* createArray(const TArray* src = nullptr) const = 0; @@ -86,7 +89,8 @@ class StepTHnT : public StepTHn { public: StepTHnT() : StepTHn() {} - StepTHnT(const Char_t* name, const Char_t* title, const Int_t nSteps, const Int_t nAxis, Int_t* nBins, std::vector binLimits[], const char** axisTitles); + StepTHnT(const Char_t* name, const Char_t* title, const Int_t nSteps, const Int_t nAxes, Int_t* nBins, std::vector binLimits[], const char** axisTitles); + StepTHnT(const char* name, const char* title, const int nSteps, const int nAxes, const int* nBins, const double* xmin, const double* xmax); ~StepTHnT() override = default; protected: @@ -110,6 +114,10 @@ typedef StepTHnT StepTHnD; template void StepTHn::Fill(Int_t istep, const Ts&... valuesAndWeight) { + if (istep >= mNSteps) { + LOGF(FATAL, "Selected step for filling is not in range of StepTHn."); + } + constexpr int nParams = sizeof...(Ts); // TODO Find a way to avoid the array double tempArray[nParams] = {static_cast(valuesAndWeight)...}; @@ -118,7 +126,7 @@ void StepTHn::Fill(Int_t istep, const Ts&... valuesAndWeight) if (nParams == mNVars + 1) { weight = tempArray[mNVars]; } else if (nParams != mNVars) { - LOGF(fatal, "Fill called with invalid number of parameters (%d vs %d)", mNVars, nParams); + LOGF(FATAL, "Fill called with invalid number of parameters (%d vs %d)", mNVars, nParams); } // fill axis cache @@ -170,7 +178,7 @@ void StepTHn::Fill(Int_t istep, const Ts&... valuesAndWeight) LOGF(info, "Created values container for step %d", istep); } - if (weight != 1) { + if (weight != 1.) { // initialize with already filled entries (which have been filled with weight == 1), in this case mSumw2 := mValues if (!mSumw2[istep]) { mSumw2[istep] = createArray(); diff --git a/Framework/Core/src/HistogramRegistry.cxx b/Framework/Core/src/HistogramRegistry.cxx index bb19012367fe3..f833f2ec8f9ce 100644 --- a/Framework/Core/src/HistogramRegistry.cxx +++ b/Framework/Core/src/HistogramRegistry.cxx @@ -9,6 +9,7 @@ // or submit itself to any jurisdiction. #include "Framework/HistogramRegistry.h" +#include namespace o2::framework { @@ -29,15 +30,14 @@ const std::map> HistFacto 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) -}; + CALLB(StepTHnF), CALLB(StepTHnD)}; #undef CALLB // create histogram from specification and insert it into the registry void HistogramRegistry::insert(const HistogramSpec& histSpec) { - validateHash(histSpec.hash, histSpec.name.data()); + validateHistName(histSpec.name.data(), histSpec.hash); const uint32_t idx = imask(histSpec.hash); for (auto i = 0u; i < MAX_REGISTRY_SIZE; ++i) { TObject* rawPtr = nullptr; @@ -53,8 +53,10 @@ void HistogramRegistry::insert(const HistogramSpec& histSpec) LOGF(FATAL, R"(Internal array of HistogramRegistry "%s" is full.)", mName); } -void HistogramRegistry::validateHash(const uint32_t hash, const char* name) +// helper function that checks if histogram name can be used in registry +void HistogramRegistry::validateHistName(const char* name, const uint32_t hash) { + // validate that hash is unique auto it = std::find(mRegistryKey.begin(), mRegistryKey.end(), hash); if (it != mRegistryKey.end()) { auto idx = it - mRegistryKey.begin(); @@ -62,6 +64,13 @@ void HistogramRegistry::validateHash(const uint32_t hash, const char* name) std::visit([&](const auto& hist) { collidingName = hist->GetName(); }, mRegistryValue[idx]); LOGF(FATAL, R"(Hash collision in HistogramRegistry "%s"! Please rename histogram "%s" or "%s".)", mName, name, collidingName); } + // validate that name contains only allowed characters + /* + TODO: exactly determine constraints for valid histogram names + if (!std::regex_match(name, std::regex("[A-Za-z0-9\-_/]+"))) { + LOGF(FATAL, R"(Histogram name "%s" contains invalid characters.)", name); + } + */ } void HistogramRegistry::add(const HistogramSpec& histSpec) diff --git a/Framework/Core/src/StepTHn.cxx b/Framework/Core/src/StepTHn.cxx index 9bc3fe8391138..9c61ec920b18a 100644 --- a/Framework/Core/src/StepTHn.cxx +++ b/Framework/Core/src/StepTHn.cxx @@ -25,60 +25,68 @@ #include "THn.h" #include "TMath.h" -ClassImp(StepTHn) - templateClassImp(StepTHnT) - - StepTHn::StepTHn() : mNBins(0), - mNVars(0), - mNSteps(0), - mValues(nullptr), - mSumw2(nullptr), - mTarget(nullptr), - mAxisCache(nullptr), - mNbinsCache(nullptr), - mLastVars(nullptr), - mLastBins(nullptr), - mPrototype(nullptr) +ClassImp(StepTHn); +templateClassImp(StepTHnT); + +StepTHn::StepTHn() : mNBins(0), + mNVars(0), + mNSteps(0), + mValues(nullptr), + mSumw2(nullptr), + mTarget(nullptr), + mAxisCache(nullptr), + mNbinsCache(nullptr), + mLastVars(nullptr), + mLastBins(nullptr), + mPrototype(nullptr) { // Default constructor (for streaming) } -StepTHn::StepTHn(const Char_t* name, const Char_t* title, const Int_t nSteps, const Int_t nAxis, Int_t* nBins, std::vector binEdges[], const char** axisTitles) : TNamed(name, title), - mNBins(0), - mNVars(nAxis), - mNSteps(nSteps), - mValues(nullptr), - mSumw2(nullptr), - mTarget(nullptr), - mAxisCache(nullptr), - mNbinsCache(nullptr), - mLastVars(nullptr), - mLastBins(nullptr), - mPrototype(nullptr) +StepTHn::StepTHn(const Char_t* name, const Char_t* title, const Int_t nSteps, const Int_t nAxes) : TNamed(name, title), + mNBins(0), + mNVars(nAxes), + mNSteps(nSteps), + mValues(nullptr), + mSumw2(nullptr), + mTarget(nullptr), + mAxisCache(nullptr), + mNbinsCache(nullptr), + mLastVars(nullptr), + mLastBins(nullptr), + mPrototype(nullptr) { // Constructor // // This will create a container for steps. The memory for such a step is only allocated once the first value is filled. // Therefore you can easily create many steps which are only filled under certain analysis settings. - // For each step a dimensional histogram is created. + // For each step a dimensional histogram is created. // The axis have bins. The bin edges are given in . If there are only two bin edges, equidistant binning is set. init(); } +// root-like constructor template -StepTHnT::StepTHnT(const Char_t* name, const Char_t* title, const Int_t nSteps, const Int_t nAxis, - Int_t* nBins, std::vector binEdges[], const char** axisTitles) : StepTHn(name, title, nSteps, nAxis, nBins, binEdges, axisTitles) +StepTHnT::StepTHnT(const char* name, const char* title, const int nSteps, const int nAxes, const int* nBins, const double* xmin, const double* xmax) : StepTHn(name, title, nSteps, nAxes) { mNBins = 1; for (Int_t i = 0; i < mNVars; i++) { mNBins *= nBins[i]; } - if constexpr (std::is_same_v) { - mPrototype = new THnSparseD(Form("%s_sparse", name), title, nAxis, nBins); - } else { - mPrototype = new THnSparseF(Form("%s_sparse", name), title, nAxis, nBins); + mPrototype = new THnSparseT(Form("%s_sparse", name), title, nAxes, nBins, xmin, xmax); +} + +template +StepTHnT::StepTHnT(const Char_t* name, const Char_t* title, const Int_t nSteps, const Int_t nAxes, + Int_t* nBins, std::vector binEdges[], const char** axisTitles) : StepTHn(name, title, nSteps, nAxes) +{ + mNBins = 1; + for (Int_t i = 0; i < mNVars; i++) { + mNBins *= nBins[i]; } + mPrototype = new THnSparseT(Form("%s_sparse", name), title, nAxes, nBins); + for (Int_t i = 0; i < mNVars; i++) { if (nBins[i] + 1 == binEdges[i].size()) { // variable-width binning mPrototype->GetAxis(i)->Set(nBins[i], &(binEdges[i])[0]); @@ -202,15 +210,14 @@ void StepTHn::Copy(TObject& c) const } if (mPrototype) { - target.mPrototype = dynamic_cast(mPrototype->Clone()); + target.mPrototype = dynamic_cast(mPrototype->Clone()); } } template Long64_t StepTHnT::Merge(TCollection* list) { - // Merge a list of StepTHn objects with this (needed for - // PROOF). + // Merge a list of StepTHn objects with this (needed for PROOF). // Returns the number of merged objects (including this). if (!list) { diff --git a/Framework/Core/test/test_HistogramRegistry.cxx b/Framework/Core/test/test_HistogramRegistry.cxx index d59587c444ad3..2e0234ff7c1c4 100644 --- a/Framework/Core/test/test_HistogramRegistry.cxx +++ b/Framework/Core/test/test_HistogramRegistry.cxx @@ -107,3 +107,24 @@ BOOST_AUTO_TEST_CASE(HistogramRegistryExpressionFill) registry.fill(HIST("xy"), tests, test::x > 3.0f && test::y > -5.0f); BOOST_CHECK_EQUAL(registry.get(HIST("xy"))->GetEntries(), 2); } + +BOOST_AUTO_TEST_CASE(HistogramRegistryStepTHn) +{ + HistogramRegistry registry{"registry"}; + + registry.add("stepTHnF", "a", {kStepTHnF, {{100, -10.0f, 10.01f}, {100, -10.0f, 10.01f}}, 2}); + registry.add("stepTHnD", "b", {kStepTHnD, {{100, -10.0f, 10.01f}, {100, -10.0f, 10.01f}}, 3}); + registry.addClone("stepTHnD", "stepTHnD2"); + + auto& histo = registry.get(HIST("stepTHnD")); + BOOST_REQUIRE_EQUAL(histo->getNSteps(), 3); + + // fill first step at position (0,3) + registry.fill(HIST("stepTHnF"), 0, 0., 3.); + // fill second step (0,4) + registry.fill(HIST("stepTHnF"), 1, 0., 4.); + + registry.fill(HIST("stepTHnD2"), 1, 0., 4.); + + registry.print(); +}