ThePEG  1.8.0
Histogram1D.h
1 // -*- C++ -*-
2 //
3 // Histogram1D.h is a part of ThePEG - Toolkit for HEP Event Generation
4 // Copyright (C) 1999-2011 Leif Lonnblad
5 //
6 // ThePEG is licenced under version 2 of the GPL, see COPYING for details.
7 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
8 //
9 #ifndef LWH_Histogram1D_H
10 #define LWH_Histogram1D_H
11 //
12 // This is the declaration of the Histogram1D class.
13 //
14 
15 #include "AIHistogram1D.h"
16 #include "ManagedObject.h"
17 #include "Axis.h"
18 #include "VariAxis.h"
19 #include <vector>
20 #include <stdexcept>
21 
22 namespace LWH {
23 
24 using namespace AIDA;
25 
29 class Histogram1D: public IHistogram1D, public ManagedObject {
30 
31 public:
32 
34  friend class HistogramFactory;
35 
36 public:
37 
41  Histogram1D(int n, double lo, double up)
42  : fax(new Axis(n, lo, up)), vax(0),
43  sum(n + 2), sumw(n + 2), sumw2(n + 2), sumxw(n + 2), sumx2w(n + 2) {
44  ax = fax;
45  }
46 
50  Histogram1D(const std::vector<double> & edges)
51  : fax(0), vax(new VariAxis(edges)),
52  sum(edges.size() + 1), sumw(edges.size() + 1), sumw2(edges.size() + 1),
53  sumxw(edges.size() + 1), sumx2w(edges.size() + 1) {
54  ax = vax;
55  }
56 
61  : IBaseHistogram(h), IHistogram(h), IHistogram1D(h), ManagedObject(h),
62  fax(0), vax(0), sum(h.sum), sumw(h.sumw), sumw2(h.sumw2),
63  sumxw(h.sumxw), sumx2w(h.sumx2w) {
64  const VariAxis * hvax = dynamic_cast<const VariAxis *>(h.ax);
65  if ( vax ) ax = vax = new VariAxis(*hvax);
66  else ax = fax = new Axis(dynamic_cast<const Axis &>(*h.ax));
67 }
68 
70  virtual ~Histogram1D() {
71  delete ax;
72  }
73 
78  std::string title() const {
79  return theTitle;
80  }
81 
86  std::string name() const {
87  return theTitle;
88  }
89 
95  bool setTitle(const std::string & title) {
96  theTitle = title;
97  return true;
98  }
99 
103  IAnnotation & annotation() {
104  throw std::runtime_error("LWH cannot handle annotations");
105  }
106 
110  const IAnnotation & annotation() const {
111  throw std::runtime_error("LWH cannot handle annotations");
112  }
113 
118  int dimension() const {
119  return 1;
120  }
121 
126  bool reset() {
127  sum = std::vector<int>(ax->bins() + 2);
128  sumw = std::vector<double>(ax->bins() + 2);
129  sumxw = std::vector<double>(ax->bins() + 2);
130  sumx2w = std::vector<double>(ax->bins() + 2);
131  sumw2 = std::vector<double>(ax->bins() + 2);
132  return true;
133  }
134 
140  int entries() const {
141  int si = 0;
142  for ( int i = 2; i < ax->bins() + 2; ++i ) si += sum[i];
143  return si;
144  }
145 
153  int allEntries() const {
154  return entries() + extraEntries();
155  }
156 
161  int extraEntries() const {
162  return sum[0] + sum[1];
163  }
164 
170  double equivalentBinEntries() const {
171  double sw = 0.0;
172  double sw2 = 0.0;
173  for ( int i = 2; i < ax->bins() + 2; ++i ) {
174  sw += sumw[i];
175  sw2 += sumw2[i];
176  }
177  return sw2/(sw*sw);
178  }
179 
186  double sumBinHeights() const {
187  double sw = 0.0;
188  for ( int i = 2; i < ax->bins() + 2; ++i ) sw += sumw[i];
189  return sw;
190  }
191 
197  double sumAllBinHeights() const {
198  return sumBinHeights() + sumExtraBinHeights();
199  }
200 
205  double sumExtraBinHeights() const {
206  return sumw[0] + sumw[1];
207  }
208 
214  double minBinHeight() const {
215  double minw = sumw[2];
216  for ( int i = 3; i < ax->bins() + 2; ++i ) minw = std::min(minw, sumw[i]);
217  return minw;
218  }
219 
225  double maxBinHeight() const{
226  double maxw = sumw[2];
227  for ( int i = 3; i < ax->bins() + 2; ++i ) maxw = std::max(maxw, sumw[i]);
228  return maxw;
229  }
230 
238  bool fill(double x, double weight = 1.) {
239  int i = ax->coordToIndex(x) + 2;
240  ++sum[i];
241  sumw[i] += weight;
242  sumxw[i] += x*weight;
243  sumx2w[i] += x*x*weight;
244  sumw2[i] += weight*weight;
245  return weight >= 0 && weight <= 1;
246  }
247 
253  double binMean(int index) const {
254  int i = index + 2;
255  return sumw[i] != 0.0? sumxw[i]/sumw[i]:
256  ( vax? vax->binMidPoint(index): fax->binMidPoint(index) );
257  };
258 
264  double binRms(int index) const {
265  int i = index + 2;
266  return sumw[i] == 0.0 || sum[i] < 2? ax->binWidth(index):
267  std::sqrt(std::max(sumw[i]*sumx2w[i] - sumxw[i]*sumxw[i], 0.0))/sumw[i];
268  };
269 
276  int binEntries(int index) const {
277  return sum[index + 2];
278  }
279 
286  double binHeight(int index) const {
287  return sumw[index + 2];
288  }
289 
296  double binError(int index) const {
297  return std::sqrt(sumw2[index + 2]);
298  }
299 
304  double mean() const {
305  double s = 0.0;
306  double sx = 0.0;
307  for ( int i = 2; i < ax->bins() + 2; ++i ) {
308  s += sumw[i];
309  sx += sumxw[i];
310  }
311  return s != 0.0? sx/s: 0.0;
312  }
313 
318  double rms() const {
319  double s = 0.0;
320  double sx = 0.0;
321  double sx2 = 0.0;
322  for ( int i = 2; i < ax->bins() + 2; ++i ) {
323  s += sumw[i];
324  sx += sumxw[i];
325  sx2 += sumx2w[i];
326  }
327  return s != 0.0? std::sqrt(std::max(s*sx2 - sx*sx, 0.0))/s:
328  ax->upperEdge() - ax->lowerEdge();
329  }
330 
335  const IAxis & axis() const {
336  return *ax;
337  }
338 
346  int coordToIndex(double coord) const {
347  return ax->coordToIndex(coord);
348  }
349 
355  bool add(const Histogram1D & h) {
356  if ( ax->upperEdge() != h.ax->upperEdge() ||
357  ax->lowerEdge() != h.ax->lowerEdge() ||
358  ax->bins() != h.ax->bins() ) return false;
359  for ( int i = 0; i < ax->bins() + 2; ++i ) {
360  sum[i] += h.sum[i];
361  sumw[i] += h.sumw[i];
362  sumxw[i] += h.sumxw[i];
363  sumx2w[i] += h.sumx2w[i];
364  sumw2[i] += h.sumw2[i];
365  }
366  return true;
367  }
368 
374  bool add(const IHistogram1D & hist) {
375  return add(dynamic_cast<const Histogram1D &>(hist));
376  }
377 
382  bool scale(double s) {
383  for ( int i = 0; i < ax->bins() + 2; ++i ) {
384  sumw[i] *= s;
385  sumxw[i] *= s;
386  sumx2w[i] *= s;
387  sumw2[i] *= s*s;
388  }
389  return true;
390  }
391 
399  void normalize(double intg) {
400  double oldintg = sumAllBinHeights();
401  if ( oldintg == 0.0 ) return;
402  for ( int i = 0; i < ax->bins() + 2; ++i ) {
403  double fac = intg/oldintg;
404  if ( i >= 2 ) fac /= (ax->binUpperEdge(i - 2) - ax->binLowerEdge(i - 2));
405  sumw[i] *= fac;
406  sumxw[i] *= fac;
407  sumx2w[i] *= fac;
408  sumw2[i] *= fac*fac;
409  }
410  }
411 
416  double integral() const {
417  double intg = sumw[0] + sumw[1];
418  for ( int i = 2; i < ax->bins() + 2; ++i )
419  intg += sumw[i]*(ax->binUpperEdge(i - 2) - ax->binLowerEdge(i - 2));
420  return intg;
421  }
422 
427  void * cast(const std::string &) const {
428  return 0;
429  }
430 
434  bool writeXML(std::ostream & os, std::string path, std::string name) {
435  os << " <histogram1d name=\"" << name
436  << "\"\n title=\"" << title()
437  << "\" path=\"" << path
438  << "\">\n <axis max=\"" << ax->upperEdge()
439  << "\" numberOfBins=\"" << ax->bins()
440  << "\" min=\"" << ax->lowerEdge()
441  << "\" direction=\"x\"";
442  if ( vax ) {
443  os << ">\n";
444  for ( int i = 0, N = ax->bins() - 1; i < N; ++i )
445  os << " <binBorder value=\"" << ax->binUpperEdge(i) << "\"/>\n";
446  os << " </axis>\n";
447  } else {
448  os << "/>\n";
449  }
450  os << " <statistics entries=\"" << entries()
451  << "\">\n <statistic mean=\"" << mean()
452  << "\" direction=\"x\"\n rms=\"" << rms()
453  << "\"/>\n </statistics>\n <data1d>\n";
454  for ( int i = 0; i < ax->bins() + 2; ++i ) if ( sum[i] ) {
455  os << " <bin1d binNum=\"";
456  if ( i == 0 ) os << "UNDERFLOW";
457  else if ( i == 1 ) os << "OVERFLOW";
458  else os << i - 2;
459  os << "\" entries=\"" << sum[i]
460  << "\" height=\"" << sumw[i]
461  << "\"\n error=\"" << std::sqrt(sumw2[i])
462  << "\" error2=\"" << sumw2[i]
463  << "\"\n weightedMean=\"" << binMean(i - 2)
464  << "\" weightedRms=\"" << binRms(i - 2)
465  << "\"/>\n";
466  }
467  os << " </data1d>\n </histogram1d>" << std::endl;
468  return true;
469  }
470 
471 
476  bool writeFLAT(std::ostream & os, std::string path, std::string name) {
477  os << "# " << path << "/" << name << " " << ax->lowerEdge()
478  << " " << ax->bins() << " " << ax->upperEdge()
479  << " \"" << title() << " \"" << std::endl;
480  for ( int i = 2; i < ax->bins() + 2; ++i )
481  os << 0.5*(ax->binLowerEdge(i - 2) + ax->binUpperEdge(i - 2)) << " "
482  << sumw[i] << " " << sqrt(sumw2[i]) << " " << sum[i] << std::endl;
483  os << std::endl;
484  return true;
485  }
486 
487 private:
488 
490  std::string theTitle;
491 
493  IAxis * ax;
494 
497 
500 
502  std::vector<int> sum;
503 
505  std::vector<double> sumw;
506 
508  std::vector<double> sumw2;
509 
511  std::vector<double> sumxw;
512 
514  std::vector<double> sumx2w;
515 
517  IAnnotation * anno;
518 
519 };
520 
521 }
522 
523 #endif /* LWH_Histogram1D_H */