Migrated observable plaq.

diff --git a/u1casc-ordinary/obs_plaq.hpp b/u1casc-ordinary/obs_plaq.hpp
new file mode 100644 (file)
index 0000000..5dacf68
--- /dev/null
+++ b/u1casc-ordinary/obs_plaq.hpp
@@ -0,0 +1,97 @@
+#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 plaq;
+  };
+  obs_plaq(o815 *_O815);
+  
+private:
+  void _start();
+  void _meas(bool loadedobs, const int& nthmeas);
+  void _finish();
+ 
+  obsmem* OM;
+
+  void plaqCompute();
+  static double plaqSus(vector< vector<double> > *vals, void *para);
+
+  sim *Sim;
+  obstat<double,double> oPlaq;
+};
+
+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)
+    plaqCompute();
+
+  oPlaq.addMeas(OM->plaq);
+};
+
+void obs_plaq::_finish() {
+  *out << O815->paraQ->getParaVals();
+
+  int compid_plaq, compid_plaqsus;
+
+  compid_plaq = oPlaq.computeEasy();
+  compid_plaqsus = oPlaq.computeJack(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::plaqCompute()
+{
+  OM->plaq = 0;
+  
+  for( int ix=0; ix<Sim->lsize4; ix++ )
+    for( int inu=0; inu<3; inu++ )
+      for( int isigma=inu+1; isigma<4; isigma++ ) {
+       OM->plaq += real( Sim->U[ ix*4 + inu ] * Sim->U[ (*Sim->nb)[ix*8+inu]*4 + isigma ]  
+                         * conj(Sim->U[ (*Sim->nb)[ix*8+isigma]*4 + inu ]) * conj(Sim->U[ ix*4 + isigma ]) );
+      }
+  
+  OM->plaq /= Sim->lsize4*6;
+}
+
+double obs_plaq::plaqSus(vector< vector<double> > *vals, void *para) {
+  double mean=0, mean2=0;
+  for(vector< vector<double> >::iterator valIt = vals->begin(); valIt != vals->end(); ++valIt)  
+    {
+      mean += (*valIt)[0];
+      mean2 += pow((*valIt)[0],2);
+    }
+  mean /= vals->size();
+  mean2 /= vals->size();
+  
+  return ( mean2 - pow(mean,2) ) * *(int*)para * 6;
+}
+
+#endif

diff --git a/u1casc-ordinary/u1casc-ordinary.cpp b/u1casc-ordinary/u1casc-ordinary.cpp
index 4956ee1b304e1fd2d9366eb24a0557570f1a505b..554f7e4281d1ec271614e25f46d340b130ed6d1b 100644 (file)
--- a/u1casc-ordinary/u1casc-ordinary.cpp
+++ b/u1casc-ordinary/u1casc-ordinary.cpp
@@ -11,6 +11,7 @@ sim *Sim;
 //const complex<double> _i_ = complex<double>(0.0,1.0);
 
 #include "obs_phi2.hpp"
+#include "obs_plaq.hpp"
 
 o815::comoption specOps[] = {
   { "kappaone", required_argument, NULL, 'r', "set inverse mass kappa_1", "min:max:inc" },
@@ -65,16 +66,10 @@ void parseLonelyArgs()
       *O815->out->log << "MASTER: registered observable: phi2" << endl << flush;
       O815->observables.push_back(new obs_phi2(O815));
     }
-    /*
-    else if ( strcmp(*lonit, "phi2_hist") == 0 ) {
-      *O815->out->log << "MASTER: registered observable: phi2_hist" << endl << flush;
-      O815->observables.push_back(new obs_phi2_hist(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, "phip2_hist") == 0 ) {
-      *O815->out->log << "MASTER: registered observable: phip2_hist" << endl << flush;
-      O815->observables.push_back(new obs_phip2_hist(O815));
-    }
-    */
   }
 }