Skip to content

DPL Analysis: make HistogramRegistry support StepTHn #5266

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jgrosseo merged 5 commits into from
Jan 27, 2021
Merged

DPL Analysis: make HistogramRegistry support StepTHn #5266

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
4b0ddd0
DPL Analysis: add sanity check for histogram names before inserting i…
mario-krueger Dec 16, 2020
964fe3d
DPL Analysis: make HistogramRegistry support StepTHn
mario-krueger Jan 22, 2021
9fc9aa8
properly set name in Clone()
mario-krueger Jan 23, 2021
e991033
change order of ctor arguments
mario-krueger Jan 26, 2021
492beb7
remove unnecessary clone
mario-krueger Jan 27, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 7 additions & 0 deletions Analysis/Tutorials/src/histogramRegistry.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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());
}
}
};
Expand Down
36 changes: 24 additions & 12 deletions Framework/Core/include/Framework/HistogramRegistry.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<THn>, std::shared_ptr<THnSparse>, std::shared_ptr<TH3>, std::shared_ptr<TH2>, std::shared_ptr<TH1>, std::shared_ptr<TProfile3D>, std::shared_ptr<TProfile2D>, std::shared_ptr<TProfile>>;
using HistPtr = std::variant<std::shared_ptr<THn>, std::shared_ptr<THnSparse>, std::shared_ptr<TH3>, std::shared_ptr<TH2>, std::shared_ptr<TH1>, std::shared_ptr<TProfile3D>, std::shared_ptr<TProfile2D>, std::shared_ptr<TProfile>, std::shared_ptr<StepTHn>>;

//**************************************************************************************************
/**
Expand Down Expand Up @@ -117,9 +117,10 @@ struct AxisSpec {
*/
//**************************************************************************************************
struct HistogramConfigSpec {
HistogramConfigSpec(HistType type_, std::vector<AxisSpec> axes_)
HistogramConfigSpec(HistType type_, std::vector<AxisSpec> axes_, uint8_t nSteps_ = 1)
: type(type_),
axes(axes_)
axes(axes_),
nSteps(nSteps_)
{
}
HistogramConfigSpec() = default;
Expand Down Expand Up @@ -154,6 +155,7 @@ struct HistogramConfigSpec {

HistType type{HistType::kUndefinedHist};
std::vector<AxisSpec> axes{};
uint32_t nSteps{1}; // variable used only in StepTHn
};

//**************************************************************************************************
Expand Down Expand Up @@ -218,7 +220,7 @@ struct HistFactory {
}

// create histogram
std::shared_ptr<T> hist{generateHist<T>(histSpec.name, histSpec.title, nAxes, nBins, lowerBounds, upperBounds)};
std::shared_ptr<T> hist{generateHist<T>(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;
Expand Down Expand Up @@ -280,7 +282,7 @@ struct HistFactory {
template <typename T>
static TAxis* getAxis(const int i, std::shared_ptr<T>& hist)
{
if constexpr (std::is_base_of_v<THnBase, T>) {
if constexpr (std::is_base_of_v<THnBase, T> || std::is_base_of_v<StepTHn, T>) {
return hist->GetAxis(i);
} else {
return (i == 0) ? hist->GetXaxis()
Expand All @@ -296,9 +298,11 @@ struct HistFactory {
// helper function to generate the actual histograms
template <typename T>
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<THnBase, T>) {
if constexpr (std::is_base_of_v<StepTHn, T>) {
return new T(name.data(), title.data(), nSteps, nDim, nBins, lowerBounds, upperBounds);
} else if constexpr (std::is_base_of_v<THnBase, T>) {
return new T(name.data(), title.data(), nDim, nBins, lowerBounds, upperBounds);
} else if constexpr (std::is_base_of_v<TH3, T>) {
return (nDim != 3) ? nullptr
Expand Down Expand Up @@ -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<THn, THnSparse, TProfile3D, TH3, TProfile2D, TH2, TProfile, TH1>(obj);
return castToVariant<THn, THnSparse, TProfile3D, TH3, TProfile2D, TH2, TProfile, TH1, StepTHn>(obj);
}
return std::nullopt;
}
Expand Down Expand Up @@ -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<THnBase, T>;
constexpr bool validComplexFillStep = std::is_base_of_v<StepTHn, T>;

if constexpr (validSimpleFill) {
hist->Fill(static_cast<double>(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<double>(positionAndWeight)...};
double weight{1.};
Expand All @@ -391,6 +398,10 @@ struct HistFiller {
template <typename... Cs, typename R, typename T>
static void fillHistAny(std::shared_ptr<R>& hist, const T& table, const o2::framework::expressions::Filter& filter)
{
if constexpr (std::is_base_of_v<StepTHn, T>) {
LOGF(FATAL, "Table filling is not (yet?) supported for StepTHn.");
return;
}
auto filtered = o2::soa::Filtered<T>{{table.asArrowTable()}, o2::framework::expressions::createSelection(table.asArrowTable(), filter)};
for (auto& t : filtered) {
fillHistAny(hist, (*(static_cast<Cs>(t).getIterator()))...);
Expand Down Expand Up @@ -615,7 +626,7 @@ class HistogramRegistry
template <typename T>
void insertClone(const HistName& histName, const std::shared_ptr<T>& 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)]);
Expand All @@ -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 <typename T>
Expand Down
16 changes: 12 additions & 4 deletions Framework/Core/include/Framework/StepTHn.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<Double_t> binLimits[], const char** axisTitles);
StepTHn(const Char_t* name, const Char_t* title, const Int_t nSteps, const Int_t nAxes);
~StepTHn() override;

template <typename... Ts>
Expand All @@ -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;
Expand Down Expand Up @@ -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<Double_t> 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<Double_t> 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:
Expand All @@ -110,6 +114,10 @@ typedef StepTHnT<TArrayD> StepTHnD;
template <typename... Ts>
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<double>(valuesAndWeight)...};
Expand All @@ -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
Expand Down Expand Up @@ -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();
Expand Down
17 changes: 13 additions & 4 deletions Framework/Core/src/HistogramRegistry.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
// or submit itself to any jurisdiction.

#include "Framework/HistogramRegistry.h"
#include <regex>

namespace o2::framework
{
Expand All @@ -29,15 +30,14 @@ const std::map<HistType, std::function<HistPtr(const HistogramSpec&)>> 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;
Expand All @@ -53,15 +53,24 @@ 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();
std::string collidingName{};
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)
Expand Down
79 changes: 43 additions & 36 deletions Framework/Core/src/StepTHn.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<Double_t> 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 <nSteps> 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 <nAxis> dimensional histogram is created.
// For each step a <nAxes> dimensional histogram is created.
// The axis have <nBins[i]> bins. The bin edges are given in <binEdges[i]>. If there are only two bin edges, equidistant binning is set.

init();
}

// root-like constructor
template <class TemplateArray>
StepTHnT<TemplateArray>::StepTHnT(const Char_t* name, const Char_t* title, const Int_t nSteps, const Int_t nAxis,
Int_t* nBins, std::vector<Double_t> binEdges[], const char** axisTitles) : StepTHn(name, title, nSteps, nAxis, nBins, binEdges, axisTitles)
StepTHnT<TemplateArray>::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<TemplateArray, TArrayD>) {
mPrototype = new THnSparseD(Form("%s_sparse", name), title, nAxis, nBins);
} else {
mPrototype = new THnSparseF(Form("%s_sparse", name), title, nAxis, nBins);
mPrototype = new THnSparseT<TemplateArray>(Form("%s_sparse", name), title, nAxes, nBins, xmin, xmax);
}

template <class TemplateArray>
StepTHnT<TemplateArray>::StepTHnT(const Char_t* name, const Char_t* title, const Int_t nSteps, const Int_t nAxes,
Int_t* nBins, std::vector<Double_t> 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<TemplateArray>(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]);
Expand Down Expand Up @@ -202,15 +210,14 @@ void StepTHn::Copy(TObject& c) const
}

if (mPrototype) {
target.mPrototype = dynamic_cast<THnSparseF*>(mPrototype->Clone());
target.mPrototype = dynamic_cast<THnSparse*>(mPrototype->Clone());
}
}

template <class TemplateArray>
Long64_t StepTHnT<TemplateArray>::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) {
Expand Down
Loading