Bug Summary

File:/volatile/halld/gluex/nightly/2024-03-18/Linux_CentOS7.7-x86_64-gcc4.8.5/halld_recon/src/libraries/DAQ/Df250EmulatorAlgorithm_v2.cc
Location:line 255, column 10
Description:Division by zero

Annotated Source Code

1#include <DAQ/Df250EmulatorAlgorithm_v2.h>
2
3// corresponds to version 0x0C0C or 0x0C0D of the fADC250 firmwarer
4
5Df250EmulatorAlgorithm_v2::Df250EmulatorAlgorithm_v2(JEventLoop *loop){
6 // Enables forced use of default values
7 FORCE_DEFAULT = 0;
8 // Default values for the essential parameters
9 NSA_DEF = 20;
10 NSB_DEF = 5;
11 THR_DEF = 120;
12 NPED_DEF = 4;
13 MAXPED_DEF = 512;
14 NSAT_DEF = 2;
15
16 // DEBUG
17 NSA_DEF = 15;
18 NSB_DEF = 1;
19 THR_DEF = 108;
20 NPED_DEF = 4;
21 MAXPED_DEF = 512;
22 NSAT_DEF = 2;
23
24
25 // Set verbosity
26 VERBOSE = 0;
27
28 if(gPARMS){
29 gPARMS->SetDefaultParameter("EMULATION250:FORCE_DEFAULT", FORCE_DEFAULT,"Set to >0 to force use of default values");
30 gPARMS->SetDefaultParameter("EMULATION250:NSA", NSA_DEF,"Set NSA for firmware emulation, will be overwritten by BORConfig if present");
31 gPARMS->SetDefaultParameter("EMULATION250:NSB", NSB_DEF,"Set NSB for firmware emulation, will be overwritten by BORConfig if present");
32 gPARMS->SetDefaultParameter("EMULATION250:THR", THR_DEF,"Set threshold for firmware emulation, will be overwritten by BORConfig if present");
33 gPARMS->SetDefaultParameter("EMULATION250:NPED", NPED_DEF,"Set NPED for firmware emulation, will be overwritten by BORConfig if present");
34 gPARMS->SetDefaultParameter("EMULATION250:MAXPED", MAXPED_DEF,"Set MAXPED for firmware emulation, will be overwritten by BORConfig if present");
35 gPARMS->SetDefaultParameter("EMULATION250:NSAT", NSAT_DEF,"Set NSAT for firmware emulation, will be overwritten by BORConfig if present");
36 gPARMS->SetDefaultParameter("EMULATION250:VERBOSE", VERBOSE,"Set verbosity for f250 emulation");
37 }
38}
39
40void Df250EmulatorAlgorithm_v2::EmulateFirmware(const Df250WindowRawData* rawData,
41 std::vector<Df250PulseData*> &pdat_objs)
42{
43 // This is the main routine called by JEventSource_EVIO::GetObjects() and serves as the entry point for the code.
44 if (VERBOSE > 0) {
1
Taking false branch
→
45 jout << " Df250EmulatorAlgorithm_v2::EmulateFirmware ==> Starting emulation <==" << endl;
46 jout << "rocid : " << rawData->rocid << " slot: " << rawData->slot << " channel: " << rawData->channel << endl;
47 }
48
49 // First check that we have window raw data available
50 if (rawData == NULL__null) {
2
←
Taking false branch
→
51 jerr << " ERROR: Df250EmulatorAlgorithm_v2::EmulateFirmware - raw sample data is missing" << endl;
52 jerr << " Contact mstaib@jlab.org" << endl;
53 return;
54 }
55
56 // We need the channel number to get the threshold
57 uint32_t channel = rawData->channel;
58
59 // First grab the config objects from the raw data and get the quantities we need from them
60 // The only things we need for this version of the f250 firmware are NSB, NSA, and the threshold.
61 // These are all stored in the BOR config. We can grab this from the raw data since that was already associated in JEventSource_EVIO::GetObjects.
62 const Df250BORConfig *f250BORConfig = NULL__null;
63 rawData->GetSingle(f250BORConfig);
64
65 uint32_t NSA;
66 int32_t NSB;
67 uint32_t NPED, MAXPED;
68 uint16_t THR;
69 uint16_t NSAT;
70 //If this does not exist, or we force it, use the default values
71 if (f250BORConfig == NULL__null || FORCE_DEFAULT){
72 static int counter = 0;
73 NSA = NSA_DEF;
74 NSB = NSB_DEF;
75 THR = THR_DEF;
76 NPED = NPED_DEF;
3
←
Value assigned to 'NPED'
→
77 MAXPED = MAXPED_DEF;
78 NSAT = NSAT_DEF;
79 if (counter < 10){
4
←
Taking true branch
→
80 counter++;
81 if (counter == 10) jout << " WARNING Df250EmulatorAlgorithm_v2::EmulateFirmware No Df250BORConfig == Using default values == LAST WARNING" << endl;
5
←
Taking false branch
→
82 else jout << " WARNING Df250EmulatorAlgorithm_v2::EmulateFirmware No Df250BORConfig == Using default values " << endl;
83 //<< rawData->rocid << "/" << rawData->slot << "/" << rawData->channel << endl;
84 }
85 }
86 else{
87 NSA = f250BORConfig->NSA;
88 NSB = f250BORConfig->NSB;
89 THR = f250BORConfig->adc_thres[channel];
90 NPED = f250BORConfig->NPED;
91 MAXPED = f250BORConfig->MaxPed;
92 NSAT = f250BORConfig->NSAT;
93 //if (VERBOSE > 0) jout << "Df250EmulatorAlgorithm_v2::EmulateFirmware NSA: " << NSA << " NSB: " << NSB << " THR: " << THR << endl;
94 }
95
96 if (VERBOSE > 0) jout << "Df250EmulatorAlgorithm_v2::EmulateFirmware NSA: " << NSA << " NSB: " << NSB << " THR: " << THR << endl;
6
←
Taking false branch
→
97
98 // Note that in principle we could get this information from the Df250Config objects as well, but generally only NPED and the value of NSA+NSB are saved
99 // not the individual NSA and NSB values
100
101 /*
102 // TEST
103 if( (rawData->rocid >= 31) || (rawData->rocid <= 46) ) {
104 NSA = NSA_DEF;
105 NSB = NSB_DEF;
106 }
107 */
108
109 // quality bits
110 bool bad_pedestal = false;
111 bool bad_timing_pedestal = false;
112 bool no_timing_calculation = false;
113
114 // Now we can start to loop over the raw data
115 // This requires a few passes due to some features in the way the quantities are calculated...
116 // The first step is to scan the samples for TC (threshold crossing sample) and compute the
117 // integrals of all pulses found.
118
119 vector<uint16_t> samples = rawData->samples;
120 uint16_t NW = samples.size();
121 uint32_t npulses = 0;
122 const int max_pulses = 3;
123 uint32_t TC[max_pulses] = {};
124 uint32_t TMIN[max_pulses] = {3};
125 //uint32_t TNSAT[max_pulses] = {};
126 uint32_t pulse_integral[max_pulses] = {};
127 bool has_overflow_samples[max_pulses] = {false};
128 bool has_underflow_samples[max_pulses] = {false};
129 uint32_t number_samples_above_threshold[max_pulses] = {0};
130 bool NSA_beyond_PTW[max_pulses] = {false};
131 bool vpeak_beyond_NSA[max_pulses] = {false};
132 bool vpeak_not_found[max_pulses] = {false};
133
134 // some verbose debugging output
135 if(VERBOSE > 0) {
7
←
Taking false branch
→
136 for (unsigned int i=0; i < NW; i++) {
137 if(VERBOSE > 2) {
138 if(samples[i] == 0x1fff)
139 jout << "Overflow at sample " << i << endl;
140 if(samples[i] == 0x1000)
141 jout << "Underflow at sample " << i << endl;
142 }
143 if (VERBOSE > 5) jout << "Df250EmulatorAlgorithm_v2::EmulateFirmware samples[" << i << "]: " << samples[i] << endl;
144 }
145 }
146
147
148 // look for the threhold crossings and compute the integrals
149 //unsigned int MAX_SAMPLE = (NSB>0) ? (NW-NSAT) : (NW-NSAT+NSB-1)); // check this
150 unsigned int MAX_SAMPLE = NW-NSAT;
151 for (unsigned int i=0; i < MAX_SAMPLE; i++) {
8
←
Assuming 'i' is >= 'MAX_SAMPLE'
→
9
←
Loop condition is false. Execution continues on line 231
→
152 if ((samples[i] & 0xfff) > THR) {
153 if (VERBOSE > 1) {
154 jout << "threshold crossing at " << i << endl;
155 }
156
157 // save threshold crossing - could be overwritten
158 TC[npulses] = i+1;
159
160 // check that we have more than NSAT samples over threshold
161 if( NSAT>1 ){
162 int samples_over_threshold = 1;
163
164 if(i==0) {
165 // the algorithm only terminates if we dip below threshold...
166 for(unsigned int j=i+1; ((samples[j]&0xfff)>=THR) && (j<MAX_SAMPLE+1); j++) {
167 // only count samples actually above threshold
168 if ((samples[j] & 0xfff) > THR)
169 samples_over_threshold++;
170
171 if( samples_over_threshold == NSAT ) {
172 //TC[npulses] = j+1;
173 //i=j;
174 break;
175 }
176
177 }
178 } else {
179 for(unsigned int j=i+1; ((samples[j]&0xfff)>THR) && (j<MAX_SAMPLE+1); j++) {
180 samples_over_threshold++;
181
182 if( samples_over_threshold == NSAT )
183 break;
184 }
185 }
186
187 // if we couldn't find NSAT samples above threshold, move on...
188 if( samples_over_threshold != NSAT )
189 continue;
190 }
191 //else {
192 //TNSAT[npulses] = TC[npulses];
193 //}
194
195
196 // calculate integral
197 unsigned int ibegin;
198 if(NSB > 0)
199 ibegin = i > uint32_t(NSB) ? (i - NSB) : 0; // Set to beginning of window if too early
200 else {
201 ibegin = i - NSB;
202 if(ibegin > uint32_t(NW)) // make sure we don't start looking outside the window
203 break;
204 }
205 unsigned int iend = (i + NSA) < uint32_t(NW) ? (i + NSA) : NW; // Set to last sample if too late
206 // check to see if NSA extends beyond the end of the window
207 NSA_beyond_PTW[npulses] = (i + NSA - 1) >= uint32_t(NW);
208 for (i = ibegin; i < iend; ++i) {
209 pulse_integral[npulses] += (samples[i] & 0xfff);
210 // quality monitoring
211 if(samples[i] == 0x1fff) {
212 has_overflow_samples[npulses] = true;
213 }
214 if(samples[i] == 0x1000) {
215 has_underflow_samples[npulses] = true;
216 }
217 // count number of samples within NSA that are above thresholds
218 if( (i+1>=TC[npulses]) && ((samples[i] & 0xfff) > THR) )
219 number_samples_above_threshold[npulses]++;
220 }
221 for (; i < NW && (samples[i] & 0xfff) >= THR; ++i) {}
222 if (++npulses == max_pulses)
223 break;
224 TMIN[npulses] = i;
225 }
226 }
227
228 // That concludes the first pass over the data.
229 // Now we can head into the fine timing pass over the data.
230
231 uint32_t VPEAK[max_pulses] = {};
232 uint32_t TPEAK[max_pulses] = {};
233 uint16_t TMID[max_pulses] = {};
234 uint16_t VMID[max_pulses] = {};
235 uint16_t TFINE[max_pulses] = {};
236 uint32_t pulse_time[max_pulses] = {};
237
238 // The pulse pedestal is the sum of NPED (4-15) samples at the beginning of the window
239 uint32_t pedestal = 0;
240 uint32_t VMIN = 0; // VMIN is just the average of the first 4 samples, needed for timing algorithm
241 for (unsigned int i=0; i < NPED; i++) {
10
←
Assuming 'i' is >= 'NPED'
→
11
←
Loop condition is false. Execution continues on line 255
→
242 pedestal += (samples[i] & 0xfff);
243 if(i<4)
244 VMIN += (samples[i] & 0xfff);
245 // error conditions
246 // sample larger than MaxPed
247 if ((samples[i] & 0xfff) > MAXPED) {
248 bad_pedestal = true;
249 }
250 // samples with underflow/overflow
251 if( (samples[i] == 0x1fff) || (samples[i] == 0x1000) ) {
252 bad_pedestal = true;
253 }
254 }
255 VMIN /= NPED; // compute average
12
←
Division by zero
256
257 // error conditions for timing algorithm
258 //bool pedestal_underflow = false;
259 //for (unsigned int i=0; i < 5; i++) {
260 for (unsigned int i=0; i < 4; i++) {
261 // We set the "Time Quality bit 0" to 1 if any of the first 5 samples is greated than MaxPed or TET...
262 if ( ((samples[i] & 0xfff) > MAXPED) || ((samples[i] & 0xfff) > THR) ) {
263 bad_timing_pedestal = true;
264 }
265 // ... or is overflow or underflow
266 if ( (samples[i] == 0x1000) || (samples[i] == 0x1fff) ) {
267 bad_timing_pedestal = true;
268 }
269 //}
270
271 // TEST
272 //for (unsigned int i=0; i < 4; i++) {
273
274 // "If any of the first 5 samples is greater than TET the TDC will NOT proceed..."
275 // Waiit for iiit...
276 if( (samples[i] & 0xfff) > THR ) {
277 no_timing_calculation = true;
278 }
279 }
280
281
282 for (unsigned int p=0; p < npulses; ++p) {
283 // TEST
284 //TC[p] = TNSAT[p];
285
286 // "If any of the first 5 samples is greater than TET or underflow the TDC will NOT proceed
287 // 1. pulse time is set to TC
288 // 2. pulse peak is set to zero
289 // 3. Time quality bits 0 and 1 are set to 1"
290 if(no_timing_calculation) {
291 TMID[p] = TC[p];
292 TFINE[p] = 0;
293 VPEAK[p] = 0;
294 vpeak_not_found[p] = true; // this is "time quality bit 1"
295 // "Time Quality bit 0" should already be set
296 } // should just put an else here...
297
298 // we set up a loop so that we can break out of it at appropriate times...
299 // note that currently the timing algorithm is run when the pedestal has underflow samples,
300 // but according to the documentation, it shouldn't...
301 //while ( (!no_timing_calculation || pedestal_underflow) && true) {
302 while ( (!no_timing_calculation) && true) {
303 //if (VMIN == 99999) {
304 // VPEAK[p] = 0;
305 // reportTC[p] = true;
306 // pulse_time[p] = (TC[p] << 6);
307 // break;
308 // }
309
310 // search for the peak of the pulse
311 // has to be after the threshold crossing (NO?)
312 // has to be before the last sample
313 unsigned int ipeak;
314 for (ipeak = TC[p]; (int)ipeak < NW-1; ++ipeak) {
315 //for (ipeak = TC[p]+1; ipeak < NW-1; ++ipeak) {
316 if ((samples[ipeak] & 0xfff) < (samples[ipeak-1] & 0xfff)) {
317 VPEAK[p] = (samples[ipeak-1] & 0xfff);
318 TPEAK[p] = ipeak-1;
319 break;
320 }
321 }
322
323 // check to see if the peak is beyond the NSA
324 if(ipeak > TC[p]+NSA)
325 vpeak_beyond_NSA[p] = true;
326
327 if (VERBOSE > 1) {
328 jout << " pulse " << p << ": VMIN: " << VMIN
329 << " TC: " << TC[p] << " VPEAK: " << VPEAK[p] << endl;
330 }
331
332 // set error conditions in case we didn't find the peak
333 if (VPEAK[p] == 0) {
334 TMID[p] = TC[p];
335 TFINE[p] = 0;
336 VPEAK[p] = 0;
337 vpeak_beyond_NSA[p] = true;
338 vpeak_not_found[p] = true;
339 break;
340 }
341
342 // VMID is the half amplitude
343 VMID[p] = (VMIN + VPEAK[p]) >> 1;
344
345 /*
346 for (unsigned int i = TMIN[p] + 1; i < (uint32_t)ipeak; ++i) { // old
347 if ((samples[i] & 0xfff) > VMID[p]) {
348 TMID[p] = i;
349 break;
350 }
351 }
352 */
353
354 // look down the leading edge for the sample that satisfies V(N1) <= VMID < V(N+1)
355 // N1 is then the coarse time
356 // note that when analyzing pulses after the first, we could end up with a time for the
357 // second pulse that is before the first one! this is a little crazy, but is how
358 // the algorithm is currently implemented
359 for (unsigned int i = TPEAK[p]; i >= 1; --i) {
360 if ( ((samples[i-1] & 0xfff) <= VMID[p]) && ((samples[i] & 0xfff) > VMID[p]) ) { // V(N1) <= VMID < V(N+1)
361 // if ( ((samples[i-1] & 0xfff) <= VMID[p]) // V(N1) <= VMID < V(N+1)
362 // || ( (samples[i-1] & 0xfff) > (samples[i] & 0xfff) ) ) { // we aren't on the leading edge anymore
363 TMID[p] = i;
364 break;
365 }
366 }
367
368 if (TMID[p] == 0) { // handle the case where we couldn't find a coarse time - redundant?
369 TFINE[p] = 0;
370 }
371 else {
372 // fine timing algorithm (see documentation)
373 int Vnext = (samples[TMID[p]] & 0xfff);
374 int Vlast = (samples[TMID[p]-1] & 0xfff);
375 if (VERBOSE > 2) {
376 jout << " TMIN = " << TMIN[p] << " TMID = " << TMID[p] << " TPEAK = " << TPEAK[p] << endl
377 << " VMID = " << VMID[p] << " Vnext = " << Vnext << " Vlast = " << Vlast << endl;
378 }
379 if (Vnext > Vlast && VMID[p] >= Vlast)
380 TFINE[p] = 64 * (VMID[p] - Vlast) / (Vnext - Vlast);
381 else
382 TFINE[p] = 62;
383 if(TFINE[p] == 64)
384 TFINE[p] = 0;
385 }
386 pulse_time[p] = ((TMID[p]-1) << 6) + TFINE[p];
387 break;
388 }
389 VMIN = (VMIN < 99999)? VMIN : 0; // deprecated?
390
391 if (VERBOSE > 1) {
392 jout << " pulse " << p << ": VMID: " << VMID[p] << " TMID: " << TMID[p]
393 << " TFINE: " << TFINE[p] << " time: " << pulse_time[p]
394 << " integral: " << pulse_integral[p] << endl;
395 if (VERBOSE > 2) {
396 jout << " TMIN = " << TMIN[p] << " TMID = " << TMID[p] << " TPEAK = " << TPEAK[p] << endl;
397 //<< " VMID = " << VMID[p] << " Vnext = " << Vnext << " Vlast = " << Vlast << endl;
398 }
399 }
400
401 // algorithm is finished, fill the information
402 Df250PulseData* f250PulseData;
403 if( p < pdat_objs.size() ) {
404 f250PulseData = pdat_objs[p];
405
406 if(f250PulseData == NULL__null) {
407 jerr << " NULL f250PulseData object!" << endl;
408 continue;
409 }
410 } else {
411 // make a fresh object if one does not exist
412 f250PulseData = new Df250PulseData;
413
414 f250PulseData->rocid = rawData->rocid;
415 f250PulseData->slot = rawData->slot;
416 f250PulseData->channel = rawData->channel;
417 f250PulseData->itrigger = rawData->itrigger;
418 // word 1
419 f250PulseData->event_within_block = 1;
420 f250PulseData->QF_pedestal = bad_pedestal;
421 f250PulseData->pedestal = pedestal;
422 // word 2
423 f250PulseData->integral = pulse_integral[p];
424 f250PulseData->QF_NSA_beyond_PTW = NSA_beyond_PTW[p];
425 f250PulseData->QF_overflow = has_overflow_samples[p];
426 f250PulseData->QF_underflow = has_underflow_samples[p];
427 f250PulseData->nsamples_over_threshold = number_samples_above_threshold[p];
428 // word 3
429 f250PulseData->course_time = TMID[p];
430 f250PulseData->fine_time = TFINE[p];
431 f250PulseData->QF_vpeak_beyond_NSA = vpeak_beyond_NSA[p];
432 f250PulseData->QF_vpeak_not_found = vpeak_not_found[p];
433 f250PulseData->QF_bad_pedestal = bad_timing_pedestal;
434 // other information
435 f250PulseData->pulse_number = p;
436 f250PulseData->nsamples_integral = NSA + NSB;
437 f250PulseData->nsamples_pedestal = NPED;
438 f250PulseData->emulated = true;
439
440 f250PulseData->AddAssociatedObject(rawData);
441 const_cast<Df250WindowRawData*>(rawData)->AddAssociatedObject(f250PulseData);
442 pdat_objs.push_back(f250PulseData);
443 }
444
445 // copy over emulated values
446 f250PulseData->integral_emulated = pulse_integral[p];
447 f250PulseData->pedestal_emulated = pedestal;
448 f250PulseData->pulse_peak_emulated = VPEAK[p];
449 f250PulseData->course_time_emulated = TMID[p];
450 f250PulseData->fine_time_emulated = TFINE[p];
451
452 // check the emulated quality factors as well
453 uint32_t QF = 0; // make single quality factor number for compactness
454 if( bad_pedestal ) QF |= (1<<0);
455 if( NSA_beyond_PTW[p] ) QF |= (1<<1);
456 if( has_overflow_samples[p] ) QF |= (1<<2);
457 if( has_underflow_samples[p] ) QF |= (1<<3);
458 if( vpeak_beyond_NSA[p] ) QF |= (1<<4);
459 if( vpeak_not_found[p] ) QF |= (1<<5);
460 if( bad_timing_pedestal ) QF |= (1<<6);
461 f250PulseData->QF_emulated = QF;
462
463 if(VERBOSE > 3) {
464 cout << boolalpha;
465 cout << "bad_pedestal = " << bad_pedestal << endl;
466 cout << "NSA_beyond_PTW = " << NSA_beyond_PTW[p] << endl;
467 cout << "has_overflow_samples = " << has_overflow_samples[p] << endl;
468 cout << "has_underflow_samples = " << has_underflow_samples[p] << endl;
469 cout << "vpeak_beyond_NSA = " << vpeak_beyond_NSA[p] << endl;
470 cout << "vpeak_not_found = " << vpeak_not_found[p] << endl;
471 cout << "bad_timing_pedestal = " << bad_timing_pedestal << endl;
472 cout << "total QF = " << QF << endl;
473 }
474
475 // if we are using the emulated values, copy them
476 if( f250PulseData->emulated ) {
477 f250PulseData->integral = f250PulseData->integral_emulated;
478 f250PulseData->pedestal = f250PulseData->pedestal_emulated;
479 f250PulseData->pulse_peak = f250PulseData->pulse_peak_emulated;
480 f250PulseData->course_time = f250PulseData->course_time_emulated;
481 f250PulseData->fine_time = f250PulseData->fine_time_emulated;
482
483 /*
484 if( (rawData->rocid >= 31) || (rawData->rocid <= 46) ) {
485 f250PulseData->nsamples_integral = NSA + NSB;
486 }
487 */
488
489 }
490
491 }
492
493 if (VERBOSE > 0) jout << " Df250EmulatorAlgorithm_v2::EmulateFirmware ==> Emulation complete <==" << endl;
494 return;
495}