Analysis/LWH/Histogram2D.h

0001 // -*- C++ -*-
0002 #ifndef LWH_Histogram2D_H
0003 #define LWH_Histogram2D_H
0004 //
0005 // This is the declaration of the Histogram1D class.
0006 //
0007
0008 #include "AIHistogram2D.h"
0009 #include "ManagedObject.h"
0010 #include "Axis.h"
0011 #include "VariAxis.h"
0012 #include <vector>
0013 #include <stdexcept>
0014
0015 #include <iostream>
0016
0017 namespace LWH {
0018
0019   using namespace AIDA;
0020
0021
0022   /**
0023    * User level interface to 1D Histogram.
0024    */
0025   class Histogram2D: public IHistogram2D, public ManagedObject {
0026
0027   public:
0028
0029     /** HistFactory is a friend. */
0030     friend class HistogramFactory;
0031
0032   public:
0033
0034     /**
0035      * Standard constructor.
0036      */
0037     Histogram2D(int nx, double lox, double upx,
0038         int ny, double loy, double upy)
0039       : xfax(new Axis(nx, lox, upx)), xvax(0), yfax(new Axis(ny, loy, upy)),
0040     sum(nx + 2, std::vector<int>(ny + 2)),
0041     sumw(nx + 2, std::vector<double>(ny + 2)),
0042     sumw2(nx + 2, std::vector<double>(ny + 2)),
0043     sumxw(nx + 2, std::vector<double>(ny + 2)),
0044     sumx2w(nx + 2, std::vector<double>(ny + 2)),
0045     sumyw(nx + 2, std::vector<double>(ny + 2)),
0046     sumy2w(nx + 2, std::vector<double>(ny + 2)) {
0047       xax = xfax;
0048       yax = yfax;
0049     }
0050
0051     /**
0052      * Standard constructor for variable bin width.
0053      */
0054     Histogram2D(const std::vector<double> & xedges,
0055         const std::vector<double> & yedges)
0056       : xfax(0), xvax(new VariAxis(xedges)),
0057     yfax(0), yvax(new VariAxis(xedges)),
0058         sum(xedges.size() + 1, std::vector<int>(yedges.size() + 1)),
0059     sumw(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
0060     sumw2(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
0061         sumxw(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
0062     sumx2w(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
0063         sumyw(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
0064     sumy2w(xedges.size() + 1, std::vector<double>(yedges.size() + 1)) {
0065       xax = xvax;
0066       yax = yvax;
0067     }
0068
0069     /**
0070      * Copy constructor.
0071      */
0072     Histogram2D(const Histogram2D & h)
0073       : IBaseHistogram(h), IHistogram(h), IHistogram2D(h), ManagedObject(h),
0074         xfax(0), xvax(0),  yfax(0), yvax(0),
0075     sum(h.sum), sumw(h.sumw), sumw2(h.sumw2),
0076         sumxw(h.sumxw), sumx2w(h.sumx2w) ,
0077         sumyw(h.sumyw), sumy2w(h.sumy2w){
0078       const VariAxis * hxvax = dynamic_cast<const VariAxis *>(h.xax);
0079       if ( hxvax ) xax = xvax = new VariAxis(*hxvax);
0080       else xax = xfax = new Axis(dynamic_cast<const Axis &>(*h.xax));
0081       const VariAxis * hyvax = dynamic_cast<const VariAxis *>(h.yax);
0082       if ( hyvax ) yax = yvax = new VariAxis(*hyvax);
0083       else yax = yfax = new Axis(dynamic_cast<const Axis &>(*h.yax));
0084   }
0085
0086     /// Destructor.
0087     virtual ~Histogram2D() {
0088       delete xax;
0089       delete yax;
0090     }
0091
0092     /**
0093      * Get the Histogram's title.
0094      * @return The Histogram's title.
0095      */
0096     std::string title() const {
0097       return theTitle;
0098     }
0099
0100     /**
0101      * Get the Histogram's name.
0102      * @return The Histogram's name
0103      */
0104     std::string name() const {
0105       return title();
0106     }
0107
0108     /**
0109      * Set the histogram title.
0110      * @param title The title.
0111      * @return false If title cannot be changed.
0112      */
0113     bool setTitle(const std::string & title) {
0114       theTitle = title;
0115       return true;
0116     }
0117
0118     /**
0119      * Not implemented in LWH. will throw an exception.
0120      */
0121     IAnnotation & annotation() {
0122       throw std::runtime_error("LWH cannot handle annotations");
0123       return *anno;
0124     }
0125
0126     /**
0127      * Not implemented in LWH. will throw an exception.
0128      */
0129     const IAnnotation & annotation() const {
0130       throw std::runtime_error("LWH cannot handle annotations");
0131       return *anno;
0132     }
0133
0134     /**
0135      * Get the Histogram's dimension.
0136      * @return The Histogram's dimension.
0137      */
0138     int dimension() const {
0139       return 2;
0140     }
0141
0142     /**
0143      * Reset the Histogram; as if just created.
0144      * @return false If something goes wrong.
0145      */
0146     bool reset() {
0147       const int nx = xax->bins() + 2;
0148       const int ny = yax->bins() + 2;
0149       sum = std::vector< std::vector<int> >(nx, std::vector<int>(ny));
0150       sumw = std::vector< std::vector<double> >(nx, std::vector<double>(ny));
0151       sumw2 = sumw;
0152       sumxw = sumw;
0153       sumx2w = sumw;
0154       sumyw = sumw;
0155       sumy2w = sumw;
0156       return true;
0157     }
0158
0159     /**
0160      * Get the number of in-range entries in the Histogram.
0161      * @return The number of in-range entries.
0162      *
0163      */
0164     int entries() const {
0165       int si = 0;
0166       for ( int ix = 2; ix < xax->bins() + 2; ++ix )
0167     for ( int iy = 2; iy < yax->bins() + 2; ++iy ) si += sum[ix][iy];
0168       return si;
0169     }
0170
0171     /**
0172      * Sum of the entries in all the IHistogram's bins,
0173      * i.e in-range bins, UNDERFLOW and OVERFLOW.
0174      * This is equivalent to the number of times the
0175      * method fill was invoked.
0176      * @return The sum of all the entries.
0177      */
0178     int allEntries() const {
0179       return entries() + extraEntries();
0180     }
0181
0182     /**
0183      * Number of entries in the UNDERFLOW and OVERFLOW bins.
0184      * @return The number of entries outside the range of the IHistogram.
0185      */
0186     int extraEntries() const {
0187       int esum = sum[0][0] + sum[1][0] + sum[0][1] + sum[1][1];
0188       for ( int ix = 2; ix < xax->bins() + 2; ++ix )
0189     esum += sum[ix][0] + sum[ix][1];
0190       for ( int iy = 2; iy < yax->bins() + 2; ++iy )
0191     esum += sum[0][iy] + sum[1][iy];
0192       return esum;
0193     }
0194
0195     /**
0196      * Number of equivalent entries,
0197      * i.e. <tt>SUM[ weight ] ^ 2 / SUM[ weight^2 ]</tt>
0198      * @return The number of equivalent entries.
0199      */
0200     double equivalentBinEntries() const {
0201       double sw = 0.0;
0202       double sw2 = 0.0;
0203       for ( int ix = 2; ix < xax->bins() + 2; ++ix )
0204     for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
0205       sw += sumw[ix][iy];
0206       sw2 += sumw2[ix][iy];
0207     }
0208       return sw2/(sw*sw);
0209     }
0210
0211     /**
0212      * Sum of in-range bin heights in the IHistogram,
0213      * UNDERFLOW and OVERFLOW bins are excluded.
0214      * @return The sum of the in-range bins heights.
0215      *
0216      */
0217     double sumBinHeights() const {
0218       double sw = 0.0;
0219       for ( int ix = 2; ix < xax->bins() + 2; ++ix )
0220     for ( int iy = 2; iy < yax->bins() + 2; ++iy ) sw += sumw[ix][iy];
0221       return sw;
0222     }
0223
0224     /**
0225      * Sum of the heights of all the IHistogram's bins,
0226      * i.e in-range bins, UNDERFLOW and OVERFLOW.
0227      * @return The sum of all the bins heights.
0228      */
0229     double sumAllBinHeights() const {
0230       return sumBinHeights() + sumExtraBinHeights();
0231     }
0232
0233     /**
0234      * Sum of heights in the UNDERFLOW and OVERFLOW bins.
0235      * @return The sum of the heights of the out-of-range bins.
0236      */
0237     double sumExtraBinHeights() const {
0238       int esum = sumw[0][0] + sumw[1][0] + sumw[0][1] + sumw[1][1];
0239       for ( int ix = 2; ix < xax->bins() + 2; ++ix )
0240     esum += sumw[ix][0] + sumw[ix][1];
0241       for ( int iy = 2; iy < yax->bins() + 2; ++iy )
0242     esum += sumw[0][iy] + sumw[1][iy];
0243       return esum;
0244     }
0245
0246     /**
0247      * Minimum height of the in-range bins,
0248      * i.e. not considering the UNDERFLOW and OVERFLOW bins.
0249      * @return The minimum height among the in-range bins.
0250      */
0251     double minBinHeight() const {
0252       double minw = sumw[2][2];
0253       for ( int ix = 2; ix < xax->bins() + 2; ++ix )
0254     for ( int iy = 2; iy < yax->bins() + 2; ++iy )
0255       minw = std::min(minw, sumw[ix][iy]);
0256       return minw;
0257     }
0258
0259     /**
0260      * Maximum height of the in-range bins,
0261      * i.e. not considering the UNDERFLOW and OVERFLOW bins.
0262      * @return The maximum height among the in-range bins.
0263      */
0264     double maxBinHeight() const{
0265       double maxw = sumw[2][2];
0266       for ( int ix = 2; ix < xax->bins() + 2; ++ix )
0267     for ( int iy = 2; iy < yax->bins() + 2; ++iy )
0268       maxw = std::max(maxw, sumw[ix][iy]);
0269       return maxw;
0270     }
0271
0272     /**
0273      * Fill the IHistogram1D with a value and the
0274      * corresponding weight.
0275      * @param x      The value to be filled in.
0276      * @param weight The corresponding weight (by default 1).
0277      * @return false If the weight is <0 or >1 (?).
0278      */
0279     bool fill(double x, double y, double weight = 1.) {
0280       int ix = xax->coordToIndex(x) + 2;
0281       int iy = yax->coordToIndex(y) + 2;
0282       ++sum[ix][iy];
0283       sumw[ix][iy] += weight;
0284       sumxw[ix][iy] += x*weight;
0285       sumx2w[ix][iy] += x*x*weight;
0286       sumyw[ix][iy] += y*weight;
0287       sumy2w[ix][iy] += y*y*weight;
0288       sumw2[ix][iy] += weight*weight;
0289       return weight >= 0 && weight <= 1;
0290     }
0291
0292     /**
0293      * The weighted mean along the x-axis of a bin.
0294      * @param xindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0295      * @param yindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0296      * @return      The mean of the corresponding bin.
0297      */
0298     double binMeanX(int xindex, int yindex) const {
0299       int ix = xindex + 2;
0300       int iy = yindex + 2;
0301       return sumw[ix][iy] != 0.0? sumxw[ix][iy]/sumw[ix][iy]:
0302         ( xvax? xvax->binMidPoint(xindex): xfax->binMidPoint(xindex) );
0303     };
0304
0305     /**
0306      * The weighted mean along the y-axis of a bin.
0307      * @param xindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0308      * @param yindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0309      * @return      The mean of the corresponding bin.
0310      */
0311     double binMeanY(int xindex, int yindex) const {
0312       int ix = xindex + 2;
0313       int iy = yindex + 2;
0314       return sumw[ix][iy] != 0.0? sumyw[ix][iy]/sumw[ix][iy]:
0315         ( yvax? yvax->binMidPoint(yindex): xfax->binMidPoint(yindex) );
0316     };
0317
0318     /**
0319      * The weighted x-RMS of a bin.
0320      * @param xindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0321      * @param yindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0322      * @return      The RMS of the corresponding bin.
0323      */
0324     double binRmsX(int xindex, int yindex) const {
0325       int ix = xindex + 2;
0326       int iy = yindex + 2;
0327       return sumw[ix][iy] == 0.0 || sum[ix][iy] < 2? xax->binWidth(xindex):
0328         std::sqrt(std::max(sumw[ix][iy]*sumx2w[ix][iy] -
0329                sumxw[ix][iy]*sumxw[ix][iy], 0.0))/sumw[ix][iy];
0330     };
0331
0332     /**
0333      * The weighted y-RMS of a bin.
0334      * @param xindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0335      * @param yindex The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0336      * @return      The RMS of the corresponding bin.
0337      */
0338     double binRmsY(int xindex, int yindex) const {
0339       int ix = xindex + 2;
0340       int iy = yindex + 2;
0341       return sumw[ix][iy] == 0.0 || sum[ix][iy] < 2? yax->binWidth(yindex):
0342         std::sqrt(std::max(sumw[ix][iy]*sumy2w[ix][iy] -
0343                sumyw[ix][iy]*sumyw[ix][iy], 0.0))/sumw[ix][iy];
0344     };
0345
0346     /**
0347      * Number of entries in the corresponding bin (ie the number of
0348      * times fill was called for this bin).
0349      * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0350      * @return      The number of entries in the corresponding bin.
0351      */
0352     int binEntries(int xindex, int yindex) const {
0353       return sum[xindex + 2][yindex + 2];
0354     }
0355
0356     /**
0357      * Sum of all the entries of the bins along a given x bin.
0358      * This is equivalent to <tt>projectionX().binEntries(index)</tt>.
0359      * @param index The x bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0360      * @return      The number of entries in the corresponding set of bins.
0361      *
0362      */
0363     virtual int binEntriesX(int index) const {
0364       int ret = 0;
0365       for ( int iy = 2; iy < yax->bins() + 2; ++iy )
0366     ret += sum[index + 2][iy];
0367       return ret;
0368     }
0369
0370     /**
0371      * Sum of all the entries of the bins along a given y bin.
0372      * This is equivalent to <tt>projectionY().binEntries(index)</tt>.
0373      * @param index The y bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0374      * @return      The number of entries in the corresponding set of bins.
0375      *
0376      */
0377     virtual int binEntriesY(int index) const {
0378       int ret = 0;
0379       for ( int ix = 2; ix < xax->bins() + 2; ++ix )
0380     ret += sum[ix][index + 2];
0381       return ret;
0382     }
0383
0384     /**
0385      * Total height of the corresponding bin (ie the sum of the weights
0386      * in this bin).
0387      * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0388      * @return      The height of the corresponding bin.
0389      */
0390     double binHeight(int xindex, int yindex) const {
0391       /// @todo While this is compatible with the reference AIDA
0392       /// implementation, it is not the bin height!
0393       return sumw[xindex + 2][yindex + 2];
0394     }
0395
0396     /**
0397      * Sum of all the heights of the bins along a given x bin.
0398      * This is equivalent to <tt>projectionX().binHeight(index)</tt>.
0399      * @param index The x bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0400      * @return      The sum of the heights in the corresponding set of bins.
0401      *
0402      */
0403     virtual double binHeightX(int index) const {
0404       double ret = 0;
0405       for ( int iy = 2; iy < yax->bins() + 2; ++iy )
0406     ret += sumw[index + 2][iy];
0407       return ret;
0408     }
0409
0410     /**
0411      * Sum of all the heights of the bins along a given y bin.
0412      * This is equivalent to <tt>projectionY().binHeight(index)</tt>.
0413      * @param index The y bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0414      * @return      The sum of the heights in the corresponding set of bins.
0415      *
0416      */
0417     virtual double binHeightY(int index) const {
0418       double ret = 0;
0419       for ( int ix = 2; ix < xax->bins() + 2; ++ix )
0420     ret += sumw[ix][index + 2];
0421       return ret;
0422     }
0423
0424     /**
0425      * The error of a given bin.
0426      * @param index The bin number (0...N-1) or OVERFLOW or UNDERFLOW.
0427      * @return      The error on the corresponding bin.
0428      *
0429      */
0430     double binError(int xindex, int yindex) const {
0431       return std::sqrt(sumw2[xindex + 2][yindex + 2]);
0432     }
0433
0434     /**
0435      * The mean of the IHistogram2D along the x axis.
0436      * @return The mean of the IHistogram2D along the x axis.
0437      *
0438      */
0439     double meanX() const {
0440       double s = 0.0;
0441       double sx = 0.0;
0442       for ( int ix = 2; ix < xax->bins() + 2; ++ix )
0443     for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
0444         s += sumw[ix][iy];
0445         sx += sumxw[ix][iy];
0446       }
0447       return s != 0.0? sx/s: 0.0;
0448     }
0449
0450     /**
0451      * The mean of the IHistogram2D along the y axis.
0452      * @return The mean of the IHistogram2D along the y axis.
0453      *
0454      */
0455     double meanY() const {
0456       double s = 0.0;
0457       double sy = 0.0;
0458       for ( int ix = 2; ix < xax->bins() + 2; ++ix )
0459     for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
0460         s += sumw[ix][iy];
0461         sy += sumyw[ix][iy];
0462       }
0463       return s != 0.0? sy/s: 0.0;
0464     }
0465
0466     /**
0467      * The RMS of the IHistogram2D along the x axis.
0468      * @return The RMS if the IHistogram2D along the x axis.
0469      *
0470      */
0471     double rmsX() const {
0472       double s = 0.0;
0473       double sx = 0.0;
0474       double sx2 = 0.0;
0475       for ( int ix = 2; ix < xax->bins() + 2; ++ix )
0476     for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
0477         s += sumw[ix][iy];
0478         sx += sumxw[ix][iy];
0479         sx2 += sumx2w[ix][iy];
0480       }
0481       return s != 0.0? std::sqrt(std::max(s*sx2 - sx*sx, 0.0))/s:
0482         xax->upperEdge() - xax->lowerEdge();
0483     }
0484
0485     /**
0486      * The RMS of the IHistogram2D along the x axis.
0487      * @return The RMS if the IHistogram2D along the x axis.
0488      *
0489      */
0490     double rmsY() const {
0491       double s = 0.0;
0492       double sy = 0.0;
0493       double sy2 = 0.0;
0494       for ( int ix = 2; ix < xax->bins() + 2; ++ix )
0495     for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
0496         s += sumw[ix][iy];
0497         sy += sumyw[ix][iy];
0498         sy2 += sumy2w[ix][iy];
0499       }
0500       return s != 0.0? std::sqrt(std::max(s*sy2 - sy*sy, 0.0))/s:
0501         yax->upperEdge() - yax->lowerEdge();
0502     }
0503
0504     /** The weights. */
0505     double getSumW(int xindex, int yindex) const {
0506         return sumw[xindex + 2][yindex + 2];
0507     }
0508
0509     /** The squared weights. */
0510     double getSumW2(int xindex, int yindex) const {
0511         return sumw2[xindex + 2][yindex + 2];
0512     }
0513
0514     /** The weighted x-values. */
0515     double getSumXW(int xindex, int yindex) const {
0516         return sumxw[xindex + 2][yindex + 2];
0517     }
0518
0519     /** The weighted x-square-values. */
0520     double getSumX2W(int xindex, int yindex) const {
0521         return sumx2w[xindex + 2][yindex + 2];
0522     }
0523
0524     /** The weighted x-values. */
0525     double getSumYW(int xindex, int yindex) const {
0526         return sumyw[xindex + 2][yindex + 2];
0527     }
0528
0529     /** The weighted x-square-values. */
0530     double getSumY2W(int xindex, int yindex) const {
0531         return sumy2w[xindex + 2][yindex + 2];
0532     }
0533
0534     /**
0535      * Get the x axis of the IHistogram2D.
0536      * @return The x coordinate IAxis.
0537      */
0538     const IAxis & xAxis() const {
0539       return *xax;
0540     }
0541
0542     /**
0543      * Get the y axis of the IHistogram2D.
0544      * @return The y coordinate IAxis.
0545      */
0546     const IAxis & yAxis() const {
0547       return *yax;
0548     }
0549
0550     /**
0551      * Get the bin number corresponding to a given coordinate along the
0552      * x axis.  This is a convenience method, equivalent to
0553      * <tt>axis().coordToIndex(coord)</tt>.
0554      * @param coord The coordinalte along the x axis.
0555      * @return      The corresponding bin number.
0556      */
0557     int coordToIndexX(double coord) const {
0558       return xax->coordToIndex(coord);
0559     }
0560
0561     /**
0562      * Get the bin number corresponding to a given coordinate along the
0563      * y axis.  This is a convenience method, equivalent to
0564      * <tt>axis().coordToIndex(coord)</tt>.
0565      * @param coord The coordinalte along the y axis.
0566      * @return      The corresponding bin number.
0567      */
0568     int coordToIndexY(double coord) const {
0569       return yax->coordToIndex(coord);
0570     }
0571
0572     /**
0573      * Add to this Histogram2D the contents of another IHistogram2D.
0574      * @param h The Histogram2D to be added to this IHistogram2D.
0575      * @return false If the IHistogram1Ds binnings are incompatible.
0576      */
0577     bool add(const Histogram2D & h) {
0578       if ( xax->upperEdge() != h.xax->upperEdge() ||
0579        xax->lowerEdge() != h.xax->lowerEdge() ||
0580        xax->bins() != h.xax->bins() ) return false;
0581       if ( yax->upperEdge() != h.yax->upperEdge() ||
0582        yax->lowerEdge() != h.yax->lowerEdge() ||
0583        yax->bins() != h.yax->bins() ) return false;
0584       for ( int ix = 0; ix < xax->bins() + 2; ++ix )
0585     for ( int iy = 0; iy < yax->bins() + 2; ++iy ) {
0586       sum[ix][iy] += h.sum[ix][iy];
0587       sumw[ix][iy] += h.sumw[ix][iy];
0588       sumxw[ix][iy] += h.sumxw[ix][iy];
0589       sumx2w[ix][iy] += h.sumx2w[ix][iy];
0590       sumyw[ix][iy] += h.sumyw[ix][iy];
0591       sumy2w[ix][iy] += h.sumy2w[ix][iy];
0592       sumw2[ix][iy] += h.sumw2[ix][iy];
0593     }
0594       return true;
0595     }
0596
0597     /**
0598      * Add to this IHistogram1D the contents of another IHistogram1D.
0599      * @param hist The IHistogram1D to be added to this IHistogram1D.
0600      * @return false If the IHistogram1Ds binnings are incompatible.
0601      */
0602     bool add(const IHistogram2D & hist) {
0603       return add(dynamic_cast<const Histogram2D &>(hist));
0604     }
0605
0606     /**
0607      * Scale the contents of this histogram with the given factor.
0608      * @param s the scaling factor to use.
0609      */
0610     bool scale(double s) {
0611       for ( int ix = 0; ix < xax->bins() + 2; ++ix )
0612     for ( int iy = 0; iy < yax->bins() + 2; ++iy ) {
0613       sumw[ix][iy] *= s;
0614       sumxw[ix][iy] *= s;
0615       sumx2w[ix][iy] *= s;
0616       sumyw[ix][iy] *= s;
0617       sumy2w[ix][iy] *= s;
0618       sumw2[ix][iy] *= s*s;
0619       }
0620       return true;
0621     }
0622
0623     /**
0624      * Scale the given histogram so that the integral over all bins
0625      * (including overflow) gives \a intg. This function also corrects
0626      * for the bin-widths, which means that it should only be run once
0627      * for each histogram. Further rescaling must be done with the
0628      * scale(double) function.
0629      */
0630     void normalize(double intg) {
0631       double oldintg = sumAllBinHeights();
0632       if ( oldintg == 0.0 ) return;
0633       for ( int ix = 0; ix < xax->bins() + 2; ++ix )
0634     for ( int iy = 0; iy < yax->bins() + 2; ++iy ) {
0635       double fac = intg/oldintg;
0636       if ( ix >= 2 && iy >= 2 )
0637         fac /= (xax->binUpperEdge(ix - 2) - xax->binLowerEdge(ix - 2))*
0638           (yax->binUpperEdge(iy - 2) - yax->binLowerEdge(iy - 2));
0639         sumw[ix][iy] *= fac;
0640         sumxw[ix][iy] *= fac;
0641         sumx2w[ix][iy] *= fac;
0642         sumyw[ix][iy] *= fac;
0643         sumy2w[ix][iy] *= fac;
0644         sumw2[ix][iy] *= fac*fac;
0645       }
0646     }
0647
0648     /**
0649      * Return the integral over the histogram bins assuming it has been
0650      * normalize()d.
0651      */
0652     // double integral() const {
0653     //   double intg = sumw[0] + sumw[1];
0654     //   for ( int i = 2; i < ax->bins() + 2; ++i )
0655
0656     // is this right? Leave out bin width factor?
0657
0658     //     intg += sumw[ix][iy]*(ax->binUpperEdge(i - 2) - ax->binLowerEdge(i - 2));
0659     //   return intg;
0660     // }
0661
0662     /**
0663      * Not implemented in LWH.
0664      * @return null pointer always.
0665      */
0666     void * cast(const std::string &) const {
0667       return 0;
0668     }
0669
0670     /**
0671      * Write out the histogram in the AIDA xml format.
0672      */
0673     bool writeXML(std::ostream & os, std::string path, std::string name) {
0674       //std::cout << "Writing out histogram " << name << " in AIDA file format!" << std::endl;
0675       os << "  <histogram2d name=\"" << name
0676          << "\"\n    title=\"" << title()
0677          << "\" path=\"" << path
0678          << "\">\n    <axis max=\"" << xax->upperEdge()
0679          << "\" numberOfBins=\"" << xax->bins()
0680          << "\" min=\"" << xax->lowerEdge()
0681          << "\" direction=\"x\"";
0682       if ( xvax ) {
0683         os << ">\n";
0684         for ( int i = 0, N = xax->bins() - 1; i < N; ++i )
0685           os << "      <binBorder value=\"" << xax->binUpperEdge(i) << "\"/>\n";
0686         os << "    </axis>\n";
0687       } else {
0688         os << "/>\n";
0689       }
0690       os << "    <axis max=\"" << yax->upperEdge()
0691          << "\" numberOfBins=\"" << yax->bins()
0692          << "\" min=\"" << yax->lowerEdge()
0693          << "\" direction=\"y\"";
0694       if ( yvax ) {
0695         os << ">\n";
0696         for ( int i = 0, N = yax->bins() - 1; i < N; ++i )
0697           os << "      <binBorder value=\"" << yax->binUpperEdge(i) << "\"/>\n";
0698         os << "    </axis>\n";
0699       } else {
0700         os << "/>\n";
0701       }
0702       os << "    <statistics entries=\"" << entries()
0703          << "\">\n      <statistic mean=\"" << meanX()
0704          << "\" direction=\"x\"\n        rms=\"" << rmsX()
0705          << "\"/>\n    </statistics>\n      <statistic mean=\"" << meanY()
0706          << "\" direction=\"y\"\n        rms=\"" << rmsY()
0707          << "\"/>\n    </statistics>\n    <data2d>\n";
0708       for ( int ix = 0; ix < xax->bins() + 2; ++ix )
0709     for ( int iy = 0; iy < yax->bins() + 2; ++iy )
0710       if ( sum[ix][iy] ) {
0711         os << "      <bin2d binNumX=\"";
0712         if ( ix == 0 ) os << "UNDERFLOW";
0713         else if ( ix == 1 ) os << "OVERFLOW";
0714         else os << ix - 2;
0715         os << "\" binNumY=\"";
0716         if ( iy == 0 ) os << "UNDERFLOW";
0717         else if ( iy == 1 ) os << "OVERFLOW";
0718         else os << iy - 2;
0719         os << "\" entries=\"" << sum[ix][iy]
0720            << "\" height=\"" << sumw[ix][iy]
0721            << "\"\n        error=\"" << std::sqrt(sumw2[ix][iy])
0722            << "\" error2=\"" << sumw2[ix][iy]
0723            << "\"\n        weightedMeanX=\"" << binMeanX(ix - 2, iy - 2)
0724            << "\" weightedRmsX=\"" << binRmsX(ix - 2, iy - 2)
0725            << "\"\n        weightedMeanY=\"" << binMeanY(ix - 2, iy - 2)
0726            << "\" weightedRmsY=\"" << binRmsY(ix - 2, iy - 2)
0727            << "\"/>\n";
0728         }
0729       os << "    </data2d>\n  </histogram2d>" << std::endl;
0730       return true;
0731     }
0732
0733
0734     /**
0735      * Write out the histogram in a flat text file suitable for
0736      * eg. gnuplot to read. The coloums are layed out as 'x w w2 n'.
0737      */
0738     bool writeFLAT(std::ostream & os, std::string path, std::string name) {
0739       os << "#2D " << path << "/" << name << " " << xax->lowerEdge()
0740          << " " << xax->bins() << " " << xax->upperEdge() << " "
0741      << yax->lowerEdge() << " " << yax->bins() << " " << yax->upperEdge()
0742          << " \"" << title() << "\"" << std::endl;
0743       for ( int ix = 2; ix < xax->bins() + 2; ++ix ) {
0744     for ( int iy = 2; iy < yax->bins() + 2; ++iy )
0745       os << 0.5*(xax->binLowerEdge(ix - 2)+xax->binUpperEdge(ix - 2)) << " "
0746          << 0.5*(yax->binLowerEdge(iy - 2)+yax->binUpperEdge(iy - 2))
0747          << " " << sumw[ix][iy] << " " << sqrt(sumw2[ix][iy])
0748          << " " << sum[ix][iy] << std::endl;
0749     os << std::endl;
0750       }
0751       os << std::endl;
0752       return true;
0753     }
0754
0755
0756
0757    #ifdef HAVE_ROOT
0758     /**
0759      * Write out the histogram in Root file format.
0760      */
0761     //bool writeROOT(std::ostream & os, std::string path, std::string name) {
0762     bool writeROOT(TFile* file, std::string path, std::string name) {
0763
0764       //std::cout << "Writing out histogram " << name.c_str() << " in ROOT file format" << std::endl;
0765
0766       TH1D* hist1d;
0767       int nbins;
0768       if (!vax || vax->isFixedBinning() ) {//equidistant binning (easier case)
0769         nbins = ax->bins();
0770         hist1d = new TH1D(name.c_str(), title().c_str(), nbins, ax->lowerEdge(), ax->upperEdge());
0771       }
0772       else {
0773         nbins = vax->bins();
0774         double* bins = new double[nbins+1];
0775         for (int i=0; i<nbins; ++i) {
0776       bins[ix][iy] = vax->binEdges(i).first;
0777         }
0778         bins[nbins] = vax->binEdges(nbins-1).second; //take last bin right border
0779         hist1d = new TH1D(name.c_str(), title().c_str(), nbins, bins);
0780         delete [] bins;
0781       }
0782
0783
0784       double entries = 0;
0785       for ( int i = 0; i < nbins + 2; ++i ) {
0786         if ( sum[ix][iy] ) {
0787           //i==0: underflow->RootBin(0), i==1: overflow->RootBin(NBins+1)
0788           entries = entries + sum[ix][iy];
0789           int j=i;
0790           if (i==0) j=0; //underflow
0791           else if (i==1) j=nbins+1; //overflow
0792           if (i>=2) j=i-1; //normal bin entries
0793           hist1d->SetBinContent(j, sumw[ix][iy]);
0794           hist1d->SetBinError(j, sqrt(sumw2[ix][iy]));
0795           //hist1d->Fill(binMean(i), sumw[ix][iy]);
0796         }
0797       }
0798
0799       hist1d->Sumw2();
0800       hist1d->SetEntries(entries);
0801
0802       std::string DirName; //remove preceding slash from directory name, else ROOT error
0803       for (unsigned int i=1; i<path.size(); ++i) DirName += path[i];
0804       if (!file->Get(DirName.c_str())) file->mkdir(DirName.c_str());
0805       file->cd(DirName.c_str());
0806       hist1d->Write();
0807
0808       delete hist1d;
0809
0810       return true;
0811     }
0812
0813    #endif
0814
0815
0816
0817   private:
0818
0819     /** The title */
0820     std::string theTitle;
0821
0822     /** The axis. */
0823     IAxis * xax;
0824
0825     /** Pointer (possibly null) to a axis with fixed bin width. */
0826     Axis * xfax;
0827
0828     /** Pointer (possibly null) to a axis with fixed bin width. */
0829     VariAxis * xvax;
0830
0831     /** The axis. */
0832     IAxis * yax;
0833
0834     /** Pointer (possibly null) to a axis with fixed bin width. */
0835     Axis * yfax;
0836
0837     /** Pointer (possibly null) to a axis with fixed bin width. */
0838     VariAxis * yvax;
0839
0840     /** The counts. */
0841     std::vector< std::vector<int> > sum;
0842
0843     /** The weights. */
0844     std::vector< std::vector<double> > sumw;
0845
0846     /** The squared weights. */
0847     std::vector< std::vector<double> > sumw2;
0848
0849     /** The weighted x-values. */
0850     std::vector< std::vector<double> > sumxw;
0851
0852     /** The weighted x-square-values. */
0853     std::vector< std::vector<double> > sumx2w;
0854
0855     /** The weighted y-values. */
0856     std::vector< std::vector<double> > sumyw;
0857
0858     /** The weighted y-square-values. */
0859     std::vector< std::vector<double> > sumy2w;
0860
0861     /** dummy pointer to non-existen annotation. */
0862     IAnnotation * anno;
0863
0864   };
0865
0866 }
0867
0868 #endif /* LWH_Histogram1D_H */