u1casc-flux/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 A;
  19     double B;
  20   };
  21   obs_plaq(o815 *_O815);
  22
  23 private:
  24   void _start();
  25   void _meas(bool loadedobs, const int& nthmeas);
  26   void _finish();
  27
  28   obsmem* OM;
  29
  30   static void prePlaqSus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para);
  31   static double plaqSus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para);
  32
  33   sim *Sim;
  34   obstat<double,double> oPlaq;
  35
  36   double A(const double& beta);
  37   double B(const double& beta);
  38 };
  39
  40 obs_plaq::obs_plaq(o815 *_O815) : o815::obs("plaq",
  41                                             _O815->paraQ->getParaNames() + "plaq:plaq_err:plaqsus:plaqsus_err",
  42                                             _O815, sizeof(obsmem) ) {
  43   OM = (obsmem*)obsMem;
  44   Sim = (sim*)O815->Sim;
  45 }
  46
  47 void obs_plaq::_start() {
  48   //*out << "OBS_test: start" << endl;
  49 };
  50
  51 void obs_plaq::_meas(bool loadedobs, const int& nthmeas) {
  52   if (!loadedobs) {
  53     OM->A = A((*paraQ)["beta"]);
  54     OM->B = B((*paraQ)["beta"]);
  55   }
  56
  57   double tmpMeas[2] = {OM->A, OM->B};
  58   oPlaq.addMeas(tmpMeas, 2);
  59 };
  60
  61 void obs_plaq::_finish() {
  62   *out << O815->paraQ->getParaVals();
  63
  64   int compid_plaq, compid_plaqsus;
  65
  66   compid_plaq = oPlaq.computeEasy();
  67   compid_plaqsus = oPlaq.computeJack(obs_plaq::prePlaqSus, obs_plaq::plaqSus, &(Sim->lsize4));
  68
  69   *out << "\t" << oPlaq.getMean(compid_plaq) << "\t" << oPlaq.getErr(compid_plaq);
  70   *out << "\t" << oPlaq.getMean(compid_plaqsus) << "\t" << oPlaq.getErr(compid_plaqsus);
  71
  72   *out << endl;
  73
  74   oPlaq.reset();
  75 };
  76
  77 void obs_plaq::prePlaqSus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para) {
  78   preCalculated->push_back(0);
  79   preCalculated->push_back(0);
  80
  81   for(vector< vector<double> >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) {
  82     (*preCalculated)[0] += (*valIt)[1] + pow((*valIt)[0], 2)*(*(int*)para*6);
  83     (*preCalculated)[1] += (*valIt)[0];
  84   }
  85 }
  86
  87 double obs_plaq::plaqSus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para) {
  88   return ( (*preCalculated)[0] - pow((*excludedMeas)[0], 2) * *(int*)para * 6 - (*excludedMeas)[1] ) / (nmeas-1)
  89     - pow( ( (*preCalculated)[1] - (*excludedMeas)[0] ) / (nmeas-1), 2 )
  90     * *(int*)para * 6;
  91 }
  92
  93 double obs_plaq::A(const double& beta)
  94 {
  95   double A = 0;
  96
  97   for( int x=0; x<Sim->lsize4; x++ )
  98     for( int sigma=0; sigma<3; sigma++ )
  99       for( int tau=sigma+1; tau<4; tau++ ) {
 100         A += Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1 ) / Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] );
 101         A += Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] / beta;
 102       }
 103
 104   return A/(6*Sim->lsize4);
 105 }
 106
 107 double obs_plaq::B(const double& beta)
 108 {
 109   double B=0;
 110
 111   for( int x=0; x<Sim->lsize4; x++ )
 112     for( int sigma=0; sigma<3; sigma++ )
 113       for( int tau=sigma+1; tau<4; tau++ ) {
 114         B += Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 2 ) / Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] );
 115         B += ( (Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1) * Sim->I(Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1) ) /
 116           ( beta *  Sim->I(Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ]) );
 117         B -= pow(Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1 ) / Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] ), 2);
 118         B -= ( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] * Sim->I(Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1) ) /
 119           ( beta *  Sim->I(Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ]) );
 120         B -= Sim->p[ x*6 + Sim->iPlaq(sigma,tau)] / pow(beta,2);
 121       }
 122
 123   return B/(6*Sim->lsize4);
 124 }
 125
 126 #endif