Belle II Software development
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
13namespace 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
41namespace 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