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

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