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

Last change on this file since 11969 was 11969, checked in by tbretz, 13 years ago
Added simple function to just return the maxima within a pre-defined range.
File size: 11.6 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
210 /*
211 scalegain *= scaleabs;
212 scalerel *= scaleabs;
213
214 for (size_t ch=0; ch<1440; ch++)
215 {
216 const size_t pos = ch*roi;
217 const int16_t spos = start[ch];
218 if (spos<0)
219 {
220 memset(vec+pos, 0, roi);
221 continue;
222 }
223
224 for (size_t i=0; i<roi; i++)
225 {
226 // Value is relative to trigger
227 // Offset is relative to DRS pipeline
228 // Abs is corresponding index relative to DRS pipeline
229 const size_t rel = pos + i;
230 const size_t abs = pos + (spos+i)%1024;
231
232 const int64_t v =
233 + int64_t(val[rel]) *scaleabs*scalerel
234 - int64_t(trgoff[rel])*scaleabs
235 - int64_t(offset[abs])*scalerel
236 ;
237
238 const int64_t div = int64_t(gain[abs])*scaleabs*scalerel;
239
240 vec[rel] = double(v)*scalegain/div;
241 }
242 }*/
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 {
257 // Value is relative to trigger
258 // Offset is relative to DRS pipeline
259 // Abs is corresponding index relative to DRS pipeline
260 const size_t rel = pos + i;
261 const size_t abs = pos + (spos+i)%1024;
262
263 const int64_t v =
264 + (int64_t(val[rel])*scaleabs-offset[abs])*scalerel
265 - trgoff[rel]
266 ;
267
268 const int64_t div = gain[abs]*scalerel;
269 vec[rel] = double(v)*scalegain/div;
270 }
271 }
272 }
273
274 /*
275 static void Apply(float *vec, const int16_t *val, const int16_t *start, uint32_t roi)
276 {
277 for (size_t ch=0; ch<1440; ch++)
278 {
279 const size_t pos = ch*roi;
280
281 const int16_t spos = start[ch];
282 if (spos<0)
283 {
284 memset(vec+pos, 0, roi);
285 continue;
286 }
287
288 for (size_t i=0; i<roi; i++)
289 vec[pos+i] = float(val[pos+i])/2;
290 }
291 }*/
292
293 pair<vector<double>,vector<double> > GetSampleStats() const
294 {
295 if (fNumEntries==0)
296 return make_pair(vector<double>(),vector<double>());
297
298 vector<double> mean(fSum.size());
299 vector<double> error(fSum.size());
300
301 vector<int64_t>::const_iterator it = fSum.begin();
302 vector<int64_t>::const_iterator i2 = fSum2.begin();
303 vector<double>::iterator im = mean.begin();
304 vector<double>::iterator ie = error.begin();
305
306 while (it!=fSum.end())
307 {
308 *im = /*cnt<fResult.size() ? fResult[cnt] :*/ double(*it)/fNumEntries;
309 *ie = sqrt(double(*i2*int64_t(fNumEntries) - *it * *it))/fNumEntries;
310
311 im++;
312 ie++;
313 it++;
314 i2++;
315 }
316
317
318 /*
319 valarray<double> ...
320
321 mean /= fNumEntries;
322 error = sqrt(error/fNumEntries - mean*mean);
323 */
324
325 return make_pair(mean, error);
326 }
327
328 void GetSampleStats(float *ptr, float scale) const
329 {
330 if (fNumEntries==0)
331 {
332 memset(ptr, 0, sizeof(float)*1024*1440*2);
333 return;
334 }
335
336 vector<int64_t>::const_iterator it = fSum.begin();
337 vector<int64_t>::const_iterator i2 = fSum2.begin();
338
339 while (it!=fSum.end())
340 {
341 *ptr = scale*double(*it)/fNumEntries;
342 *(ptr+1024*1440) = scale*sqrt(double(*i2*int64_t(fNumEntries) - *it * *it))/fNumEntries;
343
344 ptr++;
345 it++;
346 i2++;
347 }
348 }
349
350 static void GetPixelStats(float *ptr, const float *data, uint16_t roi)
351 {
352 if (roi==0)
353 return;
354
355 for (int i=0; i<1440; i++)
356 {
357 const float *vec = data+i*roi;
358
359 int pos = 0;
360 double sum = vec[0];
361 double sum2 = vec[0]*vec[0];
362 for (int j=1; j<roi; j++)
363 {
364 sum += vec[j];
365 sum2 += vec[j]*vec[j];
366
367 if (vec[j]>vec[pos])
368 pos = j;
369 }
370 sum /= roi;
371 sum2 /= roi;
372
373 *(ptr+0*1440+i) = sum;
374 *(ptr+1*1440+i) = sqrt(sum2 - sum * sum);
375 *(ptr+2*1440+i) = vec[pos];
376 *(ptr+3*1440+i) = pos;
377 }
378 }
379
380 static void GetPixelMax(float *ptr, const float *data, uint16_t roi, uint16_t first, uint16_t last)
381 {
382 if (roi==0 || first>=roi || last>=roi || last<first)
383 return;
384
385
386
387 for (int i=0; i<1440; i++)
388 {
389 const float *vec = data+i*roi+first;
390
391 int pos = 0;
392 for (int j=first+1; j<=last; j++)
393 if (vec[j]>vec[pos])
394 pos = j;
395
396 *(ptr+1440+i) = vec[pos];
397 }
398 }
399
400 const vector<int64_t> &GetSum() const { return fSum; }
401
402 uint64_t GetNumEntries() const { return fNumEntries; }
403};
404
Note: See TracBrowser for help on using the repository browser.