Skip to content

Commit fb9f2d1

Browse files
mario-kruegermario-krueger
authored
add some histogram helpers (#4239)
1 parent e2c74a2 commit fb9f2d1

3 files changed

Lines changed: 426 additions & 0 deletions

File tree

‎Analysis/Core/include/Analysis/HistHelpers.h‎

Lines changed: 274 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,274 @@
1+
// Copyright CERN and copyright holders of ALICE O2. This software is
2+
// distributed under the terms of the GNU General Public License v3 (GPL
3+
// Version 3), copied verbatim in the file "COPYING".
4+
//
5+
// See http://alice-o2.web.cern.ch/license for full licensing information.
6+
//
7+
// In applying this license CERN does not waive the privileges and immunities
8+
// granted to it by virtue of its status as an Intergovernmental Organization
9+
// or submit itself to any jurisdiction.
10+
11+
/// \author Mario Krueger <mario.kruger@cern.ch>
12+
//
13+
// Some helper templates to simplify working with histograms
14+
//
15+
16+
#ifndef HistHelpers_H
17+
#define HistHelpers_H
18+
19+
#include <TH1.h>
20+
#include <TH2.h>
21+
#include <TH3.h>
22+
#include <THn.h>
23+
#include <THnSparse.h>
24+
#include <TObjArray.h>
25+
26+
#include "Framework/Logger.h"
27+
28+
namespace o2
29+
{
30+
namespace experimental
31+
{
32+
namespace histhelpers
33+
{
34+
35+
template <typename T>
36+
struct is_shared_ptr : std::false_type {
37+
};
38+
template <typename T>
39+
struct is_shared_ptr<std::shared_ptr<T>> : std::true_type {
40+
};
41+
42+
//**************************************************************************************************
43+
/**
44+
* Container class for storing and saving root histograms of any type.
45+
* TObjArray inheritance (and concomitant raw pointer gynmnastics) is used to interface with existing O2 file writing functionality.
46+
*/
47+
//**************************************************************************************************
48+
class HistContainer : public TObjArray
49+
{
50+
public:
51+
HistContainer(const std::string& name) : TObjArray()
52+
{
53+
SetBit(TObject::kSingleKey, false); // not working; seems like WriteObjectAny ignores this...
54+
SetOwner(false); // let container handle object deletion
55+
SetName(name.data());
56+
}
57+
58+
using HistType = std::variant<std::shared_ptr<TH3>, std::shared_ptr<TH2>, std::shared_ptr<TH1>, std::shared_ptr<THn>, std::shared_ptr<THnSparse>>;
59+
60+
template <typename T>
61+
void Add(uint8_t histID, T&& hist)
62+
{
63+
if (mHistos.find(histID) != mHistos.end()) {
64+
LOGF(WARNING, "HistContainer %s already holds a histogram at histID = %d. Overriding it now...", GetName(), histID);
65+
TObject* oldPtr = nullptr;
66+
std::visit([&](const auto& sharedPtr) { oldPtr = sharedPtr.get(); }, mHistos[histID]);
67+
TObjArray::Remove(oldPtr);
68+
}
69+
// if shared poiner is provided as argument, object itself is used, otherwise copied
70+
if constexpr (is_shared_ptr<T>::value)
71+
mHistos[histID] = DownCast(hist);
72+
else {
73+
mHistos[histID] = DownCast(std::make_shared<T>(hist));
74+
}
75+
TObject* rawPtr = nullptr;
76+
std::visit([&](const auto& sharedPtr) { rawPtr = sharedPtr.get(); }, mHistos[histID]);
77+
TObjArray::Add(rawPtr);
78+
}
79+
80+
// gets the underlying histogram pointer
81+
// unfortunately we cannot automatically infer type here
82+
// so one has to use Get<TH1>(), Get<TH2>(), Get<TH3>(), Get<THn>(), Get<THnSparse>()
83+
template <typename T>
84+
std::shared_ptr<T> Get(uint8_t histID)
85+
{
86+
return *std::get_if<std::shared_ptr<T>>(&mHistos[histID]);
87+
}
88+
89+
template <typename... Ts>
90+
void Fill(uint8_t histID, Ts&&... position)
91+
{
92+
std::visit([this, &position...](auto&& hist) { GenericFill(hist, std::forward<Ts>(position)...); }, mHistos[histID]);
93+
}
94+
95+
private:
96+
std::map<uint8_t, HistType> mHistos;
97+
98+
// disallow user to call TObjArray::Add()
99+
using TObjArray::Add;
100+
101+
template <typename T, typename... Ts>
102+
void GenericFill(std::shared_ptr<T> hist, const Ts&... position)
103+
{
104+
constexpr bool isValidTH3 = (std::is_same_v<TH3, T> && sizeof...(Ts) == 3);
105+
constexpr bool isValidTH2 = (std::is_same_v<TH2, T> && sizeof...(Ts) == 2);
106+
constexpr bool isValidTH1 = (std::is_same_v<TH1, T> && sizeof...(Ts) == 1);
107+
constexpr bool isValidPrimitive = isValidTH1 || isValidTH2 || isValidTH3;
108+
// unfortunately we dont know at compile the dimension of THn(Sparse)
109+
constexpr bool isValidComplex = std::is_base_of<THnBase, T>::value;
110+
111+
if constexpr (isValidPrimitive) {
112+
hist->Fill(static_cast<double>(position)...);
113+
} else if constexpr (isValidComplex) {
114+
// savety check for n dimensional histograms (runtime overhead)
115+
// if(hist->GetNdimensions() != sizeof...(position)) return;
116+
double tempArray[] = {static_cast<double>(position)...};
117+
hist->Fill(tempArray);
118+
}
119+
};
120+
121+
template <typename T>
122+
HistType DownCast(std::shared_ptr<T> hist)
123+
{
124+
// since the variant HistType only knows the interface classes (TH1,TH2,TH3)
125+
// assigning an actual derived type of TH2 and TH3 (e.g. TH2F, TH3I)
126+
// will confuse the compiler since they are both TH2(3) and TH1 at the same time
127+
// it will not know which alternative of HistType to select
128+
// therefore, in these cases we have to explicitly cast to the correct interface type first
129+
if constexpr (std::is_base_of_v<TH3, T>) {
130+
return std::static_pointer_cast<TH3>(hist);
131+
} else if constexpr (std::is_base_of_v<TH2, T>) {
132+
return std::static_pointer_cast<TH2>(hist);
133+
} else {
134+
// all other cases can be left untouched
135+
return hist;
136+
}
137+
}
138+
};
139+
140+
//**************************************************************************************************
141+
/**
142+
* Helper class to build all kinds of root histograms with a streamlined user interface.
143+
*/
144+
//**************************************************************************************************
145+
146+
struct Axis {
147+
std::string name{};
148+
std::string title{};
149+
std::vector<double> binEdges{};
150+
int nBins{}; // 0 when bin edges are specified directly FIXME: make this std::optional
151+
};
152+
153+
template <typename RootHist_t>
154+
class HistBuilder
155+
{
156+
public:
157+
HistBuilder() : mAxes{}, mHist{nullptr} {}
158+
HistBuilder(const HistBuilder&) = delete; // non construction-copyable
159+
HistBuilder& operator=(const HistBuilder&) = delete; // non copyable
160+
161+
void AddAxis(const Axis& axis) { mAxes.push_back(axis); }
162+
void AddAxes(const std::vector<Axis>& axes)
163+
{
164+
mAxes.insert(mAxes.end(), axes.begin(), axes.end());
165+
}
166+
void AddAxis(const std::string& name, const std::string& title, const int& nBins,
167+
const double& lowerEdge, const double& upperEdge)
168+
{
169+
mAxes.push_back({name, title, {lowerEdge, upperEdge}, nBins});
170+
}
171+
void AddAxis(const std::string& name, const std::string& title,
172+
const std::vector<double>& binEdges)
173+
{
174+
mAxes.push_back({name, title, binEdges, 0});
175+
}
176+
177+
std::shared_ptr<RootHist_t> GenerateHist(const std::string& name, bool hasWeights = false)
178+
{
179+
const std::size_t MAX_DIM{10};
180+
const std::size_t nAxes{mAxes.size()};
181+
if (nAxes == 0 || nAxes > MAX_DIM)
182+
return nullptr;
183+
184+
int nBins[MAX_DIM]{0};
185+
double lowerBounds[MAX_DIM]{0.0};
186+
double upperBounds[MAX_DIM]{0.0};
187+
188+
// first figure out number of bins and dimensions
189+
std::string title = "[ ";
190+
for (std::size_t i = 0; i < nAxes; i++) {
191+
nBins[i] = (mAxes[i].nBins) ? mAxes[i].nBins : mAxes[i].binEdges.size() - 1;
192+
lowerBounds[i] = mAxes[i].binEdges.front();
193+
upperBounds[i] = mAxes[i].binEdges.back();
194+
title += mAxes[i].name;
195+
if (i < nAxes - 1)
196+
title += " : ";
197+
else
198+
title += " ]";
199+
}
200+
201+
// create histogram
202+
mHist.reset(HistFactory(name, title, nAxes, nBins, lowerBounds, upperBounds));
203+
204+
if (!mHist) {
205+
LOGF(WARNING, "ERROR: The number of specified dimensions does not match the type.");
206+
return nullptr;
207+
}
208+
209+
// set axis properties
210+
for (std::size_t i = 0; i < nAxes; i++) {
211+
TAxis* axis{GetAxis(i)};
212+
if (axis) {
213+
axis->SetTitle(mAxes[i].title.data());
214+
if constexpr (std::is_base_of_v<THnBase, RootHist_t>)
215+
axis->SetName((std::to_string(i) + "-" + mAxes[i].name).data());
216+
217+
// move the bin edges in case a variable binnining was requested
218+
if (!mAxes[i].nBins) {
219+
if (!std::is_sorted(std::begin(mAxes[i].binEdges), std::end(mAxes[i].binEdges))) {
220+
LOGF(WARNING, "ERROR: The bin edges specified for axis %s in histogram %s are not in increasing order!", mAxes[i].name, name);
221+
return nullptr;
222+
}
223+
axis->Set(nBins[i], mAxes[i].binEdges.data());
224+
}
225+
}
226+
}
227+
if (hasWeights)
228+
mHist->Sumw2();
229+
mAxes.clear(); // clean up after creating the root histogram
230+
return mHist;
231+
}
232+
233+
private:
234+
std::vector<Axis> mAxes;
235+
std::shared_ptr<RootHist_t> mHist;
236+
237+
RootHist_t* HistFactory(const std::string& name, const std::string& title, const std::size_t nDim,
238+
const int nBins[], const double lowerBounds[], const double upperBounds[])
239+
{
240+
if constexpr (std::is_base_of_v<THnBase, RootHist_t>) {
241+
return new RootHist_t(name.data(), title.data(), nDim, nBins, lowerBounds, upperBounds);
242+
} else if constexpr (std::is_base_of_v<TH3, RootHist_t>) {
243+
return (nDim != 3) ? nullptr
244+
: new RootHist_t(name.data(), title.data(), nBins[0], lowerBounds[0],
245+
upperBounds[0], nBins[1], lowerBounds[1], upperBounds[1],
246+
nBins[2], lowerBounds[2], upperBounds[2]);
247+
} else if constexpr (std::is_base_of_v<TH2, RootHist_t>) {
248+
return (nDim != 2) ? nullptr
249+
: new RootHist_t(name.data(), title.data(), nBins[0], lowerBounds[0],
250+
upperBounds[0], nBins[1], lowerBounds[1], upperBounds[1]);
251+
} else if constexpr (std::is_base_of_v<TH1, RootHist_t>) {
252+
return (nDim != 1)
253+
? nullptr
254+
: new RootHist_t(name.data(), title.data(), nBins[0], lowerBounds[0], upperBounds[0]);
255+
}
256+
return nullptr;
257+
}
258+
259+
TAxis* GetAxis(const int i)
260+
{
261+
if constexpr (std::is_base_of_v<THnBase, RootHist_t>) {
262+
return mHist->GetAxis(i);
263+
} else {
264+
return (i == 0) ? mHist->GetXaxis()
265+
: (i == 1) ? mHist->GetYaxis() : (i == 2) ? mHist->GetZaxis() : nullptr;
266+
}
267+
}
268+
};
269+
270+
} // namespace histhelpers
271+
} // namespace experimental
272+
} // namespace o2
273+
274+
#endif

‎Analysis/Tutorials/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -119,3 +119,8 @@ o2_add_dpl_workflow(extended-columns
119119
o2_add_dpl_workflow(zdc-vzero-iteration
120120
SOURCES src/ZDCVZeroIteration.cxx
121121
COMPONENT_NAME AnalysisTutorial)
122+
123+
o2_add_dpl_workflow(hist-helpers-test
124+
SOURCES src/histHelpersTest.cxx
125+
PUBLIC_LINK_LIBRARIES O2::AnalysisDataModel O2::AnalysisCore
126+
COMPONENT_NAME AnalysisTutorial)

0 commit comments

Comments
 (0)