#include "analysis.h"
#include <cmath>
#include <algorithm>
#include <iterator>
#include <fstream>

using namespace school;
using namespace std;

vector<double> analysis::regular_bin_edges(double lower, double upper, size_t n_bins) {
  vector<double> res;
  double step = (upper - lower)/n_bins;
  for(size_t k = 0; k <= n_bins; ++k )
    res.push_back(lower + k*step);
  return res;
}

vector<double> analysis::log_bin_edges(double lower, double upper, size_t n_bins) {
  vector<double> res;
  double step = log10(upper/lower)/n_bins;
  for (size_t k = 0; k <= n_bins; ++k )
    res.push_back(lower*pow(10.0,k*step));
  return res;
}

analysis::analysis()
: _M_name("unknown") {}

analysis::~analysis() {}

analysis::analysis(const string& name, const vector<double>& edges)
: _M_name(name) {
  for ( vector<double>::const_iterator e = edges.begin() + 1; e != edges.end(); ++e )
    _M_bins[*e] = bin(*(e-1),*e);
}

void analysis::book(double observable, double w) {
  if ( observable < _M_bins.begin()->first ) return;
  map<double,bin>::iterator b = _M_bins.upper_bound(observable);
  if ( b == _M_bins.end() ) return;
  b->second.count(w);
}

analysis::dsigma_point::dsigma_point(const analysis::bin& b, unsigned long n_points) {
  value = 0.5*(b.lower_bound + b.upper_bound);
  double width = b.upper_bound-b.lower_bound;
  cross_section = b.sum_of_weights/width/n_points;
  cross_section_error =
    b.sum_of_squared_weights/n_points -
    b.sum_of_weights*b.sum_of_weights/n_points/n_points;
  cross_section_error /= n_points - 1;
  cross_section_error = sqrt(cross_section_error);
}

list<analysis::dsigma_point> analysis::make_dsigma(unsigned long n_points) const {
  list<dsigma_point> res;
  for ( map<double,bin>::const_iterator b = _M_bins.begin();
	b != _M_bins.end(); ++b ) {
    res.push_back(dsigma_point(b->second,n_points));
  }
  return res;
}

void analysis::output_dsigma(unsigned long n_points) const {
  list<dsigma_point> points = make_dsigma(n_points);
  ofstream out((_M_name+".dat").c_str());
  copy(points.begin(),points.end(),ostream_iterator<dsigma_point>(out,"\n"));
}