YODA is hosted by Hepforge, IPPP Durham
YODA - Yet more Objects for Data Analysis  1.3.1
ROOTCnv.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // This file is part of YODA -- Yet more Objects for Data Analysis
4 // Copyright (C) 2008-2015 The YODA collaboration (see AUTHORS for details)
5 //
6 #ifndef YODA_ROOTCnv_h
7 #define YODA_ROOTCnv_h
8 
9 #include "YODA/Histo1D.h"
10 #include "YODA/Histo2D.h"
11 #include "YODA/Profile1D.h"
12 
13 #include "TH1.h"
14 #include "TH2.h"
15 #include "TProfile.h"
16 #include "TGraphAsymmErrors.h"
17 #include "TVectorF.h"
18 
19 namespace YODA {
20 
21 
23 
24 
26 
28 
30 
31 
32  // /// @brief Convert a ROOT 1D histogram to a YODA Histo1D
33  // ///
34  // /// Note that ROOT's histograms do not contain enough information to properly rebuild
35  // /// @a x distributions within bins, in underflow and overflow bins, or across the whole histogram.
36  // inline Histo1D toHisto1D(const TH1& th1) {
37  // std::vector<HistoBin1D> bins;
38  // TArrayD sumw2s = *th1.GetSumw2();
39  // Dbn1D dbn_uflow, dbn_oflow;
40  // double sumWtot(0), sumW2tot(0)
41  // for (int i = 0; i =< th1.GetNbinsX()+1; ++i) {
42  // Dbn1D dbn(static_cast<unsigned long>(th1.GetBinContent(i)), th1.GetBinContent(i), sumw2s[i], 0, 0);
43  // // th1.GetBinContent(i)*th1.GetBinCenter(i), th1.GetBinContent(i)*sqr(th1.GetBinCenter(i)));
44  // if (i == 0) dbn_uflow = dbn;
45  // else if (i == th1.GetNbinsX()+1) dbn_oflow = dbn;
46  // else bins.push_back(HistoBin1D(std::make_pair(th1.GetBinLowEdge(i), th1.GetBinLowEdge(i+1)), dbn));
47  // sumWtot += th1.GetBinContent(i);
48  // sumW2tot += sumw2s[i];
49  // }
50  // Dbn1D dbn_tot(static_cast<unsigned long>(th1.GetEntries()), sumWtot, sumW2tot, 0, 0);
51 
52  // Histo1D rtn(bins, dbn_tot, dbn_uflow, const Dbn1D& dbn_oflow, th1.GetName(), th1.GetTitle());
53  // rtn.addAnnotation("XLabel", th1.GetXaxis->GetTitle());
54  // rtn.addAnnotation("YLabel", th1.GetYaxis->GetTitle());
55  // return rtn;
56  // }
57 
58 
59  // /// @brief Convert a ROOT 1D histogram to a YODA Histo1D
60  // ///
61  // /// Note that ROOT's histograms do not contain enough information to properly rebuild
62  // /// @a x distributions within bins, in underflow and overflow bins, or across the whole histogram.
63  // inline Histo1D toHisto1D(const TH1* th1) {
64  // return toHisto1D(*th1);
65  // }
66 
67 
69 
70 
74  inline Scatter2D toScatter2D(const TH1& th1) {
75  Scatter2D rtn;
76  for (int i = 1; i <= th1.GetNbinsX(); ++i) {
77  const double x = th1.GetBinCenter(i);
78  const double exminus = x - th1.GetBinLowEdge(i);
79  const double explus = th1.GetBinLowEdge(i+1) - x;
80  const double width = exminus + explus;
81  rtn.addPoint(x, th1.GetBinContent(i)/width,
82  exminus, explus,
83  th1.GetBinErrorLow(i)/width, th1.GetBinErrorUp(i)/width);
84  }
85  rtn.addAnnotation("XLabel", th1.GetXaxis()->GetTitle());
86  rtn.addAnnotation("YLabel", th1.GetYaxis()->GetTitle());
87  return rtn;
88  }
89 
90 
92  inline Scatter2D toScatter2D(const TH1* th1) {
93  return toScatter2D(*th1);
94  }
95 
97 
98 
99 
101 
102 
103 
105 
106 
110  inline TH1D toTH1D(const Histo1D& h) {
111  // Work out bin edges first
112  std::vector<double> edges;
113  edges.reserve(h.numBins()+1);
114  edges.push_back(h.bin(0).xMin());
115  for (size_t i = 0; i < h.numBins(); ++i) {
116  const HistoBin1D& b = h.bin(i);
117  if (!fuzzyEquals(edges.back(), b.xMin())) edges.push_back(b.xMin());
118  if (!fuzzyEquals(edges.back(), b.xMax())) edges.push_back(b.xMax());
119  }
120  // Book ROOT histogram
121  TH1D rtn(h.path().c_str(), h.title().c_str(), edges.size()-1, &edges[0]);
122  rtn.Sumw2();
123  TArrayD& sumw2s = *rtn.GetSumw2();
124  for (int i = 1; i <= rtn.GetNbinsX(); ++i) {
125  try {
126  const HistoBin1D& b = h.binAt(rtn.GetBinCenter(i)); // throws if in a gap
127  rtn.SetBinContent(i, b.sumW());
128  sumw2s[i] = b.sumW2();
129  } catch (const Exception& e) { }
130  }
131  // Overflows
132  rtn.SetBinContent(0, h.underflow().sumW());
133  rtn.SetBinContent(rtn.GetNbinsX()+1, h.overflow().sumW());
134  sumw2s[0] = h.underflow().sumW2();
135  sumw2s[rtn.GetNbinsX()+1] = h.overflow().sumW2();
136  // Labels
137  if (h.hasAnnotation("XLabel")) rtn.SetXTitle(h.annotation("XLabel").c_str());
138  if (h.hasAnnotation("YLabel")) rtn.SetYTitle(h.annotation("YLabel").c_str());
139  return rtn;
140  }
141 
142 
143 
147  inline TH2D toTH2D(const Histo2D& h) {
148  // Work out bin edges first
149  std::vector<double> xedges, yedges;
150  xedges.reserve(h.numBins()+1);
151  yedges.reserve(h.numBins()+1);
152  xedges.push_back(h.bin(0).xMin());
153  yedges.push_back(h.bin(0).yMin());
154  for (size_t i = 0; i < h.numBins(); ++i) {
155  const HistoBin2D& b = h.bin(i);
156  xedges.push_back(b.xMin());
157  xedges.push_back(b.xMax());
158  yedges.push_back(b.yMin());
159  yedges.push_back(b.yMax());
160  }
161  // Sort and remove (fuzzy) duplicate edges
162  std::sort(xedges.begin(), xedges.end());
163  std::sort(yedges.begin(), yedges.end());
164  const std::vector<double>::iterator xlast = std::unique(xedges.begin(), xedges.end());
165  const std::vector<double>::iterator ylast = std::unique(yedges.begin(), yedges.end());
166  xedges.erase(xlast, xedges.end());
167  yedges.erase(ylast, yedges.end());
168 
169  // Book ROOT histogram
170  TH2D rtn(h.path().c_str(), h.title().c_str(), xedges.size()-1, &xedges[0], yedges.size()-1, &yedges[0]);
171  rtn.Sumw2();
172  TArrayD& sumw2s = *rtn.GetSumw2();
173  for (int ix = 1; ix <= rtn.GetNbinsX(); ++ix) {
174  for (int iy = 1; iy <= rtn.GetNbinsY(); ++iy) {
175  const int i = rtn.GetBin(ix, iy);
176  try {
177  const HistoBin2D& b = h.binAt(rtn.GetBinCenter(ix), rtn.GetBinCenter(iy)); // throws if in a gap
178  rtn.SetBinContent(i, b.sumW());
179  sumw2s[i] = b.sumW2();
180  } catch (const Exception& e) { }
181  }
182  }
183  // Overflows
185  // rtn.SetBinContent(0, h.underflow().sumW());
186  // rtn.SetBinContent(rtn.GetNbinsX()+1, h.overflow().sumW());
187  // sumw2s[0] = h.underflow().sumW2();
188  // sumw2s[rtn.GetNbinsX()+1] = h.overflow().sumW2();
189  // Labels
190  if (h.hasAnnotation("XLabel")) rtn.SetXTitle(h.annotation("XLabel").c_str());
191  if (h.hasAnnotation("YLabel")) rtn.SetYTitle(h.annotation("YLabel").c_str());
192  if (h.hasAnnotation("ZLabel")) rtn.SetZTitle(h.annotation("ZLabel").c_str());
193  return rtn;
194  }
195 
196 
197 
199 
200 
204  inline TProfile toTProfile(const Profile1D& p) {
205  // Work out bin edges first
206  std::vector<double> edges;
207  edges.reserve(p.numBins()+1);
208  edges.push_back(p.bin(0).xMin());
209  for (size_t i = 0; i < p.numBins(); ++i) {
210  const ProfileBin1D& b = p.bin(i);
211  if (!fuzzyEquals(edges.back(), b.xMin())) edges.push_back(b.xMin());
212  if (!fuzzyEquals(edges.back(), b.xMax())) edges.push_back(b.xMax());
213  }
214  // Book ROOT histogram
215  TProfile rtn(p.path().c_str(), p.title().c_str(), edges.size()-1, &edges[0]);
216  rtn.Sumw2();
217  Double_t* sumwys = rtn.GetArray(); //< YUCK!!!
218  TArrayD& sumwy2s = *rtn.GetSumw2(); //< YUCK!!!
219  for (int i = 1; i <= rtn.GetNbinsX(); ++i) {
220  try {
221  const ProfileBin1D& b = p.binAt(rtn.GetBinCenter(i)); // throws if in a gap
222  // rtn.SetBinContent(i, b.mean());
223  // rtn.SetBinError(i, b.stdErr());
225  // - for sum(y*y): TProfile::GetSumw2()
226  // - for sum(y): TProfile::GetArray()
227  // - for sum(w): TProfile::SetBinEntries(bin, w)
228  // Clearly, the names of accessors/methods are confusing...
229  rtn.SetBinEntries(i, b.sumW());
230  sumwys[i] = b.sumWY();
231  sumwy2s[i] = b.sumWY2();
232  } catch (const Exception& e) { }
233  }
234  // Overflows
235  rtn.SetBinEntries(0, p.underflow().sumW());
236  rtn.SetBinEntries(rtn.GetNbinsX()+1, p.overflow().sumW());
237  sumwys[0] = p.underflow().sumWY();
238  sumwys[0] = p.underflow().sumWY();
239  sumwy2s[rtn.GetNbinsX()+1] = p.overflow().sumWY2();
240  sumwy2s[rtn.GetNbinsX()+1] = p.overflow().sumWY2();
241  // Labels
242  if (p.hasAnnotation("XLabel")) rtn.SetXTitle(p.annotation("XLabel").c_str());
243  if (p.hasAnnotation("YLabel")) rtn.SetYTitle(p.annotation("YLabel").c_str());
244  return rtn;
245  }
246 
247 
249  inline TH1D toTH1D(const Profile1D& p) {
250  // Work out bin edges first
251  std::vector<double> edges;
252  edges.reserve(p.numBins()+1);
253  edges.push_back(p.bin(0).xMin());
254  for (size_t i = 0; i < p.numBins(); ++i) {
255  const ProfileBin1D& b = p.bin(i);
256  if (!fuzzyEquals(edges.back(), b.xMin())) edges.push_back(b.xMin());
257  if (!fuzzyEquals(edges.back(), b.xMax())) edges.push_back(b.xMax());
258  }
259  // Book ROOT histogram
260  TH1D rtn(p.path().c_str(), p.title().c_str(), edges.size()-1, &edges[0]);
261  for (int i = 1; i <= rtn.GetNbinsX(); ++i) {
262  try {
263  const ProfileBin1D& b = p.binAt(rtn.GetBinCenter(i)); // throws if in a gap
264  rtn.SetBinContent(i, b.mean());
265  rtn.SetBinError(i, b.stdErr());
266  } catch (const Exception& e) { }
267  }
268  // Overflows
269  rtn.SetBinContent(0, p.underflow().yMean());
270  rtn.SetBinContent(rtn.GetNbinsX()+1, p.overflow().yMean());
271  rtn.SetBinError(0, p.underflow().yStdErr());
272  rtn.SetBinError(rtn.GetNbinsX()+1, p.overflow().yStdErr());
273  // Labels
274  if (p.hasAnnotation("XLabel")) rtn.SetXTitle(p.annotation("XLabel").c_str());
275  if (p.hasAnnotation("YLabel")) rtn.SetYTitle(p.annotation("YLabel").c_str());
276  return rtn;
277  }
278 
279 
281 
282 
283 
284 
285 
289  inline TGraphAsymmErrors toTGraph(const Scatter2D& s) {
290  TVectorF xs(s.numPoints()), ys(s.numPoints());
291  TVectorF exls(s.numPoints()), exhs(s.numPoints());
292  TVectorF eyls(s.numPoints()), eyhs(s.numPoints());
293  for (size_t i = 0; i < s.numPoints(); ++i) {
294  const Point2D& p = s.point(i);
295  xs[i] = p.x();
296  ys[i] = p.y();
297  exls[i] = p.xErrMinus();
298  exhs[i] = p.xErrPlus();
299  eyls[i] = p.yErrMinus();
300  eyhs[i] = p.yErrPlus();
301  }
302  // Make the ROOT object... mm, the constructors don't take name+title, unlike all this histos!
303  TGraphAsymmErrors rtn(xs, ys, exls, exhs, eyls, eyhs);
304  rtn.SetName(s.path().c_str());
305  rtn.SetTitle(s.title().c_str());
306  // Labels
307  if (s.hasAnnotation("XLabel")) rtn.GetXaxis()->SetTitle(s.annotation("XLabel").c_str());
308  if (s.hasAnnotation("YLabel")) rtn.GetYaxis()->SetTitle(s.annotation("YLabel").c_str());
309  return rtn;
310  }
311 
312 
315  inline TGraphAsymmErrors toTGraph(const Histo1D& h) {
316  return toTGraph(mkScatter(h));
317  }
318 
319 
322  inline TGraphAsymmErrors toTGraph(const Profile1D& p) {
323  return toTGraph(mkScatter(p));
324  }
325 
326 
328 
329 
331 
332 
333 }
334 
335 #endif
double sumW() const
The sum of weights.
Definition: Dbn1D.h:152
TProfile toTProfile(const Profile1D &p)
Convert a YODA Scatter2D to a ROOT TH1D.
Definition: ROOTCnv.h:204
TH1D toTH1D(const Histo1D &h)
Convert a YODA Histo1D to a ROOT 1D histogram.
Definition: ROOTCnv.h:110
double xMin() const
Lower limit of the bin (inclusive).
Definition: Bin1D.h:120
size_t numBins() const
Number of bins on this axis (not counting under/overflow)
Definition: Profile1D.h:196
void addPoint(const Point2D &pt)
Insert a new point.
Definition: Scatter2D.h:204
const HistoBin1D & binAt(double x) const
Access a bin by coordinate (const version)
Definition: Histo1D.h:203
A very generic data type which is just a collection of 2D data points with errors.
Definition: Scatter2D.h:23
Generic unspecialised YODA runtime error.
Definition: Exceptions.h:20
size_t numBins() const
Number of bins.
Definition: Histo2D.h:260
double xMax() const
Upper x limit of the bin (exclusive).
Definition: Bin2D.h:141
double xMax() const
Upper limit of the bin (exclusive).
Definition: Bin1D.h:125
double sumW() const
The sum of weights.
Definition: Bin2D.h:296
double sumW() const
The sum of weights.
Definition: Dbn2D.h:176
Dbn1D & underflow()
Access underflow (non-const version)
Definition: Histo1D.h:213
bool fuzzyEquals(double a, double b, double tolerance=1E-5)
Compare two floating point numbers for equality with a degree of fuzziness.
Definition: MathUtils.h:68
Dbn2D & overflow()
Access overflow (non-const version)
Definition: Profile1D.h:234
double xMin() const
Lower x limit of the bin (inclusive).
Definition: Bin2D.h:136
Dbn1D & overflow()
Access overflow (non-const version)
Definition: Histo1D.h:219
double sumW2() const
The sum of weights squared.
Definition: Dbn1D.h:157
A one-dimensional histogram.
Definition: Histo1D.h:27
double yErrMinus() const
Get negative y-error value.
Definition: Point2D.h:212
const ProfileBin1D & binAt(double x) const
Access a bin by x-coordinate (const version)
Definition: Profile1D.h:221
A Bin1D specialised for handling histogram-type information.
Definition: HistoBin1D.h:21
size_t numPoints() const
Number of points in the scatter.
Definition: Scatter2D.h:167
bool hasAnnotation(const std::string &name) const
Check if an annotation is defined.
double sumW2() const
The sum of weights squared.
Definition: Bin2D.h:301
const std::string & annotation(const std::string &name) const
Get an annotation by name (as a string)
Point2D & point(size_t index)
Get a reference to the point with index index (non-const)
Definition: Scatter2D.h:185
double mean() const
The mean of the y distribution.
Definition: ProfileBin1D.h:111
double yMin() const
Lower y limit of the bin (inclusive).
Definition: Bin2D.h:152
double sumWY2() const
The sum of y^2*weight.
Definition: Dbn2D.h:201
double stdErr() const
The standard error on the mean.
Definition: ProfileBin1D.h:126
double xErrPlus() const
Get positive x-error value.
Definition: Point2D.h:138
void addAnnotation(const std::string &name, const T &value)
Add or set an annotation by name.
size_t numBins() const
Number of bins on this axis (not counting under/overflow)
Definition: Histo1D.h:176
Scatter2D toScatter2D(const TH1 &th1)
Definition: ROOTCnv.h:74
const HistoBin2D & binAt(double x, double y) const
Access a bin by coordinate (const version)
Definition: Histo2D.h:256
double sumWY() const
The sum of y*weight.
Definition: Dbn2D.h:196
TGraphAsymmErrors toTGraph(const Scatter2D &s)
Convert a YODA Scatter2D to a ROOT TGraphAsymmErrors.
Definition: ROOTCnv.h:289
double x() const
Get x value.
Definition: Point2D.h:99
double sumW() const
The sum of weights.
Definition: Bin1D.h:212
Scatter1D mkScatter(const Counter &c)
Make a Scatter1D representation of a Histo1D.
Definition: Scatter1D.cc:8
double y() const
Get y value.
Definition: Point2D.h:105
A Bin1D specialised for handling profile-type information.
Definition: ProfileBin1D.h:23
A Bin2D specialised for handling histogram-type information.
Definition: HistoBin2D.h:21
HistoBin1D & bin(size_t index)
Access a bin by index (non-const version)
Definition: Histo1D.h:192
double sumWY() const
The sum of y*weight.
Definition: ProfileBin1D.h:147
HistoBin2D & bin(size_t index)
Access a bin by index (non-const version)
Definition: Histo2D.h:247
TH2D toTH2D(const Histo2D &h)
Convert a YODA Histo2D to a ROOT 2D histogram.
Definition: ROOTCnv.h:147
double sumWY2() const
The sum of y^2 * weight.
Definition: ProfileBin1D.h:152
double yMean() const
Weighted mean, , of distribution.
Definition: Dbn2D.h:133
const std::string path() const
double yErrPlus() const
Get positive y-error value.
Definition: Point2D.h:217
A two-dimensional histogram.
Definition: Histo2D.h:30
const std::string title() const
double yMax() const
Upper y limit of the bin (exclusive).
Definition: Bin2D.h:157
A one-dimensional profile histogram.
Definition: Profile1D.h:32
double yStdErr() const
Weighted standard error on the mean, , of distribution.
Definition: Dbn2D.h:151
double xErrMinus() const
Get negative x-error value.
Definition: Point2D.h:133
ProfileBin1D & bin(size_t index)
Access a bin by index (non-const version)
Definition: Profile1D.h:213
Dbn2D & underflow()
Access underflow (non-const version)
Definition: Profile1D.h:229
double sumW2() const
The sum of weights squared.
Definition: Bin1D.h:217
A 2D data point to be contained in a Scatter2D.
Definition: Point2D.h:17