]> git.treefish.org Git - phys/latlib.git/blob - obs.cpp
added outfile-flush
[phys/latlib.git] / 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 }