Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Df250EmulatorAlgorithm_v3.cc
Go to the documentation of this file.
2 
3 // corresponds to version 0x0C12 of the fADC250 firmware
4 
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 
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) {
45  jout << " Df250EmulatorAlgorithm_v3::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) {
51  jerr << " ERROR: Df250EmulatorAlgorithm_v3::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;
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 || FORCE_DEFAULT){
72  static int counter = 0;
73  NSA = NSA_DEF;
74  NSB = NSB_DEF;
75  THR = THR_DEF;
76  NPED = NPED_DEF;
77  MAXPED = MAXPED_DEF;
78  NSAT = NSAT_DEF;
79  if (counter < 10){
80  counter++;
81  if (counter == 10) jout << " WARNING Df250EmulatorAlgorithm_v3::EmulateFirmware No Df250BORConfig == Using default values == LAST WARNING" << endl;
82  else jout << " WARNING Df250EmulatorAlgorithm_v3::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_v3::EmulateFirmware NSA: " << NSA << " NSB: " << NSB << " THR: " << THR << endl;
94  }
95 
96  if (VERBOSE > 0) jout << "Df250EmulatorAlgorithm_v3::EmulateFirmware NSA: " << NSA << " NSB: " << NSB << " THR: " << THR << endl;
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) {
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_v3::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++) {
152  for (unsigned int i=0; i < MAX_SAMPLE; i++) {
153  if ((samples[i] & 0xfff) > THR) {
154  if (VERBOSE > 1) {
155  jout << "threshold crossing at " << i << endl;
156  }
157 
158  // save threshold crossing - could be overwritten
159  TC[npulses] = i+1;
160 
161  // check that we have more than NSAT samples over threshold
162  if( NSAT>1 ){
163  int samples_over_threshold = 1;
164 
165  if(i==0) {
166  // the algorithm only terminates if we dip below threshold...
167  //for(unsigned int j=i+1; ((samples[j]&0xfff)>=THR) && (j<MAX_SAMPLE+1); j++) {
168  for(unsigned int j=i+1; ((samples[j]&0xfff)>THR) && (j<MAX_SAMPLE+1); j++) {
169  // only count samples actually above threshold
170  if ((samples[j] & 0xfff) > THR)
171  samples_over_threshold++;
172 
173  if( samples_over_threshold == NSAT ) {
174  //TC[npulses] = j+1;
175  //i=j;
176  break;
177  }
178 
179  }
180  } else {
181  for(unsigned int j=i+1; ((samples[j]&0xfff)>THR) && (j<MAX_SAMPLE+1); j++) {
182  samples_over_threshold++;
183 
184  if( samples_over_threshold == NSAT )
185  break;
186  }
187  }
188 
189  // if we couldn't find NSAT samples above threshold, move on...
190  if( samples_over_threshold != NSAT )
191  continue;
192  }
193  //else {
194  //TNSAT[npulses] = TC[npulses];
195  //}
196 
197 
198  // calculate integral
199  unsigned int ibegin;
200  if(NSB > 0)
201  ibegin = i > uint32_t(NSB) ? (i - NSB) : 0; // Set to beginning of window if too early
202  else {
203  ibegin = i - NSB;
204  if(ibegin > uint32_t(NW)) // make sure we don't start looking outside the window
205  break;
206  }
207  unsigned int iend = (i + NSA) < uint32_t(NW) ? (i + NSA) : NW; // Set to last sample if too late
208  // check to see if NSA extends beyond the end of the window
209  NSA_beyond_PTW[npulses] = (i + NSA - 1) >= uint32_t(NW);
210  for (i = ibegin; i < iend; ++i) {
211  pulse_integral[npulses] += (samples[i] & 0xfff);
212  // quality monitoring
213  if(samples[i] == 0x1fff) {
214  has_overflow_samples[npulses] = true;
215  }
216  if(samples[i] == 0x1000) {
217  has_underflow_samples[npulses] = true;
218  }
219  // count number of samples within NSA that are above thresholds
220  if( (i+1>=TC[npulses]) && ((samples[i] & 0xfff) > THR) )
221  number_samples_above_threshold[npulses]++;
222  }
223  for (; i < NW && (samples[i] & 0xfff) >= THR; ++i) {}
224  if (++npulses == max_pulses)
225  break;
226  TMIN[npulses] = i;
227  }
228  }
229 
230  // That concludes the first pass over the data.
231  // Now we can head into the fine timing pass over the data.
232 
233  uint32_t VPEAK[max_pulses] = {};
234  uint32_t TPEAK[max_pulses] = {};
235  uint16_t TMID[max_pulses] = {};
236  uint16_t VMID[max_pulses] = {};
237  uint16_t TFINE[max_pulses] = {};
238  uint32_t pulse_time[max_pulses] = {};
239 
240  // The pulse pedestal is the sum of NPED (4-15) samples at the beginning of the window
241  uint32_t pedestal = 0;
242  uint32_t VMIN = 0; // VMIN is just the average of the first 4 samples, needed for timing algorithm
243  for (unsigned int i=0; i < NPED; i++) {
244  pedestal += (samples[i] & 0xfff);
245  if(i<4)
246  VMIN += (samples[i] & 0xfff);
247  // error conditions
248  // sample larger than MaxPed
249  if ((samples[i] & 0xfff) > MAXPED) {
250  bad_pedestal = true;
251  }
252  // samples with underflow/overflow
253  if( (samples[i] == 0x1fff) || (samples[i] == 0x1000) ) {
254  bad_pedestal = true;
255  }
256  }
257  VMIN /= NPED; // compute average
258 
259  // error conditions for timing algorithm
260  //bool pedestal_underflow = false;
261  for (unsigned int i=0; i < 4; i++) {
262  // We set the "Time Quality bit 0" to 1 if any of the first 4 samples is greated than MaxPed or TET...
263  if ( ((samples[i] & 0xfff) > MAXPED) || ((samples[i] & 0xfff) > THR) ) {
264  bad_timing_pedestal = true;
265  }
266  // ... or is overflow or underflow
267  if ( (samples[i] == 0x1000) || (samples[i] == 0x1fff) ) {
268  bad_timing_pedestal = true;
269  }
270  //}
271 
272  // "If any of the first 4 samples is greater than TET the TDC will NOT proceed..."
273  // Waiit for iiit...
274  if( (samples[i] & 0xfff) > THR ) {
275  no_timing_calculation = true;
276  }
277  }
278 
279 
280  for (unsigned int p=0; p < npulses; ++p) {
281 
282  // "If any of the first 4 samples is greater than TET or underflow the TDC will NOT proceed
283  // 1. pulse time is set to TC
284  // 2. pulse peak is set to zero - not anymore!
285  // 3. Time quality bits 0 and 1 are set to 1"
286  if(no_timing_calculation) {
287  TMID[p] = TC[p];
288  TFINE[p] = 0;
289  VPEAK[p] = 0;
290  //vpeak_not_found[p] = true; // this is "time quality bit 1"
291  // "Time Quality bit 0" should already be set
292  } // should just put an else here...
293 
294  // we set up a loop so that we can break out of it at appropriate times...
295  // note that currently the timing algorithm is run when the pedestal has underflow samples,
296  // but according to the documentation, it shouldn't...
297  // NOTE that we should always run the calculation since we always need to search for the
298  // peak position, just in some edge cases we don't need to run the timing algorithm
299 
300  //while ( (!no_timing_calculation || pedestal_underflow) && true) {
301  while (true) {
302  //if (VMIN == 99999) {
303  // VPEAK[p] = 0;
304  // reportTC[p] = true;
305  // pulse_time[p] = (TC[p] << 6);
306  // break;
307  // }
308 
309  // search for the peak of the pulse
310  // has to be after the threshold crossing (NO?)
311  // has to be before the last sample
312  unsigned int ipeak;
313  for (ipeak = TC[p]; (int)ipeak < NW-1; ++ipeak) {
314  //for (ipeak = TC[p]+1; ipeak < NW-1; ++ipeak) {
315  if ((samples[ipeak] & 0xfff) < (samples[ipeak-1] & 0xfff)) {
316  VPEAK[p] = (samples[ipeak-1] & 0xfff);
317  TPEAK[p] = ipeak-1;
318  break;
319  }
320  }
321 
322  // check to see if the peak is beyond the NSA
323  if(ipeak > TC[p]+NSA)
324  vpeak_beyond_NSA[p] = true;
325 
326  if (VERBOSE > 1) {
327  jout << " pulse " << p << ": VMIN: " << VMIN
328  << " TC: " << TC[p] << " VPEAK: " << VPEAK[p] << endl;
329  }
330 
331  // set error conditions in case we didn't find the peak
332  if (VPEAK[p] == 0) {
333  TMID[p] = TC[p];
334  TFINE[p] = 0;
335  VPEAK[p] = 0;
336  vpeak_beyond_NSA[p] = true;
337  vpeak_not_found[p] = true;
338  break;
339  }
340 
341  // we have found the peak position, now ignore the timing calculation if need be...
342  if(no_timing_calculation)
343  break;
344 
345  // VMID is the half amplitude
346  VMID[p] = (VMIN + VPEAK[p]) >> 1;
347 
348 
349  // look down the leading edge for the sample that satisfies V(N1) <= VMID < V(N+1)
350  // N1 is then the coarse time
351  // note that when analyzing pulses after the first, we could end up with a time for the
352  // second pulse that is before the first one! this is a little crazy, but is how
353  // the algorithm is currently implemented
354  for (unsigned int i = TPEAK[p]; i >= 1; --i) {
355  if ( ((samples[i-1] & 0xfff) <= VMID[p]) && ((samples[i] & 0xfff) > VMID[p]) ) { // V(N1) <= VMID < V(N+1)
356  // if ( ((samples[i-1] & 0xfff) <= VMID[p]) // V(N1) <= VMID < V(N+1)
357  // || ( (samples[i-1] & 0xfff) > (samples[i] & 0xfff) ) ) { // we aren't on the leading edge anymore
358  TMID[p] = i;
359  break;
360  }
361  }
362 
363  if (TMID[p] == 0) { // handle the case where we couldn't find a coarse time - redundant?
364  TFINE[p] = 0;
365  }
366  else {
367  // fine timing algorithm (see documentation)
368  int Vnext = (samples[TMID[p]] & 0xfff);
369  int Vlast = (samples[TMID[p]-1] & 0xfff);
370  if (VERBOSE > 2) {
371  jout << " TMIN = " << TMIN[p] << " TMID = " << TMID[p] << " TPEAK = " << TPEAK[p] << endl
372  << " VMID = " << VMID[p] << " Vnext = " << Vnext << " Vlast = " << Vlast << endl;
373  }
374  if (Vnext > Vlast && VMID[p] >= Vlast)
375  TFINE[p] = 64 * (VMID[p] - Vlast) / (Vnext - Vlast);
376  else
377  TFINE[p] = 62;
378  if(TFINE[p] == 64)
379  TFINE[p] = 0;
380  }
381  pulse_time[p] = ((TMID[p]-1) << 6) + TFINE[p];
382  break;
383  }
384  VMIN = (VMIN < 99999)? VMIN : 0; // deprecated?
385 
386  if (VERBOSE > 1) {
387  jout << " pulse " << p << ": VMID: " << VMID[p] << " TMID: " << TMID[p]
388  << " TFINE: " << TFINE[p] << " time: " << pulse_time[p]
389  << " integral: " << pulse_integral[p] << endl;
390  if (VERBOSE > 2) {
391  jout << " TMIN = " << TMIN[p] << " TMID = " << TMID[p] << " TPEAK = " << TPEAK[p] << endl;
392  //<< " VMID = " << VMID[p] << " Vnext = " << Vnext << " Vlast = " << Vlast << endl;
393  }
394  }
395 
396  // algorithm is finished, fill the information
397  Df250PulseData* f250PulseData;
398  if( p < pdat_objs.size() ) {
399  f250PulseData = pdat_objs[p];
400 
401  if(f250PulseData == NULL) {
402  jerr << " NULL f250PulseData object!" << endl;
403  continue;
404  }
405  } else {
406  // make a fresh object if one does not exist
407  f250PulseData = new Df250PulseData;
408 
409  f250PulseData->rocid = rawData->rocid;
410  f250PulseData->slot = rawData->slot;
411  f250PulseData->channel = rawData->channel;
412  f250PulseData->itrigger = rawData->itrigger;
413  // word 1
414  f250PulseData->event_within_block = 1;
415  f250PulseData->QF_pedestal = bad_pedestal;
416  f250PulseData->pedestal = pedestal;
417  // word 2
418  f250PulseData->integral = pulse_integral[p];
419  f250PulseData->QF_NSA_beyond_PTW = NSA_beyond_PTW[p];
420  f250PulseData->QF_overflow = has_overflow_samples[p];
421  f250PulseData->QF_underflow = has_underflow_samples[p];
422  f250PulseData->nsamples_over_threshold = number_samples_above_threshold[p];
423  // word 3
424  f250PulseData->course_time = TMID[p];
425  f250PulseData->fine_time = TFINE[p];
426  f250PulseData->QF_vpeak_beyond_NSA = vpeak_beyond_NSA[p];
427  f250PulseData->QF_vpeak_not_found = vpeak_not_found[p];
428  f250PulseData->QF_bad_pedestal = bad_timing_pedestal;
429  // other information
430  f250PulseData->pulse_number = p;
431  f250PulseData->nsamples_integral = NSA + NSB;
432  f250PulseData->nsamples_pedestal = NPED;
433  f250PulseData->emulated = true;
434 
435  f250PulseData->AddAssociatedObject(rawData);
436  const_cast<Df250WindowRawData*>(rawData)->AddAssociatedObject(f250PulseData);
437  pdat_objs.push_back(f250PulseData);
438  }
439 
440  // copy over emulated values
441  f250PulseData->integral_emulated = pulse_integral[p];
442  f250PulseData->pedestal_emulated = pedestal;
443  f250PulseData->pulse_peak_emulated = VPEAK[p];
444  f250PulseData->course_time_emulated = TMID[p];
445  f250PulseData->fine_time_emulated = TFINE[p];
446 
447  // check the emulated quality factors as well
448  uint32_t QF = 0; // make single quality factor number for compactness
449  if( bad_pedestal ) QF |= (1<<0);
450  if( NSA_beyond_PTW[p] ) QF |= (1<<1);
451  if( has_overflow_samples[p] ) QF |= (1<<2);
452  if( has_underflow_samples[p] ) QF |= (1<<3);
453  if( vpeak_beyond_NSA[p] ) QF |= (1<<4);
454  if( vpeak_not_found[p] ) QF |= (1<<5);
455  if( bad_timing_pedestal ) QF |= (1<<6);
456  f250PulseData->QF_emulated = QF;
457 
458  if(VERBOSE > 3) {
459  cout << boolalpha;
460  cout << "bad_pedestal = " << bad_pedestal << endl;
461  cout << "NSA_beyond_PTW = " << NSA_beyond_PTW[p] << endl;
462  cout << "has_overflow_samples = " << has_overflow_samples[p] << endl;
463  cout << "has_underflow_samples = " << has_underflow_samples[p] << endl;
464  cout << "vpeak_beyond_NSA = " << vpeak_beyond_NSA[p] << endl;
465  cout << "vpeak_not_found = " << vpeak_not_found[p] << endl;
466  cout << "bad_timing_pedestal = " << bad_timing_pedestal << endl;
467  cout << "total QF = " << QF << endl;
468  }
469 
470  // if we are using the emulated values, copy them
471  if( f250PulseData->emulated ) {
472  f250PulseData->integral = f250PulseData->integral_emulated;
473  f250PulseData->pedestal = f250PulseData->pedestal_emulated;
474  f250PulseData->pulse_peak = f250PulseData->pulse_peak_emulated;
475  f250PulseData->course_time = f250PulseData->course_time_emulated;
476  f250PulseData->fine_time = f250PulseData->fine_time_emulated;
477 
478  /*
479  if( (rawData->rocid >= 31) || (rawData->rocid <= 46) ) {
480  f250PulseData->nsamples_integral = NSA + NSB;
481  }
482  */
483 
484  }
485 
486  }
487 
488  if (VERBOSE > 0) jout << " Df250EmulatorAlgorithm_v3::EmulateFirmware ==> Emulation complete <==" << endl;
489  return;
490 }
uint32_t fine_time
uint32_t pedestal_emulated
Value calculated from raw data (if available)
vector< uint16_t > samples
double counter
Definition: FitGains.C:151
uint32_t course_time_emulated
Value calculated from raw data (if available) - debug.
bool emulated
true if made from Window Raw Data
uint32_t integral
uint32_t nsamples_integral
number of samples used in integral
uint32_t pedestal
uint32_t fine_time_emulated
Value calculated from raw data (if available) - debug.
uint32_t pulse_peak_emulated
Value calculated from raw data (if available)
uint32_t integral_emulated
Value calculated from raw data (if available)
uint32_t nsamples_over_threshold
void EmulateFirmware(const Df250WindowRawData *rawData, std::vector< Df250PulseData * > &pdatt_objs)
uint32_t channel
Definition: DDAQAddress.h:34
uint32_t pulse_peak
uint32_t nsamples_pedestal
number of samples used in pedestal
const double TMIN
static TH1I * pedestal[nChan]
uint32_t rocid
Definition: DDAQAddress.h:32
uint16_t adc_thres[16]
Definition: bor_roc.h:91
uint32_t course_time
uint32_t QF_emulated
uint32_t itrigger
Definition: DDAQAddress.h:35
uint32_t event_within_block
uint32_t pulse_number
pulse number for this channel, this event starting from 0
uint32_t slot
Definition: DDAQAddress.h:33