38 fNumberIterations = 1;
43 fStatisticType = kFitOptimChiCounts;
44 fAlphaOptim = kFitAlphaHalving;
46 fFitTaylor = kFitTaylorOrderFirst;
152 if (numberPeaks <= 0){
153 Error (
"TSpectrum2Fit",
"Invalid number of peaks, must be > than 0");
300 Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap =
313 c = c * t * (da1 + t * (da2 + t * da3));
325 Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap =
337 c = (-1.) * dap * c * t * t * (da1 + t * (2. * da2 + t * 3. * da3)) -
355 if (pw > 10) c *= a2;
356 if (pw > 12) c *= a2;
375 Double_t sk = 0,
b, lambdak, normk, normk_old = 0;
381 for (i = 0; i < size; i++) {
382 a[i][size + 2] = -a[i][size];
383 for (j = 0; j < size; j++) {
384 a[i][size + 2] += a[i][j] * a[j][size + 1];
386 normk += a[i][size + 2] * a[i][size + 2];
391 sk = normk / normk_old;
395 for (i = 0; i < size; i++) {
396 a[i][size + 3] = -a[i][size + 2] + sk * a[i][size + 3];
401 for (i = 0; i < size; i++) {
402 for (j = 0,
b = 0; j < size; j++) {
403 b += a[i][j] * a[j][size + 3];
405 lambdak +=
b * a[i][size + 3];
408 lambdak = normk / lambdak;
412 for (i = 0; i < size; i++)
413 a[i][size + 1] += lambdak * a[i][size + 3];
443 Double_t r, p,
r1,
e,
ex,
ey, vx, s2, px, py, rx, ry, erx, ery;
446 for (j = 0; j < numOfFittedPeaks; j++) {
447 p = (x - parameter[7 * j + 1]) / sigmax;
448 r = (y - parameter[7 * j + 2]) / sigmay;
450 e = (p * p - 2 * ro * p * r + r *
r) / (2 * (1 - ro * ro));
459 erx =
Erfc(p / s2 + 1 / (2 * bx)), ery =
460 Erfc(r / s2 + 1 / (2 * by));
461 ex = p / (s2 * bx), ey = r / (s2 * by);
463 px =
exp(ex) * erx, py =
exp(ey) * ery;
465 r1 += 0.5 * txy * px * py;
468 rx =
Erfc(p / s2), ry =
Erfc(r / s2);
469 r1 += 0.5 * sxy * rx * ry;
471 vx = vx + parameter[7 * j] *
r1;
473 p = (x - parameter[7 * j + 5]) / sigmax;
484 erx =
Erfc(p / s2 + 1 / (2 * bx));
495 vx = vx + parameter[7 * j + 3] *
r1;
497 r = (y - parameter[7 * j + 6]) / sigmay;
508 ery =
Erfc(r / s2 + 1 / (2 * by));
519 vx = vx + parameter[7 * j + 4] *
r1;
522 vx = vx + a0 + ax * x + ay *
y;
545 Double_t p,
r,
r1 = 0,
e,
ex,
ey, px, py, rx, ry, erx, ery, s2;
546 p = (x - x0) / sigmax;
547 r = (y - y0) / sigmay;
550 e = (p * p - 2 * ro * p * r + r *
r) / (2 * (1 - ro * ro));
559 erx =
Erfc(p / s2 + 1 / (2 * bx)), ery =
560 Erfc(r / s2 + 1 / (2 * by));
561 ex = p / (s2 * bx), ey = r / (s2 * by);
563 px =
exp(ex) * erx, py =
exp(ey) * ery;
565 r1 += 0.5 * txy * px * py;
568 rx =
Erfc(p / s2), ry =
Erfc(r / s2);
569 r1 += 0.5 * sxy * rx * ry;
592 p = (x - x0) / sigmax;
604 erx =
Erfc(p / s2 + 1 / (2 * bx));
640 Double_t p,
r,
r1 = 0,
e,
ex,
ey, px, py, rx, ry, erx, ery, s2;
641 p = (x - x0) / sigmax;
642 r = (y - y0) / sigmay;
645 e = (p * p - 2 * ro * p * r + r *
r) / (2 * (1 - ro * ro));
652 e = -(ro * r - p) / sigmax;
653 e =
e / (1 - ro * ro);
658 (-
Erfc(p / s2 + 1 / (2 * bx)) / (s2 * bx * sigmax) -
659 Derfc(p / s2 + 1 / (2 * bx)) / (s2 * sigmax)), ery =
660 Erfc(r / s2 + 1 / (2 * by));
661 ex = p / (s2 * bx), ey = r / (s2 * by);
663 px =
exp(ex) * erx, py =
exp(ey) * ery;
665 r1 += 0.5 * txy * px * py;
668 rx = -
Derfc(p / s2) / (s2 * sigmax), ry =
Erfc(r / s2);
669 r1 += 0.5 * sxy * rx * ry;
695 p = (x - x0) / sigmax;
696 r = (y - y0) / sigmay;
698 e = (p * p - 2 * ro * p * r + r *
r) / (2 * (1 - ro * ro));
705 e = -(ro * r - p) / sigmax;
706 e =
e / (1 - ro * ro);
707 r1 = r1 * (
e *
e - 1 / ((1 - ro * ro) * sigmax * sigmax));
735 Double_t p,
r,
r1 = 0,
e,
ex,
ey, px, py, rx, ry, erx, ery, s2;
736 p = (x - x0) / sigmax;
737 r = (y - y0) / sigmay;
740 e = (p * p - 2 * ro * p * r + r *
r) / (2 * (1 - ro * ro));
747 e = -(ro * p -
r) / sigmay;
748 e =
e / (1 - ro * ro);
753 (-
Erfc(r / s2 + 1 / (2 * by)) / (s2 * by * sigmay) -
754 Derfc(r / s2 + 1 / (2 * by)) / (s2 * sigmay)), erx =
755 Erfc(p / s2 + 1 / (2 * bx));
756 ex = p / (s2 * bx), ey = r / (s2 * by);
758 px =
exp(ex) * erx, py =
exp(ey) * ery;
760 r1 += 0.5 * txy * px * py;
763 ry = -
Derfc(r / s2) / (s2 * sigmay), rx =
Erfc(p / s2);
764 r1 += 0.5 * sxy * rx * ry;
790 p = (x - x0) / sigmax;
791 r = (y - y0) / sigmay;
793 e = (p * p - 2 * ro * p * r + r *
r) / (2 * (1 - ro * ro));
800 e = -(ro * p -
r) / sigmay;
801 e =
e / (1 - ro * ro);
802 r1 = r1 * (
e *
e - 1 / ((1 - ro * ro) * sigmay * sigmay));
825 p = (x - x0) / sigmax;
835 r1 = r1 * p / sigmax;
839 (-
Erfc(p / s2 + 1 / (2 * bx)) / (s2 * bx * sigmax) -
840 Derfc(p / s2 + 1 / (2 * bx)) / (s2 * sigmax));
847 rx = -
Derfc(p / s2) / (s2 * sigmax);
869 p = (x - x0) / sigmax;
878 r1 = r1 * (p * p / (sigmax * sigmax) - 1 / (sigmax * sigmax));
905 0,
e,
a,
b, x0, y0, s2, px, py, rx, ry, erx, ery,
ex,
ey;
908 for (j = 0; j < numOfFittedPeaks; j++) {
909 a = parameter[7 * j];
910 x0 = parameter[7 * j + 1];
911 y0 = parameter[7 * j + 2];
912 p = (x - x0) / sigmax;
913 r = (y - y0) / sigmay;
915 e = (p * p - 2 * ro * p * r + r *
r) / (2 * (1 - ro * ro));
922 b = -(ro * p * r - p * p) / sigmax;
923 e =
e * b / (1 - ro * ro);
927 -
Erfc(p / s2 + 1 / (2 * bx)) * p / (s2 * bx * sigmax) -
928 Derfc(p / s2 + 1 / (2 * bx)) * p / (s2 * sigmax), ery =
929 Erfc(r / s2 + 1 / (2 * by));
930 ex = p / (s2 * bx), ey = r / (s2 * by);
932 px =
exp(ex) * erx, py =
exp(ey) * ery;
934 e += 0.5 * txy * px * py;
937 rx = -
Derfc(p / s2) * p / (s2 * sigmax), ry =
Erfc(r / s2);
938 e += 0.5 * sxy * rx * ry;
943 x0 = parameter[7 * j + 5];
944 p = (x - x0) / sigmax;
952 e = 2 * b *
e / sigmax;
956 (-
Erfc(p / s2 + 1 / (2 * bx)) * p / (s2 * bx * sigmax) -
957 Derfc(p / s2 + 1 / (2 * bx)) * p / (s2 * sigmax));
964 rx = -
Derfc(p / s2) * p / (s2 * sigmax);
967 r1 += parameter[7 * j + 3] *
e;
992 for (j = 0; j < numOfFittedPeaks; j++) {
993 a = parameter[7 * j];
994 x0 = parameter[7 * j + 1];
995 y0 = parameter[7 * j + 2];
996 p = (x - x0) / sigmax;
997 r = (y - y0) / sigmay;
999 e = (p * p - 2 * ro * p * r + r *
r) / (2 * (1 - ro * ro));
1006 b = -(ro * p * r - p * p) / sigmax;
1007 e =
e * (b * b / (1 - ro * ro) -
1008 (3 * p * p - 2 * ro * p * r) / (sigmax * sigmax)) / (1 -
1015 x0 = parameter[7 * j + 5];
1016 p = (x - x0) / sigmax;
1024 e =
e * (4 * b * b - 6 *
b) / (sigmax * sigmax);
1025 r1 += parameter[7 * j + 3] *
e;
1052 0,
e,
a,
b, x0, y0, s2, px, py, rx, ry, erx, ery,
ex,
ey;
1055 for (j = 0; j < numOfFittedPeaks; j++) {
1056 a = parameter[7 * j];
1057 x0 = parameter[7 * j + 1];
1058 y0 = parameter[7 * j + 2];
1059 p = (x - x0) / sigmax;
1060 r = (y - y0) / sigmay;
1062 e = (p * p - 2 * ro * p * r + r *
r) / (2 * (1 - ro * ro));
1069 b = -(ro * p * r - r *
r) / sigmay;
1070 e =
e * b / (1 - ro * ro);
1074 -
Erfc(r / s2 + 1 / (2 * by)) * r / (s2 * by * sigmay) -
1075 Derfc(r / s2 + 1 / (2 * by)) * r / (s2 * sigmay), erx =
1076 Erfc(p / s2 + 1 / (2 * bx));
1077 ex = p / (s2 * bx), ey = r / (s2 * by);
1079 px =
exp(ex) * erx, py =
exp(ey) * ery;
1081 e += 0.5 * txy * px * py;
1084 ry = -
Derfc(r / s2) * r / (s2 * sigmay), rx =
Erfc(p / s2);
1085 e += 0.5 * sxy * rx * ry;
1090 y0 = parameter[7 * j + 6];
1091 r = (y - y0) / sigmay;
1099 e = 2 * b *
e / sigmay;
1103 (-
Erfc(r / s2 + 1 / (2 * by)) * r / (s2 * by * sigmay) -
1104 Derfc(r / s2 + 1 / (2 * by)) * r / (s2 * sigmay));
1111 ry = -
Derfc(r / s2) * r / (s2 * sigmay);
1114 r1 += parameter[7 * j + 4] *
e;
1139 for (j = 0; j < numOfFittedPeaks; j++) {
1140 a = parameter[7 * j];
1141 x0 = parameter[7 * j + 1];
1142 y0 = parameter[7 * j + 2];
1143 p = (x - x0) / sigmax;
1144 r = (y - y0) / sigmay;
1146 e = (p * p - 2 * ro * p * r + r *
r) / (2 * (1 - ro * ro));
1153 b = -(ro * p * r - r *
r) / sigmay;
1154 e =
e * (b * b / (1 - ro * ro) -
1155 (3 * r * r - 2 * ro * r * p) / (sigmay * sigmay)) / (1 -
1162 y0 = parameter[7 * j + 6];
1163 r = (y - y0) / sigmay;
1171 e =
e * (4 * b * b - 6 *
b) / (sigmay * sigmay);
1172 r1 += parameter[7 * j + 4] *
e;
1197 for (j = 0; j < numOfFittedPeaks; j++) {
1198 a = parameter[7 * j];
1199 x0 = parameter[7 * j + 1];
1200 y0 = parameter[7 * j + 2];
1204 rx = (px * px - 2 * r * px * qx + qx * qx);
1205 ex = rx / (2 * (1 - r *
r));
1212 tx = px * qx / (1 - r *
r);
1213 tx = tx - r * rx / ((1 - r *
r) * (1 - r * r));
1214 vx = vx + a * ex * tx;
1236 Double_t p,
r,
r1 = 0,
ex,
ey, px, py, erx, ery, s2, x0, y0,
a;
1239 for (j = 0; j < numOfFittedPeaks; j++) {
1240 a = parameter[7 * j];
1241 x0 = parameter[7 * j + 1];
1242 y0 = parameter[7 * j + 2];
1243 p = (x - x0) / sigmax;
1244 r = (y - y0) / sigmay;
1246 erx =
Erfc(p / s2 + 1 / (2 * bx)), ery =
1247 Erfc(r / s2 + 1 / (2 * by));
1248 ex = p / (s2 * bx), ey = r / (s2 * by);
1250 px =
exp(
ex) * erx, py =
exp(ey) * ery;
1252 r1 += 0.5 * a * px * py;
1275 for (j = 0; j < numOfFittedPeaks; j++) {
1276 a = parameter[7 * j];
1277 x0 = parameter[7 * j + 1];
1278 y0 = parameter[7 * j + 2];
1279 p = (x - x0) / sigmax;
1280 r = (y - y0) / sigmay;
1281 rx =
Erfc(p / s2), ry =
Erfc(r / s2);
1282 r1 += 0.5 * a * rx * ry;
1305 for (j = 0; j < numOfFittedPeaks; j++) {
1306 ax = parameter[7 * j + 3];
1307 x0 = parameter[7 * j + 5];
1308 p = (x - x0) / sigmax;
1310 erx =
Erfc(p / s2 + 1 / (2 * bx));
1315 r1 += 0.5 * ax * px;
1338 for (j = 0; j < numOfFittedPeaks; j++) {
1339 ax = parameter[7 * j + 4];
1340 x0 = parameter[7 * j + 6];
1341 p = (x - x0) / sigmax;
1343 erx =
Erfc(p / s2 + 1 / (2 * bx));
1348 r1 += 0.5 * ax * px;
1369 for (j = 0; j < numOfFittedPeaks; j++) {
1370 ax = parameter[7 * j + 3];
1371 x0 = parameter[7 * j + 5];
1372 p = (x - x0) / sigmax;
1375 r1 += 0.5 * ax * rx;
1396 for (j = 0; j < numOfFittedPeaks; j++) {
1397 ax = parameter[7 * j + 4];
1398 x0 = parameter[7 * j + 6];
1399 p = (x - x0) / sigmax;
1402 r1 += 0.5 * ax * rx;
1425 Double_t p,
r,
r1 = 0,
a, x0, y0, s2, px, py, erx, ery,
ex,
ey;
1428 for (j = 0; j < numOfFittedPeaks; j++) {
1429 a = parameter[7 * j];
1430 x0 = parameter[7 * j + 1];
1431 y0 = parameter[7 * j + 2];
1432 p = (x - x0) / sigmax;
1433 r = (y - y0) / sigmay;
1437 -
Erfc(p / s2 + 1 / (2 * bx)) * p / (s2 * bx * bx) -
1438 Derfc(p / s2 + 1 / (2 * bx)) / (s2 * bx * bx), ery =
1439 Erfc(r / s2 + 1 / (2 * by));
1440 ex = p / (s2 * bx), ey = r / (s2 * by);
1442 px =
exp(ex) * erx, py =
exp(ey) * ery;
1444 r1 += 0.5 *
a * txy * px * py;
1446 a = parameter[7 * j + 3];
1447 x0 = parameter[7 * j + 5];
1448 p = (x - x0) / sigmax;
1452 (-
Erfc(p / s2 + 1 / (2 * bx)) * p / (s2 * bx * bx) -
1453 Derfc(p / s2 + 1 / (2 * bx)) / (s2 * bx * bx));
1457 r1 += 0.5 *
a * tx * px;
1481 Double_t p,
r,
r1 = 0,
a, x0, y0, s2, px, py, erx, ery,
ex,
ey;
1484 for (j = 0; j < numOfFittedPeaks; j++) {
1485 a = parameter[7 * j];
1486 x0 = parameter[7 * j + 1];
1487 y0 = parameter[7 * j + 2];
1488 p = (x - x0) / sigmax;
1489 r = (y - y0) / sigmay;
1493 -
Erfc(r / s2 + 1 / (2 * by)) * r / (s2 * by * by) -
1494 Derfc(r / s2 + 1 / (2 * by)) / (s2 * by * by), erx =
1495 Erfc(p / s2 + 1 / (2 * bx));
1496 ex = p / (s2 * bx), ey = r / (s2 * by);
1498 px =
exp(ex) * erx, py =
exp(ey) * ery;
1500 r1 += 0.5 *
a * txy * px * py;
1502 a = parameter[7 * j + 4];
1503 y0 = parameter[7 * j + 6];
1504 r = (y - y0) / sigmay;
1508 (-
Erfc(r / s2 + 1 / (2 * by)) * r / (s2 * by * by) -
1509 Derfc(r / s2 + 1 / (2 * by)) / (s2 * by * by));
1513 r1 += 0.5 *
a * ty * py;
1537 r = 2 * a * pi * sx * sy *
r;
1559 r = 2 * pi * sx * sy *
r;
1582 r = a * 2 * pi * sy *
r;
1605 r = a * 2 * pi * sx *
r;
1628 r = -a * 2 * pi * sx * sy * ro /
r;
1851 Int_t i, i1, i2, j, k, shift =
1852 7 *
fNPeaks + 14, peak_vel, size, iter, pw,
1854 Double_t a,
b,
c, d = 0, alpha, chi_opt, yw, ywm,
f, chi2, chi_min, chi =
1855 0,
pi, pmin = 0, chi_cel = 0, chi_er;
1857 for (i = 0, j = 0; i <
fNPeaks; i++) {
1858 working_space[7 * i] =
fAmpInit[i];
1860 working_space[shift + j] =
fAmpInit[i];
1905 working_space[7 * i + 2] =
fRoInit;
1907 working_space[shift + j] =
fRoInit;
1910 working_space[7 * i + 3] =
fA0Init;
1912 working_space[shift + j] =
fA0Init;
1915 working_space[7 * i + 4] =
fAxInit;
1917 working_space[shift + j] =
fAxInit;
1920 working_space[7 * i + 5] =
fAyInit;
1922 working_space[shift + j] =
fAyInit;
1925 working_space[7 * i + 6] =
fTxyInit;
1927 working_space[shift + j] =
fTxyInit;
1930 working_space[7 * i + 7] =
fSxyInit;
1932 working_space[shift + j] =
fSxyInit;
1935 working_space[7 * i + 8] =
fTxInit;
1937 working_space[shift + j] =
fTxInit;
1940 working_space[7 * i + 9] =
fTyInit;
1942 working_space[shift + j] =
fTyInit;
1945 working_space[7 * i + 10] =
fSxyInit;
1947 working_space[shift + j] =
fSxInit;
1950 working_space[7 * i + 11] =
fSyInit;
1952 working_space[shift + j] =
fSyInit;
1955 working_space[7 * i + 12] =
fBxInit;
1957 working_space[shift + j] =
fBxInit;
1960 working_space[7 * i + 13] =
fByInit;
1962 working_space[shift + j] =
fByInit;
1967 for (j = 0; j < size; j++) {
1968 working_space[2 * shift + j] = 0, working_space[3 * shift + j] = 0;
1973 chi_opt = 0, pw =
fPower - 2;
1976 yw = source[i1][i2];
1978 f =
Shape2(fNPeaks, i1, i2,
1979 working_space, working_space[peak_vel],
1980 working_space[peak_vel + 1],
1981 working_space[peak_vel + 2],
1982 working_space[peak_vel + 3],
1983 working_space[peak_vel + 4],
1984 working_space[peak_vel + 5],
1985 working_space[peak_vel + 6],
1986 working_space[peak_vel + 7],
1987 working_space[peak_vel + 8],
1988 working_space[peak_vel + 9],
1989 working_space[peak_vel + 10],
1990 working_space[peak_vel + 11],
1991 working_space[peak_vel + 12],
1992 working_space[peak_vel + 13]);
2000 chi_opt += (yw -
f) * (yw - f) / ywm;
2020 for (j = 0, k = 0; j <
fNPeaks; j++) {
2023 working_space[7 * j + 1],
2024 working_space[7 * j + 2],
2025 working_space[peak_vel],
2026 working_space[peak_vel + 1],
2027 working_space[peak_vel + 2],
2028 working_space[peak_vel + 6],
2029 working_space[peak_vel + 7],
2030 working_space[peak_vel + 12],
2031 working_space[peak_vel + 13]);
2035 b = a * (yw * yw - f *
f) / (ywm * ywm);
2036 working_space[2 * shift + k] += b *
c;
2037 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2038 working_space[3 * shift + k] += b *
c;
2042 b = a * (yw -
f) / ywm;
2043 working_space[2 * shift + k] += b *
c;
2045 working_space[3 * shift + k] += b *
c;
2052 working_space[7 * j],
2053 working_space[7 * j + 1],
2054 working_space[7 * j + 2],
2055 working_space[peak_vel],
2056 working_space[peak_vel + 1],
2057 working_space[peak_vel + 2],
2058 working_space[peak_vel + 6],
2059 working_space[peak_vel + 7],
2060 working_space[peak_vel + 12],
2061 working_space[peak_vel + 13]);
2064 working_space[7 * j],
2065 working_space[7 * j + 1],
2066 working_space[7 * j + 2],
2067 working_space[peak_vel],
2068 working_space[peak_vel + 1],
2069 working_space[peak_vel + 2]);
2075 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
2084 b = a * (yw * yw - f *
f) / (ywm * ywm);
2085 working_space[2 * shift + k] += b *
c;
2086 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2087 working_space[3 * shift + k] += b *
c;
2091 b = a * (yw -
f) / ywm;
2092 working_space[2 * shift + k] += b *
c;
2094 working_space[3 * shift + k] += b *
c;
2101 working_space[7 * j],
2102 working_space[7 * j + 1],
2103 working_space[7 * j + 2],
2104 working_space[peak_vel],
2105 working_space[peak_vel + 1],
2106 working_space[peak_vel + 2],
2107 working_space[peak_vel + 6],
2108 working_space[peak_vel + 7],
2109 working_space[peak_vel + 12],
2110 working_space[peak_vel + 13]);
2113 working_space[7 * j],
2114 working_space[7 * j + 1],
2115 working_space[7 * j + 2],
2116 working_space[peak_vel],
2117 working_space[peak_vel + 1],
2118 working_space[peak_vel + 2]);
2124 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
2133 b = a * (yw * yw - f *
f) / (ywm * ywm);
2134 working_space[2 * shift + k] += b *
c;
2135 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2136 working_space[3 * shift + k] += b *
c;
2140 b = a * (yw -
f) / ywm;
2141 working_space[2 * shift + k] += b *
c;
2143 working_space[3 * shift + k] += b *
c;
2149 a =
Derampx(i1, working_space[7 * j + 5],
2150 working_space[peak_vel],
2151 working_space[peak_vel + 8],
2152 working_space[peak_vel + 10],
2153 working_space[peak_vel + 12]);
2157 b = a * (yw * yw - f *
f) / (ywm * ywm);
2158 working_space[2 * shift + k] += b *
c;
2159 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2160 working_space[3 * shift + k] += b *
c;
2164 b = a * (yw -
f) / ywm;
2165 working_space[2 * shift + k] += b *
c;
2167 working_space[3 * shift + k] += b *
c;
2173 a =
Derampx(i2, working_space[7 * j + 6],
2174 working_space[peak_vel + 1],
2175 working_space[peak_vel + 9],
2176 working_space[peak_vel + 11],
2177 working_space[peak_vel + 13]);
2181 b = a * (yw * yw - f *
f) / (ywm * ywm);
2182 working_space[2 * shift + k] += b *
c;
2183 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2184 working_space[3 * shift + k] += b *
c;
2188 b = a * (yw -
f) / ywm;
2189 working_space[2 * shift + k] += b *
c;
2191 working_space[3 * shift + k] += b *
c;
2197 a =
Deri01(i1, working_space[7 * j + 3],
2198 working_space[7 * j + 5],
2199 working_space[peak_vel],
2200 working_space[peak_vel + 8],
2201 working_space[peak_vel + 10],
2202 working_space[peak_vel + 12]);
2204 d =
Derderi01(i1, working_space[7 * j + 3],
2205 working_space[7 * j + 5],
2206 working_space[peak_vel]);
2212 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
2221 b = a * (yw * yw - f *
f) / (ywm * ywm);
2222 working_space[2 * shift + k] += b *
c;
2223 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2224 working_space[3 * shift + k] += b *
c;
2228 b = a * (yw -
f) / ywm;
2229 working_space[2 * shift + k] += b *
c;
2231 working_space[3 * shift + k] += b *
c;
2237 a =
Deri01(i2, working_space[7 * j + 4],
2238 working_space[7 * j + 6],
2239 working_space[peak_vel + 1],
2240 working_space[peak_vel + 9],
2241 working_space[peak_vel + 11],
2242 working_space[peak_vel + 13]);
2244 d =
Derderi01(i2, working_space[7 * j + 4],
2245 working_space[7 * j + 6],
2246 working_space[peak_vel + 1]);
2252 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
2261 b = a * (yw * yw - f *
f) / (ywm * ywm);
2262 working_space[2 * shift + k] += b *
c;
2263 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2264 working_space[3 * shift + k] += b *
c;
2268 b = a * (yw -
f) / ywm;
2269 working_space[2 * shift + k] += b *
c;
2271 working_space[3 * shift + k] += b *
c;
2279 working_space, working_space[peak_vel],
2280 working_space[peak_vel + 1],
2281 working_space[peak_vel + 2],
2282 working_space[peak_vel + 6],
2283 working_space[peak_vel + 7],
2284 working_space[peak_vel + 8],
2285 working_space[peak_vel + 10],
2286 working_space[peak_vel + 12],
2287 working_space[peak_vel + 13]);
2291 working_space[peak_vel],
2292 working_space[peak_vel + 1],
2293 working_space[peak_vel + 2]);
2299 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
2307 b = a * (yw * yw - f *
f) / (ywm * ywm);
2308 working_space[2 * shift + k] += b *
c;
2309 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2310 working_space[3 * shift + k] += b *
c;
2314 b = a * (yw -
f) / ywm;
2315 working_space[2 * shift + k] += b *
c;
2317 working_space[3 * shift + k] += b *
c;
2324 working_space, working_space[peak_vel],
2325 working_space[peak_vel + 1],
2326 working_space[peak_vel + 2],
2327 working_space[peak_vel + 6],
2328 working_space[peak_vel + 7],
2329 working_space[peak_vel + 9],
2330 working_space[peak_vel + 11],
2331 working_space[peak_vel + 12],
2332 working_space[peak_vel + 13]);
2336 working_space[peak_vel],
2337 working_space[peak_vel + 1],
2338 working_space[peak_vel + 2]);
2344 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
2352 b = a * (yw * yw - f *
f) / (ywm * ywm);
2353 working_space[2 * shift + k] += b *
c;
2354 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2355 working_space[3 * shift + k] += b *
c;
2359 b = a * (yw -
f) / ywm;
2360 working_space[2 * shift + k] += b *
c;
2362 working_space[3 * shift + k] += b *
c;
2368 a =
Derro(fNPeaks, i1, i2,
2369 working_space, working_space[peak_vel],
2370 working_space[peak_vel + 1],
2371 working_space[peak_vel + 2]);
2377 if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
2385 b = a * (yw * yw - f *
f) / (ywm * ywm);
2386 working_space[2 * shift + k] += b *
c;
2387 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2388 working_space[3 * shift + k] += b *
c;
2392 b = a * (yw -
f) / ywm;
2393 working_space[2 * shift + k] += b *
c;
2395 working_space[3 * shift + k] += b *
c;
2405 b = a * (yw * yw - f *
f) / (ywm * ywm);
2406 working_space[2 * shift + k] += b *
c;
2407 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2408 working_space[3 * shift + k] += b *
c;
2412 b = a * (yw -
f) / ywm;
2413 working_space[2 * shift + k] += b *
c;
2415 working_space[3 * shift + k] += b *
c;
2425 b = a * (yw * yw - f *
f) / (ywm * ywm);
2426 working_space[2 * shift + k] += b *
c;
2427 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2428 working_space[3 * shift + k] += b *
c;
2432 b = a * (yw -
f) / ywm;
2433 working_space[2 * shift + k] += b *
c;
2435 working_space[3 * shift + k] += b *
c;
2445 b = a * (yw * yw - f *
f) / (ywm * ywm);
2446 working_space[2 * shift + k] += b *
c;
2447 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2448 working_space[3 * shift + k] += b *
c;
2452 b = a * (yw -
f) / ywm;
2453 working_space[2 * shift + k] += b *
c;
2455 working_space[3 * shift + k] += b *
c;
2461 a =
Dertxy(fNPeaks, i1, i2,
2462 working_space, working_space[peak_vel],
2463 working_space[peak_vel + 1],
2464 working_space[peak_vel + 12],
2465 working_space[peak_vel + 13]);
2469 b = a * (yw * yw - f *
f) / (ywm * ywm);
2470 working_space[2 * shift + k] += b *
c;
2471 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2472 working_space[3 * shift + k] += b *
c;
2476 b = a * (yw -
f) / ywm;
2477 working_space[2 * shift + k] += b *
c;
2479 working_space[3 * shift + k] += b *
c;
2485 a =
Dersxy(fNPeaks, i1, i2,
2486 working_space, working_space[peak_vel],
2487 working_space[peak_vel + 1]);
2491 b = a * (yw * yw - f *
f) / (ywm * ywm);
2492 working_space[2 * shift + k] += b *
c;
2493 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2494 working_space[3 * shift + k] += b *
c;
2498 b = a * (yw -
f) / ywm;
2499 working_space[2 * shift + k] += b *
c;
2501 working_space[3 * shift + k] += b *
c;
2507 a =
Dertx(fNPeaks, i1, working_space,
2508 working_space[peak_vel],
2509 working_space[peak_vel + 12]);
2513 b = a * (yw * yw - f *
f) / (ywm * ywm);
2514 working_space[2 * shift + k] += b *
c;
2515 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2516 working_space[3 * shift + k] += b *
c;
2520 b = a * (yw -
f) / ywm;
2521 working_space[2 * shift + k] += b *
c;
2523 working_space[3 * shift + k] += b *
c;
2529 a =
Derty(fNPeaks, i2, working_space,
2530 working_space[peak_vel + 1],
2531 working_space[peak_vel + 13]);
2535 b = a * (yw * yw - f *
f) / (ywm * ywm);
2536 working_space[2 * shift + k] += b *
c;
2537 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2538 working_space[3 * shift + k] += b *
c;
2542 b = a * (yw -
f) / ywm;
2543 working_space[2 * shift + k] += b *
c;
2545 working_space[3 * shift + k] += b *
c;
2551 a =
Dersx(fNPeaks, i1, working_space,
2552 working_space[peak_vel]);
2556 b = a * (yw * yw - f *
f) / (ywm * ywm);
2557 working_space[2 * shift + k] += b *
c;
2558 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2559 working_space[3 * shift + k] += b *
c;
2563 b = a * (yw -
f) / ywm;
2564 working_space[2 * shift + k] += b *
c;
2566 working_space[3 * shift + k] += b *
c;
2572 a =
Dersy(fNPeaks, i2, working_space,
2573 working_space[peak_vel + 1]);
2577 b = a * (yw * yw - f *
f) / (ywm * ywm);
2578 working_space[2 * shift + k] += b *
c;
2579 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2580 working_space[3 * shift + k] += b *
c;
2584 b = a * (yw -
f) / ywm;
2585 working_space[2 * shift + k] += b *
c;
2587 working_space[3 * shift + k] += b *
c;
2593 a =
Derbx(fNPeaks, i1, i2,
2594 working_space, working_space[peak_vel],
2595 working_space[peak_vel + 1],
2596 working_space[peak_vel + 6],
2597 working_space[peak_vel + 8],
2598 working_space[peak_vel + 12],
2599 working_space[peak_vel + 13]);
2603 b = a * (yw * yw - f *
f) / (ywm * ywm);
2604 working_space[2 * shift + k] += b *
c;
2605 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2606 working_space[3 * shift + k] += b *
c;
2610 b = a * (yw -
f) / ywm;
2611 working_space[2 * shift + k] += b *
c;
2613 working_space[3 * shift + k] += b *
c;
2619 a =
Derby(fNPeaks, i1, i2,
2620 working_space, working_space[peak_vel],
2621 working_space[peak_vel + 1],
2622 working_space[peak_vel + 6],
2623 working_space[peak_vel + 8],
2624 working_space[peak_vel + 12],
2625 working_space[peak_vel + 13]);
2629 b = a * (yw * yw - f *
f) / (ywm * ywm);
2630 working_space[2 * shift + k] += b *
c;
2631 b = a * a * (4 * yw - 2 *
f) / (ywm * ywm);
2632 working_space[3 * shift + k] += b *
c;
2636 b = a * (yw -
f) / ywm;
2637 working_space[2 * shift + k] += b *
c;
2639 working_space[3 * shift + k] += b *
c;
2646 for (j = 0; j < size; j++) {
2647 if (
TMath::Abs(working_space[3 * shift + j]) > 0.000001)
2648 working_space[2 * shift + j] = working_space[2 * shift + j] /
TMath::Abs(working_space[3 * shift + j]);
2650 working_space[2 * shift + j] = 0;
2659 for (j = 0; j < size; j++) {
2660 working_space[4 * shift + j] = working_space[shift + j];
2666 chi_min = 10000 * chi2;
2669 chi_min = 0.1 * chi2;
2671 for (
pi = 0.1; flag == 0 &&
pi <= 100;
pi += 0.1) {
2672 for (j = 0; j < size; j++) {
2673 working_space[shift + j] = working_space[4 * shift + j] +
pi * alpha * working_space[2 * shift + j];
2675 for (i = 0, j = 0; i <
fNPeaks; i++) {
2677 if (working_space[shift + j] < 0)
2678 working_space[shift + j] = 0;
2679 working_space[7 * i] = working_space[shift + j];
2683 if (working_space[shift + j] <
fXmin)
2684 working_space[shift + j] =
fXmin;
2685 if (working_space[shift + j] >
fXmax)
2686 working_space[shift + j] =
fXmax;
2687 working_space[7 * i + 1] = working_space[shift + j];
2691 if (working_space[shift + j] <
fYmin)
2692 working_space[shift + j] =
fYmin;
2693 if (working_space[shift + j] >
fYmax)
2694 working_space[shift + j] =
fYmax;
2695 working_space[7 * i + 2] = working_space[shift + j];
2699 if (working_space[shift + j] < 0)
2700 working_space[shift + j] = 0;
2701 working_space[7 * i + 3] = working_space[shift + j];
2705 if (working_space[shift + j] < 0)
2706 working_space[shift + j] = 0;
2707 working_space[7 * i + 4] = working_space[shift + j];
2711 if (working_space[shift + j] <
fXmin)
2712 working_space[shift + j] =
fXmin;
2713 if (working_space[shift + j] >
fXmax)
2714 working_space[shift + j] =
fXmax;
2715 working_space[7 * i + 5] = working_space[shift + j];
2719 if (working_space[shift + j] <
fYmin)
2720 working_space[shift + j] =
fYmin;
2721 if (working_space[shift + j] >
fYmax)
2722 working_space[shift + j] =
fYmax;
2723 working_space[7 * i + 6] = working_space[shift + j];
2728 if (working_space[shift + j] < 0.001) {
2729 working_space[shift + j] = 0.001;
2731 working_space[peak_vel] = working_space[shift + j];
2735 if (working_space[shift + j] < 0.001) {
2736 working_space[shift + j] = 0.001;
2738 working_space[peak_vel + 1] = working_space[shift + j];
2742 if (working_space[shift + j] < -1) {
2743 working_space[shift + j] = -1;
2745 if (working_space[shift + j] > 1) {
2746 working_space[shift + j] = 1;
2748 working_space[peak_vel + 2] = working_space[shift + j];
2752 working_space[peak_vel + 3] = working_space[shift + j];
2756 working_space[peak_vel + 4] = working_space[shift + j];
2760 working_space[peak_vel + 5] = working_space[shift + j];
2764 working_space[peak_vel + 6] = working_space[shift + j];
2768 working_space[peak_vel + 7] = working_space[shift + j];
2772 working_space[peak_vel + 8] = working_space[shift + j];
2776 working_space[peak_vel + 9] = working_space[shift + j];
2780 working_space[peak_vel + 10] = working_space[shift + j];
2784 working_space[peak_vel + 11] = working_space[shift + j];
2788 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
2789 if (working_space[shift + j] < 0)
2790 working_space[shift + j] = -0.001;
2792 working_space[shift + j] = 0.001;
2794 working_space[peak_vel + 12] = working_space[shift + j];
2798 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
2799 if (working_space[shift + j] < 0)
2800 working_space[shift + j] = -0.001;
2802 working_space[shift + j] = 0.001;
2804 working_space[peak_vel + 13] = working_space[shift + j];
2810 yw = source[i1][i2];
2814 working_space[peak_vel],
2815 working_space[peak_vel + 1],
2816 working_space[peak_vel + 2],
2817 working_space[peak_vel + 3],
2818 working_space[peak_vel + 4],
2819 working_space[peak_vel + 5],
2820 working_space[peak_vel + 6],
2821 working_space[peak_vel + 7],
2822 working_space[peak_vel + 8],
2823 working_space[peak_vel + 9],
2824 working_space[peak_vel + 10],
2825 working_space[peak_vel + 11],
2826 working_space[peak_vel + 12],
2827 working_space[peak_vel + 13]);
2840 chi2 += (yw -
f) * (yw - f) / ywm;
2848 pmin =
pi, chi_min = chi2;
2858 for (j = 0; j < size; j++) {
2859 working_space[shift + j] = working_space[4 * shift + j] + pmin * alpha * working_space[2 * shift + j];
2861 for (i = 0, j = 0; i <
fNPeaks; i++) {
2863 if (working_space[shift + j] < 0)
2864 working_space[shift + j] = 0;
2865 working_space[7 * i] = working_space[shift + j];
2869 if (working_space[shift + j] <
fXmin)
2870 working_space[shift + j] =
fXmin;
2871 if (working_space[shift + j] >
fXmax)
2872 working_space[shift + j] =
fXmax;
2873 working_space[7 * i + 1] = working_space[shift + j];
2877 if (working_space[shift + j] <
fYmin)
2878 working_space[shift + j] =
fYmin;
2879 if (working_space[shift + j] >
fYmax)
2880 working_space[shift + j] =
fYmax;
2881 working_space[7 * i + 2] = working_space[shift + j];
2885 if (working_space[shift + j] < 0)
2886 working_space[shift + j] = 0;
2887 working_space[7 * i + 3] = working_space[shift + j];
2891 if (working_space[shift + j] < 0)
2892 working_space[shift + j] = 0;
2893 working_space[7 * i + 4] = working_space[shift + j];
2897 if (working_space[shift + j] <
fXmin)
2898 working_space[shift + j] =
fXmin;
2899 if (working_space[shift + j] >
fXmax)
2900 working_space[shift + j] =
fXmax;
2901 working_space[7 * i + 5] = working_space[shift + j];
2905 if (working_space[shift + j] <
fYmin)
2906 working_space[shift + j] =
fYmin;
2907 if (working_space[shift + j] >
fYmax)
2908 working_space[shift + j] =
fYmax;
2909 working_space[7 * i + 6] = working_space[shift + j];
2914 if (working_space[shift + j] < 0.001) {
2915 working_space[shift + j] = 0.001;
2917 working_space[peak_vel] = working_space[shift + j];
2921 if (working_space[shift + j] < 0.001) {
2922 working_space[shift + j] = 0.001;
2924 working_space[peak_vel + 1] = working_space[shift + j];
2928 if (working_space[shift + j] < -1) {
2929 working_space[shift + j] = -1;
2931 if (working_space[shift + j] > 1) {
2932 working_space[shift + j] = 1;
2934 working_space[peak_vel + 2] = working_space[shift + j];
2938 working_space[peak_vel + 3] = working_space[shift + j];
2942 working_space[peak_vel + 4] = working_space[shift + j];
2946 working_space[peak_vel + 5] = working_space[shift + j];
2950 working_space[peak_vel + 6] = working_space[shift + j];
2954 working_space[peak_vel + 7] = working_space[shift + j];
2958 working_space[peak_vel + 8] = working_space[shift + j];
2962 working_space[peak_vel + 9] = working_space[shift + j];
2966 working_space[peak_vel + 10] = working_space[shift + j];
2970 working_space[peak_vel + 11] = working_space[shift + j];
2974 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
2975 if (working_space[shift + j] < 0)
2976 working_space[shift + j] = -0.001;
2978 working_space[shift + j] = 0.001;
2980 working_space[peak_vel + 12] = working_space[shift + j];
2984 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
2985 if (working_space[shift + j] < 0)
2986 working_space[shift + j] = -0.001;
2988 working_space[shift + j] = 0.001;
2990 working_space[peak_vel + 13] = working_space[shift + j];
2998 for (j = 0; j < size; j++) {
2999 working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j];
3001 for (i = 0, j = 0; i <
fNPeaks; i++) {
3003 if (working_space[shift + j] < 0)
3004 working_space[shift + j] = 0;
3005 working_space[7 * i] = working_space[shift + j];
3009 if (working_space[shift + j] <
fXmin)
3010 working_space[shift + j] =
fXmin;
3011 if (working_space[shift + j] >
fXmax)
3012 working_space[shift + j] =
fXmax;
3013 working_space[7 * i + 1] = working_space[shift + j];
3017 if (working_space[shift + j] <
fYmin)
3018 working_space[shift + j] =
fYmin;
3019 if (working_space[shift + j] >
fYmax)
3020 working_space[shift + j] =
fYmax;
3021 working_space[7 * i + 2] = working_space[shift + j];
3025 if (working_space[shift + j] < 0)
3026 working_space[shift + j] = 0;
3027 working_space[7 * i + 3] = working_space[shift + j];
3031 if (working_space[shift + j] < 0)
3032 working_space[shift + j] = 0;
3033 working_space[7 * i + 4] = working_space[shift + j];
3037 if (working_space[shift + j] <
fXmin)
3038 working_space[shift + j] =
fXmin;
3039 if (working_space[shift + j] >
fXmax)
3040 working_space[shift + j] =
fXmax;
3041 working_space[7 * i + 5] = working_space[shift + j];
3045 if (working_space[shift + j] <
fYmin)
3046 working_space[shift + j] =
fYmin;
3047 if (working_space[shift + j] >
fYmax)
3048 working_space[shift + j] =
fYmax;
3049 working_space[7 * i + 6] = working_space[shift + j];
3054 if (working_space[shift + j] < 0.001) {
3055 working_space[shift + j] = 0.001;
3057 working_space[peak_vel] = working_space[shift + j];
3061 if (working_space[shift + j] < 0.001) {
3062 working_space[shift + j] = 0.001;
3064 working_space[peak_vel + 1] = working_space[shift + j];
3068 if (working_space[shift + j] < -1) {
3069 working_space[shift + j] = -1;
3071 if (working_space[shift + j] > 1) {
3072 working_space[shift + j] = 1;
3074 working_space[peak_vel + 2] = working_space[shift + j];
3078 working_space[peak_vel + 3] = working_space[shift + j];
3082 working_space[peak_vel + 4] = working_space[shift + j];
3086 working_space[peak_vel + 5] = working_space[shift + j];
3090 working_space[peak_vel + 6] = working_space[shift + j];
3094 working_space[peak_vel + 7] = working_space[shift + j];
3098 working_space[peak_vel + 8] = working_space[shift + j];
3102 working_space[peak_vel + 9] = working_space[shift + j];
3106 working_space[peak_vel + 10] = working_space[shift + j];
3110 working_space[peak_vel + 11] = working_space[shift + j];
3114 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
3115 if (working_space[shift + j] < 0)
3116 working_space[shift + j] = -0.001;
3118 working_space[shift + j] = 0.001;
3120 working_space[peak_vel + 12] = working_space[shift + j];
3124 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
3125 if (working_space[shift + j] < 0)
3126 working_space[shift + j] = -0.001;
3128 working_space[shift + j] = 0.001;
3130 working_space[peak_vel + 13] = working_space[shift + j];
3136 yw = source[i1][i2];
3138 f =
Shape2(fNPeaks, i1, i2,
3139 working_space, working_space[peak_vel],
3140 working_space[peak_vel + 1],
3141 working_space[peak_vel + 2],
3142 working_space[peak_vel + 3],
3143 working_space[peak_vel + 4],
3144 working_space[peak_vel + 5],
3145 working_space[peak_vel + 6],
3146 working_space[peak_vel + 7],
3147 working_space[peak_vel + 8],
3148 working_space[peak_vel + 9],
3149 working_space[peak_vel + 10],
3150 working_space[peak_vel + 11],
3151 working_space[peak_vel + 12],
3152 working_space[peak_vel + 13]);
3165 chi += (yw -
f) * (yw - f) / ywm;
3173 alpha = alpha * chi_opt / (2 * chi);
3176 alpha = alpha / 10.0;
3179 }
while (((chi > chi_opt
3184 for (j = 0; j < size; j++) {
3185 working_space[4 * shift + j] = 0;
3186 working_space[2 * shift + j] = 0;
3188 for (i1 =
fXmin, chi_cel = 0; i1 <=
fXmax; i1++) {
3190 yw = source[i1][i2];
3193 f =
Shape2(fNPeaks, i1, i2,
3194 working_space, working_space[peak_vel],
3195 working_space[peak_vel + 1],
3196 working_space[peak_vel + 2],
3197 working_space[peak_vel + 3],
3198 working_space[peak_vel + 4],
3199 working_space[peak_vel + 5],
3200 working_space[peak_vel + 6],
3201 working_space[peak_vel + 7],
3202 working_space[peak_vel + 8],
3203 working_space[peak_vel + 9],
3204 working_space[peak_vel + 10],
3205 working_space[peak_vel + 11],
3206 working_space[peak_vel + 12],
3207 working_space[peak_vel + 13]);
3208 chi_opt = (yw -
f) * (yw - f) / yw;
3209 chi_cel += (yw -
f) * (yw - f) / yw;
3212 for (j = 0, k = 0; j <
fNPeaks; j++) {
3215 working_space[7 * j + 1],
3216 working_space[7 * j + 2],
3217 working_space[peak_vel],
3218 working_space[peak_vel + 1],
3219 working_space[peak_vel + 2],
3220 working_space[peak_vel + 6],
3221 working_space[peak_vel + 7],
3222 working_space[peak_vel + 12],
3223 working_space[peak_vel + 13]);
3226 working_space[2 * shift + k] += chi_opt *
c;
3228 working_space[4 * shift + k] += b *
c;
3234 working_space[7 * j],
3235 working_space[7 * j + 1],
3236 working_space[7 * j + 2],
3237 working_space[peak_vel],
3238 working_space[peak_vel + 1],
3239 working_space[peak_vel + 2],
3240 working_space[peak_vel + 6],
3241 working_space[peak_vel + 7],
3242 working_space[peak_vel + 12],
3243 working_space[peak_vel + 13]);
3246 working_space[2 * shift + k] += chi_opt *
c;
3248 working_space[4 * shift + k] += b *
c;
3254 working_space[7 * j],
3255 working_space[7 * j + 1],
3256 working_space[7 * j + 2],
3257 working_space[peak_vel],
3258 working_space[peak_vel + 1],
3259 working_space[peak_vel + 2],
3260 working_space[peak_vel + 6],
3261 working_space[peak_vel + 7],
3262 working_space[peak_vel + 12],
3263 working_space[peak_vel + 13]);
3266 working_space[2 * shift + k] += chi_opt *
c;
3268 working_space[4 * shift + k] += b *
c;
3273 a =
Derampx(i1, working_space[7 * j + 5],
3274 working_space[peak_vel],
3275 working_space[peak_vel + 8],
3276 working_space[peak_vel + 10],
3277 working_space[peak_vel + 12]);
3280 working_space[2 * shift + k] += chi_opt *
c;
3282 working_space[4 * shift + k] += b *
c;
3287 a =
Derampx(i2, working_space[7 * j + 6],
3288 working_space[peak_vel + 1],
3289 working_space[peak_vel + 9],
3290 working_space[peak_vel + 11],
3291 working_space[peak_vel + 13]);
3294 working_space[2 * shift + k] += chi_opt *
c;
3296 working_space[4 * shift + k] += b *
c;
3301 a =
Deri01(i1, working_space[7 * j + 3],
3302 working_space[7 * j + 5],
3303 working_space[peak_vel],
3304 working_space[peak_vel + 8],
3305 working_space[peak_vel + 10],
3306 working_space[peak_vel + 12]);
3309 working_space[2 * shift + k] += chi_opt *
c;
3311 working_space[4 * shift + k] += b *
c;
3316 a =
Deri01(i2, working_space[7 * j + 4],
3317 working_space[7 * j + 6],
3318 working_space[peak_vel + 1],
3319 working_space[peak_vel + 9],
3320 working_space[peak_vel + 11],
3321 working_space[peak_vel + 13]);
3324 working_space[2 * shift + k] += chi_opt *
c;
3326 working_space[4 * shift + k] += b *
c;
3333 working_space, working_space[peak_vel],
3334 working_space[peak_vel + 1],
3335 working_space[peak_vel + 2],
3336 working_space[peak_vel + 6],
3337 working_space[peak_vel + 7],
3338 working_space[peak_vel + 8],
3339 working_space[peak_vel + 10],
3340 working_space[peak_vel + 12],
3341 working_space[peak_vel + 13]);
3344 working_space[2 * shift + k] += chi_opt *
c;
3346 working_space[4 * shift + k] += b *
c;
3352 working_space, working_space[peak_vel],
3353 working_space[peak_vel + 1],
3354 working_space[peak_vel + 2],
3355 working_space[peak_vel + 6],
3356 working_space[peak_vel + 7],
3357 working_space[peak_vel + 9],
3358 working_space[peak_vel + 11],
3359 working_space[peak_vel + 12],
3360 working_space[peak_vel + 13]);
3363 working_space[2 * shift + k] += chi_opt *
c;
3365 working_space[4 * shift + k] += b *
c;
3370 a =
Derro(fNPeaks, i1, i2,
3371 working_space, working_space[peak_vel],
3372 working_space[peak_vel + 1],
3373 working_space[peak_vel + 2]);
3376 working_space[2 * shift + k] += chi_opt *
c;
3378 working_space[4 * shift + k] += b *
c;
3386 working_space[2 * shift + k] += chi_opt *
c;
3388 working_space[4 * shift + k] += b *
c;
3396 working_space[2 * shift + k] += chi_opt *
c;
3398 working_space[4 * shift + k] += b *
c;
3406 working_space[2 * shift + k] += chi_opt *
c;
3408 working_space[4 * shift + k] += b *
c;
3413 a =
Dertxy(fNPeaks, i1, i2,
3414 working_space, working_space[peak_vel],
3415 working_space[peak_vel + 1],
3416 working_space[peak_vel + 12],
3417 working_space[peak_vel + 13]);
3420 working_space[2 * shift + k] += chi_opt *
c;
3422 working_space[4 * shift + k] += b *
c;
3427 a =
Dersxy(fNPeaks, i1, i2,
3428 working_space, working_space[peak_vel],
3429 working_space[peak_vel + 1]);
3432 working_space[2 * shift + k] += chi_opt *
c;
3434 working_space[4 * shift + k] += b *
c;
3439 a =
Dertx(fNPeaks, i1, working_space,
3440 working_space[peak_vel],
3441 working_space[peak_vel + 12]);
3444 working_space[2 * shift + k] += chi_opt *
c;
3446 working_space[4 * shift + k] += b *
c;
3451 a =
Derty(fNPeaks, i2, working_space,
3452 working_space[peak_vel + 1],
3453 working_space[peak_vel + 13]);
3456 working_space[2 * shift + k] += chi_opt *
c;
3458 working_space[4 * shift + k] += b *
c;
3463 a =
Dersx(fNPeaks, i1, working_space,
3464 working_space[peak_vel]);
3467 working_space[2 * shift + k] += chi_opt *
c;
3469 working_space[4 * shift + k] += b *
c;
3474 a =
Dersy(fNPeaks, i2, working_space,
3475 working_space[peak_vel + 1]);
3478 working_space[2 * shift + k] += chi_opt *
c;
3480 working_space[4 * shift + k] += b *
c;
3485 a =
Derbx(fNPeaks, i1, i2,
3486 working_space, working_space[peak_vel],
3487 working_space[peak_vel + 1],
3488 working_space[peak_vel + 6],
3489 working_space[peak_vel + 8],
3490 working_space[peak_vel + 12],
3491 working_space[peak_vel + 13]);
3494 working_space[2 * shift + k] += chi_opt *
c;
3496 working_space[4 * shift + k] += b *
c;
3501 a =
Derby(fNPeaks, i1, i2,
3502 working_space, working_space[peak_vel],
3503 working_space[peak_vel + 1],
3504 working_space[peak_vel + 6],
3505 working_space[peak_vel + 8],
3506 working_space[peak_vel + 12],
3507 working_space[peak_vel + 13]);
3510 working_space[2 * shift + k] += chi_opt *
c;
3512 working_space[4 * shift + k] += b *
c;
3520 chi_er = chi_cel /
b;
3521 for (i = 0, j = 0; i <
fNPeaks; i++) {
3523 Volume(working_space[7 * i], working_space[peak_vel],
3524 working_space[peak_vel + 1], working_space[peak_vel + 2]);
3528 a =
Derpa2(working_space[peak_vel],
3529 working_space[peak_vel + 1],
3530 working_space[peak_vel + 2]);
3531 b = working_space[4 * shift + j];
3541 working_space[peak_vel + 1],
3542 working_space[peak_vel + 2]);
3543 b = working_space[4 * shift + peak_vel];
3553 working_space[peak_vel],
3554 working_space[peak_vel + 2]);
3555 b = working_space[4 * shift + peak_vel + 1];
3564 a =
Derpro(working_space[shift + j], working_space[peak_vel],
3565 working_space[peak_vel + 1],
3566 working_space[peak_vel + 2]);
3567 b = working_space[4 * shift + peak_vel + 2];
3582 fAmpCalc[i] = working_space[shift + j];
3583 if (working_space[3 * shift + j] != 0)
3594 if (working_space[3 * shift + j] != 0)
3605 if (working_space[3 * shift + j] != 0)
3616 if (working_space[3 * shift + j] != 0)
3627 if (working_space[3 * shift + j] != 0)
3638 if (working_space[3 * shift + j] != 0)
3649 if (working_space[3 * shift + j] != 0)
3661 if (working_space[3 * shift + j] != 0)
3672 if (working_space[3 * shift + j] != 0)
3682 fRoCalc = working_space[shift + j];
3683 if (working_space[3 * shift + j] != 0)
3693 fA0Calc = working_space[shift + j];
3694 if (working_space[3 * shift + j] != 0)
3704 fAxCalc = working_space[shift + j];
3705 if (working_space[3 * shift + j] != 0)
3715 fAyCalc = working_space[shift + j];
3716 if (working_space[3 * shift + j] != 0)
3726 fTxyCalc = working_space[shift + j];
3727 if (working_space[3 * shift + j] != 0)
3737 fSxyCalc = working_space[shift + j];
3738 if (working_space[3 * shift + j] != 0)
3748 fTxCalc = working_space[shift + j];
3749 if (working_space[3 * shift + j] != 0)
3759 fTyCalc = working_space[shift + j];
3760 if (working_space[3 * shift + j] != 0)
3770 fSxCalc = working_space[shift + j];
3771 if (working_space[3 * shift + j] != 0)
3781 fSyCalc = working_space[shift + j];
3782 if (working_space[3 * shift + j] != 0)
3792 fBxCalc = working_space[shift + j];
3793 if (working_space[3 * shift + j] != 0)
3803 fByCalc = working_space[shift + j];
3804 if (working_space[3 * shift + j] != 0)
3817 f =
Shape2(fNPeaks, i1, i2,
3818 working_space, working_space[peak_vel],
3819 working_space[peak_vel + 1],
3820 working_space[peak_vel + 2],
3821 working_space[peak_vel + 3],
3822 working_space[peak_vel + 4],
3823 working_space[peak_vel + 5],
3824 working_space[peak_vel + 6],
3825 working_space[peak_vel + 7],
3826 working_space[peak_vel + 8],
3827 working_space[peak_vel + 9],
3828 working_space[peak_vel + 10],
3829 working_space[peak_vel + 11],
3830 working_space[peak_vel + 12],
3831 working_space[peak_vel + 13]);
3835 delete [] working_space;
3945 Int_t i, i1, i2, j, k, shift =
3946 7 *
fNPeaks + 14, peak_vel, size, iter, regul_cycle,
3948 Double_t a,
b,
c, alpha, chi_opt, yw, ywm,
f, chi2, chi_min, chi = 0
3949 ,
pi, pmin = 0, chi_cel = 0, chi_er;
3951 for (i = 0, j = 0; i <
fNPeaks; i++) {
3952 working_space[7 * i] =
fAmpInit[i];
3954 working_space[shift + j] =
fAmpInit[i];
3999 working_space[7 * i + 2] =
fRoInit;
4001 working_space[shift + j] =
fRoInit;
4004 working_space[7 * i + 3] =
fA0Init;
4006 working_space[shift + j] =
fA0Init;
4009 working_space[7 * i + 4] =
fAxInit;
4011 working_space[shift + j] =
fAxInit;
4014 working_space[7 * i + 5] =
fAyInit;
4016 working_space[shift + j] =
fAyInit;
4019 working_space[7 * i + 6] =
fTxyInit;
4021 working_space[shift + j] =
fTxyInit;
4024 working_space[7 * i + 7] =
fSxyInit;
4026 working_space[shift + j] =
fSxyInit;
4029 working_space[7 * i + 8] =
fTxInit;
4031 working_space[shift + j] =
fTxInit;
4034 working_space[7 * i + 9] =
fTyInit;
4036 working_space[shift + j] =
fTyInit;
4039 working_space[7 * i + 10] =
fSxyInit;
4041 working_space[shift + j] =
fSxInit;
4044 working_space[7 * i + 11] =
fSyInit;
4046 working_space[shift + j] =
fSyInit;
4049 working_space[7 * i + 12] =
fBxInit;
4051 working_space[shift + j] =
fBxInit;
4054 working_space[7 * i + 13] =
fByInit;
4056 working_space[shift + j] =
fByInit;
4061 for (i = 0; i < size; i++)
4062 working_matrix[i] =
new Double_t[size + 4];
4064 for (j = 0; j < size; j++) {
4065 working_space[3 * shift + j] = 0;
4066 for (k = 0; k <= size; k++) {
4067 working_matrix[j][k] = 0;
4077 for (j = 0, k = 0; j <
fNPeaks; j++) {
4079 working_space[2 * shift + k] =
4081 working_space[7 * j + 1],
4082 working_space[7 * j + 2],
4083 working_space[peak_vel],
4084 working_space[peak_vel + 1],
4085 working_space[peak_vel + 2],
4086 working_space[peak_vel + 6],
4087 working_space[peak_vel + 7],
4088 working_space[peak_vel + 12],
4089 working_space[peak_vel + 13]);
4093 working_space[2 * shift + k] =
4095 working_space[7 * j],
4096 working_space[7 * j + 1],
4097 working_space[7 * j + 2],
4098 working_space[peak_vel],
4099 working_space[peak_vel + 1],
4100 working_space[peak_vel + 2],
4101 working_space[peak_vel + 6],
4102 working_space[peak_vel + 7],
4103 working_space[peak_vel + 12],
4104 working_space[peak_vel + 13]);
4108 working_space[2 * shift + k] =
4110 working_space[7 * j],
4111 working_space[7 * j + 1],
4112 working_space[7 * j + 2],
4113 working_space[peak_vel],
4114 working_space[peak_vel + 1],
4115 working_space[peak_vel + 2],
4116 working_space[peak_vel + 6],
4117 working_space[peak_vel + 7],
4118 working_space[peak_vel + 12],
4119 working_space[peak_vel + 13]);
4123 working_space[2 * shift + k] =
4124 Derampx(i1, working_space[7 * j + 5],
4125 working_space[peak_vel],
4126 working_space[peak_vel + 8],
4127 working_space[peak_vel + 10],
4128 working_space[peak_vel + 12]);
4132 working_space[2 * shift + k] =
4133 Derampx(i2, working_space[7 * j + 6],
4134 working_space[peak_vel + 1],
4135 working_space[peak_vel + 9],
4136 working_space[peak_vel + 11],
4137 working_space[peak_vel + 13]);
4141 working_space[2 * shift + k] =
4142 Deri01(i1, working_space[7 * j + 3],
4143 working_space[7 * j + 5],
4144 working_space[peak_vel],
4145 working_space[peak_vel + 8],
4146 working_space[peak_vel + 10],
4147 working_space[peak_vel + 12]);
4151 working_space[2 * shift + k] =
4152 Deri01(i2, working_space[7 * j + 4],
4153 working_space[7 * j + 6],
4154 working_space[peak_vel + 1],
4155 working_space[peak_vel + 9],
4156 working_space[peak_vel + 11],
4157 working_space[peak_vel + 13]);
4161 working_space[2 * shift + k] =
4163 working_space, working_space[peak_vel],
4164 working_space[peak_vel + 1],
4165 working_space[peak_vel + 2],
4166 working_space[peak_vel + 6],
4167 working_space[peak_vel + 7],
4168 working_space[peak_vel + 8],
4169 working_space[peak_vel + 10],
4170 working_space[peak_vel + 12],
4171 working_space[peak_vel + 13]);
4175 working_space[2 * shift + k] =
4177 working_space, working_space[peak_vel],
4178 working_space[peak_vel + 1],
4179 working_space[peak_vel + 2],
4180 working_space[peak_vel + 6],
4181 working_space[peak_vel + 7],
4182 working_space[peak_vel + 9],
4183 working_space[peak_vel + 11],
4184 working_space[peak_vel + 12],
4185 working_space[peak_vel + 13]);
4189 working_space[2 * shift + k] =
4190 Derro(fNPeaks, i1, i2,
4191 working_space, working_space[peak_vel],
4192 working_space[peak_vel + 1],
4193 working_space[peak_vel + 2]);
4197 working_space[2 * shift + k] = 1.;
4201 working_space[2 * shift + k] = i1;
4205 working_space[2 * shift + k] = i2;
4209 working_space[2 * shift + k] =
4211 working_space, working_space[peak_vel],
4212 working_space[peak_vel + 1],
4213 working_space[peak_vel + 12],
4214 working_space[peak_vel + 13]);
4218 working_space[2 * shift + k] =
4220 working_space, working_space[peak_vel],
4221 working_space[peak_vel + 1]);
4225 working_space[2 * shift + k] =
4226 Dertx(fNPeaks, i1, working_space,
4227 working_space[peak_vel],
4228 working_space[peak_vel + 12]);
4232 working_space[2 * shift + k] =
4233 Derty(fNPeaks, i2, working_space,
4234 working_space[peak_vel + 1],
4235 working_space[peak_vel + 13]);
4239 working_space[2 * shift + k] =
4240 Dersx(fNPeaks, i1, working_space,
4241 working_space[peak_vel]);
4245 working_space[2 * shift + k] =
4246 Dersy(fNPeaks, i2, working_space,
4247 working_space[peak_vel + 1]);
4251 working_space[2 * shift + k] =
4252 Derbx(fNPeaks, i1, i2,
4253 working_space, working_space[peak_vel],
4254 working_space[peak_vel + 1],
4255 working_space[peak_vel + 6],
4256 working_space[peak_vel + 8],
4257 working_space[peak_vel + 12],
4258 working_space[peak_vel + 13]);
4262 working_space[2 * shift + k] =
4263 Derby(fNPeaks, i1, i2,
4264 working_space, working_space[peak_vel],
4265 working_space[peak_vel + 1],
4266 working_space[peak_vel + 6],
4267 working_space[peak_vel + 8],
4268 working_space[peak_vel + 12],
4269 working_space[peak_vel + 13]);
4272 yw = source[i1][i2];
4274 f =
Shape2(fNPeaks, i1, i2,
4275 working_space, working_space[peak_vel],
4276 working_space[peak_vel + 1],
4277 working_space[peak_vel + 2],
4278 working_space[peak_vel + 3],
4279 working_space[peak_vel + 4],
4280 working_space[peak_vel + 5],
4281 working_space[peak_vel + 6],
4282 working_space[peak_vel + 7],
4283 working_space[peak_vel + 8],
4284 working_space[peak_vel + 9],
4285 working_space[peak_vel + 10],
4286 working_space[peak_vel + 11],
4287 working_space[peak_vel + 12],
4288 working_space[peak_vel + 13]);
4296 chi_opt += (yw -
f) * (yw - f) / ywm;
4314 for (j = 0; j < size; j++) {
4315 for (k = 0; k < size; k++) {
4316 b = working_space[2 * shift +
4317 j] * working_space[2 * shift +
4320 b = b * (4 * yw - 2 *
f) / ywm;
4321 working_matrix[j][k] +=
b;
4323 working_space[3 * shift + j] +=
b;
4327 b = (f * f - yw * yw) / (ywm * ywm);
4331 for (j = 0; j < size; j++) {
4332 working_matrix[j][size] -=
4333 b * working_space[2 * shift + j];
4337 for (i = 0; i < size; i++) {
4338 working_matrix[i][size + 1] = 0;
4341 for (i = 0; i < size; i++) {
4342 working_space[2 * shift + i] = working_matrix[i][size + 1];
4351 for (j = 0; j < size; j++) {
4352 working_space[4 * shift + j] = working_space[shift + j];
4358 chi_min = 10000 * chi2;
4361 chi_min = 0.1 * chi2;
4363 for (
pi = 0.1; flag == 0 &&
pi <= 100;
pi += 0.1) {
4364 for (j = 0; j < size; j++) {
4365 working_space[shift + j] = working_space[4 * shift + j] +
pi * alpha * working_space[2 * shift + j];
4367 for (i = 0, j = 0; i <
fNPeaks; i++) {
4369 if (working_space[shift + j] < 0)
4370 working_space[shift + j] = 0;
4371 working_space[7 * i] = working_space[shift + j];
4375 if (working_space[shift + j] <
fXmin)
4376 working_space[shift + j] =
fXmin;
4377 if (working_space[shift + j] >
fXmax)
4378 working_space[shift + j] =
fXmax;
4379 working_space[7 * i + 1] = working_space[shift + j];
4383 if (working_space[shift + j] <
fYmin)
4384 working_space[shift + j] =
fYmin;
4385 if (working_space[shift + j] >
fYmax)
4386 working_space[shift + j] =
fYmax;
4387 working_space[7 * i + 2] = working_space[shift + j];
4391 if (working_space[shift + j] < 0)
4392 working_space[shift + j] = 0;
4393 working_space[7 * i + 3] = working_space[shift + j];
4397 if (working_space[shift + j] < 0)
4398 working_space[shift + j] = 0;
4399 working_space[7 * i + 4] = working_space[shift + j];
4403 if (working_space[shift + j] <
fXmin)
4404 working_space[shift + j] =
fXmin;
4405 if (working_space[shift + j] >
fXmax)
4406 working_space[shift + j] =
fXmax;
4407 working_space[7 * i + 5] = working_space[shift + j];
4411 if (working_space[shift + j] <
fYmin)
4412 working_space[shift + j] =
fYmin;
4413 if (working_space[shift + j] >
fYmax)
4414 working_space[shift + j] =
fYmax;
4415 working_space[7 * i + 6] = working_space[shift + j];
4420 if (working_space[shift + j] < 0.001) {
4421 working_space[shift + j] = 0.001;
4423 working_space[peak_vel] = working_space[shift + j];
4427 if (working_space[shift + j] < 0.001) {
4428 working_space[shift + j] = 0.001;
4430 working_space[peak_vel + 1] = working_space[shift + j];
4434 if (working_space[shift + j] < -1) {
4435 working_space[shift + j] = -1;
4437 if (working_space[shift + j] > 1) {
4438 working_space[shift + j] = 1;
4440 working_space[peak_vel + 2] = working_space[shift + j];
4444 working_space[peak_vel + 3] = working_space[shift + j];
4448 working_space[peak_vel + 4] = working_space[shift + j];
4452 working_space[peak_vel + 5] = working_space[shift + j];
4456 working_space[peak_vel + 6] = working_space[shift + j];
4460 working_space[peak_vel + 7] = working_space[shift + j];
4464 working_space[peak_vel + 8] = working_space[shift + j];
4468 working_space[peak_vel + 9] = working_space[shift + j];
4472 working_space[peak_vel + 10] = working_space[shift + j];
4476 working_space[peak_vel + 11] = working_space[shift + j];
4480 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
4481 if (working_space[shift + j] < 0)
4482 working_space[shift + j] = -0.001;
4484 working_space[shift + j] = 0.001;
4486 working_space[peak_vel + 12] = working_space[shift + j];
4490 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
4491 if (working_space[shift + j] < 0)
4492 working_space[shift + j] = -0.001;
4494 working_space[shift + j] = 0.001;
4496 working_space[peak_vel + 13] = working_space[shift + j];
4502 yw = source[i1][i2];
4506 working_space[peak_vel],
4507 working_space[peak_vel + 1],
4508 working_space[peak_vel + 2],
4509 working_space[peak_vel + 3],
4510 working_space[peak_vel + 4],
4511 working_space[peak_vel + 5],
4512 working_space[peak_vel + 6],
4513 working_space[peak_vel + 7],
4514 working_space[peak_vel + 8],
4515 working_space[peak_vel + 9],
4516 working_space[peak_vel + 10],
4517 working_space[peak_vel + 11],
4518 working_space[peak_vel + 12],
4519 working_space[peak_vel + 13]);
4532 chi2 += (yw -
f) * (yw - f) / ywm;
4540 pmin =
pi, chi_min = chi2;
4550 for (j = 0; j < size; j++) {
4551 working_space[shift + j] = working_space[4 * shift + j] + pmin * alpha * working_space[2 * shift + j];
4553 for (i = 0, j = 0; i <
fNPeaks; i++) {
4555 if (working_space[shift + j] < 0)
4556 working_space[shift + j] = 0;
4557 working_space[7 * i] = working_space[shift + j];
4561 if (working_space[shift + j] <
fXmin)
4562 working_space[shift + j] =
fXmin;
4563 if (working_space[shift + j] >
fXmax)
4564 working_space[shift + j] =
fXmax;
4565 working_space[7 * i + 1] = working_space[shift + j];
4569 if (working_space[shift + j] <
fYmin)
4570 working_space[shift + j] =
fYmin;
4571 if (working_space[shift + j] >
fYmax)
4572 working_space[shift + j] =
fYmax;
4573 working_space[7 * i + 2] = working_space[shift + j];
4577 if (working_space[shift + j] < 0)
4578 working_space[shift + j] = 0;
4579 working_space[7 * i + 3] = working_space[shift + j];
4583 if (working_space[shift + j] < 0)
4584 working_space[shift + j] = 0;
4585 working_space[7 * i + 4] = working_space[shift + j];
4589 if (working_space[shift + j] <
fXmin)
4590 working_space[shift + j] =
fXmin;
4591 if (working_space[shift + j] >
fXmax)
4592 working_space[shift + j] =
fXmax;
4593 working_space[7 * i + 5] = working_space[shift + j];
4597 if (working_space[shift + j] <
fYmin)
4598 working_space[shift + j] =
fYmin;
4599 if (working_space[shift + j] >
fYmax)
4600 working_space[shift + j] =
fYmax;
4601 working_space[7 * i + 6] = working_space[shift + j];
4606 if (working_space[shift + j] < 0.001) {
4607 working_space[shift + j] = 0.001;
4609 working_space[peak_vel] = working_space[shift + j];
4613 if (working_space[shift + j] < 0.001) {
4614 working_space[shift + j] = 0.001;
4616 working_space[peak_vel + 1] = working_space[shift + j];
4620 if (working_space[shift + j] < -1) {
4621 working_space[shift + j] = -1;
4623 if (working_space[shift + j] > 1) {
4624 working_space[shift + j] = 1;
4626 working_space[peak_vel + 2] = working_space[shift + j];
4630 working_space[peak_vel + 3] = working_space[shift + j];
4634 working_space[peak_vel + 4] = working_space[shift + j];
4638 working_space[peak_vel + 5] = working_space[shift + j];
4642 working_space[peak_vel + 6] = working_space[shift + j];
4646 working_space[peak_vel + 7] = working_space[shift + j];
4650 working_space[peak_vel + 8] = working_space[shift + j];
4654 working_space[peak_vel + 9] = working_space[shift + j];
4658 working_space[peak_vel + 10] = working_space[shift + j];
4662 working_space[peak_vel + 11] = working_space[shift + j];
4666 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
4667 if (working_space[shift + j] < 0)
4668 working_space[shift + j] = -0.001;
4670 working_space[shift + j] = 0.001;
4672 working_space[peak_vel + 12] = working_space[shift + j];
4676 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
4677 if (working_space[shift + j] < 0)
4678 working_space[shift + j] = -0.001;
4680 working_space[shift + j] = 0.001;
4682 working_space[peak_vel + 13] = working_space[shift + j];
4690 for (j = 0; j < size; j++) {
4691 working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j];
4693 for (i = 0, j = 0; i <
fNPeaks; i++) {
4695 if (working_space[shift + j] < 0)
4696 working_space[shift + j] = 0;
4697 working_space[7 * i] = working_space[shift + j];
4701 if (working_space[shift + j] <
fXmin)
4702 working_space[shift + j] =
fXmin;
4703 if (working_space[shift + j] >
fXmax)
4704 working_space[shift + j] =
fXmax;
4705 working_space[7 * i + 1] = working_space[shift + j];
4709 if (working_space[shift + j] <
fYmin)
4710 working_space[shift + j] =
fYmin;
4711 if (working_space[shift + j] >
fYmax)
4712 working_space[shift + j] =
fYmax;
4713 working_space[7 * i + 2] = working_space[shift + j];
4717 if (working_space[shift + j] < 0)
4718 working_space[shift + j] = 0;
4719 working_space[7 * i + 3] = working_space[shift + j];
4723 if (working_space[shift + j] < 0)
4724 working_space[shift + j] = 0;
4725 working_space[7 * i + 4] = working_space[shift + j];
4729 if (working_space[shift + j] <
fXmin)
4730 working_space[shift + j] =
fXmin;
4731 if (working_space[shift + j] >
fXmax)
4732 working_space[shift + j] =
fXmax;
4733 working_space[7 * i + 5] = working_space[shift + j];
4737 if (working_space[shift + j] <
fYmin)
4738 working_space[shift + j] =
fYmin;
4739 if (working_space[shift + j] >
fYmax)
4740 working_space[shift + j] =
fYmax;
4741 working_space[7 * i + 6] = working_space[shift + j];
4746 if (working_space[shift + j] < 0.001) {
4747 working_space[shift + j] = 0.001;
4749 working_space[peak_vel] = working_space[shift + j];
4753 if (working_space[shift + j] < 0.001) {
4754 working_space[shift + j] = 0.001;
4756 working_space[peak_vel + 1] = working_space[shift + j];
4760 if (working_space[shift + j] < -1) {
4761 working_space[shift + j] = -1;
4763 if (working_space[shift + j] > 1) {
4764 working_space[shift + j] = 1;
4766 working_space[peak_vel + 2] = working_space[shift + j];
4770 working_space[peak_vel + 3] = working_space[shift + j];
4774 working_space[peak_vel + 4] = working_space[shift + j];
4778 working_space[peak_vel + 5] = working_space[shift + j];
4782 working_space[peak_vel + 6] = working_space[shift + j];
4786 working_space[peak_vel + 7] = working_space[shift + j];
4790 working_space[peak_vel + 8] = working_space[shift + j];
4794 working_space[peak_vel + 9] = working_space[shift + j];
4798 working_space[peak_vel + 10] = working_space[shift + j];
4802 working_space[peak_vel + 11] = working_space[shift + j];
4806 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
4807 if (working_space[shift + j] < 0)
4808 working_space[shift + j] = -0.001;
4810 working_space[shift + j] = 0.001;
4812 working_space[peak_vel + 12] = working_space[shift + j];
4816 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
4817 if (working_space[shift + j] < 0)
4818 working_space[shift + j] = -0.001;
4820 working_space[shift + j] = 0.001;
4822 working_space[peak_vel + 13] = working_space[shift + j];
4828 yw = source[i1][i2];
4830 f =
Shape2(fNPeaks, i1, i2,
4831 working_space, working_space[peak_vel],
4832 working_space[peak_vel + 1],
4833 working_space[peak_vel + 2],
4834 working_space[peak_vel + 3],
4835 working_space[peak_vel + 4],
4836 working_space[peak_vel + 5],
4837 working_space[peak_vel + 6],
4838 working_space[peak_vel + 7],
4839 working_space[peak_vel + 8],
4840 working_space[peak_vel + 9],
4841 working_space[peak_vel + 10],
4842 working_space[peak_vel + 11],
4843 working_space[peak_vel + 12],
4844 working_space[peak_vel + 13]);
4857 chi += (yw -
f) * (yw - f) / ywm;
4865 alpha = alpha * chi_opt / (2 * chi);
4868 alpha = alpha / 10.0;
4871 }
while (((chi > chi_opt
4876 for (j = 0; j < size; j++) {
4877 working_space[4 * shift + j] = 0;
4878 working_space[2 * shift + j] = 0;
4880 for (i1 =
fXmin, chi_cel = 0; i1 <=
fXmax; i1++) {
4882 yw = source[i1][i2];
4885 f =
Shape2(fNPeaks, i1, i2,
4886 working_space, working_space[peak_vel],
4887 working_space[peak_vel + 1],
4888 working_space[peak_vel + 2],
4889 working_space[peak_vel + 3],
4890 working_space[peak_vel + 4],
4891 working_space[peak_vel + 5],
4892 working_space[peak_vel + 6],
4893 working_space[peak_vel + 7],
4894 working_space[peak_vel + 8],
4895 working_space[peak_vel + 9],
4896 working_space[peak_vel + 10],
4897 working_space[peak_vel + 11],
4898 working_space[peak_vel + 12],
4899 working_space[peak_vel + 13]);
4900 chi_opt = (yw -
f) * (yw - f) / yw;
4901 chi_cel += (yw -
f) * (yw - f) / yw;
4904 for (j = 0, k = 0; j <
fNPeaks; j++) {
4907 working_space[7 * j + 1],
4908 working_space[7 * j + 2],
4909 working_space[peak_vel],
4910 working_space[peak_vel + 1],
4911 working_space[peak_vel + 2],
4912 working_space[peak_vel + 6],
4913 working_space[peak_vel + 7],
4914 working_space[peak_vel + 12],
4915 working_space[peak_vel + 13]);
4917 working_space[2 * shift + k] += chi_opt;
4919 working_space[4 * shift + k] +=
b;
4925 working_space[7 * j],
4926 working_space[7 * j + 1],
4927 working_space[7 * j + 2],
4928 working_space[peak_vel],
4929 working_space[peak_vel + 1],
4930 working_space[peak_vel + 2],
4931 working_space[peak_vel + 6],
4932 working_space[peak_vel + 7],
4933 working_space[peak_vel + 12],
4934 working_space[peak_vel + 13]);
4936 working_space[2 * shift + k] += chi_opt;
4938 working_space[4 * shift + k] +=
b;
4944 working_space[7 * j],
4945 working_space[7 * j + 1],
4946 working_space[7 * j + 2],
4947 working_space[peak_vel],
4948 working_space[peak_vel + 1],
4949 working_space[peak_vel + 2],
4950 working_space[peak_vel + 6],
4951 working_space[peak_vel + 7],
4952 working_space[peak_vel + 12],
4953 working_space[peak_vel + 13]);
4955 working_space[2 * shift + k] += chi_opt;
4957 working_space[4 * shift + k] +=
b;
4962 a =
Derampx(i1, working_space[7 * j + 5],
4963 working_space[peak_vel],
4964 working_space[peak_vel + 8],
4965 working_space[peak_vel + 10],
4966 working_space[peak_vel + 12]);
4968 working_space[2 * shift + k] += chi_opt;
4970 working_space[4 * shift + k] +=
b;
4975 a =
Derampx(i2, working_space[7 * j + 6],
4976 working_space[peak_vel + 1],
4977 working_space[peak_vel + 9],
4978 working_space[peak_vel + 11],
4979 working_space[peak_vel + 13]);
4981 working_space[2 * shift + k] += chi_opt;
4983 working_space[4 * shift + k] +=
b;
4988 a =
Deri01(i1, working_space[7 * j + 3],
4989 working_space[7 * j + 5],
4990 working_space[peak_vel],
4991 working_space[peak_vel + 8],
4992 working_space[peak_vel + 10],
4993 working_space[peak_vel + 12]);
4995 working_space[2 * shift + k] += chi_opt;
4997 working_space[4 * shift + k] +=
b;
5002 a =
Deri01(i2, working_space[7 * j + 4],
5003 working_space[7 * j + 6],
5004 working_space[peak_vel + 1],
5005 working_space[peak_vel + 9],
5006 working_space[peak_vel + 11],
5007 working_space[peak_vel + 13]);
5009 working_space[2 * shift + k] += chi_opt;
5011 working_space[4 * shift + k] +=
b;
5018 working_space, working_space[peak_vel],
5019 working_space[peak_vel + 1],
5020 working_space[peak_vel + 2],
5021 working_space[peak_vel + 6],
5022 working_space[peak_vel + 7],
5023 working_space[peak_vel + 8],
5024 working_space[peak_vel + 10],
5025 working_space[peak_vel + 12],
5026 working_space[peak_vel + 13]);
5028 working_space[2 * shift + k] += chi_opt;
5030 working_space[4 * shift + k] +=
b;
5036 working_space, working_space[peak_vel],
5037 working_space[peak_vel + 1],
5038 working_space[peak_vel + 2],
5039 working_space[peak_vel + 6],
5040 working_space[peak_vel + 7],
5041 working_space[peak_vel + 9],
5042 working_space[peak_vel + 11],
5043 working_space[peak_vel + 12],
5044 working_space[peak_vel + 13]);
5046 working_space[2 * shift + k] += chi_opt;
5048 working_space[4 * shift + k] +=
b;
5053 a =
Derro(fNPeaks, i1, i2,
5054 working_space, working_space[peak_vel],
5055 working_space[peak_vel + 1],
5056 working_space[peak_vel + 2]);
5058 working_space[2 * shift + k] += chi_opt;
5060 working_space[4 * shift + k] +=
b;
5067 working_space[2 * shift + k] += chi_opt;
5069 working_space[4 * shift + k] +=
b;
5076 working_space[2 * shift + k] += chi_opt;
5078 working_space[4 * shift + k] +=
b;
5085 working_space[2 * shift + k] += chi_opt;
5087 working_space[4 * shift + k] +=
b;
5092 a =
Dertxy(fNPeaks, i1, i2,
5093 working_space, working_space[peak_vel],
5094 working_space[peak_vel + 1],
5095 working_space[peak_vel + 12],
5096 working_space[peak_vel + 13]);
5098 working_space[2 * shift + k] += chi_opt;
5100 working_space[4 * shift + k] +=
b;
5105 a =
Dersxy(fNPeaks, i1, i2,
5106 working_space, working_space[peak_vel],
5107 working_space[peak_vel + 1]);
5109 working_space[2 * shift + k] += chi_opt;
5111 working_space[4 * shift + k] +=
b;
5116 a =
Dertx(fNPeaks, i1, working_space,
5117 working_space[peak_vel],
5118 working_space[peak_vel + 12]);
5120 working_space[2 * shift + k] += chi_opt;
5122 working_space[4 * shift + k] +=
b;
5127 a =
Derty(fNPeaks, i2, working_space,
5128 working_space[peak_vel + 1],
5129 working_space[peak_vel + 13]);
5131 working_space[2 * shift + k] += chi_opt;
5133 working_space[4 * shift + k] +=
b;
5138 a =
Dersx(fNPeaks, i1, working_space,
5139 working_space[peak_vel]);
5141 working_space[2 * shift + k] += chi_opt;
5143 working_space[4 * shift + k] +=
b;
5148 a =
Dersy(fNPeaks, i2, working_space,
5149 working_space[peak_vel + 1]);
5151 working_space[2 * shift + k] += chi_opt;
5153 working_space[4 * shift + k] +=
b;
5158 a =
Derbx(fNPeaks, i1, i2,
5159 working_space, working_space[peak_vel],
5160 working_space[peak_vel + 1],
5161 working_space[peak_vel + 6],
5162 working_space[peak_vel + 8],
5163 working_space[peak_vel + 12],
5164 working_space[peak_vel + 13]);
5166 working_space[2 * shift + k] += chi_opt;
5168 working_space[4 * shift + k] +=
b;
5173 a =
Derby(fNPeaks, i1, i2,
5174 working_space, working_space[peak_vel],
5175 working_space[peak_vel + 1],
5176 working_space[peak_vel + 6],
5177 working_space[peak_vel + 8],
5178 working_space[peak_vel + 12],
5179 working_space[peak_vel + 13]);
5181 working_space[2 * shift + k] += chi_opt;
5183 working_space[4 * shift + k] +=
b;
5191 chi_er = chi_cel /
b;
5192 for (i = 0, j = 0; i <
fNPeaks; i++) {
5194 Volume(working_space[7 * i], working_space[peak_vel],
5195 working_space[peak_vel + 1], working_space[peak_vel + 2]);
5199 a =
Derpa2(working_space[peak_vel],
5200 working_space[peak_vel + 1],
5201 working_space[peak_vel + 2]);
5202 b = working_space[4 * shift + j];
5212 working_space[peak_vel + 1],
5213 working_space[peak_vel + 2]);
5214 b = working_space[4 * shift + peak_vel];
5224 working_space[peak_vel],
5225 working_space[peak_vel + 2]);
5226 b = working_space[4 * shift + peak_vel + 1];
5235 a =
Derpro(working_space[shift + j], working_space[peak_vel],
5236 working_space[peak_vel + 1],
5237 working_space[peak_vel + 2]);
5238 b = working_space[4 * shift + peak_vel + 2];
5253 fAmpCalc[i] = working_space[shift + j];
5254 if (working_space[3 * shift + j] != 0)
5265 if (working_space[3 * shift + j] != 0)
5276 if (working_space[3 * shift + j] != 0)
5287 if (working_space[3 * shift + j] != 0)
5298 if (working_space[3 * shift + j] != 0)
5309 if (working_space[3 * shift + j] != 0)
5320 if (working_space[3 * shift + j] != 0)
5332 if (working_space[3 * shift + j] != 0)
5343 if (working_space[3 * shift + j] != 0)
5353 fRoCalc = working_space[shift + j];
5354 if (working_space[3 * shift + j] != 0)
5364 fA0Calc = working_space[shift + j];
5365 if (working_space[3 * shift + j] != 0)
5375 fAxCalc = working_space[shift + j];
5376 if (working_space[3 * shift + j] != 0)
5386 fAyCalc = working_space[shift + j];
5387 if (working_space[3 * shift + j] != 0)
5397 fTxyCalc = working_space[shift + j];
5398 if (working_space[3 * shift + j] != 0)
5408 fSxyCalc = working_space[shift + j];
5409 if (working_space[3 * shift + j] != 0)
5419 fTxCalc = working_space[shift + j];
5420 if (working_space[3 * shift + j] != 0)
5430 fTyCalc = working_space[shift + j];
5431 if (working_space[3 * shift + j] != 0)
5441 fSxCalc = working_space[shift + j];
5442 if (working_space[3 * shift + j] != 0)
5452 fSyCalc = working_space[shift + j];
5453 if (working_space[3 * shift + j] != 0)
5463 fBxCalc = working_space[shift + j];
5464 if (working_space[3 * shift + j] != 0)
5474 fByCalc = working_space[shift + j];
5475 if (working_space[3 * shift + j] != 0)
5488 f =
Shape2(fNPeaks, i1, i2,
5489 working_space, working_space[peak_vel],
5490 working_space[peak_vel + 1],
5491 working_space[peak_vel + 2],
5492 working_space[peak_vel + 3],
5493 working_space[peak_vel + 4],
5494 working_space[peak_vel + 5],
5495 working_space[peak_vel + 6],
5496 working_space[peak_vel + 7],
5497 working_space[peak_vel + 8],
5498 working_space[peak_vel + 9],
5499 working_space[peak_vel + 10],
5500 working_space[peak_vel + 11],
5501 working_space[peak_vel + 12],
5502 working_space[peak_vel + 13]);
5507 for (i = 0; i < size; i++)
delete [] working_matrix[i];
5508 delete [] working_matrix;
5509 delete [] working_space;
5525 if(xmin<0 || xmax <= xmin || ymin<0 || ymax <= ymin){
5526 Error(
"SetFitParameters",
"Wrong range");
5529 if (numberIterations <= 0){
5530 Error(
"SetFitParameters",
"Invalid number of iterations, must be positive");
5533 if (alpha <= 0 || alpha > 1){
5534 Error (
"SetFitParameters",
"Invalid step coefficient alpha, must be > than 0 and <=1");
5540 Error(
"SetFitParameters",
"Wrong type of statistic");
5545 Error(
"SetFitParameters",
"Wrong optimization algorithm");
5551 Error(
"SetFitParameters",
"Wrong power");
5556 Error(
"SetFitParameters",
"Wrong order of Taylor development");
5581 void TSpectrum2Fit::SetPeakParameters(
Double_t sigmaX,
Bool_t fixSigmaX,
Double_t sigmaY,
Bool_t fixSigmaY,
Double_t ro,
Bool_t fixRo,
const Double_t *positionInitX,
const Bool_t *fixPositionX,
const Double_t *positionInitY,
const Bool_t *fixPositionY,
const Double_t *positionInitX1,
const Bool_t *fixPositionX1,
const Double_t *positionInitY1,
const Bool_t *fixPositionY1,
const Double_t *ampInit,
const Bool_t *fixAmp,
const Double_t *ampInitX1,
const Bool_t *fixAmpX1,
const Double_t *ampInitY1,
const Bool_t *fixAmpY1)
5583 if (sigmaX <= 0 || sigmaY <= 0){
5584 Error (
"SetPeakParameters",
"Invalid sigma, must be > than 0");
5587 if (ro < -1 || ro > 1){
5588 Error (
"SetPeakParameters",
"Invalid ro, must be from region <-1,1>");
5593 if(positionInitX[i] <
fXmin || positionInitX[i] >
fXmax){
5594 Error (
"SetPeakParameters",
"Invalid peak position, must be in the range fXmin, fXmax");
5597 if(positionInitY[i] <
fYmin || positionInitY[i] >
fYmax){
5598 Error (
"SetPeakParameters",
"Invalid peak position, must be in the range fYmin, fYmax");
5601 if(positionInitX1[i] <
fXmin || positionInitX1[i] >
fXmax){
5602 Error (
"SetPeakParameters",
"Invalid ridge position, must be in the range fXmin, fXmax");
5605 if(positionInitY1[i] <
fYmin || positionInitY1[i] >
fYmax){
5606 Error (
"SetPeakParameters",
"Invalid ridge position, must be in the range fYmin, fYmax");
5610 Error (
"SetPeakParameters",
"Invalid peak amplitude, must be > than 0");
5613 if(ampInitX1[i] < 0){
5614 Error (
"SetPeakParameters",
"Invalid x ridge amplitude, must be > than 0");
5617 if(ampInitY1[i] < 0){
5618 Error (
"SetPeakParameters",
"Invalid y ridge amplitude, must be > than 0");
5679 void TSpectrum2Fit::SetTailParameters(
Double_t tInitXY,
Bool_t fixTxy,
Double_t tInitX,
Bool_t fixTx,
Double_t tInitY,
Bool_t fixTy,
Double_t bInitX,
Bool_t fixBx,
Double_t bInitY,
Bool_t fixBy,
Double_t sInitXY,
Bool_t fixSxy,
Double_t sInitX,
Bool_t fixSx,
Double_t sInitY,
Bool_t fixSy)
5757 amplitudeErrors[i] =
fAmpErr[i];
5856 void TSpectrum2Fit::GetTailParameters(
Double_t &txy,
Double_t &txyErr,
Double_t &tx,
Double_t &txErr,
Double_t &ty,
Double_t &tyErr,
Double_t &bx,
Double_t &bxErr,
Double_t &by,
Double_t &byErr,
Double_t &sxy,
Double_t &sxyErr,
Double_t &sx,
Double_t &sxErr,
Double_t &sy,
Double_t &syErr)
Bool_t * fFixAmp
[fNPeaks] array of logical values which allow to fix appropriate amplitudes of 2D peaks (not fit)...
Double_t Derpsigmay(Double_t a, Double_t sx, Double_t ro)
This function calculates derivative of the volume of a peak according to sigmay.
Double_t Derpa2(Double_t sx, Double_t sy, Double_t ro)
This function calculates derivative of the volume of a peak according to amplitude.
Double_t fByCalc
calculated value of b parameter for 1D ridges in y direction
Double_t Derampx(Double_t x, Double_t x0, Double_t sigmax, Double_t tx, Double_t sx, Double_t bx)
This function calculates derivative of 2D peaks shape function (see manual) according to amplitude of...
Bool_t fFixSx
logical value of s parameter for 1D ridges in x direction, which allows to fix the parameter (not to ...
Double_t Dersigmay(Int_t numOfFittedPeaks, Double_t x, Double_t y, const Double_t *parameter, Double_t sigmax, Double_t sigmay, Double_t ro, Double_t txy, Double_t sxy, Double_t ty, Double_t sy, Double_t bx, Double_t by)
This function calculates derivative of peaks shape function (see manual) according to sigmax of peaks...
Double_t * fAmpCalcY1
[fNPeaks] array of calculated values of amplitudes of 1D ridges in y direction, output parameters ...
Int_t fYmax
last fitted channel in y direction
Double_t Dertx(Int_t numOfFittedPeaks, Double_t x, const Double_t *parameter, Double_t sigmax, Double_t bx)
This function calculates derivative of peaks shape function (see manual) according to relative amplit...
Double_t * fPositionErrY1
[fNPeaks] array of y positions errors of 1D ridges, output parameters
Double_t Deri01(Double_t x, Double_t ax, Double_t x0, Double_t sigmax, Double_t tx, Double_t sx, Double_t bx)
This function calculates derivative of 2D peaks shape function (see manual) according to x position o...
Advanced 2-dimensional spectra fitting functions.
Double_t * fPositionCalcX1
[fNPeaks] array of calculated x positions of 1D ridges, output parameters
Bool_t * fFixPositionY1
[fNPeaks] array of logical values which allow to fix appropriate y positions of 1D ridges (not fit)...
Bool_t * fFixAmpY1
[fNPeaks] array of logical values which allow to fix appropriate amplitudes of 1D ridges in y directi...
Double_t fAxErr
error value of background ax parameter
Double_t fSxInit
initial value of s parameter for 1D ridges in x direction (relative amplitude of step), for details see html manual and references
Bool_t fFixSxy
logical value of s parameter for 2D peaks, which allows to fix the parameter (not to fit)...
Double_t fTxInit
initial value of t parameter for 1D ridges in x direction (relative amplitude of tail), for details see html manual and references
Double_t * fPositionInitY
[fNPeaks] array of initial values of y positions of 2D peaks, input parameters
Int_t fXmin
first fitted channel in x direction
void GetBackgroundParameters(Double_t &a0, Double_t &a0Err, Double_t &ax, Double_t &axErr, Double_t &ay, Double_t &ayErr)
This function gets the background parameters and their errors.
Double_t Derpsigmax(Double_t a, Double_t sy, Double_t ro)
This function calculates derivative of the volume of a peak according to sigmax.
Int_t fYmin
first fitted channel in y direction
Bool_t * fFixAmpX1
[fNPeaks] array of logical values which allow to fix appropriate amplitudes of 1D ridges in x directi...
Double_t fSyErr
error value of s parameter for 1D ridges in y direction
void GetTailParameters(Double_t &txy, Double_t &txyErr, Double_t &tx, Double_t &txErr, Double_t &ty, Double_t &tyErr, Double_t &bx, Double_t &bxErr, Double_t &by, Double_t &byErr, Double_t &sxy, Double_t &sxyErr, Double_t &sx, Double_t &sxErr, Double_t &sy, Double_t &syErr)
This function gets the tail parameters and their errors.
Bool_t * fFixPositionX
[fNPeaks] array of logical values which allow to fix appropriate x positions of 2D peaks (not fit)...
Double_t fByInit
initial value of b parameter for 1D ridges in y direction (slope), for details see html manual and re...
Bool_t fFixA0
logical value of a0 parameter, which allows to fix the parameter (not to fit).
Double_t Derfc(Double_t x)
This function calculates derivative of error function of x.
Double_t fA0Init
initial value of background a0 parameter(backgroud is estimated as a0+ax*x+ay*y)
Bool_t fFixSy
logical value of s parameter for 1D ridges in y direction, which allows to fix the parameter (not to ...
Double_t Derpro(Double_t a, Double_t sx, Double_t sy, Double_t ro)
This function calculates derivative of the volume of a peak according to ro.
Double_t Deri02(Double_t x, Double_t y, Double_t a, Double_t x0, Double_t y0, Double_t sigmax, Double_t sigmay, Double_t ro, Double_t txy, Double_t sxy, Double_t bx, Double_t by)
This function calculates derivative of 2D peaks shape function (see manual) according to x position o...
Double_t fSigmaCalcY
calculated value of sigma y parameter
Double_t * fPositionCalcY1
[fNPeaks] array of calculated y positions of 1D ridges, output parameters
Double_t * fPositionInitY1
[fNPeaks] array of initial y positions of 1D ridges, input parameters
void GetAmplitudes(Double_t *amplitudes, Double_t *amplitudesX1, Double_t *amplitudesY1)
This function gets the amplitudes of fitted 2D peaks and 1D ridges.
Double_t * fPositionErrY
[fNPeaks] array of error values of y positions of 2D peaks, output parameters
Double_t Derty(Int_t numOfFittedPeaks, Double_t x, const Double_t *parameter, Double_t sigmax, Double_t bx)
This function calculates derivative of peaks shape function (see manual) according to relative amplit...
TSpectrum2Fit(void)
Default constructor.
Double_t fTxyErr
error value of t parameter for 2D peaks
Double_t * fPositionInitX
[fNPeaks] array of initial values of x positions of 2D peaks, input parameters
void FitStiefel(Double_t **source)
This function fits the source spectrum.
void GetPositions(Double_t *positionsX, Double_t *positionsY, Double_t *positionsX1, Double_t *positionsY1)
This function gets the positions of fitted 2D peaks and 1D ridges.
Double_t Derdersigmay(Int_t numOfFittedPeaks, Double_t x, Double_t y, const Double_t *parameter, Double_t sigmax, Double_t sigmay, Double_t ro)
This function calculates second derivative of peaks shape function (see manual) according to sigmay o...
Double_t * fAmpCalc
[fNPeaks] array of calculated values of amplitudes of 2D peaks, output parameters ...
Double_t fTxErr
error value of t parameter for 1D ridges in x direction
Double_t * fAmpCalcX1
[fNPeaks] array of calculated values of amplitudes of 1D ridges in x direction, output parameters ...
Double_t Derby(Int_t numOfFittedPeaks, Double_t x, Double_t y, const Double_t *parameter, Double_t sigmax, Double_t sigmay, Double_t txy, Double_t ty, Double_t bx, Double_t by)
This function calculates derivative of peaks shape function (see manual) according to slope by...
Double_t fSigmaInitY
initial value of sigma y parameter
Double_t Derderi01(Double_t x, Double_t ax, Double_t x0, Double_t sigmax)
This function calculates second derivative of 2D peaks shape function (see manual) according to x pos...
Double_t fAyInit
initial value of background ay parameter(backgroud is estimated as a0+ax*x+ay*y)
Double_t fRoErr
error value of correlation coefficient
Double_t Dersx(Int_t numOfFittedPeaks, Double_t x, const Double_t *parameter, Double_t sigmax)
This function calculates derivative of peaks shape function (see manual) according to relative amplit...
Bool_t fFixTy
logical value of t parameter for 1D ridges in y direction, which allows to fix the parameter (not to ...
Double_t Shape2(Int_t numOfFittedPeaks, Double_t x, Double_t y, const Double_t *parameter, Double_t sigmax, Double_t sigmay, Double_t ro, Double_t a0, Double_t ax, Double_t ay, Double_t txy, Double_t sxy, Double_t tx, Double_t ty, Double_t sx, Double_t sy, Double_t bx, Double_t by)
This function calculates 2D peaks shape function (see manual)
Double_t fAxCalc
calculated value of background ax parameter
Double_t * fPositionErrX1
[fNPeaks] array of x positions errors of 1D ridges, output parameters
Double_t * fAmpInitX1
[fNPeaks] array of initial values of amplitudes of 1D ridges in x direction, input parameters ...
Double_t * fAmpErrY1
[fNPeaks] array of amplitudes errors of 1D ridges in y direction, output parameters ...
void FitAwmi(Double_t **source)
This function fits the source spectrum.
The TNamed class is the base class for all named ROOT classes.
Bool_t fFixTxy
logical value of t parameter for 2D peaks, which allows to fix the parameter (not to fit)...
Int_t fFitTaylor
order of Taylor expansion, possible values kFitTaylorOrderFirst, kFitTaylorOrderSecond. It applies only for Awmi fitting function.
Double_t * fAmpInit
[fNPeaks] array of initial values of amplitudes of 2D peaks, input parameters
Double_t Deramp2(Double_t x, Double_t y, Double_t x0, Double_t y0, Double_t sigmax, Double_t sigmay, Double_t ro, Double_t txy, Double_t sxy, Double_t bx, Double_t by)
This function calculates derivative of 2D peaks shape function (see manual) according to amplitude of...
Double_t Ourpowl(Double_t a, Int_t pw)
power function
Int_t fXmax
last fitted channel in x direction
Double_t fSxCalc
calculated value of s parameter for 1D ridges in x direction
Bool_t fFixRo
logical value of correlation coefficient, which allows to fix the parameter (not to fit)...
Double_t fSyInit
initial value of s parameter for 1D ridges in y direction (relative amplitude of step), for details see html manual and references
Double_t * fAmpErrX1
[fNPeaks] array of amplitudes errors of 1D ridges in x direction, output parameters ...
Double_t fBxErr
error value of b parameter for 1D ridges in x direction
void GetSigmaX(Double_t &sigmaX, Double_t &sigmaErrX)
This function gets the sigma x parameter and its error.
void SetBackgroundParameters(Double_t a0Init, Bool_t fixA0, Double_t axInit, Bool_t fixAx, Double_t ayInit, Bool_t fixAy)
This function sets the following fitting parameters of background:
Int_t fAlphaOptim
optimization of convergence algorithm, possible values kFitAlphaHalving, kFitAlphaOptimal ...
Bool_t fFixAy
logical value of ay parameter, which allows to fix the parameter (not to fit).
void GetVolumes(Double_t *volumes)
This function gets the volumes of fitted 2D peaks.
Double_t * fPositionCalcY
[fNPeaks] array of calculated values of y positions of 2D peaks, output parameters ...
void SetPeakParameters(Double_t sigmaX, Bool_t fixSigmaX, Double_t sigmaY, Bool_t fixSigmaY, Double_t ro, Bool_t fixRo, const Double_t *positionInitX, const Bool_t *fixPositionX, const Double_t *positionInitY, const Bool_t *fixPositionY, const Double_t *positionInitX1, const Bool_t *fixPositionX1, const Double_t *positionInitY1, const Bool_t *fixPositionY1, const Double_t *ampInit, const Bool_t *fixAmp, const Double_t *ampInitX1, const Bool_t *fixAmpX1, const Double_t *ampInitY1, const Bool_t *fixAmpY1)
This function sets the following fitting parameters of peaks:
Double_t fTxyInit
initial value of t parameter for 2D peaks (relative amplitude of tail), for details see html manual a...
Double_t fAyCalc
calculated value of background ay parameter
void GetPositionErrors(Double_t *positionErrorsX, Double_t *positionErrorsY, Double_t *positionErrorsX1, Double_t *positionErrorsY1)
This function gets the errors of positions of fitted 2D peaks and 1D ridges.
Double_t * fVolumeErr
[fNPeaks] array of volumes errors of 2D peaks, output parameters
Bool_t fFixSigmaX
logical value of sigma x parameter, which allows to fix the parameter (not to fit).
unsigned int r1[N_CITIES]
Bool_t fFixSigmaY
logical value of sigma y parameter, which allows to fix the parameter (not to fit).
Int_t fNumberIterations
number of iterations in fitting procedure, input parameter, it should be > 0
Bool_t * fFixPositionX1
[fNPeaks] array of logical values which allow to fix appropriate x positions of 1D ridges (not fit)...
Double_t Derj02(Double_t x, Double_t y, Double_t a, Double_t x0, Double_t y0, Double_t sigmax, Double_t sigmay, Double_t ro, Double_t txy, Double_t sxy, Double_t bx, Double_t by)
This function calculates derivative of 2D peaks shape function (see manual) according to y position o...
Double_t fTyCalc
calculated value of t parameter for 1D ridges in y direction
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
void GetAmplitudeErrors(Double_t *amplitudeErrors, Double_t *amplitudeErrorsX1, Double_t *amplitudeErrorsY1)
This function gets the amplitudes of fitted 2D peaks and 1D ridges.
Bool_t fFixBy
logical value of b parameter for 1D ridges in y direction, which allows to fix the parameter (not to ...
Int_t fPower
possible values kFitPower2,4,6,8,10,12, for details see references. It applies only for Awmi fitting ...
void SetTailParameters(Double_t tInitXY, Bool_t fixTxy, Double_t tInitX, Bool_t fixTx, Double_t tInitY, Bool_t fixTy, Double_t bInitX, Bool_t fixBx, Double_t bInitY, Bool_t fixBy, Double_t sInitXY, Bool_t fixSxy, Double_t sInitX, Bool_t fixSx, Double_t sInitY, Bool_t fixSy)
This function sets the following fitting parameters of tails of peaks.
void GetSigmaY(Double_t &sigmaY, Double_t &sigmaErrY)
This function gets the sigma y parameter and its error.
Double_t fSxyErr
error value of s parameter for 2D peaks
void GetRo(Double_t &ro, Double_t &roErr)
This function gets the ro parameter and its error.
Int_t fStatisticType
type of statistics, possible values kFitOptimChiCounts (chi square statistics with counts as weightin...
Double_t Derdersigmax(Int_t numOfFittedPeaks, Double_t x, Double_t y, const Double_t *parameter, Double_t sigmax, Double_t sigmay, Double_t ro)
This function calculates second derivative of peaks shape function (see manual) according to sigmax o...
Double_t fByErr
error value of b parameter for 1D ridges in y direction
Double_t Dertxy(Int_t numOfFittedPeaks, Double_t x, Double_t y, const Double_t *parameter, Double_t sigmax, Double_t sigmay, Double_t bx, Double_t by)
This function calculates derivative of peaks shape function (see manual) according to relative amplit...
Double_t fA0Calc
calculated value of background a0 parameter
Bool_t fFixAx
logical value of ax parameter, which allows to fix the parameter (not to fit).
Double_t fTxyCalc
calculated value of t parameter for 2D peaks
Double_t Derderi02(Double_t x, Double_t y, Double_t a, Double_t x0, Double_t y0, Double_t sigmax, Double_t sigmay, Double_t ro)
This function calculates second derivative of 2D peaks shape function (see manual) according to x pos...
void GetVolumeErrors(Double_t *volumeErrors)
This function gets errors of the volumes of fitted 2D peaks.
Double_t fSigmaCalcX
calculated value of sigma x parameter
Double_t fAyErr
error value of background ay parameter
Double_t fTxCalc
calculated value of t parameter for 1D ridges in x direction
Double_t fSxyCalc
calculated value of s parameter for 2D peaks
Double_t fA0Err
error value of background a0 parameter
Double_t * fAmpErr
[fNPeaks] array of amplitudes errors of 2D peaks, output parameters
Double_t * fVolume
[fNPeaks] array of calculated volumes of 2D peaks, output parameters
Bool_t fFixTx
logical value of t parameter for 1D ridges in x direction, which allows to fix the parameter (not to ...
Double_t fBxInit
initial value of b parameter for 1D ridges in x direction (slope), for details see html manual and re...
Double_t fSyCalc
calculated value of s parameter for 1D ridges in y direction
Double_t Derro(Int_t numOfFittedPeaks, Double_t x, Double_t y, const Double_t *parameter, Double_t sx, Double_t sy, Double_t r)
This function calculates derivative of peaks shape function (see manual) according to correlation coe...
Bool_t * fFixPositionY
[fNPeaks] array of logical values which allow to fix appropriate y positions of 2D peaks (not fit)...
Double_t Dersxy(Int_t numOfFittedPeaks, Double_t x, Double_t y, const Double_t *parameter, Double_t sigmax, Double_t sigmay)
This function calculates derivative of peaks shape function (see manual) according to relative amplit...
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Double_t fRoCalc
calculated value of correlation coefficient
Double_t Erfc(Double_t x)
This function calculates error function of x.
Double_t Derderj02(Double_t x, Double_t y, Double_t a, Double_t x0, Double_t y0, Double_t sigmax, Double_t sigmay, Double_t ro)
This function calculates second derivative of 2D peaks shape function (see manual) according to y pos...
Double_t * fPositionErrX
[fNPeaks] array of error values of x positions of 2D peaks, output parameters
Double_t fSigmaErrX
error value of sigma x parameter
void StiefelInversion(Double_t **a, Int_t size)
This function calculates solution of the system of linear equations.
Double_t * fPositionCalcX
[fNPeaks] array of calculated values of x positions of 2D peaks, output parameters ...
Int_t fNPeaks
number of peaks present in fit, input parameter, it should be > 0
void SetFitParameters(Int_t xmin, Int_t xmax, Int_t ymin, Int_t ymax, Int_t numberIterations, Double_t alpha, Int_t statisticType, Int_t alphaOptim, Int_t power, Int_t fitTaylor)
This function sets the following fitting parameters:
Double_t fBxCalc
calculated value of b parameter for 1D ridges in x direction
Double_t fSigmaInitX
initial value of sigma x parameter
Double_t fSxyInit
initial value of s parameter for 2D peaks (relative amplitude of step), for details see html manual a...
Double_t fAxInit
initial value of background ax parameter(backgroud is estimated as a0+ax*x+ay*y)
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Double_t Volume(Double_t a, Double_t sx, Double_t sy, Double_t ro)
This function calculates volume of a peak.
Double_t * fPositionInitX1
[fNPeaks] array of initial x positions of 1D ridges, input parameters
Double_t fSigmaErrY
error value of sigma y parameter
Double_t fAlpha
convergence coefficient, input parameter, it should be positive number and <=1, for details see refer...
Double_t Sqrt(Double_t x)
virtual ~TSpectrum2Fit()
Destructor.
Double_t fTyInit
initial value of t parameter for 1D ridges in y direction (relative amplitude of tail), for details see html manual and references
Double_t Dersigmax(Int_t numOfFittedPeaks, Double_t x, Double_t y, const Double_t *parameter, Double_t sigmax, Double_t sigmay, Double_t ro, Double_t txy, Double_t sxy, Double_t tx, Double_t sx, Double_t bx, Double_t by)
This function calculates derivative of peaks shape function (see manual) according to sigmax of peaks...
Double_t fRoInit
initial value of correlation coefficient
Double_t fTyErr
error value of t parameter for 1D ridges in y direction
Double_t Derbx(Int_t numOfFittedPeaks, Double_t x, Double_t y, const Double_t *parameter, Double_t sigmax, Double_t sigmay, Double_t txy, Double_t tx, Double_t bx, Double_t by)
This function calculates derivative of peaks shape function (see manual) according to slope bx...
Double_t * fAmpInitY1
[fNPeaks] array of initial values of amplitudes of 1D ridges in y direction, input parameters ...
Double_t fChi
here the fitting functions return resulting chi square
Double_t fSxErr
error value of s parameter for 1D ridges in x direction
Bool_t fFixBx
logical value of b parameter for 1D ridges in x direction, which allows to fix the parameter (not to ...
Double_t Dersy(Int_t numOfFittedPeaks, Double_t x, const Double_t *parameter, Double_t sigmax)
This function calculates derivative of peaks shape function (see manual) according to relative amplit...