Histogram2D.h

// -*- C++ -*-
#ifndef LWH_Histogram2D_H
#define LWH_Histogram2D_H
//
// This is the declaration of the Histogram1D class.
//

#include "AIHistogram2D.h"
#include "ManagedObject.h"
#include "Axis.h"
#include "VariAxis.h"
#include <vector>
#include <stdexcept>

#include <iostream>

namespace LWH {

  using namespace AIDA;


  class Histogram2D: public IHistogram2D, public ManagedObject {
    
  public:

    friend class HistogramFactory;

  public:

    Histogram2D(int nx, double lox, double upx,
                int ny, double loy, double upy)
      : xfax(new Axis(nx, lox, upx)), xvax(0), yfax(new Axis(ny, loy, upy)),
        sum(nx + 2, std::vector<int>(ny + 2)),
        sumw(nx + 2, std::vector<double>(ny + 2)),
        sumw2(nx + 2, std::vector<double>(ny + 2)),
        sumxw(nx + 2, std::vector<double>(ny + 2)),
        sumx2w(nx + 2, std::vector<double>(ny + 2)),
        sumyw(nx + 2, std::vector<double>(ny + 2)),
        sumy2w(nx + 2, std::vector<double>(ny + 2)) {
      xax = xfax;
      yax = yfax;
    }

    Histogram2D(const std::vector<double> & xedges,
                const std::vector<double> & yedges)
      : xfax(0), xvax(new VariAxis(xedges)),
        yfax(0), yvax(new VariAxis(xedges)),
        sum(xedges.size() + 1, std::vector<int>(yedges.size() + 1)),
        sumw(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
        sumw2(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
        sumxw(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
        sumx2w(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
        sumyw(xedges.size() + 1, std::vector<double>(yedges.size() + 1)),
        sumy2w(xedges.size() + 1, std::vector<double>(yedges.size() + 1)) {
      xax = xvax;
      yax = yvax;
    }

    Histogram2D(const Histogram2D & h)
      : IBaseHistogram(h), IHistogram(h), IHistogram2D(h), ManagedObject(h),
        xfax(0), xvax(0),  yfax(0), yvax(0),
        sum(h.sum), sumw(h.sumw), sumw2(h.sumw2),
        sumxw(h.sumxw), sumx2w(h.sumx2w) ,
        sumyw(h.sumyw), sumy2w(h.sumy2w){
      const VariAxis * hxvax = dynamic_cast<const VariAxis *>(h.xax);
      if ( hxvax ) xax = xvax = new VariAxis(*hxvax);
      else xax = xfax = new Axis(dynamic_cast<const Axis &>(*h.xax));
      const VariAxis * hyvax = dynamic_cast<const VariAxis *>(h.yax);
      if ( hyvax ) yax = yvax = new VariAxis(*hyvax);
      else yax = yfax = new Axis(dynamic_cast<const Axis &>(*h.yax));
  }

    virtual ~Histogram2D() {
      delete xax;
      delete yax;
    }

    std::string title() const {
      return theTitle;
    }

    std::string name() const {
      return title();
    }

    bool setTitle(const std::string & title) {
      theTitle = title;
      return true;
    }

    IAnnotation & annotation() {
      throw std::runtime_error("LWH cannot handle annotations");
      return *anno;
    }

    const IAnnotation & annotation() const {
      throw std::runtime_error("LWH cannot handle annotations");
      return *anno;
    }

    int dimension() const {
      return 2;
    }

    bool reset() {
      const int nx = xax->bins() + 2;
      const int ny = yax->bins() + 2;
      sum = std::vector< std::vector<int> >(nx, std::vector<int>(ny));
      sumw = std::vector< std::vector<double> >(nx, std::vector<double>(ny));
      sumw2 = sumw;
      sumxw = sumw;
      sumx2w = sumw;
      sumyw = sumw;
      sumy2w = sumw;
      return true;
    }

    int entries() const {
      int si = 0;
      for ( int ix = 2; ix < xax->bins() + 2; ++ix )
        for ( int iy = 2; iy < yax->bins() + 2; ++iy ) si += sum[ix][iy];
      return si;
    }

    int allEntries() const {
      return entries() + extraEntries();
    }

    int extraEntries() const {
      int esum = sum[0][0] + sum[1][0] + sum[0][1] + sum[1][1];
      for ( int ix = 2; ix < xax->bins() + 2; ++ix )
        esum += sum[ix][0] + sum[ix][1];
      for ( int iy = 2; iy < yax->bins() + 2; ++iy )
        esum += sum[0][iy] + sum[1][iy];
      return esum;
    }

    double equivalentBinEntries() const {
      double sw = 0.0;
      double sw2 = 0.0;
      for ( int ix = 2; ix < xax->bins() + 2; ++ix )
        for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
          sw += sumw[ix][iy];
          sw2 += sumw2[ix][iy];
        }
      return sw2/(sw*sw);
    }

    double sumBinHeights() const {
      double sw = 0.0;
      for ( int ix = 2; ix < xax->bins() + 2; ++ix )
        for ( int iy = 2; iy < yax->bins() + 2; ++iy ) sw += sumw[ix][iy];
      return sw;
    }

    double sumAllBinHeights() const {
      return sumBinHeights() + sumExtraBinHeights();
    }

    double sumExtraBinHeights() const {
      int esum = sumw[0][0] + sumw[1][0] + sumw[0][1] + sumw[1][1];
      for ( int ix = 2; ix < xax->bins() + 2; ++ix )
        esum += sumw[ix][0] + sumw[ix][1];
      for ( int iy = 2; iy < yax->bins() + 2; ++iy )
        esum += sumw[0][iy] + sumw[1][iy];
      return esum;
    }

    double minBinHeight() const {
      double minw = sumw[2][2];
      for ( int ix = 2; ix < xax->bins() + 2; ++ix )
        for ( int iy = 2; iy < yax->bins() + 2; ++iy )
          minw = std::min(minw, sumw[ix][iy]);
      return minw;
    }

    double maxBinHeight() const{
      double maxw = sumw[2][2];
      for ( int ix = 2; ix < xax->bins() + 2; ++ix )
        for ( int iy = 2; iy < yax->bins() + 2; ++iy )
          maxw = std::max(maxw, sumw[ix][iy]);
      return maxw;
    }

    bool fill(double x, double y, double weight = 1.) {
      int ix = xax->coordToIndex(x) + 2;
      int iy = yax->coordToIndex(y) + 2;
      ++sum[ix][iy];
      sumw[ix][iy] += weight;
      sumxw[ix][iy] += x*weight;
      sumx2w[ix][iy] += x*x*weight;
      sumyw[ix][iy] += y*weight;
      sumy2w[ix][iy] += y*y*weight;
      sumw2[ix][iy] += weight*weight;
      return weight >= 0 && weight <= 1;
    }

    double binMeanX(int xindex, int yindex) const {
      int ix = xindex + 2;
      int iy = yindex + 2;
      return sumw[ix][iy] != 0.0? sumxw[ix][iy]/sumw[ix][iy]:
        ( xvax? xvax->binMidPoint(xindex): xfax->binMidPoint(xindex) );
    };

    double binMeanY(int xindex, int yindex) const {
      int ix = xindex + 2;
      int iy = yindex + 2;
      return sumw[ix][iy] != 0.0? sumyw[ix][iy]/sumw[ix][iy]:
        ( yvax? yvax->binMidPoint(yindex): xfax->binMidPoint(yindex) );
    };

    double binRmsX(int xindex, int yindex) const {
      int ix = xindex + 2;
      int iy = yindex + 2;
      return sumw[ix][iy] == 0.0 || sum[ix][iy] < 2? xax->binWidth(xindex):
        std::sqrt(std::max(sumw[ix][iy]*sumx2w[ix][iy] -
                           sumxw[ix][iy]*sumxw[ix][iy], 0.0))/sumw[ix][iy];
    };

    double binRmsY(int xindex, int yindex) const {
      int ix = xindex + 2;
      int iy = yindex + 2;
      return sumw[ix][iy] == 0.0 || sum[ix][iy] < 2? yax->binWidth(yindex):
        std::sqrt(std::max(sumw[ix][iy]*sumy2w[ix][iy] -
                           sumyw[ix][iy]*sumyw[ix][iy], 0.0))/sumw[ix][iy];
    };

    int binEntries(int xindex, int yindex) const {
      return sum[xindex + 2][yindex + 2];
    }

    virtual int binEntriesX(int index) const {
      int ret = 0;
      for ( int iy = 2; iy < yax->bins() + 2; ++iy )
        ret += sum[index + 2][iy];
      return ret;
    }

    virtual int binEntriesY(int index) const {
      int ret = 0;
      for ( int ix = 2; ix < xax->bins() + 2; ++ix )
        ret += sum[ix][index + 2];
      return ret;
    }

    double binHeight(int xindex, int yindex) const {
      return sumw[xindex + 2][yindex + 2];
    }

    virtual double binHeightX(int index) const {
      double ret = 0;
      for ( int iy = 2; iy < yax->bins() + 2; ++iy )
        ret += sumw[index + 2][iy];
      return ret;
    }

    virtual double binHeightY(int index) const {
      double ret = 0;
      for ( int ix = 2; ix < xax->bins() + 2; ++ix )
        ret += sumw[ix][index + 2];
      return ret;
    }

    double binError(int xindex, int yindex) const {
      return std::sqrt(sumw2[xindex + 2][yindex + 2]);
    }

    double meanX() const {
      double s = 0.0;
      double sx = 0.0;
      for ( int ix = 2; ix < xax->bins() + 2; ++ix )
        for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
        s += sumw[ix][iy];
        sx += sumxw[ix][iy];
      }
      return s != 0.0? sx/s: 0.0;
    }

    double meanY() const {
      double s = 0.0;
      double sy = 0.0;
      for ( int ix = 2; ix < xax->bins() + 2; ++ix )
        for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
        s += sumw[ix][iy];
        sy += sumyw[ix][iy];
      }
      return s != 0.0? sy/s: 0.0;
    }

    double rmsX() const {
      double s = 0.0;
      double sx = 0.0;
      double sx2 = 0.0;
      for ( int ix = 2; ix < xax->bins() + 2; ++ix )
        for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
        s += sumw[ix][iy];
        sx += sumxw[ix][iy];
        sx2 += sumx2w[ix][iy];
      }
      return s != 0.0? std::sqrt(std::max(s*sx2 - sx*sx, 0.0))/s:
        xax->upperEdge() - xax->lowerEdge();
    }

    double rmsY() const {
      double s = 0.0;
      double sy = 0.0;
      double sy2 = 0.0;
      for ( int ix = 2; ix < xax->bins() + 2; ++ix )
        for ( int iy = 2; iy < yax->bins() + 2; ++iy ) {
        s += sumw[ix][iy];
        sy += sumyw[ix][iy];
        sy2 += sumy2w[ix][iy];
      }
      return s != 0.0? std::sqrt(std::max(s*sy2 - sy*sy, 0.0))/s:
        yax->upperEdge() - yax->lowerEdge();
    }

    double getSumW(int xindex, int yindex) const {
        return sumw[xindex + 2][yindex + 2];
    }

    double getSumW2(int xindex, int yindex) const {
        return sumw2[xindex + 2][yindex + 2];
    }

    double getSumXW(int xindex, int yindex) const {
        return sumxw[xindex + 2][yindex + 2];
    }

    double getSumX2W(int xindex, int yindex) const {
        return sumx2w[xindex + 2][yindex + 2];
    }
    
    double getSumYW(int xindex, int yindex) const {
        return sumyw[xindex + 2][yindex + 2];
    }

    double getSumY2W(int xindex, int yindex) const {
        return sumy2w[xindex + 2][yindex + 2];
    }
    
    const IAxis & xAxis() const {
      return *xax;
    }

    const IAxis & yAxis() const {
      return *yax;
    }

    int coordToIndexX(double coord) const {
      return xax->coordToIndex(coord);
    }

    int coordToIndexY(double coord) const {
      return yax->coordToIndex(coord);
    }

    bool add(const Histogram2D & h) {
      if ( xax->upperEdge() != h.xax->upperEdge() ||
           xax->lowerEdge() != h.xax->lowerEdge() ||
           xax->bins() != h.xax->bins() ) return false;
      if ( yax->upperEdge() != h.yax->upperEdge() ||
           yax->lowerEdge() != h.yax->lowerEdge() ||
           yax->bins() != h.yax->bins() ) return false;
      for ( int ix = 0; ix < xax->bins() + 2; ++ix )
        for ( int iy = 0; iy < yax->bins() + 2; ++iy ) {
          sum[ix][iy] += h.sum[ix][iy];
          sumw[ix][iy] += h.sumw[ix][iy];
          sumxw[ix][iy] += h.sumxw[ix][iy];
          sumx2w[ix][iy] += h.sumx2w[ix][iy];
          sumyw[ix][iy] += h.sumyw[ix][iy];
          sumy2w[ix][iy] += h.sumy2w[ix][iy];
          sumw2[ix][iy] += h.sumw2[ix][iy];
        }
      return true;
    }

    bool add(const IHistogram2D & hist) {
      return add(dynamic_cast<const Histogram2D &>(hist));
    }

    bool scale(double s) {
      for ( int ix = 0; ix < xax->bins() + 2; ++ix )
        for ( int iy = 0; iy < yax->bins() + 2; ++iy ) {
          sumw[ix][iy] *= s;
          sumxw[ix][iy] *= s;
          sumx2w[ix][iy] *= s;
          sumyw[ix][iy] *= s;
          sumy2w[ix][iy] *= s;
          sumw2[ix][iy] *= s*s;
      }
      return true;
    }

    void normalize(double intg) {
      double oldintg = sumAllBinHeights();
      if ( oldintg == 0.0 ) return;
      for ( int ix = 0; ix < xax->bins() + 2; ++ix )
        for ( int iy = 0; iy < yax->bins() + 2; ++iy ) {
          double fac = intg/oldintg;
          if ( ix >= 2 && iy >= 2 )
            fac /= (xax->binUpperEdge(ix - 2) - xax->binLowerEdge(ix - 2))*
              (yax->binUpperEdge(iy - 2) - yax->binLowerEdge(iy - 2));
        sumw[ix][iy] *= fac;
        sumxw[ix][iy] *= fac;
        sumx2w[ix][iy] *= fac;
        sumyw[ix][iy] *= fac;
        sumy2w[ix][iy] *= fac;
        sumw2[ix][iy] *= fac*fac;
      }
    }

    // double integral() const {
    //   double intg = sumw[0] + sumw[1];
    //   for ( int i = 2; i < ax->bins() + 2; ++i )

    // is this right? Leave out bin width factor?

    //     intg += sumw[ix][iy]*(ax->binUpperEdge(i - 2) - ax->binLowerEdge(i - 2));
    //   return intg;
    // }

    void * cast(const std::string &) const {
      return 0;
    }

    bool writeXML(std::ostream & os, std::string path, std::string name) {
      //std::cout << "Writing out histogram " << name << " in AIDA file format!" << std::endl;
      os << "  <histogram2d name=\"" << name
         << "\"\n    title=\"" << title()
         << "\" path=\"" << path
         << "\">\n    <axis max=\"" << xax->upperEdge()
         << "\" numberOfBins=\"" << xax->bins()
         << "\" min=\"" << xax->lowerEdge()
         << "\" direction=\"x\"";
      if ( xvax ) {
        os << ">\n";
        for ( int i = 0, N = xax->bins() - 1; i < N; ++i )
          os << "      <binBorder value=\"" << xax->binUpperEdge(i) << "\"/>\n";
        os << "    </axis>\n";
      } else {
        os << "/>\n";
      }
      os << "    <axis max=\"" << yax->upperEdge()
         << "\" numberOfBins=\"" << yax->bins()
         << "\" min=\"" << yax->lowerEdge()
         << "\" direction=\"y\"";
      if ( yvax ) {
        os << ">\n";
        for ( int i = 0, N = yax->bins() - 1; i < N; ++i )
          os << "      <binBorder value=\"" << yax->binUpperEdge(i) << "\"/>\n";
        os << "    </axis>\n";
      } else {
        os << "/>\n";
      }
      os << "    <statistics entries=\"" << entries()
         << "\">\n      <statistic mean=\"" << meanX()
         << "\" direction=\"x\"\n        rms=\"" << rmsX()
         << "\"/>\n    </statistics>\n      <statistic mean=\"" << meanY()
         << "\" direction=\"y\"\n        rms=\"" << rmsY()
         << "\"/>\n    </statistics>\n    <data2d>\n";
      for ( int ix = 0; ix < xax->bins() + 2; ++ix )
        for ( int iy = 0; iy < yax->bins() + 2; ++iy )
          if ( sum[ix][iy] ) {
            os << "      <bin2d binNumX=\"";
            if ( ix == 0 ) os << "UNDERFLOW";
            else if ( ix == 1 ) os << "OVERFLOW";
            else os << ix - 2;
            os << "\" binNumY=\"";
            if ( iy == 0 ) os << "UNDERFLOW";
            else if ( iy == 1 ) os << "OVERFLOW";
            else os << iy - 2;
            os << "\" entries=\"" << sum[ix][iy]
               << "\" height=\"" << sumw[ix][iy]
               << "\"\n        error=\"" << std::sqrt(sumw2[ix][iy])
               << "\" error2=\"" << sumw2[ix][iy]
               << "\"\n        weightedMeanX=\"" << binMeanX(ix - 2, iy - 2)
               << "\" weightedRmsX=\"" << binRmsX(ix - 2, iy - 2)
               << "\"\n        weightedMeanY=\"" << binMeanY(ix - 2, iy - 2)
               << "\" weightedRmsY=\"" << binRmsY(ix - 2, iy - 2)
               << "\"/>\n";
        }
      os << "    </data2d>\n  </histogram2d>" << std::endl;
      return true;
    }


    bool writeFLAT(std::ostream & os, std::string path, std::string name) {
      os << "#2D " << path << "/" << name << " " << xax->lowerEdge()
         << " " << xax->bins() << " " << xax->upperEdge() << " "
         << yax->lowerEdge() << " " << yax->bins() << " " << yax->upperEdge()
         << " \"" << title() << "\"" << std::endl;
      for ( int ix = 2; ix < xax->bins() + 2; ++ix ) {
        for ( int iy = 2; iy < yax->bins() + 2; ++iy )
          os << 0.5*(xax->binLowerEdge(ix - 2)+xax->binUpperEdge(ix - 2)) << " "
             << 0.5*(yax->binLowerEdge(iy - 2)+yax->binUpperEdge(iy - 2))
             << " " << sumw[ix][iy] << " " << sqrt(sumw2[ix][iy])
             << " " << sum[ix][iy] << std::endl;
        os << std::endl;
      }
      os << std::endl;
      return true;
    }



   #ifdef HAVE_ROOT

    //bool writeROOT(std::ostream & os, std::string path, std::string name) {
    bool writeROOT(TFile* file, std::string path, std::string name) {

      //std::cout << "Writing out histogram " << name.c_str() << " in ROOT file format" << std::endl;

      TH1D* hist1d;
      int nbins;
      if (!vax || vax->isFixedBinning() ) {//equidistant binning (easier case)
        nbins = ax->bins();
        hist1d = new TH1D(name.c_str(), title().c_str(), nbins, ax->lowerEdge(), ax->upperEdge());
      }
      else {
        nbins = vax->bins();
        double* bins = new double[nbins+1];
        for (int i=0; i<nbins; ++i) {
      bins[ix][iy] = vax->binEdges(i).first;
        }
        bins[nbins] = vax->binEdges(nbins-1).second; //take last bin right border
        hist1d = new TH1D(name.c_str(), title().c_str(), nbins, bins);
        delete [] bins;
      }


      double entries = 0;
      for ( int i = 0; i < nbins + 2; ++i ) {
        if ( sum[ix][iy] ) {
          //i==0: underflow->RootBin(0), i==1: overflow->RootBin(NBins+1)
          entries = entries + sum[ix][iy];
          int j=i;
          if (i==0) j=0; //underflow
          else if (i==1) j=nbins+1; //overflow
          if (i>=2) j=i-1; //normal bin entries
          hist1d->SetBinContent(j, sumw[ix][iy]);
          hist1d->SetBinError(j, sqrt(sumw2[ix][iy]));
          //hist1d->Fill(binMean(i), sumw[ix][iy]);
        }
      }

      hist1d->Sumw2();
      hist1d->SetEntries(entries);

      std::string DirName; //remove preceding slash from directory name, else ROOT error
      for (unsigned int i=1; i<path.size(); ++i) DirName += path[i];
      if (!file->Get(DirName.c_str())) file->mkdir(DirName.c_str());
      file->cd(DirName.c_str());
      hist1d->Write();

      delete hist1d;

      return true;
    }

   #endif



  private:

    std::string theTitle;

    IAxis * xax;

    Axis * xfax;

    VariAxis * xvax;

    IAxis * yax;

    Axis * yfax;

    VariAxis * yvax;

    std::vector< std::vector<int> > sum;

    std::vector< std::vector<double> > sumw;

    std::vector< std::vector<double> > sumw2;

    std::vector< std::vector<double> > sumxw;

    std::vector< std::vector<double> > sumx2w;

    std::vector< std::vector<double> > sumyw;

    std::vector< std::vector<double> > sumy2w;

    IAnnotation * anno;

  };

}

#endif /* LWH_Histogram1D_H */