Histogram1D.h

// -*- C++ -*-
//
// Histogram1D.h is a part of ThePEG - Toolkit for HEP Event Generation
// Copyright (C) 1999-2019 Leif Lonnblad
//
// ThePEG is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef LWH_Histogram1D_H
#define LWH_Histogram1D_H
//
// This is the declaration of the Histogram1D class.
//

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

namespace LWH {

using namespace AIDA;

class Histogram1D: public IHistogram1D, public ManagedObject {

public:

  friend class HistogramFactory;

public:

  Histogram1D(int n, double lo, double up)
    : fax(new Axis(n, lo, up)), vax(0),
      sum(n + 2), sumw(n + 2), sumw2(n + 2), sumxw(n + 2), sumx2w(n + 2) {
    ax = fax;
  }

  Histogram1D(const std::vector<double> & edges)
    : fax(0), vax(new VariAxis(edges)),
      sum(edges.size() + 1), sumw(edges.size() + 1), sumw2(edges.size() + 1),
      sumxw(edges.size() + 1), sumx2w(edges.size() + 1) {
    ax = vax;
  }

  Histogram1D(const Histogram1D & h)
    : IBaseHistogram(h), IHistogram(h), IHistogram1D(h), ManagedObject(h),
      fax(0), vax(0), sum(h.sum), sumw(h.sumw), sumw2(h.sumw2),
      sumxw(h.sumxw), sumx2w(h.sumx2w) {
    const VariAxis * hvax = dynamic_cast<const VariAxis *>(h.ax);
    if ( vax ) ax = vax = new VariAxis(*hvax);
    else ax = fax = new Axis(dynamic_cast<const Axis &>(*h.ax));
}

  virtual ~Histogram1D() {
    delete ax;
  }

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

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

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

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

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

  int dimension() const {
    return 1;
  }

  bool reset() {
    sum = std::vector<int>(ax->bins() + 2);
    sumw = std::vector<double>(ax->bins() + 2);
    sumxw = std::vector<double>(ax->bins() + 2);
    sumx2w = std::vector<double>(ax->bins() + 2);
    sumw2 = std::vector<double>(ax->bins() + 2);
    return true;
  }

  int entries() const {
    int si = 0;
    for ( int i = 2; i < ax->bins() + 2; ++i ) si += sum[i];
    return si;
  }

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

  int extraEntries() const {
    return sum[0] + sum[1];
  }

  double equivalentBinEntries() const {
    double sw = 0.0;
    double sw2 = 0.0;
    for ( int i = 2; i < ax->bins() + 2; ++i ) {
      sw += sumw[i];
      sw2 += sumw2[i];
    }
    return sw2/(sw*sw);
  }
    
  double sumBinHeights() const {
    double sw = 0.0;
    for ( int i = 2; i < ax->bins() + 2; ++i ) sw += sumw[i];
    return sw;
  }
    
  double sumAllBinHeights() const {
    return sumBinHeights() + sumExtraBinHeights();
  }

  double sumExtraBinHeights() const {
    return sumw[0] + sumw[1];
  }

  double minBinHeight() const {
    double minw = sumw[2];
    for ( int i = 3; i < ax->bins() + 2; ++i ) minw = std::min(minw, sumw[i]);
    return minw;
  }

  double maxBinHeight() const{
    double maxw = sumw[2];
    for ( int i = 3; i < ax->bins() + 2; ++i ) maxw = std::max(maxw, sumw[i]);
    return maxw;
  }

  bool fill(double x, double weight = 1.) {
    int i = ax->coordToIndex(x) + 2;
    ++sum[i];
    sumw[i] += weight;
    sumxw[i] += x*weight;
    sumx2w[i] += x*x*weight;
    sumw2[i] += weight*weight;
    return weight >= 0 && weight <= 1;
  }

  double binMean(int index) const {
    int i = index + 2;
    return sumw[i] != 0.0? sumxw[i]/sumw[i]:
      ( vax? vax->binMidPoint(index): fax->binMidPoint(index) );
  };

  double binRms(int index) const {
    int i = index + 2;
    return sumw[i] == 0.0 || sum[i] < 2? ax->binWidth(index):
      std::sqrt(std::max(sumw[i]*sumx2w[i] - sumxw[i]*sumxw[i], 0.0))/sumw[i];
  };

  int binEntries(int index) const {
    return sum[index + 2];
  }

  double binHeight(int index) const {
    return sumw[index + 2];
  }

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

  double mean() const {
    double s = 0.0;
    double sx = 0.0;
    for ( int i = 2; i < ax->bins() + 2; ++i ) {
      s += sumw[i];
      sx += sumxw[i];
    }
    return s != 0.0? sx/s: 0.0;
  }

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

  const IAxis & axis() const {
    return *ax;
  }

  int coordToIndex(double coord) const {
    return ax->coordToIndex(coord);
  }
  
  bool add(const Histogram1D & h) {
    if ( ax->upperEdge() != h.ax->upperEdge() ||
         ax->lowerEdge() != h.ax->lowerEdge() ||
         ax->bins() != h.ax->bins() ) return false;
    for ( int i = 0; i < ax->bins() + 2; ++i ) {
      sum[i] += h.sum[i];
      sumw[i] += h.sumw[i];
      sumxw[i] += h.sumxw[i];
      sumx2w[i] += h.sumx2w[i];
      sumw2[i] += h.sumw2[i];
    }
    return true;
  }

  bool add(const IHistogram1D & hist) {
    return add(dynamic_cast<const Histogram1D &>(hist));
  }

  bool scale(double s) {
    for ( int i = 0; i < ax->bins() + 2; ++i ) {
      sumw[i] *= s;
      sumxw[i] *= s;
      sumx2w[i] *= s;
      sumw2[i] *= s*s;
    }
    return true;
  }

  void normalize(double intg) {
    double oldintg = sumAllBinHeights();
    if ( oldintg == 0.0 ) return;
    for ( int i = 0; i < ax->bins() + 2; ++i ) {
      double fac = intg/oldintg;
      if ( i >= 2 ) fac /= (ax->binUpperEdge(i - 2) - ax->binLowerEdge(i - 2));
      sumw[i] *= fac;
      sumxw[i] *= fac;
      sumx2w[i] *= fac;
      sumw2[i] *= fac*fac;
    }
  }

  double integral() const {
    double intg = sumw[0] + sumw[1];
    for ( int i = 2; i < ax->bins() + 2; ++i )
      intg += sumw[i]*(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) {
    os << "  <histogram1d name=\"" << name
       << "\"\n    title=\"" << title()
       << "\" path=\"" << path
       << "\">\n    <axis max=\"" << ax->upperEdge()
       << "\" numberOfBins=\"" << ax->bins()
       << "\" min=\"" << ax->lowerEdge()
       << "\" direction=\"x\"";
    if ( vax ) {
      os << ">\n";
      for ( int i = 0, N = ax->bins() - 1; i < N; ++i )
        os << "      <binBorder value=\"" << ax->binUpperEdge(i) << "\"/>\n";
      os << "    </axis>\n";
    } else {
      os << "/>\n";
    }
    os << "    <statistics entries=\"" << entries()
       << "\">\n      <statistic mean=\"" << mean()
       << "\" direction=\"x\"\n        rms=\"" << rms()
       << "\"/>\n    </statistics>\n    <data1d>\n";
    for ( int i = 0; i < ax->bins() + 2; ++i ) if ( sum[i] ) {
      os << "      <bin1d binNum=\"";
      if ( i == 0 ) os << "UNDERFLOW";
      else if ( i == 1 ) os << "OVERFLOW";
      else os << i - 2;
      os << "\" entries=\"" << sum[i]
         << "\" height=\"" << sumw[i]
         << "\"\n        error=\"" << std::sqrt(sumw2[i])
         << "\" error2=\"" << sumw2[i]
         << "\"\n        weightedMean=\"" << binMean(i - 2)
         << "\" weightedRms=\"" << binRms(i - 2)
         << "\"/>\n";
    }
    os << "    </data1d>\n  </histogram1d>" << std::endl;
    return true;
  }


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

private:

  std::string theTitle;

  IAxis * ax;

  Axis * fax;

  VariAxis * vax;

  std::vector<int> sum;

  std::vector<double> sumw;

  std::vector<double> sumw2;

  std::vector<double> sumxw;

  std::vector<double> sumx2w;

};

}

#endif /* LWH_Histogram1D_H */