#ifndef OBS_PHI2_HPP
#define OBS_PHI2_HPP

#include "latlib/o815/o815.h"

#include "latlib/writeout.h"

#include "latlib/obstat.hpp"

#include <iostream>

using namespace std;

class obs_phi2 : public o815::obs {

public:
  struct obsmem {
    double phi2;
  };
  obs_phi2(o815 *_O815);
  
private:
  void _start();
  void _meas(bool loadedobs, const int& nthmeas);
  void _finish();
 
  obsmem* OM;

  void phi2Compute();
  static double phi2Sus(vector< vector<double> > *vals, void *para);

  sim *Sim;
  obstat<double,double> oPhi2;
};

obs_phi2::obs_phi2(o815 *_O815) : o815::obs("phi2", 
					    _O815->paraQ->getParaNames() + "phi2:phi2_err:phi2sus:phi2sus_err", 
					    _O815, sizeof(obsmem) ) {
  OM = (obsmem*)obsMem;
  Sim = (sim*)O815->Sim;
}

void obs_phi2::_start() {
  //*out << "OBS_test: start" << endl;
};

void obs_phi2::_meas(bool loadedobs, const int& nthmeas) {
  if (!loadedobs)
    phi2Compute();

  oPhi2.addMeas(OM->phi2);
};

void obs_phi2::_finish() {
  *out << O815->paraQ->getParaVals();

  int compid_phi2, compid_phi2sus;

  compid_phi2 = oPhi2.computeEasy();
  compid_phi2sus = oPhi2.computeJack(obs_phi2::phi2Sus, &(Sim->LSIZE2));

  *out << "\t" << oPhi2.getMean(compid_phi2) << "\t" << oPhi2.getErr(compid_phi2);
  *out << "\t" << oPhi2.getMean(compid_phi2sus) << "\t" << oPhi2.getErr(compid_phi2sus);

  *out << endl;

  oPhi2.reset();
};

void obs_phi2::phi2Compute()
{
  OM->phi2 = 0;

  for (int ix = 0; ix < Sim->LSIZE2; ix++)
    OM->phi2 += norm( Sim->conf[ix].phi );

  OM->phi2 /= Sim->LSIZE2;
}

double obs_phi2::phi2Sus(vector< vector<double> > *vals, void *para) {
  double mean=0, mean2=0;
  for(vector< vector<double> >::iterator valIt = vals->begin(); valIt != vals->end(); ++valIt)  
    {
      mean += (*valIt)[0];
      mean2 += pow((*valIt)[0],2);
    }
  mean /= vals->size();
  mean2 /= vals->size();
  
  return ( mean2 - pow(mean,2) ) * *(int*)para;
}

#endif