CbmRoot
PronyFitter.cxx
Go to the documentation of this file.
1 #include "PronyFitter.h"
2 #include "PolynomComplexRoots.h"
3 #include "PolynomRealRoots.h"
4 
5 namespace PsdSignalFitting {
6 
7 
8  PronyFitter::PronyFitter(int model_order,
9  int exponent_number,
10  int gate_beg,
11  int gate_end) {
12  Initialize(model_order, exponent_number, gate_beg, gate_end);
13  }
14 
15  void PronyFitter::Initialize(int model_order,
16  int exponent_number,
17  int gate_beg,
18  int gate_end) {
19  fModelOrder = model_order;
20  fExpNumber = exponent_number;
21  fGateBeg = gate_beg;
22  fGateEnd = gate_end;
23  AllocData();
24  }
25 
27  fz = new std::complex<float>[fExpNumber + 1];
28  fh = new std::complex<float>[fExpNumber + 1];
29  for (int i = 0; i < fExpNumber + 1; i++) {
30  fz[i] = {0., 0.};
31  fh[i] = {0., 0.};
32  }
33  }
34 
35  void PronyFitter::SetWaveform(std::vector<uint16_t>& uWfm,
36  uint16_t uZeroLevel) {
37  fuWfm = uWfm;
38  fuZeroLevel = uZeroLevel;
39  fSampleTotal = fuWfm.size();
40  fuFitWfm.resize(fSampleTotal);
41  }
42 
43  int PronyFitter::CalcSignalBegin(float front_time_beg_03,
44  float front_time_end) {
45  return std::ceil((3 * front_time_beg_03 - front_time_end) / 2.);
46  }
47 
48  void PronyFitter::SetSignalBegin(int SignalBeg) {
49  fSignalBegin = SignalBeg;
50  if (fIsDebug)
51  printf("\nsignal begin %i zero level %u\n", fSignalBegin, fuZeroLevel);
52  }
53 
55  float rho_f = 999;
56  float rho_b = 999;
57  std::vector<float> a_f;
58  std::vector<float> a_b;
59 
60  CovarianceDirect(rho_f, a_f, rho_b, a_b);
61  //CovarianceQRmod(rho_f, a_f, rho_b, a_b);
62 
63  if (fIsDebug) {
64  printf("LSerr %e, SLE roots forward ", rho_f);
65  for (uint8_t i = 0; i < a_f.size(); i++)
66  printf(" %e ", a_f.at(i));
67  printf("\n");
68  printf("LSerr %e, SLE roots backward ", rho_b);
69  for (uint8_t i = 0; i < a_b.size(); i++)
70  printf(" %e ", a_b.at(i));
71  printf("\n");
72  }
73 
74  float* a_arr = new float[fModelOrder + 1];
75  for (int i = 0; i < fModelOrder + 1; i++)
76  a_arr[i] = 0.;
77 
78  //for(uint8_t i = 0; i < a_f.size(); i++) a_arr[i+1] = 0.5*(a_f.at(i) + a_b.at(i));
79  for (uint8_t i = 0; i < a_f.size(); i++)
80  a_arr[i + 1] = a_f.at(i);
81  a_arr[0] = 1.;
82 
83  float* zr = new float[fModelOrder];
84  float* zi = new float[fModelOrder];
85  for (int i = 0; i < fModelOrder; i++) {
86  zr[i] = 0.;
87  zi[i] = 0.;
88  }
89 
90  int total_roots;
91  //polynomRealRoots(zr, fModelOrder, a_arr, total_roots);
92  polynomComplexRoots(zr, zi, fModelOrder, a_arr, total_roots);
93  if (fIsDebug) {
94  printf("forward polinom roots ");
95  for (int i = 0; i < fModelOrder; i++)
96  printf("%.5f%+.5fi ", zr[i], zi[i]);
97  printf("\n");
98 
99  printf("forward freqs ");
100  for (int i = 0; i < fModelOrder; i++)
101  printf("%.5f ", log(zr[i]));
102  printf("\n");
103  }
104 
105  std::complex<float>* z_arr = new std::complex<float>[fExpNumber + 1];
106  for (int i = 0; i < fExpNumber; i++) {
107  if (std::isfinite(zr[i]))
108  z_arr[i + 1] = {zr[i], zi[i]};
109  else
110  z_arr[i + 1] = 0.;
111  }
112  z_arr[0] = {1., 0.};
113  SetHarmonics(z_arr);
114  fTotalPolRoots = total_roots;
115 
116  delete[] a_arr;
117  delete[] zr;
118  delete[] zi;
119  delete[] z_arr;
120  }
121 
122  void PronyFitter::CovarianceDirect(float& /*rho_f*/,
123  std::vector<float>& a_f,
124  float& /*rho_b*/,
125  std::vector<float>& /*a_b*/) {
126  std::vector<int32_t> kiWfmSignal;
127  //filtering constant fraction
128  for (int sample_curr = fSignalBegin; sample_curr <= fGateEnd; sample_curr++)
129  kiWfmSignal.push_back(fuWfm.at(sample_curr) - fuZeroLevel);
130  int n = kiWfmSignal.size();
131 
132  float** Rkm_arr = new float*[fModelOrder];
133  for (int i = 0; i < fModelOrder; i++) {
134  Rkm_arr[i] = new float[fModelOrder];
135  for (int j = 0; j < fModelOrder; j++)
136  Rkm_arr[i][j] = 0.;
137  }
138 
139  float* R0k_arr = new float[fModelOrder];
140  for (int i = 0; i < fModelOrder; i++)
141  R0k_arr[i] = 0.;
142 
143  //Regression via linear prediction forward
144  for (int i = 1; i <= fModelOrder; i++)
145  for (int j = 1; j <= fModelOrder; j++)
146  for (int sample_curr = fModelOrder; sample_curr < n; sample_curr++)
147  Rkm_arr[i - 1][j - 1] += (float) (kiWfmSignal.at(sample_curr - i)
148  * kiWfmSignal.at(sample_curr - j));
149 
150  for (int i = 1; i <= fModelOrder; i++)
151  for (int sample_curr = fModelOrder; sample_curr < n; sample_curr++)
152  R0k_arr[i - 1] -= (float) (kiWfmSignal.at(sample_curr)
153  * kiWfmSignal.at(sample_curr - i));
154 
155  if (fIsDebug) {
156  printf("system forward\n");
157  for (int i = 0; i < fModelOrder; i++) {
158  for (int j = 0; j < fModelOrder; j++)
159  printf("%e ", Rkm_arr[i][j]);
160  printf("%e\n", R0k_arr[i]);
161  }
162  }
163 
164  float* a = new float[fModelOrder];
165  for (int i = 0; i < fModelOrder; i++)
166  a[i] = 0.;
167 
168  SolveSLEGauss(a, Rkm_arr, R0k_arr, fModelOrder);
169  //SolveSLECholesky(a, Rkm_arr, R0k_arr, fModelOrder);
170  if (fIsDebug) {
171  printf("SLE roots ");
172  for (int i = 0; i < fModelOrder; i++)
173  printf(" %e ", a[i]);
174  printf("\n");
175  }
176 
177  a_f.resize(fModelOrder);
178  for (int i = 0; i < fModelOrder; i++)
179  a_f.at(i) = a[i];
180 
181  for (int i = 0; i < fModelOrder; i++)
182  delete[] Rkm_arr[i];
183  delete[] Rkm_arr;
184  delete[] R0k_arr;
185  delete[] a;
186  }
187 
188  void PronyFitter::CovarianceQRmod(float& rho_f,
189  std::vector<float>& a_f,
190  float& rho_b,
191  std::vector<float>& a_b) {
192 
193  /*
194 % Copyright (c) 2019 by S. Lawrence Marple Jr.
195 %
196 %
197 p
198 --
199 order of linear prediction/autoregressive filter
200 %
201 x
202 --
203 vector of data samples
204 %
205 rho_f
206 --
207 least squares estimate of forward linear prediction variance
208 %
209 a_f
210 --
211 vector of forward linear prediction/autoregressive
212 parameters
213 %
214 rho_b
215 --
216 least squares estimate of backward linear prediction
217 variance
218 %
219 a_b
220 --
221 vector of backward linear prediction/autoregressive
222 parameters
223 
224 */
225 
226 
227  //******** Initialization ********
228 
229  std::vector<int32_t> kiWfmSignal;
230  //filtering constant fraction
231  for (int sample_curr = fSignalBegin; sample_curr <= fGateEnd; sample_curr++)
232  kiWfmSignal.push_back(fuWfm.at(sample_curr) - fuZeroLevel);
233  int n = kiWfmSignal.size();
234  if (2 * fModelOrder + 1 > n) {
235  if (fIsDebug)
236  printf("ERROR: Order too high; will make solution singular\n");
237  return;
238  }
239 
240  float r1 = 0.;
241  for (int k = 1; k <= n - 2; k++)
242  r1 += std::pow(kiWfmSignal.at(k), 2);
243 
244  float r2 = std::pow(kiWfmSignal.front(), 2);
245  float r3 = std::pow(kiWfmSignal.back(), 2);
246 
247  rho_f = r1 + r3;
248  rho_b = r1 + r2;
249  r1 = 1. / (rho_b + r3);
250 
251  std::vector<float> c, d;
252  c.push_back(kiWfmSignal.back() * r1);
253  d.push_back(kiWfmSignal.front() * r1);
254 
255  float gam = 0.;
256  float del = 0.;
257  std::vector<float> ef, eb, ec, ed;
258  std::vector<float> coeffs;
259  coeffs.resize(6);
260 
261  for (int v_iter = 0; v_iter < n; v_iter++) {
262  ef.push_back(kiWfmSignal.at(v_iter));
263  eb.push_back(kiWfmSignal.at(v_iter));
264  ec.push_back(c.at(0) * kiWfmSignal.at(v_iter));
265  ed.push_back(d.at(0) * kiWfmSignal.at(v_iter));
266  }
267 
268  //******** Main Recursion ********
269 
270  for (int k = 1; k <= fModelOrder; k++) {
271  if (rho_f <= 0 || rho_b <= 0) {
272  if (fIsDebug)
273  printf("PsdPronyFitter::ERROR: prediction squared error was less "
274  "than or equal to zero\n");
275  return;
276  }
277 
278  gam = 1. - ec.at(n - k);
279  del = 1. - ed.front();
280  if (gam <= 0 || gam > 1 || del <= 0 || del > 1) {
281  if (fIsDebug)
282  printf("PsdPronyFitter::ERROR: GAM or DEL gain factor not in "
283  "expected range 0 to 1\n");
284  return;
285  }
286 
287  // computation for k-th order reflection coefficients
288  std::vector<float> eff, ebb;
289  eff.resize(n);
290  ebb.resize(n);
291  float delta = 0.;
292  for (int v_iter = 0; v_iter < n - 1; v_iter++) {
293  eff.at(v_iter) = ef.at(v_iter + 1);
294  ebb.at(v_iter) = eb.at(v_iter);
295  delta += eff.at(v_iter) * ebb.at(v_iter);
296  }
297 
298  float k_f = -delta / rho_b;
299  float k_b = -delta / rho_f;
300 
301  //% order updates for squared prediction errors rho_f and rho_b
302  rho_f = rho_f * (1 - k_f * k_b);
303  rho_b = rho_b * (1 - k_f * k_b);
304 
305  //% order updates for linear prediction parameter arrays a_f and a_b
306  std::vector<float> temp_af = a_f;
307  std::reverse(std::begin(a_b), std::end(a_b));
308  for (uint8_t i = 0; i < a_f.size(); i++)
309  a_f.at(i) += k_f * a_b.at(i);
310  a_f.push_back(k_f);
311 
312  std::reverse(std::begin(a_b), std::end(a_b));
313  std::reverse(std::begin(temp_af), std::end(temp_af));
314  for (uint8_t i = 0; i < a_b.size(); i++)
315  a_b.at(i) += k_b * temp_af.at(i);
316  a_b.push_back(k_b);
317 
318  //% check if maximum order has been reached
319  if (k == fModelOrder) {
320  rho_f = rho_f / (n - fModelOrder);
321  rho_b = rho_b / (n - fModelOrder);
322  return;
323  }
324 
325  //% order updates for prediction error arrays ef and eb
326  for (int v_iter = 0; v_iter < n; v_iter++) {
327  eb.at(v_iter) = ebb.at(v_iter) + k_b * eff.at(v_iter);
328  ef.at(v_iter) = eff.at(v_iter) + k_f * ebb.at(v_iter);
329  }
330 
331  //% coefficients for next set of updates
332  coeffs.at(0) = ec.front();
333  coeffs.at(1) = coeffs.at(0) / del;
334  coeffs.at(2) = coeffs.at(0) / gam;
335 
336  //% time updates for gain arrays c' and d"
337  std::vector<float> temp_c = c;
338  for (uint8_t v_iter = 0; v_iter < c.size(); v_iter++) {
339  c.at(v_iter) += coeffs.at(1) * d.at(v_iter);
340  d.at(v_iter) += coeffs.at(2) * temp_c.at(v_iter);
341  }
342 
343  //% time updates for ec' and ed"
344  std::vector<float> temp_ec = ec;
345  for (int v_iter = 0; v_iter < n; v_iter++) {
346  ec.at(v_iter) += coeffs.at(1) * ed.at(v_iter);
347  ed.at(v_iter) += coeffs.at(2) * temp_ec.at(v_iter);
348  }
349 
350  std::vector<float> ecc, edd;
351  ecc.resize(n);
352  edd.resize(n);
353  for (int v_iter = 0; v_iter < n - 1; v_iter++) {
354  ecc.at(v_iter) = ec.at(v_iter + 1);
355  edd.at(v_iter) = ed.at(v_iter);
356  }
357 
358  if (rho_f <= 0 || rho_b <= 0) {
359  if (fIsDebug)
360  printf("PsdPronyFitter::ERROR2: prediction squared error was less "
361  "than or equal to zero\n");
362  return;
363  }
364  gam = 1 - ecc.at(n - k - 1); //n-k?
365  del = 1 - edd.front();
366  if (gam <= 0 || gam > 1 || del <= 0 || del > 1) {
367  if (fIsDebug)
368  printf("PsdPronyFitter::ERROR2: GAM or DEL gain factor not in "
369  "expected range 0 to 1\n");
370  return;
371  }
372 
373 
374  //% coefficients for next set of updates
375  coeffs.at(0) = ef.front();
376  coeffs.at(1) = eb.at(n - k - 1); //n-k?
377  coeffs.at(2) = coeffs.at(1) / rho_b;
378  coeffs.at(3) = coeffs.at(0) / rho_f;
379  coeffs.at(4) = coeffs.at(0) / del;
380  coeffs.at(5) = coeffs.at(1) / gam;
381 
382  //% order updates for c and d; time updates for a_f' and a_b"
383  std::vector<float> temp_ab = a_b;
384  std::reverse(std::begin(temp_ab), std::end(temp_ab));
385  std::reverse(std::begin(c), std::end(c));
386 
387  for (uint8_t i = 0; i < a_b.size(); i++)
388  a_b.at(i) += coeffs.at(5) * c.at(i);
389  std::reverse(std::begin(c), std::end(c));
390 
391  for (uint8_t i = 0; i < c.size(); i++)
392  c.at(i) += coeffs.at(2) * temp_ab.at(i);
393  c.push_back(coeffs.at(2));
394 
395  std::vector<float> temp_af2 = a_f;
396  for (uint8_t i = 0; i < a_f.size(); i++)
397  a_f.at(i) += coeffs.at(4) * d.at(i);
398 
399  for (uint8_t i = 0; i < d.size(); i++)
400  d.at(i) += coeffs.at(3) * temp_af2.at(i);
401  d.insert(d.begin(), coeffs.at(3));
402 
403  //% time updates for rho_f' and rho_b"
404  rho_f = rho_f - coeffs.at(4) * coeffs.at(0);
405  rho_b = rho_b - coeffs.at(5) * coeffs.at(1);
406 
407  if (rho_f <= 0 || rho_b <= 0) {
408  if (fIsDebug)
409  printf("PsdPronyFitter::ERROR3: prediction squared error was less "
410  "than or equal to zero\n");
411  return;
412  }
413 
414  //% order updates for ec and ed; time updates for ef' and eb"
415  for (int v_iter = 0; v_iter < n; v_iter++) {
416  ec.at(v_iter) = ecc.at(v_iter) + coeffs.at(2) * eb.at(v_iter);
417  eb.at(v_iter) = eb.at(v_iter) + coeffs.at(5) * ecc.at(v_iter);
418  ed.at(v_iter) = edd.at(v_iter) + coeffs.at(3) * ef.at(v_iter);
419  ef.at(v_iter) = ef.at(v_iter) + coeffs.at(4) * edd.at(v_iter);
420  }
421  }
422  }
423 
424  void PronyFitter::SetHarmonics(std::complex<float>* z) {
425  std::memcpy(fz, z, (fExpNumber + 1) * sizeof(std::complex<float>));
426  }
427 
428  void PronyFitter::SetExternalHarmonics(std::complex<float> z1,
429  std::complex<float> z2) {
430  std::complex<float>* z_arr = new std::complex<float>[fExpNumber + 1];
431  for (int i = 0; i <= fExpNumber; i++)
432  z_arr[i] = {0., 0.};
433  z_arr[0] = {1., 0.};
434  z_arr[1] = z1;
435  z_arr[2] = z2;
436  SetHarmonics(z_arr);
437  delete[] z_arr;
438  }
439 
440  std::complex<float>* PronyFitter::GetHarmonics() { return fz; }
441 
443 
444  void PronyFitter::MakePileUpRejection(int /*time_max*/) {}
445 
447  std::complex<float>** Zik_arr = new std::complex<float>*[fExpNumber + 1];
448  for (int i = 0; i < fExpNumber + 1; i++) {
449  Zik_arr[i] = new std::complex<float>[fExpNumber + 1];
450  for (int j = 0; j < fExpNumber + 1; j++)
451  Zik_arr[i][j] = {0., 0.};
452  }
453 
454  std::complex<float>* Zyk_arr = new std::complex<float>[fExpNumber + 1];
455  for (int i = 0; i < fExpNumber + 1; i++)
456  Zyk_arr[i] = {0., 0.};
457 
458  int samples_in_gate = fGateEnd - fSignalBegin + 1;
459  const std::complex<float> unit = {1., 0.};
460 
461  for (int i = 0; i <= fExpNumber; i++) {
462  for (int j = 0; j <= fExpNumber; j++) {
463  std::complex<float> temp = std::conj(fz[i]) * fz[j];
464  if (std::abs(temp - unit) > 1e-3) {
465  Zik_arr[i][j] = static_cast<std::complex<float>>(
466  (std::pow(temp, static_cast<float>(samples_in_gate)) - unit)
467  / (temp - unit));
468  } else {
469  Zik_arr[i][j] = static_cast<std::complex<float>>(samples_in_gate);
470  }
471  }
472  }
473 
474  std::complex<float>* z_power = new std::complex<float>[fExpNumber + 1];
475  for (int i = 0; i < fExpNumber + 1; i++)
476  z_power[i] = unit;
477 
478  for (int i = 0; i <= fExpNumber; i++) {
479  for (int sample_curr = fSignalBegin; sample_curr <= fGateEnd;
480  sample_curr++) {
481  Zyk_arr[i] += (std::complex<float>) (std::conj(z_power[i])
482  * (float) fuWfm.at(sample_curr));
483  z_power[i] *= fz[i];
484  }
485  }
486 
487  if (fIsDebug) {
488  printf("\nampl calculation\n");
489  for (int i = 0; i <= fExpNumber; i++) {
490  for (int j = 0; j <= fExpNumber; j++)
491  printf(
492  "%e%+ei ", std::real(Zik_arr[i][j]), std::imag(Zik_arr[i][j]));
493  printf(
494  " %e%+ei\n", std::real(Zyk_arr[i]), std::imag(Zyk_arr[i]));
495  }
496  }
497 
498  SolveSLEGauss(fh, Zik_arr, Zyk_arr, fExpNumber + 1);
499 
500  if (fIsDebug) {
501  printf("amplitudes\n%.0f%+.0fi ", std::real(fh[0]), std::imag(fh[0]));
502  for (int i = 1; i < fExpNumber + 1; i++)
503  printf("%e%+ei ", std::real(fh[i]), std::imag(fh[i]));
504  printf("\n\n");
505  }
506 
507  for (int i = 0; i < fExpNumber + 1; i++)
508  z_power[i] = unit;
509 
510  std::complex<float> fit_ampl_in_sample = {0., 0.};
511  fuFitZeroLevel = (uint16_t) std::real(fh[0]);
512  for (int sample_curr = 0; sample_curr < fSampleTotal; sample_curr++) {
513  fit_ampl_in_sample = {0., 0.};
514  if ((sample_curr >= fSignalBegin) && (sample_curr <= fGateEnd)) {
515  for (int i = 0; i < fExpNumber + 1; i++) {
516  fit_ampl_in_sample += fh[i] * z_power[i];
517  z_power[i] *= fz[i];
518  }
519  fuFitWfm.at(sample_curr) = (uint16_t) std::real(fit_ampl_in_sample);
520  } else
521  fuFitWfm.at(sample_curr) = fuFitZeroLevel;
522  }
523 
524  if (fIsDebug) {
525  printf("waveform:\n");
526  for (uint8_t i = 0; i < fuWfm.size(); i++)
527  printf("%u ", fuWfm.at(i));
528 
529  printf("\nfit waveform:\n");
530  for (uint8_t i = 0; i < fuFitWfm.size(); i++)
531  printf("%u ", fuFitWfm.at(i));
532 
533  printf("\nzero level %u\n\n", fuZeroLevel);
534  }
535 
536  for (int i = 0; i < fExpNumber + 1; i++)
537  delete[] Zik_arr[i];
538  delete[] Zik_arr;
539  delete[] Zyk_arr;
540  delete[] z_power;
541  }
542 
543 
544  std::complex<float>* PronyFitter::GetAmplitudes() { return fh; }
545 
546  float PronyFitter::GetIntegral(int gate_beg, int gate_end) {
547  float integral = 0.;
548  for (int sample_curr = gate_beg; sample_curr <= gate_end; sample_curr++)
549  integral += (float) fuFitWfm.at(sample_curr) - fuFitZeroLevel;
550 
551  if (std::isfinite(integral)) return integral;
552  return 0;
553  }
554 
555  uint16_t PronyFitter::GetFitValue(int sample_number) {
556  if (std::isfinite(fuFitWfm.at(sample_number)))
557  return fuFitWfm.at(sample_number);
558  return 0;
559  }
560 
562  std::complex<float> amplitude = {0., 0.};
563  if (x < GetSignalBeginFromPhase()) return std::real(fh[0]);
564  amplitude += fh[0];
565  for (int i = 1; i < fExpNumber + 1; i++)
566  amplitude += fh[i] * std::pow(fz[i], x - fSignalBegin);
567 
568  if (std::isfinite(std::real(amplitude))) return std::real(amplitude);
569  return 0;
570  }
571 
572  float PronyFitter::GetZeroLevel() { return (float) fuFitZeroLevel; }
573 
575  return fSignalBegin
576  + std::real(
577  (std::log(-fh[2] * std::log(fz[2])) - std::log(fh[1] * log(fz[1])))
578  / (std::log(fz[1]) - std::log(fz[2])));
579  }
580 
582  if (std::real(fh[2] / fh[1]) < 0)
583  return fSignalBegin
584  + std::real(std::log(-fh[2] / fh[1]) / std::log(fz[1] / fz[2]));
585  return -999.;
586  }
587 
589  return GetFitValue(GetSignalMaxTime());
590  }
591 
592  float PronyFitter::GetX(float level, int first_sample, int last_sample) {
593  int step = 0;
594  if (first_sample < last_sample)
595  step = 1;
596  else
597  step = -1;
598  float result_sample = 0.;
599  int sample_to_check = first_sample;
600  float amplitude = 0.;
601  float amplitude_prev = GetFitValue(sample_to_check - step);
602  while ((first_sample - sample_to_check) * (last_sample - sample_to_check)
603  <= 0) {
604  amplitude = GetFitValue(sample_to_check);
605  if ((level - amplitude) * (level - amplitude_prev) <= 0) {
606  result_sample = LevelBy2Points(sample_to_check,
607  amplitude,
608  sample_to_check - step,
609  amplitude_prev,
610  level);
611  return result_sample;
612  }
613  amplitude_prev = amplitude;
614  sample_to_check += step;
615  }
616 
617  return 0;
618  }
619 
620  float PronyFitter::GetX(float level,
621  int first_sample,
622  int last_sample,
623  float step) {
624  float result_sample = 0.;
625  float sample_to_check = (float) first_sample;
626  std::complex<float> amplitude = {0., 0.};
627  std::complex<float> amplitude_prev = GetFitValue(sample_to_check - step);
628  while ((first_sample - sample_to_check) * (last_sample - sample_to_check)
629  <= 0) {
630  amplitude = GetFitValue(sample_to_check);
631  if ((level - std::real(amplitude)) * (level - std::real(amplitude_prev))
632  <= 0) {
633  if (amplitude != amplitude_prev)
634  result_sample = LevelBy2Points(sample_to_check,
635  std::real(amplitude),
636  sample_to_check - step,
637  std::real(amplitude_prev),
638  level);
639  return result_sample;
640  }
641  amplitude_prev = amplitude;
642  sample_to_check += step;
643  }
644 
645  return 0;
646  }
647 
649  float Y1,
650  float X2,
651  float Y2,
652  float Y0) {
653  return (X1 * Y0 - X1 * Y2 - X2 * Y0 + X2 * Y1) / (Y1 - Y2);
654  }
655 
656  float PronyFitter::GetRSquare(int gate_beg, int gate_end) {
657  float R2 = 0.;
658  float RSS = 0.;
659  float TSS = 0.;
660  int m = gate_end - gate_beg + 1;
661  int params_number = 1 + 2 * fModelOrder;
662  if (m <= params_number) return 999;
663  float average = 0.;
664  for (int sample_curr = gate_beg; sample_curr <= gate_end; sample_curr++)
665  average += fuWfm.at(sample_curr);
666  average /= m;
667 
668  for (int sample_curr = gate_beg; sample_curr <= gate_end; sample_curr++) {
669  RSS += std::pow(fuFitWfm.at(sample_curr) - fuWfm.at(sample_curr), 2);
670  TSS += std::pow(fuWfm.at(sample_curr) - average, 2);
671  }
672  if (TSS == 0) return 999;
673  R2 =
674  RSS
675  / TSS; // correct definition is R2=1.-RSS/TSS, but R2=RSS/TSS is more convenient
676 
677  float R2_adj = R2 * (m - 1) / (m - params_number);
678  return R2_adj;
679  }
680 
683  }
684 
685  float PronyFitter::GetChiSquare(int gate_beg, int gate_end, int time_max) {
686  float chi2 = 0.;
687  int freedom_counter = 0;
688  int regions_number = 10;
689  float amplitude_max = std::abs(fuWfm.at(time_max) - fuZeroLevel);
690  if (amplitude_max == 0) return 999;
691 
692  int* probability_exp = new int[regions_number];
693  int* probability_theor = new int[regions_number];
694  float* amplitude_regions = new float[regions_number + 1];
695  amplitude_regions[0] = 0.;
696  for (int i = 0; i < regions_number; i++) {
697  probability_exp[i] = 0;
698  probability_theor[i] = 0;
699  amplitude_regions[i + 1] = (i + 1) * amplitude_max / regions_number;
700  }
701 
702  for (int sample_curr = gate_beg; sample_curr <= gate_end; sample_curr++) {
703  for (int i = 0; i < regions_number; i++) {
704  if ((std::abs(fuWfm.at(sample_curr) - fuZeroLevel)
705  > amplitude_regions[i])
706  && (std::abs(fuWfm.at(sample_curr) - fuZeroLevel)
707  <= amplitude_regions[i + 1]))
708  probability_exp[i]++;
709  if ((std::abs(fuFitWfm.at(sample_curr) - fuFitZeroLevel)
710  > amplitude_regions[i])
711  && (std::abs(fuFitWfm.at(sample_curr) - fuFitZeroLevel)
712  <= amplitude_regions[i + 1]))
713  probability_theor[i]++;
714  }
715  }
716 
717  for (int i = 0; i < regions_number; i++) {
718  if (probability_exp[i] > 0) {
719  chi2 += std::pow(probability_exp[i] - probability_theor[i], 2.)
720  / (probability_exp[i]);
721  freedom_counter++;
722  }
723  }
724 
725  if (freedom_counter > 0) chi2 /= freedom_counter;
726  delete[] probability_exp;
727  delete[] probability_theor;
728  delete[] amplitude_regions;
729 
730  return chi2;
731  }
732 
733  float PronyFitter::GetDeltaInSample(int sample) {
734  return fuFitWfm.at(sample) - fuWfm.at(sample);
735  }
736 
737  /*
738 void PronyFitter::DrawFit(TObjArray *check_fit_arr, TString hist_title)
739 {
740  float *sample_arr = new float[fSampleTotal];
741  for(int i = 0; i < fSampleTotal; i++)
742  sample_arr[i] = (float) i;
743 
744  TGraph* tgr_ptr = new TGraph( fSampleTotal, sample_arr, fuWfm);
745  TGraph* tgr_ptr_fit = new TGraph( fSampleTotal, sample_arr, fuFitWfm);
746  TCanvas *canv_ptr = new TCanvas(hist_title.Data());
747  tgr_ptr->SetTitle(hist_title.Data());
748  tgr_ptr->Draw();
749 
750  tgr_ptr_fit->SetLineColor(kRed);
751  tgr_ptr_fit->SetLineWidth(2);
752  tgr_ptr_fit->Draw("same");
753 
754  check_fit_arr->Add(canv_ptr);
755 
756  delete[] sample_arr;
757 }
758 */
759 
761  int last_sample) {
762  float best_R2 = 0.;
763  int best_signal_begin = 0;
764  bool IsReasonableRoot;
765  bool IsGoodFit = false;
766  int good_fit_counter = 0;
767 
768  for (int signal_begin = first_sample; signal_begin <= last_sample;
769  signal_begin++) {
770  SetSignalBegin(signal_begin);
772  IsReasonableRoot = true;
773  for (int j = 0; j < fExpNumber; j++)
774  IsReasonableRoot = IsReasonableRoot && (std::abs(fz[j + 1]) > 1e-6)
775  && (std::abs(fz[j + 1]) < 1e1);
776  IsGoodFit = (fTotalPolRoots > 0) && (IsReasonableRoot);
777 
778  if (IsGoodFit) {
779  if (fIsDebug)
780  printf("good fit candidate at signal begin %i\n", signal_begin);
781  good_fit_counter++;
783  float R2 = GetRSquare(fGateBeg, fGateEnd);
784  if (good_fit_counter == 1) {
785  best_R2 = R2;
786  best_signal_begin = signal_begin;
787  }
788  if (R2 < best_R2) {
789  best_R2 = R2;
790  best_signal_begin = signal_begin;
791  }
792  }
793  }
794 
795  return best_signal_begin;
796  }
797 
798  int PronyFitter::ChooseBestSignalBegin(int first_sample, int last_sample) {
799  float best_R2 = 0.;
800  int best_signal_begin = first_sample;
801 
802  for (int signal_begin = first_sample; signal_begin <= last_sample;
803  signal_begin++) {
804  SetSignalBegin(signal_begin);
806  float R2 = GetRSquare(fGateBeg, fGateEnd);
807  if (signal_begin == first_sample) best_R2 = R2;
808  if (R2 < best_R2) {
809  best_R2 = R2;
810  best_signal_begin = signal_begin;
811  }
812  }
813 
814  return best_signal_begin;
815  }
816 
817  void PronyFitter::SolveSLEGauss(float* x, float** r, float* b, int n) {
818  bool solvable = true;
819  int maxRow;
820  float maxEl, tmp, c;
821  float** a = new float*[n];
822  for (int i = 0; i < n; i++) {
823  a[i] = new float[n + 1];
824  for (int j = 0; j < n + 1; j++)
825  a[i][j] = 0.;
826  }
827 
828  for (int i = 0; i < n; i++) {
829  for (int j = 0; j < n; j++)
830  a[i][j] = r[i][j];
831  a[i][n] = b[i];
832  }
833 
834  for (int i = 0; i < n; i++) {
835  maxEl = std::abs(a[i][i]);
836  maxRow = i;
837  for (int k = i + 1; k < n; k++)
838  if (abs(a[k][i]) > maxEl) {
839  maxEl = std::abs(a[k][i]);
840  maxRow = k;
841  }
842 
843  if (maxEl == 0) {
844  solvable = false;
845  if (fIsDebug) printf("SLE has not solution\n");
846  }
847 
848  for (int k = i; k < n + 1; k++) {
849  tmp = a[maxRow][k];
850  a[maxRow][k] = a[i][k];
851  a[i][k] = tmp;
852  }
853 
854  for (int k = i + 1; k < n; k++) {
855  c = -a[k][i] / a[i][i];
856  for (int j = i; j < n + 1; j++) {
857  if (i == j)
858  a[k][j] = 0.;
859  else
860  a[k][j] += c * a[i][j];
861  }
862  }
863  }
864 
865  for (int i = n - 1; i >= 0; i--) {
866  x[i] = a[i][n] / a[i][i];
867  for (int k = i - 1; k >= 0; k--)
868  a[k][n] -= a[k][i] * x[i];
869  }
870 
871  if (!solvable) {
872  for (int i = n - 1; i >= 0; i--)
873  x[i] = 0.;
874  }
875 
876  for (int i = 0; i < n; i++)
877  delete[] a[i];
878  delete[] a;
879  }
880 
881  void PronyFitter::SolveSLEGauss(std::complex<float>* x,
882  std::complex<float>** r,
883  std::complex<float>* b,
884  int n) {
885  bool solvable = true;
886  int maxRow;
887  float maxEl;
888  std::complex<float> tmp, c;
889  std::complex<float>** a = new std::complex<float>*[n];
890  for (int i = 0; i < n; i++) {
891  a[i] = new std::complex<float>[n + 1];
892  for (int j = 0; j < n + 1; j++)
893  a[i][j] = {0., 0.};
894  }
895 
896  for (int i = 0; i < n; i++) {
897  for (int j = 0; j < n; j++)
898  a[i][j] = r[i][j];
899  a[i][n] = b[i];
900  }
901 
902  for (int i = 0; i < n; i++) {
903  maxEl = std::abs(a[i][i]);
904  maxRow = i;
905  for (int k = i + 1; k < n; k++)
906  if (std::abs(a[k][i]) > maxEl) {
907  maxEl = std::abs(a[k][i]);
908  maxRow = k;
909  }
910 
911  if (maxEl == 0) {
912  solvable = false;
913  if (fIsDebug) printf("PsdPronyFitter:: SLE has not solution\n");
914  }
915 
916  for (int k = i; k < n + 1; k++) {
917  tmp = a[maxRow][k];
918  a[maxRow][k] = a[i][k];
919  a[i][k] = tmp;
920  }
921 
922  for (int k = i + 1; k < n; k++) {
923  c = -a[k][i] / a[i][i];
924  for (int j = i; j < n + 1; j++) {
925  if (i == j)
926  a[k][j] = 0.;
927  else
928  a[k][j] += c * a[i][j];
929  }
930  }
931  }
932 
933  for (int i = n - 1; i >= 0; i--) {
934  x[i] = a[i][n] / a[i][i];
935  for (int k = i - 1; k >= 0; k--)
936  a[k][n] -= a[k][i] * x[i];
937  }
938 
939  if (!solvable) {
940  for (int i = n - 1; i >= 0; i--)
941  x[i] = {0., 0.};
942  }
943 
944  for (int i = 0; i < n; i++)
945  delete[] a[i];
946  delete[] a;
947  }
948 
949  void PronyFitter::SolveSLECholesky(float* x, float** a, float* b, int n) {
950  float temp;
951  float** u = new float*[n];
952  for (int i = 0; i < n; i++) {
953  u[i] = new float[n];
954  for (int j = 0; j < n; j++)
955  u[i][j] = 0.;
956  }
957 
958  float* y = new float[n];
959  for (int i = 0; i < n; i++)
960  y[i] = 0.;
961 
962  for (int i = 0; i < n; i++) {
963  temp = 0.;
964  for (int k = 0; k < i; k++)
965  temp = temp + u[k][i] * u[k][i];
966  u[i][i] = std::sqrt(a[i][i] - temp);
967  for (int j = i; j < n; j++) {
968  temp = 0.;
969  for (int k = 0; k < i; k++)
970  temp = temp + u[k][i] * u[k][j];
971  u[i][j] = (a[i][j] - temp) / u[i][i];
972  }
973  }
974 
975  for (int i = 0; i < n; i++) {
976  temp = 0.;
977  for (int k = 0; k < i; k++)
978  temp = temp + u[k][i] * y[k];
979  y[i] = (b[i] - temp) / u[i][i];
980  }
981 
982  for (int i = n - 1; i >= 0; i--) {
983  temp = 0.;
984  for (int k = i + 1; k < n; k++)
985  temp = temp + u[i][k] * x[k];
986  x[i] = (y[i] - temp) / u[i][i];
987  }
988 
989  for (int i = 0; i < n; i++)
990  delete[] u[i];
991  delete[] u;
992  delete[] y;
993  }
994 
996  delete[] fz;
997  delete[] fh;
998  }
999 
1001  fModelOrder = 0;
1002  fExpNumber = 0;
1003  fGateBeg = 0;
1004  fGateEnd = 0;
1005  fSampleTotal = 0;
1006  fuZeroLevel = 0.;
1007  fSignalBegin = 0;
1008  fTotalPolRoots = 0;
1009  fuFitWfm.clear();
1010  DeleteData();
1011  }
1012 
1013 
1014 } // namespace PsdSignalFitting
PsdSignalFitting::PronyFitter::fSignalBegin
int fSignalBegin
Definition: PronyFitter.h:98
PsdSignalFitting::PronyFitter::SetSignalBegin
void SetSignalBegin(int SignalBeg)
Definition: PronyFitter.cxx:48
PsdSignalFitting::PronyFitter::DeleteData
void DeleteData()
Definition: PronyFitter.cxx:995
PsdSignalFitting::PronyFitter::CovarianceDirect
void CovarianceDirect(float &rho_f, std::vector< float > &a_f, float &rho_b, std::vector< float > &a_b)
Definition: PronyFitter.cxx:122
PsdSignalFitting::PronyFitter::fuZeroLevel
uint16_t fuZeroLevel
Definition: PronyFitter.h:102
PsdSignalFitting::PronyFitter::SetHarmonics
void SetHarmonics(std::complex< float > *z)
Definition: PronyFitter.cxx:424
PronyFitter.h
PsdSignalFitting::PronyFitter::fGateBeg
int fGateBeg
Definition: PronyFitter.h:94
PsdSignalFitting::PronyFitter::fGateEnd
int fGateEnd
Definition: PronyFitter.h:95
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
PsdSignalFitting::PronyFitter::GetIntegral
float GetIntegral(int gate_beg, int gate_end)
Definition: PronyFitter.cxx:546
PsdSignalFitting::PronyFitter::GetRSquareSignal
float GetRSquareSignal()
Definition: PronyFitter.cxx:681
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
PsdSignalFitting::PronyFitter::ChooseBestSignalBeginHarmonics
int ChooseBestSignalBeginHarmonics(int first_sample, int last_sample)
Definition: PronyFitter.cxx:760
PsdSignalFitting::PronyFitter::GetDeltaInSample
float GetDeltaInSample(int sample)
Definition: PronyFitter.cxx:733
PsdSignalFitting::PronyFitter::Clear
void Clear()
Definition: PronyFitter.cxx:1000
PolynomComplexRoots.h
PsdSignalFitting::PronyFitter::GetNumberPolRoots
int GetNumberPolRoots()
Definition: PronyFitter.cxx:442
PsdSignalFitting::PronyFitter::fuFitZeroLevel
uint16_t fuFitZeroLevel
Definition: PronyFitter.h:106
PsdSignalFitting::PronyFitter::GetSignalMaxTime
float GetSignalMaxTime()
Definition: PronyFitter.cxx:574
polynomComplexRoots
void polynomComplexRoots(float *wr, float *wi, int n, float *a, int &numr)
Definition: PolynomComplexRoots.h:262
PsdSignalFitting::PronyFitter::fuFitWfm
std::vector< uint16_t > fuFitWfm
Definition: PronyFitter.h:105
PsdSignalFitting::PronyFitter::GetAmplitudes
std::complex< float > * GetAmplitudes()
Definition: PronyFitter.cxx:544
PsdSignalFitting::PronyFitter::GetFitValue
uint16_t GetFitValue(int sample_number)
Definition: PronyFitter.cxx:555
PsdSignalFitting::PronyFitter::SetWaveform
void SetWaveform(std::vector< uint16_t > &uWfm, uint16_t uZeroLevel)
Definition: PronyFitter.cxx:35
PsdSignalFitting::PronyFitter::LevelBy2Points
float LevelBy2Points(float X1, float Y1, float X2, float Y2, float Y0)
Definition: PronyFitter.cxx:648
z2
Double_t z2[nSects2]
Definition: pipe_v16a_mvdsts100.h:11
PsdSignalFitting::PronyFitter::GetRSquare
float GetRSquare(int gate_beg, int gate_end)
Definition: PronyFitter.cxx:656
PsdSignalFitting::PronyFitter::fTotalPolRoots
int fTotalPolRoots
Definition: PronyFitter.h:99
d
double d
Definition: P4_F64vec2.h:24
PsdSignalFitting::PronyFitter::CovarianceQRmod
void CovarianceQRmod(float &rho_f, std::vector< float > &a_f, float &rho_b, std::vector< float > &a_b)
Definition: PronyFitter.cxx:188
log
friend F32vec4 log(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:135
PsdSignalFitting::PronyFitter::fz
std::complex< float > * fz
Definition: PronyFitter.h:103
PsdSignalFitting::PronyFitter::fuWfm
std::vector< uint16_t > fuWfm
Definition: PronyFitter.h:101
PsdSignalFitting::PronyFitter::CalcSignalBegin
int CalcSignalBegin(float front_time_beg_03, float front_time_end)
Definition: PronyFitter.cxx:43
PsdSignalFitting::PronyFitter::fIsDebug
bool fIsDebug
Definition: PronyFitter.h:91
PsdSignalFitting::PronyFitter::GetZeroLevel
float GetZeroLevel()
Definition: PronyFitter.cxx:572
PsdSignalFitting::PronyFitter::fh
std::complex< float > * fh
Definition: PronyFitter.h:104
PsdSignalFitting::PronyFitter::SolveSLEGauss
void SolveSLEGauss(float *x, float **r, float *b, int n)
Definition: PronyFitter.cxx:817
PsdSignalFitting::PronyFitter::GetX
float GetX(float level, int first_sample, int last_sample)
Definition: PronyFitter.cxx:592
PsdSignalFitting::PronyFitter::Initialize
void Initialize(int model_order, int exponent_number, int gate_beg, int gate_end)
Definition: PronyFitter.cxx:15
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
PsdSignalFitting::PronyFitter::fExpNumber
int fExpNumber
Definition: PronyFitter.h:93
PolynomRealRoots.h
PsdSignalFitting::PronyFitter::ChooseBestSignalBegin
int ChooseBestSignalBegin(int first_sample, int last_sample)
Definition: PronyFitter.cxx:798
PsdSignalFitting::PronyFitter::fSampleTotal
int fSampleTotal
Definition: PronyFitter.h:96
PsdSignalFitting::PronyFitter::SolveSLECholesky
void SolveSLECholesky(float *x, float **a, float *b, int n)
Definition: PronyFitter.cxx:949
PsdSignalFitting::PronyFitter::CalculateFitAmplitudes
void CalculateFitAmplitudes()
Definition: PronyFitter.cxx:446
PsdSignalFitting::PronyFitter::AllocData
void AllocData()
Definition: PronyFitter.cxx:26
PsdSignalFitting::PronyFitter::GetHarmonics
std::complex< float > * GetHarmonics()
Definition: PronyFitter.cxx:440
PsdSignalFitting::PronyFitter::GetMaxAmplitude
float GetMaxAmplitude()
Definition: PronyFitter.cxx:588
z1
Double_t z1[nSects1]
Definition: pipe_v16a_mvdsts100.h:6
PsdSignalFitting::PronyFitter::CalculateFitHarmonics
void CalculateFitHarmonics()
Definition: PronyFitter.cxx:54
PsdSignalFitting::PronyFitter::PronyFitter
PronyFitter()
Definition: PronyFitter.h:23
PsdSignalFitting::PronyFitter::fModelOrder
int fModelOrder
Definition: PronyFitter.h:92
PsdSignalFitting::PronyFitter::SetExternalHarmonics
void SetExternalHarmonics(std::complex< float > z1, std::complex< float > z2)
Definition: PronyFitter.cxx:428
PsdSignalFitting
Definition: PronyFitter.cxx:5
PsdSignalFitting::PronyFitter::GetChiSquare
float GetChiSquare(int gate_beg, int gate_end, int time_max)
Definition: PronyFitter.cxx:685
PsdSignalFitting::PronyFitter::MakePileUpRejection
void MakePileUpRejection(int time_max)
Definition: PronyFitter.cxx:444
PsdSignalFitting::PronyFitter::GetSignalBeginFromPhase
float GetSignalBeginFromPhase()
Definition: PronyFitter.cxx:581