Bug Summary

File:libraries/DAQ/fa125algo.cc
Location:line 195, column 22
Description:Division by zero

Annotated Source Code

1//
2// $Id$
3//
4//
5// (See comments in fa125algo16.h)
6//
7// The #defines below might normally be placed in the header file, but
8// these have some generic names and are intended only for this CDC
9// timing algorithm. They may also be replaced with JANA configuration
10// parameters and/or CCDB constants.
11//
12// The main entry point here is the fa125_algos routine. The samples
13// are passed in along with references to many different return values.
14//
15
16#include <stdlib.h>
17
18#include <fa125algo.h>
19
20#ifndef _DBG_cout<<"libraries/DAQ/fa125algo.cc"<<":"<<20
<<" "
21#define _DBG_cout<<"libraries/DAQ/fa125algo.cc"<<":"<<21
<<" " cout<<__FILE__"libraries/DAQ/fa125algo.cc"<<":"<<__LINE__21<<" "
22#define _DBG__cout<<"libraries/DAQ/fa125algo.cc"<<":"<<22
<<endl cout<<__FILE__"libraries/DAQ/fa125algo.cc"<<":"<<__LINE__22<<endl
23#endif
24
25// FA125 emulation functions cdc_hit, cdc_time etc for NPK=1 (only the first peak is identified)
26
27// The ADC data buffer has NW samples in total, comprising
28// [unused samples][NPED samples for pedestal][samples WINDOW_START to WINDOW_END][20 or more extra samples]
29//
30// WINDOW_START is the sample immediately after the NPED samples in the pedestal window.
31// The hit search starts with sample WINDOW_START+PG and ends including sample WINDOW_END. This is for ease of use with older data.
32// In the new firmware, there are no unused samples at the start, and 20 samples at the end.
33// The first sample is numbered sample 0.
34// The pedestal returned is the mean of the NPED2 samples ending PED samples before the sample containing the hit threshold crossing
35//
36// This requires that:
37// WINDOW_START >= NPED
38// WINDOW_END+20 < NW (number of samples supplied)
39// NPED2 >= NPED
40//
41//
42// Peak amplitude and integral are returned without pedestal subtraction
43
44 // cdc_time q_code values:
45
46 // q_code Time returned Condition
47 // 0 Leading edge time Good
48 // 1 X*10 - 29 Sample value of 0 found
49 // 1 X*10 - 28 Sample value greater than PED_MAX found in adc[0 to PED]
50 // 1 X*10 - 27 Sample values lie below the high timing threshold
51 // 1 TCL*10 + 4 Low timing threshold crossing sample TCL occurs too late to upsample
52 // 1 TCL*10 + 5 One or more upsampled values are negative
53 // 1 TCL*10 + 9 The upsampled values are too low
54
55
56
57
58// Wrapper for routine below
59//void fa125_algos(int rocid, vector<uint16_t> samples, fa125_algos_data_t &fa125_algos_data)
60
61void fa125_algos(int rocid, vector<uint16_t> samples, fa125_algos_data_t &fa125_algos_data, uint32_t CDC_WS, uint32_t CDC_WE, uint32_t CDC_IE, uint32_t CDC_NP, uint32_t CDC_NP2, uint32_t CDC_PG, uint32_t CDC_H, uint32_t CDC_TH, uint32_t CDC_TL, uint32_t FDC_WS, uint32_t FDC_WE, uint32_t FDC_IE, uint32_t FDC_NP, uint32_t FDC_NP2, uint32_t FDC_PG, uint32_t FDC_H, uint32_t FDC_TH, uint32_t FDC_TL)
62{
63 // Since we would have to write a lot of "cdc_algos_data."'s, make another
64 // reference that's smaller so we just have to write "d." instead
65 fa125_algos_data_t &d = fa125_algos_data;
66
67 bool cdchit = kFALSE;
68 if ((rocid > 24)&&(rocid < 29)) cdchit = kTRUE; //CDC has rocid 25 to 28
1
Assuming 'rocid' is <= 24
→
69
70 if (cdchit) {
2
←
Taking false branch
→
71
72 d.WINDOW_START = CDC_WS;
73 d.WINDOW_END = CDC_WE;
74 d.INT_END = CDC_IE;
75
76 d.NPED = CDC_NP;
77 d.NPED2 = CDC_NP2;
78 d.PG = CDC_PG;
79
80 d.HIT_THRES = CDC_H;
81 d.HIGH_THRESHOLD = CDC_TH;
82 d.LOW_THRESHOLD = CDC_TL;
83
84 } else {
85
86 d.WINDOW_START = FDC_WS;
87 d.WINDOW_END = FDC_WE;
88 d.INT_END = FDC_IE;
89
90 d.NPED = FDC_NP;
91 d.NPED2 = FDC_NP2;
92 d.PG = FDC_PG;
93
94 d.HIT_THRES = FDC_H;
95 d.HIGH_THRESHOLD = FDC_TH;
96 d.LOW_THRESHOLD = FDC_TL;
97
98 }
99
100
101 if (samples.size()<=(uint32_t)d.WINDOW_END + 20) {
3
←
Taking false branch
→
102 cout << "The number of samples passed into the fa125_algos routine (" << samples.size() << ") is less than the" << endl;
103 cout << "minimum required by the parameters in use (" << d.WINDOW_END+21 << "). " << endl;
104 cout << "Parameter WE (" << d.WINDOW_END << ") should be decreased to " << samples.size()-21 << " or less." << endl;
105 exit(-1);
106 }
107
108 if (d.NPED2 > d.NPED) {
4
←
Taking false branch
→
109 cout << "Parameter NPED is too small or NPED2 is too large. " << endl;
110 cout << "NPED (" << d.NPED << ") should be increased or NPED2 (" << d.NPED2 << ") decreased until NPED >= NPED2." << endl;
111 exit(-1);
112 }
113
114 if (d.WINDOW_START < d.NPED) {
5
←
Taking false branch
→
115 cout << "Parameter WS (" << d.WINDOW_START << ") is too small or NPED (" << d.NPED << ") too large." << endl;
116 cout << "WS should be >= NPED." << endl;
117 exit(-1);
118 }
119
120 // Copy uint16_t samples into Int_t type array so we can pass it into the cdc_algos2
121 // routine that does the actual work
122 Int_t adc[d.WINDOW_END+21];
123 for(uint32_t i=0; i<=(uint32_t)d.WINDOW_END+20; i++) adc[i] = (Int_t)samples[i];
6
←
Loop condition is true. Entering loop body
→
7
←
Loop condition is false. Execution continues on line 126
→
124
125 // Call the actual routine that does the heavy lifting
126 fa125_algos(d.time, d.q_code, d.pedestal, d.integral, d.overflows, d.maxamp, adc, d.WINDOW_START, d.WINDOW_END, d.INT_END, d.NPED, d.NPED2, d.PG, d.HIT_THRES, d.HIGH_THRESHOLD, d.LOW_THRESHOLD);
8
←
Calling 'fa125_algos'
→
127
128}
129
130
131void fa125_algos(Int_t &time, Int_t &q_code, Int_t &pedestal, Long_t &integral, Int_t &overflows, Int_t &maxamp, Int_t adc[], Int_t WINDOW_START, Int_t WINDOW_END, Int_t INT_END, Int_t NPED, Int_t NPED2, Int_t PG, Int_t HIT_THRES, Int_t HIGH_THRESHOLD, Int_t LOW_THRESHOLD) {
132
133
134 const Int_t NU = 20; //number of samples sent to time algo
135 const Int_t PED = 5; //sample to be used as pedestal for timing is in place 5
136
137 const Int_t XTHR_SAMPLE = PED + PG;
138
139
140 Int_t adc_subset[NU];
141
142 Int_t hitfound=0; //hit found or not (1=found,0=not)
143 Int_t hitsample=-1; // if hit found, sample number of threshold crossing
144 Int_t timesample=0;
145
146 Int_t i=0;
147
148 time=0; // hit time in 0.1xsamples since start of buffer passed to cdc_time
149 q_code=-1; // quality code, 0=good, 1=returned rough estimate
150 pedestal=0; // pedestal just before hit
151 integral=0; // signal integral, total
152 overflows=0; // count of samples with overflow bit set (need raw data, not possible from my root files)
153 maxamp=0; // signal amplitude at first max after hit
154
155
156 // look for hit using mean pedestal of NPED samples before trigger
157 cdc_hit(hitfound, hitsample, pedestal, adc, WINDOW_START, WINDOW_END, HIT_THRES, NPED, NPED2, PG);
9
←
Passing value via 8th parameter 'NPED'
→
10
←
Calling 'cdc_hit'
→
158
159
160 if (hitfound==1) {
161
162 for (i=0; i<NU; i++) {
163 adc_subset[i] = adc[hitsample+i-XTHR_SAMPLE];
164 }
165
166 cdc_time(time, q_code, adc_subset, NU, PG, HIGH_THRESHOLD, LOW_THRESHOLD);
167
168 timesample = hitsample-XTHR_SAMPLE + (Int_t)(0.1*time); //sample number containing leading edge sample
169
170 cdc_integral(integral, overflows, timesample, adc, WINDOW_END, INT_END);
171
172 cdc_max(maxamp, hitsample, adc, WINDOW_END);
173
174 time = 10*(hitsample-XTHR_SAMPLE) + time; // integer number * 0.1 samples
175
176 }
177
178}
179
180
181
182
183void cdc_hit(Int_t &hitfound, Int_t &hitsample, Int_t &pedestal, Int_t adc[], Int_t WINDOW_START, Int_t WINDOW_END, Int_t HIT_THRES, Int_t NPED, Int_t NPED2, Int_t PG) {
184
185
186 pedestal=0; //pedestal
187 Int_t threshold=0;
188
189 Int_t i=0;
190
191 // calc pedestal as mean of NPED samples before trigger
192
193 for (i=0; i<NPED; i++) pedestal += adc[WINDOW_START-1-NPED+i];
11
←
Assuming 'i' is >= 'NPED'
→
12
←
Loop condition is false. Execution continues on line 195
→
194
195 pedestal = pedestal/NPED; // Integer div is ok as fpga will do 2 rightshifts
13
←
Division by zero
196
197 threshold = pedestal + HIT_THRES;
198
199 // look for threshold crossing
200 i = WINDOW_START - 1 + PG;
201 hitfound = 0;
202
203 while ((hitfound==0) && (i<WINDOW_END-1)) {
204
205 i++;
206
207 if (adc[i] >= threshold) {
208 if (adc[i+1] >= threshold) {
209 hitfound = 1;
210 hitsample = i;
211 }
212 }
213 }
214
215 if (hitfound == 1) {
216
217 //calculate new pedestal ending just before the hit
218
219 pedestal = 0;
220
221 for (i=0; i<NPED2; i++) {
222 pedestal += adc[hitsample-PG-i];
223 }
224
225 pedestal = pedestal/NPED2;
226 }
227
228
229}
230
231
232
233
234void cdc_integral(Long_t& integral, Int_t& overflows, Int_t timesample, Int_t adc[], Int_t WINDOW_END, Int_t INT_END) {
235
236 Int_t i=0;
237
238 integral = 0;
239 overflows = 0;
240
241 Int_t lastsample = timesample + INT_END - 1;
242
243 if (lastsample > WINDOW_END) lastsample = WINDOW_END;
244
245 for (i = timesample; i <= lastsample; i++ ) {
246
247 integral += (Long_t)adc[i];
248 if (adc[i]==(Long_t)4095) overflows++; // only a placeholder at present. need to test on sample's overflow bit
249
250 }
251
252
253}
254
255
256
257
258void cdc_max(Int_t& maxamp, Int_t hitsample, Int_t adc[], Int_t WINDOW_END) {
259
260 maxamp = 0;
261
262
263
264 int nd = 0; // num decreasing
265
266 int i = hitsample; //start value for max amp sample
267 int imax = i;
268 int pfound = 1;
269 int debug=0;
270 int peakfound = 0;
271
272 if (debug>4) printf(" Starting with potential peak at xthr adc[%i] %i\n",i,adc[i]);
273
274 while (i<WINDOW_END && !peakfound) {
275
276 i++;
277
278 if (debug>4) printf(" Inspecting adc[%i] %i\n",i,adc[i]);
279
280 if (adc[i] > adc[i-1]) {
281
282 pfound = 1;
283 imax = i;
284 nd = 0; //reset count of descending values
285 if (debug>4) printf(" Increase - this could be a peak\n");
286
287 } else if (adc[i] < adc[i-1]) {
288
289 if (pfound) {
290 nd++;
291 if (nd>1) peakfound = 1;
292 if (debug>4) printf(" Decrease after potential peak - increase decrease counter to %i\n",nd);
293 }
294
295 }
296 }
297
298 maxamp = adc[imax];
299
300
301}
302
303
304
305
306
307
308void cdc_time(Int_t &le_time, Int_t &q_code, Int_t adc[], Int_t NU, Int_t PG, Int_t THRES_HIGH, Int_t THRES_LOW) {
309
310
311 // adc[NU] array of samples
312 // NU=20 size of array
313 // PG pedestal gap - hit threshold crossing is PG samples after adc[PED] - the single sample to be used as pedestal here
314 // THRES_HIGH high timing threshold (eg 4 x pedestal-width )
315 // THRES_LOW high timing threshold (eg 1 x pedestal-width )
316 //
317 // le_time leading edge time as 0.1x number of samples since first sample supplied
318 // q_code quality code, 0=good, >0=rough estimate (firmware returns 0=good, 1=not so good)
319 //
320 // q_code Time returned Condition
321 // 0 Leading edge time Good
322 // 1 X*10 - 29 Sample value of 0 found
323 // 1 X*10 - 28 Sample value greater than PED_MAX found in adc[0 to PED]
324 // 1 X*10 - 27 Sample values lie below the high timing threshold
325 // 1 TCL*10 + 4 Low timing threshold crossing sample TCL occurs too late to upsample
326 // 1 TCL*10 + 5 One or more upsampled values are negative
327 // 1 TCL*10 + 9 The upsampled values are too low
328 //
329
330
331
332 const Int_t NUPSAMPLED = 6; // number of upsampled values to calculate, minimum is 6
333 const Int_t SET_ADC_MIN = 20; // adjust adc values so that the min is at 20
334 const Int_t LIMIT_PED_MAX = 511; // max acceptable value in adc[0 to PED]
335 const Int_t PED = 5; // take local pedestal to be adc[PED]
336
337 const Int_t START_SEARCH = PED+1; // -- start looking for hi threshold xing with this sample
338
339 const Int_t X = PED + PG; // hit threshold crossing sample is adc[X]
340 const Int_t ROUGH_TIME = (X*10)-30; // -- add onto this to return rough time estimates
341
342 Int_t iubuf[NUPSAMPLED] = {0}; // array of upsampled values; iubuf[0] maps to low thres xing sample
343
344 Int_t adc_thres_hi = 0; // high threshold
345 Int_t adc_thres_lo = 0; // low threshold
346
347 // -- contributions to hit time, these are summed together eventually, units of sample/10
348 Int_t itime1 = 0; // which sample
349 Int_t itime2 = 0; // which minisample
350 Int_t itime3 = 0; // correction from interpolation
351
352 // -- search vars
353 Int_t adc_sample_hi = 0; // integer range 0 to NU := 0; --sample number for adc val at or above hi thres
354 Int_t adc_sample_lo = 0; // integer range 0 to NU := 0; -- sample num for adc val at or below lo thres
355 Int_t adc_sample_lo2 = 0; // integer range 0 to 12:= 0; -- minisample num for adc val at or below lo thres
356
357 Bool_t over_threshold = kFALSE;
358 Bool_t below_threshold = kFALSE;
359
360
361 // upsampling checks
362 Int_t ups_adjust = 0;
363
364 Int_t i = 0;
365
366
367
368 //check all samples are >0
369
370 Bool_t adczero = kFALSE;
371
372 i = 0;
373 while ((!adczero)&&(i<NU)) {
374
375 if (adc[i] == 0) {
376 adczero = kTRUE;
377 }
378
379 i++;
380 }
381
382
383 if (adczero) {
384
385 le_time = ROUGH_TIME + 1;
386 q_code = 1;
387 return;
388
389 }
390
391
392
393
394 //check all samples from 0 to pedestal are <= LIMIT_PED_MAX
395
396 Bool_t pedlimit = kFALSE;
397
398 i = 0;
399 while ((!pedlimit)&&(i<PED+1)) {
400
401 if (adc[i] > LIMIT_PED_MAX) {
402 pedlimit = kTRUE;
403 }
404
405 i++;
406 }
407
408
409
410 if (pedlimit) {
411
412 le_time = ROUGH_TIME + 2;
413 q_code = 2;
414 return;
415
416 }
417
418 // add offset to move min val in subset equal to SET_ADC_MIN
419 // this is to move samples away from 0 to avoid upsampled pts going -ve (on a curve betw 2 samples)
420
421 Int_t adcmin = 4095;
422
423 i=0;
424
425 while (i<NU) {
426
427 if (adc[i] < adcmin) {
428 adcmin = adc[i];
429 }
430
431 i++;
432 }
433
434 Int_t adcoffset = SET_ADC_MIN - adcmin;
435
436 i=0;
437
438 while (i<NU) {
439 adc[i] = adc[i] + adcoffset;
440 i++;
441 }
442
443 // eg if adcmin is 100, setmin is 30, adcoffset = 30 - 100 = -70, move adc down by 70
444
445
446 //////////////////////////////
447
448 // calc thresholds
449
450 adc_thres_hi = adc[PED] + THRES_HIGH;
451 adc_thres_lo = adc[PED] + THRES_LOW;
452
453 // search for high threshold crossing
454
455 over_threshold = kFALSE;
456 i = START_SEARCH;
457
458 while ((!over_threshold)&&(i<NU)) {
459
460 if (adc[i] >= adc_thres_hi) {
461 adc_sample_hi = i;
462 over_threshold = kTRUE;
463 }
464
465 i++;
466 }
467
468
469 if (!over_threshold) {
470
471 le_time = ROUGH_TIME + 3;
472 q_code = 3;
473 return;
474
475 }
476
477
478 // search for low threshold crossing
479
480 below_threshold = kFALSE;
481 i = adc_sample_hi-1;
482
483 while ((!below_threshold) && (i>=PED)) {
484
485 if (adc[i] <= adc_thres_lo) {
486 adc_sample_lo = i;
487 itime1 = i*10;
488 below_threshold = kTRUE;
489 }
490
491 i--;
492 }
493
494
495
496 if (adc[adc_sample_lo] == adc_thres_lo) { // no need to upsample
497
498 le_time = itime1;
499 q_code = 0;
500 return;
501
502 }
503
504
505 if (adc_sample_lo > NU-7) { // too late to upsample
506
507 le_time = itime1 + 4;
508 q_code = 4;
509 return;
510
511 }
512
513
514
515 //upsample values from adc_sample_lo to adc_sample_lo + 1
516
517 upsamplei(adc, adc_sample_lo, iubuf, NUPSAMPLED);
518
519
520
521 //check upsampled values are >0
522
523 Bool_t negups = kFALSE;
524
525 i=0;
526 while ((!negups)&&(i<NUPSAMPLED)) {
527
528 if (iubuf[i] < 0 ) {
529 negups = kTRUE;
530 }
531
532 i++;
533 }
534
535
536 if (negups) {
537
538 le_time = itime1 + 5;
539 q_code = 5;
540 return;
541
542 }
543
544
545 // correct errors
546 // iubuf[0] should be equal to adc[adc_sample_lo] and iubuf[5] should equal adc[adc_sample_lo+1]
547 // very steep pulse gradients cause errors in upsampling with larger errors in the later values
548 // match iubuf[0] to adc[adc_sample_lo] so that the threshold crossing must be at or after iubuf[0]
549
550 ups_adjust = iubuf[0] - adc[adc_sample_lo];
551
552 // move threshold correspondingly instead of correcting upsampled values
553
554 adc_thres_lo = adc_thres_lo + ups_adjust;
555
556 // check that threshold crossing lies within the range of iubuf[0 to 5]
557
558 if (iubuf[NUPSAMPLED-1]<= adc_thres_lo) { //bad upsampling
559
560 le_time = itime1 + 9; //midway
561 q_code = 6;
562 return;
563
564 }
565
566
567
568 // search through upsampled array
569
570 below_threshold = kFALSE;
571 i = NUPSAMPLED-2;
572
573 while ((!below_threshold) && (i>=0)) {
574
575 if (iubuf[i] <= adc_thres_lo) {
576 adc_sample_lo2 = i;
577 below_threshold = kTRUE;
578 }
579
580 i--;
581 }
582
583
584
585 if (!below_threshold) { //upsampled points did not go below thres
586
587 printf("upsampled points did not go below threshold - should be impossible\n");
588 le_time = 0;
589 q_code = 9;
590 return;
591 }
592
593
594
595
596 itime2 = adc_sample_lo2*2; // convert from sample/5 to sample/10
597
598
599
600 //interpolate
601
602 itime3 = 0;
603
604
605 if (iubuf[adc_sample_lo2] != adc_thres_lo) {
606
607 if (2*adc_thres_lo >= iubuf[adc_sample_lo2] + iubuf[adc_sample_lo2+1]) itime3 = 1;
608
609 }
610
611
612 le_time = itime1 + itime2 + itime3; // -- this is time from first sample point, in 1/10ths of samples
613 q_code = 0;
614
615
616}
617
618
619void upsamplei(Int_t x[], Int_t startpos, Int_t z[], const Int_t NUPSAMPLED) {
620
621 // x is array of samples
622 // z is array of upsampled data
623 // startpos is where to start upsampling in array x, only need to upsample a small region
624
625
626 const Int_t nz = NUPSAMPLED;
627
628 // const Int_t Kscale = 32768;
629 // const Int_t K[43]={-8,-18,-27,-21,10,75,165,249,279,205,-2,-323,-673,-911,-873,-425,482,1773,3247,4618,5591,5943,5591,4618,3247,1773,482,-425,-873,-911,-673,-323,-2,205,279,249,165,75,10,-21,-27,-18,-8}; //32768
630
631 Int_t k,j,dk;
632
633
634 const Int_t Kscale = 16384;
635 const Int_t K[43] = {-4, -9, -13, -10, 5, 37, 82, 124, 139, 102, -1, -161, -336, -455, -436, -212, 241, 886, 1623, 2309, 2795, 2971, 2795, 2309, 1623, 886, 241, -212, -436, -455, -336, -161, -1, 102, 139, 124, 82, 37, 5, -10, -13, -9, -4};
636
637
638 //don't need to calculate whole range of k possible
639 //earliest value k=42 corresponds to sample 4.2
640 // k=43 sample 4.4
641 // k=46 sample 5.0
642
643
644 // sample 4 (if possible) would be at k=41
645 // sample 4.2 k=42
646 // sample 5 k=46
647
648 // sample x k=41 + (x-4)*5
649 // sample x-0.2 k=40 + (x-4)*5
650
651
652 Int_t firstk = 41 + (startpos-4)*5;
653
654
655 for (k=firstk; k<firstk+nz; k++) {
656
657 dk = k - firstk;
658
659 z[dk]=0.0;
660
661 for (j=k%5;j<43;j+=5) {
662
663 z[dk] += x[(k-j)/5]*K[j];
664
665 }
666
667 // printf("dk %i z %i 5z %i 5z/scale %i\n",dk,z[dk],5.0*z[dk],5.0*z[dk]/Kscale);
668
669 z[dk] = (Int_t)(5*z[dk])/Kscale;
670
671 }
672
673}
674