yoda is hosted by Hepforge, IPPP Durham
YODA - Yet more Objects for Data Analysis  1.7.2
Profile2D.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-2017 The YODA collaboration (see AUTHORS for details)
5 //
6 #include "YODA/Profile2D.h"
7 #include "YODA/Scatter3D.h"
8 #include "YODA/Histo2D.h"
9 
10 using namespace std;
11 
12 namespace YODA {
13 
14 
15  void Profile2D::fill(double x, double y, double z, double weight, double fraction) {
16  if ( std::isnan(x) ) throw RangeError("X is NaN");
17  if ( std::isnan(y) ) throw RangeError("Y is NaN");
18  if ( std::isnan(z) ) throw RangeError("Z is NaN");
19 
20  // Fill the overall distribution
21  _axis.totalDbn().fill(x, y, z, weight, fraction);
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, z, weight, fraction);
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, z, weight, fraction);
37  // }
38 
39  // Lock the axis now that a fill has happened
40  _axis._setLock(true);
41  }
42 
43 
44  void Profile2D::fillBin(size_t i, double z, double weight, double fraction) {
45  pair<double, double> mid = bin(i).xyMid();
46  fill(mid.first, mid.second, z, weight, fraction);
47  }
48 
49 
51 
52  double Profile2D::numEntries(bool includeoverflows) const {
53  if (includeoverflows) return totalDbn().numEntries();
54  unsigned long n = 0;
55  for (const Bin& b : bins()) n += b.numEntries();
56  return n;
57  }
58 
59 
60  double Profile2D::effNumEntries(bool includeoverflows) const {
61  if (includeoverflows) return totalDbn().effNumEntries();
62  double n = 0;
63  for (const Bin& b : bins()) n += b.effNumEntries();
64  return n;
65  }
66 
67 
68  double Profile2D::sumW(bool includeoverflows) const {
69  if (includeoverflows) return _axis.totalDbn().sumW2();
70  double sumw = 0;
71  for (const Bin& b : bins()) sumw += b.sumW();
72  return sumw;
73  }
74 
75 
76  double Profile2D::sumW2(bool includeoverflows) const {
77  if (includeoverflows) return _axis.totalDbn().sumW2();
78  double sumw2 = 0;
79  for (const Bin& b : bins()) sumw2 += b.sumW2();
80  return sumw2;
81  }
82 
83  // ^^^^^^^^^^^
84 
85 
86  double Profile2D::xMean(bool includeoverflows) const {
87  if (includeoverflows) return _axis.totalDbn().xMean();
88  Dbn3D dbn;
89  for (const ProfileBin2D& b : bins()) dbn += b.dbn();
90  return dbn.xMean();
91  }
92 
93 
94  double Profile2D::yMean(bool includeoverflows) const {
95  if (includeoverflows) return _axis.totalDbn().yMean();
96  Dbn3D dbn;
97  for (const ProfileBin2D& b : bins()) dbn += b.dbn();
98  return dbn.yMean();
99  }
100 
101 
102  double Profile2D::xVariance(bool includeoverflows) const {
103  if (includeoverflows) return _axis.totalDbn().xVariance();
104  Dbn3D dbn;
105  for (const ProfileBin2D& b : bins()) dbn += b.dbn();
106  return dbn.xVariance();
107  }
108 
109 
110  double Profile2D::yVariance(bool includeoverflows) const {
111  if (includeoverflows) return _axis.totalDbn().yVariance();
112  Dbn3D dbn;
113  for (const ProfileBin2D& b : bins()) dbn += b.dbn();
114  return dbn.yVariance();
115  }
116 
117 
118  double Profile2D::xStdErr(bool includeoverflows) const {
119  if (includeoverflows) return _axis.totalDbn().xStdErr();
120  Dbn3D dbn;
121  for (const ProfileBin2D& b : bins()) dbn += b.dbn();
122  return dbn.xStdErr();
123  }
124 
125 
126  double Profile2D::yStdErr(bool includeoverflows) const {
127  if (includeoverflows) return _axis.totalDbn().yStdErr();
128  Dbn3D dbn;
129  for (const ProfileBin2D& b : bins()) dbn += b.dbn();
130  return dbn.yStdErr();
131  }
132 
133 
134  double Profile2D::xRMS(bool includeoverflows) const {
135  if (includeoverflows) return _axis.totalDbn().xRMS();
136  Dbn3D dbn;
137  for (const ProfileBin2D& b : bins()) dbn += b.dbn();
138  return dbn.xRMS();
139  }
140 
141 
142  double Profile2D::yRMS(bool includeoverflows) const {
143  if (includeoverflows) return _axis.totalDbn().yRMS();
144  Dbn3D dbn;
145  for (const ProfileBin2D& b : bins()) dbn += b.dbn();
146  return dbn.yRMS();
147  }
148 
149 
151 
152 
154  Profile2D::Profile2D(const Profile2D& p, const std::string& path)
155  : AnalysisObject("Profile2D",
156  (path.size() == 0) ? p.path() : path,
157  p, p.title()),
158  _axis(p._axis)
159  { }
160 
161 
163  Profile2D::Profile2D(const Scatter3D& s, const std::string& path)
164  : AnalysisObject("Profile2D",
165  (path.size() == 0) ? s.path() : path,
166  s, s.title())
167  {
168  Bins bins;
169  for (const Scatter3D::Point& p : s.points()) {
170  bins.push_back(ProfileBin2D(p.xMin(), p.yMin(), p.xMax(), p.yMax()));
171  }
172  _axis = Profile2DAxis(bins);
173  }
174 
175 
177  Profile2D::Profile2D(const Histo2D& h, const std::string& path)
178  : AnalysisObject("Profile2D", (path.size() == 0) ? h.path() : path, h, h.title())
179  {
180  Bins bins;
181  for (const HistoBin2D& b : h.bins()) {
182  bins.push_back(ProfileBin2D(b.xMin(), b.yMin(), b.xMax(), b.yMax()));
183  }
184  _axis = Profile2DAxis(bins);
185  }
186 
187 
189  Scatter3D divide(const Profile2D& numer, const Profile2D& denom) {
190  Scatter3D rtn;
191 
192  for (size_t i = 0; i < numer.numBins(); ++i) {
193  const ProfileBin2D& b1 = numer.bin(i);
194  const ProfileBin2D& b2 = denom.bin(i);
195 
197  if (!fuzzyEquals(b1.xMin(), b2.xMin()) || !fuzzyEquals(b1.xMax(), b2.xMax()))
198  throw BinningError("x binnings are not equivalent in " + numer.path() + " / " + denom.path());
199  if (!fuzzyEquals(b1.yMin(), b2.yMin()) || !fuzzyEquals(b1.yMax(), b2.yMax()))
200  throw BinningError("y binnings are not equivalent in " + numer.path() + " / " + denom.path());
201 
202  // Assemble the x value and error
203  // Use the midpoint of the "bin" for the new central x value, in the absence of better information
204  const double x = b1.xMid();
205  const double exminus = x - b1.xMin();
206  const double explus = b1.xMax() - x;
207 
208  // Assemble the y value and error
209  // Use the midpoint of the "bin" for the new central y value, in the absence of better information
210  const double y = b1.yMid();
211  const double eyminus = y - b1.yMin();
212  const double eyplus = b1.yMax() - y;
213 
214  // Assemble the z value and error
215  double z = std::numeric_limits<double>::quiet_NaN();
216  double ez = std::numeric_limits<double>::quiet_NaN();
217  try {
218  if (b2.mean() == 0 || (b1.mean() == 0 && b1.stdErr() != 0)) {
219  // z = std::numeric_limits<double>::quiet_NaN();
220  // ez = std::numeric_limits<double>::quiet_NaN();
221  // throw LowStatsError("Requested division of empty bin");
222  } else {
223  z = b1.mean() / b2.mean();
225  const double relerr_1 = b1.stdErr() != 0 ? b1.stdErr()/b1.mean() : 0;
226  const double relerr_2 = b2.stdErr() != 0 ? b2.stdErr()/b2.mean() : 0;
227  ez = z * sqrt(sqr(relerr_1) + sqr(relerr_2));
228  }
229  } catch (const LowStatsError& e) {
230  // z = std::numeric_limits<double>::quiet_NaN();
231  // ez = std::numeric_limits<double>::quiet_NaN();
232  }
233 
236  //const double eyplus = y * sqrt( sqr(p1.yErrPlus()/p1.y()) + sqr(p2.yErrMinus()/p2.y()) );
237  //const double eyminus = y * sqrt( sqr(p1.yErrMinus()/p1.y()) + sqr(p2.yErrPlus()/p2.y()) );
238  rtn.addPoint(x, y, z, exminus, explus, eyminus, eyplus, ez, ez);
239  }
240 
241  assert(rtn.numPoints() == numer.numBins());
242  return rtn;
243  }
244 
245 
247 
248 
249 }
double yMean() const
Weighted mean, , of distribution.
Definition: Dbn3D.h:186
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:155
double yVariance() const
Weighted variance, , of distribution.
Definition: Dbn3D.h:195
A Bin1D specialised for handling profile-type information.
Definition: ProfileBin2D.h:23
double yMax() const
Upper y limit of the bin (exclusive).
Definition: Bin2D.h:163
std::vector< YODA::HistoBin2D > & bins()
Access the bin vector (non-const version)
Definition: Histo2D.h:260
double xMin() const
Lower x limit of the bin (inclusive).
Definition: Bin2D.h:142
Scatter1D divide(const Counter &numer, const Counter &denom)
Definition: Counter.cc:25
const std::string path() const
Get the AO path.
Profile2D(const std::string &path="", const std::string &title="")
Default constructor.
Definition: Profile2D.h:49
double stdErr() const
The standard error on the mean.
Definition: ProfileBin2D.h:134
double yMin() const
Lower y limit of the bin (inclusive).
Definition: Bin2D.h:158
STL namespace.
Error for e.g. use of invalid bin ranges.
Definition: Exceptions.h:34
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:72
double xMax() const
Upper x limit of the bin (exclusive).
Definition: Bin2D.h:147
double xMean() const
Weighted mean, , of distribution.
Definition: Dbn3D.h:183
Axis::Bins Bins
Definition: Profile2D.h:36
double xMid() const
Middle of the bin in x.
Definition: Bin2D.h:169
size_t numPoints() const
Number of points in the scatter.
Definition: Scatter3D.h:196
ProfileBin2D & bin(size_t index)
Access a bin by index (non-const)
Definition: Profile2D.h:252
double xRMS() const
Weighted RMS, , of distribution.
Definition: Dbn3D.h:219
double yRMS() const
Weighted RMS, , of distribution.
Definition: Dbn3D.h:222
std::vector< YODA::ProfileBin2D > & bins()
Access the bin vector (non-const)
Definition: Profile2D.h:245
NUM sqr(NUM a)
Named number-type squaring operation.
Definition: MathUtils.h:207
A two-dimensional profile histogram.
Definition: Profile2D.h:31
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 mean() const
The mean of the z distribution.
Definition: ProfileBin2D.h:119
Error for general binning problems.
Definition: Exceptions.h:27
double xStdErr() const
Weighted standard error on the mean, , of distribution.
Definition: Dbn3D.h:210
Axis2D< ProfileBin2D, Dbn3D > Profile2DAxis
Convenience typedef.
Definition: Profile2D.h:24
size_t numBins() const
Number of bins of this axis (not counting under/over flow)
Definition: Profile2D.h:268
double yStdErr() const
Weighted standard error on the mean, , of distribution.
Definition: Dbn3D.h:213
const std::string title() const
Get the AO title.
A very generic data type which is just a collection of 3D data points with errors.
Definition: Scatter3D.h:24
A Bin2D specialised for handling histogram-type information.
Definition: HistoBin2D.h:21
AnalysisObject is the base class for histograms and scatters.
Points & points()
Get the collection of points (non-const)
Definition: Scatter3D.h:202
A 2D distribution.
Definition: Dbn3D.h:16
A two-dimensional histogram.
Definition: Histo2D.h:31
double yMid() const
Middle of the bin in y.
Definition: Bin2D.h:174
double xVariance() const
Weighted variance, , of distribution.
Definition: Dbn3D.h:192
void addPoint(const Point3D &pt)
Insert a new point.
Definition: Scatter3D.h:233