Bug Summary

File:libraries/DAQ/Df250EmulatorAlgorithm_v2.cc
Location:line 208, column 10
Description:Division by zero

Annotated Source Code

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