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