Belle II Software  release-05-02-19
ECLDspEmulator.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Mikhail Remnev *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 //ECL
11 #include <ecl/utility/ECLDspEmulator.h>
12 //Framework
13 #include <framework/logging/Logger.h>
14 
15 namespace Belle2 {
20  namespace ECL {
21  template<typename T>
22  T setInRange(T val, T min, T max)
23  {
24  if (val < min) return min;
25  if (val > max) return max;
26  return val;
27  }
28 
29 
30  namespace ShapeFitter {
31  bool amplitudeOverflow(long long amp)
32  {
33  // There are 18 bits reserved for the amplitude,
34  // in range [-128, 262015]
35  static const long long max_amp = 0x3FFFF - 128;
36  return amp > max_amp;
37  }
38  }
39  }
41 }
42 
43 namespace Belle2 {
48  namespace ECL {
49  template <typename INT>
50  ECLShapeFit lftda_(const INT* f, const INT* f1, const INT* fg41,
51  const INT* fg43, const INT* fg31, const INT* fg32,
52  const INT* fg33, int* y, int ttrig2, int la_thr,
53  int hit_thr, int skip_thr, int k_a, int k_b,
54  int k_c, int k_16, int k1_chi, int k2_chi,
55  int chi_thres, bool adjusted_timing)
56  {
57  // Typical plot of y_i (i=0..31)
58  // +-------------------------------------------------------+
59  // | |
60  // | ------------------------- |
61  // | ^ xxxx |
62  // | | xx x |
63  // | | x x |
64  // | | x xx |
65  // | A1 | xx x |
66  // | (ampl) | x x |
67  // | | x xx |
68  // | | xx x |
69  // | | x x |
70  // | | xx x |
71  // | | x xx |
72  // | | xx xx |
73  // | v x xxxx |
74  // |xxxxxxxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxx
75  // | ^ |
76  // |<-|-------------------> |
77  // | | n_16 |
78  // | | |
79  // | | C1 |
80  // | v (pedestal) |
81  // +-------------------------------------------------------+
82  //
83  // 0 10 20 30
84  // (point number)
85  //
86 
87  using namespace ShapeFitter;
88 
89  enum QualityFlag {
90  c_GoodQuality = 0,
91  c_InternalError = 1,
92  c_LowAmp = 2,
93  c_BadChi2 = 3
94  };
95 
96  /*** VARIABLE DEFINITION ***/
97 
98  const static long long int k_np[16] = {
99  65536,
100  32768,
101  21845,
102  16384,
103  13107,
104  10923,
105  9362,
106  8192,
107  7282,
108  6554,
109  5958,
110  5461,
111  5041,
112  4681,
113  4369,
114  4096
115  };
116 
117  // Trigger time 0-23
118  const int ttrig = ttrig2 > 0 ? ttrig2 / 6 : 0;
119  // Number of points for pedestal measurement
120  const int n_16 = 16;
121 
122  if (k_16 + n_16 != 16) {
123  B2ERROR("Disagreement in the number of points: "
124  << k_16 << " and " << n_16);
125  }
126 
127  // Initial time index
128  int it0;
129  if (!adjusted_timing) {
130  it0 = 48 + ((23 - ttrig) << 2);
131  } else {
132  it0 = 48 + ((143 - ttrig2) << 2) / 6;
133  }
134  // Min time index, max time index
135  int it_l = 0, it_h = 191;
136  // Time index, must be within [it_l, it_h]
137  int it = setInRange(it0, it_l, it_h);
138  // Amplitude from the fit without time correction
139  // (assuming t_0 == trigger time)
140  long long A2;
141  // Amplitude from the fit with time correction
142  long long A1;
143  // Estimation of (Ampl * delta T) value
144  long long B1;
145  // Pedestal, used to calculate chi^2
146  long long C1 = 0;
147 
148  // Quality flag (https://confluence.desy.de/display/BI/ECL+Quality+flag)
149  int validity_code = c_GoodQuality;
150  bool skip_fit = false;
151  // Set to true if amplitude is less than SKIP_THR
152  bool skip_thr_flag = false;
153 
154  //== Calculate sum of first 16 points in the waveform.
155  // This sum is used for pedestal estimation.
156 
157  long long z00 = 0;
158  for (int i = k_16; i < 16; i++)
159  z00 += y[i];
160 
161  const int kz_s = 0;
162  const long long z0 = z00 >> kz_s;
163 
164  // Struct with fit results
165  ECLShapeFit result;
166  result.fit.resize(31);
167 
168  //== Check for pedestal amplitude overflow
169 
170  if (z00 > 0x3FFFF) {
171  skip_fit = true;
172  validity_code = c_InternalError;
173  A1 = -128;
174  }
175 
176  /*** FIRST APPROXIMATION ***/
177 
178  //== First approximation without time correction
179  // (assuming t_0 == trigger time)
180 
181  A2 = fg41[ttrig * 16] * z0;
182 
183  for (int i = 1; i < 16; i++)
184  A2 += y[15 + i] * (long long)fg41[ttrig * 16 + i];
185 
186  A2 += (1 << (k_a - 1));
187  A2 >>= k_a;
188 
189  //== Check if amplitude estimation is too large.
190 
191  if (amplitudeOverflow(A2)) {
192  A1 = A2 >> 3;
193  if (amplitudeOverflow(A1)) A1 >>= 3;
194  A1 -= 112;
195  B2DEBUG(15, A1 << " 2");
196 
197  validity_code = c_InternalError;
198  skip_fit = true;
199  }
200 
201  //== Check if amplitude is less than LA_THR
202  // (https://confluence.desy.de/display/BI/Electronics+Thresholds)
203 
204  bool low_ampl = false;
205  if (A2 < la_thr) low_ampl = true;
206 
207  /*** MAIN PART ***/
208 
209  //== Main part of the algorithm, estimate amplitude, time and pedestal
210 
211  // Time estimation in ADC ticks.
212  // (See TDC time in https://confluence.desy.de/display/BI/ECL+Technical+Notes)
213  int T = 0;
214 
215  if (!skip_fit && !low_ampl) {
216  for (int iter = 1; iter <= 3; iter++) {
217 
218  //== Get amplitude and (ampl*time)
219  // at time index 'it'
220 
221  A1 = fg31[it * 16] * z0;
222  B1 = fg32[it * 16] * z0;
223 
224  for (int i = 1; i < 16; i++) {
225  A1 += fg31[it * 16 + i] * (long long)y[15 + i];
226  B1 += fg32[it * 16 + i] * (long long)y[15 + i];
227  }
228  A1 += (1 << (k_a - 1));
229  A1 >>= k_a;
230 
231  //== Check if amplitude estimation is negative.
232 
233  if (A1 < -128) {
234  if (A2 > 0) {
235  A1 = A2;
236  validity_code = c_LowAmp;
237  } else {
238  A1 = -128;
239  validity_code = c_InternalError;
240  }
241 
242  skip_fit = true;
243  break;
244  }
245 
246  //== Check if amplitude estimation is too large.
247 
248  if (amplitudeOverflow(A1)) {
249  A1 >>= 3;
250  if (amplitudeOverflow(A1)) A1 >>= 3;
251  A1 -= 112;
252  B2DEBUG(15, A1 << " 1");
253 
254  validity_code = c_InternalError;
255  skip_fit = true;
256  break;
257  }
258 
259  //== Check if amplitude is less than LA_THR
260  // (https://confluence.desy.de/display/BI/Electronics+Thresholds)
261 
262  if (A1 < la_thr) {
263  it = it0;
264  low_ampl = true;
265  break;
266  }
267 
268  //== Check if amplitude is less than SKIP_THR
269  // (https://confluence.desy.de/display/BI/Electronics+Thresholds)
270 
271  if (A1 < skip_thr) {
272  skip_thr_flag = true;
273  }
274 
275  // Internal variables for fit algo
276  long long B2, B3;
277 
278  //== Estimate t_new as t_0 - B/A
279 
280  if (k_b >= 13) {
281  B2 = B1 >> (k_b - 13);
282  } else {
283  B2 = B1 << (13 - k_b);
284  }
285  B2 += (A1 << 13);
286  B3 = (B2 / A1);
287 
288  int it_uncut = (B3 + 16) >> 5;
289  int delta_it = it_uncut - 256;
290  delta_it = setInRange(delta_it, -191, 191);
291 
292  int it_current = it;
293 
294  it += delta_it;
295  it = setInRange(it, it_l, it_h);
296 
297  if (iter == 3) {
298  int delta_T;
299 
300  if (B3 >= 0x37FE)
301  delta_T = 2047;
302  else if ((B3 >> 1) <= 0x400)
303  delta_T = -2047;
304  else
305  delta_T = ((B3 * 3 + 4) >> 3) - 3072;
306 
307  int t_uncut = it_current * 12 + delta_T;
308  t_uncut = setInRange(t_uncut, -4096, 4095);
309 
310  T = -t_uncut;
311  if (!adjusted_timing) {
312  T += ((210 - ttrig2) << 3);
313  } else {
314  T += ((215 - ttrig2) << 3) - 4;
315  }
316 
317  T = setInRange(T, -2048, 2047);
318 
319  //== Estimate pedestal
320 
321  C1 = fg33[it * 16] * z0;
322  for (int i = 1; i < 16; i++)
323  C1 += fg33[it * 16 + i] * (long long)y[15 + i];
324  C1 += (1 << (k_c - 1));
325  C1 >>= k_c;
326  }
327  } // for (int iter...)
328  } // if (!skip_fit && !low_ampl)
329 
330  //== Use initial time estimation for low amplitude
331 
332  if (!skip_fit && low_ampl) {
333  A1 = A2;
334  if (A1 < -128) {
335  skip_fit = true;
336  validity_code = c_InternalError;
337  A1 = -128;
338  } else {
339  validity_code = c_LowAmp;
340  B1 = 0;
341 
342  //== Estimate pedestal
343 
344  C1 = fg43[ttrig * 16] * z0;
345  for (int i = 1; i < 16; i++)
346  C1 += fg43[ttrig * 16 + i] * (long long)y[15 + i];
347  C1 += (1 << (k_c - 1));
348  C1 >>= k_c;
349  }
350  }
351 
352  //== Estimate chi^2
353  long long chi_sq = 0;
354 
355  if (!skip_fit) {
356  //== Get fit function values in 31 points
357 
358  int B1_chi = B1 >> k_b;
359 
360  // Points 0-15 contain identical pedestal value
361  result.fit[0] = A1 * f[it * 16] + B1_chi * f1[it * 16];
362  result.fit[0] >>= k1_chi;
363  result.fit[0] += C1;
364  for (int i = 1; i <= 15; i++) {
365  result.fit[i] = result.fit[0];
366  }
367  // Points 16-30 contain the signal
368  for (int i = 1; i < 16; i++) {
369  result.fit[i + 15] = A1 * f[it * 16 + i] + B1_chi * f1[it * 16 + i];
370  result.fit[i + 15] >>= k1_chi;
371  result.fit[i + 15] += C1;
372  }
373 
374  //== Get actual threshold for chi^2, based on amplitude fit.
375 
376  long long chi_thr;
377  chi_thr = (A1 >> 1) * (A1 >> 1);
378  chi_thr >>= (k2_chi - 2); // a1n
379  chi_thr += chi_thres; // ch2_int
380 
381  //== Get chi^2
382 
383  long long chi;
384 
385  chi = n_16 * result.fit[0] - z00;
386 
387  chi_sq = chi * chi;
388  chi_sq *= k_np[n_16 - 1];
389  chi_sq >>= 16;
390  for (int i = 1; i < 16; i++) {
391  chi = result.fit[i + 15] - y[i + 15];
392  chi_sq += chi * chi;
393  }
394 
395  if (chi_sq > chi_thr && validity_code != c_LowAmp) validity_code = c_BadChi2;
396  if (C1 < 0 && validity_code == c_GoodQuality) validity_code = c_BadChi2;
397  }
398 
399  /*** ***/
400 
401  //== Compare signal peak to HIT_THR
402  // (See https://confluence.desy.de/display/BI/Electronics+Thresholds)
403 
404  int hit_val = y[20] + y[21] - (y[12] + y[13] + y[14] + y[15]) / 2;
405  if (hit_val <= hit_thr) {
406  validity_code += 4;
407  }
408 
409  //== Compare amplitude to SKIP_THR
410  // (See https://confluence.desy.de/display/BI/Electronics+Thresholds)
411 
412  if (A1 < skip_thr && validity_code != c_InternalError) {
413  skip_thr_flag = true;
414  }
415 
416  //== Set output values
417 
418  result.amp = A1;
419  result.time = T;
420  result.quality = validity_code % 4;
421  result.hit_thr = validity_code / 4;
422  result.skip_thr = skip_thr_flag;
423  result.low_amp = low_ampl;
424  result.pedestal = result.fit[0];
425 
426  result.chi2 = chi_sq;
427 
428  return result;
429  }
430 
431  template ECLShapeFit lftda_<short>(const short* f, const short* f1, const short* fg41,
432  const short* fg43, const short* fg31, const short* fg32,
433  const short* fg33, int* y, int ttrig2, int la_thr,
434  int hit_thr, int skip_thr, int k_a, int k_b,
435  int k_c, int k_16, int k1_chi, int k2_chi,
436  int chi_thres, bool adjusted_timing);
437  template ECLShapeFit lftda_<int>(const int* f, const int* f1, const int* fg41,
438  const int* fg43, const int* fg31, const int* fg32,
439  const int* fg33, int* y, int ttrig2, int la_thr,
440  int hit_thr, int skip_thr, int k_a, int k_b,
441  int k_c, int k_16, int k1_chi, int k2_chi,
442  int chi_thres, bool adjusted_timing);
443  }
445 }
446 
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19