#include "histo.h"
#include "run.h"
#include "bin.h"
#include <string>
#include <iostream>
#include <iomanip>
#include <math.h>
#include <TMatrixD.h>
#include <TVectorD.h>
#include <TTree.h>

using std::string;
using namespace std;

Histo::Histo(string def) {
  sscanf (def.c_str(), "%d %s %d %lf %lf", &type, var_name, &no_of_bins, &min, &max);
  bins = new Bin [no_of_bins];
  out_of_bounds = 0;
}

Histo::~Histo() {
  delete [] bins;
}

void Histo::init(Run *run) {
  pointer_to_x = run->get_address(var_name);
}

void Histo::add(Run *run) {
  if (this != NULL) {
    Double_t x = *pointer_to_x;
//    stat_x.add(x);
    int i = int((x - min) / (max - min) * no_of_bins);
    if (i >= 0 && i < no_of_bins) bins[i].add(x, run);
    else out_of_bounds++;
  }
}

/*
void  Histo::out() {

  cout << "Histogram: " << var_name << endl << endl;
  cout << "OOB: " << out_of_bounds << endl;
  cout << "Stat: " << stat_x.get_mean() << " " << stat_x.get_min() << " " << stat_x.get_max() << endl;
  for (int n = 0; n < no_of_bins; n++) {
    cout << setw(3) << n;
    cout << setw(10) << bins[n].get_n();
    cout << fixed << setprecision(5);
    cout << setw(14) << bins[n].get_Px() << setw(10) << bins[n].get_dPx();
    cout << setw(14) << bins[n].get_Py() << setw(10) << bins[n].get_dPy();
    cout << setw(14) << bins[n].get_Pz() << setw(10) << bins[n].get_dPz();
    cout << setw(14) << bins[n].get_R() << setw(10) << bins[n].get_dR();
    cout << endl;
  }
  cout << endl;
  for (int n = 0; n < no_of_bins; n++) {
    cout << setw(3) << n;
    cout << setw(10) << bins[n].get_n();
    cout << fixed << setprecision(5);
    cout << setw(14) << bins[n].get_Px() << setw(10) << bins[n].get_dPx();
    cout << setw(14) << bins[n].get_Py() << setw(10) << bins[n].get_dPy();
    cout << setw(14) << bins[n].get_Pz() << setw(10) << bins[n].get_dPz();
    cout << setw(14) << bins[n].get_R() << setw(10) << bins[n].get_dR();
    cout << endl;
  }
}
*/

void Histo::init_all(Run *run) {
  if (this != NULL && run != NULL) {
    Histo *h = this;
    while (h != NULL) {
      h->init(run);
      h = h->next_histo;
    }
  }
}

void Histo::add_all(Run *run) {
  if (this != NULL && run != NULL) {
    Histo *h = this;
    while (h != NULL) {
      h->add(run);
      h = h->next_histo;
    }
  }
}

/*
void Histo::out_all() {
  if (this != NULL) {
    Histo *h = this;
    while (h != NULL) {
      h->out();
      h = h->next_histo;
    }
  }
}
*/
void Histo::set_next_histo(Histo* histo) {
  next_histo = histo;
}

Histo* Histo::get_next_histo() {
  return next_histo;
}

string Histo::get_name() {
  return var_name;
}


string Histo::get_type() {
  switch (type) {
  case 1: return "hist";
  case 2: return "evnt";
  case 3: return "ac";
  default: return "unknown_type";
  }
}

// ----------------------------------------------------------------------------

ostream& operator<< (ostream &os, const Histo &obj) {
  os << endl << "Histogram: " << obj.var_name << endl << endl;
  os << "OOB: " << obj.out_of_bounds << endl;

  bool out_calc = false;
  for (int n = 0; n < obj.no_of_bins; n++) {
    if (obj.bins[n].get_px_prime_calc() != 0.0 ||
	obj.bins[n].get_pz_prime_calc() != 0.0) {
      out_calc = true;
    }
  }

  // calculate the mean value over all bins of the form-factor ratio for
  // later use in the estimation of the analyzing power.

  Double_t R = 0;
  Double_t dR = 0;
  Double_t s1 = 0;
  Double_t s2 = 0;
  for (int n = 0; n < obj.no_of_bins; n++) {
    if (obj.bins[n].get_dR() > 0) {
      s1 += obj.bins[n].get_R() / pow(obj.bins[n].get_dR(), 2);
      s2 += 1. / pow(obj.bins[n].get_dR(), 2);
    }
  }
  if (s2 > 0) {
    R = s1 / s2;
    dR = sqrt(1 / s2);
  }

  // Phil: TODO: why is there a cout here?
  cout << "XXXX: " << R << " " << dR << endl;

  for (int n = 0; n < obj.no_of_bins; n++) {
    switch (obj.type) {
    case 1:
    default :
      os << setw(3) << n;
      os << setw(10) << obj.bins[n].get_n();
      os << setw(13) << scientific << setprecision(3) << obj.bins[n].get_x();
      os << fixed << setprecision(4);
      os << setw(11) << obj.bins[n].get_Px() << setw(8) << obj.bins[n].get_dPx();
      if (out_calc) {os << setw(8) << obj.bins[n].get_px_prime_calc();};
      os << setw(11) << obj.bins[n].get_Py() << setw(8) << obj.bins[n].get_dPy();
      os << setw(11) << obj.bins[n].get_Pz() << setw(8) << obj.bins[n].get_dPz();
      if (out_calc) {os << setw(8) << obj.bins[n].get_pz_prime_calc();};
      os << setw(11) << obj.bins[n].get_R() << setw(8) << obj.bins[n].get_dR();
      if (out_calc) {os << setw(8) << obj.bins[n].get_r_calc();};

      os << setw(11) << obj.bins[n].get_a1() << setw(8) << obj.bins[n].get_da1();
      os << setw(8) << obj.bins[n].get_b1() << setw(8) << obj.bins[n].get_db1();

      os << endl;
      break;

    case 2:
      os << setw(13) << scientific << setprecision(3) << obj.bins[n].get_x();
      os << setw(10) << obj.bins[n].get_np();
      os << setw(10) << obj.bins[n].get_nm();
      os << endl;
      break;

    case 3:
      os << setw(13) << scientific << setprecision(3) << obj.bins[n].get_x();
      os << setw(10) << obj.bins[n].get_n();
      os << fixed << setprecision(4) << setw(13) << obj.bins[n].get_Ac();
      os << fixed << setprecision(4) << setw(10) << obj.bins[n].get_dAc();
      os << fixed << setprecision(4) << setw(13) << obj.bins[n].get_Ac(R, dR);
      os << fixed << setprecision(4) << setw(10) << obj.bins[n].get_dAc(R, dR);
      os << fixed << setprecision(4) << setw(13) << obj.bins[n].get_Ac_calc();
      os << endl;
      break;
    }
  }

  os << endl;
  return os;
}