9 #ifndef LWH_HistogramFactory_H
10 #define LWH_HistogramFactory_H
15 #include "AIHistogramFactory.h"
16 #include "Histogram1D.h"
17 #include "Histogram2D.h"
53 IManagedObject * mo =
dynamic_cast<IManagedObject *
>(hist);
54 if ( !mo )
return false;
55 return tree->rm(tree->findPath(*mo));
61 ICloud1D * createCloud1D(
const std::string &,
const std::string &,
62 int = -1,
const std::string & =
"") {
63 return error<ICloud1D>(
"ICloud1D");
70 return error<ICloud1D>(
"ICloud1D");
76 ICloud1D *
createCopy(
const std::string &,
const ICloud1D &) {
77 return error<ICloud1D>(
"ICloud1D");
83 ICloud2D * createCloud2D(
const std::string &,
const std::string &,
int = -1,
84 const std::string & =
"") {
85 return error<ICloud2D>(
"ICloud2D");
93 return error<ICloud2D>(
"ICloud2D");
99 ICloud2D *
createCopy(
const std::string &,
const ICloud2D &) {
100 return error<ICloud2D>(
"ICloud2D");
106 ICloud3D * createCloud3D(
const std::string &,
const std::string &,
int = -1,
107 const std::string & =
"") {
108 return error<ICloud3D>(
"ICloud3D");
115 return error<ICloud3D>(
"ICloud3D");
121 ICloud3D *
createCopy(
const std::string &,
const ICloud3D &) {
122 return error<ICloud3D>(
"ICloud3D");
142 createHistogram1D(
const std::string & path,
const std::string & title,
143 int nBins,
double lowerEdge,
double upperEdge,
144 const std::string & =
"") {
147 if ( !tree->insert(path, hist) ) {
150 throw std::runtime_error(
"LWH could not create histogram '"
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);
189 createHistogram1D(
const std::string & path,
const std::string & title,
190 const std::vector<double> & binEdges,
191 const std::string & =
"") {
194 if ( !tree->insert(path, hist) ) {
197 throw std::runtime_error(
"LWH could not create histogram '"
214 createCopy(
const std::string & path,
const IHistogram1D & hist) {
216 h->
setTitle(path.substr(path.rfind(
'/') + 1));
217 if ( !tree->insert(path, h) ) {
220 throw std::runtime_error(
"LWH could not create a copy of histogram '"
221 + hist.title() +
"'." );
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 & =
"") {
236 if ( !tree->insert(path, hist) ) {
239 throw std::runtime_error(
"LWH could not create histogram '"
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);
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 & =
"") {
265 if ( !tree->insert(path, hist) ) {
268 throw std::runtime_error(
"LWH could not create histogram '"
278 createCopy(
const std::string & path,
const IHistogram2D & hist) {
280 h->
setTitle(path.substr(path.rfind(
'/') + 1));
281 if ( !tree->insert(path, h) ) {
284 throw std::runtime_error(
"LWH could not create a copy of histogram '"
285 + hist.title() +
"'." );
293 IHistogram3D * createHistogram3D(
const std::string &,
const std::string &,
294 int,
double,
double,
int,
double,
double,
296 const std::string & =
"") {
297 return error<IHistogram3D>(
"IHistogram3D");
303 IHistogram3D * createHistogram3D(
const std::string &,
int,
double,
double,
304 int,
double,
double,
int,
double,
double) {
305 return error<IHistogram3D>(
"IHistogram3D");
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");
322 IHistogram3D *
createCopy(
const std::string &,
const IHistogram3D &) {
323 return error<IHistogram3D>(
"IHistogram3D");
329 IProfile1D * createProfile1D(
const std::string &,
const std::string &,
330 int,
double,
double,
const std::string & =
"") {
331 return error<IProfile1D>(
"IProfile1D");
337 IProfile1D * createProfile1D(
const std::string &,
const std::string &,
338 int,
double,
double,
double,
double,
339 const std::string & =
"") {
340 return error<IProfile1D>(
"IProfile1D");
346 IProfile1D * createProfile1D(
const std::string &,
const std::string &,
347 const std::vector<double> &,
348 const std::string & =
"") {
349 return error<IProfile1D>(
"IProfile1D");
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");
365 return error<IProfile1D>(
"IProfile1D");
371 IProfile1D * createProfile1D(
const std::string &,
372 int,
double,
double,
double,
double) {
373 return error<IProfile1D>(
"IProfile1D");
379 IProfile1D *
createCopy(
const std::string &,
const IProfile1D &) {
380 return error<IProfile1D>(
"IProfile1D");
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");
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");
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");
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");
425 IProfile2D * createProfile2D(
const std::string &,
int,
double,
double,
426 int,
double,
double) {
427 return error<IProfile2D>(
"IProfile2D");
433 IProfile2D * createProfile2D(
const std::string &,
int,
double,
double,
434 int,
double,
double,
double,
double) {
435 return error<IProfile2D>(
"IProfile2D");
441 IProfile2D *
createCopy(
const std::string &,
const IProfile2D &) {
442 return error<IProfile2D>(
"IProfile2D");
458 if ( !checkBins(hist1, hist2) )
return 0;
460 h->
setTitle(path.substr(path.rfind(
'/') + 1));
462 if ( !tree->insert(path, h) )
return 0;
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));
496 if ( !checkBins(h1, h2) )
return 0;
498 h->
setTitle(path.substr(path.rfind(
'/') + 1));
499 for (
int i = 0; i < h->
ax->bins() + 2; ++i ) {
504 if ( !tree->insert(path, h) )
return 0;
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));
538 if ( !checkBins(h1, h2) )
return 0;
540 h->
setTitle(path.substr(path.rfind(
'/') + 1));
541 for (
int i = 0; i < h->
ax->bins() + 2; ++i ) {
546 if ( !tree->insert(path, h) )
return 0;
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));
580 if ( !checkBins(h1, h2) )
return 0;
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 ) {
594 if ( !tree->insert(path, h) )
return 0;
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));
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;
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;
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;
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)) )
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)) )
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));
674 if ( !checkBins(h1, h2) )
return 0;
676 h->
setTitle(path.substr(path.rfind(
'/') + 1));
678 if ( !tree->insert(path, h) ) {
690 if ( !checkBins(h1, h2) ) {
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];
706 if ( !tree->insert(path, h) ) {
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));
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));
737 if ( !checkBins(h1, h2) )
return 0;
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];
747 if ( !tree->insert(path, h) ) {
759 if ( !checkBins(h1,h2) )
return 0;
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 ) {
766 h->
sumw[ix][iy] = h->
sumw2[ix][iy] = 0.0;
774 if ( !tree->insert(path, h) ) {
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));
794 IHistogram3D * add(
const std::string &,
795 const IHistogram3D &,
const IHistogram3D &) {
796 return error<IHistogram3D>(
"3D histograms");
802 IHistogram3D * subtract(
const std::string &,
803 const IHistogram3D &,
const IHistogram3D &) {
804 return error<IHistogram3D>(
"3D histograms");
810 IHistogram3D * multiply(
const std::string &,
811 const IHistogram3D &,
const IHistogram3D &) {
812 return error<IHistogram3D>(
"3D histograms");
818 IHistogram3D * divide(
const std::string &,
819 const IHistogram3D &,
const IHistogram3D &) {
820 return error<IHistogram3D>(
"3D histograms");
827 IHistogram1D *
projectionX(
const std::string & path,
const IHistogram2D & h) {
828 return projectionX(path, dynamic_cast<const Histogram2D &>(h));
836 return sliceX(path, h, 0, h.
yax->bins() - 1);
843 IHistogram1D *
projectionY(
const std::string & path,
const IHistogram2D & h) {
844 return projectionY(path, dynamic_cast<const Histogram2D &>(h));
852 return sliceY(path, h, 0, h.
xax->bins() - 1);
860 sliceX(
const std::string & path,
const IHistogram2D & h,
int i) {
861 return sliceX(path, dynamic_cast<const Histogram2D &>(h), i, i);
870 return sliceX(path, h, i, i);
878 sliceY(
const std::string & path,
const IHistogram2D & h,
int i) {
879 return sliceY(path, dynamic_cast<const Histogram2D &>(h), i, i);
887 return sliceY(path, h, i, i);
895 sliceX(
const std::string & path,
const IHistogram2D & h,
int il,
int iu) {
896 return sliceX(path, dynamic_cast<const Histogram2D &>(h), il, iu);
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));
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];
924 if ( !tree->insert(path, h1) ) {
936 sliceY(
const std::string & path,
const IHistogram2D & h,
int il,
int iu) {
937 return sliceY(path, dynamic_cast<const Histogram2D &>(h), il, iu);
941 sliceY(
const std::string & path,
const Histogram2D & h2,
int il,
int iu) {
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));
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];
961 if ( !tree->insert(path, h1) ) {
972 IHistogram2D *
projectionXY(
const std::string &,
const IHistogram3D &) {
973 return error<IHistogram2D>(
"2D histograms");
980 IHistogram2D *
projectionXZ(
const std::string &,
const IHistogram3D &) {
981 return error<IHistogram2D>(
"2D histograms");
988 IHistogram2D *
projectionYZ(
const std::string &,
const IHistogram3D &) {
989 return error<IHistogram2D>(
"2D histograms");
997 IHistogram2D *
sliceXY(
const std::string &,
const IHistogram3D &,
int,
int) {
998 return error<IHistogram2D>(
"2D histograms");
1006 IHistogram2D *
sliceXZ(
const std::string &,
const IHistogram3D &,
int,
int) {
1007 return error<IHistogram2D>(
"2D histograms");
1015 IHistogram2D *
sliceYZ(
const std::string &,
const IHistogram3D &,
int,
int) {
1016 return error<IHistogram2D>(
"2D histograms");
1023 template <
typename T>
1025 throw std::runtime_error(
"LWH cannot handle " + feature +
".");