Migrated observable phi2.

diff --git a/u1casc-ordinary/.gitignore b/u1casc-ordinary/.gitignore
index 332a25adcf98e9a80a85fccc5d1cba74afc67aee..89fd3b4ba0e67600bccebd9c5987fddce7f7ed82 100644 (file)
--- a/u1casc-ordinary/.gitignore
+++ b/u1casc-ordinary/.gitignore
@@ -4,3 +4,5 @@ Makefile
 cmake_install.cmake
 CMakeCache.txt
 u1casc-ordinary
+data
+out

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

diff --git a/u1casc-ordinary/u1casc-ordinary.cpp b/u1casc-ordinary/u1casc-ordinary.cpp
index f6f05fd13cf91aebcdbdd806c2b92547f4590ad3..4956ee1b304e1fd2d9366eb24a0557570f1a505b 100644 (file)
--- a/u1casc-ordinary/u1casc-ordinary.cpp
+++ b/u1casc-ordinary/u1casc-ordinary.cpp
@@ -10,14 +10,16 @@ sim *Sim;
 
 //const complex<double> _i_ = complex<double>(0.0,1.0);
 
-//#include "obs_phi2.hpp"
+#include "obs_phi2.hpp"
 
 o815::comoption specOps[] = {
-  { "kappaone", required_argument, NULL, 'k', "set inverse mass kappa_1", "min:max:inc" },
-  { "kappatwo", required_argument, NULL, 'm', "set inverse mass kappa_2", "min:max:inc" },
+  { "kappaone", required_argument, NULL, 'r', "set inverse mass kappa_1", "min:max:inc" },
+  { "kappatwo", required_argument, NULL, 's', "set inverse mass kappa_2", "min:max:inc" },
   { "beta", required_argument, NULL, 'b', "set inverse gauge coupling beta", "min:max:inc" },
   { "lambdaone", required_argument, NULL, 'p', "set quartic coupling lambda_1", "min:max:inc" },
   { "lambdatwo", required_argument, NULL, 'q', "set quartic coupling lambda_2", "min:max:inc" },
+  { "kappa", required_argument, NULL, 'k', "set kappa_1 = kappa_2 = kappa", "min:max:inc" },
+  { "lambda", required_argument, NULL, 'l', "set lambda_1 = lambda_2 = lambda", "min:max:inc" },
   { "", 0, NULL, 0, "", "" }
 };
 
@@ -25,8 +27,28 @@ void parseSpecOps()
 {
   for (int isopt = 0; isopt < O815->parsedSpecOps.size(); isopt++)
     switch(O815->parsedSpecOps[isopt].first) {
-    case 'm':
-      O815->paraQ->addRange("mass", O815->parsedSpecOps[isopt].second);
+    case 'r':
+      O815->paraQ->addRange("kappaone", O815->parsedSpecOps[isopt].second);
+      break;
+    case 's':
+      O815->paraQ->addRange("kappatwo", O815->parsedSpecOps[isopt].second);
+      break;
+    case 'b':
+      O815->paraQ->addRange("beta", O815->parsedSpecOps[isopt].second);
+      break;
+    case 'p':
+      O815->paraQ->addRange("lambdaone", O815->parsedSpecOps[isopt].second);
+      break;
+    case 'q':
+      O815->paraQ->addRange("lambdatwo", O815->parsedSpecOps[isopt].second);
+      break;
+    case 'l':
+      O815->paraQ->addRange("lambdaone", O815->parsedSpecOps[isopt].second);
+      O815->paraQ->linkParas("lambdatwo", "lambdaone");
+      break;
+    case 'k':
+      O815->paraQ->addRange("kappaone", O815->parsedSpecOps[isopt].second);
+      O815->paraQ->linkParas("kappatwo", "kappaone");
       break;
     }
 }
@@ -39,11 +61,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, "phi2_hist") == 0 ) {
       *O815->out->log << "MASTER: registered observable: phi2_hist" << endl << flush;
       O815->observables.push_back(new obs_phi2_hist(O815));
@@ -63,7 +85,7 @@ int main (int argc, char *argv[])
   O815->addPara("beta", 1);
   O815->addPara("lambdaone", 1);
   O815->addPara("lambdatwo", 1);
-  O815->addPara("kappaone", 1);
+  O815->addPara("kappaone", 1);  
   O815->addPara("kappatwo", 1);
 
   parseSpecOps();