+#ifndef OBS_HPP
+#define OBS_HPP
+
+#include <list>
+#include <string>
+#include <map>
+#include <math.h>
+
+using namespace std;
+
+template <typename obstype>
+class obs
+{
+ struct result{
+ double val;
+ double err;
+ };
+
+public:
+ void addMeas(obstype val[], int valsize);
+ void addMeas(const obstype& val);
+
+ void computeMean(const string& compid) { mean(compid, &measurements); }
+ void computeJack(const string& compid, double (*func)(list<obstype> *vals, void *para), void *para=NULL);
+
+ double getMean(const string& compid) { return computations[compid].val; }
+ double getErr(const string& compid) { return computations[compid].err; }
+
+private:
+ list<obstype> measurements;
+ map<string,result> computations;
+
+ void mean(const string& compid, list<double> *meas);
+};
+
+template <typename obstype>
+void obs<obstype>::addMeas(obstype val[], int valsize)
+{
+ for(int i=0; i<valsize; i++) measurements.push_back(val[i]);
+}
+
+template <typename obstype>
+void obs<obstype>::addMeas(const obstype& val)
+{
+ measurements.push_back(val);
+}
+
+template <typename obstype>
+void obs<obstype>::mean(const string& compid, list<double> *meas)
+{
+ computations[compid].val = 0;
+ computations[compid].err = 0;
+
+ int nmeas = meas->size();
+
+ for(list<double>::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
+ computations[compid].val += *measIt;
+ computations[compid].val /= nmeas;
+
+ for(list<double>::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
+ computations[compid].err += pow( *measIt - computations[compid].val, 2 );
+ computations[compid].err /= nmeas*(nmeas-1);
+ computations[compid].err = sqrt( computations[compid].err );
+}
+
+template <typename obstype>
+void obs<obstype>::computeJack(const string& compid, double (*func)(list<obstype> *vals, void *para), void *para)
+{
+ int nmeas=measurements.size();
+ double manymeans[nmeas];
+
+ computations[compid].val = 0;
+ computations[compid].err = 0;
+
+ int imeas=0;
+ for(typename list<obstype>::iterator removedIt = measurements.begin(); removedIt != measurements.end(); imeas++)
+ {
+ obstype removed = *removedIt;
+
+ typename list<obstype>::iterator nextAfterIt = removedIt;
+ ++nextAfterIt;
+
+ measurements.erase(removedIt);
+
+ manymeans[imeas] = func(&measurements, para);
+ computations[compid].val += manymeans[imeas];
+
+ measurements.insert(nextAfterIt, removed);
+
+ removedIt = nextAfterIt;
+ }
+ computations[compid].val /= nmeas;
+
+ for(int imean=0; imean<nmeas; imean++)
+ computations[compid].err += pow( computations[compid].val - manymeans[imean], 2 );
+ computations[compid].err *= (double)(nmeas-1)/nmeas;
+ computations[compid].err = sqrt(computations[compid].err);
+}
+
+#endif