obs.hpp

   1 #ifndef OBS_HPP
   2 #define OBS_HPP
   3
   4 #include <vector>
   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)(vector< vector <obstype> > *vals, void *para)=NULL, void *para=NULL);
  21   void computeJack(double (*func)(vector< 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   void reset();
  27
  28 private:
  29   struct result{
  30     double val;
  31     double err;
  32   };
  33
  34   vector< vector<obstype> > measurements;
  35   map<string,result> computations;
  36
  37   void mean(const string& compid, vector< vector <double> > *meas, const int& ival);
  38 };
  39
  40
  41 template <typename obstype>
  42 void obs<obstype>::reset()
  43 {
  44   computations.clear();
  45   measurements.clear();
  46 }
  47
  48 template <typename obstype>
  49 void obs<obstype>::addMeas(const obstype& val)
  50 {
  51   measurements.push_back( vector<obstype>(1,val) );
  52 }
  53
  54 template <typename obstype>
  55 void obs<obstype>::addMeas(obstype val[], int valsize)
  56 {
  57   vector<obstype> tmpvec;
  58   for(int i=0; i<valsize; i++) tmpvec.push_back(val[i]);
  59   measurements.push_back(tmpvec);
  60 }
  61
  62 template <typename obstype>
  63 void obs<obstype>::mean(const string& compid, vector< vector<double> > *meas, const int& ival)
  64 {
  65   computations[compid].val = 0;
  66   computations[compid].err = 0;
  67
  68   int nmeas = meas->size();
  69
  70   for(vector< vector<double> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  71     computations[compid].val += (*measIt)[ival];
  72   computations[compid].val /= nmeas;
  73
  74   for(vector< vector<double> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  75     computations[compid].err += pow( (*measIt)[ival] - computations[compid].val, 2 );
  76   computations[compid].err /= nmeas*(nmeas-1);
  77   computations[compid].err = sqrt( computations[compid].err );
  78 }
  79
  80 template <typename obstype>
  81 void obs<obstype>::computeJack(const string& compid, double (*func)(vector< vector<obstype> > *vals, void *para), void *para)
  82 {
  83   int nmeas=measurements.size();
  84   double manymeans[nmeas];
  85
  86   computations[compid].val = 0;
  87   computations[compid].err = 0;
  88
  89   int imeas=0;
  90   for(typename vector< vector<obstype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++)
  91     {
  92       vector<obstype> removed = *removedIt;
  93
  94       *removedIt = measurements.back();
  95       measurements.pop_back();
  96
  97       manymeans[imeas] = func(&measurements, para);
  98       computations[compid].val += manymeans[imeas];
  99
 100       measurements.push_back( *removedIt );
 101       *removedIt = removed;
 102     }
 103   computations[compid].val /= nmeas;
 104
 105   for(int imean=0; imean<nmeas; imean++)
 106     computations[compid].err += pow( manymeans[imean] - computations[compid].val, 2 );
 107   computations[compid].err *= (double)(nmeas-1)/nmeas;
 108   computations[compid].err = sqrt(computations[compid].err);
 109 }
 110
 111 #endif