Migrated observable plaq.

[phys/u1casc.git] / u1casc-flux / obs_plaq.hpp

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