source: trunk/Mars/mcore/DrsCalib.h@ 11816

Last change on this file since 11816 was 11816, checked in by tbretz, 13 years ago
File size: 10.5 KB
Line 
1class CalibData
2{
3protected:
4 uint64_t fNumEntries;
5
6 size_t fNumSamples;
7 size_t fNumChannels;
8
9 vector<int64_t> fSum;
10 vector<int64_t> fSum2;
11
12 /*
13 vector<map<int16_t, uint8_t> > fMedian;
14 vector<int16_t> fResult;
15 */
16public:
17 CalibData() : fNumEntries(0), fNumSamples(0), fNumChannels(0) { }
18 void Reset()
19 {
20 fNumEntries = 0;
21 fNumSamples = 0;
22 fNumChannels = 0;
23
24 fSum.clear();
25 fSum2.clear();
26 }
27
28 void InitSize(uint16_t channels, uint16_t samples)
29 {
30 fNumChannels = channels;
31 fNumSamples = samples;
32
33 fSum.resize(samples*channels);
34 fSum2.resize(samples*channels);
35
36 //fMedian.resize(40*9*1024);
37 }
38/*
39 continue;
40
41 if (ch<40*9)
42 {
43 fMedian[abs][val[rel]]++;
44
45 int n= 8;
46 if (fNumEntries>0 && fNumEntries%(n*100)==0)
47 {
48 fResult.resize(40*9*1024);
49
50 uint16_t sum = 0;
51 for (map<int16_t, uint8_t>::const_iterator it=fMedian[abs].begin();
52 it!=fMedian[abs].end(); it++)
53 {
54 map<int16_t, uint8_t>::const_iterator is = it;
55 is++;
56 if (is==fMedian[abs].end())
57 break;
58
59 sum += it->second;
60
61 double med = 0;
62
63 med += int64_t(it->second)*int64_t(it->first);
64 med += int64_t(is->second)*int64_t(is->first);
65 med /= int64_t(it->second)+int64_t(is->second);
66
67 if (sum==n*50)
68 {
69 fResult[abs] = it->first;
70 cout << ch << " MED+(50): " << it->first << endl;
71 }
72 if (sum<n*50 && sum+is->second>n*50)
73 {
74 fResult[abs] = med;
75 cout << ch << " MED-(50): " << med << endl;
76 }
77 }
78 cout << ch << " AVG=" << float(fSum[abs])/fNumEntries << endl;
79 fMedian[abs].clear();
80 }
81 } // -2029 -2012 -2003 -1996 -1990 // -1955
82 */
83
84 void AddRel(const int16_t *val, const int16_t *start)
85 {
86 for (size_t ch=0; ch<fNumChannels; ch++)
87 {
88 const int16_t spos = start[ch];
89 if (spos<0)
90 continue;
91
92 const size_t pos = ch*fNumSamples;
93
94 for (size_t i=0; i<fNumSamples; i++)
95 {
96 // Value is relative to trigger
97 // Abs is corresponding index relative to DRS pipeline
98 const size_t rel = pos + i;
99 const size_t abs = pos + (spos+i)%fNumSamples;
100
101 const int64_t v = val[rel];
102
103 fSum[abs] += v;
104 fSum2[abs] += v*v;
105 }
106 }
107
108 fNumEntries++;
109 }
110
111 void AddRel(const int16_t *val, const int16_t *start,
112 const int32_t *offset, const uint32_t scale)
113 {
114 for (size_t ch=0; ch<fNumChannels; ch++)
115 {
116 const int16_t spos = start[ch];
117 if (spos<0)
118 continue;
119
120 const size_t pos = ch*fNumSamples;
121
122 for (size_t i=0; i<fNumSamples; i++)
123 {
124 // Value is relative to trigger
125 // Offset is relative to DRS pipeline
126 // Abs is corresponding index relative to DRS pipeline
127 const size_t rel = pos + i;
128 const size_t abs = pos + (spos+i)%fNumSamples;
129
130 const int64_t v = int64_t(val[rel])*scale-offset[abs];
131
132 fSum[abs] += v;
133 fSum2[abs] += v*v;
134 }
135 }
136
137 fNumEntries++;
138 }
139
140 void AddAbs(const int16_t *val, const int16_t *start,
141 const int32_t *offset, const uint32_t scale)
142 {
143 for (size_t ch=0; ch<fNumChannels; ch++)
144 {
145 const int16_t spos = start[ch];
146 if (spos<0)
147 continue;
148
149 const size_t pos = ch*fNumSamples;
150
151 for (size_t i=0; i<fNumSamples; i++)
152 {
153 // Value is relative to trigger
154 // Offset is relative to DRS pipeline
155 // Abs is corresponding index relative to DRS pipeline
156 const size_t rel = pos + i;
157 const size_t abs = pos + (spos+i)%fNumSamples;
158
159 const int64_t v = int64_t(val[rel])*scale-offset[abs];
160
161 fSum[rel] += v;
162 fSum2[rel] += v*v;
163 }
164 }
165
166 fNumEntries++;
167 }
168
169 static void Apply(int16_t *val, const int16_t *start, uint32_t roi,
170 const int32_t *offset, const uint32_t scaleabs,
171 const int64_t *gain, const uint32_t scalegain,
172 const int64_t *trgoff, const uint32_t scalerel)
173 {
174 return ;// "FIXME: scalegain exceeds valid range..."
175 for (size_t ch=0; ch<1440; ch++)
176 {
177 const int16_t spos = start[ch];
178 if (spos<0)
179 continue;
180
181 const size_t pos = ch*roi;
182
183 for (size_t i=0; i<roi; i++)
184 {
185 // Value is relative to trigger
186 // Offset is relative to DRS pipeline
187 // Abs is corresponding index relative to DRS pipeline
188 const size_t rel = pos + i;
189 const size_t abs = pos + (spos+i)%1024;
190
191 const int64_t v =
192 + int64_t(val[rel]) *scaleabs*scalerel
193 - int64_t(trgoff[rel])*scaleabs
194 - int64_t(offset[abs])*scalerel
195 ;
196
197 // The gain value is multiplied by fNumOffset and fNumGain
198 const int64_t div = int64_t(gain[abs])*scaleabs*scalerel;
199 val[rel] = (v*scalegain)/div;
200 }
201 }
202 }
203
204 static void Apply(float *vec, const int16_t *val, const int16_t *start, uint32_t roi,
205 const int32_t *offset, const uint32_t scaleabs,
206 const int64_t *gain, const uint64_t scalegain,
207 const int64_t *trgoff, const uint64_t scalerel)
208 {
209 for (size_t ch=0; ch<1440; ch++)
210 {
211 const size_t pos = ch*roi;
212
213 const int16_t spos = start[ch];
214 if (spos<0)
215 {
216 memset(vec+pos, 0, roi);
217 continue;
218 }
219
220 for (size_t i=0; i<roi; i++)
221 {
222 // Value is relative to trigger
223 // Offset is relative to DRS pipeline
224 // Abs is corresponding index relative to DRS pipeline
225 const size_t rel = pos + i;
226 const size_t abs = pos + (spos+i)%1024;
227
228 const int64_t v =
229 + int64_t(val[rel]) *scaleabs*scalerel
230 - int64_t(trgoff[rel])*scaleabs
231 - int64_t(offset[abs])*scalerel
232 ;
233
234 const int64_t div = int64_t(gain[abs])*scaleabs*scalerel;
235
236 vec[rel] = double(v)*scalegain/div;
237 }
238 }
239 }
240
241 /*
242 static void Apply(float *vec, const int16_t *val, const int16_t *start, uint32_t roi)
243 {
244 for (size_t ch=0; ch<1440; ch++)
245 {
246 const size_t pos = ch*roi;
247
248 const int16_t spos = start[ch];
249 if (spos<0)
250 {
251 memset(vec+pos, 0, roi);
252 continue;
253 }
254
255 for (size_t i=0; i<roi; i++)
256 vec[pos+i] = float(val[pos+i])/2;
257 }
258 }*/
259
260 pair<vector<double>,vector<double> > GetSampleStats() const
261 {
262 if (fNumEntries==0)
263 return make_pair(vector<double>(),vector<double>());
264
265 vector<double> mean(fSum.size());
266 vector<double> error(fSum.size());
267
268 vector<int64_t>::const_iterator it = fSum.begin();
269 vector<int64_t>::const_iterator i2 = fSum2.begin();
270 vector<double>::iterator im = mean.begin();
271 vector<double>::iterator ie = error.begin();
272
273 while (it!=fSum.end())
274 {
275 *im = /*cnt<fResult.size() ? fResult[cnt] :*/ double(*it)/fNumEntries;
276 *ie = sqrt(double(*i2*int64_t(fNumEntries) - *it * *it))/fNumEntries;
277
278 im++;
279 ie++;
280 it++;
281 i2++;
282 }
283
284
285 /*
286 valarray<double> ...
287
288 mean /= fNumEntries;
289 error = sqrt(error/fNumEntries - mean*mean);
290 */
291
292 return make_pair(mean, error);
293 }
294
295 void GetSampleStats(float *ptr, float scale) const
296 {
297 if (fNumEntries==0)
298 {
299 memset(ptr, 0, sizeof(float)*1024*1440*2);
300 return;
301 }
302
303 vector<int64_t>::const_iterator it = fSum.begin();
304 vector<int64_t>::const_iterator i2 = fSum2.begin();
305
306 while (it!=fSum.end())
307 {
308 *ptr = scale*double(*it)/fNumEntries;
309 *(ptr+1024*1440) = scale*sqrt(double(*i2*int64_t(fNumEntries) - *it * *it))/fNumEntries;
310
311 ptr++;
312 it++;
313 i2++;
314 }
315 }
316
317 static void GetPixelStats(float *ptr, const float *data, uint16_t roi)
318 {
319 for (int i=0; i<1440; i++)
320 {
321 const float *vec = data+i*roi;
322
323 double max = vec[0];
324 int pos = 0;
325 double sum = vec[0];
326 double sum2 = vec[0]*vec[0];
327 for (int j=1; j<roi; j++)
328 {
329 sum += vec[j];
330 sum2 += vec[j]*vec[j];
331
332 if (vec[j]>max)
333 {
334 max = vec[j];
335 pos = j;
336 }
337 }
338 sum /= roi;
339 sum2 /= roi;
340
341 *(ptr+0*1440+i) = sum;
342 *(ptr+1*1440+i) = sqrt(sum2 - sum * sum);
343 *(ptr+2*1440+i) = max;
344 *(ptr+3*1440+i) = pos;
345 }
346 }
347
348 const vector<int64_t> &GetSum() const { return fSum; }
349
350 uint64_t GetNumEntries() const { return fNumEntries; }
351};
352
353
354/*
355Int_t TH1::MyFill(Int_t x, Double_t w)
356{
357 fEntries++;
358
359 fArrayBin[x] += w;
360
361 fSumw2.fArray[x] += w*w
362 fTsumw += w;
363 fTsumw2 += w*w;
364 fTsumwx += w*x;
365 fTsumwx2 += w*x*x;
366 return bin;
367}
368*/
Note: See TracBrowser for help on using the repository browser.