...

[phys/latlib.git] / obs.hpp

diff --git a/obs.hpp b/obs.hpp
new file mode 100644 (file)
index 0000000..ae544e4
--- /dev/null
+++ b/obs.hpp
@@ -0,0 +1,100 @@
+#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