--- /dev/null
+#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();
+ }
+}