From 7664dbeaba693ea47d5fd405c20f5db109d4da13 Mon Sep 17 00:00:00 2001
From: Alexander Schmidt <alex@treefish.org>
Date: Thu, 27 Feb 2014 22:32:17 +0100
Subject: [PATCH 1/1] Migrated observable phi2.

---
 u1casc-flux/obs_phi2.hpp    | 118 ++++++++++++++++++++++++++++++++++++
 u1casc-flux/sim.hpp         |   4 +-
 u1casc-flux/u1casc-flux.cpp |  10 +--
 3 files changed, 123 insertions(+), 9 deletions(-)
 create mode 100644 u1casc-flux/obs_phi2.hpp

diff --git a/u1casc-flux/obs_phi2.hpp b/u1casc-flux/obs_phi2.hpp
new file mode 100644
index 0000000..dc499c4
--- /dev/null
+++ b/u1casc-flux/obs_phi2.hpp
@@ -0,0 +1,118 @@
+#ifndef OBS_PHI2_HPP
+#define OBS_PHI2_HPP
+
+#include "latlib/o815/o815.h"
+
+#include "latlib/writeout.h"
+
+#include "latlib/obstat.hpp"
+
+#include <iostream>
+
+using namespace std;
+
+class obs_phi2 : public o815::obs {
+
+public:
+  struct obsmem {
+    double A[2];
+    double B[2];
+    double C[2];
+  };
+  obs_phi2(o815 *_O815);
+  
+private:
+  void _start();
+  void _meas(bool loadedobs, const int& nthmeas);
+  void _finish();
+ 
+  obsmem* OM;
+
+  void phi2Compute();
+
+  static void prePhi2Sus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para);
+  static double phi2Sus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para);
+
+  sim *Sim;
+  obstat<double,double> oPhi2[2];
+};
+
+obs_phi2::obs_phi2(o815 *_O815) : o815::obs("phi2", 
+					    _O815->paraQ->getParaNames() + 
+					    "phi2one:phi2one_err:phi2onesus:phi2onesus_err"
+					    ":phi2two:phi2two_err:phi2twosus:phi2twosus_err", 
+					    _O815, sizeof(obsmem) ) {
+  OM = (obsmem*)obsMem;
+  Sim = (sim*)O815->Sim;
+}
+
+void obs_phi2::_start() {
+  //*out << "OBS_test: start" << endl;
+};
+
+void obs_phi2::_meas(bool loadedobs, const int& nthmeas) {
+  if (!loadedobs)
+    phi2Compute();
+
+  for (int iflav=0; iflav<2; iflav++) {
+    double tmpMeas[3] = {OM->A[iflav], OM->B[iflav], OM->C[iflav]};
+    oPhi2[iflav].addMeas(tmpMeas, 3);
+  }
+};
+
+void obs_phi2::_finish() {
+  *out << O815->paraQ->getParaVals();
+
+  int compid_phi2[2], compid_phi2sus[2];
+
+  for (int iflav=0; iflav<2; iflav++) {
+    compid_phi2[iflav] = oPhi2[iflav].computeEasy();
+    compid_phi2sus[iflav] = oPhi2[iflav].computeJack(obs_phi2::prePhi2Sus, obs_phi2::phi2Sus, &(Sim->lsize4));
+
+    *out << "\t" << oPhi2[iflav].getMean(compid_phi2[iflav]) / Sim->lsize4 << "\t" << oPhi2[iflav].getErr(compid_phi2[iflav]) / Sim->lsize4;
+    *out << "\t" << oPhi2[iflav].getMean(compid_phi2sus[iflav]) << "\t" << oPhi2[iflav].getErr(compid_phi2sus[iflav]);
+
+    oPhi2[iflav].reset();
+  }
+  
+  *out << endl;
+};
+
+void obs_phi2::phi2Compute()
+{
+  for(int i=0; i<2; i++) {
+    OM->A[i] = 0;
+    OM->B[i] = 0;
+    OM->C[i] = 0;
+    
+    for (int x = 0; x < Sim->lsize4; x++) {
+      double wp2w;
+      
+      wp2w = Sim->WF( Sim->wArg(x,i) + 2, Sim->wArg(x,i), i );
+
+      OM->A[i] += wp2w;
+      OM->B[i] += Sim->WF( Sim->wArg(x,i) + 4, Sim->wArg(x,i), i );
+      OM->C[i] += pow(wp2w, 2);
+    }
+  }
+}
+
+void obs_phi2::prePhi2Sus(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] += pow( (*valIt)[0], 2 ) + (*valIt)[1] - (*valIt)[2]; 
+      (*preCalculated)[1] += (*valIt)[0];
+    }
+}
+
+double obs_phi2::phi2Sus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para) {
+  return (( (*preCalculated)[0] - pow((*excludedMeas)[0], 2) - (*excludedMeas)[1] + (*excludedMeas)[2] ) / (nmeas-1)
+	  - pow( ( (*preCalculated)[1] - (*excludedMeas)[0] ) / (nmeas-1), 2 ))
+    / *(int*)para;
+}
+
+#endif
diff --git a/u1casc-flux/sim.hpp b/u1casc-flux/sim.hpp
index 9eca3d7..0da2b8a 100644
--- a/u1casc-flux/sim.hpp
+++ b/u1casc-flux/sim.hpp
@@ -20,6 +20,8 @@ public:
   neigh *nb;
   int *l, *lp, *p;
   double kappa[2], lambda[2], Mu[2], beta;
+  int wArg(int x, int i);
+  double WF(int nom, int denom, int iflav);
 
 private:
   void _makeSweep();
@@ -33,8 +35,6 @@ private:
   int iPlaq(int i, int j);
   int lpUpdate(int x, int mu, int i);
   double I(int p);
-  int wArg(int x, int i);
-  double WF(int nom, int denom, int iflav);
   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 0061116..d0d6347 100644
--- a/u1casc-flux/u1casc-flux.cpp
+++ b/u1casc-flux/u1casc-flux.cpp
@@ -8,9 +8,7 @@
 o815 *O815;
 sim *Sim;
 
-//const complex<double> _i_ = complex<double>(0.0,1.0);
-
-//#include "obs_phi2.hpp"
+#include "obs_phi2.hpp"
 //#include "obs_plaq.hpp"
 //#include "obs_corr.hpp"
 
@@ -75,12 +73,11 @@ void helpHeader()
 
 void parseLonelyArgs()
 {
-  /*
   for (vector<char*>::iterator lonit = O815->lonelyArgs.begin(); lonit != O815->lonelyArgs.end(); ++lonit) {
     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));
@@ -88,9 +85,8 @@ void parseLonelyArgs()
     else if ( strcmp(*lonit, "corr") == 0 ) {
       *O815->out->log << "MASTER: registered observable: corr" << endl << flush;
       O815->observables.push_back(new obs_corr(O815));
-    }
+      } */
   }
-  */
 }
 
 int main (int argc, char *argv[])
-- 
2.39.5