Migrated observable plaq.

diff --git a/u1casc-flux/obs_plaq.hpp b/u1casc-flux/obs_plaq.hpp
new file mode 100644 (file)
index 0000000..249f091
--- /dev/null
+++ b/u1casc-flux/obs_plaq.hpp
@@ -0,0 +1,126 @@
+#ifndef OBS_PLAQ_HPP
+#define OBS_PLAQ_HPP
+
+#include "latlib/o815/o815.h"
+
+#include "latlib/writeout.h"
+
+#include "latlib/obstat.hpp"
+
+#include <iostream>
+
+using namespace std;
+
+class obs_plaq : public o815::obs {
+
+public:
+  struct obsmem {
+    double A;
+    double B;
+  };
+  obs_plaq(o815 *_O815);
+  
+private:
+  void _start();
+  void _meas(bool loadedobs, const int& nthmeas);
+  void _finish();
+ 
+  obsmem* OM;
+
+  static void prePlaqSus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para);
+  static double plaqSus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para);
+
+  sim *Sim;
+  obstat<double,double> oPlaq;
+
+  double A(const double& beta);
+  double B(const double& beta);
+};
+
+obs_plaq::obs_plaq(o815 *_O815) : o815::obs("plaq", 
+                                           _O815->paraQ->getParaNames() + "plaq:plaq_err:plaqsus:plaqsus_err", 
+                                           _O815, sizeof(obsmem) ) {
+  OM = (obsmem*)obsMem;
+  Sim = (sim*)O815->Sim;
+}
+
+void obs_plaq::_start() {
+  //*out << "OBS_test: start" << endl;
+};
+
+void obs_plaq::_meas(bool loadedobs, const int& nthmeas) {
+  if (!loadedobs) {
+    OM->A = A((*paraQ)["beta"]);
+    OM->B = B((*paraQ)["beta"]);
+  }
+
+  double tmpMeas[2] = {OM->A, OM->B};
+  oPlaq.addMeas(tmpMeas, 2);
+};
+
+void obs_plaq::_finish() {
+  *out << O815->paraQ->getParaVals();
+
+  int compid_plaq, compid_plaqsus;
+
+  compid_plaq = oPlaq.computeEasy();
+  compid_plaqsus = oPlaq.computeJack(obs_plaq::prePlaqSus, obs_plaq::plaqSus, &(Sim->lsize4));
+
+  *out << "\t" << oPlaq.getMean(compid_plaq) << "\t" << oPlaq.getErr(compid_plaq);
+  *out << "\t" << oPlaq.getMean(compid_plaqsus) << "\t" << oPlaq.getErr(compid_plaqsus);
+
+  *out << endl;
+
+  oPlaq.reset();
+};
+
+void obs_plaq::prePlaqSus(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)[1] + pow((*valIt)[0], 2)*(*(int*)para*6);
+    (*preCalculated)[1] += (*valIt)[0];
+  }
+}
+
+double obs_plaq::plaqSus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para) {
+  return ( (*preCalculated)[0] - pow((*excludedMeas)[0], 2) * *(int*)para * 6 - (*excludedMeas)[1] ) / (nmeas-1)
+    - pow( ( (*preCalculated)[1] - (*excludedMeas)[0] ) / (nmeas-1), 2 )
+    * *(int*)para * 6;
+}
+
+double obs_plaq::A(const double& beta)
+{
+  double A = 0;
+
+  for( int x=0; x<Sim->lsize4; x++ )
+    for( int sigma=0; sigma<3; sigma++ )
+      for( int tau=sigma+1; tau<4; tau++ ) {
+       A += Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1 ) / Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] );
+       A += Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] / beta;
+      }
+  
+  return A/(6*Sim->lsize4);
+}
+
+double obs_plaq::B(const double& beta)
+{
+  double B=0;
+
+  for( int x=0; x<Sim->lsize4; x++ )
+    for( int sigma=0; sigma<3; sigma++ )
+      for( int tau=sigma+1; tau<4; tau++ ) {
+       B += Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 2 ) / Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] );
+       B += ( (Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1) * Sim->I(Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1) ) / 
+         ( beta *  Sim->I(Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ]) );
+       B -= pow(Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1 ) / Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] ), 2);
+       B -= ( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] * Sim->I(Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1) ) / 
+         ( beta *  Sim->I(Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ]) );
+       B -= Sim->p[ x*6 + Sim->iPlaq(sigma,tau)] / pow(beta,2);
+      }
+  
+  return B/(6*Sim->lsize4);
+}
+
+#endif

diff --git a/u1casc-flux/sim.hpp b/u1casc-flux/sim.hpp
index 0da2b8afc76530b93ad3ec0908a65ac089856f22..7fd59aff106cd447785b33d106fd6b1fad5bc001 100644 (file)
--- a/u1casc-flux/sim.hpp
+++ b/u1casc-flux/sim.hpp
@@ -22,6 +22,8 @@ public:
   double kappa[2], lambda[2], Mu[2], beta;
   int wArg(int x, int i);
   double WF(int nom, int denom, int iflav);
   double kappa[2], lambda[2], Mu[2], beta;
   int wArg(int x, int i);
   double WF(int nom, int denom, int iflav);

+  int iPlaq(int i, int j);
+  double I(int p);

 
 private:
   void _makeSweep();
 
 private:
   void _makeSweep();


@@ -32,9 +34,7 @@ private:
   int cubeUpdate(int x, int orient);
   void getCube(int cube[3], int orient);
   int plaqUpdate(int x, int mu, int nu, int i);
   int cubeUpdate(int x, int orient);
   void getCube(int cube[3], int orient);
   int plaqUpdate(int x, int mu, int nu, int i);

-  int iPlaq(int i, int j);

   int lpUpdate(int x, int mu, int i);
   int lpUpdate(int x, int mu, int i);

-  double I(int p);

   static double wkern(double t, void *params);
   double wGsl(int n, int iflav);
   struct bcache { 
   static double wkern(double t, void *params);
   double wGsl(int n, int iflav);
   struct bcache { 

diff --git a/u1casc-flux/u1casc-flux.cpp b/u1casc-flux/u1casc-flux.cpp
index d0d63472a78f6122835c859feacfb2ef14ccecfc..64869e43321109663342fda305d927d487d6075e 100644 (file)
--- a/u1casc-flux/u1casc-flux.cpp
+++ b/u1casc-flux/u1casc-flux.cpp
@@ -9,7 +9,7 @@ o815 *O815;
 sim *Sim;
 
 #include "obs_phi2.hpp"
 sim *Sim;
 
 #include "obs_phi2.hpp"

-//#include "obs_plaq.hpp"
+#include "obs_plaq.hpp"

 //#include "obs_corr.hpp"
 
 o815::comoption specOps[] = {
 //#include "obs_corr.hpp"
 
 o815::comoption specOps[] = {


@@ -77,15 +77,17 @@ void parseLonelyArgs()
     if ( strcmp(*lonit, "phi2") == 0 ) {
       *O815->out->log << "MASTER: registered observable: phi2" << endl << flush;
       O815->observables.push_back(new obs_phi2(O815));
     if ( strcmp(*lonit, "phi2") == 0 ) {
       *O815->out->log << "MASTER: registered observable: phi2" << endl << flush;
       O815->observables.push_back(new obs_phi2(O815));

-    } /*
+    }

     else if ( strcmp(*lonit, "plaq") == 0 ) {
       *O815->out->log << "MASTER: registered observable: plaq" << endl << flush;
       O815->observables.push_back(new obs_plaq(O815));
     }
     else if ( strcmp(*lonit, "plaq") == 0 ) {
       *O815->out->log << "MASTER: registered observable: plaq" << endl << flush;
       O815->observables.push_back(new obs_plaq(O815));
     }

-    else if ( strcmp(*lonit, "corr") == 0 ) {
+    /*
+      else if ( strcmp(*lonit, "corr") == 0 ) {

       *O815->out->log << "MASTER: registered observable: corr" << endl << flush;
       O815->observables.push_back(new obs_corr(O815));
       *O815->out->log << "MASTER: registered observable: corr" << endl << flush;
       O815->observables.push_back(new obs_corr(O815));

-      } */
+      } 
+    */

   }
 }
 
   }
 }