Added observable phi2.

[phys/heatbath.git] / obs_phi2.hpp

diff --git a/obs_phi2.hpp b/obs_phi2.hpp
new file mode 100644 (file)
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