#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[2];
  };
  obs_phi2(o815 *_O815);
  
private:
  void _start();
  void _meas(bool loadedobs, const int& nthmeas);
  void _finish();
 
  obsmem* OM;

  void phi2Compute();

  static void prePhi2Sus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para);
  static double phi2Sus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para);

  sim *Sim;
  obstat<double,double> oPhi2[2];
};

obs_phi2::obs_phi2(o815 *_O815) : o815::obs("phi2", 
					    _O815->paraQ->getParaNames() + 
					    "phi2one:phi2one_err:phi2onesus:phi2onesus_err"
					    ":phi2two:phi2two_err:phi2twosus:phi2twosus_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();

  for (int iflav=0; iflav<2; iflav++)
    oPhi2[iflav].addMeas(OM->phi2[iflav]);
};

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

  int compid_phi2[2], compid_phi2sus[2];

  for (int iflav=0; iflav<2; iflav++) {
    compid_phi2[iflav] = oPhi2[iflav].computeEasy();
    compid_phi2sus[iflav] = oPhi2[iflav].computeJack(obs_phi2::prePhi2Sus, obs_phi2::phi2Sus, &(Sim->lsize4));

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

    oPhi2[iflav].reset();
  }
  
  *out << endl;
};

void obs_phi2::phi2Compute()
{
  for (int iflav=0; iflav<2; iflav++) {
    OM->phi2[iflav] = 0;
    
    for (int ix = 0; ix < Sim->lsize4; ix++)
      OM->phi2[iflav] += norm(Sim->phi[ iflav*Sim->lsize4 + ix ]);

    OM->phi2[iflav] /= Sim->lsize4;
  }
}

void obs_phi2::prePhi2Sus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para) {
  preCalculated->push_back(0);
  preCalculated->push_back(0);

  for(vector< vector<double> >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) {
    (*preCalculated)[0] += (*valIt)[0];
    (*preCalculated)[1] += pow((*valIt)[0],2);
  }
}

double obs_phi2::phi2Sus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para) {
  return (( (*preCalculated)[1] - pow((*excludedMeas)[0], 2) ) / (nmeas-1)
	  - pow( ( (*preCalculated)[0] - (*excludedMeas)[0] ) / (nmeas-1), 2 ))
    * *(int*)para;
}

#endif