ThePEG  1.8.0
HistogramFactory.h
1 // -*- C++ -*-
2 //
3 // HistogramFactory.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_HistogramFactory_H
10 #define LWH_HistogramFactory_H
11 //
12 // This is the declaration of the HistogramFactory class.
13 //
14 
15 #include "AIHistogramFactory.h"
16 #include "Histogram1D.h"
17 #include "Histogram2D.h"
18 #include "Tree.h"
19 #include <string>
20 #include <stdexcept>
21 
22 namespace LWH {
23 
24 using namespace AIDA;
25 
32 class HistogramFactory: public IHistogramFactory {
33 
34 public:
35 
40  : tree(&t) {}
41 
45  virtual ~HistogramFactory() {}
46 
52  bool destroy(IBaseHistogram * hist) {
53  IManagedObject * mo = dynamic_cast<IManagedObject *>(hist);
54  if ( !mo ) return false;
55  return tree->rm(tree->findPath(*mo));
56  }
57 
61  ICloud1D * createCloud1D(const std::string &, const std::string &,
62  int = -1, const std::string & = "") {
63  return error<ICloud1D>("ICloud1D");
64  }
65 
69  ICloud1D * createCloud1D(const std::string &) {
70  return error<ICloud1D>("ICloud1D");
71  }
72 
76  ICloud1D * createCopy(const std::string &, const ICloud1D &) {
77  return error<ICloud1D>("ICloud1D");
78  }
79 
83  ICloud2D * createCloud2D(const std::string &, const std::string &, int = -1,
84  const std::string & = "") {
85  return error<ICloud2D>("ICloud2D");
86  }
87 
88 
92  ICloud2D * createCloud2D(const std::string &) {
93  return error<ICloud2D>("ICloud2D");
94  }
95 
99  ICloud2D * createCopy(const std::string &, const ICloud2D &) {
100  return error<ICloud2D>("ICloud2D");
101  }
102 
106  ICloud3D * createCloud3D(const std::string &, const std::string &, int = -1,
107  const std::string & = "") {
108  return error<ICloud3D>("ICloud3D");
109  }
110 
114  ICloud3D * createCloud3D(const std::string &) {
115  return error<ICloud3D>("ICloud3D");
116  }
117 
121  ICloud3D * createCopy(const std::string &, const ICloud3D &) {
122  return error<ICloud3D>("ICloud3D");
123  }
124 
141  IHistogram1D *
142  createHistogram1D(const std::string & path, const std::string & title,
143  int nBins, double lowerEdge, double upperEdge,
144  const std::string & = "") {
145  Histogram1D * hist = new Histogram1D(nBins, lowerEdge, upperEdge);
146  hist->setTitle(title);
147  if ( !tree->insert(path, hist) ) {
148  delete hist;
149  hist = 0;
150  throw std::runtime_error("LWH could not create histogram '"
151  + title + "'." );
152  }
153  return hist;
154  }
155 
168  IHistogram1D *
169  createHistogram1D(const std::string & pathAndTitle,
170  int nBins, double lowerEdge, double upperEdge) {
171  std::string title = pathAndTitle.substr(pathAndTitle.rfind('/') + 1);
172  return createHistogram1D(pathAndTitle, title, nBins, lowerEdge, upperEdge);
173  }
174 
175 
188  IHistogram1D *
189  createHistogram1D(const std::string & path, const std::string & title,
190  const std::vector<double> & binEdges,
191  const std::string & = "") {
192  Histogram1D * hist = new Histogram1D(binEdges);
193  hist->setTitle(title);
194  if ( !tree->insert(path, hist) ) {
195  delete hist;
196  hist = 0;
197  throw std::runtime_error("LWH could not create histogram '"
198  + title + "'." );
199  }
200  return hist;
201  }
202 
213  IHistogram1D *
214  createCopy(const std::string & path, const IHistogram1D & hist) {
215  Histogram1D * h = new Histogram1D(dynamic_cast<const Histogram1D &>(hist));
216  h->setTitle(path.substr(path.rfind('/') + 1));
217  if ( !tree->insert(path, h) ) {
218  delete h;
219  h = 0;
220  throw std::runtime_error("LWH could not create a copy of histogram '"
221  + hist.title() + "'." );
222  }
223  return h;
224  }
225 
229  IHistogram2D *
230  createHistogram2D(const std::string & path, const std::string & title,
231  int nx, double xlo, double xup,
232  int ny, double ylo, double yup,
233  const std::string & = "") {
234  Histogram2D * hist = new Histogram2D(nx, xlo, xup, ny, ylo, yup);
235  hist->setTitle(title);
236  if ( !tree->insert(path, hist) ) {
237  delete hist;
238  hist = 0;
239  throw std::runtime_error("LWH could not create histogram '"
240  + title + "'." );
241  }
242  return hist;
243  }
244 
248  IHistogram2D * createHistogram2D(const std::string & pathAndTitle,
249  int nx, double xlo, double xup,
250  int ny, double ylo, double yup) {
251  std::string title = pathAndTitle.substr(pathAndTitle.rfind('/') + 1);
252  return createHistogram2D(pathAndTitle, title, nx, xlo, xup, ny, ylo, yup);
253  }
254 
258  IHistogram2D *
259  createHistogram2D(const std::string & path, const std::string & title,
260  const std::vector<double> & xedges,
261  const std::vector<double> & yedges,
262  const std::string & = "") {
263  Histogram2D * hist = new Histogram2D(xedges, yedges);
264  hist->setTitle(title);
265  if ( !tree->insert(path, hist) ) {
266  delete hist;
267  hist = 0;
268  throw std::runtime_error("LWH could not create histogram '"
269  + title + "'." );
270  }
271  return hist;
272  }
273 
277  IHistogram2D *
278  createCopy(const std::string & path, const IHistogram2D & hist) {
279  Histogram2D * h = new Histogram2D(dynamic_cast<const Histogram2D &>(hist));
280  h->setTitle(path.substr(path.rfind('/') + 1));
281  if ( !tree->insert(path, h) ) {
282  delete h;
283  h = 0;
284  throw std::runtime_error("LWH could not create a copy of histogram '"
285  + hist.title() + "'." );
286  }
287  return h;
288  }
289 
293  IHistogram3D * createHistogram3D(const std::string &, const std::string &,
294  int, double, double, int, double, double,
295  int, double, double,
296  const std::string & = "") {
297  return error<IHistogram3D>("IHistogram3D");
298  }
299 
303  IHistogram3D * createHistogram3D(const std::string &, int, double, double,
304  int, double, double, int, double, double) {
305  return error<IHistogram3D>("IHistogram3D");
306  }
307 
311  IHistogram3D * createHistogram3D(const std::string &, const std::string &,
312  const std::vector<double> &,
313  const std::vector<double> &,
314  const std::vector<double> &,
315  const std::string & = "") {
316  return error<IHistogram3D>("IHistogram3D");
317  }
318 
322  IHistogram3D * createCopy(const std::string &, const IHistogram3D &) {
323  return error<IHistogram3D>("IHistogram3D");
324  }
325 
329  IProfile1D * createProfile1D(const std::string &, const std::string &,
330  int, double, double, const std::string & = "") {
331  return error<IProfile1D>("IProfile1D");
332  }
333 
337  IProfile1D * createProfile1D(const std::string &, const std::string &,
338  int, double, double, double, double,
339  const std::string & = "") {
340  return error<IProfile1D>("IProfile1D");
341  }
342 
346  IProfile1D * createProfile1D(const std::string &, const std::string &,
347  const std::vector<double> &,
348  const std::string & = "") {
349  return error<IProfile1D>("IProfile1D");
350  }
351 
355  IProfile1D * createProfile1D(const std::string &, const std::string &,
356  const std::vector<double> &, double, double,
357  const std::string & = "") {
358  return error<IProfile1D>("IProfile1D");
359  }
360 
364  IProfile1D * createProfile1D(const std::string &, int, double, double) {
365  return error<IProfile1D>("IProfile1D");
366  }
367 
371  IProfile1D * createProfile1D(const std::string &,
372  int, double, double, double, double) {
373  return error<IProfile1D>("IProfile1D");
374  }
375 
379  IProfile1D * createCopy(const std::string &, const IProfile1D &) {
380  return error<IProfile1D>("IProfile1D");
381  }
382 
386  IProfile2D * createProfile2D(const std::string &, const std::string &,
387  int, double, double, int, double, double,
388  const std::string & = "") {
389  return error<IProfile2D>("IProfile2D");
390  }
391 
395  IProfile2D * createProfile2D(const std::string &, const std::string &,
396  int, double, double, int,
397  double, double, double, double,
398  const std::string & = "") {
399  return error<IProfile2D>("IProfile2D");
400  }
401 
405  IProfile2D * createProfile2D(const std::string &, const std::string &,
406  const std::vector<double> &,
407  const std::vector<double> &,
408  const std::string & = "") {
409  return error<IProfile2D>("IProfile2D");
410  }
411 
415  IProfile2D * createProfile2D(const std::string &, const std::string &,
416  const std::vector<double> &,
417  const std::vector<double> &,
418  double, double, const std::string & = "") {
419  return error<IProfile2D>("IProfile2D");
420  }
421 
425  IProfile2D * createProfile2D(const std::string &, int, double, double,
426  int, double, double) {
427  return error<IProfile2D>("IProfile2D");
428  }
429 
433  IProfile2D * createProfile2D(const std::string &, int, double, double,
434  int, double, double, double, double) {
435  return error<IProfile2D>("IProfile2D");
436  }
437 
441  IProfile2D * createCopy(const std::string &, const IProfile2D &) {
442  return error<IProfile2D>("IProfile2D");
443  }
444 
456  Histogram1D * add(const std::string & path,
457  const Histogram1D & hist1, const Histogram1D & hist2) {
458  if ( !checkBins(hist1, hist2) ) return 0;
459  Histogram1D * h = new Histogram1D(hist1);
460  h->setTitle(path.substr(path.rfind('/') + 1));
461  h->add(hist2);
462  if ( !tree->insert(path, h) ) return 0;
463  return h;
464  }
465 
477  IHistogram1D * add(const std::string & path,
478  const IHistogram1D & hist1, const IHistogram1D & hist2) {
479  return add(path, dynamic_cast<const Histogram1D &>(hist1),
480  dynamic_cast<const Histogram1D &>(hist2));
481  }
482 
494  Histogram1D * subtract(const std::string & path,
495  const Histogram1D & h1, const Histogram1D & h2) {
496  if ( !checkBins(h1, h2) ) return 0;
497  Histogram1D * h = new Histogram1D(h1);
498  h->setTitle(path.substr(path.rfind('/') + 1));
499  for ( int i = 0; i < h->ax->bins() + 2; ++i ) {
500  h->sum[i] += h2.sum[i];
501  h->sumw[i] -= h2.sumw[i];
502  h->sumw2[i] += h2.sumw2[i];
503  }
504  if ( !tree->insert(path, h) ) return 0;
505  return h;
506  }
507 
519  IHistogram1D * subtract(const std::string & path, const IHistogram1D & hist1,
520  const IHistogram1D & hist2) {
521  return subtract(path, dynamic_cast<const Histogram1D &>(hist1),
522  dynamic_cast<const Histogram1D &>(hist2));
523  }
524 
536  Histogram1D * multiply(const std::string & path,
537  const Histogram1D & h1, const Histogram1D & h2) {
538  if ( !checkBins(h1, h2) ) return 0;
539  Histogram1D * h = new Histogram1D(h1);
540  h->setTitle(path.substr(path.rfind('/') + 1));
541  for ( int i = 0; i < h->ax->bins() + 2; ++i ) {
542  h->sumw[i] *= h2.sumw[i];
543  h->sumw2[i] += h1.sumw[i]*h1.sumw[i]*h2.sumw2[i] +
544  h2.sumw[i]*h2.sumw[i]*h1.sumw2[i];
545  }
546  if ( !tree->insert(path, h) ) return 0;
547  return h;
548  }
549 
561  IHistogram1D * multiply(const std::string & path, const IHistogram1D & hist1,
562  const IHistogram1D & hist2) {
563  return multiply(path, dynamic_cast<const Histogram1D &>(hist1),
564  dynamic_cast<const Histogram1D &>(hist2));
565  }
566 
578  Histogram1D * divide(const std::string & path,
579  const Histogram1D & h1, const Histogram1D & h2) {
580  if ( !checkBins(h1, h2) ) return 0;
581  Histogram1D * h = new Histogram1D(h1);
582  h->setTitle(path.substr(path.rfind('/') + 1));
583  for ( int i = 0; i < h->ax->bins() + 2; ++i ) {
584  if ( h2.sum[i] == 0 || h2.sumw[i] == 0.0 ) {
585  h->sum[i] = 0;
586  h->sumw[i] = h->sumw2[i] = 0.0;
587  continue;
588  }
589  h->sumw[i] /= h2.sumw[i];
590  h->sumw2[i] = h1.sumw2[i]/(h2.sumw[i]*h2.sumw[i]) +
591  h1.sumw[i]*h1.sumw[i]*h2.sumw2[i]/
592  (h2.sumw[i]*h2.sumw[i]*h2.sumw[i]*h2.sumw[i]);
593  }
594  if ( !tree->insert(path, h) ) return 0;
595  return h;
596  }
597 
609  IHistogram1D * divide(const std::string & path, const IHistogram1D & hist1,
610  const IHistogram1D & hist2) {
611  return divide(path, dynamic_cast<const Histogram1D &>(hist1),
612  dynamic_cast<const Histogram1D &>(hist2));
613  }
614 
615  inline bool _neq(double a, double b, double eps = 1e-5) const {
616  if ( a == 0 && b == 0 ) return false;
617  if ( abs(a-b) < eps*(abs(a) + abs(b)) ) return false;
618  return true;
619  }
620 
624  bool checkBins(const Histogram1D & h1, const Histogram1D & h2) const {
625  if ( _neq(h1.ax->upperEdge(), h2.ax->upperEdge()) ||
626  _neq(h1.ax->lowerEdge(), h2.ax->lowerEdge()) ||
627  _neq(h1.ax->bins(), h2.ax->bins()) ) return false;
628  if ( h1.fax && h2.fax ) return true;
629  for ( int i = 0; i < h1.ax->bins(); ++i ) {
630  if ( _neq(h1.ax->binUpperEdge(i), h2.ax->binUpperEdge(i)) ||
631  _neq(h1.ax->binLowerEdge(i), h2.ax->binLowerEdge(i)) ) return false;
632  }
633  return true;
634  }
635 
639  bool checkBins(const Histogram2D & h1, const Histogram2D & h2) const {
640  if (_neq( h1.xax->upperEdge(), h2.xax->upperEdge()) ||
641  _neq( h1.xax->lowerEdge(), h2.xax->lowerEdge()) ||
642  h1.xax->bins() != h2.xax->bins() ) return false;
643  if (_neq( h1.yax->upperEdge(), h2.yax->upperEdge()) ||
644  _neq( h1.yax->lowerEdge(), h2.yax->lowerEdge()) ||
645  h1.yax->bins() != h2.yax->bins() ) return false;
646  if ( h1.xfax && h2.xfax && h1.yfax && h2.yfax ) return true;
647  for ( int i = 0; i < h1.xax->bins(); ++i ) {
648  if ( _neq(h1.xax->binUpperEdge(i), h2.xax->binUpperEdge(i)) ||
649  _neq(h1.xax->binLowerEdge(i), h2.xax->binLowerEdge(i)) )
650  return false;
651  }
652  for ( int i = 0; i < h1.yax->bins(); ++i ) {
653  if ( _neq(h1.yax->binUpperEdge(i), h2.yax->binUpperEdge(i)) ||
654  _neq(h1.yax->binLowerEdge(i), h2.yax->binLowerEdge(i)) )
655  return false;
656  }
657  return true;
658  }
659 
663  IHistogram2D * add(const std::string & path,
664  const IHistogram2D & hist1, const IHistogram2D & hist2) {
665  return add(path, dynamic_cast<const Histogram2D &>(hist1),
666  dynamic_cast<const Histogram2D &>(hist2));
667  }
668 
672  Histogram2D * add(const std::string & path,
673  const Histogram2D & h1, const Histogram2D & h2) {
674  if ( !checkBins(h1, h2) ) return 0;
675  Histogram2D * h = new Histogram2D(h1);
676  h->setTitle(path.substr(path.rfind('/') + 1));
677  h->add(h2);
678  if ( !tree->insert(path, h) ) {
679  delete h;
680  return 0;
681  }
682  return h;
683  }
684 
688  Histogram2D * subtract(const std::string & path,
689  const Histogram2D & h1, const Histogram2D & h2) {
690  if ( !checkBins(h1, h2) ) {
691  //std::cout << "!!!!!!!" << std::endl;
692  return 0;
693  }
694  Histogram2D * h = new Histogram2D(h1);
695  h->setTitle(path.substr(path.rfind('/') + 1));
696  for ( int ix = 0; ix < h->xax->bins() + 2; ++ix )
697  for ( int iy = 0; iy < h->yax->bins() + 2; ++iy ) {
698  h->sum[ix][iy] += h2.sum[ix][iy];
699  h->sumw[ix][iy] -= h2.sumw[ix][iy];
700  h->sumw2[ix][iy] += h2.sumw2[ix][iy];
701  h->sumxw[ix][iy] -= h2.sumxw[ix][iy];
702  h->sumx2w[ix][iy] -= h2.sumx2w[ix][iy];
703  h->sumyw[ix][iy] -= h2.sumyw[ix][iy];
704  h->sumy2w[ix][iy] -= h2.sumy2w[ix][iy];
705  }
706  if ( !tree->insert(path, h) ) {
707  //std::cout << "&&&&&&&" << std::endl;
708  delete h;
709  return 0;
710  }
711  return h;
712  }
713 
717  IHistogram2D * subtract(const std::string & path,
718  const IHistogram2D & h1, const IHistogram2D & h2) {
719  return subtract(path, dynamic_cast<const Histogram2D &>(h1),
720  dynamic_cast<const Histogram2D &>(h2));
721  }
722 
726  IHistogram2D * multiply(const std::string & path,
727  const IHistogram2D & h1, const IHistogram2D & h2) {
728  return multiply(path, dynamic_cast<const Histogram2D &>(h1),
729  dynamic_cast<const Histogram2D &>(h2));
730  }
731 
735  Histogram2D * multiply(const std::string & path,
736  const Histogram2D & h1, const Histogram2D & h2) {
737  if ( !checkBins(h1, h2) ) return 0;
738  Histogram2D * h = new Histogram2D(h1);
739  h->setTitle(path.substr(path.rfind('/') + 1));
740  for ( int ix = 0; ix < h->xax->bins() + 2; ++ix )
741  for ( int iy = 0; iy < h->yax->bins() + 2; ++iy ) {
742  h->sum[ix][iy] *= h2.sum[ix][iy];
743  h->sumw[ix][iy] *= h2.sumw[ix][iy];
744  h->sumw2[ix][iy] += h1.sumw[ix][iy]*h1.sumw[ix][iy]*h2.sumw2[ix][iy] +
745  h2.sumw[ix][iy]*h2.sumw[ix][iy]*h1.sumw2[ix][iy];
746  }
747  if ( !tree->insert(path, h) ) {
748  delete h;
749  return 0;
750  }
751  return h;
752  }
753 
757  Histogram2D * divide(const std::string & path,
758  const Histogram2D & h1, const Histogram2D & h2) {
759  if ( !checkBins(h1,h2) ) return 0;
760  Histogram2D * h = new Histogram2D(h1);
761  h->setTitle(path.substr(path.rfind('/') + 1));
762  for ( int ix = 0; ix < h->xax->bins() + 2; ++ix )
763  for ( int iy = 0; iy < h->yax->bins() + 2; ++iy ) {
764  if ( h2.sum[ix][iy] == 0 || h2.sumw[ix][iy] == 0.0 ) {
765  h->sum[ix][iy] = 0;
766  h->sumw[ix][iy] = h->sumw2[ix][iy] = 0.0;
767  continue;
768  }
769  h->sumw[ix][iy] /= h2.sumw[ix][iy];
770  h->sumw2[ix][iy] = h1.sumw2[ix][iy]/(h2.sumw[ix][iy]*h2.sumw[ix][iy]) +
771  h1.sumw[ix][iy]*h1.sumw[ix][iy]*h2.sumw2[ix][iy]/
772  (h2.sumw[ix][iy]*h2.sumw[ix][iy]*h2.sumw[ix][iy]*h2.sumw[ix][iy]);
773  }
774  if ( !tree->insert(path, h) ) {
775  delete h;
776  return 0;
777  }
778  return h;
779  }
780 
781 
785  IHistogram2D * divide(const std::string & path,
786  const IHistogram2D & h1, const IHistogram2D & h2) {
787  return divide(path, dynamic_cast<const Histogram2D &>(h1),
788  dynamic_cast<const Histogram2D &>(h2));
789  }
790 
794  IHistogram3D * add(const std::string &,
795  const IHistogram3D &, const IHistogram3D &) {
796  return error<IHistogram3D>("3D histograms");
797  }
798 
802  IHistogram3D * subtract(const std::string &,
803  const IHistogram3D &, const IHistogram3D &) {
804  return error<IHistogram3D>("3D histograms");
805  }
806 
810  IHistogram3D * multiply(const std::string &,
811  const IHistogram3D &, const IHistogram3D &) {
812  return error<IHistogram3D>("3D histograms");
813  }
814 
818  IHistogram3D * divide(const std::string &,
819  const IHistogram3D &, const IHistogram3D &) {
820  return error<IHistogram3D>("3D histograms");
821  }
822 
827  IHistogram1D * projectionX(const std::string & path, const IHistogram2D & h) {
828  return projectionX(path, dynamic_cast<const Histogram2D &>(h));
829  }
830 
835  Histogram1D * projectionX(const std::string & path, const Histogram2D & h) {
836  return sliceX(path, h, 0, h.yax->bins() - 1);
837  }
838 
843  IHistogram1D * projectionY(const std::string & path, const IHistogram2D & h) {
844  return projectionY(path, dynamic_cast<const Histogram2D &>(h));
845  }
846 
851  Histogram1D * projectionY(const std::string & path, const Histogram2D & h) {
852  return sliceY(path, h, 0, h.xax->bins() - 1);
853  }
854 
859  IHistogram1D *
860  sliceX(const std::string & path, const IHistogram2D & h, int i) {
861  return sliceX(path, dynamic_cast<const Histogram2D &>(h), i, i);
862  }
863 
868  Histogram1D *
869  sliceX(const std::string & path, const Histogram2D & h, int i) {
870  return sliceX(path, h, i, i);
871  }
872 
877  IHistogram1D *
878  sliceY(const std::string & path, const IHistogram2D & h, int i) {
879  return sliceY(path, dynamic_cast<const Histogram2D &>(h), i, i);
880  }
881 
886  Histogram1D * sliceY(const std::string & path, const Histogram2D & h, int i) {
887  return sliceY(path, h, i, i);
888  }
889 
894  IHistogram1D *
895  sliceX(const std::string & path, const IHistogram2D & h, int il, int iu) {
896  return sliceX(path, dynamic_cast<const Histogram2D &>(h), il, iu);
897  }
898 
903  Histogram1D *
904  sliceX(const std::string & path, const Histogram2D & h2, int il, int iu) {
905  Histogram1D * h1;
906  if ( h2.xfax )
907  h1 = new Histogram1D(h2.xfax->bins(), h2.xfax->lowerEdge(),
908  h2.xfax->upperEdge());
909  else {
910  std::vector<double> edges(h2.xax->bins() + 1);
911  edges.push_back(h2.xax->lowerEdge());
912  for ( int i = 0; i < h2.xax->bins(); ++i )
913  edges.push_back(h2.xax->binLowerEdge(i));
914  h1 = new Histogram1D(edges);
915  }
916  for ( int ix = 0; ix < h2.xax->bins() + 2; ++ix )
917  for ( int iy = il + 2; iy <= iu + 2; ++iy ) {
918  h1->sum[ix] += h2.sum[ix][iy];
919  h1->sumw[ix] += h2.sumw[ix][iy];
920  h1->sumw2[ix] += h2.sumw2[ix][iy];
921  h1->sumxw[ix] += h2.sumxw[ix][iy];
922  h1->sumx2w[ix] += h2.sumx2w[ix][iy];
923  }
924  if ( !tree->insert(path, h1) ) {
925  delete h1;
926  return 0;
927  }
928  return h1;
929  }
930 
935  IHistogram1D *
936  sliceY(const std::string & path, const IHistogram2D & h, int il, int iu) {
937  return sliceY(path, dynamic_cast<const Histogram2D &>(h), il, iu);
938  }
939 
940  Histogram1D *
941  sliceY(const std::string & path, const Histogram2D & h2, int il, int iu) {
942  Histogram1D * h1;
943  if ( h2.yfax )
944  h1 = new Histogram1D(h2.yfax->bins(), h2.yfax->lowerEdge(),
945  h2.yfax->upperEdge());
946  else {
947  std::vector<double> edges(h2.yax->bins() + 1);
948  edges.push_back(h2.yax->lowerEdge());
949  for ( int i = 0; i < h2.yax->bins(); ++i )
950  edges.push_back(h2.yax->binLowerEdge(i));
951  h1 = new Histogram1D(edges);
952  }
953  for ( int iy = 0; iy < h2.yax->bins() + 2; ++iy )
954  for ( int ix = il + 2; ix <= iu + 2; ++ix ) {
955  h1->sum[iy] += h2.sum[ix][iy];
956  h1->sumw[iy] += h2.sumw[ix][iy];
957  h1->sumw2[iy] += h2.sumw2[ix][iy];
958  h1->sumxw[iy] += h2.sumyw[ix][iy];
959  h1->sumx2w[iy] += h2.sumy2w[ix][iy];
960  }
961  if ( !tree->insert(path, h1) ) {
962  delete h1;
963  return 0;
964  }
965  return h1;
966  }
967 
972  IHistogram2D * projectionXY(const std::string &, const IHistogram3D &) {
973  return error<IHistogram2D>("2D histograms");
974  }
975 
980  IHistogram2D * projectionXZ(const std::string &, const IHistogram3D &) {
981  return error<IHistogram2D>("2D histograms");
982  }
983 
988  IHistogram2D * projectionYZ(const std::string &, const IHistogram3D &) {
989  return error<IHistogram2D>("2D histograms");
990  }
991 
997  IHistogram2D * sliceXY(const std::string &, const IHistogram3D &, int, int) {
998  return error<IHistogram2D>("2D histograms");
999  }
1000 
1006  IHistogram2D * sliceXZ(const std::string &, const IHistogram3D &, int, int) {
1007  return error<IHistogram2D>("2D histograms");
1008  }
1009 
1015  IHistogram2D * sliceYZ(const std::string &, const IHistogram3D &, int, int) {
1016  return error<IHistogram2D>("2D histograms");
1017  }
1018 
1019 
1020 private:
1021 
1023  template <typename T>
1024  static T * error(std::string feature) {
1025  throw std::runtime_error("LWH cannot handle " + feature + ".");
1026  }
1027 
1030 
1031 };
1032 
1033 }
1034 
1035 #endif /* LWH_HistogramFactory_H */
std::vector< double > sumx2w
The weighted x-square-values.
Definition: Histogram1D.h:514
ICloud2D * createCloud2D(const std::string &)
LWH cannot create a ICloud2D, an unbinned 2-dimensional histogram.
ICloud3D * createCopy(const std::string &, const ICloud3D &)
LWH cannot create a copy of an ICloud3D.
IHistogram1D * add(const std::string &path, const IHistogram1D &hist1, const IHistogram1D &hist2)
Create an IHistogram1D by adding two IHistogram1D.
IHistogram3D * createHistogram3D(const std::string &, const std::string &, int, double, double, int, double, double, int, double, double, const std::string &="")
LWH cannot create a IHistogram3D.
IHistogram2D * createHistogram2D(const std::string &pathAndTitle, int nx, double xlo, double xup, int ny, double ylo, double yup)
Create a IHistogram2D.
virtual ~HistogramFactory()
Destructor.
Histogram1D * sliceY(const std::string &path, const Histogram2D &h, int i)
LWH cannot create an IHistogram1D by slicing an IHistogram2D parallel to the x axis at a given bin...
IProfile2D * createProfile2D(const std::string &, int, double, double, int, double, double, double, double)
LWH cannot create a IProfile2D.
Histogram1D * projectionX(const std::string &path, const Histogram2D &h)
LWH cannot create an IHistogram1D by projecting an IHistogram2D along its x axis. ...
bool setTitle(const std::string &title)
Set the histogram title.
Definition: Histogram2D.h:113
IHistogram1D * subtract(const std::string &path, const IHistogram1D &hist1, const IHistogram1D &hist2)
Create an IHistogram1D by subtracting two IHistogram1D.
bool destroy(IBaseHistogram *hist)
Destroy an IBaseHistogram object.
IHistogram3D * subtract(const std::string &, const IHistogram3D &, const IHistogram3D &)
LWH cannot create an IHistogram3D by subtracting two IHistogram3D.
double upperEdge() const
Get the upper edge of the IAxis.
Definition: Axis.h:69
IHistogram3D * createHistogram3D(const std::string &, const std::string &, const std::vector< double > &, const std::vector< double > &, const std::vector< double > &, const std::string &="")
LWH cannot create a IHistogram3D.
Histogram2D * add(const std::string &path, const Histogram2D &h1, const Histogram2D &h2)
LWH cannot create an IHistogram2D by adding two IHistogram2D.
IHistogram2D * sliceXZ(const std::string &, const IHistogram3D &, int, int)
LWH cannot create an IHistogram2D by slicing an IHistogram3D perpendicular to the Y axis...
IProfile1D * createProfile1D(const std::string &, const std::string &, int, double, double, const std::string &="")
LWH cannot create a IProfile1D.
ICloud1D * createCloud1D(const std::string &)
LWH cannot create a ICloud1D, an unbinned 1-dimensional histogram.
IHistogram3D * multiply(const std::string &, const IHistogram3D &, const IHistogram3D &)
LWH cannot create an IHistogram3D by multiplying two IHistogram3D.
std::vector< std::vector< double > > sumyw
The weighted y-values.
Definition: Histogram2D.h:856
IAxis * ax
The axis.
Definition: Histogram1D.h:493
HistogramFactory(Tree &t)
Standard constructor.
User level interface to 1D Histogram.
Definition: Histogram2D.h:25
IProfile2D * createCopy(const std::string &, const IProfile2D &)
LWH cannot create a copy of an IProfile2D.
ICloud3D * createCloud3D(const std::string &)
LWH cannot create a ICloud3D, an unbinned 3-dimensional histogram.
Tree * tree
The tree where the actual histograms are stored.
IProfile1D * createProfile1D(const std::string &, const std::string &, int, double, double, double, double, const std::string &="")
LWH cannot create a IProfile1D.
Histogram2D * multiply(const std::string &path, const Histogram2D &h1, const Histogram2D &h2)
LWH cannot create an IHistogram2D by multiplying two IHistogram2D.
IHistogram2D * createCopy(const std::string &path, const IHistogram2D &hist)
Create a copy of an IHistogram2D.
ICloud2D * createCloud2D(const std::string &, const std::string &, int=-1, const std::string &="")
LWH cannot create a ICloud2D, an unbinned 2-dimensional histogram.
Histogram1D * sliceX(const std::string &path, const Histogram2D &h, int i)
LWH cannot create an IHistogram1D by slicing an IHistogram2D parallel to the y axis at a given bin...
std::vector< double > sumw
The weights.
Definition: Histogram1D.h:505
IHistogram1D * projectionY(const std::string &path, const IHistogram2D &h)
LWH cannot create an IHistogram1D by projecting an IHistogram2D along its y axis. ...
bool add(const Histogram2D &h)
Add to this Histogram2D the contents of another IHistogram2D.
Definition: Histogram2D.h:577
Histogram1D * divide(const std::string &path, const Histogram1D &h1, const Histogram1D &h2)
Create n Histogram1D by dividing two Histogram1D.
IHistogram2D * subtract(const std::string &path, const IHistogram2D &h1, const IHistogram2D &h2)
LWH cannot create an IHistogram2D by subtracting two IHistogram2D.
IHistogram3D * add(const std::string &, const IHistogram3D &, const IHistogram3D &)
LWH cannot create an IHistogram3D by adding two IHistogram3D.
IHistogram2D * divide(const std::string &path, const IHistogram2D &h1, const IHistogram2D &h2)
LWH cannot create an IHistogram2D by dividing two IHistogram2D.
IHistogram1D * createCopy(const std::string &path, const IHistogram1D &hist)
Create a copy of an IHistogram1D.
IHistogram1D * sliceY(const std::string &path, const IHistogram2D &h, int i)
LWH cannot create an IHistogram1D by slicing an IHistogram2D parallel to the x axis at a given bin...
IProfile1D * createCopy(const std::string &, const IProfile1D &)
LWH cannot create a copy of an IProfile1D.
IProfile1D * createProfile1D(const std::string &, int, double, double)
LWH cannot create a IProfile1D.
IHistogram2D * add(const std::string &path, const IHistogram2D &hist1, const IHistogram2D &hist2)
LWH cannot create an IHistogram2D by adding two IHistogram2D.
IHistogram1D * sliceY(const std::string &path, const IHistogram2D &h, int il, int iu)
LWH cannot create an IHistogram1D by slicing an IHistogram2D parallel to the x axis between two bins ...
IHistogram2D * projectionYZ(const std::string &, const IHistogram3D &)
LWH cannot create an IHistogram2D by projecting an IHistogram3D on the y-z plane. ...
IHistogram2D * multiply(const std::string &path, const IHistogram2D &h1, const IHistogram2D &h2)
LWH cannot create an IHistogram2D by multiplying two IHistogram2D.
IProfile2D * createProfile2D(const std::string &, const std::string &, const std::vector< double > &, const std::vector< double > &, double, double, const std::string &="")
LWH cannot create a IProfile2D.
std::vector< int > sum
The counts.
Definition: Histogram1D.h:502
IHistogram2D * projectionXZ(const std::string &, const IHistogram3D &)
LWH cannot create an IHistogram2D by projecting an IHistogram3D on the x-z plane. ...
IHistogram1D * projectionX(const std::string &path, const IHistogram2D &h)
LWH cannot create an IHistogram1D by projecting an IHistogram2D along its x axis. ...
ICloud1D * createCopy(const std::string &, const ICloud1D &)
LWH cannot create a copy of an ICloud1D.
Histogram1D * subtract(const std::string &path, const Histogram1D &h1, const Histogram1D &h2)
Create a Histogram1D by subtracting two Histogram1D.
IHistogram2D * createHistogram2D(const std::string &path, const std::string &title, const std::vector< double > &xedges, const std::vector< double > &yedges, const std::string &="")
Create a IHistogram2D.
IProfile1D * createProfile1D(const std::string &, const std::string &, const std::vector< double > &, double, double, const std::string &="")
LWH cannot create a IProfile1D.
bool checkBins(const Histogram1D &h1, const Histogram1D &h2) const
Check if two histograms have the same bins.
bool checkBins(const Histogram2D &h1, const Histogram2D &h2) const
Check if two histograms have the same bins.
std::vector< std::vector< double > > sumw2
The squared weights.
Definition: Histogram2D.h:847
IHistogram1D * multiply(const std::string &path, const IHistogram1D &hist1, const IHistogram1D &hist2)
Create an IHistogram1D by multiplying two IHistogram1D.
std::vector< std::vector< double > > sumw
The weights.
Definition: Histogram2D.h:844
IHistogram3D * createHistogram3D(const std::string &, int, double, double, int, double, double, int, double, double)
LWH cannot create a IHistogram3D.
IHistogram1D * createHistogram1D(const std::string &path, const std::string &title, const std::vector< double > &binEdges, const std::string &="")
Create a IHistogram1D.
int bins() const
The number of bins (excluding underflow and overflow) on the IAxis.
Definition: Axis.h:76
IAxis * xax
The axis.
Definition: Histogram2D.h:823
IHistogram1D * divide(const std::string &path, const IHistogram1D &hist1, const IHistogram1D &hist2)
Create an IHistogram1D by dividing two IHistogram1D.
ICloud3D * createCloud3D(const std::string &, const std::string &, int=-1, const std::string &="")
LWH cannot create a ICloud3D, an unbinned 3-dimensional histogram.
Histogram2D * subtract(const std::string &path, const Histogram2D &h1, const Histogram2D &h2)
LWH cannot create an IHistogram2D by subtracting two IHistogram2D.
ICloud1D * createCloud1D(const std::string &, const std::string &, int=-1, const std::string &="")
LWH cannot create a ICloud1D, an unbinned 1-dimensional histogram.
double lowerEdge() const
Get the lower edge of the IAxis.
Definition: Axis.h:62
IHistogram1D * sliceX(const std::string &path, const IHistogram2D &h, int i)
LWH cannot create an IHistogram1D by slicing an IHistogram2D parallel to the y axis at a given bin...
Axis * yfax
Pointer (possibly null) to a axis with fixed bin width.
Definition: Histogram2D.h:835
Histogram1D * projectionY(const std::string &path, const Histogram2D &h)
LWH cannot create an IHistogram1D by projecting an IHistogram2D along its y axis. ...
User level interface for factory classes of Histograms (binned, unbinned, and profile).
User level interface to 1D Histogram.
Definition: Histogram1D.h:29
IProfile2D * createProfile2D(const std::string &, const std::string &, int, double, double, int, double, double, const std::string &="")
LWH cannot create a IProfile2D.
IHistogram2D * projectionXY(const std::string &, const IHistogram3D &)
LWH cannot create an IHistogram2D by projecting an IHistogram3D on the x-y plane. ...
IProfile1D * createProfile1D(const std::string &, const std::string &, const std::vector< double > &, const std::string &="")
LWH cannot create a IProfile1D.
IProfile2D * createProfile2D(const std::string &, int, double, double, int, double, double)
LWH cannot create a IProfile2D.
static T * error(std::string feature)
Throw a suitable error.
bool setTitle(const std::string &title)
Set the histogram title.
Definition: Histogram1D.h:95
std::vector< std::vector< double > > sumx2w
The weighted x-square-values.
Definition: Histogram2D.h:853
Histogram2D * divide(const std::string &path, const Histogram2D &h1, const Histogram2D &h2)
LWH cannot create an IHistogram2D by dividing two IHistogram2D.
IAxis * yax
The axis.
Definition: Histogram2D.h:832
IProfile2D * createProfile2D(const std::string &, const std::string &, int, double, double, int, double, double, double, double, const std::string &="")
LWH cannot create a IProfile2D.
IHistogram2D * sliceXY(const std::string &, const IHistogram3D &, int, int)
LWH cannot create an IHistogram2D by slicing an IHistogram3D perpendicular to the Z axis...
std::vector< std::vector< int > > sum
The counts.
Definition: Histogram2D.h:841
IProfile2D * createProfile2D(const std::string &, const std::string &, const std::vector< double > &, const std::vector< double > &, const std::string &="")
LWH cannot create a IProfile2D.
std::vector< std::vector< double > > sumy2w
The weighted y-square-values.
Definition: Histogram2D.h:859
Axis * fax
Pointer (possibly null) to a axis with fixed bin width.
Definition: Histogram1D.h:496
The LWH namespace contains a Light-Weight Histogram package which implements the most rudimentary his...
IHistogram1D * createHistogram1D(const std::string &pathAndTitle, int nBins, double lowerEdge, double upperEdge)
Create a IHistogram1D.
Histogram1D * add(const std::string &path, const Histogram1D &hist1, const Histogram1D &hist2)
Create a Histogram1D by adding two Histogram1D.
The Tree class is a simple implementation of the AIDA::ITree interface.
Definition: Tree.h:32
IProfile1D * createProfile1D(const std::string &, int, double, double, double, double)
LWH cannot create a IProfile1D.
IHistogram3D * divide(const std::string &, const IHistogram3D &, const IHistogram3D &)
LWH cannot create an IHistogram3D by dividing two IHistogram3D.
bool add(const Histogram1D &h)
Add to this Histogram1D the contents of another IHistogram1D.
Definition: Histogram1D.h:355
IHistogram2D * createHistogram2D(const std::string &path, const std::string &title, int nx, double xlo, double xup, int ny, double ylo, double yup, const std::string &="")
Create a IHistogram2D.
std::vector< double > sumxw
The weighted x-values.
Definition: Histogram1D.h:511
ICloud2D * createCopy(const std::string &, const ICloud2D &)
LWH cannot create a copy of an ICloud2D.
IHistogram2D * sliceYZ(const std::string &, const IHistogram3D &, int, int)
LWH cannot create an IHistogram2D by slicing an IHistogram3D perpendicular to the X axis...
Histogram1D * sliceX(const std::string &path, const Histogram2D &h2, int il, int iu)
LWH cannot create an IHistogram1D by slicing an IHistogram2D parallel to the y axis between two bins ...
IHistogram1D * sliceX(const std::string &path, const IHistogram2D &h, int il, int iu)
LWH cannot create an IHistogram1D by slicing an IHistogram2D parallel to the y axis between two bins ...
Histogram1D * multiply(const std::string &path, const Histogram1D &h1, const Histogram1D &h2)
Create a Histogram1D by multiplying two Histogram1D.
IHistogram3D * createCopy(const std::string &, const IHistogram3D &)
LWH cannot create a copy of an IHistogram3D.
IHistogram1D * createHistogram1D(const std::string &path, const std::string &title, int nBins, double lowerEdge, double upperEdge, const std::string &="")
Create a IHistogram1D.
std::vector< double > sumw2
The squared weights.
Definition: Histogram1D.h:508
std::vector< std::vector< double > > sumxw
The weighted x-values.
Definition: Histogram2D.h:850
Axis * xfax
Pointer (possibly null) to a axis with fixed bin width.
Definition: Histogram2D.h:826