2 #ifndef LWH_Histogram2D_H
3 #define LWH_Histogram2D_H
8 #include "AIHistogram2D.h"
9 #include "ManagedObject.h"
38 int ny,
double loy,
double upy)
39 : xfax(new
Axis(nx, lox, upx)), xvax(0), yfax(new
Axis(ny, loy, upy)),
40 sum(nx + 2, std::vector<int>(ny + 2)),
41 sumw(nx + 2, std::vector<double>(ny + 2)),
42 sumw2(nx + 2, std::vector<double>(ny + 2)),
43 sumxw(nx + 2, std::vector<double>(ny + 2)),
44 sumx2w(nx + 2, std::vector<double>(ny + 2)),
45 sumyw(nx + 2, std::vector<double>(ny + 2)),
46 sumy2w(nx + 2, std::vector<double>(ny + 2)) {
55 const std::vector<double> & yedges)
56 : xfax(0), xvax(new
VariAxis(xedges)),
58 sum(xedges.size() + 1, std::vector<int>(yedges.size() + 1)),
59 sumw(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
60 sumw2(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
61 sumxw(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
62 sumx2w(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
63 sumyw(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
64 sumy2w(xedges.size() + 1, std::vector<double>(yedges.size() + 1)) {
73 : IBaseHistogram(h), IHistogram(h), IHistogram2D(h),
ManagedObject(h),
74 xfax(0), xvax(0), yfax(0), yvax(0),
75 sum(h.sum), sumw(h.sumw), sumw2(h.sumw2),
76 sumxw(h.sumxw), sumx2w(h.sumx2w) ,
77 sumyw(h.sumyw), sumy2w(h.sumy2w){
79 if ( hxvax ) xax = xvax =
new VariAxis(*hxvax);
80 else xax = xfax =
new Axis(dynamic_cast<const Axis &>(*h.
xax));
82 if ( hyvax ) yax = yvax =
new VariAxis(*hyvax);
83 else yax = yfax =
new Axis(dynamic_cast<const Axis &>(*h.
yax));
122 throw std::runtime_error(
"LWH cannot handle annotations");
130 throw std::runtime_error(
"LWH cannot handle annotations");
147 const int nx = xax->bins() + 2;
148 const int ny = yax->bins() + 2;
149 sum = std::vector< std::vector<int> >(nx, std::vector<int>(ny));
150 sumw = std::vector< std::vector<double> >(nx, std::vector<double>(ny));
166 for (
int ix = 2; ix < xax->bins() + 2; ++ix )
167 for (
int iy = 2; iy < yax->bins() + 2; ++iy ) si += sum[ix][iy];
179 return entries() + extraEntries();
187 int esum = sum[0][0] + sum[1][0] + sum[0][1] + sum[1][1];
188 for (
int ix = 2; ix < xax->bins() + 2; ++ix )
189 esum += sum[ix][0] + sum[ix][1];
190 for (
int iy = 2; iy < yax->bins() + 2; ++iy )
191 esum += sum[0][iy] + sum[1][iy];
203 for (
int ix = 2; ix < xax->bins() + 2; ++ix )
204 for (
int iy = 2; iy < yax->bins() + 2; ++iy ) {
206 sw2 += sumw2[ix][iy];
219 for (
int ix = 2; ix < xax->bins() + 2; ++ix )
220 for (
int iy = 2; iy < yax->bins() + 2; ++iy ) sw += sumw[ix][iy];
230 return sumBinHeights() + sumExtraBinHeights();
238 int esum = sumw[0][0] + sumw[1][0] + sumw[0][1] + sumw[1][1];
239 for (
int ix = 2; ix < xax->bins() + 2; ++ix )
240 esum += sumw[ix][0] + sumw[ix][1];
241 for (
int iy = 2; iy < yax->bins() + 2; ++iy )
242 esum += sumw[0][iy] + sumw[1][iy];
252 double minw = sumw[2][2];
253 for (
int ix = 2; ix < xax->bins() + 2; ++ix )
254 for (
int iy = 2; iy < yax->bins() + 2; ++iy )
255 minw = std::min(minw, sumw[ix][iy]);
265 double maxw = sumw[2][2];
266 for (
int ix = 2; ix < xax->bins() + 2; ++ix )
267 for (
int iy = 2; iy < yax->bins() + 2; ++iy )
268 maxw = std::max(maxw, sumw[ix][iy]);
279 bool fill(
double x,
double y,
double weight = 1.) {
280 int ix = xax->coordToIndex(x) + 2;
281 int iy = yax->coordToIndex(y) + 2;
283 sumw[ix][iy] += weight;
284 sumxw[ix][iy] += x*weight;
285 sumx2w[ix][iy] += x*x*weight;
286 sumyw[ix][iy] += y*weight;
287 sumy2w[ix][iy] += y*y*weight;
288 sumw2[ix][iy] += weight*weight;
289 return weight >= 0 && weight <= 1;
301 return sumw[ix][iy] != 0.0? sumxw[ix][iy]/sumw[ix][iy]:
302 ( xvax? xvax->binMidPoint(xindex): xfax->binMidPoint(xindex) );
314 return sumw[ix][iy] != 0.0? sumyw[ix][iy]/sumw[ix][iy]:
315 ( yvax? yvax->binMidPoint(yindex): xfax->binMidPoint(yindex) );
324 double binRmsX(
int xindex,
int yindex)
const {
327 return sumw[ix][iy] == 0.0 || sum[ix][iy] < 2? xax->binWidth(xindex):
328 std::sqrt(std::max(sumw[ix][iy]*sumx2w[ix][iy] -
329 sumxw[ix][iy]*sumxw[ix][iy], 0.0))/sumw[ix][iy];
338 double binRmsY(
int xindex,
int yindex)
const {
341 return sumw[ix][iy] == 0.0 || sum[ix][iy] < 2? yax->binWidth(yindex):
342 std::sqrt(std::max(sumw[ix][iy]*sumy2w[ix][iy] -
343 sumyw[ix][iy]*sumyw[ix][iy], 0.0))/sumw[ix][iy];
353 return sum[xindex + 2][yindex + 2];
365 for (
int iy = 2; iy < yax->bins() + 2; ++iy )
366 ret += sum[index + 2][iy];
379 for (
int ix = 2; ix < xax->bins() + 2; ++ix )
380 ret += sum[ix][index + 2];
393 return sumw[xindex + 2][yindex + 2];
405 for (
int iy = 2; iy < yax->bins() + 2; ++iy )
406 ret += sumw[index + 2][iy];
419 for (
int ix = 2; ix < xax->bins() + 2; ++ix )
420 ret += sumw[ix][index + 2];
431 return std::sqrt(sumw2[xindex + 2][yindex + 2]);
442 for (
int ix = 2; ix < xax->bins() + 2; ++ix )
443 for (
int iy = 2; iy < yax->bins() + 2; ++iy ) {
447 return s != 0.0? sx/s: 0.0;
458 for (
int ix = 2; ix < xax->bins() + 2; ++ix )
459 for (
int iy = 2; iy < yax->bins() + 2; ++iy ) {
463 return s != 0.0? sy/s: 0.0;
475 for (
int ix = 2; ix < xax->bins() + 2; ++ix )
476 for (
int iy = 2; iy < yax->bins() + 2; ++iy ) {
479 sx2 += sumx2w[ix][iy];
481 return s != 0.0? std::sqrt(std::max(s*sx2 - sx*sx, 0.0))/s:
482 xax->upperEdge() - xax->lowerEdge();
494 for (
int ix = 2; ix < xax->bins() + 2; ++ix )
495 for (
int iy = 2; iy < yax->bins() + 2; ++iy ) {
498 sy2 += sumy2w[ix][iy];
500 return s != 0.0? std::sqrt(std::max(s*sy2 - sy*sy, 0.0))/s:
501 yax->upperEdge() - yax->lowerEdge();
505 double getSumW(
int xindex,
int yindex)
const {
506 return sumw[xindex + 2][yindex + 2];
511 return sumw2[xindex + 2][yindex + 2];
516 return sumxw[xindex + 2][yindex + 2];
521 return sumx2w[xindex + 2][yindex + 2];
526 return sumyw[xindex + 2][yindex + 2];
531 return sumy2w[xindex + 2][yindex + 2];
558 return xax->coordToIndex(coord);
569 return yax->coordToIndex(coord);
578 if ( xax->upperEdge() != h.
xax->upperEdge() ||
579 xax->lowerEdge() != h.
xax->lowerEdge() ||
580 xax->bins() != h.
xax->bins() )
return false;
581 if ( yax->upperEdge() != h.
yax->upperEdge() ||
582 yax->lowerEdge() != h.
yax->lowerEdge() ||
583 yax->bins() != h.
yax->bins() )
return false;
584 for (
int ix = 0; ix < xax->bins() + 2; ++ix )
585 for (
int iy = 0; iy < yax->bins() + 2; ++iy ) {
586 sum[ix][iy] += h.
sum[ix][iy];
587 sumw[ix][iy] += h.
sumw[ix][iy];
588 sumxw[ix][iy] += h.
sumxw[ix][iy];
589 sumx2w[ix][iy] += h.
sumx2w[ix][iy];
590 sumyw[ix][iy] += h.
sumyw[ix][iy];
591 sumy2w[ix][iy] += h.
sumy2w[ix][iy];
592 sumw2[ix][iy] += h.
sumw2[ix][iy];
602 bool add(
const IHistogram2D & hist) {
603 return add(dynamic_cast<const Histogram2D &>(hist));
611 for (
int ix = 0; ix < xax->bins() + 2; ++ix )
612 for (
int iy = 0; iy < yax->bins() + 2; ++iy ) {
618 sumw2[ix][iy] *= s*s;
631 double oldintg = sumAllBinHeights();
632 if ( oldintg == 0.0 )
return;
633 for (
int ix = 0; ix < xax->bins() + 2; ++ix )
634 for (
int iy = 0; iy < yax->bins() + 2; ++iy ) {
635 double fac = intg/oldintg;
636 if ( ix >= 2 && iy >= 2 )
637 fac /= (xax->binUpperEdge(ix - 2) - xax->binLowerEdge(ix - 2))*
638 (yax->binUpperEdge(iy - 2) - yax->binLowerEdge(iy - 2));
640 sumxw[ix][iy] *= fac;
641 sumx2w[ix][iy] *= fac;
642 sumyw[ix][iy] *= fac;
643 sumy2w[ix][iy] *= fac;
644 sumw2[ix][iy] *= fac*fac;
666 void *
cast(
const std::string &)
const {
673 bool writeXML(std::ostream & os, std::string path, std::string name) {
675 os <<
" <histogram2d name=\"" << name
676 <<
"\"\n title=\"" << title()
677 <<
"\" path=\"" << path
678 <<
"\">\n <axis max=\"" << xax->upperEdge()
679 <<
"\" numberOfBins=\"" << xax->bins()
680 <<
"\" min=\"" << xax->lowerEdge()
681 <<
"\" direction=\"x\"";
684 for (
int i = 0, N = xax->bins() - 1; i < N; ++i )
685 os <<
" <binBorder value=\"" << xax->binUpperEdge(i) <<
"\"/>\n";
690 os <<
" <axis max=\"" << yax->upperEdge()
691 <<
"\" numberOfBins=\"" << yax->bins()
692 <<
"\" min=\"" << yax->lowerEdge()
693 <<
"\" direction=\"y\"";
696 for (
int i = 0, N = yax->bins() - 1; i < N; ++i )
697 os <<
" <binBorder value=\"" << yax->binUpperEdge(i) <<
"\"/>\n";
702 os <<
" <statistics entries=\"" << entries()
703 <<
"\">\n <statistic mean=\"" << meanX()
704 <<
"\" direction=\"x\"\n rms=\"" << rmsX()
705 <<
"\"/>\n </statistics>\n <statistic mean=\"" << meanY()
706 <<
"\" direction=\"y\"\n rms=\"" << rmsY()
707 <<
"\"/>\n </statistics>\n <data2d>\n";
708 for (
int ix = 0; ix < xax->bins() + 2; ++ix )
709 for (
int iy = 0; iy < yax->bins() + 2; ++iy )
711 os <<
" <bin2d binNumX=\"";
712 if ( ix == 0 ) os <<
"UNDERFLOW";
713 else if ( ix == 1 ) os <<
"OVERFLOW";
715 os <<
"\" binNumY=\"";
716 if ( iy == 0 ) os <<
"UNDERFLOW";
717 else if ( iy == 1 ) os <<
"OVERFLOW";
719 os <<
"\" entries=\"" << sum[ix][iy]
720 <<
"\" height=\"" << sumw[ix][iy]
721 <<
"\"\n error=\"" << std::sqrt(sumw2[ix][iy])
722 <<
"\" error2=\"" << sumw2[ix][iy]
723 <<
"\"\n weightedMeanX=\"" << binMeanX(ix - 2, iy - 2)
724 <<
"\" weightedRmsX=\"" << binRmsX(ix - 2, iy - 2)
725 <<
"\"\n weightedMeanY=\"" << binMeanY(ix - 2, iy - 2)
726 <<
"\" weightedRmsY=\"" << binRmsY(ix - 2, iy - 2)
729 os <<
" </data2d>\n </histogram2d>" << std::endl;
738 bool writeFLAT(std::ostream & os, std::string path, std::string name) {
739 os <<
"#2D " << path <<
"/" << name <<
" " << xax->lowerEdge()
740 <<
" " << xax->bins() <<
" " << xax->upperEdge() <<
" "
741 << yax->lowerEdge() <<
" " << yax->bins() <<
" " << yax->upperEdge()
742 <<
" \"" << title() <<
"\"" << std::endl;
743 for (
int ix = 2; ix < xax->bins() + 2; ++ix ) {
744 for (
int iy = 2; iy < yax->bins() + 2; ++iy )
745 os << 0.5*(xax->binLowerEdge(ix - 2)+xax->binUpperEdge(ix - 2)) <<
" "
746 << 0.5*(yax->binLowerEdge(iy - 2)+yax->binUpperEdge(iy - 2))
747 <<
" " << sumw[ix][iy] <<
" " << sqrt(sumw2[ix][iy])
748 <<
" " << sum[ix][iy] << std::endl;
762 bool writeROOT(TFile* file, std::string path, std::string name) {
768 if (!vax || vax->isFixedBinning() ) {
770 hist1d =
new TH1D(name.c_str(), title().c_str(), nbins, ax->lowerEdge(), ax->upperEdge());
774 double* bins =
new double[nbins+1];
775 for (
int i=0; i<nbins; ++i) {
776 bins[ix][iy] = vax->binEdges(i).first;
778 bins[nbins] = vax->binEdges(nbins-1).second;
779 hist1d =
new TH1D(name.c_str(), title().c_str(), nbins, bins);
785 for (
int i = 0; i < nbins + 2; ++i ) {
788 entries = entries + sum[ix][iy];
791 else if (i==1) j=nbins+1;
793 hist1d->SetBinContent(j, sumw[ix][iy]);
794 hist1d->SetBinError(j, sqrt(sumw2[ix][iy]));
800 hist1d->SetEntries(entries);
803 for (
unsigned int i=1; i<path.size(); ++i) DirName += path[i];
804 if (!file->Get(DirName.c_str())) file->mkdir(DirName.c_str());
805 file->cd(DirName.c_str());
841 std::vector< std::vector<int> >
sum;
844 std::vector< std::vector<double> >
sumw;
847 std::vector< std::vector<double> >
sumw2;
850 std::vector< std::vector<double> >
sumxw;
853 std::vector< std::vector<double> >
sumx2w;
856 std::vector< std::vector<double> >
sumyw;
859 std::vector< std::vector<double> >
sumy2w;