Migrated observable phi2.

[phys/u1casc.git] / u1casc-ordinary / obs_phi2.hpp

diff --git a/u1casc-ordinary/obs_phi2.hpp b/u1casc-ordinary/obs_phi2.hpp
new file mode 100644 (file)
index 0000000..9afc390
--- /dev/null
+++ b/u1casc-ordinary/obs_phi2.hpp
@@ -0,0 +1,100 @@
+#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[2];
+  };
+  obs_phi2(o815 *_O815);
+  
+private:
+  void _start();
+  void _meas(bool loadedobs, const int& nthmeas);
+  void _finish();
+ 
+  obsmem* OM;
+
+  void phi2Compute();
+  static double phi2Sus(vector< vector<double> > *vals, 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++)
+    oPhi2[iflav].addMeas(OM->phi2[iflav]);
+};
+
+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::phi2Sus, &(Sim->lsize4));
+
+    *out << "\t" << oPhi2[iflav].getMean(compid_phi2[iflav]) << "\t" << oPhi2[iflav].getErr(compid_phi2[iflav]);
+    *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 iflav=0; iflav<2; iflav++) {
+    OM->phi2[iflav] = 0;
+    
+    for (int ix = 0; ix < Sim->lsize4; ix++)
+      OM->phi2[iflav] += norm(Sim->phi[ iflav*Sim->lsize4 + ix ]);
+
+    OM->phi2[iflav] /= 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