u1casc-ordinary/obs_plaq.hpp

   1 #ifndef OBS_PLAQ_HPP
   2 #define OBS_PLAQ_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_plaq : public o815::obs {
  15
  16 public:
  17   struct obsmem {
  18     double plaq;
  19   };
  20   obs_plaq(o815 *_O815);
  21
  22 private:
  23   void _start();
  24   void _meas(bool loadedobs, const int& nthmeas);
  25   void _finish();
  26
  27   obsmem* OM;
  28
  29   void plaqCompute();
  30
  31   static void prePlaqSus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para);
  32   static double plaqSus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para);
  33
  34   sim *Sim;
  35   obstat<double,double> oPlaq;
  36 };
  37
  38 obs_plaq::obs_plaq(o815 *_O815) : o815::obs("plaq",
  39                                             _O815->paraQ->getParaNames() + "plaq:plaq_err:plaqsus:plaqsus_err",
  40                                             _O815, sizeof(obsmem) ) {
  41   OM = (obsmem*)obsMem;
  42   Sim = (sim*)O815->Sim;
  43 }
  44
  45 void obs_plaq::_start() {
  46   //*out << "OBS_test: start" << endl;
  47 };
  48
  49 void obs_plaq::_meas(bool loadedobs, const int& nthmeas) {
  50   if (!loadedobs)
  51     plaqCompute();
  52
  53   oPlaq.addMeas(OM->plaq);
  54 };
  55
  56 void obs_plaq::_finish() {
  57   *out << O815->paraQ->getParaVals();
  58
  59   int compid_plaq, compid_plaqsus;
  60
  61   compid_plaq = oPlaq.computeEasy();
  62   compid_plaqsus = oPlaq.computeJack(obs_plaq::prePlaqSus, obs_plaq::plaqSus, &(Sim->lsize4));
  63
  64   *out << "\t" << oPlaq.getMean(compid_plaq) << "\t" << oPlaq.getErr(compid_plaq);
  65   *out << "\t" << oPlaq.getMean(compid_plaqsus) << "\t" << oPlaq.getErr(compid_plaqsus);
  66
  67   *out << endl;
  68
  69   oPlaq.reset();
  70 };
  71
  72 void obs_plaq::plaqCompute()
  73 {
  74   OM->plaq = 0;
  75
  76   for( int ix=0; ix<Sim->lsize4; ix++ )
  77     for( int inu=0; inu<3; inu++ )
  78       for( int isigma=inu+1; isigma<4; isigma++ ) {
  79         OM->plaq += real( Sim->U[ ix*4 + inu ] * Sim->U[ (*Sim->nb)[ix*8+inu]*4 + isigma ]
  80                           * conj(Sim->U[ (*Sim->nb)[ix*8+isigma]*4 + inu ]) * conj(Sim->U[ ix*4 + isigma ]) );
  81       }
  82
  83   OM->plaq /= Sim->lsize4*6;
  84 }
  85
  86 void obs_plaq::prePlaqSus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para) {
  87   preCalculated->push_back(0);
  88   preCalculated->push_back(0);
  89
  90   for(vector< vector<double> >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) {
  91     (*preCalculated)[0] += (*valIt)[0];
  92     (*preCalculated)[1] += pow((*valIt)[0],2);
  93   }
  94 }
  95
  96 double obs_plaq::plaqSus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para) {
  97   return (( (*preCalculated)[1] - pow((*excludedMeas)[0], 2) ) / (nmeas-1)
  98           - pow( ( (*preCalculated)[0] - (*excludedMeas)[0] ) / (nmeas-1), 2 ))
  99     * *(int*)para * 6;
 100 }
 101
 102 #endif