Migrated observable n.

diff --git a/u1casc-flux/obs_n.hpp b/u1casc-flux/obs_n.hpp
new file mode 100644 (file)
index 0000000..ffc1148
--- /dev/null
+++ b/u1casc-flux/obs_n.hpp
@@ -0,0 +1,108 @@
+#ifndef OBS_N_HPP
+#define OBS_N_HPP
+
+#include "latlib/o815/o815.h"
+
+#include "latlib/writeout.h"
+
+#include "latlib/obstat.hpp"
+
+#include <iostream>
+
+using namespace std;
+
+class obs_n : public o815::obs {
+
+public:
+  struct obsmem {
+    int lx4sum[2];
+  };
+  obs_n(o815 *_O815);
+  
+private:
+  void _start();
+  void _meas(bool loadedobs, const int& nthmeas);
+  void _finish();
+ 
+  obsmem* OM;
+
+  static void preNSus(vector< vector <double> > *allVals, vector <double> *preCalculated, void *para);
+  static double nSus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para);
+
+  sim *Sim;
+  obstat<double,double> oN[2];
+
+  int lx4sum(int i);
+};
+
+obs_n::obs_n(o815 *_O815) : o815::obs("n", 
+                                           _O815->paraQ->getParaNames() + 
+                                           "none:none_err:nonesus:nonesus_err"
+                                           ":ntwo:ntwo_err:ntwosus:ntwosus_err", 
+                                           _O815, sizeof(obsmem) ) 
+{
+  OM = (obsmem*)obsMem;
+  Sim = (sim*)O815->Sim;
+}
+
+void obs_n::_start() {
+  //*out << "OBS_test: start" << endl;
+};
+
+void obs_n::_meas(bool loadedobs, const int& nthmeas) {
+  if (!loadedobs) {
+    OM->lx4sum[0] = lx4sum(0);
+    OM->lx4sum[1] = lx4sum(1);
+  }
+
+  oN[0].addMeas(OM->lx4sum[0]);
+  oN[1].addMeas(OM->lx4sum[1]);
+};
+
+void obs_n::_finish() {
+  *out << O815->paraQ->getParaVals();
+
+  int compid_n[2], compid_nsus[2];
+
+  for (int iflav=0; iflav<2; iflav++) {
+    compid_n[iflav] = oN[iflav].computeEasy();
+    compid_nsus[iflav] = oN[iflav].computeJack(obs_n::preNSus, obs_n::nSus, &(Sim->lsize4));
+
+    *out << "\t" << oN[iflav].getMean(compid_n[iflav]) / Sim->lsize4 << "\t" << oN[iflav].getErr(compid_n[iflav]) / Sim->lsize4;
+    *out << "\t" << oN[iflav].getMean(compid_nsus[iflav]) << "\t" << oN[iflav].getErr(compid_nsus[iflav]);
+
+    oN[iflav].reset();
+  }
+  
+  *out << endl;
+};
+
+int obs_n::lx4sum(int i)
+{
+  int lx4sum = 0;
+  
+  for( int ix = 0; ix < Sim->lsize4; ix++ )
+    lx4sum += Sim->l[ i*Sim->lsize4*4 + ix*4 + 3 ];
+
+  return lx4sum;
+}
+
+void obs_n::preNSus(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] += (*valIt)[0];
+      (*preCalculated)[1] += pow( (*valIt)[0], 2 );
+    }
+}
+
+double obs_n::nSus(vector <double> *preCalculated, vector<double> *excludedMeas, int nmeas, void *para) {
+  return (( (*preCalculated)[1] - pow((*excludedMeas)[0], 2) ) / (nmeas-1)
+         - pow( ( (*preCalculated)[0] - (*excludedMeas)[0] ) / (nmeas-1), 2 ))
+    / *(int*)para;
+}
+
+#endif

diff --git a/u1casc-flux/u1casc-flux.cpp b/u1casc-flux/u1casc-flux.cpp
index 64869e43321109663342fda305d927d487d6075e..639caa5681544b47809f8d04c2f362b29985f0c3 100644 (file)
--- a/u1casc-flux/u1casc-flux.cpp
+++ b/u1casc-flux/u1casc-flux.cpp
@@ -10,7 +10,7 @@ sim *Sim;
 
 #include "obs_phi2.hpp"
 #include "obs_plaq.hpp"
 
 #include "obs_phi2.hpp"
 #include "obs_plaq.hpp"

-//#include "obs_corr.hpp"
+#include "obs_n.hpp"

 
 o815::comoption specOps[] = {
   { "kappaone", required_argument, NULL, 'r', "set inverse mass kappa_1", "min:max:inc" },
 
 o815::comoption specOps[] = {
   { "kappaone", required_argument, NULL, 'r', "set inverse mass kappa_1", "min:max:inc" },


@@ -82,12 +82,10 @@ void parseLonelyArgs()
       *O815->out->log << "MASTER: registered observable: plaq" << endl << flush;
       O815->observables.push_back(new obs_plaq(O815));
     }
       *O815->out->log << "MASTER: registered observable: plaq" << endl << flush;
       O815->observables.push_back(new obs_plaq(O815));
     }

-    /*
-      else if ( strcmp(*lonit, "corr") == 0 ) {
-      *O815->out->log << "MASTER: registered observable: corr" << endl << flush;
-      O815->observables.push_back(new obs_corr(O815));
-      } 
-    */
+    else if ( strcmp(*lonit, "n") == 0 ) {
+      *O815->out->log << "MASTER: registered observable: n" << endl << flush;
+      O815->observables.push_back(new obs_n(O815));
+    }

   }
 }
 
   }
 }