thepeg is hosted by Hepforge, IPPP Durham
ThePEG  2.2.1
HistogramFactory.h
1 // -*- C++ -*-
2 //
3 // HistogramFactory.h is a part of ThePEG - Toolkit for HEP Event Generation
4 // Copyright (C) 1999-2019 Leif Lonnblad
5 //
6 // ThePEG is licenced under version 3 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  using std::abs;
617  if ( a == 0 && b == 0 ) return false;
618  if ( abs(a-b) < eps*(abs(a) + abs(b)) ) return false;
619  return true;
620  }
621 
625  bool checkBins(const Histogram1D & h1, const Histogram1D & h2) const {
626  if ( _neq(h1.ax->upperEdge(), h2.ax->upperEdge()) ||
627  _neq(h1.ax->lowerEdge(), h2.ax->lowerEdge()) ||
628  _neq(h1.ax->bins(), h2.ax->bins()) ) return false;
629  if ( h1.fax && h2.fax ) return true;
630  for ( int i = 0; i < h1.ax->bins(); ++i ) {
631  if ( _neq(h1.ax->binUpperEdge(i), h2.ax->binUpperEdge(i)) ||
632  _neq(h1.ax->binLowerEdge(i), h2.ax->binLowerEdge(i)) ) return false;
633  }
634  return true;
635  }
636 
640  bool checkBins(const Histogram2D & h1, const Histogram2D & h2) const {
641  if (_neq( h1.xax->upperEdge(), h2.xax->upperEdge()) ||
642  _neq( h1.xax->lowerEdge(), h2.xax->lowerEdge()) ||
643  h1.xax->bins() != h2.xax->bins() ) return false;
644  if (_neq( h1.yax->upperEdge(), h2.yax->upperEdge()) ||
645  _neq( h1.yax->lowerEdge(), h2.yax->lowerEdge()) ||
646  h1.yax->bins() != h2.yax->bins() ) return false;
647  if ( h1.xfax && h2.xfax && h1.yfax && h2.yfax ) return true;
648  for ( int i = 0; i < h1.xax->bins(); ++i ) {
649  if ( _neq(h1.xax->binUpperEdge(i), h2.xax->binUpperEdge(i)) ||
650  _neq(h1.xax->binLowerEdge(i), h2.xax->binLowerEdge(i)) )
651  return false;
652  }
653  for ( int i = 0; i < h1.yax->bins(); ++i ) {
654  if ( _neq(h1.yax->binUpperEdge(i), h2.yax->binUpperEdge(i)) ||
655  _neq(h1.yax->binLowerEdge(i), h2.yax->binLowerEdge(i)) )
656  return false;
657  }
658  return true;
659  }
660 
664  IHistogram2D * add(const std::string & path,
665  const IHistogram2D & hist1, const IHistogram2D & hist2) {
666  return add(path, dynamic_cast<const Histogram2D &>(hist1),
667  dynamic_cast<const Histogram2D &>(hist2));
668  }
669 
673  Histogram2D * add(const std::string & path,
674  const Histogram2D & h1, const Histogram2D & h2) {
675  if ( !checkBins(h1, h2) ) return 0;
676  Histogram2D * h = new Histogram2D(h1);
677  h->setTitle(path.substr(path.rfind('/') + 1));
678  h->add(h2);
679  if ( !tree->insert(path, h) ) {
680  delete h;
681  return 0;
682  }
683  return h;
684  }
685 
689  Histogram2D * subtract(const std::string & path,
690  const Histogram2D & h1, const Histogram2D & h2) {
691  if ( !checkBins(h1, h2) ) {
692  //std::cout << "!!!!!!!" << std::endl;
693  return 0;
694  }
695  Histogram2D * h = new Histogram2D(h1);
696  h->setTitle(path.substr(path.rfind('/') + 1));
697  for ( int ix = 0; ix < h->xax->bins() + 2; ++ix )
698  for ( int iy = 0; iy < h->yax->bins() + 2; ++iy ) {
699  h->sum[ix][iy] += h2.sum[ix][iy];
700  h->sumw[ix][iy] -= h2.sumw[ix][iy];
701  h->sumw2[ix][iy] += h2.sumw2[ix][iy];
702  h->sumxw[ix][iy] -= h2.sumxw[ix][iy];
703  h->sumx2w[ix][iy] -= h2.sumx2w[ix][iy];
704  h->sumyw[ix][iy] -= h2.sumyw[ix][iy];
705  h->sumy2w[ix][iy] -= h2.sumy2w[ix][iy];
706  }
707  if ( !tree->insert(path, h) ) {
708  //std::cout << "&&&&&&&" << std::endl;
709  delete h;
710  return 0;
711  }
712  return h;
713  }
714 
718  IHistogram2D * subtract(const std::string & path,
719  const IHistogram2D & h1, const IHistogram2D & h2) {
720  return subtract(path, dynamic_cast<const Histogram2D &>(h1),
721  dynamic_cast<const Histogram2D &>(h2));
722  }
723 
727  IHistogram2D * multiply(const std::string & path,
728  const IHistogram2D & h1, const IHistogram2D & h2) {
729  return multiply(path, dynamic_cast<const Histogram2D &>(h1),
730  dynamic_cast<const Histogram2D &>(h2));
731  }
732 
736  Histogram2D * multiply(const std::string & path,
737  const Histogram2D & h1, const Histogram2D & h2) {
738  if ( !checkBins(h1, h2) ) return 0;
739  Histogram2D * h = new Histogram2D(h1);
740  h->setTitle(path.substr(path.rfind('/') + 1));
741  for ( int ix = 0; ix < h->xax->bins() + 2; ++ix )
742  for ( int iy = 0; iy < h->yax->bins() + 2; ++iy ) {
743  h->sum[ix][iy] *= h2.sum[ix][iy];
744  h->sumw[ix][iy] *= h2.sumw[ix][iy];
745  h->sumw2[ix][iy] += h1.sumw[ix][iy]*h1.sumw[ix][iy]*h2.sumw2[ix][iy] +
746  h2.sumw[ix][iy]*h2.sumw[ix][iy]*h1.sumw2[ix][iy];
747  }
748  if ( !tree->insert(path, h) ) {
749  delete h;
750  return 0;
751  }
752  return h;
753  }
754 
758  Histogram2D * divide(const std::string & path,
759  const Histogram2D & h1, const Histogram2D & h2) {
760  if ( !checkBins(h1,h2) ) return 0;
761  Histogram2D * h = new Histogram2D(h1);
762  h->setTitle(path.substr(path.rfind('/') + 1));
763  for ( int ix = 0; ix < h->xax->bins() + 2; ++ix )
764  for ( int iy = 0; iy < h->yax->bins() + 2; ++iy ) {
765  if ( h2.sum[ix][iy] == 0 || h2.sumw[ix][iy] == 0.0 ) {
766  h->sum[ix][iy] = 0;
767  h->sumw[ix][iy] = h->sumw2[ix][iy] = 0.0;
768  continue;
769  }
770  h->sumw[ix][iy] /= h2.sumw[ix][iy];
771  h->sumw2[ix][iy] = h1.sumw2[ix][iy]/(h2.sumw[ix][iy]*h2.sumw[ix][iy]) +
772  h1.sumw[ix][iy]*h1.sumw[ix][iy]*h2.sumw2[ix][iy]/
773  (h2.sumw[ix][iy]*h2.sumw[ix][iy]*h2.sumw[ix][iy]*h2.sumw[ix][iy]);
774  }
775  if ( !tree->insert(path, h) ) {
776  delete h;
777  return 0;
778  }
779  return h;
780  }
781 
782 
786  IHistogram2D * divide(const std::string & path,
787  const IHistogram2D & h1, const IHistogram2D & h2) {
788  return divide(path, dynamic_cast<const Histogram2D &>(h1),
789  dynamic_cast<const Histogram2D &>(h2));
790  }
791 
795  IHistogram3D * add(const std::string &,
796  const IHistogram3D &, const IHistogram3D &) {
797  return error<IHistogram3D>("3D histograms");
798  }
799 
803  IHistogram3D * subtract(const std::string &,
804  const IHistogram3D &, const IHistogram3D &) {
805  return error<IHistogram3D>("3D histograms");
806  }
807 
811  IHistogram3D * multiply(const std::string &,
812  const IHistogram3D &, const IHistogram3D &) {
813  return error<IHistogram3D>("3D histograms");
814  }
815 
819  IHistogram3D * divide(const std::string &,
820  const IHistogram3D &, const IHistogram3D &) {
821  return error<IHistogram3D>("3D histograms");
822  }
823 
828  IHistogram1D * projectionX(const std::string & path, const IHistogram2D & h) {
829  return projectionX(path, dynamic_cast<const Histogram2D &>(h));
830  }
831 
836  Histogram1D * projectionX(const std::string & path, const Histogram2D & h) {
837  return sliceX(path, h, 0, h.yax->bins() - 1);
838  }
839 
844  IHistogram1D * projectionY(const std::string & path, const IHistogram2D & h) {
845  return projectionY(path, dynamic_cast<const Histogram2D &>(h));
846  }
847 
852  Histogram1D * projectionY(const std::string & path, const Histogram2D & h) {
853  return sliceY(path, h, 0, h.xax->bins() - 1);
854  }
855 
860  IHistogram1D *
861  sliceX(const std::string & path, const IHistogram2D & h, int i) {
862  return sliceX(path, dynamic_cast<const Histogram2D &>(h), i, i);
863  }
864 
869  Histogram1D *
870  sliceX(const std::string & path, const Histogram2D & h, int i) {
871  return sliceX(path, h, i, i);
872  }
873 
878  IHistogram1D *
879  sliceY(const std::string & path, const IHistogram2D & h, int i) {
880  return sliceY(path, dynamic_cast<const Histogram2D &>(h), i, i);
881  }
882 
887  Histogram1D * sliceY(const std::string & path, const Histogram2D & h, int i) {
888  return sliceY(path, h, i, i);
889  }
890 
895  IHistogram1D *
896  sliceX(const std::string & path, const IHistogram2D & h, int il, int iu) {
897  return sliceX(path, dynamic_cast<const Histogram2D &>(h), il, iu);
898  }
899 
904  Histogram1D *
905  sliceX(const std::string & path, const Histogram2D & h2, int il, int iu) {
906  Histogram1D * h1;
907  if ( h2.xfax )
908  h1 = new Histogram1D(h2.xfax->bins(), h2.xfax->lowerEdge(),
909  h2.xfax->upperEdge());
910  else {
911  std::vector<double> edges(h2.xax->bins() + 1);
912  edges.push_back(h2.xax->lowerEdge());
913  for ( int i = 0; i < h2.xax->bins(); ++i )
914  edges.push_back(h2.xax->binLowerEdge(i));
915  h1 = new Histogram1D(edges);
916  }
917  for ( int ix = 0; ix < h2.xax->bins() + 2; ++ix )
918  for ( int iy = il + 2; iy <= iu + 2; ++iy ) {
919  h1->sum[ix] += h2.sum[ix][iy];
920  h1->sumw[ix] += h2.sumw[ix][iy];
921  h1->sumw2[ix] += h2.sumw2[ix][iy];
922  h1->sumxw[ix] += h2.sumxw[ix][iy];
923  h1->sumx2w[ix] += h2.sumx2w[ix][iy];
924  }
925  if ( !tree->insert(path, h1) ) {
926  delete h1;
927  return 0;
928  }
929  return h1;
930  }
931 
936  IHistogram1D *
937  sliceY(const std::string & path, const IHistogram2D & h, int il, int iu) {
938  return sliceY(path, dynamic_cast<const Histogram2D &>(h), il, iu);
939  }
940 
941  Histogram1D *
942  sliceY(const std::string & path, const Histogram2D & h2, int il, int iu) {
943  Histogram1D * h1;
944  if ( h2.yfax )
945  h1 = new Histogram1D(h2.yfax->bins(), h2.yfax->lowerEdge(),
946  h2.yfax->upperEdge());
947  else {
948  std::vector<double> edges(h2.yax->bins() + 1);
949  edges.push_back(h2.yax->lowerEdge());
950  for ( int i = 0; i < h2.yax->bins(); ++i )
951  edges.push_back(h2.yax->binLowerEdge(i));
952  h1 = new Histogram1D(edges);
953  }
954  for ( int iy = 0; iy < h2.yax->bins() + 2; ++iy )
955  for ( int ix = il + 2; ix <= iu + 2; ++ix ) {
956  h1->sum[iy] += h2.sum[ix][iy];
957  h1->sumw[iy] += h2.sumw[ix][iy];
958  h1->sumw2[iy] += h2.sumw2[ix][iy];
959  h1->sumxw[iy] += h2.sumyw[ix][iy];
960  h1->sumx2w[iy] += h2.sumy2w[ix][iy];
961  }
962  if ( !tree->insert(path, h1) ) {
963  delete h1;
964  return 0;
965  }
966  return h1;
967  }
968 
973  IHistogram2D * projectionXY(const std::string &, const IHistogram3D &) {
974  return error<IHistogram2D>("2D histograms");
975  }
976 
981  IHistogram2D * projectionXZ(const std::string &, const IHistogram3D &) {
982  return error<IHistogram2D>("2D histograms");
983  }
984 
989  IHistogram2D * projectionYZ(const std::string &, const IHistogram3D &) {
990  return error<IHistogram2D>("2D histograms");
991  }
992 
998  IHistogram2D * sliceXY(const std::string &, const IHistogram3D &, int, int) {
999  return error<IHistogram2D>("2D histograms");
1000  }
1001 
1007  IHistogram2D * sliceXZ(const std::string &, const IHistogram3D &, int, int) {
1008  return error<IHistogram2D>("2D histograms");
1009  }
1010 
1016  IHistogram2D * sliceYZ(const std::string &, const IHistogram3D &, int, int) {
1017  return error<IHistogram2D>("2D histograms");
1018  }
1019 
1020 
1021 private:
1022 
1024  template <typename T>
1025  static T * error(std::string feature) {
1026  throw std::runtime_error("LWH cannot handle " + feature + ".");
1027  }
1028 
1031 
1032 };
1033 
1034 }
1035 
1036 #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