YODA is hosted by Hepforge, IPPP Durham
 YODA - Yet more Objects for Data Analysis  1.3.1
Histo2D.cc
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 #include "YODA/Histo2D.h"
7 #include "YODA/Profile2D.h"
8 #include "YODA/Scatter3D.h"
9 #include <cmath>
10
11 using namespace std;
12
13 namespace YODA {
14
15
16  void Histo2D::fill(double x, double y, double weight) {
17  if ( std::isnan(x) ) throw RangeError("X is NaN");
18  if ( std::isnan(y) ) throw RangeError("Y is NaN");
19
20  // Fill the overall distribution
21  _axis.totalDbn().fill(x, y, weight);
22
23  // Fill the bins and overflows
25  if (inRange(x, _axis.xMin(), _axis.xMax()) && inRange(y, _axis.yMin(), _axis.yMax())) {
26  try {
28  binAt(x, y).fill(x, y, weight);
29  } catch (const RangeError& re) { }
30  }
32  // else {
33  // size_t ix(0), iy(0);
34  // if (x < _axis.xMin()) ix = -1; else if (x >= _axis.xMax()) ix = 1;
35  // if (y < _axis.yMin()) iy = -1; else if (y >= _axis.yMax()) iy = 1;
36  // _axis.outflow(ix, iy).fill(x, y, weight);
37  // }
38
39  // Lock the axis now that a fill has happened
40  _axis._setLock(true);
41  }
42
43
44  void Histo2D::fillBin(size_t i, double weight) {
45  pair<double, double> mid = bin(i).xyMid();
46  fill(mid.first, mid.second, weight);
47  }
48
49
50  double Histo2D::sumW(bool includeoverflows) const {
51  if (includeoverflows) return _axis.totalDbn().sumW();
52  double sumw = 0;
53  BOOST_FOREACH (const HistoBin2D& b, bins()) sumw += b.sumW();
54  return sumw;
55  }
56
57
58  double Histo2D::sumW2(bool includeoverflows) const {
59  if (includeoverflows) return _axis.totalDbn().sumW2();
60  double sumw2 = 0;
61  BOOST_FOREACH (const HistoBin2D& b, bins()) sumw2 += b.sumW2();
62  return sumw2;
63  }
64
65
66  double Histo2D::xMean(bool includeoverflows) const {
67  if (includeoverflows) return _axis.totalDbn().xMean();
68  double sumwx = 0;
69  BOOST_FOREACH (const HistoBin2D& b, bins()) sumwx += b.sumWX();
70  return sumwx/sumW();
71  }
72
73
74  double Histo2D::yMean(bool includeoverflows) const {
75  if (includeoverflows) return _axis.totalDbn().yMean();
76  double sumwy = 0;
77  BOOST_FOREACH (const HistoBin2D& b, bins()) sumwy += b.sumWY();
78  return sumwy/sumW();
79  }
80
81
82  double Histo2D::xVariance(bool includeoverflows) const {
83  if (includeoverflows) return _axis.totalDbn().xVariance();
85  double sigma2 = 0;
86  const double xMean = this->xMean();
87  for (size_t i = 0; i < bins().size(); ++i) {
88  const double diff = bin(i).xFocus() - xMean;
89  sigma2 += diff * diff * bin(i).sumW();
90  }
91  return sigma2/sumW();
92  }
93
94
95  double Histo2D::yVariance(bool includeoverflows) const {
96  if (includeoverflows) return _axis.totalDbn().yVariance();
98  double sigma2 = 0;
99  const double yMean = this->yMean();
100  for (size_t i = 0; i < bins().size(); ++i) {
101  const double diff = bin(i).yFocus() - yMean;
102  sigma2 += diff * diff * bin(i).sumW();
103  }
104  return sigma2/sumW();
105  }
106
107
108  double Histo2D::xStdErr(bool includeoverflows) const {
109  if (includeoverflows) return _axis.totalDbn().xStdErr();
110  const double effNumEntries = sumW(false)*sumW(false)/sumW2(false);
111  return std::sqrt(xVariance(false) / effNumEntries);
112  }
113
114
115  double Histo2D::yStdErr(bool includeoverflows) const {
116  if (includeoverflows) return _axis.totalDbn().yStdErr();
117  const double effNumEntries = sumW(false)*sumW(false)/sumW2(false);
118  return std::sqrt(yVariance(false) / effNumEntries);
119  }
120
121
122  // double Profile2D::xRMS(bool includeoverflows) const {
123  // if (includeoverflows) return _axis.totalDbn().xRMS();
124  // /// @todo Finish
125  // }
126
127
128  // double Profile2D::yRMS(bool includeoverflows) const {
129  // if (includeoverflows) return _axis.totalDbn().yRMS();
130  // /// @todo Finish
131  // }
132
133
135
136
138  Histo2D::Histo2D(const Histo2D& h, const std::string& path)
139  : AnalysisObject("Histo2D", (path.size() == 0) ? h.path() : path, h, h.title()),
140  _axis(h._axis)
141  { }
142
143
145  Histo2D::Histo2D(const Scatter3D& s, const std::string& path)
146  : AnalysisObject("Histo2D", (path.size() == 0) ? s.path() : path, s, s.title())
147  {
148  std::vector<HistoBin2D> bins;
149  BOOST_FOREACH (const Scatter3D::Point& p, s.points()) {
150  bins.push_back(HistoBin2D(p.xMin(), p.xMax(), p.yMin(), p.yMax()));
151  }
152  _axis = Histo2DAxis(bins);
153  }
154
155
157  Histo2D::Histo2D(const Profile2D& p, const std::string& path)
158  : AnalysisObject("Histo2D", (path.size() == 0) ? p.path() : path, p, p.title())
159  {
160  std::vector<HistoBin2D> bins;
161  BOOST_FOREACH (const ProfileBin2D& b, p.bins()) {
162  bins.push_back(HistoBin2D(b.xMin(), b.xMax(), b.yMin(), b.yMax()));
163  }
164  _axis = Histo2DAxis(bins);
165  }
166
167
169
170
171  // Histo1D Histo2D::cutterX(double atY, const std::string& path, const std::string& title) {
172  // if (!_axis.isGrid()) throw GridError("Attempt to cut a Histo2D that is not a grid!");
173
174  // if (atY < yMin() || atY > highEdgeY()) throw RangeError("Y is outside the grid");
175  // vector<HistoBin1D> tempBins;
176
177  // for (double i = binByCoord(xMin(), atY).xMin(); i < highEdgeX(); i += binByCoord(i, atY).widthX()) {
178  // const HistoBin2D& b2 = binByCoord(i, atY);
179  // const Dbn1D dbn2(b2.numEntries(), b2.sumW(), b2.sumW2(), b2.sumWX(), b2.sumWX2());
180  // tempBins.push_back(HistoBin1D(b2.xMin(), b2.highEdgeX(), dbn2));
181  // }
182
183  // /// Setting under/over flows
184  // Dbn2D underflow;
185  // underflow += _axis.outflows()[7][_axis.getBinRow(_axis.getBinIndex(xMin(), atY))];
186
187  // Dbn2D overflow;
188  // overflow += _axis.outflows()[3][_axis.getBinRow(_axis.getBinIndex(xMin(), atY))];
189
190  // return Histo1D(tempBins, _axis.totalDbn().transformX(), underflow.transformX(), overflow.transformX(), path, title);
191
192  // }
193
194
195  // Histo1D Histo2D::cutterY(double atX, const std::string& path, const std::string& title) {
196  // if (!_axis.isGrid()) throw GridError("Attempt to cut a Histo2D that is not a grid!");
197
198  // if (atX < xMin() || atX > highEdgeX()) throw RangeError("X is outside the grid");
199  // vector<HistoBin1D> tempBins;
200
201  // for (double i = binByCoord(atX, yMin()).yMin(); i < highEdgeY(); i += binByCoord(atX, i).widthY()) {
202  // const HistoBin2D& b2 = binByCoord(atX, i);
203  // const Dbn1D dbn2(b2.numEntries(), b2.sumW(), b2.sumW2(), b2.sumWX(), b2.sumWX2());
204  // tempBins.push_back(HistoBin1D(b2.yMin(), b2.highEdgeY(), dbn2));
205  // }
206
207  // // Setting under/over flows
208  // Dbn2D underflow;
209  // underflow += _axis.outflows()[1][_axis.getBinColumn(_axis.getBinIndex(atX, yMin()))];
210
211  // Dbn2D overflow;
212  // overflow += _axis.outflows()[5][_axis.getBinColumn(_axis.getBinIndex(atX, yMin()))];
213  // Dbn2D total = _axis.totalDbn();
214
215  // // Making sure that we rotate our distributions, as we are cutting parallel to Y axis now
216  // total.flipXY();
217  // underflow.flipXY();
218  // overflow.flipXY();
219
220  // return Histo1D(tempBins, total.transformX(), underflow.transformX(), overflow.transformX(), path, title);
221  // }
222
223
224  // Profile1D Histo2D::mkProfileX() {
225  // if (!_axis.isGrid()) throw GridError("Profile1D cannot be made from a histogram that is not a grid!");
226
227  // vector<ProfileBin1D> prof;
228  // for(int i = xMin() + _axis.bin(0).xMid(); i < highEdgeX(); i+= _axis.bin(0).widthX()) {
229  // HistoBin2D& bin(_axis.binByCoord(i, yMin()));
230  // HistoBin2D composite(bin.xMin(), bin.xMax(), bin.yMin(), bin.yMax()) ;
231  // for(int j = yMin() + _axis.bin(0).yMid(); j < highEdgeY(); j += _axis.bin(0).widthY()) {
232  // composite += _axis.binByCoord(i, j);
233  // }
234  // prof.push_back(composite.transformX());
235  // }
236
237  // vector<vector<Dbn2D> >& outflows = _axis.outflows();
238
239  // /// Properly setting an underflow
240  // Dbn2D underflow;
241  // underflow += outflows[0][0]; underflow += outflows[6][0];
242  // for(size_t i = 0; i < outflows[7].size(); ++i) {
243  // underflow += outflows[7][i];
244  // }
245
246  // /// Setting an overflow
247  // Dbn2D overflow;
248  // overflow += outflows[2][0]; overflow += outflows[4][0];
249  // for(size_t i = 0; i < outflows[3].size(); ++i) {
250  // overflow += outflows[3][i];
251  // }
252
253  // /// And constructing a profile 1D from all this data
254  // Profile1D ret(prof, _axis.totalDbn(), underflow, overflow);
255  // return ret;
256
257  // }
258
259  // Profile1D Histo2D::mkProfileY() {
260  // if (!_axis.isGrid()) throw GridError("Profile1D cannot be made from a histogram that is not a grid!");
261
262  // vector<ProfileBin1D> prof;
263  // for(int i = yMin() + _axis.bin(0).yMid(); i < highEdgeY(); i+= _axis.bin(0).widthY()) {
264  // HistoBin2D& bin(_axis.binByCoord(i, yMin()));
265  // HistoBin2D composite(bin.xMin(), bin.xMax(), bin.yMin(), bin.yMax()) ;
266  // for(int j = xMin() + _axis.bin(0).xMid(); j < highEdgeX(); j += _axis.bin(0).widthX()) {
267  // composite += _axis.binByCoord(i, j);
268  // }
269  // prof.push_back(composite.transformY());
270  // }
271
272  // vector<vector<Dbn2D> >& outflows = _axis.outflows();
273
274  // /// Properly setting an underflow
275  // Dbn2D underflow;
276  // underflow += outflows[0][0]; underflow += outflows[2][0];
277  // for(size_t i = 0; i < outflows[1].size(); ++i) {
278  // underflow += outflows[1][i];
279  // }
280
281  // /// Setting an overflow
282  // Dbn2D overflow;
283  // overflow += outflows[6][0]; overflow += outflows[4][0];
284  // for(size_t i = 0; i < outflows[5].size(); ++i) {
285  // overflow += outflows[5][i];
286  // }
287
288  // /// Setting a flipped total distribution
289  // Dbn2D td = _axis.totalDbn();
290  // td.flipXY();
291
292  // /// And constructing a profile 1D from all this data
293  // Profile1D ret(prof, td, underflow, overflow);
294  // return ret;
295  // }
296
297
298  Scatter3D divide(const Histo2D& numer, const Histo2D& denom) {
300  return divide(mkScatter(numer), mkScatter(denom));
301  }
302
303
304  Scatter3D efficiency(const Histo2D& accepted, const Histo2D& total) {
305  Scatter3D tmp = divide(accepted, total);
306  for (size_t i = 0; i < accepted.numBins(); ++i) {
307  const HistoBin2D& b_acc = accepted.bin(i);
308  const HistoBin2D& b_tot = total.bin(i);
309  Point3D& point = tmp.point(i);
310
312
313  // Check that the numerator is consistent with being a subset of the denominator (NOT effNumEntries here!)
314  if (b_acc.numEntries() > b_tot.numEntries() || b_acc.sumW() > b_tot.sumW())
315  throw UserError("Attempt to calculate an efficiency when the numerator is not a subset of the denominator");
316
317  // If no entries on the denominator, set eff = err = 0 and move to the next bin
320  double eff = 0, err = 0;
321  if (b_tot.sumW() != 0) {
322  // Calculate the values and errors
323  // const double eff = b_acc.effNumEntries() / b_tot.effNumEntries();
324  // const double ey = sqrt( b_acc.effNumEntries() * (1 - b_acc.effNumEntries()/b_tot.effNumEntries()) ) / b_tot.effNumEntries();
325  eff = b_acc.sumW() / b_tot.sumW(); //< Actually this is already calculated by the division...
326  err = sqrt(abs( ((1-2*eff)*b_acc.sumW2() + sqr(eff)*b_tot.sumW2()) / sqr(b_tot.sumW()) ));
327  // assert(point.y() == eff); //< @todo Correct? So we don't need to reset the eff on the next line?
328  }
329
331
332  point.setZ(eff, err);
333  }
334  return tmp;
335
336  }
337
338
340
341
342 }
bool inRange(NUM value, NUM low, NUM high, RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN)
Determine if value is in the range low to high, for floating point numbers.
Definition: MathUtils.h:151
A Bin1D specialised for handling profile-type information.
Definition: ProfileBin2D.h:23
std::vector< YODA::HistoBin2D > & bins()
Access the bin vector (non-const version)
Definition: Histo2D.h:241
unsigned long numEntries() const
The number of entries.
Definition: Bin2D.h:286
Scatter1D divide(const Counter &numer, const Counter &denom)
Definition: Counter.cc:24
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 sumW() const
The sum of weights.
Definition: Bin2D.h:296
Error for e.g. use of invalid bin ranges.
Definition: Exceptions.h:34
double xMin() const
Lower x limit of the bin (inclusive).
Definition: Bin2D.h:136
Error for problems introduced outside YODA, to put it nicely.
Definition: Exceptions.h:93
std::vector< YODA::ProfileBin2D > & bins()
Access the bin vector (non-const)
Definition: Profile2D.h:236
NUM sqr(NUM a)
Named number-type squaring operation.
Definition: MathUtils.h:203
Point3D & point(size_t index)
Get a reference to the point with index index.
Definition: Scatter3D.h:199
A two-dimensional profile histogram.
Definition: Profile2D.h:30
double sumW2() const
The sum of weights squared.
Definition: Bin2D.h:301
A 3D data point to be contained in a Scatter3D.
Definition: Point3D.h:17
double yMin() const
Lower y limit of the bin (inclusive).
Definition: Bin2D.h:152
Axis2D< HistoBin2D, Dbn2D > Histo2DAxis
Convenience typedef.
Definition: Histo2D.h:23
Scatter1D mkScatter(const Counter &c)
Make a Scatter1D representation of a Histo1D.
Definition: Scatter1D.cc:8
double yMax() const
Get value plus positive y-error.
Definition: Point3D.h:243
A very generic data type which is just a collection of 3D data points with errors.
Definition: Scatter3D.h:23
A Bin2D specialised for handling histogram-type information.
Definition: HistoBin2D.h:21
AnalysisObject is the base class for histograms and scatters.
HistoBin2D & bin(size_t index)
Access a bin by index (non-const version)
Definition: Histo2D.h:247
double sumWY() const
The sum of y*weight.
Definition: Bin2D.h:311
Points & points()
Get the collection of points (non-const)
Definition: Scatter3D.h:187
double yMin() const
Get value minus negative y-error.
Definition: Point3D.h:238
A two-dimensional histogram.
Definition: Histo2D.h:30
double yMax() const
Upper y limit of the bin (exclusive).
Definition: Bin2D.h:157
Histo2D(const std::string &path="", const std::string &title="")
Default constructor.
Definition: Histo2D.h:44
void setZ(double z)
Set z value.
Definition: Point3D.h:108
double xMin() const
Get value minus negative x-error.
Definition: Point3D.h:175
double sumWX() const
The sum of x*weight.
Definition: Bin2D.h:306
Scatter1D efficiency(const Counter &accepted, const Counter &total)
Calculate an efficiency ratio of two counters.
Definition: Counter.cc:38
double xMax() const
Get value plus positive x-error.
Definition: Point3D.h:180