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;
void normalize(double intg)
Scale the given histogram so that the integral over all bins (including overflow) gives intg...
double binRmsY(int xindex, int yindex) const
The weighted y-RMS of a bin.
VariAxis * yvax
Pointer (possibly null) to a axis with fixed bin width.
double getSumXW(int xindex, int yindex) const
The weighted x-values.
Histogram2D(int nx, double lox, double upx, int ny, double loy, double upy)
Standard constructor.
virtual ~Histogram2D()
Destructor.
virtual double binHeightX(int index) const
Sum of all the heights of the bins along a given x bin.
VariAxis * xvax
Pointer (possibly null) to a axis with fixed bin width.
double binRmsX(int xindex, int yindex) const
The weighted x-RMS of a bin.
int binEntries(int xindex, int yindex) const
Number of entries in the corresponding bin (ie the number of times fill was called for this bin)...
const IAxis & yAxis() const
Get the y axis of the IHistogram2D.
double equivalentBinEntries() const
Number of equivalent entries, i.e.
bool setTitle(const std::string &title)
Set the histogram title.
std::string name() const
Get the Histogram's name.
std::vector< std::vector< double > > sumyw
The weighted y-values.
User level interface to 1D Histogram.
const IAnnotation & annotation() const
Not implemented in LWH.
std::string title() const
Get the Histogram's title.
double sumAllBinHeights() const
Sum of the heights of all the IHistogram's bins, i.e in-range bins, UNDERFLOW and OVERFLOW...
double getSumYW(int xindex, int yindex) const
The weighted x-values.
std::string theTitle
The title.
double rmsX() const
The RMS of the IHistogram2D along the x axis.
const IAxis & xAxis() const
Get the x axis of the IHistogram2D.
bool add(const Histogram2D &h)
Add to this Histogram2D the contents of another IHistogram2D.
double binError(int xindex, int yindex) const
The error of a given bin.
double getSumY2W(int xindex, int yindex) const
The weighted x-square-values.
int extraEntries() const
Number of entries in the UNDERFLOW and OVERFLOW bins.
double binMeanX(int xindex, int yindex) const
The weighted mean along the x-axis of a bin.
double binHeight(int xindex, int yindex) const
Total height of the corresponding bin (ie the sum of the weights in this bin).
double meanX() const
The mean of the IHistogram2D along the x axis.
std::vector< std::vector< double > > sumw2
The squared weights.
double binMeanY(int xindex, int yindex) const
The weighted mean along the y-axis of a bin.
Histogram2D(const Histogram2D &h)
Copy constructor.
std::vector< std::vector< double > > sumw
The weights.
double getSumX2W(int xindex, int yindex) const
The weighted x-square-values.
double sumExtraBinHeights() const
Sum of heights in the UNDERFLOW and OVERFLOW bins.
int entries() const
Get the number of in-range entries in the Histogram.
virtual int binEntriesY(int index) const
Sum of all the entries of the bins along a given y bin.
int coordToIndexY(double coord) const
Get the bin number corresponding to a given coordinate along the y axis.
bool add(const IHistogram2D &hist)
Add to this IHistogram1D the contents of another IHistogram1D.
bool writeFLAT(std::ostream &os, std::string path, std::string name)
Write out the histogram in a flat text file suitable for eg.
double minBinHeight() const
Minimum height of the in-range bins, i.e.
double rmsY() const
The RMS of the IHistogram2D along the x axis.
Histogram2D(const std::vector< double > &xedges, const std::vector< double > &yedges)
Standard constructor for variable bin width.
An VariAxis represents a binned histogram axis.
Axis * yfax
Pointer (possibly null) to a axis with fixed bin width.
IAnnotation & annotation()
Not implemented in LWH.
User level interface for factory classes of Histograms (binned, unbinned, and profile).
int coordToIndexX(double coord) const
Get the bin number corresponding to a given coordinate along the x axis.
std::vector< std::vector< double > > sumx2w
The weighted x-square-values.
int dimension() const
Get the Histogram's dimension.
An Axis represents a binned histogram axis.
std::vector< std::vector< int > > sum
The counts.
virtual double binHeightY(int index) const
Sum of all the heights of the bins along a given y bin.
std::vector< std::vector< double > > sumy2w
The weighted y-square-values.
bool scale(double s)
Scale the contents of this histogram with the given factor.
The LWH namespace contains a Light-Weight Histogram package which implements the most rudimentary his...
double getSumW(int xindex, int yindex) const
The weights.
bool writeXML(std::ostream &os, std::string path, std::string name)
Write out the histogram in the AIDA xml format.
double maxBinHeight() const
Maximum height of the in-range bins, i.e.
void * cast(const std::string &) const
Return the integral over the histogram bins assuming it has been normalize()d.
IAnnotation * anno
dummy pointer to non-existen annotation.
bool reset()
Reset the Histogram; as if just created.
bool fill(double x, double y, double weight=1.)
Fill the IHistogram1D with a value and the corresponding weight.
double sumBinHeights() const
Sum of in-range bin heights in the IHistogram, UNDERFLOW and OVERFLOW bins are excluded.
int allEntries() const
Sum of the entries in all the IHistogram's bins, i.e in-range bins, UNDERFLOW and OVERFLOW...
double meanY() const
The mean of the IHistogram2D along the y axis.
std::vector< std::vector< double > > sumxw
The weighted x-values.
virtual int binEntriesX(int index) const
Sum of all the entries of the bins along a given x bin.
double getSumW2(int xindex, int yindex) const
The squared weights.
Axis * xfax
Pointer (possibly null) to a axis with fixed bin width.