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