#ifndef OBS_PHIP2_HIST_HPP
#define OBS_PHIP2_HIST_HPP

#include "latlib/o815/o815.h"

#include "latlib/writeout.h"

#include "latlib/obstat.hpp"

#include <iostream>
#include <stdlib.h>

using namespace std;

class obs_phip2_hist : public o815::obs {

public:
  struct obsmem {
    double phip2[4];
  };
  obs_phip2_hist(o815 *_O815, const int& avgmany=2);
  
private:
  void _start() {}
  void _meas(bool loadedobs, const int& nthmeas);
  void _finish() {}
  obsmem *OM;
  sim *Sim;
  double avgphip2[4];
  static string intToStr(const int& number);
  int avgmany;
};

string obs_phip2_hist::intToStr(const int& number)
{
  stringstream ss;
  ss << number;
  return ss.str();
}

obs_phip2_hist::obs_phip2_hist(o815 *_O815, const int& _avgmany) : o815::obs("phip2", 
									     _O815->paraQ->getParaNames() + "nthstep:imode:kt:kx:phip2absdiffequi", _O815, 
									     sizeof(obsmem), 
									     "_hist_avg" + intToStr(_avgmany) ) {
  OM = (obsmem*)obsMem;
  Sim = (sim*)O815->Sim;
  avgmany = _avgmany;
}

void obs_phip2_hist::_meas(bool loadedobs, const int& nthmeas) {
  for (int kpimode=0; kpimode<4; kpimode++) {
    const double ppt = 2.*M_PI/O815->comargs.lsize[0] * ( 0 + int( kpimode/4. * O815->comargs.lsize[0] ) ) ;
    const double ppx = 2.*M_PI/O815->comargs.lsize[1] * ( 0 + int( kpimode/4. * O815->comargs.lsize[1] ) );

    if (!loadedobs) {
      complex<double> phitildepp = 0;
      
      for (int ixpt = 0; ixpt < O815->comargs.lsize[0]; ixpt++)
	for (int ixpx = 0; ixpx < O815->comargs.lsize[1]; ixpx++)
	  phitildepp += Sim->conf[ ixpt*O815->comargs.lsize[1] + ixpx ].phi 
	    * exp ( - _i_*(double)ppx*(double)ixpx - _i_*(double)ppt*(double)ixpt );
      
      OM->phip2[ kpimode ] = norm( phitildepp) / Sim->LSIZE2;
    }

    if ( nthmeas%avgmany == 0 )
      avgphip2[kpimode] = 0;
    
    avgphip2[kpimode] += OM->phip2[ kpimode ];

    if ( (nthmeas+1)%avgmany == 0 ) {
      *out << O815->paraQ->getParaVals();
      *out << "\t" << ( Sim->nequi + (nthmeas+1)*Sim->nskip );
      *out << "\t" << kpimode;
      *out << "\t" << 0 + int( kpimode/4. * O815->comargs.lsize[0] );
      *out << "\t" << 0 + int( kpimode/4. * O815->comargs.lsize[1] );
      *out << "\t" << abs( avgphip2[kpimode] / (double)avgmany - 1./(ppt*ppt + ppx*ppx + Sim->m*Sim->m) ) << endl;
    }
  }
};

#endif