obstat.hpp

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