source: branches/Mars_IncreaseNsb/mcore/DrsCalib.h@ 19904

Last change on this file since 19904 was 18461, checked in by tbretz, 9 years ago
Fixed a problem with the intiatlization of the arrays.
File size: 49.2 KB
Line 
1#ifndef MARS_DrsCalib
2#define MARS_DrsCalib
3
4#include <math.h> // fabs
5#include <errno.h> // errno
6
7#ifndef MARS_fits
8#include "fits.h"
9#endif
10
11#ifndef MARS_ofits
12#include "ofits.h"
13#endif
14
15class DrsCalibrate
16{
17protected:
18 int64_t fNumEntries;
19
20 size_t fNumSamples;
21 size_t fNumChannels;
22
23 std::vector<int64_t> fSum;
24 std::vector<int64_t> fSum2;
25
26public:
27 DrsCalibrate() : fNumEntries(0), fNumSamples(0), fNumChannels(0)
28 {
29 fSum.reserve(1024*1440);
30 fSum2.reserve(1024*1440);
31 }
32
33 void Reset()
34 {
35 fNumEntries = 0;
36 fNumSamples = 0;
37 fNumChannels = 0;
38
39 fSum.clear();
40 fSum2.clear();
41 }
42
43 void InitSize(uint16_t channels, uint16_t samples)
44 {
45 fNumChannels = channels;
46 fNumSamples = samples;
47
48 fSum.assign(samples*channels, 0);
49 fSum2.assign(samples*channels, 0);
50 }
51
52 void AddRel(const int16_t *val, const int16_t *start)
53 {
54 /*
55 for (size_t ch=0; ch<fNumChannels; ch++)
56 {
57 const int16_t &spos = start[ch];
58 if (spos<0)
59 continue;
60
61 const size_t pos = ch*1024;
62 for (size_t i=0; i<1024; i++)
63 {
64 // Value is relative to trigger
65 // Abs is corresponding index relative to DRS pipeline
66 const size_t rel = pos + i;
67 const size_t abs = pos + (spos+i)%1024;
68
69 const int64_t v = val[rel];
70
71 fSum[abs] += v;
72 fSum2[abs] += v*v;
73 }
74 }*/
75
76 // This version is 2.5 times faster because the compilers optimization
77 // is not biased by the evaluation of %1024
78 for (size_t ch=0; ch<fNumChannels; ch++)
79 {
80 const int16_t &spos = start[ch];
81 if (spos<0)
82 continue;
83
84 const size_t pos = ch*1024;
85
86 const int16_t *beg_val = val + pos;
87 int64_t *beg_sum = fSum.data() + pos;
88 int64_t *beg_sum2 = fSum2.data() + pos;
89
90 const int16_t *pval = beg_val; // val[rel]
91 int64_t *psum = beg_sum + spos; // fSum[abs]
92 int64_t *psum2 = beg_sum2 + spos; // fSum2[abs]
93
94 while (psum<beg_sum+1024)
95 {
96 const int64_t v = *pval++;
97
98 *psum++ += v;
99 *psum2++ += v*v;
100 }
101
102 psum = beg_sum;
103 psum2 = beg_sum2;
104
105 while (pval<beg_val+1024)
106 {
107 const int64_t v = *pval++;
108
109 *psum++ += v;
110 *psum2++ += v*v;
111 }
112 }
113
114 fNumEntries++;
115 }
116
117 void AddRel(const int16_t *val, const int16_t *start,
118 const int32_t *offset, const int64_t scale)
119 {
120 /*
121 for (size_t ch=0; ch<fNumChannels; ch++)
122 {
123 const int16_t spos = start[ch];
124 if (spos<0)
125 continue;
126
127 const size_t pos = ch*1024;
128
129 for (size_t i=0; i<fNumSamples; i++)
130 {
131 // Value is relative to trigger
132 // Offset is relative to DRS pipeline
133 // Abs is corresponding index relative to DRS pipeline
134 const size_t rel = pos + i;
135 const size_t abs = pos + (spos+i)%1024;
136
137 const int64_t v = int64_t(val[rel])*scale-offset[abs];
138
139 fSum[abs] += v;
140 fSum2[abs] += v*v;
141 }
142 }*/
143
144 // This version is 2.5 times faster because the compilers optimization
145 // is not biased by the evaluation of %1024
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*1024;
153
154 const int16_t *beg_val = val + pos;
155 const int32_t *beg_offset = offset + pos;
156 int64_t *beg_sum = fSum.data() + pos;
157 int64_t *beg_sum2 = fSum2.data() + pos;
158
159
160 const int16_t *pval = beg_val; // val[rel]
161 const int32_t *poffset = beg_offset + spos; // offset[abs]
162 int64_t *psum = beg_sum + spos; // fSum[abs]
163 int64_t *psum2 = beg_sum2 + spos; // fSum2[abs]
164
165 while (psum<beg_sum+1024)
166 {
167 const int64_t v = int64_t(*pval++)*scale - *poffset++;
168
169 *psum++ += v;
170 *psum2++ += v*v;
171 }
172
173 psum = beg_sum;
174 psum2 = beg_sum2;
175 poffset = beg_offset;
176
177 while (pval<beg_val+1024)
178 {
179 const int64_t v = int64_t(*pval++)*scale - *poffset++;
180
181 *psum++ += v;
182 *psum2++ += v*v;
183 }
184 }
185
186 fNumEntries++;
187 }
188 /*
189 void AddAbs(const int16_t *val, const int16_t *start,
190 const int32_t *offset, const uint32_t scale)
191 {
192 for (size_t ch=0; ch<fNumChannels; ch++)
193 {
194 const int16_t spos = start[ch];
195 if (spos<0)
196 continue;
197
198 const size_t pos = ch*fNumSamples;
199
200 for (size_t i=0; i<fNumSamples; i++)
201 {
202 // Value is relative to trigger
203 // Offset is relative to DRS pipeline
204 // Abs is corresponding index relative to DRS pipeline
205 const size_t rel = pos + i;
206 const size_t abs = pos + (spos+i)%1024;
207
208 const int64_t v = int64_t(val[rel])*scale-offset[abs];
209
210 fSum[rel] += v;
211 fSum2[rel] += v*v;
212 }
213 }
214
215 fNumEntries++;
216 }*/
217
218 void AddAbs(const int16_t *val, const int16_t *start,
219 const int32_t *offset, const int64_t scale)
220 {
221 /*
222 // 1440 without tm, 1600 with tm
223 for (size_t ch=0; ch<fNumChannels; ch++)
224 {
225 const int16_t spos = start[ch];
226 if (spos<0)
227 continue;
228
229 const size_t pos = ch*fNumSamples;
230 const size_t drs = ch>1439 ? ((ch-1440)*9+8)*1024 : ch*1024;
231
232 for (size_t i=0; i<fNumSamples; i++)
233 {
234 // Value is relative to trigger
235 // Offset is relative to DRS pipeline
236 // Abs is corresponding index relative to DRS pipeline
237 const size_t rel = pos + i;
238 const size_t abs = drs + (spos+i)%1024;
239
240 const int64_t v = int64_t(val[rel])*scale-offset[abs];
241
242 fSum[rel] += v;
243 fSum2[rel] += v*v;
244 }
245 }*/
246
247 // This version is 1.5 times faster because the compilers optimization
248 // is not biased by the evaluation of %1024
249 for (size_t ch=0; ch<fNumChannels; ch++)
250 {
251 const int16_t &spos = start[ch];
252 if (spos<0)
253 continue;
254
255 const size_t pos = ch*fNumSamples;
256
257 const int32_t *beg_offset = offset + ch*1024;
258 const int16_t *beg_val = val + pos;
259 int64_t *beg_sum = fSum.data() + pos;
260 int64_t *beg_sum2 = fSum2.data() + pos;
261
262
263 const int16_t *pval = beg_val; // val[rel]
264 const int32_t *poffset = beg_offset + spos; // offset[abs]
265 int64_t *psum = beg_sum; // fSum[rel]
266 int64_t *psum2 = beg_sum2; // fSum2[rel]
267
268 if (spos+fNumSamples>1024)
269 {
270 while (poffset<beg_offset+1024)
271 {
272 const int64_t v = int64_t(*pval++)*scale - *poffset++;
273
274 *psum++ += v;
275 *psum2++ += v*v;
276 }
277
278 poffset = beg_offset;
279 }
280
281 while (psum<beg_sum+fNumSamples)
282 {
283 const int64_t v = int64_t(*pval++)*scale - *poffset++;
284
285 *psum++ += v;
286 *psum2++ += v*v;
287 }
288 }
289
290 fNumEntries++;
291 }
292
293
294 static void ApplyCh(float *vec, const int16_t *val, int16_t start, uint32_t roi,
295 const int32_t *offset, const int64_t scaleabs,
296 const int64_t *gain, const int64_t scalegain)
297 {
298 if (start<0)
299 {
300 memset(vec, 0, roi);
301 return;
302 }
303 /*
304 for (size_t i=0; i<roi; i++)
305 {
306 // Value is relative to trigger
307 // Offset is relative to DRS pipeline
308 // Abs is corresponding index relative to DRS pipeline
309 const size_t abs = (start+i)%1024;
310
311 const int64_t v =
312 + int64_t(val[i])*scaleabs-offset[abs]
313 ;
314
315 const int64_t div = gain[abs];
316 vec[i] = div==0 ? 0 : double(v)*scalegain/div;
317 }
318 */
319
320 // This version is faster because the compilers optimization
321 // is not biased by the evaluation of %1024
322 // (Here we are dominated by numerics... improvement ~10%)
323 const int32_t *poffset = offset + start; // offset[abs]
324 const int64_t *pgain = gain + start; // gain[abs]
325 const int16_t *pval = val; // val[rel]
326 float *pvec = vec; // vec[rel]
327
328 if (start+roi>1024)
329 {
330 while (poffset<offset+1024)
331 {
332 const int64_t v =
333 + int64_t(*pval++)*scaleabs - *poffset++
334 ;
335
336 *pvec++ = *pgain==0 ? 0 : double(v)*scalegain / *pgain;
337
338 pgain++;
339 }
340
341 poffset = offset;
342 pgain = gain;
343 }
344
345 while (pvec<vec+roi)
346 {
347 const int64_t v =
348 + int64_t(*pval++)*scaleabs - *poffset++
349 ;
350
351 *pvec++ = *pgain==0 ? 0 : double(v)*scalegain / *pgain;
352
353 pgain++;
354 }
355 }
356
357 static void ApplyCh(float *vec, const int16_t *val, int16_t start, uint32_t roi,
358 const int32_t *offset, const int64_t scaleabs,
359 const int64_t *gain, const int64_t scalegain,
360 const int64_t *trgoff, const int64_t scalerel)
361 {
362 if (start<0)
363 {
364 memset(vec, 0, roi);
365 return;
366 }
367 /*
368 for (size_t i=0; i<roi; i++)
369 {
370 // Value is relative to trigger
371 // Offset is relative to DRS pipeline
372 // Abs is corresponding index relative to DRS pipeline
373 const size_t abs = (start+i)%1024;
374
375 const int64_t v =
376 + (int64_t(val[i])*scaleabs-offset[abs])*scalerel
377 - trgoff[i]
378 ;
379
380 const int64_t div = gain[abs]*scalerel;
381 vec[i] = div==0 ? 0 : double(v)*scalegain/div;
382 }
383 */
384 // (Here we are dominated by numerics... improvement ~10%)
385 const int32_t *poffset = offset + start; // offset[abs]
386 const int64_t *pgain = gain + start; // gain[abs]
387 const int16_t *pval = val; // val[rel]
388 const int64_t *ptrgoff = trgoff; // trgoff[rel]
389 float *pvec = vec; // vec[rel]
390
391 if (start+roi>1024)
392 {
393 while (poffset<offset+1024)
394 {
395 const int64_t v =
396 + (int64_t(*pval++)*scaleabs - *poffset++)*scalerel
397 - *ptrgoff++;
398 ;
399
400 const int64_t div = *pgain * scalerel;
401 *pvec++ = div==0 ? 0 : double(v)*scalegain / div;
402
403 pgain++;
404 }
405
406 poffset = offset;
407 pgain = gain;
408 }
409
410 while (pvec<vec+roi)
411 {
412 const int64_t v =
413 + (int64_t(*pval++)*scaleabs - *poffset++)*scalerel
414 - *ptrgoff++;
415 ;
416
417 const int64_t div = *pgain * scalerel;
418 *pvec++ = div==0 ? 0 : double(v)*scalegain / div;
419
420 pgain++;
421 }
422 }
423
424 static double FindStep(const size_t ch0, const float *vec, int16_t roi, const int16_t pos, const uint16_t *map=NULL)
425 {
426 // We have about 1% of all cases which are not ahndled here,
427 // because the baseline jumps up just before the readout window
428 // and down just after it. In this cases we could determine the jump
429 // from the board time difference or throw the event away.
430 if (pos==0 || pos>=roi)
431 return 0;
432
433 double step = 0; // before
434 double rms = 0; // before
435 int cnt = 0;
436
437 // Exclude TM channel
438 for (int p=0; p<8; p++)
439 {
440 const size_t hw = ch0+p;
441 const size_t sw = (map?map[hw]:hw)*roi + pos;
442
443 const double diff = vec[sw]-vec[sw-1];
444
445 step += diff;
446 rms += (vec[sw]-vec[sw-1])*(vec[sw]-vec[sw-1]);
447
448 cnt++;
449 }
450
451 return cnt==0 ? 0 : step/cnt;
452 }
453
454 static void SubtractStep(const size_t ch0, const double avg, float *vec, int16_t roi, int32_t pos, const uint16_t *map=NULL)
455 {
456 if (pos==0 || pos>=roi)
457 return;
458
459 const int begin = avg>0 ? pos : 0;
460 const int end = avg>0 ? roi : pos;
461
462 const double sub = fabs(avg);
463
464 for (int p=0; p<9; p++)
465 {
466 for (int j=begin; j<end; j++)
467 {
468 const size_t hw = ch0+p;
469 const size_t sw = (map?map[hw]:hw)*roi + j;
470
471 vec[sw] -= sub;
472 }
473 }
474 }
475
476 struct Step
477 {
478 Step() : avg(0), rms(0), pos(0), cnt(0) { }
479 double avg;
480 double rms;
481 double pos;
482 uint16_t cnt;
483
484 static bool sort(const Step &s, const Step &r) { return s.avg<r.avg; }
485 };
486
487 static Step AverageSteps(const std::vector<Step>::iterator beg, const std::vector<Step>::iterator end)
488 {
489 Step rc;
490 for (auto it=beg; it!=end; it++)
491 {
492 rc.pos += it->pos;
493 rc.avg += it->avg;
494 rc.rms += it->avg*it->avg;
495 }
496
497 rc.cnt = end-beg;
498
499 rc.pos /= rc.cnt;
500 rc.avg /= rc.cnt;
501 rc.rms /= rc.cnt;
502 rc.rms -= rc.avg*rc.avg;
503 rc.rms = rc.rms<0 ? 0 : sqrt(rc.rms);
504
505 return rc;
506 }
507
508
509 static Step CorrectStep(float *vec, uint16_t nch, uint16_t roi,
510 const int16_t *prev, const int16_t *start,
511 const int16_t offset, const uint16_t *map=NULL)
512 {
513
514 std::vector<Step> list;
515 list.reserve(nch);
516
517 // Fill steps into array
518 // Exclude broken pixels?
519 // Remove maximum and minimum patches (4max and 4min)?
520 for (size_t ch=0; ch<nch; ch += 9)
521 {
522 if (prev[ch]<0 || start[ch]<0)
523 continue;
524
525 const int16_t dist = (prev[ch]-start[ch]+1024+offset)%1024;
526 const double step = FindStep(ch, vec, roi, dist, map);
527 if (step==0)
528 continue;
529
530 Step rc;
531 rc.pos = dist;
532 rc.avg = step;
533 list.push_back(rc);
534 }
535
536 if (list.empty())
537 return Step();
538
539 Step rc = AverageSteps(list.begin(), list.begin()+list.size());;
540
541 if (rc.avg==0)
542 return Step();
543
544 // std::cout << " A0=" << rc.avg << " rms=" << rc.rms << std::endl;
545 if (rc.rms>5)
546 {
547 sort(list.begin(), list.end(), Step::sort);
548
549 //for (auto it=list.begin(); it!=list.end(); it++)
550 // std::cout << " " << it->avg << std::endl;
551
552 const size_t skip = list.size()/10;
553 rc = AverageSteps(list.begin()+skip, list.begin()+list.size()-skip);
554
555 // std::cout << " A1=" << rc.avg << " rms=" << rc.rms << std::endl;
556 }
557
558 for (size_t ch=0; ch<nch; ch += 9)
559 {
560 const int16_t dist = (prev[ch]-start[ch]+1024+offset)%1024;
561 SubtractStep(ch, rc.avg, vec, roi, dist, map);
562 }
563
564 return rc;
565 }
566
567 static void RemoveSpikes(float *p, uint32_t roi)
568 {
569 if (roi<4)
570 return;
571
572 for (size_t i=1; i<roi-2; i++)
573 {
574 if (p[i]-p[i-1]>25 && p[i]-p[i+1]>25)
575 {
576 p[i] = (p[i-1]+p[i+1])/2;
577 }
578
579 if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22)
580 {
581 p[i] = (p[i-1]+p[i+2])/2;
582 p[i+1] = p[i];
583 }
584 }
585 }
586
587 static void RemoveSpikes2(float *p, uint32_t roi)//from Werner
588 {
589 if (roi<4)
590 return;
591
592 std::vector<float> Ameas(p, p+roi);
593
594 std::vector<float> diff(roi);
595 for (size_t i=1; i<roi-1; i++)
596 diff[i] = (p[i-1] + p[i+1])/2 - p[i];
597
598 //std::vector<float> N1mean(roi);
599 //for (size_t i=1; i<roi-1; i++)
600 // N1mean[i] = (p[i-1] + p[i+1])/2;
601
602 const float fract = 0.8;
603
604 for (size_t i=0; i<roi-3; i++)
605 {
606 if (diff[i]<5)
607 continue;
608
609 if (Ameas[i+2] - (Ameas[i] + Ameas[i+3])/2 > 10)
610 {
611 p[i+1]= (Ameas[i+3] - Ameas[i])/3 + Ameas[i];
612 p[i+2]= 2*(Ameas[i+3] - Ameas[i])/3 + Ameas[i];
613
614 i += 3;
615
616 continue;
617 }
618
619 if ( (diff[i+1]<-diff[i]*fract*2) && (diff[i+2]>10) )
620 {
621 p[i+1] = (Ameas[i]+Ameas[i+2])/2;
622 diff[i+2] = (p[i+1] + Ameas[i+3])/2 - Ameas[i+2];
623
624 i += 2;
625 }
626
627 // const float x = Ameas[i] - N1mean[i];
628 // if (x > -5.)
629 // continue;
630
631 // if (Ameas[i+2] - (Ameas[i] + Ameas[i+3])/2. > 10.)
632 // {
633 // p[i+1]= (Ameas[i+3] - Ameas[i])/3 + Ameas[i];
634 // p[i+2]= 2*(Ameas[i+3] - Ameas[i])/3 + Ameas[i];
635 // i += 3;
636 // continue;
637 // }
638
639 // const float xp = Ameas[i+1] - N1mean[i+1];
640 // const float xpp = Ameas[i+2] - N1mean[i+2];
641
642 // if ( (xp > -2.*x*fract) && (xpp < -10.) )
643 // {
644 // p[i+1] = N1mean[i+1];
645 // N1mean[i+2] = Ameas[i+1] - Ameas[i+3]/2;
646 //
647 // i += 2;
648 // }
649 }
650 }
651
652 static void RemoveSpikes3(float *vec, uint32_t roi)//from Werner
653 {
654 if (roi<4)
655 return;
656
657 const float SingleCandidateTHR = -10.;
658 const float DoubleCandidateTHR = -5.;
659
660 const std::vector<float> src(vec, vec+roi);
661
662 std::vector<float> diff(roi);
663 for (size_t i=1; i<roi-1; i++)
664 diff[i] = src[i] - (src[i-1] + src[i+1])/2;
665
666 // find the spike and replace it by mean value of neighbors
667 for (unsigned int i=1; i<roi-3; i++)
668 {
669 // Speed up (no leading edge)
670 if (diff[i]>=DoubleCandidateTHR)
671 continue;
672
673 //bool checkDouble = false;
674
675 // a single spike candidate
676 if (diff[i]<SingleCandidateTHR)
677 {
678 // check consistency with a single channel spike
679 if (diff[i+1] > -1.6*diff[i])
680 {
681 vec[i+1] = (src[i] + src[i+2]) / 2;
682
683 i += 2;
684
685 /*** NEW ***/
686 continue;
687 /*** NEW ***/
688 }
689 /*
690 else
691 {
692 // do nothing - not really a single spike,
693 // but check if it is a double
694 checkDouble = true;
695 }*/
696 }
697
698 // a double spike candidate
699 //if (diff[i]>DoubleCandidateTHR || checkDouble == 1)
700 {
701 // check the consistency with a double spike
702 if ((diff[i+1] > -DoubleCandidateTHR) &&
703 (diff[i+2] > -DoubleCandidateTHR))
704 {
705 vec[i+1] = (src[i+3] - src[i])/3 + src[i];
706 vec[i+2] = 2*(src[i+3] - src[i])/3 + src[i];
707
708 //vec[i] = (src[i-1] + src[i+2]) / 2.;
709 //vec[i+1] = (src[i-1] + src[i+2]) / 2.;
710
711 //do not care about the next sample it was the spike
712 i += 3;
713 }
714 }
715 }
716 }
717
718 static void RemoveSpikes4(float *ptr, uint32_t roi)
719 {
720 if (roi<7)
721 return;
722
723 for (uint32_t i=0; i<roi-6; i++)
724 {
725 double d10, d21, d32, d43, d54;
726
727 // ============================================
728 d43 = ptr[i+4]-ptr[i+3];
729 d54 = ptr[i+5]-ptr[i+4];
730
731 if ((d43>35 && -d54>35) || (d43<-35 && -d54<-35))
732 {
733 ptr[i+4] = (ptr[i+3]+ptr[i+5])/2;
734 }
735
736 // ============================================
737 d32 = ptr[i+3]-ptr[i+2];
738 d54 = ptr[i+5]-ptr[i+4];
739
740 if ((d32>9 && -d54>13 && d32-d54>31)/* || (d32<-13 && -d54<-13 && d32+d54<-63)*/)
741 {
742 double avg0 = (ptr[i+2]+ptr[i+5])/2;
743 double avg1 = (ptr[i+3]+ptr[i+4])/2;
744
745 ptr[i+3] = ptr[i+3] - avg1+avg0;
746 ptr[i+4] = ptr[i+4] - avg1+avg0;
747 }
748
749 // ============================================
750 d21 = ptr[i+2]-ptr[i+1];
751 d54 = ptr[i+5]-ptr[i+4];
752
753 if (d21>15 && -d54>17)
754 {
755 double avg0 = (ptr[i+1]+ptr[i+5])/2;
756 double avg1 = (ptr[i+2]+ptr[i+3]+ptr[i+4])/3;
757
758 ptr[i+2] = ptr[i+2] - avg1+avg0;
759 ptr[i+3] = ptr[i+3] - avg1+avg0;
760 ptr[i+4] = ptr[i+4] - avg1+avg0;
761 }
762
763 // ============================================
764 d10 = ptr[i+1]-ptr[i];
765 d54 = ptr[i+5]-ptr[i+4];
766
767 if (d10>18 && -d54>20)
768 {
769 double avg0 = (ptr[i]+ptr[i+5])/2;
770 double avg1 = (ptr[i+1]+ptr[i+2]+ptr[i+3]+ptr[i+4])/4;
771
772 ptr[i+1] = ptr[i+1] - avg1+avg0;
773 ptr[i+2] = ptr[i+2] - avg1+avg0;
774 ptr[i+3] = ptr[i+3] - avg1+avg0;
775 ptr[i+4] = ptr[i+4] - avg1+avg0;
776 }
777 }
778 }
779
780 static void SlidingAverage(float *const vec, const uint32_t roi, const uint16_t w)
781 {
782 if (w==0 || w>roi)
783 return;
784
785 for (float *pix=vec; pix<vec+1440*roi; pix += roi)
786 {
787 for (float *ptr=pix; ptr<pix+roi-w; ptr++)
788 {
789 for (float *p=ptr+1; p<ptr+w; p++)
790 *ptr += *p;
791 *ptr /= w;
792 }
793 }
794 }
795
796 std::pair<std::vector<double>,std::vector<double> > GetSampleStats() const
797 {
798 if (fNumEntries==0)
799 return make_pair(std::vector<double>(),std::vector<double>());
800
801 std::vector<double> mean(fSum.size());
802 std::vector<double> error(fSum.size());
803
804 std::vector<int64_t>::const_iterator it = fSum.begin();
805 std::vector<int64_t>::const_iterator i2 = fSum2.begin();
806 std::vector<double>::iterator im = mean.begin();
807 std::vector<double>::iterator ie = error.begin();
808
809 while (it!=fSum.end())
810 {
811 *im = /*cnt<fResult.size() ? fResult[cnt] :*/ double(*it)/fNumEntries;
812 *ie = sqrt(double(*i2*int64_t(fNumEntries) - *it * *it))/fNumEntries;
813
814 im++;
815 ie++;
816 it++;
817 i2++;
818 }
819
820
821 /*
822 valarray<double> ...
823
824 mean /= fNumEntries;
825 error = sqrt(error/fNumEntries - mean*mean);
826 */
827
828 return make_pair(mean, error);
829 }
830
831 void GetSampleStats(float *ptr, float scale) const
832 {
833 const size_t sz = fNumSamples*fNumChannels;
834
835 if (fNumEntries==0)
836 {
837 memset(ptr, 0, sizeof(float)*sz*2);
838 return;
839 }
840
841 std::vector<int64_t>::const_iterator it = fSum.begin();
842 std::vector<int64_t>::const_iterator i2 = fSum2.begin();
843
844 while (it!=fSum.end())
845 {
846 *ptr = scale*double(*it)/fNumEntries;
847 *(ptr+sz) = scale*sqrt(double(*i2*fNumEntries - *it * *it))/fNumEntries;
848
849 ptr++;
850 it++;
851 i2++;
852 }
853 }
854
855 static double GetPixelStats(float *ptr, const float *data, uint16_t roi, uint16_t begskip=0, uint16_t endskip=0)
856 {
857 if (roi==0)
858 return -1;
859
860 // Skip first 10 samples
861 const uint beg = roi>begskip ? begskip : 0;
862 const uint end = roi-beg>endskip ? roi-endskip : roi;
863 const uint len = end-beg;
864
865 double max = 0;
866 double patch = 0;
867 for (uint i=0; i<1440; i++)
868 {
869 const float *vec = data+i*roi;
870
871 uint pos = beg;
872 double sum = vec[beg];
873 double sum2 = vec[beg]*vec[beg];
874
875 for (uint j=beg+1; j<end; j++)
876 {
877 sum += vec[j];
878 sum2 += vec[j]*vec[j];
879
880 if (vec[j]>vec[pos])
881 pos = j;
882 }
883 sum /= len;
884 sum2 /= len;
885 sum2 -= sum*sum;
886
887 if (i%9!=8)
888 patch += vec[pos];
889 else
890 {
891 if (patch > max)
892 max = patch;
893 patch = 0;
894 }
895
896 *(ptr+0*1440+i) = sum;
897 *(ptr+1*1440+i) = sum2<0 ? 0 : sqrt(sum2);
898 *(ptr+2*1440+i) = vec[pos];
899 *(ptr+3*1440+i) = pos;
900 }
901
902 return max/8;
903 }
904
905 static void GetPixelMax(float *max, const float *data, uint16_t roi, int32_t first, int32_t last)
906 {
907 if (roi==0 || first<0 || last<0 || first>=roi || last>=roi || last<first)
908 return;
909
910 for (int i=0; i<1440; i++)
911 {
912 const float *beg = data+i*roi+first;
913 const float *end = data+i*roi+last;
914
915 const float *pmax = beg;
916
917 for (const float *ptr=beg+1; ptr<=end; ptr++)
918 if (*ptr>*pmax)
919 pmax = ptr;
920
921 max[i] = *pmax;
922 }
923 }
924
925 const std::vector<int64_t> &GetSum() const { return fSum; }
926
927 int64_t GetNumEntries() const { return fNumEntries; }
928};
929
930class DrsCalibrateTime
931{
932public:
933 int64_t fNumEntries;
934
935 size_t fNumSamples;
936 size_t fNumChannels;
937
938 std::vector<std::pair<double, double>> fStat;
939
940public:
941 DrsCalibrateTime() : fNumEntries(0), fNumSamples(0), fNumChannels(0)
942 {
943 InitSize(160, 1024);
944 }
945
946 DrsCalibrateTime(const DrsCalibrateTime &p) : fNumEntries(p.fNumEntries), fNumSamples(p.fNumSamples), fNumChannels(p.fNumChannels), fStat(p.fStat)
947 {
948 }
949 virtual ~DrsCalibrateTime()
950 {
951 }
952
953 double Sum(uint32_t i) const { return fStat[i].first; }
954 double W(uint32_t i) const { return fStat[i].second; }
955
956 virtual void InitSize(uint16_t channels, uint16_t samples)
957 {
958 fNumChannels = channels;
959 fNumSamples = samples;
960
961 fNumEntries = 0;
962
963 fStat.clear();
964
965 fStat.resize(samples*channels);
966 }
967
968 void Reset()
969 {
970 for (auto it=fStat.begin(); it!=fStat.end(); it++)
971 {
972 it->first = 0;
973 it->second = 0;
974 }
975 }
976
977 void AddT(const float *val, const int16_t *start, signed char edge=0)
978 {
979 if (fNumSamples!=1024 || fNumChannels!=160)
980 return;
981
982 // Rising or falling edge detection has the advantage that
983 // we are much less sensitive to baseline shifts
984
985 for (size_t ch=0; ch<160; ch++)
986 {
987 const size_t tm = ch*9+8;
988
989 const int16_t spos = start[tm];
990 if (spos<0)
991 continue;
992
993 const size_t pos = ch*1024;
994
995 double p_prev = 0;
996 int32_t i_prev = -1;
997
998 for (size_t i=0; i<1024-1; i++)
999 {
1000 const size_t rel = tm*1024 + i;
1001
1002 const float &v0 = val[rel]; //-avg;
1003 const float &v1 = val[rel+1];//-avg;
1004
1005 // If edge is positive ignore all falling edges
1006 if (edge>0 && v0>0)
1007 continue;
1008
1009 // If edge is negative ignore all falling edges
1010 if (edge<0 && v0<0)
1011 continue;
1012
1013 // Check if there is a zero crossing
1014 if ((v0<0 && v1<0) || (v0>0 && v1>0))
1015 continue;
1016
1017 // Calculate the position p of the zero-crossing
1018 // within the interval [rel, rel+1] relative to rel
1019 // by linear interpolation.
1020 const double p = v0==v1 ? 0.5 : v0/(v0-v1);
1021
1022 // If this was at least the second zero-crossing detected
1023 if (i_prev>=0)
1024 {
1025 // Calculate the distance l between the
1026 // current and the last zero-crossing
1027 const double l = i+p - (i_prev+p_prev);
1028
1029 // By summation, the average length of each
1030 // cell is calculated. For the first and last
1031 // fraction of a cell, the fraction is applied
1032 // as a weight.
1033 const double w0 = 1-p_prev;
1034 fStat[pos+(spos+i_prev)%1024].first += w0*l;
1035 fStat[pos+(spos+i_prev)%1024].second += w0;
1036
1037 for (size_t k=i_prev+1; k<i; k++)
1038 {
1039 fStat[pos+(spos+k)%1024].first += l;
1040 fStat[pos+(spos+k)%1024].second += 1;
1041 }
1042
1043 const double w1 = p;
1044 fStat[pos+(spos+i)%1024].first += w1*l;
1045 fStat[pos+(spos+i)%1024].second += w1;
1046 }
1047
1048 // Remember this zero-crossing position
1049 p_prev = p;
1050 i_prev = i;
1051 }
1052 }
1053 fNumEntries++;
1054 }
1055
1056 void FillEmptyBins()
1057 {
1058 for (int ch=0; ch<160; ch++)
1059 {
1060 const auto beg = fStat.begin() + ch*1024;
1061 const auto end = beg + 1024;
1062
1063 double avg = 0;
1064 uint32_t num = 0;
1065 for (auto it=beg; it!=end; it++)
1066 {
1067 if (it->second<fNumEntries-0.5)
1068 continue;
1069
1070 avg += it->first / it->second;
1071 num++;
1072 }
1073 avg /= num;
1074
1075 for (auto it=beg; it!=end; it++)
1076 {
1077 if (it->second>=fNumEntries-0.5)
1078 continue;
1079
1080 // {
1081 // result[i+1].first = *is2;
1082 // result[i+1].second = *iw2;
1083 // }
1084 // else
1085 // {
1086 it->first = avg*fNumEntries;
1087 it->second = fNumEntries;
1088 // }
1089 }
1090 }
1091 }
1092
1093 DrsCalibrateTime GetComplete() const
1094 {
1095 DrsCalibrateTime rc(*this);
1096 rc.FillEmptyBins();
1097 return rc;
1098 }
1099
1100 void CalcResult()
1101 {
1102 for (int ch=0; ch<160; ch++)
1103 {
1104 const auto beg = fStat.begin() + ch*1024;
1105 const auto end = beg + 1024;
1106
1107 // First calculate the average length s of a single
1108 // zero-crossing interval in the whole range [0;1023]
1109 // (which is identical to the/ wavelength of the
1110 // calibration signal)
1111 double s = 0;
1112 double w = 0;
1113 for (auto it=beg; it!=end; it++)
1114 {
1115 s += it->first;
1116 w += it->second;
1117 }
1118 s /= w;
1119
1120 // Dividing the average length s of the zero-crossing
1121 // interval in the range [0;1023] by the average length
1122 // in the interval [0;n] yields the relative size of
1123 // the interval in the range [0;n].
1124 //
1125 // Example:
1126 // Average [0;1023]: 10.00 (global interval size in samples)
1127 // Average [0;512]: 8.00 (local interval size in samples)
1128 //
1129 // Globally, on average one interval is sampled by 10 samples.
1130 // In the sub-range [0;512] one interval is sampled on average
1131 // by 8 samples.
1132 // That means that the interval contains 64 periods, while
1133 // in the ideal case (each sample has the same length), it
1134 // should contain 51.2 periods.
1135 // So, the sampling position 512 corresponds to a time 640,
1136 // while in the ideal case with equally spaces samples,
1137 // it would correspond to a time 512.
1138 //
1139 // The offset (defined as 'ideal - real') is then calculated
1140 // as 512*(1-10/8) = -128, so that the time is calculated as
1141 // 'sampling position minus offset'
1142 //
1143 double sumw = 0;
1144 double sumv = 0;
1145 int n = 0;
1146
1147 // Sums about many values are numerically less stable than
1148 // just sums over less. So we do the exercise from both sides.
1149 // First from the left
1150 for (auto it=beg; it!=end-512; it++, n++)
1151 {
1152 const double valv = it->first;
1153 const double valw = it->second;
1154
1155 it->first = sumv>0 ? n*(1-s*sumw/sumv) : 0;
1156
1157 sumv += valv;
1158 sumw += valw;
1159 }
1160
1161 sumw = 0;
1162 sumv = 0;
1163 n = 1;
1164
1165 // Second from the right
1166 for (auto it=end-1; it!=beg-1+512; it--, n++)
1167 {
1168 const double valv = it->first;
1169 const double valw = it->second;
1170
1171 sumv += valv;
1172 sumw += valw;
1173
1174 it->first = sumv>0 ? n*(s*sumw/sumv-1) : 0;
1175 }
1176
1177 // A crosscheck has shown, that the values from the left
1178 // and right perfectly agree over the whole range. This means
1179 // the a calculation from just one side would be enough, but
1180 // doing it from both sides might still make the numerics
1181 // a bit more stable.
1182 }
1183 }
1184
1185 DrsCalibrateTime GetResult() const
1186 {
1187 DrsCalibrateTime rc(*this);
1188 rc.CalcResult();
1189 return rc;
1190 }
1191
1192 double Offset(uint32_t ch, double pos) const
1193 {
1194 const auto p = fStat.begin() + ch*1024;
1195
1196 const uint32_t f = floor(pos);
1197
1198 const double v0 = p[f].first;
1199 const double v1 = p[(f+1)%1024].first;
1200
1201 return v0 + fmod(pos, 1)*(v1-v0);
1202 }
1203
1204 double Calib(uint32_t ch, double pos) const
1205 {
1206 return pos-Offset(ch, pos);
1207 }
1208
1209 /*
1210 // See MDrsCalibrationTime
1211 std::string ReadFitsImp(const std::string &str)
1212 {
1213 fits file(str);
1214 if (!file)
1215 {
1216 std::ostringstream msg;
1217 msg << "Could not open file '" << str << "': " << strerror(errno);
1218 return msg.str();
1219 }
1220
1221 if (file.GetStr("TELESCOP")!="FACT")
1222 {
1223 std::ostringstream msg;
1224 msg << "Reading '" << str << "' failed: Not a valid FACT file (TELESCOP not FACT in header)";
1225 return msg.str();
1226 }
1227
1228 if (file.GetNumRows()!=1)
1229 {
1230 std::ostringstream msg;
1231 msg << "Reading '" << str << "' failed: Number of rows in table is not 1.";
1232 return msg.str();
1233 }
1234
1235 fNumSamples = file.GetUInt("NROI");
1236 fNumChannels = file.GetUInt("NCH");
1237 fNumEntries = file.GetUInt("NBTIME");
1238
1239 fStat.resize(fNumSamples*fNumChannels);
1240
1241 double *f = reinterpret_cast<double*>(file.SetPtrAddress("CellWidthMean"));
1242 double *s = reinterpret_cast<double*>(file.SetPtrAddress("CellWidthRms"));
1243
1244 if (!file.GetNextRow())
1245 {
1246 std::ostringstream msg;
1247 msg << "Reading data from " << str << " failed.";
1248 return msg.str();
1249 }
1250
1251 for (uint32_t i=0; i<fNumSamples*fNumChannels; i++)
1252 {
1253 fStat[i].first = f[i];
1254 fStat[i].second = s[i];
1255 }
1256
1257 return std::string();
1258 }
1259
1260 std::string WriteFitsImp(const std::string &filename, uint32_t night=0) const
1261 {
1262 ofits file(filename.c_str());
1263 if (!file)
1264 {
1265 std::ostringstream msg;
1266 msg << "Could not open file '" << filename << "': " << strerror(errno);
1267 return msg.str();
1268 }
1269
1270 file.SetDefaultKeys();
1271 file.AddColumnDouble(fStat.size(), "CellWidthMean", "ratio", "Relative cell width mean");
1272 file.AddColumnDouble(fStat.size(), "CellWidthRms", "ratio", "Relative cell width rms");
1273
1274 file.SetInt("ADCRANGE", 2000, "Dynamic range of the ADC in mV");
1275 file.SetInt("DACRANGE", 2500, "Dynamic range of the DAC in mV");
1276 file.SetInt("ADC", 12, "Resolution of ADC in bits");
1277 file.SetInt("DAC", 16, "Resolution of DAC in bits");
1278 file.SetInt("NPIX", 1440, "Number of channels in the camera");
1279 file.SetInt("NTM", 0, "Number of time marker channels");
1280 file.SetInt("NROI", fNumSamples, "Region of interest");
1281 file.SetInt("NCH", fNumChannels, "Number of chips");
1282 file.SetInt("NBTIME", fNumEntries, "Num of entries for time calibration");
1283
1284 file.WriteTableHeader("DrsCellTimes");
1285 if (night>0)
1286 file.SetInt("NIGHT", night, "Night as int");
1287
1288 std::vector<double> data(fNumSamples*fNumChannels*2);
1289
1290 for (uint32_t i=0; i<fNumSamples*fNumChannels; i++)
1291 {
1292 data[i] = fStat[i].first;
1293 data[fNumSamples*fNumChannels+i] = fStat[i].second;
1294 }
1295
1296 if (!file.WriteRow(data.data(), data.size()*sizeof(double)))
1297 {
1298 std::ostringstream msg;
1299 msg << "Writing data to " << filename << " failed.";
1300 return msg.str();
1301 }
1302
1303 return std::string();
1304 }*/
1305};
1306
1307struct DrsCalibration
1308{
1309 std::vector<int32_t> fOffset;
1310 std::vector<int64_t> fGain;
1311 std::vector<int64_t> fTrgOff;
1312
1313 int64_t fNumOffset;
1314 int64_t fNumGain;
1315 int64_t fNumTrgOff;
1316
1317 uint32_t fStep;
1318 uint16_t fRoi; // Region of interest for trgoff
1319 uint16_t fNumTm; // Number of time marker channels in trgoff
1320
1321 std::string fDateObs;
1322 std::string fDateRunBeg[3];
1323 std::string fDateRunEnd[3];
1324 std::string fDateEnd;
1325
1326// uint16_t fDAC[8];
1327
1328 DrsCalibration() :
1329 fOffset(1440*1024, 0),
1330 fGain(1440*1024, 4096),
1331 fTrgOff (1600*1024, 0),
1332 fNumOffset(1),
1333 fNumGain(2000),
1334 fNumTrgOff(1),
1335 fStep(0),
1336 fDateObs("1970-01-01T00:00:00"),
1337 fDateEnd("1970-01-01T00:00:00")
1338 {
1339 for (int i=0; i<3; i++)
1340 {
1341 fDateRunBeg[i] = "1970-01-01T00:00:00";
1342 fDateRunEnd[i] = "1970-01-01T00:00:00";
1343 }
1344 }
1345
1346 DrsCalibration(const DrsCalibration &cpy) :
1347 fOffset(cpy.fOffset),
1348 fGain(cpy.fGain),
1349 fTrgOff(cpy.fTrgOff),
1350 fNumOffset(cpy.fNumOffset),
1351 fNumGain(cpy.fNumGain),
1352 fNumTrgOff(cpy.fNumTrgOff),
1353 fStep(cpy.fStep),
1354 fRoi(cpy.fRoi),
1355 fNumTm(cpy.fNumTm),
1356 fDateObs(cpy.fDateObs),
1357 fDateEnd(cpy.fDateEnd)
1358 {
1359 for (int i=0; i<3; i++)
1360 {
1361 fDateRunBeg[i] = cpy.fDateRunBeg[i];
1362 fDateRunEnd[i] = cpy.fDateRunEnd[i];
1363 }
1364 }
1365
1366 void Clear()
1367 {
1368 // Default gain:
1369 // 0.575*[45590]*2.5V / 2^16 = 0.99999 V
1370 fOffset.assign(1440*1024, 0);
1371 fGain.assign (1440*1024, 4096);
1372 fTrgOff.assign(1600*1024, 0);
1373
1374 fNumOffset = 1;
1375 fNumGain = 2000;
1376 fNumTrgOff = 1;
1377
1378 fStep = 0;
1379
1380 fDateObs = "1970-01-01T00:00:00";
1381 fDateEnd = "1970-01-01T00:00:00";
1382
1383 for (int i=0; i<3; i++)
1384 {
1385 fDateRunBeg[i] = "1970-01-01T00:00:00";
1386 fDateRunEnd[i] = "1970-01-01T00:00:00";
1387 }
1388 }
1389
1390 std::string ReadFitsImp(const std::string &str, std::vector<float> &vec)
1391 {
1392 fits file(str);
1393 if (!file)
1394 {
1395 std::ostringstream msg;
1396 msg << "Could not open file '" << str << "': " << strerror(errno);
1397 return msg.str();
1398 }
1399
1400 if (file.GetStr("TELESCOP")!="FACT")
1401 {
1402 std::ostringstream msg;
1403 msg << "Reading '" << str << "' failed: Not a valid FACT file (TELESCOP not FACT in header)";
1404 return msg.str();
1405 }
1406
1407 if (!file.HasKey("STEP"))
1408 {
1409 std::ostringstream msg;
1410 msg << "Reading '" << str << "' failed: Is not a DRS calib file (STEP not found in header)";
1411 return msg.str();
1412 }
1413
1414 if (file.GetNumRows()!=1)
1415 {
1416 std::ostringstream msg;
1417 msg << "Reading '" << str << "' failed: Number of rows in table is not 1.";
1418 return msg.str();
1419 }
1420
1421 fStep = file.GetUInt("STEP");
1422 fNumOffset = file.GetInt("NBOFFSET");
1423 fNumGain = file.GetInt("NBGAIN");
1424 fNumTrgOff = file.GetInt("NBTRGOFF");
1425 fRoi = file.GetUInt("NROI");
1426 fNumTm = file.HasKey("NTM") ? file.GetUInt("NTM") : 0;
1427
1428 if (file.HasKey("DATE-OBS"))
1429 fDateObs = file.GetStr("DATE-OBS");
1430 if (file.HasKey("DATE-END"))
1431 fDateEnd = file.GetStr("DATE-END");
1432
1433 if (file.HasKey("RUN0-BEG"))
1434 fDateRunBeg[0]= file.GetStr("RUN0-BEG");
1435 if (file.HasKey("RUN1-BEG"))
1436 fDateRunBeg[1]= file.GetStr("RUN1-BEG");
1437 if (file.HasKey("RUN2-BEG"))
1438 fDateRunBeg[2]= file.GetStr("RUN2-BEG");
1439 if (file.HasKey("RUN0-END"))
1440 fDateRunEnd[0]= file.GetStr("RUN0-END");
1441 if (file.HasKey("RUN1-END"))
1442 fDateRunEnd[1]= file.GetStr("RUN1-END");
1443 if (file.HasKey("RUN2-END"))
1444 fDateRunEnd[2]= file.GetStr("RUN2-END");
1445/*
1446 fDAC[0] = file.GetUInt("DAC_A");
1447 fDAC[1] = file.GetUInt("DAC_B");
1448 fDAC[4] = file.GetUInt("DAC_C");
1449*/
1450 vec.resize(1440*1024*4 + (1440+fNumTm)*fRoi*2 + 4);
1451
1452 float *base = vec.data();
1453
1454 reinterpret_cast<uint32_t*>(base)[0] = fRoi;
1455
1456 file.SetPtrAddress("RunNumberBaseline", base+1, 1);
1457 file.SetPtrAddress("RunNumberGain", base+2, 1);
1458 file.SetPtrAddress("RunNumberTriggerOffset", base+3, 1);
1459 file.SetPtrAddress("BaselineMean", base+4+0*1024*1440, 1024*1440);
1460 file.SetPtrAddress("BaselineRms", base+4+1*1024*1440, 1024*1440);
1461 file.SetPtrAddress("GainMean", base+4+2*1024*1440, 1024*1440);
1462 file.SetPtrAddress("GainRms", base+4+3*1024*1440, 1024*1440);
1463 file.SetPtrAddress("TriggerOffsetMean", base+4+4*1024*1440, fRoi*1440);
1464 file.SetPtrAddress("TriggerOffsetRms", base+4+4*1024*1440+fRoi*1440, fRoi*1440);
1465 if (fNumTm>0)
1466 {
1467 file.SetPtrAddress("TriggerOffsetTMMean", base+4+4*1024*1440+ 2*fRoi*1440, fRoi*fNumTm);
1468 file.SetPtrAddress("TriggerOffsetTMRms", base+4+4*1024*1440+ 2*fRoi*1440+ fRoi*fNumTm, fRoi*fNumTm);
1469 }
1470
1471 if (!file.GetNextRow())
1472 {
1473 std::ostringstream msg;
1474 msg << "Reading data from " << str << " failed.";
1475 return msg.str();
1476 }
1477/*
1478 fDAC[2] = fDAC[1];
1479 fDAC[4] = fDAC[1];
1480
1481 fDAC[5] = fDAC[4];
1482 fDAC[6] = fDAC[4];
1483 fDAC[7] = fDAC[4];
1484*/
1485 fOffset.resize(1024*1440);
1486 fGain.resize(1024*1440);
1487
1488 fTrgOff.resize(fRoi*(1440+fNumTm));
1489
1490 // Convert back to ADC counts: 256/125 = 4096/2000
1491 // Convert back to sum (mean * num_entries)
1492 for (int i=0; i<1024*1440; i++)
1493 {
1494 fOffset[i] = fNumOffset *256*base[i+1024*1440*0+4]/125;
1495 fGain[i] = fNumOffset*fNumGain*256*base[i+1024*1440*2+4]/125;
1496 }
1497
1498 for (int i=0; i<fRoi*1440; i++)
1499 fTrgOff[i] = fNumOffset*fNumTrgOff*256*base[i+1024*1440*4+4]/125;
1500
1501 for (int i=0; i<fRoi*fNumTm; i++)
1502 fTrgOff[i+1440*fRoi] = fNumOffset*fNumTrgOff*256*base[i+1024*1440*4+2*fRoi*1440+4]/125;
1503
1504
1505 // DAC: 0..2.5V == 0..65535
1506 // V-mV: 1000
1507 //fNumGain *= 2500*50000;
1508 //for (int i=0; i<1024*1440; i++)
1509 // fGain[i] *= 65536;
1510 if (fStep==0)
1511 {
1512 for (int i=0; i<1024*1440; i++)
1513 fGain[i] = fNumOffset*4096;
1514 }
1515 else
1516 {
1517 fNumGain *= 1953125;
1518 for (int i=0; i<1024*1440; i++)
1519 fGain[i] *= 1024;
1520 }
1521
1522 // Now mark the stored DRS data as "officially valid"
1523 // However, this is not thread safe. It only ensures that
1524 // this data is not used before it is completely and correctly
1525 // read.
1526 fStep++;
1527
1528 return std::string();
1529 }
1530
1531 std::string WriteFitsImp(const std::string &filename, const std::vector<float> &vec, uint32_t night=0) const
1532 {
1533 const size_t n = 1440*1024*4 + 1440*fRoi*2 + fNumTm*fRoi*2 + 3;
1534
1535 ofits file(filename.c_str());
1536 if (!file)
1537 {
1538 std::ostringstream msg;
1539 msg << "Could not open file '" << filename << "': " << strerror(errno);
1540 return msg.str();
1541 }
1542
1543 file.AddColumnInt("RunNumberBaseline");
1544 file.AddColumnInt("RunNumberGain");
1545 file.AddColumnInt("RunNumberTriggerOffset");
1546
1547 file.AddColumnFloat(1024*1440, "BaselineMean", "mV");
1548 file.AddColumnFloat(1024*1440, "BaselineRms", "mV");
1549 file.AddColumnFloat(1024*1440, "GainMean", "mV");
1550 file.AddColumnFloat(1024*1440, "GainRms", "mV");
1551 file.AddColumnFloat(fRoi*1440, "TriggerOffsetMean", "mV");
1552 file.AddColumnFloat(fRoi*1440, "TriggerOffsetRms", "mV");
1553 file.AddColumnFloat(fRoi*fNumTm, "TriggerOffsetTMMean", "mV");
1554 file.AddColumnFloat(fRoi*fNumTm, "TriggerOffsetTMRms", "mV");
1555
1556 file.SetDefaultKeys();
1557 if (night>0)
1558 file.SetInt("NIGHT", night, "Night as int");
1559
1560 file.SetStr("DATE-OBS", fDateObs, "First event of whole DRS calibration");
1561 file.SetStr("DATE-END", fDateEnd, "Last event of whole DRS calibration");
1562 file.SetStr("RUN0-BEG", fDateRunBeg[0], "First event of run 0");
1563 file.SetStr("RUN1-BEG", fDateRunBeg[1], "First event of run 1");
1564 file.SetStr("RUN2-BEG", fDateRunBeg[2], "First event of run 2");
1565 file.SetStr("RUN0-END", fDateRunEnd[0], "Last event of run 0");
1566 file.SetStr("RUN1-END", fDateRunEnd[1], "Last event of run 1");
1567 file.SetStr("RUN2-END", fDateRunEnd[2], "Last event of run 2");
1568
1569 file.SetInt("STEP", fStep, "");
1570
1571 file.SetInt("ADCRANGE", 2000, "Dynamic range of the ADC in mV");
1572 file.SetInt("DACRANGE", 2500, "Dynamic range of the DAC in mV");
1573 file.SetInt("ADC", 12, "Resolution of ADC in bits");
1574 file.SetInt("DAC", 16, "Resolution of DAC in bits");
1575 file.SetInt("NPIX", 1440, "Number of channels in the camera");
1576 file.SetInt("NTM", fNumTm, "Number of time marker channels");
1577 file.SetInt("NROI", fRoi, "Region of interest");
1578
1579 file.SetInt("NBOFFSET", fNumOffset, "Num of entries for offset calibration");
1580 file.SetInt("NBGAIN", fNumGain/1953125, "Num of entries for gain calibration");
1581 file.SetInt("NBTRGOFF", fNumTrgOff, "Num of entries for trigger offset calibration");
1582
1583 // file.WriteKeyNT("DAC_A", fData.fDAC[0], "Level of DAC 0 in DAC counts") ||
1584 // file.WriteKeyNT("DAC_B", fData.fDAC[1], "Leval of DAC 1-3 in DAC counts") ||
1585 // file.WriteKeyNT("DAC_C", fData.fDAC[4], "Leval of DAC 4-7 in DAC counts") ||
1586
1587 file.WriteTableHeader("DrsCalibration");
1588
1589 if (!file.WriteRow(vec.data()+1, n*sizeof(float)))
1590 {
1591 std::ostringstream msg;
1592 msg << "Writing data to " << filename << " failed.";
1593 return msg.str();
1594 }
1595
1596 return std::string();
1597 }
1598
1599
1600 std::string ReadFitsImp(const std::string &str)
1601 {
1602 std::vector<float> vec;
1603 return ReadFitsImp(str, vec);
1604 }
1605
1606 bool IsValid() { return fStep>2; }
1607
1608 bool Apply(float *vec, const int16_t *val, const int16_t *start, uint32_t roi)
1609 {
1610 if (roi!=fRoi)
1611 {
1612 for (size_t ch=0; ch<1440; ch++)
1613 {
1614 const size_t pos = ch*roi;
1615 const size_t drs = ch*1024;
1616
1617 DrsCalibrate::ApplyCh(vec+pos, val+pos, start[ch], roi,
1618 fOffset.data()+drs, fNumOffset,
1619 fGain.data() +drs, fNumGain);
1620 }
1621
1622 return false;
1623 }
1624
1625 for (size_t ch=0; ch<1440; ch++)
1626 {
1627 const size_t pos = ch*fRoi;
1628 const size_t drs = ch*1024;
1629
1630 DrsCalibrate::ApplyCh(vec+pos, val+pos, start[ch], roi,
1631 fOffset.data()+drs, fNumOffset,
1632 fGain.data() +drs, fNumGain,
1633 fTrgOff.data()+pos, fNumTrgOff);
1634 }
1635
1636 for (size_t ch=0; ch<fNumTm; ch++)
1637 {
1638 const size_t pos = (ch+1440)*fRoi;
1639 const size_t drs = (ch*9+8)*1024;
1640
1641 DrsCalibrate::ApplyCh(vec+pos, val+pos, start[ch], roi,
1642 fOffset.data()+drs, fNumOffset,
1643 fGain.data() +drs, fNumGain,
1644 fTrgOff.data()+pos, fNumTrgOff);
1645 }
1646
1647 return true;
1648 }
1649};
1650
1651#endif
Note: See TracBrowser for help on using the repository browser.