X-Git-Url: http://git.treefish.org/~alex/phys/latlib.git/blobdiff_plain/7b7c33502031b3d9fd9429768c29c53ed2c5d6f2..9f6e969663a4126813ae8a37d87acdb698ef8e62:/obs.cpp

diff --git a/obs.cpp b/obs.cpp
new file mode 100644
index 0000000..f9a7dbf
--- /dev/null
+++ b/obs.cpp
@@ -0,0 +1,165 @@
+#include "obs.h"
+
+#include <math.h>
+
+using namespace std;
+
+namespace obs
+{
+  vector<meas> mymeas;
+
+  double addMeas(int id, double val)
+  {
+    meas tmpmeas;
+    tmpmeas.id = id;
+    tmpmeas.val.push_back(val);
+    mymeas.push_back(tmpmeas);
+    return val;
+  }
+
+  double* addMeas(int id, double *val, int valsize)
+  {
+    meas tmpmeas;
+    tmpmeas.id = id;
+    for( int ival=0; ival<valsize; ival++ ) tmpmeas.val.push_back( val[ival] );
+    mymeas.push_back(tmpmeas);
+    return val;
+  }
+
+  double meanval(int id)
+  {
+    double meanval = 0;
+    int vecsize = mymeas.size();
+    int ivals=0;
+    for (int i=0; i<vecsize; i++)
+      {
+	if ( mymeas.at(i).id == id )
+	  {
+	    meanval += mymeas.at(i).val[0];
+	    ivals++;
+	  }
+      }
+    return meanval/ivals;
+  }
+
+  double meanerr(int id)
+  {
+    int ivals = 0;
+    int vecsize = mymeas.size();
+    double mymean = meanval(id);
+    double meanerr = 0;
+
+    for (int i=0; i<vecsize; i++)
+      {
+	if ( mymeas.at(i).id == id )
+	  {
+	    meanerr += pow(mymeas.at(i).val[0]-mymean, 2);
+	    ivals++;
+	  }
+      }
+
+    meanerr = meanerr / ( ivals*(ivals-1) );
+    return sqrt(meanerr);
+  }
+
+  obsa normObs(int id)
+  {
+    obsa tmp;
+    tmp.mean = meanval(id);
+    tmp.err = meanerr(id);
+    return tmp;
+  }
+
+  obsa jackObs(int id, double (*func)(vector<double> vals))
+  {
+    int nmeas = 0;
+    vector<double> mymeans;
+    int vecsize = mymeas.size();
+    obsa jackobs;
+
+    jackobs.mean = 0;
+
+    for (int ijack=0; ijack < nmeas || nmeas == 0; ijack++)
+      {
+
+	nmeas = 0;
+	vector<double> myvals;
+
+	for (int i=0; i<vecsize; i++)
+	  {
+	    if ( mymeas.at(i).id == id )
+	      {
+		nmeas++;
+		if ( i == ijack ) continue;
+		myvals.push_back( mymeas.at(i).val[0] );
+	      }
+	  }
+
+	mymeans.push_back( func(myvals) );
+	jackobs.mean += mymeans.at(mymeans.size()-1);
+      }
+
+    jackobs.mean /= nmeas;
+    
+    /* compute jack error */
+    vecsize = mymeans.size();
+    jackobs.err = 0;
+    for (int i=0; i<vecsize; i++)
+      {
+	jackobs.err += pow(mymeans.at(i) - jackobs.mean, 2);
+      }
+    jackobs.err *= (double)(nmeas-1)/nmeas;
+    jackobs.err = sqrt(jackobs.err);
+
+    return jackobs;
+  }
+
+  obsa jackObs(int id, double (*func)(vector< vector<double> > vals))
+  {
+    int nmeas = 0;
+    vector<double> mymeans;
+    int vecsize = mymeas.size();
+    obsa jackobs;
+
+    jackobs.mean = 0;
+
+    for (int ijack=0; ijack < nmeas || nmeas == 0; ijack++)
+      {
+
+	nmeas = 0;
+	vector< vector<double> > myvals;
+
+	for (int i=0; i<vecsize; i++)
+	  {
+	    if ( mymeas.at(i).id == id )
+	      {
+		nmeas++;
+		if ( i == ijack ) continue;
+		myvals.push_back( mymeas[i].val );
+	      }
+	  }
+
+	mymeans.push_back( func(myvals) );
+	jackobs.mean += mymeans.at(mymeans.size()-1);
+      }
+
+    jackobs.mean /= nmeas;
+    
+    /* compute jack error */
+    vecsize = mymeans.size();
+    jackobs.err = 0;
+    for (int i=0; i<vecsize; i++)
+      {
+	jackobs.err += pow(mymeans.at(i) - jackobs.mean, 2);
+      }
+    jackobs.err *= (double)(nmeas-1)/nmeas;
+    jackobs.err = sqrt(jackobs.err);
+
+    return jackobs;
+  }
+
+  void reset()
+  {
+    mymeas.clear();
+  }
+}