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(void (*preMeasFunc)(vector< vector <meastype> > *allVals, vector <meastype> *preCalculated, void *para),
  21                   restype (*measFunc)(vector <meastype> *preCalculated, vector <meastype> *excludedmeas, int nmeas, void *para), void *para=NULL);
  22
  23   restype getMean(int compid) { return computations[compid].val; }
  24   restype getErr(int compid) { return computations[compid].err; }
  25
  26   void reset();
  27
  28 private:
  29   struct result{
  30     restype val;
  31     restype err;
  32   };
  33
  34   vector< vector<meastype> > measurements;
  35   vector<result> computations;
  36
  37   int mean(vector< vector <meastype> > *meas, const int& ival);
  38 };
  39
  40 template <typename meastype, typename restype>
  41 void obstat<meastype,restype>::reset()
  42 {
  43   computations.clear();
  44   measurements.clear();
  45 }
  46
  47 template <typename meastype, typename restype>
  48 void obstat<meastype,restype>::addMeas(const meastype& val)
  49 {
  50   measurements.push_back( vector<meastype>(1,val) );
  51 }
  52
  53 template <typename meastype, typename restype>
  54 void obstat<meastype,restype>::addMeas(meastype val[], int valsize)
  55 {
  56   vector<meastype> tmpvec;
  57   for (int i=0; i<valsize; i++)
  58     tmpvec.push_back(val[i]);
  59   measurements.push_back(tmpvec);
  60 }
  61
  62 template <typename meastype, typename restype>
  63 int obstat<meastype,restype>::mean(vector< vector<meastype> > *meas, const int& ival)
  64 {
  65   result meanres;
  66
  67   meanres.val = 0;
  68   meanres.err = 0;
  69
  70   int nmeas = meas->size();
  71
  72   for (typename vector< vector<meastype> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  73     meanres.val += (*measIt)[ival];
  74   meanres.val /= nmeas;
  75
  76   for (typename vector< vector<meastype> >::iterator measIt = meas->begin(); measIt != meas->end(); ++measIt)
  77     meanres.err += pow( (*measIt)[ival] - meanres.val, 2 );
  78   meanres.err = sqrt( meanres.err ) / (restype)nmeas;
  79
  80   computations.push_back(meanres);
  81
  82   return computations.size()-1;
  83 }
  84
  85 template <typename meastype, typename restype>
  86 int obstat<meastype,restype>::computeJack(void (*preMeasFunc)(vector< vector <meastype> > *allVals, vector<meastype> *preCalculated, void *para),
  87                                           restype (*measFunc)(vector<meastype> *preCalculated, vector<meastype> *excludedmeas, int nmeas, void *para),
  88                                           void *para) {
  89   int nmeas=measurements.size();
  90   restype *manymeans = new restype[nmeas];
  91   result jackres;
  92   vector<meastype> preCalculated;
  93
  94   jackres.val = 0;
  95   jackres.err = 0;
  96
  97   preMeasFunc(&measurements, &preCalculated, para);
  98
  99   int imeas=0;
 100   for(typename vector< vector<meastype> >::iterator removedIt = measurements.begin(); removedIt != measurements.end(); ++removedIt, imeas++) {
 101     manymeans[imeas] = measFunc(&preCalculated, &(*removedIt), measurements.size(), para);
 102     jackres.val += manymeans[imeas];
 103   }
 104   jackres.val /= nmeas;
 105
 106   for(int imean=0; imean<nmeas; imean++)
 107     jackres.err += pow( manymeans[imean] - jackres.val, 2 );
 108   jackres.err *= (double)(nmeas-1)/nmeas;
 109   jackres.err = sqrt(jackres.err);
 110
 111   computations.push_back(jackres);
 112
 113   delete [] manymeans;
 114
 115   return computations.size()-1;
 116 }
 117
 118 #endif