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