11 #include <ecl/utility/ECLDspEmulator.h>
13 #include <framework/logging/Logger.h>
22 T setInRange(T val, T min, T max)
24 if (val < min)
return min;
25 if (val > max)
return max;
30 namespace ShapeFitter {
31 bool amplitudeOverflow(
long long amp)
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)
84 using namespace ShapeFitter;
95 const static long long int k_np[16] = {
115 const int ttrig = ttrig2 > 0 ? ttrig2 / 6 : 0;
119 if (k_16 + n_16 != 16) {
120 B2ERROR(
"Disagreement in the number of points: "
121 << k_16 <<
" and " << n_16);
126 if (!adjusted_timing) {
127 it0 = 48 + ((23 - ttrig) << 2);
129 it0 = 48 + ((143 - ttrig2) << 2) / 6;
132 int it_l = 0, it_h = 191;
134 int it = setInRange(it0, it_l, it_h);
146 int validity_code = c_GoodQuality;
147 bool skip_fit =
false;
149 bool skip_thr_flag =
false;
155 for (
int i = k_16; i < 16; i++)
159 const long long z0 = z00 >> kz_s;
163 result.fit.resize(31);
169 validity_code = c_InternalError;
178 A2 = fg41[ttrig * 16] * z0;
180 for (
int i = 1; i < 16; i++)
181 A2 += y[15 + i] * (
long long)fg41[ttrig * 16 + i];
183 A2 += (1 << (k_a - 1));
188 if (amplitudeOverflow(A2)) {
190 if (amplitudeOverflow(A1)) A1 >>= 3;
192 B2DEBUG(15, A1 <<
" 2");
194 validity_code = c_InternalError;
201 bool low_ampl =
false;
202 if (A2 < la_thr) low_ampl =
true;
212 if (!skip_fit && !low_ampl) {
213 for (
int iter = 1; iter <= 3; iter++) {
218 A1 = fg31[it * 16] * z0;
219 B1 = fg32[it * 16] * z0;
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];
225 A1 += (1 << (k_a - 1));
233 validity_code = c_LowAmp;
236 validity_code = c_InternalError;
245 if (amplitudeOverflow(A1)) {
247 if (amplitudeOverflow(A1)) A1 >>= 3;
249 B2DEBUG(15, A1 <<
" 1");
251 validity_code = c_InternalError;
269 skip_thr_flag =
true;
278 B2 = B1 >> (k_b - 13);
280 B2 = B1 << (13 - k_b);
285 int it_uncut = (B3 + 16) >> 5;
286 int delta_it = it_uncut - 256;
287 delta_it = setInRange(delta_it, -191, 191);
292 it = setInRange(it, it_l, it_h);
299 else if ((B3 >> 1) <= 0x400)
302 delta_T = ((B3 * 3 + 4) >> 3) - 3072;
304 int t_uncut = it_current * 12 + delta_T;
305 t_uncut = setInRange(t_uncut, -4096, 4095);
308 if (!adjusted_timing) {
309 T += ((210 - ttrig2) << 3);
311 T += ((215 - ttrig2) << 3) - 4;
314 T = setInRange(T, -2048, 2047);
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));
329 if (!skip_fit && low_ampl) {
333 validity_code = c_InternalError;
336 validity_code = c_LowAmp;
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));
350 long long chi_sq = 0;
355 int B1_chi = B1 >> k_b;
358 result.fit[0] = A1 * f[it * 16] + B1_chi * f1[it * 16];
359 result.fit[0] >>= k1_chi;
361 for (
int i = 1; i <= 15; i++) {
362 result.fit[i] = result.fit[0];
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;
374 chi_thr = (A1 >> 1) * (A1 >> 1);
375 chi_thr >>= (k2_chi - 2);
376 chi_thr += chi_thres;
382 chi = n_16 * result.fit[0] - z00;
385 chi_sq *= k_np[n_16 - 1];
387 for (
int i = 1; i < 16; i++) {
388 chi = result.fit[i + 15] - y[i + 15];
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;
401 int hit_val = y[20] + y[21] - (y[12] + y[13] + y[14] + y[15]) / 2;
402 if (hit_val <= hit_thr) {
409 if (A1 < skip_thr && validity_code != c_InternalError) {
410 skip_thr_flag =
true;
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];
423 result.chi2 = chi_sq;
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);