obs.cpp

   1 #include "obs.h"
   2
   3 #include <math.h>
   4
   5 using namespace std;
   6
   7 namespace obs
   8 {
   9   vector<meas> mymeas;
  10
  11   double addMeas(int id, double val)
  12   {
  13     meas tmpmeas;
  14     tmpmeas.id = id;
  15     tmpmeas.val.push_back(val);
  16     mymeas.push_back(tmpmeas);
  17     return val;
  18   }
  19
  20   double* addMeas(int id, double *val, int valsize)
  21   {
  22     meas tmpmeas;
  23     tmpmeas.id = id;
  24     for( int ival=0; ival<valsize; ival++ ) tmpmeas.val.push_back( val[ival] );
  25     mymeas.push_back(tmpmeas);
  26     return val;
  27   }
  28
  29   double meanval(int id)
  30   {
  31     double meanval = 0;
  32     int vecsize = mymeas.size();
  33     int ivals=0;
  34     for (int i=0; i<vecsize; i++)
  35       {
  36         if ( mymeas.at(i).id == id )
  37           {
  38             meanval += mymeas.at(i).val[0];
  39             ivals++;
  40           }
  41       }
  42     return meanval/ivals;
  43   }
  44
  45   double meanerr(int id)
  46   {
  47     int ivals = 0;
  48     int vecsize = mymeas.size();
  49     double mymean = meanval(id);
  50     double meanerr = 0;
  51
  52     for (int i=0; i<vecsize; i++)
  53       {
  54         if ( mymeas.at(i).id == id )
  55           {
  56             meanerr += pow(mymeas.at(i).val[0]-mymean, 2);
  57             ivals++;
  58           }
  59       }
  60
  61     meanerr = meanerr / ( ivals*(ivals-1) );
  62     return sqrt(meanerr);
  63   }
  64
  65   obsa normObs(int id)
  66   {
  67     obsa tmp;
  68     tmp.mean = meanval(id);
  69     tmp.err = meanerr(id);
  70     return tmp;
  71   }
  72
  73   obsa jackObs(int id, double (*func)(vector<double> vals))
  74   {
  75     int nmeas = 0;
  76     vector<double> mymeans;
  77     int vecsize = mymeas.size();
  78     obsa jackobs;
  79
  80     jackobs.mean = 0;
  81
  82     for (int ijack=0; ijack < nmeas || nmeas == 0; ijack++)
  83       {
  84
  85         nmeas = 0;
  86         vector<double> myvals;
  87
  88         for (int i=0; i<vecsize; i++)
  89           {
  90             if ( mymeas.at(i).id == id )
  91               {
  92                 nmeas++;
  93                 if ( i == ijack ) continue;
  94                 myvals.push_back( mymeas.at(i).val[0] );
  95               }
  96           }
  97
  98         mymeans.push_back( func(myvals) );
  99         jackobs.mean += mymeans.at(mymeans.size()-1);
 100       }
 101
 102     jackobs.mean /= nmeas;
 103
 104     /* compute jack error */
 105     vecsize = mymeans.size();
 106     jackobs.err = 0;
 107     for (int i=0; i<vecsize; i++)
 108       {
 109         jackobs.err += pow(mymeans.at(i) - jackobs.mean, 2);
 110       }
 111     jackobs.err *= (double)(nmeas-1)/nmeas;
 112     jackobs.err = sqrt(jackobs.err);
 113
 114     return jackobs;
 115   }
 116
 117   obsa jackObs(int id, double (*func)(vector< vector<double> > vals))
 118   {
 119     int nmeas = 0;
 120     vector<double> mymeans;
 121     int vecsize = mymeas.size();
 122     obsa jackobs;
 123
 124     jackobs.mean = 0;
 125
 126     for (int ijack=0; ijack < nmeas || nmeas == 0; ijack++)
 127       {
 128
 129         nmeas = 0;
 130         vector< vector<double> > myvals;
 131
 132         for (int i=0; i<vecsize; i++)
 133           {
 134             if ( mymeas.at(i).id == id )
 135               {
 136                 nmeas++;
 137                 if ( i == ijack ) continue;
 138                 myvals.push_back( mymeas[i].val );
 139               }
 140           }
 141
 142         mymeans.push_back( func(myvals) );
 143         jackobs.mean += mymeans.at(mymeans.size()-1);
 144       }
 145
 146     jackobs.mean /= nmeas;
 147
 148     /* compute jack error */
 149     vecsize = mymeans.size();
 150     jackobs.err = 0;
 151     for (int i=0; i<vecsize; i++)
 152       {
 153         jackobs.err += pow(mymeans.at(i) - jackobs.mean, 2);
 154       }
 155     jackobs.err *= (double)(nmeas-1)/nmeas;
 156     jackobs.err = sqrt(jackobs.err);
 157
 158     return jackobs;
 159   }
 160
 161   void reset()
 162   {
 163     mymeas.clear();
 164   }
 165 }