Histogram.h

// -*- C++ -*-
//
// Histogram.h is a part of Herwig - A multi-purpose Monte Carlo event generator
// Copyright (C) 2002-2019 The Herwig Collaboration
//
// Herwig is licenced under version 3 of the GPL, see COPYING for details.
// Please respect the MCnet academic guidelines, see GUIDELINES for details.
//
#ifndef HERWIG_Histogram_H
#define HERWIG_Histogram_H
//
// This is the declaration of the Histogram class.
//
#include "Histogram.fh"
#include "ThePEG/Interface/Interfaced.h"
#include "Statistic.h"
#include <string>
 
namespace Herwig {
 
using namespace ThePEG;
 
  namespace HistogramOptions {
    const unsigned int None      = 0;      
    const unsigned int Frame     = 1;      
    const unsigned int Errorbars = 1 << 1; 
    const unsigned int Xlog      = 1 << 2; 
    const unsigned int Ylog      = 1 << 3; 
    const unsigned int Smooth    = 1 << 4; 
    const unsigned int Rawcount  = 1 << 5; 
  }
 
class Histogram: public Interfaced {
 
public:
 
  Histogram(double lower=0., double upper=0., unsigned int nbin=0)
  : _globalStats(), _havedata(false), _bins(nbin+2),_prefactor(1.),_total(0.) {
    if (upper<lower) swap(upper,lower);
    _bins[0].limit=-1.e100;
    double limit(lower);
    double width((upper-lower)/nbin);
    for(unsigned int ix=1; ix <= nbin; ++ix) {
      _bins[ix].limit=limit;
      limit += width;
    }
    _bins.back().limit=limit;
  }
 
  Histogram(vector<double> limits) 
    : _globalStats(), _havedata(false), _bins(limits.size()+1), _prefactor(1.),_total(0.) {
    _bins[0].limit=-1.e100;
    for (size_t i=1; i<=limits.size(); ++i)
      _bins[i].limit=limits[i-1];
  }
 
  Histogram(vector<double> limits, vector<double> data, vector<double> dataerror)
    : _globalStats(), _havedata(true), _bins(limits.size()+1), _prefactor(1.),_total(0.) {
    _bins[0].limit=-1.e100;
    for (size_t i=1; i<=limits.size(); ++i)
      _bins[i].limit=limits[i-1];
    
    // no data goes into _bins[0] or _bins.back()!
    for (size_t i=1; i<=min(limits.size()-1,data.size()); ++i)
      _bins[i].data=data[i-1];
    
    for (size_t i=1; i<=min(limits.size()-1,dataerror.size()); ++i)
      _bins[i].dataerror=dataerror[i-1];
  }
  
 
public:
 
  void operator += (double input) {
    addWeighted(input,1.0);
  }
 
  void addWeighted(double input, double weight) {
    if(std::isnan(input)) return;
    unsigned int ibin;
    for(ibin=1; ibin<_bins.size(); ++ibin) {
      if(input<_bins[ibin].limit)
        break;
    }
    _bins[ibin-1].contents   += weight;
    _bins[ibin-1].contentsSq += sqr(weight);
    _globalStats += weight * input;
    _total += weight;
  }
 
  unsigned int numberOfBins() const { 
    return _bins.size()-2;
  }
 
 
  double prefactor() const {
    return _prefactor;
  }
 
  void   prefactor(double in ) {
    _prefactor=in;
  }
 
  const Statistic & globalStatistics() const {
    return _globalStats;
  }
 
  void normaliseToData();
 
  void normaliseToCrossSection();
 
  void chiSquared(double & chisq, 
                  unsigned int & ndegrees, double minfrac=0.) const;
 
  void rivetOutput(ostream & out,
                   string histogramname = string("default"),
                   string analysisname = string("default"),
                   string title  = string(),
                   string xlabel = string(),
                   string ylabel = string(),
                   bool rawcount = false,
                   double multiplicator = 1.0) const;
 
  void topdrawOutput(ostream & out,
                     unsigned int flags = 0,
                     string colour = string("BLACK"),
                     string title = string(),
                     string titlecase = string(),
                     string left = string(),
                     string leftcase = string(),
                     string bottom = string(),
                     string bottomcase = string()
                     ) const;
 
  void topdrawMCatNLO(ostream & out,
                      unsigned int flags =0 ,
                      string colour = string("BLACK"),
                      string title = string()
                      ) const;
 
  void topdrawOutputAverage(ostream & out,
                            bool frame,
                            bool errorbars,
                            bool xlog, bool ylog,
                            string colour=string("BLACK"),
                            string title=string(),
                            string titlecase =string(),
                            string left=string(),
                            string leftcase =string(),
                            string bottom=string(),
                            string bottomcase =string()) const;
 
  unsigned int visibleEntries() const;
 
  double dataNorm() const;
 
  void simpleOutput(ostream & out, bool errorbars, bool normdata=false);
 
  vector<double> dumpBins() const;
 
  Histogram ratioWith(const Histogram & h2) const;
 
 
  static
  vector<double> LogBins(double xmin, unsigned nbins, double base = 10.0);
 
 
public:
 
  static void Init();
 
protected:
 
  virtual IBPtr clone() const;
 
  virtual IBPtr fullclone() const;
 
private:
 
  Histogram & operator=(const Histogram &) = delete;
 
private:
 
  Statistic _globalStats;
 
  bool _havedata;
 
  struct Bin {
    Bin() : contents(0.0), contentsSq(0.0), 
            limit(0.0), data(0.0), dataerror(0.0), points(0) {}
    double contents;
 
    double contentsSq;
 
    double limit;
 
    double data;
 
    double dataerror;
 
    long points;
  };
 
  vector<Bin> _bins;
 
  double _prefactor;
 
  double _total;
 
 
public:
 
  vector<Bin> bins() const { return _bins; }
 
};
 
}
 
#endif /* HERWIG_Histogram_H */