obs_phip2_hist.hpp

   1 #ifndef OBS_PHIP2_HIST_HPP
   2 #define OBS_PHIP2_HIST_HPP
   3
   4 #include "latlib/o815/o815.h"
   5
   6 #include "latlib/writeout.h"
   7
   8 #include "latlib/obstat.hpp"
   9
  10 #include <iostream>
  11
  12 using namespace std;
  13
  14 class obs_phip2_hist : public o815::obs {
  15
  16 public:
  17   struct obsmem {
  18     double *phip2;
  19   };
  20   obs_phip2_hist(o815 *_O815);
  21
  22 private:
  23   void _start() {}
  24   void _meas(bool loadedobs, const int& nthmeas);
  25   void _finish() {}
  26   obsmem OM;
  27   sim *Sim;
  28 };
  29
  30 obs_phip2_hist::obs_phip2_hist(o815 *_O815) : o815::obs("phip2",
  31                                                         _O815->paraQ->getParaNames() + "nthstep:kpt:kpx:phip2", _O815,
  32                                                         sizeof(double)*_O815->comargs.lsize[0]*_O815->comargs.lsize[1],
  33                                                         "_hist") {
  34   OM.phip2 = (double*)obsMem;
  35   Sim = (sim*)O815->Sim;
  36 }
  37
  38 void obs_phip2_hist::_meas(bool loadedobs, const int& nthmeas) {
  39   if (!loadedobs) {
  40     for (int kpt = 0; kpt < O815->comargs.lsize[0]; kpt++)
  41       for (int kpx = 0; kpx < O815->comargs.lsize[1]; kpx++) {
  42         double ppt = 2.*M_PI/O815->comargs.lsize[0] * kpt;
  43         double ppx = 2.*M_PI/O815->comargs.lsize[1] * kpx;
  44
  45         OM.phip2[ kpt * O815->comargs.lsize[1] + kpx ] = 0;
  46
  47         for (int ixt = 0; ixt < O815->comargs.lsize[0]; ixt++)
  48           for (int ixx = 0; ixx < O815->comargs.lsize[1]; ixx++)
  49             for (int ixpt = 0; ixpt < O815->comargs.lsize[0]; ixpt++)
  50               for (int ixpx = 0; ixpx < O815->comargs.lsize[1]; ixpx++)
  51                 OM.phip2[ kpt * O815->comargs.lsize[1] + kpx ] += real( conj( Sim->conf[ ixt*O815->comargs.lsize[1] + ixx ].phi )
  52                                                                         * Sim->conf[ ixpt*O815->comargs.lsize[1] + ixpx ].phi
  53                                                                         * exp ( _i_*(double)ppx*(double)(ixx-ixpx) + _i_*(double)ppt*(double)(ixt-ixpt) ) );
  54
  55         OM.phip2[ kpt * O815->comargs.lsize[1] + kpx ] /= Sim->LSIZE2;
  56       }
  57   }
  58
  59   for (int kpt = 0; kpt <= O815->comargs.lsize[0]; kpt++)
  60     for (int kpx = 0; kpx <= O815->comargs.lsize[1]; kpx++) {
  61       double ppt = 2.*M_PI/O815->comargs.lsize[0] * kpt;
  62       double ppx = 2.*M_PI/O815->comargs.lsize[1] * kpx;
  63
  64       *out << O815->paraQ->getParaVals();
  65       *out << "\t" << ( Sim->nequi + (nthmeas+1)*Sim->nskip );
  66       *out << "\t" << kpt << "\t" << kpx;
  67       *out << "\t" << OM.phip2[ kpt * O815->comargs.lsize[1] + kpx ] - 1./(ppt*ppt + ppx*ppx + Sim->m*Sim->m) << endl;
  68     }
  69 };
  70
  71 #endif