yoda is hosted by Hepforge, IPPP Durham
YODA - Yet more Objects for Data Analysis  1.9.0
Scatter2D.cc
Go to the documentation of this file.
1 #include "YODA/Scatter2D.h"
2 #include "YODA/Histo1D.h"
3 #include "YODA/Profile1D.h"
4 #include <sstream>
5 #include "yaml-cpp/yaml.h"
6 #ifdef YAML_NAMESPACE
7 #define YAML YAML_NAMESPACE
8 #endif
9 
10 namespace YODA {
11 
12 
14  Scatter2D mkScatter(const Histo1D& h, bool usefocus, bool binwidthdiv) {
15  Scatter2D rtn;
16  for (const std::string& a : h.annotations()) rtn.setAnnotation(a, h.annotation(a));
17  rtn.setAnnotation("Type", h.type()); // might override the copied ones
18 
19  for (const HistoBin1D& b : h.bins()) {
20  const double x = usefocus ? b.xFocus() : b.xMid();
21  const double ex_m = x - b.xMin();
22  const double ex_p = b.xMax() - x;
23 
24  double y;
25  try {
26  y = b.sumW();
27  } catch (const Exception&) { // LowStatsError or WeightError
28  y = std::numeric_limits<double>::quiet_NaN();
29  }
30  if (binwidthdiv) y /= b.xWidth();
31  const double ey = b.relErr() * y;
32 
33  // Attach the point to its parent
34  Point2D pt(x, y, ex_m, ex_p, ey, ey);
35  pt.setParent(&rtn);
36  rtn.addPoint(pt);
37  }
38 
39  assert(h.numBins() == rtn.numPoints());
40  return rtn;
41  }
42 
43 
45  Scatter2D mkScatter(const Profile1D& p, bool usefocus, bool usestddev) {
46  Scatter2D rtn;
47  for (const std::string& a : p.annotations())
48  rtn.setAnnotation(a, p.annotation(a));
49  rtn.setAnnotation("Type", p.type());
50  for (const ProfileBin1D& b : p.bins()) {
51  const double x = usefocus ? b.xFocus() : b.xMid();
52  const double ex_m = x - b.xMin();
53  const double ex_p = b.xMax() - x;
54 
55  double y;
56  try {
57  y = b.mean();
58  } catch (const Exception&) { // LowStatsError or WeightError
59  y = std::numeric_limits<double>::quiet_NaN();
60  }
61  double ey;
62  try {
63  ey = usestddev ? b.stdDev() : b.stdErr();
64  } catch (const Exception&) { // LowStatsError or WeightError
65  ey = std::numeric_limits<double>::quiet_NaN();
66  }
67 
68  //const Point2D pt(x, y, ex_m, ex_p, ey, ey);
69  Point2D pt(x, y, ex_m, ex_p, ey, ey);
70  pt.setParent(&rtn);
71  rtn.addPoint(pt);
72  }
73  assert(p.numBins() == rtn.numPoints());
74  return rtn;
75  }
76 
77 
78  // Retrieve variations from annotation, parse them as YAML, and update the points
80  if (this->_variationsParsed) { return; }
81  if (!(this->hasAnnotation("ErrorBreakdown"))) { return; }
82  YAML::Node errorBreakdown;
83  errorBreakdown = YAML::Load(this->annotation("ErrorBreakdown"));
84 
85  if (errorBreakdown.size()) {
86  for (size_t thisPointIndex = 0; thisPointIndex < this->numPoints(); ++thisPointIndex) {
87  Point2D& thispoint = this->_points[thisPointIndex];
88  YAML::Node variations = errorBreakdown[thisPointIndex];
89  for (const auto& variation : variations) {
90  const std::string variationName = variation.first.as<std::string>();
91  double eyp = variation.second["up"].as<double>();
92  double eym = variation.second["dn"].as<double>();
93  thispoint.setYErrs(eym,eyp,variationName);
94  }
95  }
96  this-> _variationsParsed =true;
97  }
98  }
99 
100 
102  std::vector<std::string> Scatter2D::variations() const {
103  std::vector<std::string> vecVariations;
104  for (auto& point : this->_points) {
105  for (auto& it : point.errMap()) {
106  // if the variation is not already in the vector, add it!
107  if (std::find(vecVariations.begin(), vecVariations.end(), it.first) == vecVariations.end()) {
108  vecVariations.push_back(it.first);
109  }
110  }
111  }
112  return vecVariations;
113  }
114 
115 
116  std::vector<std::vector<double> > Scatter2D::covarianceMatrix(bool ignoreOffDiagonalTerms) {
117  int nPoints = this->numPoints();
118  //double covM[nPoints][nPoints] = {};
119  std::vector<std::vector<double> > covM;
120 
121 
122  // initialise cov matrix to be the right shape
123  for (int i = 0; i < nPoints; i++) {
124  std::vector<double> row;
125  row.resize(nPoints);
126  covM.push_back(row);
127  }
128 
129  // case where only have nominal, ie total uncertainty, labelled "" (empty string)
130  if (this->variations().size() == 1) {
131  for (int i = 0; i < nPoints; i++) {
132  covM[i][i] = sqr(((this->_points[i].yErrs().first+this->_points[i].yErrs().second)/2));
133  if (covM[i][i] == 0 ) covM[i][i] = 1;
134  }
135  return covM;
136  }
137  // more interesting case where we actually have some uncertainty breakdown!
138  auto systList= this->variations();
139  for (auto sname : systList){
140  if (sname.length() == 0) continue;
141  std::vector<double> systErrs;
142  systErrs.resize(nPoints);
143  for (int i = 0; i < nPoints; i++) {
144  auto point = this->_points[i];
145  try {
146  auto variations = point.errMap().at(sname);
147  // up/dn are symmetrized since this method can't handle asymmetric errors
148  systErrs[i] = (fabs(variations.first)+fabs(variations.second))*0.5;
149  } catch (const std::exception& e) { // missing bin
150  systErrs[i] = 0.0;
151  }
152  }
153  if (ignoreOffDiagonalTerms || sname.find("stat") != std::string::npos || sname.find("uncor") != std::string::npos) {
154  for (int i = 0; i < nPoints; i++) {
155  covM[i][i] += systErrs[i]*systErrs[i]; // just the diagonal; bins are considered uncorrelated
156  }
157  } else {
158  for (int i = 0; i < nPoints; i++) {
159  for (int j = 0; j < nPoints; j++) {
160  covM[i][j] += systErrs[i]*systErrs[j];
161  }
162  }
163  }
164  }
165  return covM;
166  }
167 
168 
169 }
void parseVariations()
Definition: Scatter2D.cc:79
void setYErrs(double ey, std::string source="")
Set symmetric y error (alias)
Definition: Point2D.h:263
void addPoint(const Point2D &pt)
Insert a new point.
Definition: Scatter2D.h:225
A very generic data type which is just a collection of 2D data points with errors.
Definition: Scatter2D.h:25
size_t numBins() const
Number of bins (not counting under/overflow)
Definition: Profile1D.h:239
Generic unspecialised YODA runtime error.
Definition: Exceptions.h:20
std::vector< YODA::HistoBin1D > & bins()
Access the bin vector.
Definition: Histo1D.h:224
A one-dimensional histogram.
Definition: Histo1D.h:28
void setParent(Scatter *parent)
Definition: Point.h:114
A Bin1D specialised for handling histogram-type information.
Definition: HistoBin1D.h:21
const std::map< std::string, std::pair< double, double > > & errMap() const
Get error map for direction i.
Definition: Point2D.cc:8
NUM sqr(NUM a)
Named number-type squaring operation.
Definition: MathUtils.h:208
Point2D & point(size_t index)
Get a reference to the point with index index (non-const)
Definition: Scatter2D.h:206
size_t numPoints() const
Number of points in the scatter.
Definition: Scatter2D.h:188
const std::string & annotation(const std::string &name) const
Get an annotation by name (as a string)
std::vector< YODA::ProfileBin1D > & bins()
Access the bin vector.
Definition: Profile1D.h:258
Scatter1D mkScatter(const Counter &c)
Make a Scatter1D representation of a Histo1D.
Definition: Scatter1D.cc:13
A Bin1D specialised for handling profile-type information.
Definition: ProfileBin1D.h:23
size_t numBins() const
Number of bins (not counting under/overflow)
Definition: Histo1D.h:200
virtual std::string type() const
Get name of the analysis object type.
bool hasAnnotation(const std::string &name) const
Check if an annotation is defined.
std::vector< std::string > variations() const
Get the list of variations stored in the points.
Definition: Scatter2D.cc:102
std::vector< std::vector< double > > covarianceMatrix(bool ignoreOffDiagonalTerms=false)
Definition: Scatter2D.cc:116
A one-dimensional profile histogram.
Definition: Profile1D.h:35
void setAnnotation(const std::string &name, const std::string &value)
Add or set a string-valued annotation by name.
std::vector< std::string > annotations() const
A 2D data point to be contained in a Scatter2D.
Definition: Point2D.h:19