yoda is hosted by Hepforge, IPPP Durham
YODA - Yet more Objects for Data Analysis  1.7.2
Profile1D.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/Profile1D.h"
7 #include "YODA/Histo1D.h"
8 #include "YODA/Scatter2D.h"
9 
10 #include <cmath>
11 using namespace std;
12 
13 namespace YODA {
14 
15 
16  void Profile1D::fill(double x, double y, double weight, double fraction) {
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, fraction);
22 
23  // Fill the bins and overflows
25  if (inRange(x, _axis.xMin(), _axis.xMax())) {
26  try {
28  _binAt(x).fill(x, y, weight, fraction);
29  } catch (const RangeError& re) { }
30  } else if (x < _axis.xMin()) {
31  _axis.underflow().fill(x, y, weight, fraction);
32  } else if (x >= _axis.xMax()) {
33  _axis.overflow().fill(x, y, weight, fraction);
34  }
35 
36  // Lock the axis now that a fill has happened
37  _axis._setLock(true);
38  }
39 
40 
41  void Profile1D::fillBin(size_t i, double y, double weight, double fraction) {
42  fill(bin(i).xMid(), y, weight, fraction);
43  }
44 
45 
46 
48 
49  double Profile1D::numEntries(bool includeoverflows) const {
50  if (includeoverflows) return totalDbn().numEntries();
51  unsigned long n = 0;
52  for (const Bin& b : bins()) n += b.numEntries();
53  return n;
54  }
55 
56 
57  double Profile1D::effNumEntries(bool includeoverflows) const {
58  if (includeoverflows) return totalDbn().effNumEntries();
59  double n = 0;
60  for (const Bin& b : bins()) n += b.effNumEntries();
61  return n;
62  }
63 
64 
65  double Profile1D::sumW(bool includeoverflows) const {
66  if (includeoverflows) return _axis.totalDbn().sumW();
67  double sumw = 0;
68  for (const Bin& b : bins()) sumw += b.sumW();
69  return sumw;
70  }
71 
72 
73  double Profile1D::sumW2(bool includeoverflows) const {
74  if (includeoverflows) return _axis.totalDbn().sumW2();
75  double sumw2 = 0;
76  for (const Bin& b : bins()) sumw2 += b.sumW2();
77  return sumw2;
78  }
79 
80  // ^^^^^^^^^^^^^^^^^^
81 
82 
83  double Profile1D::xMean(bool includeoverflows) const {
84  if (includeoverflows) return _axis.totalDbn().xMean();
85  Dbn2D dbn;
86  for (const ProfileBin1D& b : bins()) dbn += b.dbn();
87  return dbn.xMean();
88  }
89 
90 
91  double Profile1D::xVariance(bool includeoverflows) const {
92  if (includeoverflows) return _axis.totalDbn().xVariance();
93  Dbn2D dbn;
94  for (const ProfileBin1D& b : bins()) dbn += b.dbn();
95  return dbn.xVariance();
96  }
97 
98 
99  double Profile1D::xStdErr(bool includeoverflows) const {
100  if (includeoverflows) return _axis.totalDbn().xStdErr();
101  Dbn2D dbn;
102  for (const ProfileBin1D& b : bins()) dbn += b.dbn();
103  return dbn.xStdErr();
104  }
105 
106 
107  double Profile1D::xRMS(bool includeoverflows) const {
108  if (includeoverflows) return _axis.totalDbn().xRMS();
109  Dbn2D dbn;
110  for (const ProfileBin1D& b : bins()) dbn += b.dbn();
111  return dbn.xRMS();
112  }
113 
114 
116 
117 
119  Profile1D::Profile1D(const Profile1D& p, const std::string& path)
120  : AnalysisObject("Profile1D", (path.size() == 0) ? p.path() : path, p, p.title())
121  {
122  _axis = p._axis;
123  }
124 
125 
127  Profile1D::Profile1D(const Scatter2D& s, const std::string& path)
128  : AnalysisObject("Profile1D", (path.size() == 0) ? s.path() : path, s, s.title())
129  {
130  std::vector<ProfileBin1D> bins;
131  for (const Scatter2D::Point& p : s.points()) {
132  bins.push_back(ProfileBin1D(p.xMin(), p.xMax()));
133  }
134  _axis = Profile1DAxis(bins);
135  }
136 
137 
139  Profile1D::Profile1D(const Histo1D& h, const std::string& path)
140  : AnalysisObject("Profile1D", (path.size() == 0) ? h.path() : path, h, h.title())
141  {
142  Bins bins;
143  for (const Histo1D::Bin& b : h.bins()) {
144  bins.push_back(ProfileBin1D(b.xMin(), b.xMax()));
145  }
146  _axis = Profile1DAxis(bins);
147 
148  }
149 
150 
152 
153 
155  Scatter2D divide(const Profile1D& numer, const Profile1D& denom) {
156  Scatter2D rtn;
157 
158  for (size_t i = 0; i < numer.numBins(); ++i) {
159  const ProfileBin1D& b1 = numer.bin(i);
160  const ProfileBin1D& b2 = denom.bin(i);
161 
163  if (!fuzzyEquals(b1.xMin(), b2.xMin()) || !fuzzyEquals(b1.xMax(), b2.xMax()))
164  throw BinningError("x binnings are not equivalent in " + numer.path() + " / " + denom.path());
165 
166  // Assemble the x value and error
167  // Use the midpoint of the "bin" for the new central x value, in the absence of better information
168  const double x = b1.xMid();
169  const double exminus = x - b1.xMin();
170  const double explus = b1.xMax() - x;
171 
172  // Assemble the y value and error
173  double y = std::numeric_limits<double>::quiet_NaN();
174  double ey = std::numeric_limits<double>::quiet_NaN();
175  try {
176  if (b2.mean() == 0 || (b1.mean() == 0 && b1.stdErr() != 0)) {
177  // y = std::numeric_limits<double>::quiet_NaN();
178  // ey = std::numeric_limits<double>::quiet_NaN();
179  // throw LowStatsError("Requested division of empty bin");
180  } else {
181  y = b1.mean() / b2.mean();
183  const double relerr_1 = b1.stdErr() != 0 ? b1.stdErr()/b1.mean() : 0;
184  const double relerr_2 = b2.stdErr() != 0 ? b2.stdErr()/b2.mean() : 0;
185  ey = y * sqrt(sqr(relerr_1) + sqr(relerr_2));
186  }
187  } catch (const LowStatsError& e) {
188  // y = std::numeric_limits<double>::quiet_NaN();
189  // ey = std::numeric_limits<double>::quiet_NaN();
190  }
191 
194  //const double eyplus = y * sqrt( sqr(p1.yErrPlus()/p1.y()) + sqr(p2.yErrMinus()/p2.y()) );
195  //const double eyminus = y * sqrt( sqr(p1.yErrMinus()/p1.y()) + sqr(p2.yErrPlus()/p2.y()) );
196  rtn.addPoint(x, y, exminus, explus, ey, ey);
197  }
198 
199  assert(rtn.numPoints() == numer.numBins());
200  return rtn;
201  }
202 
203 
205 
206 
207 }
Axis1D< ProfileBin1D, Dbn2D > Profile1DAxis
Convenience typedef.
Definition: Profile1D.h:25
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 xVariance() const
Weighted variance, , of distribution.
Definition: Dbn2D.h:136
Scatter1D divide(const Counter &numer, const Counter &denom)
Definition: Counter.cc:25
const std::string path() const
Get the AO path.
void addPoint(const Point2D &pt)
Insert a new point.
Definition: Scatter2D.h:216
A very generic data type which is just a collection of 2D data points with errors.
Definition: Scatter2D.h:24
double xMean() const
Weighted mean, , of distribution.
Definition: Dbn2D.h:130
size_t numBins() const
Number of bins on this axis (not counting under/overflow)
Definition: Profile1D.h:219
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:72
std::vector< YODA::HistoBin1D > & bins()
Access the bin vector.
Definition: Histo1D.h:211
A one-dimensional histogram.
Definition: Histo1D.h:26
A Bin1D specialised for handling histogram-type information.
Definition: HistoBin1D.h:21
A 2D distribution.
Definition: Dbn2D.h:16
double xMid() const
Geometric centre of the bin, i.e. high+low/2.0.
Definition: Bin1D.h:136
NUM sqr(NUM a)
Named number-type squaring operation.
Definition: MathUtils.h:207
size_t numPoints() const
Number of points in the scatter.
Definition: Scatter2D.h:179
Errors relating to insufficient (effective) statistics.
Definition: Exceptions.h:72
Error for general binning problems.
Definition: Exceptions.h:27
std::vector< YODA::ProfileBin1D > & bins()
Access the bin vector.
Definition: Profile1D.h:235
double stdErr() const
The standard error on the mean.
Definition: ProfileBin1D.h:126
double xRMS() const
Weighted RMS, , of distribution.
Definition: Dbn2D.h:154
const std::string title() const
Get the AO title.
A Bin1D specialised for handling profile-type information.
Definition: ProfileBin1D.h:23
AnalysisObject is the base class for histograms and scatters.
double xMin() const
Lower limit of the bin (inclusive).
Definition: Bin1D.h:126
double mean() const
The mean of the y distribution.
Definition: ProfileBin1D.h:111
Profile1D(const std::string &path="", const std::string &title="")
Default constructor.
Definition: Profile1D.h:50
A one-dimensional profile histogram.
Definition: Profile1D.h:33
double xMax() const
Upper limit of the bin (exclusive).
Definition: Bin1D.h:131
ProfileBin1D & bin(size_t index)
Access a bin by index (non-const version)
Definition: Profile1D.h:242
Axis::Bins Bins
Definition: Profile1D.h:38
Points & points()
Get the collection of points (non-const)
Definition: Scatter2D.h:185
A 2D data point to be contained in a Scatter2D.
Definition: Point2D.h:18