obs.hpp

   1 #ifndef OBS_HPP
   2 #define OBS_HPP
   3
   4 #include <list>
   5 #include <string>
   6 #include <map>
   7 #include <math.h>
   8
   9 using namespace std;
  10
  11 template <typename obstype>
  12 class obs
  13 {
  14 public:
  15   void addMeas(const obstype& val);
  16   void addMeas(obstype val[], int valsize);
  17
  18   void computeEasy(const string& compid="one", const int& ival=0) { mean(compid, &measurements, ival); }
  19
  20   void computeJack(const string& compid, double (*func)(list< vector <obstype> > *vals, void *para)=NULL, void *para=NULL);
  21   void computeJack(double (*func)(list< vector <obstype> > *vals, void *para)=NULL, void *para=NULL) { computeJack("one",func,para); }
  22
  23   double getMean(const string& compid="one") { return computations[compid].val; }
  24   double getErr(const string& compid="one") { return computations[compid].err; }
  25
  26 private:
  27   struct result{
  28     double val;
  29     double err;
  30   };
  31
  32   list< vector<obstype> > measurements;
  33   map<string,result> computations;
  34
  35   void mean(const string& compid, list< vector <double> > *meas, const int& ival);
  36 };
  37
  38 template <typename obstype>
  39 void obs<obstype>::addMeas(const obstype& val)
  40 {
  41   measurements.push_back( vector<obstype>(1,val) );
  42 }
  43
  44 template <typename obstype>
  45 void obs<obstype>::addMeas(obstype val[], int valsize)
  46 {
  47   vector<obstype> tmpvec;
  48   for(int i=0; i<valsize; i++) tmpvec.push_back(val[i]);
  49   measurements.push_back(tmpvec);
  50 }
  51
  52 template <typename obstype>
  53 void obs<obstype>::mean(const string& compid, list< vector<double> > *meas, const int& ival)
  54 {
  55   computations[compid].val = 0;
  56   computations[compid].err = 0;
  57
  58   int nmeas = meas->size();
  59
  60   for(list< vector<double> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  61     computations[compid].val += (*measIt)[ival];
  62   computations[compid].val /= nmeas;
  63
  64   for(list< vector<double> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  65     computations[compid].err += pow( (*measIt)[ival] - computations[compid].val, 2 );
  66   computations[compid].err /= nmeas*(nmeas-1);
  67   computations[compid].err = sqrt( computations[compid].err );
  68 }
  69
  70 template <typename obstype>
  71 void obs<obstype>::computeJack(const string& compid, double (*func)(list< vector<obstype> > *vals, void *para), void *para)
  72 {
  73   int nmeas=measurements.size();
  74   double manymeans[nmeas];
  75
  76   computations[compid].val = 0;
  77   computations[compid].err = 0;
  78
  79   int imeas=0;
  80   for(typename list< vector<obstype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); imeas++)
  81     {
  82       vector<obstype> removed = *removedIt;
  83
  84       typename list< vector<obstype> >::iterator nextAfterIt = removedIt;
  85       ++nextAfterIt;
  86
  87       measurements.erase(removedIt);
  88
  89       manymeans[imeas] = func(&measurements, para);
  90       computations[compid].val += manymeans[imeas];
  91
  92       measurements.insert(nextAfterIt, removed);
  93
  94       removedIt = nextAfterIt;
  95     }
  96   computations[compid].val /= nmeas;
  97
  98   for(int imean=0; imean<nmeas; imean++)
  99     computations[compid].err += pow( computations[compid].val - manymeans[imean], 2 );
 100   computations[compid].err *= (double)(nmeas-1)/nmeas;
 101   computations[compid].err = sqrt(computations[compid].err);
 102 }
 103
 104 #endif