98 ring->
SetXYABP(-1., -1., -1., -1., -1.);
105 <<
"-E- CbmRichRingFitterEllipseTau::DoFit(), too many hits in the ring:"
107 ring->
SetXYABP(-1., -1., -1., -1., -1.);
130 TMatrixD PQ(5, 5, PQa);
147 TMatrixDEigen eig(PQ);
148 TMatrixD evc = eig.GetEigenVectors();
150 Double_t AlgParF = 0.;
151 AlgParF -= evc(0, 0) *
fM[
GA05];
154 AlgParF -= evc(1, 0) *
fM[
GA15];
157 AlgParF -= evc(2, 0) *
fM[
GA25];
160 AlgParF -= evc(3, 0) *
fM[
GA35];
163 AlgParF -= evc(4, 0) *
fM[
GA45];
170 const unsigned int numHits = ring->
GetNofHits();
171 const unsigned int numHits2 = 2 * numHits;
172 const unsigned int numHits3 = 3 * numHits;
173 const unsigned int numHits4 = 4 * numHits;
174 const unsigned int numHits5 = 5 * numHits;
175 const unsigned int numHits6 = 6 * numHits;
176 const double oneOverNumHits = 1. / numHits;
178 for (
unsigned int i = 0;
i < numHits;
i++) {
183 fZ[i6 + 1] =
fZT[
i + numHits] =
x *
y;
184 fZ[i6 + 2] =
fZT[
i + numHits2] =
y *
y;
185 fZ[i6 + 3] =
fZT[
i + numHits3] =
x;
186 fZ[i6 + 4] =
fZT[
i + numHits4] =
y;
187 fZ[i6 + 5] =
fZT[
i + numHits5] = 1.;
190 for (
unsigned int i = 0;
i < numHits6;
i++)
191 fM[
i] = oneOverNumHits *
fM[
i];
193 for (
int i = 0;
i < 25;
i++)
217 for (
int i = 0;
i < 25;
i++)
238 CalcChi2(Pxx, Pxy, Pyy, Px, Py, P, ring);
239 ring->
SetABCDEF(Pxx, Pxy, Pyy, Px, Py, P);
242 double QQx, QQy, Qxx, Qyy, Qx, Qy, Q;
243 double cosa, sina, cca, ssa, sin2a;
245 if (
fabs(Pxx - Pyy) > 0.1e-10) {
246 alpha = atan(Pxy / (Pxx - Pyy));
255 sin2a =
sin(2. * alpha);
256 Pxy = Pxy * sin2a / 2.;
257 Qx = Px * cosa + Py * sina;
258 Qy = -Px * sina + Py * cosa;
261 Qxx = 1. / (Pxx * cca + Pxy + Pyy * ssa);
262 Qyy = 1. / (Pxx * ssa - Pxy + Pyy * cca);
263 Q = -P + Qxx * QQx / 4. + Qyy * QQy / 4.;
267 double axisA = TMath::Sqrt(Q * Qxx);
268 double axisB = TMath::Sqrt(Q * Qyy);
269 double centerX = -xc * cosa / 2. + yc * sina / 2.;
270 double centerY = -xc * sina / 2. - yc * cosa / 2.;
272 ring->
SetXYABP(centerX, centerY, axisA, axisB, alpha);
322 const double det3_123_012 =
324 const double det3_123_013 =
326 const double det3_123_014 =
328 const double det3_123_023 =
330 const double det3_123_024 =
332 const double det3_123_034 =
334 const double det3_123_123 =
336 const double det3_123_124 =
338 const double det3_123_134 =
340 const double det3_123_234 =
342 const double det3_124_012 =
344 const double det3_124_013 =
346 const double det3_124_014 =
348 const double det3_124_023 =
350 const double det3_124_024 =
352 const double det3_124_034 =
354 const double det3_124_123 =
356 const double det3_124_124 =
358 const double det3_124_134 =
360 const double det3_124_234 =
362 const double det3_134_012 =
364 const double det3_134_013 =
366 const double det3_134_014 =
368 const double det3_134_023 =
370 const double det3_134_024 =
372 const double det3_134_034 =
374 const double det3_134_123 =
376 const double det3_134_124 =
378 const double det3_134_134 =
380 const double det3_134_234 =
382 const double det3_234_012 =
384 const double det3_234_013 =
386 const double det3_234_014 =
388 const double det3_234_023 =
390 const double det3_234_024 =
392 const double det3_234_034 =
394 const double det3_234_123 =
396 const double det3_234_124 =
398 const double det3_234_134 =
400 const double det3_234_234 =
404 const double det4_0123_0123 =
406 -
fP[
GM03] * det3_123_012;
407 const double det4_0123_0124 =
409 -
fP[
GM04] * det3_123_012;
410 const double det4_0123_0134 =
412 -
fP[
GM04] * det3_123_013;
413 const double det4_0123_0234 =
415 -
fP[
GM04] * det3_123_023;
416 const double det4_0123_1234 =
418 -
fP[
GM04] * det3_123_123;
419 const double det4_0124_0123 =
421 -
fP[
GM03] * det3_124_012;
422 const double det4_0124_0124 =
424 -
fP[
GM04] * det3_124_012;
425 const double det4_0124_0134 =
427 -
fP[
GM04] * det3_124_013;
428 const double det4_0124_0234 =
430 -
fP[
GM04] * det3_124_023;
431 const double det4_0124_1234 =
433 -
fP[
GM04] * det3_124_123;
434 const double det4_0134_0123 =
436 -
fP[
GM03] * det3_134_012;
437 const double det4_0134_0124 =
439 -
fP[
GM04] * det3_134_012;
440 const double det4_0134_0134 =
442 -
fP[
GM04] * det3_134_013;
443 const double det4_0134_0234 =
445 -
fP[
GM04] * det3_134_023;
446 const double det4_0134_1234 =
448 -
fP[
GM04] * det3_134_123;
449 const double det4_0234_0123 =
451 -
fP[
GM03] * det3_234_012;
452 const double det4_0234_0124 =
454 -
fP[
GM04] * det3_234_012;
455 const double det4_0234_0134 =
457 -
fP[
GM04] * det3_234_013;
458 const double det4_0234_0234 =
460 -
fP[
GM04] * det3_234_023;
461 const double det4_0234_1234 =
463 -
fP[
GM04] * det3_234_123;
464 const double det4_1234_0123 =
466 -
fP[
GM13] * det3_234_012;
467 const double det4_1234_0124 =
469 -
fP[
GM14] * det3_234_012;
470 const double det4_1234_0134 =
472 -
fP[
GM14] * det3_234_013;
473 const double det4_1234_0234 =
475 -
fP[
GM14] * det3_234_023;
476 const double det4_1234_1234 =
478 -
fP[
GM14] * det3_234_123;
481 const double det =
fP[
GM00] * det4_1234_1234 -
fP[
GM01] * det4_1234_0234
483 +
fP[
GM04] * det4_1234_0123;
491 const double oneOverDet = 1.0 / det;
492 const double mn1OverDet = -oneOverDet;
494 fP[
GM00] = det4_1234_1234 * oneOverDet;
495 fP[
GM01] = det4_0234_1234 * mn1OverDet;
496 fP[
GM02] = det4_0134_1234 * oneOverDet;
497 fP[
GM03] = det4_0124_1234 * mn1OverDet;
498 fP[
GM04] = det4_0123_1234 * oneOverDet;
500 fP[
GM10] = det4_1234_0234 * mn1OverDet;
501 fP[
GM11] = det4_0234_0234 * oneOverDet;
502 fP[
GM12] = det4_0134_0234 * mn1OverDet;
503 fP[
GM13] = det4_0124_0234 * oneOverDet;
504 fP[
GM14] = det4_0123_0234 * mn1OverDet;
506 fP[
GM20] = det4_1234_0134 * oneOverDet;
507 fP[
GM21] = det4_0234_0134 * mn1OverDet;
508 fP[
GM22] = det4_0134_0134 * oneOverDet;
509 fP[
GM23] = det4_0124_0134 * mn1OverDet;
510 fP[
GM24] = det4_0123_0134 * oneOverDet;
512 fP[
GM30] = det4_1234_0124 * mn1OverDet;
513 fP[
GM31] = det4_0234_0124 * oneOverDet;
514 fP[
GM32] = det4_0134_0124 * mn1OverDet;
515 fP[
GM33] = det4_0124_0124 * oneOverDet;
516 fP[
GM34] = det4_0123_0124 * mn1OverDet;
518 fP[
GM40] = det4_1234_0123 * oneOverDet;
519 fP[
GM41] = det4_0234_0123 * mn1OverDet;
520 fP[
GM42] = det4_0134_0123 * oneOverDet;
521 fP[
GM43] = det4_0124_0123 * mn1OverDet;
522 fP[
GM44] = det4_0123_0123 * oneOverDet;
528 const double*
const bp,
534 const double* arp0 = ap;
535 while (arp0 < ap + na) {
537 const double* bcp = bp;
543 while (bcp < bp + nb) {
544 cij += *arp++ * *bcp;
554 #define ROTATE(a, i, j, k, l) \
557 a[i][j] = g - s * (h + g * tau); \
558 a[k][l] = h + s * (g - h * tau)
563 double tresh, theta, tau, t, sm, s,
h, g, c;
566 unsigned int ip, iq,
i, j;
567 for (ip = 0; ip < 5; ip++) {
568 for (iq = 0; iq < 5; iq++)
573 for (ip = 0; ip < 5; ip++) {
582 for (sm = 0., ip = 0; ip < 5; ip++)
583 for (iq = ip + 1; iq < 5; iq++)
584 sm +=
fabs(a[ip][iq]);
585 if (sm == 0.) {
return; }
587 tresh = (
i < 4 ? 0.2 * sm / (5 * 5) : 0.);
589 for (ip = 0; ip < 4; ip++)
590 for (iq = ip + 1; iq < 5; iq++) {
592 g = 100. *
fabs(a[ip][iq]);
593 if (
i > 4 && (
float)
fabs(
d[ip] + g) == (
float)
fabs(
d[ip])
594 && (
float)
fabs(
d[iq] + g) == (
float)
fabs(
d[iq]))
597 else if (
fabs(a[ip][iq]) > tresh) {
603 theta = 0.5 *
h / a[ip][iq];
604 t = 1. / (
fabs(theta) +
sqrt(1. + theta * theta));
605 if (theta < 0.) t = -t;
607 c = 1. /
sqrt(1 + t * t);
616 for (j = 0; j < ip; j++) {
619 for (j = ip + 1; j < iq; j++) {
622 for (j = iq + 1; j < 5; j++) {
625 for (j = 0; j < 5; j++) {
631 for (ip = 0; ip < 5; ip++) {
642 for (
i = 0;
i < 5;
i++) {
644 for (j =
i + 1; j < 5; j++)
645 if (
d[j] >= p) p =
d[k = j];
649 for (j = 0; j < 5; j++) {