yoda is hosted by Hepforge, IPPP Durham
YODA - Yet more Objects for Data Analysis  1.9.0
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-2021 The YODA collaboration (see AUTHORS for details)
5 //
6 #include "YODA/Histo2D.h"
7 #include "YODA/Profile2D.h"
8 #include "YODA/Scatter3D.h"
10 
11 using namespace std;
12 
13 namespace YODA {
14 
15 
17  Histo2D::Histo2D(const Histo2D& h, const std::string& path)
18  : AnalysisObject("Histo2D", (path.size() == 0) ? h.path() : path, h, h.title()),
19  _axis(h._axis)
20  { }
21 
22 
24  Histo2D::Histo2D(const Scatter3D& s, const std::string& path)
25  : AnalysisObject("Histo2D", (path.size() == 0) ? s.path() : path, s, s.title())
26  {
27  std::vector<HistoBin2D> bins;
28  for (const Scatter3D::Point& p : s.points()) {
29  bins.push_back(HistoBin2D(p.xMin(), p.xMax(), p.yMin(), p.yMax()));
30  }
31  _axis = Histo2DAxis(bins);
32  }
33 
34 
36  Histo2D::Histo2D(const Profile2D& p, const std::string& path)
37  : AnalysisObject("Histo2D", (path.size() == 0) ? p.path() : path, p, p.title())
38  {
39  std::vector<HistoBin2D> bins;
40  for (const ProfileBin2D& b : p.bins()) {
41  bins.push_back(HistoBin2D(b.xMin(), b.xMax(), b.yMin(), b.yMax()));
42  }
43  _axis = Histo2DAxis(bins);
44  }
45 
46 
48 
49 
50  void Histo2D::fill(double x, double y, double weight, double fraction) {
51  if ( std::isnan(x) ) throw RangeError("X is NaN");
52  if ( std::isnan(y) ) throw RangeError("Y is NaN");
53 
54  // Fill the overall distribution
55  _axis.totalDbn().fill(x, y, weight, fraction);
56 
57  // Fill the bins and overflows
59  if (inRange(x, _axis.xMin(), _axis.xMax()) && inRange(y, _axis.yMin(), _axis.yMax())) {
60  try {
62  _binAt(x, y).fill(x, y, weight, fraction);
63  } catch (const RangeError& re) { }
64  }
66  // else {
67  // size_t ix(0), iy(0);
68  // if (x < _axis.xMin()) ix = -1; else if (x >= _axis.xMax()) ix = 1;
69  // if (y < _axis.yMin()) iy = -1; else if (y >= _axis.yMax()) iy = 1;
70  // _axis.outflow(ix, iy).fill(x, y, weight, fraction);
71  // }
72 
73  // Lock the axis now that a fill has happened
74  _axis._setLock(true);
75  }
76 
77 
78  void Histo2D::fillBin(size_t i, double weight, double fraction) {
79  pair<double, double> mid = bin(i).xyMid();
80  fill(mid.first, mid.second, weight, fraction);
81  }
82 
83 
84 
86 
87  double Histo2D::numEntries(bool includeoverflows) const {
88  if (includeoverflows) return totalDbn().numEntries();
89  unsigned long n = 0;
90  for (const Bin& b : bins()) n += b.numEntries();
91  return n;
92  }
93 
94 
95  double Histo2D::effNumEntries(bool includeoverflows) const {
96  if (includeoverflows) return totalDbn().effNumEntries();
97  double n = 0;
98  for (const Bin& b : bins()) n += b.effNumEntries();
99  return n;
100  }
101 
102 
103  double Histo2D::sumW(bool includeoverflows) const {
104  if (includeoverflows) return _axis.totalDbn().sumW();
105  double sumw = 0;
106  for (const Bin& b : bins()) sumw += b.sumW();
107  return sumw;
108  }
109 
110 
111  double Histo2D::sumW2(bool includeoverflows) const {
112  if (includeoverflows) return _axis.totalDbn().sumW2();
113  double sumw2 = 0;
114  for (const Bin& b : bins()) sumw2 += b.sumW2();
115  return sumw2;
116  }
117 
118 
120 
121 
122  double Histo2D::xMean(bool includeoverflows) const {
123  if (includeoverflows) return _axis.totalDbn().xMean();
124  Dbn2D dbn;
125  for (const HistoBin2D& b : bins()) dbn += b.dbn();
126  return dbn.xMean();
127  }
128 
129 
130  double Histo2D::yMean(bool includeoverflows) const {
131  if (includeoverflows) return _axis.totalDbn().yMean();
132  Dbn2D dbn;
133  for (const HistoBin2D& b : bins()) dbn += b.dbn();
134  return dbn.yMean();
135  }
136 
137 
138  double Histo2D::xVariance(bool includeoverflows) const {
139  if (includeoverflows) return _axis.totalDbn().xVariance();
140  Dbn2D dbn;
141  for (const HistoBin2D& b : bins()) dbn += b.dbn();
142  return dbn.xVariance();
143  }
144 
145 
146  double Histo2D::yVariance(bool includeoverflows) const {
147  if (includeoverflows) return _axis.totalDbn().yVariance();
148  Dbn2D dbn;
149  for (const HistoBin2D& b : bins()) dbn += b.dbn();
150  return dbn.yVariance();
151  }
152 
153 
154  double Histo2D::xStdErr(bool includeoverflows) const {
155  if (includeoverflows) return _axis.totalDbn().xStdErr();
156  Dbn2D dbn;
157  for (const HistoBin2D& b : bins()) dbn += b.dbn();
158  return dbn.xStdErr();
159  }
160 
161 
162  double Histo2D::yStdErr(bool includeoverflows) const {
163  if (includeoverflows) return _axis.totalDbn().yStdErr();
164  Dbn2D dbn;
165  for (const HistoBin2D& b : bins()) dbn += b.dbn();
166  return dbn.yStdErr();
167  }
168 
169 
170  double Histo2D::xRMS(bool includeoverflows) const {
171  if (includeoverflows) return _axis.totalDbn().xRMS();
172  Dbn2D dbn;
173  for (const HistoBin2D& b : bins()) dbn += b.dbn();
174  return dbn.xRMS();
175  }
176 
177 
178  double Histo2D::yRMS(bool includeoverflows) const {
179  if (includeoverflows) return _axis.totalDbn().yRMS();
180  Dbn2D dbn;
181  for (const HistoBin2D& b : bins()) dbn += b.dbn();
182  return dbn.yRMS();
183  }
184 
185 
187 
188 
189  // Histo1D Histo2D::cutterX(double atY, const std::string& path, const std::string& title) {
190  // if (!_axis.isGrid()) throw GridError("Attempt to cut a Histo2D that is not a grid!");
191 
192  // if (atY < yMin() || atY > highEdgeY()) throw RangeError("Y is outside the grid");
193  // vector<HistoBin1D> tempBins;
194 
195  // for (double i = binByCoord(xMin(), atY).xMin(); i < highEdgeX(); i += binByCoord(i, atY).widthX()) {
196  // const HistoBin2D& b2 = binByCoord(i, atY);
197  // const Dbn1D dbn2(b2.numEntries(), b2.sumW(), b2.sumW2(), b2.sumWX(), b2.sumWX2());
198  // tempBins.push_back(HistoBin1D(b2.xMin(), b2.highEdgeX(), dbn2));
199  // }
200 
201  // /// Setting under/over flows
202  // Dbn2D underflow;
203  // underflow += _axis.outflows()[7][_axis.getBinRow(_axis.getBinIndex(xMin(), atY))];
204 
205  // Dbn2D overflow;
206  // overflow += _axis.outflows()[3][_axis.getBinRow(_axis.getBinIndex(xMin(), atY))];
207 
208  // return Histo1D(tempBins, _axis.totalDbn().transformX(), underflow.transformX(), overflow.transformX(), path, title);
209 
210  // }
211 
212 
213  // Histo1D Histo2D::cutterY(double atX, const std::string& path, const std::string& title) {
214  // if (!_axis.isGrid()) throw GridError("Attempt to cut a Histo2D that is not a grid!");
215 
216  // if (atX < xMin() || atX > highEdgeX()) throw RangeError("X is outside the grid");
217  // vector<HistoBin1D> tempBins;
218 
219  // for (double i = binByCoord(atX, yMin()).yMin(); i < highEdgeY(); i += binByCoord(atX, i).widthY()) {
220  // const HistoBin2D& b2 = binByCoord(atX, i);
221  // const Dbn1D dbn2(b2.numEntries(), b2.sumW(), b2.sumW2(), b2.sumWX(), b2.sumWX2());
222  // tempBins.push_back(HistoBin1D(b2.yMin(), b2.highEdgeY(), dbn2));
223  // }
224 
225  // // Setting under/over flows
226  // Dbn2D underflow;
227  // underflow += _axis.outflows()[1][_axis.getBinColumn(_axis.getBinIndex(atX, yMin()))];
228 
229  // Dbn2D overflow;
230  // overflow += _axis.outflows()[5][_axis.getBinColumn(_axis.getBinIndex(atX, yMin()))];
231  // Dbn2D total = _axis.totalDbn();
232 
233  // // Making sure that we rotate our distributions, as we are cutting parallel to Y axis now
234  // total.flipXY();
235  // underflow.flipXY();
236  // overflow.flipXY();
237 
238  // return Histo1D(tempBins, total.transformX(), underflow.transformX(), overflow.transformX(), path, title);
239  // }
240 
241 
242  // Profile1D Histo2D::mkProfileX() {
243  // if (!_axis.isGrid()) throw GridError("Profile1D cannot be made from a histogram that is not a grid!");
244 
245  // vector<ProfileBin1D> prof;
246  // for(int i = xMin() + _axis.bin(0).xMid(); i < highEdgeX(); i+= _axis.bin(0).widthX()) {
247  // HistoBin2D& bin(_axis.binByCoord(i, yMin()));
248  // HistoBin2D composite(bin.xMin(), bin.xMax(), bin.yMin(), bin.yMax()) ;
249  // for(int j = yMin() + _axis.bin(0).yMid(); j < highEdgeY(); j += _axis.bin(0).widthY()) {
250  // composite += _axis.binByCoord(i, j);
251  // }
252  // prof.push_back(composite.transformX());
253  // }
254 
255  // vector<vector<Dbn2D> >& outflows = _axis.outflows();
256 
257  // /// Properly setting an underflow
258  // Dbn2D underflow;
259  // underflow += outflows[0][0]; underflow += outflows[6][0];
260  // for(size_t i = 0; i < outflows[7].size(); ++i) {
261  // underflow += outflows[7][i];
262  // }
263 
264  // /// Setting an overflow
265  // Dbn2D overflow;
266  // overflow += outflows[2][0]; overflow += outflows[4][0];
267  // for(size_t i = 0; i < outflows[3].size(); ++i) {
268  // overflow += outflows[3][i];
269  // }
270 
271  // /// And constructing a profile 1D from all this data
272  // Profile1D ret(prof, _axis.totalDbn(), underflow, overflow);
273  // return ret;
274 
275  // }
276 
277  // Profile1D Histo2D::mkProfileY() {
278  // if (!_axis.isGrid()) throw GridError("Profile1D cannot be made from a histogram that is not a grid!");
279 
280  // vector<ProfileBin1D> prof;
281  // for(int i = yMin() + _axis.bin(0).yMid(); i < highEdgeY(); i+= _axis.bin(0).widthY()) {
282  // HistoBin2D& bin(_axis.binByCoord(i, yMin()));
283  // HistoBin2D composite(bin.xMin(), bin.xMax(), bin.yMin(), bin.yMax()) ;
284  // for(int j = xMin() + _axis.bin(0).xMid(); j < highEdgeX(); j += _axis.bin(0).widthX()) {
285  // composite += _axis.binByCoord(i, j);
286  // }
287  // prof.push_back(composite.transformY());
288  // }
289 
290  // vector<vector<Dbn2D> >& outflows = _axis.outflows();
291 
292  // /// Properly setting an underflow
293  // Dbn2D underflow;
294  // underflow += outflows[0][0]; underflow += outflows[2][0];
295  // for(size_t i = 0; i < outflows[1].size(); ++i) {
296  // underflow += outflows[1][i];
297  // }
298 
299  // /// Setting an overflow
300  // Dbn2D overflow;
301  // overflow += outflows[6][0]; overflow += outflows[4][0];
302  // for(size_t i = 0; i < outflows[5].size(); ++i) {
303  // overflow += outflows[5][i];
304  // }
305 
306  // /// Setting a flipped total distribution
307  // Dbn2D td = _axis.totalDbn();
308  // td.flipXY();
309 
310  // /// And constructing a profile 1D from all this data
311  // Profile1D ret(prof, td, underflow, overflow);
312  // return ret;
313  // }
314 
315 
316  Scatter3D divide(const Histo2D& numer, const Histo2D& denom) {
317  Scatter3D rtn;
318 
319  for (size_t i = 0; i < numer.numBins(); ++i) {
320  const HistoBin2D& b1 = numer.bin(i);
321  const HistoBin2D& b2 = denom.bin(i);
322 
324  if (!fuzzyEquals(b1.xMin(), b2.xMin()) || !fuzzyEquals(b1.xMax(), b2.xMax()))
325  throw BinningError("x binnings are not equivalent in " + numer.path() + " / " + denom.path());
326  if (!fuzzyEquals(b1.yMin(), b2.yMin()) || !fuzzyEquals(b1.yMax(), b2.yMax()))
327  throw BinningError("y binnings are not equivalent in " + numer.path() + " / " + denom.path());
328 
329  // Assemble the x value and error
330  // Use the midpoint of the "bin" for the new central x value, in the absence of better information
331  const double x = b1.xMid();
332  const double exminus = x - b1.xMin();
333  const double explus = b1.xMax() - x;
334 
335  // Assemble the y value and error
336  // Use the midpoint of the "bin" for the new central y value, in the absence of better information
337  const double y = b1.yMid();
338  const double eyminus = y - b1.yMin();
339  const double eyplus = b1.yMax() - y;
340 
341  // Assemble the z value and error
342  double z = std::numeric_limits<double>::quiet_NaN();
343  double ez = std::numeric_limits<double>::quiet_NaN();
344  if (b2.height() == 0 || (b1.height() == 0 && b1.heightErr() != 0)) {
345  // z = std::numeric_limits<double>::quiet_NaN();
346  // ez = std::numeric_limits<double>::quiet_NaN();
347  // throw LowStatsError("Requested division of empty bin");
348  } else {
349  z = b1.height() / b2.height();
351  const double relerr_1 = b1.heightErr() != 0 ? b1.relErr() : 0;
352  const double relerr_2 = b2.heightErr() != 0 ? b2.relErr() : 0;
353  ez = z * sqrt(sqr(relerr_1) + sqr(relerr_2));
354  }
357  //const double eyplus = y * sqrt( sqr(p1.yErrPlus()/p1.y()) + sqr(p2.yErrMinus()/p2.y()) );
358  //const double eyminus = y * sqrt( sqr(p1.yErrMinus()/p1.y()) + sqr(p2.yErrPlus()/p2.y()) );
359  rtn.addPoint(x, y, z, exminus, explus, eyminus, eyplus, ez, ez);
360  }
361 
362  assert(rtn.numPoints() == numer.numBins());
363  return rtn;
364  }
365 
366 
367  Scatter3D efficiency(const Histo2D& accepted, const Histo2D& total) {
368  Scatter3D tmp = divide(accepted, total);
369  for (size_t i = 0; i < accepted.numBins(); ++i) {
370  const HistoBin2D& b_acc = accepted.bin(i);
371  const HistoBin2D& b_tot = total.bin(i);
372  Point3D& point = tmp.point(i);
373 
375 
376  // Check that the numerator is consistent with being a subset of the denominator
378  if (b_acc.numEntries() > b_tot.numEntries())
379  throw UserError("Attempt to calculate an efficiency when the numerator is not a subset of the denominator: "
380  + Utils::toStr(b_acc.numEntries()) + " entries / " + Utils::toStr(b_tot.numEntries()) + " entries");
381 
382  // If no entries on the denominator, set eff = err = 0 and move to the next bin
383  double eff = std::numeric_limits<double>::quiet_NaN();
384  double err = std::numeric_limits<double>::quiet_NaN();
385  try {
386  if (b_tot.sumW() != 0) {
387  eff = b_acc.sumW() / b_tot.sumW(); //< Actually this is already calculated by the division...
388  err = sqrt(abs( ((1-2*eff)*b_acc.sumW2() + sqr(eff)*b_tot.sumW2()) / sqr(b_tot.sumW()) ));
389  }
390  } catch (const LowStatsError& e) {
391  //
392  }
393 
395 
396  point.setZ(eff, err);
397  }
398  return tmp;
399 
400  }
401 
402 
404 
405 
406 }
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:156
double sumW2() const
The sum of weights squared.
Definition: Bin2D.h:321
A Bin1D specialised for handling profile-type information.
Definition: ProfileBin2D.h:23
double xVariance() const
Weighted variance, , of distribution.
Definition: Dbn2D.h:136
double effNumEntries(bool includeoverflows=true) const
Get the effective number of fills.
Definition: Histo2D.cc:95
double yMax() const
Upper y limit of the bin (exclusive).
Definition: Bin2D.h:171
std::vector< YODA::HistoBin2D > & bins()
Access the bin vector (non-const version)
Definition: Histo2D.h:288
double xMin() const
Lower x limit of the bin (inclusive).
Definition: Bin2D.h:150
Scatter1D divide(const Counter &numer, const Counter &denom)
Definition: Counter.cc:25
const std::string path() const
Get the AO path.
double xMean() const
Weighted mean, , of distribution.
Definition: Dbn2D.h:130
double yMin() const
Lower y limit of the bin (inclusive).
Definition: Bin2D.h:166
double yVariance() const
Weighted variance, , of distribution.
Definition: Dbn2D.h:139
double xVariance(bool includeoverflows=true) const
Get the variance in x.
Definition: Histo2D.cc:138
STL namespace.
Error for e.g. use of invalid bin ranges.
Definition: Exceptions.h:34
double xStdErr() const
Weighted standard error on the mean, , of distribution.
Definition: Dbn2D.h:148
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:73
double xMax() const
Upper x limit of the bin (exclusive).
Definition: Bin2D.h:155
std::pair< double, double > xyMid() const
The geometric centre of the bin.
Definition: Bin2D.h:187
virtual void fillBin(size_t i, double weight=1.0, double fraction=1.0)
Fill histo x-y bin i with the given weight.
Definition: Histo2D.cc:78
double xMid() const
Middle of the bin in x.
Definition: Bin2D.h:177
size_t numPoints() const
Number of points in the scatter.
Definition: Scatter3D.h:203
double yStdErr(bool includeoverflows=true) const
Get the standard error in y.
Definition: Histo2D.cc:162
double yMean() const
Weighted mean, , of distribution.
Definition: Dbn2D.h:133
double effNumEntries() const
Effective number of entries .
Definition: Dbn2D.h:171
A 2D distribution.
Definition: Dbn2D.h:16
Error for problems introduced outside YODA, to put it nicely.
Definition: Exceptions.h:100
std::vector< YODA::ProfileBin2D > & bins()
Access the bin vector (non-const)
Definition: Profile2D.h:275
NUM sqr(NUM a)
Named number-type squaring operation.
Definition: MathUtils.h:208
Point3D & point(size_t index)
Get a reference to the point with index index.
Definition: Scatter3D.h:221
A two-dimensional profile histogram.
Definition: Profile2D.h:32
A 3D data point to be contained in a Scatter3D.
Definition: Point3D.h:18
Errors relating to insufficient (effective) statistics.
Definition: Exceptions.h:72
double yVariance(bool includeoverflows=true) const
Get the variance in y.
Definition: Histo2D.cc:146
Error for general binning problems.
Definition: Exceptions.h:27
double yStdErr() const
Weighted standard error on the mean, , of distribution.
Definition: Dbn2D.h:151
double numEntries() const
The number of entries.
Definition: Bin2D.h:306
Axis2D< HistoBin2D, Dbn2D > Histo2DAxis
Convenience typedef.
Definition: Histo2D.h:25
virtual void fill(double x, double y, double weight=1.0, double fraction=1.0)
Fill histo with weight at (x,y)
Definition: Histo2D.cc:50
double numEntries(bool includeoverflows=true) const
Get the number of fills (fractional fills are possible)
Definition: Histo2D.cc:87
double yRMS(bool includeoverflows=true) const
Get the RMS in y.
Definition: Histo2D.cc:178
double xRMS() const
Weighted RMS, , of distribution.
Definition: Dbn2D.h:154
double height() const
The height of a bin.
Definition: HistoBin2D.h:106
double yMean(bool includeoverflows=true) const
Get the mean y.
Definition: Histo2D.cc:130
double xStdErr(bool includeoverflows=true) const
Get the standard error in x.
Definition: Histo2D.cc:154
double xMean(bool includeoverflows=true) const
Get the mean x.
Definition: Histo2D.cc:122
const std::string title() const
Get the AO title.
Dbn2D & totalDbn()
Access summary distribution, including gaps and overflows (non-const version)
Definition: Histo2D.h:321
A very generic data type which is just a collection of 3D data points with errors.
Definition: Scatter3D.h:25
A Bin2D specialised for handling histogram-type information.
Definition: HistoBin2D.h:21
double sumW2(bool includeoverflows=true) const
Get the sum of squared weights in histo.
Definition: Histo2D.cc:111
double yRMS() const
Weighted RMS, , of distribution.
Definition: Dbn2D.h:157
AnalysisObject is the base class for histograms and scatters.
double sumW() const
The sum of weights.
Definition: Bin2D.h:316
HistoBin2D & bin(size_t index)
Access a bin by index (non-const version)
Definition: Histo2D.h:294
double heightErr() const
Error on height.
Definition: HistoBin2D.h:111
double relErr() const
The relative size of the error (same for either volume or height errors)
Definition: HistoBin2D.h:116
Points & points()
Get the collection of points (non-const)
Definition: Scatter3D.h:209
A two-dimensional histogram.
Definition: Histo2D.h:32
Histo2D(const std::string &path="", const std::string &title="")
Default constructor.
Definition: Histo2D.h:50
double xRMS(bool includeoverflows=true) const
Get the RMS in x.
Definition: Histo2D.cc:170
double sumW(bool includeoverflows=true) const
Get the sum of weights in histo.
Definition: Histo2D.cc:103
size_t numBins() const
Number of bins.
Definition: Histo2D.h:311
double yMid() const
Middle of the bin in y.
Definition: Bin2D.h:182
void setZ(double z)
Set z value.
Definition: Point3D.h:115
void addPoint(const Point3D &pt)
Insert a new point.
Definition: Scatter3D.h:240
double numEntries() const
Number of entries (number of times fill was called, ignoring weights)
Definition: Dbn2D.h:166
Scatter1D efficiency(const Counter &accepted, const Counter &total)
Calculate an efficiency ratio of two counters.
Definition: Counter.cc:40