BmnRoot
Loading...
Searching...
No Matches
FitWLSQ.cxx
Go to the documentation of this file.
1#include <iomanip>
2#include <iostream>
3#include <float.h>
4#include <math.h>
5#include <assert.h>
6#include "FitWLSQ.h"
7
8using namespace std;
9
10FitWLSQ::FitWLSQ(Double_t *xx, Double_t cSigma1, Double_t cSigma2, Double_t cP1, Int_t nNPoints,
11 Int_t nNParams, Bool_t RBF, Bool_t LHf, Int_t RBfn) {
12 for (int i = 0; i < 3; i++) {
13 param[i] = 0;
14 for (int j = 0; j < 3; j++) {
15 cov[i][j] = 0;
16 }
17 }
18 RB_TYPE = RBfn;
19 if (RB_TYPE == 4 || RB_TYPE == 5) {
20 // Approximation-1
21 resolution = cSigma1;
22 } else {
23 // Approximation-2, Tukey weights, optimal weights
24 resolution = cSigma2;
25 }
26 Sigma1 = cSigma1;
27 Sigma2 = cSigma2;
28 P1 = cP1;
29 NParams = nNParams;
30 NPoints = nNPoints;
31 RB_FIT = RBF;
32 LH_INI = LHf;
33 numItera = 0;
34 // allocate internal dimensions
35 w = new Double_t[nNPoints];
36 wrb = new Double_t[nNPoints];
37 x = new Double_t[nNPoints];
38 // fill dimension with z coordinats
39 for (int i = 0; i < NPoints; i++) {
40 x[i] = xx[i];
41 w[i] = 0;
42 wrb[i] = 1;
43 }
44}
45
46FitWLSQ::~FitWLSQ() {
47 // delete internal dimensions
48 delete[] w;
49 delete[] wrb;
50 delete[] x;
51}
52
53inline Double_t FitWLSQ::det3(Double_t a[3][3]) {
54 Double_t val;
55 val = a[0][0] * a[1][1] * a[2][2]
56 + a[0][1] * a[1][2] * a[2][0]
57 + a[0][2] * a[1][0] * a[2][1]
58 - a[0][1] * a[1][0] * a[2][2]
59 - a[0][0] * a[1][2] * a[2][1]
60 - a[0][2] * a[1][1] * a[2][0];
61 return val;
62}
63
64Double_t FitWLSQ::alg_add(Double_t a[3][3], Int_t row, Int_t column) {
65 Int_t i, j, i_r, i_c;
66 Double_t alg_value = 0;
67 Double_t b[2][2];
68 Double_t bb;
69
70 if (NParams == 3) {
71 for (i = 0; i < 2; i++) {
72 for (j = 0; j < 2; j++) {
73 if (i < row)
74 i_r = i;
75 else
76 i_r = i + 1;
77 if (j < column)
78 i_c = j;
79 else
80 i_c = j + 1;
81
82 b[i][j] = a[i_r][i_c];
83 }
84 }
85 alg_value = pow(-1, row + column + 2) * (b[0][0] * b[1][1] - b[0][1] * b[1][0]);
86 }
87 if (NParams == 2) {
88 if (row < 1)
89 i_r = row + 1;
90 else
91 i_r = row - 1;
92
93 if (column < 1)
94 i_c = column + 1;
95 else
96 i_c = column - 1;
97
98 bb = a[i_r][i_c];
99 alg_value = pow(-1, row + column + 2) * bb;
100 }
101 return alg_value;
102}
103
104void FitWLSQ::cov_matrix(Double_t a[3][3]) {
105 Double_t dfact;
106 Int_t i, j;
107
108 // clean all matrixes
109 for (i = 0; i < 3; i++) {
110 for (j = 0; j < 3; j++) {
111 cov[i][j] = 0;
112 }
113 }
114
115 if (NParams == 3) {
116 dfact = 1 / det3(a);
117 for (i = 0; i < 3; i++) {
118 for (j = 0; j < 3; j++) {
119 cov[i][j] = dfact * alg_add(a, i, j);
120 }
121 }
122 }
123 if (NParams == 2) {
124 dfact = 1 / (a[0][0] * a[1][1] - a[0][1] * a[1][0]);
125 for (i = 0; i < 2; i++) {
126 for (j = 0; j < 2; j++) {
127 cov[i][j] = dfact * alg_add(a, i, j);
128 }
129 }
130 }
131 return;
132}
133
134Double_t FitWLSQ::f_lh(Double_t *yy, Double_t *zz, Double_t *wlh, Double_t *par) {
135 Double_t dcoef1, dcoef2, darg, deksp;
136 Double_t sum1, sum2, sum3;
137 Int_t npoint;
138 Int_t i;
139
140 dcoef1 = (Sigma2 * Sigma2 - Sigma1 * Sigma1);
141 dcoef1 = dcoef1 / (2. * Sigma1 * Sigma1 * Sigma2 * Sigma2);
142 dcoef2 = (1. - P1) * Sigma1 / (P1 * Sigma2);
143 sum1 = 0.;
144 sum2 = 0.;
145 npoint = 0;
146 for (i = 0; i < NPoints; i++) {
147 if (wlh[i]) {
148 if (NParams == 2)
149 darg = (par[1] * zz[i] + par[0] - yy[i]) * (par[1] * zz[i] + par[0] - yy[i]);
150
151 if (NParams == 3)
152 darg = (par[2] * zz[i] * zz[i] + par[1] * zz[i] + par[0] - yy[i]) *
153 (par[2] * zz[i] * zz[i] + par[1] * zz[i] + par[0] - yy[i]);
154
155 if ((dcoef1 * darg) < 600.) {
156 deksp = exp(dcoef1 * darg);
157 sum1 = sum1 - darg / (2. * Sigma1 * Sigma1);
158 sum2 = sum2 + log(1. + dcoef2 * deksp);
159 } else {
160 sum1 = sum1 - darg / (2. * Sigma1 * Sigma1);
161 sum2 = sum2 + log(dcoef2) + dcoef1 * darg;
162 }
163 }
164 npoint++;
165 }
166 sum3 = npoint * log(P1 / (Sigma1 * sqrt(2. * 3.141593)));
167
168 return -sum1 - sum2 - sum3;
169}
170
171Bool_t FitWLSQ::Initial(Double_t *trk) {
172 //unsigned char itera = 0;
173 Double_t fmin = 999999.;
174 Double_t f, rk1, rk2, rk3;
175 Double_t axmin[3] = {0};
176 Double_t axx[3] = {0};
177
178 // find initial approximation and fill dimesions to call fit iterations
179 for (int i = 0; i < NPoints; i++) {
180 w[i] = 1;
181 }
182 if (NPoints < 3) {
183 cout << "!!!Wrong number of points(<3) -> exit " << NPoints << endl;
184 return kFALSE;
185 }
186
187 if (NParams < 2 || NParams > 3) {
188 cout << "!!!Wrong number of parameters -> exit" << NParams << endl;
189 return kFALSE;
190 }
191
192 // LINE
193 if (NParams == 2) {
194 for (int i = 1; i < NPoints; i++) {
195 for (int j = 0; j < i; j++) {
196 if (w[i] && w[j]) {
197 axx[1] = (trk[i] - trk[j]) / (x[i] - x[j]);
198 axx[0] = (x[i] * trk[j] - x[j] * trk[i]) / (x[i] - x[j]);
199 f = f_lh(trk, x, w, axx);
200 if (f < fmin) {
201 fmin = f;
202 for (int k = 0; k < 3; k++) axmin[k] = axx[k];
203 }
204 }
205 }
206 }
207 }
208
209 // PARABOLA
210 if (NParams == 3) {
211 for (int i = 1; i < NPoints; i++) {
212 for (int j = 0; j < i; j++) {
213 for (int k = 0; k < j; k++) {
214 if (w[i] && w[j] && w[k]) {
215 rk1 = trk[i] / ((x[i] - x[j]) * (x[i] - x[k]));
216 rk2 = trk[j] / ((x[j] - x[i]) * (x[j] - x[k]));
217 rk3 = trk[k] / ((x[k] - x[i]) * (x[k] - x[j]));
218 axx[2] = rk1 + rk2 + rk3;
219 axx[1] = -rk1 * (x[j] + x[k]) - rk2 * (x[i] + x[k]) - rk3 * (x[i] + x[j]);
220 axx[0] = rk1 * x[j] * x[k] + rk2 * x[i] * x[k] + rk3 * x[i] * x[j];
221 f = f_lh(trk, x, w, axx);
222 if (f < fmin) {
223 fmin = f;
224 for (int n = 0; n < 3; n++) axmin[n] = axx[n];
225 }
226 }
227 }
228 }
229 }
230 }
231 for (int k = 0; k < 3; k++) param[k] = axmin[k];
232 if (!param[0] && !param[1] && !param[2]) {
233 cout << "!!!NOT VALID TRACK PARAMETERS!!!!!" << endl;
234 return kFALSE;
235 }
236 return kTRUE;
237}
238
239void FitWLSQ::CheckMatrix(Double_t *pdA, Double_t *pdB, TString sz) {
240 assert(pdA && pdB && NParams && sz);
241
242 //long l = cout.flags(ios::left);
243 cout << '\n' << sz;
244 cout << "\nMatrix A[n,n]:\n";
245 for (Int_t i = 0; i < NParams; i++) {
246 for (Int_t j = 0; j < NParams; j++)
247 cout << setw(12) << pdA[i + j * NParams];
248 cout << '\n';
249 }
250 cout << "Matrix B[n]:\n";
251 for (Int_t i = 0; i < NParams; i++)
252 cout << setw(12) << pdB[i];
253 //if (param) {
254 cout << "\nMatrix P[n]:\n";
255 for (Int_t i = 0; i < NParams; i++)
256 cout << setw(12) << param[i];
257 //}
258 cout << '\n' << flush;
259 //cout.flags( l );
260}
261
262Bool_t FitWLSQ::SymMatrix(Double_t *pdA, Double_t *pdB) {
263 assert(pdA && pdB && param && NParams);
264 Int_t n1 = NParams - 1;
265 for (Int_t i = 0u; i < n1; i++) {
266 Int_t in = i * NParams;
267 Double_t dPivot = -pdA[i + in];
268 Double_t dB = pdB[i];
269 Int_t i1 = i + 1;
270 for (Int_t j = i1; j < NParams; j++) {
271 Double_t dA = pdA[j + in];
272 Int_t jn = j * NParams;
273 pdA[jn] = dA / dPivot;
274 if (dA == 0.0)
275 continue;
276 for (Int_t k = i1; k <= j; k++) {
277 Int_t kn = k * NParams;
278 pdA[j + kn] += dA * pdA[kn];
279 }
280 pdB[j] += dB * pdA[jn];
281 }
282 }
283 Double_t dA = pdA[NParams * NParams - 1];
284 if (dA == 0.0) {
285 CheckMatrix(pdA, pdB, "matrix singularity!");
286 return kFALSE;
287 }
288 param[n1] = pdB[n1] /= dA;
289 Int_t ii = n1;
290 while (ii--) {
291 Double_t dB = -pdB[ii];
292 Int_t iin = ii * NParams;
293 for (Int_t j = ii + 1u; j < NParams; j++)
294 dB += pdA[j + iin] * param[j];
295 param[ii] = -dB / pdA[ii + iin];
296 }
297 return kTRUE;
298}
299
300Bool_t FitWLSQ::WLSQFit(Double_t *pdY) {
301 Double_t left[3][3];
302 if (NPoints < 2 || !NParams)
303 return kFALSE;
304 Int_t nn = NParams * NParams;
305 Int_t nnn = nn + NParams;
306 Double_t *pdA = new Double_t[nnn];
307 assert(pdA);
308
309 Double_t *pdB = pdA + nn;
310 for (Int_t i = 0u; i < nnn; i++)
311 pdA[i] = 0.0;
312 // calculating coefficients of the normal equation system AxP=B
313 for (Int_t n = 0u; n < NPoints; n++) {
314 Double_t dX = x[n];
315 Double_t dY = pdY[n];
316 dY *= wrb[n];
317 Double_t dW = 1;
318
319 for (Int_t i = 0u; i < NParams; i++) {
320 pdB[i] += dW * dY * pow(dX, i);
321 for (Int_t j = 0u; j <= i; j++) {
322 pdA[i + j * NParams] += dW * pow(dX, j + i);
323 }
324 }
325 }
326
327 // covariance matrix calculation
328 for (Int_t i = 0; i < NParams; i++) {
329 for (Int_t j = 0; j < NParams; j++) {
330 left[i][j] = pdA[i + j * NParams];
331 }
332 }
333
334 left[0][1] = left[1][0];
335 left[0][2] = left[2][0];
336 left[1][2] = left[2][1];
337
338 cov_matrix(left);
339
340 // solve equation
341 Bool_t b = SymMatrix(pdA, pdB);
342 delete[] pdA;
343 return b;
344}
345
346Double_t FitWLSQ::WLSQGet(Double_t dX) {
347 assert(param && NParams);
348
349 Double_t dY = param[0];
350 Double_t d = dX;
351 for (Int_t i = 1u; i < NParams; i++) {
352 dY += param[i] * d;
353 d *= dX;
354 }
355 return dY;
356}
357
358Double_t FitWLSQ::WLSQRms(Double_t *pdY) {
359 Double_t dChi2 = 0.0;
360 Double_t dSW = 0.0;
361 for (Int_t i = 0u; i < NPoints; i++) {
362 Double_t dX = x[i];
363 Double_t dY = pdY[i];
364 Double_t dW = wrb[i];
365
366 dSW += dW;
367 dY -= WLSQGet(dX);
368
369 dChi2 += dY * dY * dW;
370 }
371 return sqrt(dChi2 / dSW);
372}
373
374Bool_t FitWLSQ::par_check(Double_t *par, Double_t *par1) {
375 Int_t i;
376 Double_t val;
377 Bool_t value;
378
379 value = kTRUE;
380 for (i = 0; i < NParams; i++) {
381 val = (par1[i] - par[i]) / par1[i];
382 if (val < 0) val *= -1;
383 if (val > 0.001) value = kFALSE;
384 }
385
386 return value;
387}
388
389void FitWLSQ::RB(Double_t *yy, Int_t itera_rb) {
390 Int_t i;
391 Double_t di;
392 Double_t fitval;
393 Double_t sum_w = 0;
394 Double_t sum_wd = 0;
395 Double_t sigma_param;
396 Double_t dval;
397
398 //without annealing
399 if (RB_TYPE == 2)
400 sigma_param = 3;
401 // annealing
402 if (RB_TYPE == 3) {
403 sigma_param = 5 - itera_rb * 0.5;
404 if (sigma_param < 3) sigma_param = 3;
405 }
406
407 if (!itera_rb)
408 sigma_k = resolution;
409 if (sigma_k < resolution)
410 sigma_k = resolution;
411
412 // insert here switches
413 if (!RB_FIT || (!LH_INI && !itera_rb)) {
414 for (i = 0; i < NPoints; i++) {
415 wrb[i] = 1;
416 }
417 }
418 if ((LH_INI && RB_FIT) || (!LH_INI && itera_rb && RB_FIT)) {
419 for (i = 0; i < NPoints; i++) {
420 if (NParams == 2) fitval = param[0] + param[1] * x[i];
421 if (NParams == 3) fitval = param[0] + param[1] * x[i] + param[2] * x[i] * x[i];
422 di = yy[i] - fitval;
423 if (di < 0) di *= -1;
424
425 dval = 1 - di * di / (sigma_param * sigma_param * sigma_k * sigma_k);
426 if (di > sigma_param * sigma_k) wrb[i] = 0;
427 if (di < sigma_param * sigma_k) wrb[i] = dval * dval;
428 sum_w += wrb[i];
429 sum_wd += wrb[i] * di * di;
430 }
431 }
432
433 if ((LH_INI && RB_FIT) || (!LH_INI && itera_rb && RB_FIT))
434 sigma_k = sqrt((Double_t) sum_wd / sum_w);
435
436 return;
437}
438
439void FitWLSQ::RBN(Double_t *yy, Int_t itera_rb) {
440 Int_t i;
441 Double_t di;
442 Double_t fitval;
443 Double_t sum_w = 0;
444 Double_t sum_wd = 0;
445 Double_t sigma_p2;
446 Double_t sigma_p4;
447 Double_t sigma_pw;
448 Double_t dval;
449
450 // without annealing
451 if (RB_TYPE == 4) {
452 sigma_p2 = 21;
453 sigma_p4 = 3.5;
454 sigma_pw = 3;
455 }
456 //annealing
457 if (RB_TYPE == 5) {
458 //nastroit' po dannym dlya analiza!!!???
459 sigma_p2 = 25 - itera_rb * 0.2;
460 sigma_p4 = 5 - itera_rb * 0.2;
461 sigma_pw = 4 - itera_rb * 0.2;
462 if (sigma_p2 < 21) sigma_p2 = 21;
463 if (sigma_p4 < 3.5) sigma_p4 = 3.5;
464 if (sigma_pw < 3) sigma_pw = 3;
465 }
466
467 if (!itera_rb)
468 sigma_k = resolution;
469 if (sigma_k < resolution)
470 sigma_k = resolution;
471
472 if (!RB_FIT || (!LH_INI && !itera_rb)) {
473 for (i = 0; i < NPoints; i++) {
474 wrb[i] = 1;
475 }
476 }
477 if ((LH_INI && RB_FIT) || (!LH_INI && itera_rb && RB_FIT)) {
478 for (i = 0; i < NPoints; i++) {
479 if (NParams == 2) fitval = param[0] + param[1] * x[i];
480 if (NParams == 3) fitval = param[0] + param[1] * x[i] + param[2] * x[i] * x[i];
481 di = yy[i] - fitval;
482 if (di < 0) di *= -1;
483
484 if (di > sigma_p2 * sigma_k) wrb[i] = 0;
485 if (di < sigma_p2 * sigma_k && di > sigma_pw * sigma_k) {
486 dval = 1 - pow((di / (sigma_p2 * sigma_k)), 2);
487 wrb[i] = dval * dval * Sigma1 / (Sigma2 / P1);
488 }
489 if (di < sigma_pw * sigma_k) {
490 dval = 1 - pow((di / (sigma_p4 * sigma_k)), 4);
491 wrb[i] = dval * dval;
492 }
493
494 sum_w += wrb[i];
495 sum_wd += wrb[i] * di * di;
496 }
497 }
498
499 if ((LH_INI && RB_FIT) || (!LH_INI && itera_rb && RB_FIT))
500 sigma_k = sqrt((Double_t) sum_wd / sum_w);
501
502 return;
503}
504
505void FitWLSQ::RBM(Double_t *yy, Int_t itera_rb) {
506 // new weights definition
507 Int_t i;
508 Double_t di;
509 Double_t fitval;
510 Double_t sum_w = 0;
511 Double_t sum_wd = 0;
512 Double_t sigma_p;
513 Double_t dval;
514
515 // without annealing
516 if (RB_TYPE == 6)
517 sigma_p = 3.3;
518 // annealing
519 if (RB_TYPE == 7) {
520 sigma_p = 5. - itera_rb * 0.2;
521 if (sigma_p < 3.3) sigma_p = 3.3;
522 }
523
524 if (!itera_rb)
525 sigma_k = resolution;
526 if (sigma_k < resolution)
527 sigma_k = resolution;
528
529 if (!RB_FIT || (!LH_INI && !itera_rb)) {
530 for (i = 0; i < NPoints; i++) {
531 wrb[i] = 1;
532 }
533 }
534 if ((LH_INI && RB_FIT) || (!LH_INI && itera_rb && RB_FIT)) {
535 for (i = 0; i < NPoints; i++) {
536 if (NParams == 2) fitval = param[0] + param[1] * x[i];
537 if (NParams == 3) fitval = param[0] + param[1] * x[i] + param[2] * x[i] * x[i];
538 di = yy[i] - fitval;
539 if (di < 0) di *= -1;
540
541 if (di > sigma_p * sigma_k) wrb[i] = 0;
542 if (di < sigma_p * Sigma2 && di > sigma_p * Sigma1) {
543 dval = 1 - pow((di / (sigma_p * sigma_k)), 2);
544 wrb[i] = dval * dval * Sigma1 / (Sigma2 / P1);
545 }
546 if (di < sigma_p * Sigma1) {
547 dval = 1 - pow((di / (sigma_p * sigma_k)), 8);
548 wrb[i] = dval * dval;
549 }
550
551 sum_w += wrb[i];
552 sum_wd += wrb[i] * di * di;
553 }
554 }
555
556 if ((LH_INI && RB_FIT) || (!LH_INI && itera_rb && RB_FIT))
557 sigma_k = sqrt((Double_t) sum_wd / sum_w);
558
559 return;
560}
561
562void FitWLSQ::RB_OPT(Double_t *yy, Int_t itera_rb) {
563 Int_t i;
564 Double_t di;
565 Double_t fitval;
566 Double_t sum_w = 0;
567 Double_t sum_wd = 0;
568 Double_t sigma_param;
569 Double_t exp1, exp2, a, b, c, d;
570
571 //annealing
572 // without annealing
573 if (RB_TYPE == 0)
574 sigma_param = 3;
575 if (RB_TYPE == 1) {
576 sigma_param = 5 - itera_rb;
577 if (sigma_param < 3) sigma_param = 3;
578 }
579
580 if (!itera_rb)
581 sigma_k = resolution;
582 if (sigma_k < resolution)
583 sigma_k = resolution;
584
585 // insert here switches
586 if (!RB_FIT || (!LH_INI && !itera_rb)) {
587 for (i = 0; i < NPoints; i++) {
588 wrb[i] = 1;
589 }
590 }
591 if ((LH_INI && RB_FIT) || (!LH_INI && itera_rb && RB_FIT)) {
592 for (i = 0; i < NPoints; i++) {
593 if (NParams == 2) fitval = param[0] + param[1] * x[i];
594 if (NParams == 3) fitval = param[0] + param[1] * x[i] + param[2] * x[i] * x[i];
595 di = yy[i] - fitval;
596 if (di < 0) di *= -1;
597
598 if (di > sigma_param * sigma_k) wrb[i] = 0;
599 if (di < sigma_param * sigma_k) {
600 a = P1 * Sigma2 * Sigma2 * Sigma2;
601 b = (1 - P1) * Sigma1 * Sigma1 * Sigma1;
602 c = P1 * Sigma2;
603 d = (1 - P1) * Sigma1;
604 exp1 = exp((-1. * di * di) / (2. * Sigma1 * Sigma1));
605 exp2 = exp((-1. * di * di) / (2. * Sigma2 * Sigma2));
606 wrb[i] = ((a * exp1 + b * exp2) * (c + d)) / ((c * exp1 + d * exp2) * (a + b));
607 if (exp1 < 0.0001) {
608 wrb[i] = 0.01 - fabs(di * di * di);
609 }
610 }
611 sum_w += wrb[i];
612 sum_wd += wrb[i] * di * di;
613 }
614 }
615
616 if ((LH_INI && RB_FIT) || (!LH_INI && itera_rb && RB_FIT))
617 sigma_k = sqrt((Double_t) sum_wd / sum_w);
618
619 return;
620}
621
622Bool_t FitWLSQ::Fit(Double_t *y) {
623 Int_t itera;
624 Bool_t end_fit = kFALSE;
625 Bool_t stat;
626 Double_t par_c[3] = {0};
627 Double_t sumwrb;
628 Double_t *yyf;
629
630 Status=false;
631 yyf = new Double_t[NPoints];
632 if (LH_INI) {
633 stat = Initial(y);
634 if (!stat)
635 return kFALSE;
636 }
637 sigma_k = resolution;
638 itera = 0;
639 while (!end_fit) {
640 for (int i = 0; i < NParams; i++) par_c[i] = param[i];
641 if (RB_TYPE == 0 || RB_TYPE == 1)
642 RB_OPT(y, itera);
643 if (RB_TYPE == 2 || RB_TYPE == 3)
644 RB(y, itera);
645 if (RB_TYPE == 4 || RB_TYPE == 5)
646 RBN(y, itera);
647 if (RB_TYPE == 6 || RB_TYPE == 7)
648 RBM(y, itera);
649
650 // check number of points with non zero weights
651 //=============================================
652 sumwrb = 0;
653 for (int i = 0; i < NPoints; i++) {
654 if (wrb[i]) sumwrb++;
655 }
656 if (sumwrb < 4) {
657 //cout << "Problems with rb fit number of points with non zero weight: " << sumwrb << " " << itera << endl;
658 return kFALSE;
659 }
660 //=============================================
661
662 for (int i = 0; i < NPoints; i++) {
663 if (NParams == 3)
664 yyf[i] = y[i] - param[0] - param[1] * x[i] - param[2] * x[i] * x[i];
665 else
666 yyf[i] = y[i] - param[0] - param[1] * x[i];
667 }
668 stat = WLSQFit(yyf);
669
670 for (int i = 0; i < NParams; i++)
671 param[i] += par_c[i];
672 end_fit = par_check(param, par_c);
673 if (!stat) return kFALSE;
674 itera++;
675 if (itera > 50) {
676 //cout << " To many robust fit iterations -> exit" << endl;
677 return kFALSE;
678 }
679 }
680 numItera = itera;
681 //CovarianceMatrixPrint();
682 delete[] yyf;
683 Status=kTRUE;
684 return kTRUE;
685}
686
687void FitWLSQ::CovarianceMatrixPrint(void) {
688 //Double_t rmsval;
689 cout << "WLSQ Covariance Matrix" << endl;
690 for (int i = 0; i < 3; i++) {
691 for (int j = 0; j < 3; j++) {
692 cout << " " << cov[i][j];
693 }
694 cout << endl;
695 }
696 return;
697}
d
const Float_t d
Z-ccordinate of the first GEM-station.
Definition BmnMwpcHit.cxx:7
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition P4_F32vec4.h:34
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
i
int i
Definition P4_F32vec4.h:22
log
friend F32vec4 log(const F32vec4 &a)
Definition P4_F32vec4.h:123
exp
friend F32vec4 exp(const F32vec4 &a)
Definition P4_F32vec4.h:122
f
float f
Definition P4_F32vec4.h:21
FitWLSQ::wrb
Double_t * wrb
Definition FitWLSQ.h:19
FitWLSQ::FitWLSQ
FitWLSQ(Double_t *xx, Double_t cSigma1=0.01, Double_t cSigma2=0.07, Double_t cP1=0.9, Int_t nNPoints=6, Int_t nNParams=2, Bool_t RBF=kTRUE, Bool_t LHf=kTRUE, Int_t RBfn=7)
Definition FitWLSQ.cxx:10
FitWLSQ::Status
Bool_t Status
Definition FitWLSQ.h:20
FitWLSQ::WLSQRms
Double_t WLSQRms(Double_t *pdY)
Definition FitWLSQ.cxx:358
FitWLSQ::param
Double_t param[3]
Definition FitWLSQ.h:16
FitWLSQ::Fit
Bool_t Fit(Double_t *y)
Definition FitWLSQ.cxx:622
FitWLSQ::cov
Double_t cov[3][3]
Definition FitWLSQ.h:17
FitWLSQ::~FitWLSQ
~FitWLSQ()
Definition FitWLSQ.cxx:46
FitWLSQ::numItera
Int_t numItera
Definition FitWLSQ.h:18
std
STL namespace.