From 532d407e33ab23eb1430c7f42987b1711b00ac4e Mon Sep 17 00:00:00 2001
From: Alexander Schmidt <alex@treefish.org>
Date: Wed, 13 Nov 2013 11:08:39 +0100
Subject: [PATCH] Added observable phi2.

---
 heatbath.cpp | 14 +++++++-
 obs_phi2.hpp | 93 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 106 insertions(+), 1 deletion(-)
 create mode 100644 obs_phi2.hpp

diff --git a/heatbath.cpp b/heatbath.cpp
index 70e2140..de99650 100644
--- a/heatbath.cpp
+++ b/heatbath.cpp
@@ -2,6 +2,8 @@
 
 #include "sim.hpp"
 
+#include "obs_phi2.hpp"
+
 o815 *O815;
 sim *Sim;
 
@@ -25,6 +27,16 @@ void helpHeader()
   cout << "Usage: ./heatbath [OPTIONS] [obs1] ... [obsN]" << endl << endl;
 }
 
+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));
+    } 
+  }
+}
+
 int main (int argc, char *argv[])
 {
   O815 = new o815(argc, argv, "heatbath", specOps, &helpHeader);
@@ -37,7 +49,7 @@ int main (int argc, char *argv[])
 
   O815->Sim = new sim(O815);
 
-  //parseLonelyArgs();
+  parseLonelyArgs();
 
   O815->mainLoop();
 
diff --git a/obs_phi2.hpp b/obs_phi2.hpp
new file mode 100644
index 0000000..23a84e4
--- /dev/null
+++ b/obs_phi2.hpp
@@ -0,0 +1,93 @@
+#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 phi2;
+  };
+  obs_phi2(o815 *_O815);
+  
+private:
+  void _start();
+  void _meas(bool loadedobs);
+  void _finish();
+ 
+  obsmem* OM;
+
+  void phi2Compute();
+  static double phi2Sus(vector< vector<double> > *vals, void *para);
+
+  sim *Sim;
+  obstat<double,double> oPhi2;
+};
+
+obs_phi2::obs_phi2(o815 *_O815) : o815::obs("phi2", 
+					    _O815->paraQ->getParaNames() + "phi2:phi2_err:phi2sus:phi2sus_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) {
+  if (!loadedobs)
+    phi2Compute();
+
+  oPhi2.addMeas(OM->phi2);
+};
+
+void obs_phi2::_finish() {
+  *out << O815->paraQ->getParaVals();
+
+  int compid_phi2, compid_phi2sus;
+
+  compid_phi2 = oPhi2.computeEasy();
+  compid_phi2sus = oPhi2.computeJack(obs_phi2::phi2Sus, &(Sim->LSIZE4));
+
+  *out << "\t" << oPhi2.getMean(compid_phi2) << "\t" << oPhi2.getErr(compid_phi2);
+  *out << "\t" << oPhi2.getMean(compid_phi2sus) << "\t" << oPhi2.getErr(compid_phi2sus);
+
+  *out << endl;
+
+  oPhi2.reset();
+};
+
+void obs_phi2::phi2Compute()
+{
+  OM->phi2 = 0;
+
+  for (int ix = 0; ix < Sim->LSIZE4; ix++)
+    OM->phi2 += norm( Sim->conf[ix].phi );
+
+  OM->phi2 /= Sim->LSIZE4;
+}
+
+double obs_phi2::phi2Sus(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;
+}
+
+#endif
-- 
2.39.5