154 if (numberPeaks <= 0){
155 Error (
"TSpectrum2Fit",
"Invalid number of peaks, must be > than 0");
302 Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap =
315 c =
c * t * (da1 + t * (da2 + t * da3));
327 Double_t da1 = 0.1740121, da2 = -0.0479399, da3 = 0.3739278, dap =
339 c = (-1.) * dap *
c * t * t * (da1 + t * (2. * da2 + t * 3. * da3)) -
357 if (pw > 10)
c *= a2;
358 if (pw > 12)
c *= a2;
377 Double_t sk = 0,
b, lambdak, normk, normk_old = 0;
383 for (i = 0; i <
size; i++) {
385 for (j = 0; j <
size; j++) {
393 sk = normk / normk_old;
397 for (i = 0; i <
size; i++) {
403 for (i = 0; i <
size; i++) {
404 for (j = 0,
b = 0; j <
size; j++) {
407 lambdak +=
b *
a[i][
size + 3];
410 lambdak = normk / lambdak;
414 for (i = 0; i <
size; i++)
445 Double_t r,
p, r1,
e,
ex,
ey, vx, s2, px, py, rx, ry, erx, ery;
448 for (j = 0; j < numOfFittedPeaks; j++) {
449 p = (
x - parameter[7 * j + 1]) / sigmax;
450 r = (
y - parameter[7 * j + 2]) / sigmay;
452 e = (
p *
p - 2 * ro *
p *
r +
r *
r) / (2 * (1 - ro * ro));
461 erx =
Erfc(
p / s2 + 1 / (2 * bx)), ery =
462 Erfc(
r / s2 + 1 / (2 * by));
463 ex =
p / (s2 * bx),
ey =
r / (s2 * by);
465 px = exp(
ex) * erx, py = exp(
ey) * ery;
467 r1 += 0.5 * txy * px * py;
471 r1 += 0.5 * sxy * rx * ry;
473 vx = vx + parameter[7 * j] * r1;
475 p = (
x - parameter[7 * j + 5]) / sigmax;
486 erx =
Erfc(
p / s2 + 1 / (2 * bx));
497 vx = vx + parameter[7 * j + 3] * r1;
499 r = (
y - parameter[7 * j + 6]) / sigmay;
510 ery =
Erfc(
r / s2 + 1 / (2 * by));
521 vx = vx + parameter[7 * j + 4] * r1;
524 vx = vx + a0 + ax *
x + ay *
y;
547 Double_t p,
r, r1 = 0,
e,
ex,
ey, px, py, rx, ry, erx, ery, s2;
548 p = (
x - x0) / sigmax;
549 r = (
y - y0) / sigmay;
552 e = (
p *
p - 2 * ro *
p *
r +
r *
r) / (2 * (1 - ro * ro));
561 erx =
Erfc(
p / s2 + 1 / (2 * bx)), ery =
562 Erfc(
r / s2 + 1 / (2 * by));
563 ex =
p / (s2 * bx),
ey =
r / (s2 * by);
565 px = exp(
ex) * erx, py = exp(
ey) * ery;
567 r1 += 0.5 * txy * px * py;
571 r1 += 0.5 * sxy * rx * ry;
594 p = (
x - x0) / sigmax;
606 erx =
Erfc(
p / s2 + 1 / (2 * bx));
642 Double_t p,
r, r1 = 0,
e,
ex,
ey, px, py, rx, ry, erx, ery, s2;
643 p = (
x - x0) / sigmax;
644 r = (
y - y0) / sigmay;
647 e = (
p *
p - 2 * ro *
p *
r +
r *
r) / (2 * (1 - ro * ro));
654 e = -(ro *
r -
p) / sigmax;
655 e =
e / (1 - ro * ro);
660 (-
Erfc(
p / s2 + 1 / (2 * bx)) / (s2 * bx * sigmax) -
661 Derfc(
p / s2 + 1 / (2 * bx)) / (s2 * sigmax)), ery =
662 Erfc(
r / s2 + 1 / (2 * by));
663 ex =
p / (s2 * bx),
ey =
r / (s2 * by);
665 px = exp(
ex) * erx, py = exp(
ey) * ery;
667 r1 += 0.5 * txy * px * py;
670 rx = -
Derfc(
p / s2) / (s2 * sigmax), ry =
Erfc(
r / s2);
671 r1 += 0.5 * sxy * rx * ry;
697 p = (
x - x0) / sigmax;
698 r = (
y - y0) / sigmay;
700 e = (
p *
p - 2 * ro *
p *
r +
r *
r) / (2 * (1 - ro * ro));
707 e = -(ro *
r -
p) / sigmax;
708 e =
e / (1 - ro * ro);
709 r1 = r1 * (
e *
e - 1 / ((1 - ro * ro) * sigmax * sigmax));
737 Double_t p,
r, r1 = 0,
e,
ex,
ey, px, py, rx, ry, erx, ery, s2;
738 p = (
x - x0) / sigmax;
739 r = (
y - y0) / sigmay;
742 e = (
p *
p - 2 * ro *
p *
r +
r *
r) / (2 * (1 - ro * ro));
749 e = -(ro *
p -
r) / sigmay;
750 e =
e / (1 - ro * ro);
755 (-
Erfc(
r / s2 + 1 / (2 * by)) / (s2 * by * sigmay) -
756 Derfc(
r / s2 + 1 / (2 * by)) / (s2 * sigmay)), erx =
757 Erfc(
p / s2 + 1 / (2 * bx));
758 ex =
p / (s2 * bx),
ey =
r / (s2 * by);
760 px = exp(
ex) * erx, py = exp(
ey) * ery;
762 r1 += 0.5 * txy * px * py;
765 ry = -
Derfc(
r / s2) / (s2 * sigmay), rx =
Erfc(
p / s2);
766 r1 += 0.5 * sxy * rx * ry;
792 p = (
x - x0) / sigmax;
793 r = (
y - y0) / sigmay;
795 e = (
p *
p - 2 * ro *
p *
r +
r *
r) / (2 * (1 - ro * ro));
802 e = -(ro *
p -
r) / sigmay;
803 e =
e / (1 - ro * ro);
804 r1 = r1 * (
e *
e - 1 / ((1 - ro * ro) * sigmay * sigmay));
827 p = (
x - x0) / sigmax;
837 r1 = r1 *
p / sigmax;
841 (-
Erfc(
p / s2 + 1 / (2 * bx)) / (s2 * bx * sigmax) -
842 Derfc(
p / s2 + 1 / (2 * bx)) / (s2 * sigmax));
849 rx = -
Derfc(
p / s2) / (s2 * sigmax);
871 p = (
x - x0) / sigmax;
880 r1 = r1 * (
p *
p / (sigmax * sigmax) - 1 / (sigmax * sigmax));
907 0,
e,
a,
b, x0, y0, s2, px, py, rx, ry, erx, ery,
ex,
ey;
910 for (j = 0; j < numOfFittedPeaks; j++) {
911 a = parameter[7 * j];
912 x0 = parameter[7 * j + 1];
913 y0 = parameter[7 * j + 2];
914 p = (
x - x0) / sigmax;
915 r = (
y - y0) / sigmay;
917 e = (
p *
p - 2 * ro *
p *
r +
r *
r) / (2 * (1 - ro * ro));
924 b = -(ro *
p *
r -
p *
p) / sigmax;
925 e =
e *
b / (1 - ro * ro);
929 -
Erfc(
p / s2 + 1 / (2 * bx)) *
p / (s2 * bx * sigmax) -
930 Derfc(
p / s2 + 1 / (2 * bx)) *
p / (s2 * sigmax), ery =
931 Erfc(
r / s2 + 1 / (2 * by));
932 ex =
p / (s2 * bx),
ey =
r / (s2 * by);
934 px = exp(
ex) * erx, py = exp(
ey) * ery;
936 e += 0.5 * txy * px * py;
939 rx = -
Derfc(
p / s2) *
p / (s2 * sigmax), ry =
Erfc(
r / s2);
940 e += 0.5 * sxy * rx * ry;
945 x0 = parameter[7 * j + 5];
946 p = (
x - x0) / sigmax;
954 e = 2 *
b *
e / sigmax;
958 (-
Erfc(
p / s2 + 1 / (2 * bx)) *
p / (s2 * bx * sigmax) -
959 Derfc(
p / s2 + 1 / (2 * bx)) *
p / (s2 * sigmax));
966 rx = -
Derfc(
p / s2) *
p / (s2 * sigmax);
969 r1 += parameter[7 * j + 3] *
e;
994 for (j = 0; j < numOfFittedPeaks; j++) {
995 a = parameter[7 * j];
996 x0 = parameter[7 * j + 1];
997 y0 = parameter[7 * j + 2];
998 p = (
x - x0) / sigmax;
999 r = (
y - y0) / sigmay;
1001 e = (
p *
p - 2 * ro *
p *
r +
r *
r) / (2 * (1 - ro * ro));
1008 b = -(ro *
p *
r -
p *
p) / sigmax;
1009 e =
e * (
b *
b / (1 - ro * ro) -
1010 (3 *
p *
p - 2 * ro *
p *
r) / (sigmax * sigmax)) / (1 -
1017 x0 = parameter[7 * j + 5];
1018 p = (
x - x0) / sigmax;
1026 e =
e * (4 *
b *
b - 6 *
b) / (sigmax * sigmax);
1027 r1 += parameter[7 * j + 3] *
e;
1054 0,
e,
a,
b, x0, y0, s2, px, py, rx, ry, erx, ery,
ex,
ey;
1057 for (j = 0; j < numOfFittedPeaks; j++) {
1058 a = parameter[7 * j];
1059 x0 = parameter[7 * j + 1];
1060 y0 = parameter[7 * j + 2];
1061 p = (
x - x0) / sigmax;
1062 r = (
y - y0) / sigmay;
1064 e = (
p *
p - 2 * ro *
p *
r +
r *
r) / (2 * (1 - ro * ro));
1071 b = -(ro *
p *
r -
r *
r) / sigmay;
1072 e =
e *
b / (1 - ro * ro);
1076 -
Erfc(
r / s2 + 1 / (2 * by)) *
r / (s2 * by * sigmay) -
1077 Derfc(
r / s2 + 1 / (2 * by)) *
r / (s2 * sigmay), erx =
1078 Erfc(
p / s2 + 1 / (2 * bx));
1079 ex =
p / (s2 * bx),
ey =
r / (s2 * by);
1081 px = exp(
ex) * erx, py = exp(
ey) * ery;
1083 e += 0.5 * txy * px * py;
1086 ry = -
Derfc(
r / s2) *
r / (s2 * sigmay), rx =
Erfc(
p / s2);
1087 e += 0.5 * sxy * rx * ry;
1092 y0 = parameter[7 * j + 6];
1093 r = (
y - y0) / sigmay;
1101 e = 2 *
b *
e / sigmay;
1105 (-
Erfc(
r / s2 + 1 / (2 * by)) *
r / (s2 * by * sigmay) -
1106 Derfc(
r / s2 + 1 / (2 * by)) *
r / (s2 * sigmay));
1113 ry = -
Derfc(
r / s2) *
r / (s2 * sigmay);
1116 r1 += parameter[7 * j + 4] *
e;
1141 for (j = 0; j < numOfFittedPeaks; j++) {
1142 a = parameter[7 * j];
1143 x0 = parameter[7 * j + 1];
1144 y0 = parameter[7 * j + 2];
1145 p = (
x - x0) / sigmax;
1146 r = (
y - y0) / sigmay;
1148 e = (
p *
p - 2 * ro *
p *
r +
r *
r) / (2 * (1 - ro * ro));
1155 b = -(ro *
p *
r -
r *
r) / sigmay;
1156 e =
e * (
b *
b / (1 - ro * ro) -
1157 (3 *
r *
r - 2 * ro *
r *
p) / (sigmay * sigmay)) / (1 -
1164 y0 = parameter[7 * j + 6];
1165 r = (
y - y0) / sigmay;
1173 e =
e * (4 *
b *
b - 6 *
b) / (sigmay * sigmay);
1174 r1 += parameter[7 * j + 4] *
e;
1199 for (j = 0; j < numOfFittedPeaks; j++) {
1200 a = parameter[7 * j];
1201 x0 = parameter[7 * j + 1];
1202 y0 = parameter[7 * j + 2];
1206 rx = (px * px - 2 *
r * px * qx + qx * qx);
1207 ex = rx / (2 * (1 -
r *
r));
1214 tx = px * qx / (1 -
r *
r);
1215 tx = tx -
r * rx / ((1 -
r *
r) * (1 -
r *
r));
1216 vx = vx +
a *
ex * tx;
1238 Double_t p,
r, r1 = 0,
ex,
ey, px, py, erx, ery, s2, x0, y0,
a;
1241 for (j = 0; j < numOfFittedPeaks; j++) {
1242 a = parameter[7 * j];
1243 x0 = parameter[7 * j + 1];
1244 y0 = parameter[7 * j + 2];
1245 p = (
x - x0) / sigmax;
1246 r = (
y - y0) / sigmay;
1248 erx =
Erfc(
p / s2 + 1 / (2 * bx)), ery =
1249 Erfc(
r / s2 + 1 / (2 * by));
1250 ex =
p / (s2 * bx),
ey =
r / (s2 * by);
1252 px = exp(
ex) * erx, py = exp(
ey) * ery;
1254 r1 += 0.5 *
a * px * py;
1277 for (j = 0; j < numOfFittedPeaks; j++) {
1278 a = parameter[7 * j];
1279 x0 = parameter[7 * j + 1];
1280 y0 = parameter[7 * j + 2];
1281 p = (
x - x0) / sigmax;
1282 r = (
y - y0) / sigmay;
1284 r1 += 0.5 *
a * rx * ry;
1307 for (j = 0; j < numOfFittedPeaks; j++) {
1308 ax = parameter[7 * j + 3];
1309 x0 = parameter[7 * j + 5];
1310 p = (
x - x0) / sigmax;
1312 erx =
Erfc(
p / s2 + 1 / (2 * bx));
1317 r1 += 0.5 * ax * px;
1340 for (j = 0; j < numOfFittedPeaks; j++) {
1341 ax = parameter[7 * j + 4];
1342 x0 = parameter[7 * j + 6];
1343 p = (
x - x0) / sigmax;
1345 erx =
Erfc(
p / s2 + 1 / (2 * bx));
1350 r1 += 0.5 * ax * px;
1371 for (j = 0; j < numOfFittedPeaks; j++) {
1372 ax = parameter[7 * j + 3];
1373 x0 = parameter[7 * j + 5];
1374 p = (
x - x0) / sigmax;
1377 r1 += 0.5 * ax * rx;
1398 for (j = 0; j < numOfFittedPeaks; j++) {
1399 ax = parameter[7 * j + 4];
1400 x0 = parameter[7 * j + 6];
1401 p = (
x - x0) / sigmax;
1404 r1 += 0.5 * ax * rx;
1427 Double_t p,
r, r1 = 0,
a, x0, y0, s2, px, py, erx, ery,
ex,
ey;
1430 for (j = 0; j < numOfFittedPeaks; j++) {
1431 a = parameter[7 * j];
1432 x0 = parameter[7 * j + 1];
1433 y0 = parameter[7 * j + 2];
1434 p = (
x - x0) / sigmax;
1435 r = (
y - y0) / sigmay;
1439 -
Erfc(
p / s2 + 1 / (2 * bx)) *
p / (s2 * bx * bx) -
1440 Derfc(
p / s2 + 1 / (2 * bx)) / (s2 * bx * bx), ery =
1441 Erfc(
r / s2 + 1 / (2 * by));
1442 ex =
p / (s2 * bx),
ey =
r / (s2 * by);
1444 px = exp(
ex) * erx, py = exp(
ey) * ery;
1446 r1 += 0.5 *
a * txy * px * py;
1448 a = parameter[7 * j + 3];
1449 x0 = parameter[7 * j + 5];
1450 p = (
x - x0) / sigmax;
1454 (-
Erfc(
p / s2 + 1 / (2 * bx)) *
p / (s2 * bx * bx) -
1455 Derfc(
p / s2 + 1 / (2 * bx)) / (s2 * bx * bx));
1459 r1 += 0.5 *
a * tx * px;
1483 Double_t p,
r, r1 = 0,
a, x0, y0, s2, px, py, erx, ery,
ex,
ey;
1486 for (j = 0; j < numOfFittedPeaks; j++) {
1487 a = parameter[7 * j];
1488 x0 = parameter[7 * j + 1];
1489 y0 = parameter[7 * j + 2];
1490 p = (
x - x0) / sigmax;
1491 r = (
y - y0) / sigmay;
1495 -
Erfc(
r / s2 + 1 / (2 * by)) *
r / (s2 * by * by) -
1496 Derfc(
r / s2 + 1 / (2 * by)) / (s2 * by * by), erx =
1497 Erfc(
p / s2 + 1 / (2 * bx));
1498 ex =
p / (s2 * bx),
ey =
r / (s2 * by);
1500 px = exp(
ex) * erx, py = exp(
ey) * ery;
1502 r1 += 0.5 *
a * txy * px * py;
1504 a = parameter[7 * j + 4];
1505 y0 = parameter[7 * j + 6];
1506 r = (
y - y0) / sigmay;
1510 (-
Erfc(
r / s2 + 1 / (2 * by)) *
r / (s2 * by * by) -
1511 Derfc(
r / s2 + 1 / (2 * by)) / (s2 * by * by));
1515 r1 += 0.5 *
a * ty * py;
1539 r = 2 *
a * pi * sx * sy *
r;
1561 r = 2 * pi * sx * sy *
r;
1584 r =
a * 2 * pi * sy *
r;
1607 r =
a * 2 * pi * sx *
r;
1630 r = -
a * 2 * pi * sx * sy * ro /
r;
1853 Int_t i, i1, i2, j, k, shift =
1856 Double_t a,
b,
c,
d = 0, alpha, chi_opt, yw, ywm,
f, chi2, chi_min, chi =
1857 0, pi, pmin = 0, chi_cel = 0, chi_er;
1859 for (i = 0, j = 0; i <
fNPeaks; i++) {
1860 working_space[7 * i] =
fAmpInit[i];
1862 working_space[shift + j] =
fAmpInit[i];
1907 working_space[7 * i + 2] =
fRoInit;
1909 working_space[shift + j] =
fRoInit;
1912 working_space[7 * i + 3] =
fA0Init;
1914 working_space[shift + j] =
fA0Init;
1917 working_space[7 * i + 4] =
fAxInit;
1919 working_space[shift + j] =
fAxInit;
1922 working_space[7 * i + 5] =
fAyInit;
1924 working_space[shift + j] =
fAyInit;
1927 working_space[7 * i + 6] =
fTxyInit;
1929 working_space[shift + j] =
fTxyInit;
1932 working_space[7 * i + 7] =
fSxyInit;
1934 working_space[shift + j] =
fSxyInit;
1937 working_space[7 * i + 8] =
fTxInit;
1939 working_space[shift + j] =
fTxInit;
1942 working_space[7 * i + 9] =
fTyInit;
1944 working_space[shift + j] =
fTyInit;
1947 working_space[7 * i + 10] =
fSxyInit;
1949 working_space[shift + j] =
fSxInit;
1952 working_space[7 * i + 11] =
fSyInit;
1954 working_space[shift + j] =
fSyInit;
1957 working_space[7 * i + 12] =
fBxInit;
1959 working_space[shift + j] =
fBxInit;
1962 working_space[7 * i + 13] =
fByInit;
1964 working_space[shift + j] =
fByInit;
1969 for (j = 0; j <
size; j++) {
1970 working_space[2 * shift + j] = 0, working_space[3 * shift + j] = 0;
1975 chi_opt = 0, pw =
fPower - 2;
1978 yw = source[i1][i2];
1981 working_space, working_space[peak_vel],
1982 working_space[peak_vel + 1],
1983 working_space[peak_vel + 2],
1984 working_space[peak_vel + 3],
1985 working_space[peak_vel + 4],
1986 working_space[peak_vel + 5],
1987 working_space[peak_vel + 6],
1988 working_space[peak_vel + 7],
1989 working_space[peak_vel + 8],
1990 working_space[peak_vel + 9],
1991 working_space[peak_vel + 10],
1992 working_space[peak_vel + 11],
1993 working_space[peak_vel + 12],
1994 working_space[peak_vel + 13]);
2002 chi_opt += (yw -
f) * (yw -
f) / ywm;
2022 for (j = 0, k = 0; j <
fNPeaks; j++) {
2025 working_space[7 * j + 1],
2026 working_space[7 * j + 2],
2027 working_space[peak_vel],
2028 working_space[peak_vel + 1],
2029 working_space[peak_vel + 2],
2030 working_space[peak_vel + 6],
2031 working_space[peak_vel + 7],
2032 working_space[peak_vel + 12],
2033 working_space[peak_vel + 13]);
2037 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2038 working_space[2 * shift + k] +=
b *
c;
2039 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2040 working_space[3 * shift + k] +=
b *
c;
2044 b =
a * (yw -
f) / ywm;
2045 working_space[2 * shift + k] +=
b *
c;
2047 working_space[3 * shift + k] +=
b *
c;
2054 working_space[7 * j],
2055 working_space[7 * j + 1],
2056 working_space[7 * j + 2],
2057 working_space[peak_vel],
2058 working_space[peak_vel + 1],
2059 working_space[peak_vel + 2],
2060 working_space[peak_vel + 6],
2061 working_space[peak_vel + 7],
2062 working_space[peak_vel + 12],
2063 working_space[peak_vel + 13]);
2066 working_space[7 * j],
2067 working_space[7 * j + 1],
2068 working_space[7 * j + 2],
2069 working_space[peak_vel],
2070 working_space[peak_vel + 1],
2071 working_space[peak_vel + 2]);
2077 if (((
a +
d) <= 0 &&
a >= 0) || ((
a +
d) >= 0
2086 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2087 working_space[2 * shift + k] +=
b *
c;
2088 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2089 working_space[3 * shift + k] +=
b *
c;
2093 b =
a * (yw -
f) / ywm;
2094 working_space[2 * shift + k] +=
b *
c;
2096 working_space[3 * shift + k] +=
b *
c;
2103 working_space[7 * j],
2104 working_space[7 * j + 1],
2105 working_space[7 * j + 2],
2106 working_space[peak_vel],
2107 working_space[peak_vel + 1],
2108 working_space[peak_vel + 2],
2109 working_space[peak_vel + 6],
2110 working_space[peak_vel + 7],
2111 working_space[peak_vel + 12],
2112 working_space[peak_vel + 13]);
2115 working_space[7 * j],
2116 working_space[7 * j + 1],
2117 working_space[7 * j + 2],
2118 working_space[peak_vel],
2119 working_space[peak_vel + 1],
2120 working_space[peak_vel + 2]);
2126 if (((
a +
d) <= 0 &&
a >= 0) || ((
a +
d) >= 0
2135 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2136 working_space[2 * shift + k] +=
b *
c;
2137 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2138 working_space[3 * shift + k] +=
b *
c;
2142 b =
a * (yw -
f) / ywm;
2143 working_space[2 * shift + k] +=
b *
c;
2145 working_space[3 * shift + k] +=
b *
c;
2151 a =
Derampx(i1, working_space[7 * j + 5],
2152 working_space[peak_vel],
2153 working_space[peak_vel + 8],
2154 working_space[peak_vel + 10],
2155 working_space[peak_vel + 12]);
2159 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2160 working_space[2 * shift + k] +=
b *
c;
2161 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2162 working_space[3 * shift + k] +=
b *
c;
2166 b =
a * (yw -
f) / ywm;
2167 working_space[2 * shift + k] +=
b *
c;
2169 working_space[3 * shift + k] +=
b *
c;
2175 a =
Derampx(i2, working_space[7 * j + 6],
2176 working_space[peak_vel + 1],
2177 working_space[peak_vel + 9],
2178 working_space[peak_vel + 11],
2179 working_space[peak_vel + 13]);
2183 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2184 working_space[2 * shift + k] +=
b *
c;
2185 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2186 working_space[3 * shift + k] +=
b *
c;
2190 b =
a * (yw -
f) / ywm;
2191 working_space[2 * shift + k] +=
b *
c;
2193 working_space[3 * shift + k] +=
b *
c;
2199 a =
Deri01(i1, working_space[7 * j + 3],
2200 working_space[7 * j + 5],
2201 working_space[peak_vel],
2202 working_space[peak_vel + 8],
2203 working_space[peak_vel + 10],
2204 working_space[peak_vel + 12]);
2207 working_space[7 * j + 5],
2208 working_space[peak_vel]);
2214 if (((
a +
d) <= 0 &&
a >= 0) || ((
a +
d) >= 0
2223 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2224 working_space[2 * shift + k] +=
b *
c;
2225 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2226 working_space[3 * shift + k] +=
b *
c;
2230 b =
a * (yw -
f) / ywm;
2231 working_space[2 * shift + k] +=
b *
c;
2233 working_space[3 * shift + k] +=
b *
c;
2239 a =
Deri01(i2, working_space[7 * j + 4],
2240 working_space[7 * j + 6],
2241 working_space[peak_vel + 1],
2242 working_space[peak_vel + 9],
2243 working_space[peak_vel + 11],
2244 working_space[peak_vel + 13]);
2247 working_space[7 * j + 6],
2248 working_space[peak_vel + 1]);
2254 if (((
a +
d) <= 0 &&
a >= 0) || ((
a +
d) >= 0
2263 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2264 working_space[2 * shift + k] +=
b *
c;
2265 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2266 working_space[3 * shift + k] +=
b *
c;
2270 b =
a * (yw -
f) / ywm;
2271 working_space[2 * shift + k] +=
b *
c;
2273 working_space[3 * shift + k] +=
b *
c;
2281 working_space, working_space[peak_vel],
2282 working_space[peak_vel + 1],
2283 working_space[peak_vel + 2],
2284 working_space[peak_vel + 6],
2285 working_space[peak_vel + 7],
2286 working_space[peak_vel + 8],
2287 working_space[peak_vel + 10],
2288 working_space[peak_vel + 12],
2289 working_space[peak_vel + 13]);
2293 working_space[peak_vel],
2294 working_space[peak_vel + 1],
2295 working_space[peak_vel + 2]);
2301 if (((
a +
d) <= 0 &&
a >= 0) || ((
a +
d) >= 0 &&
a <= 0))
2309 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2310 working_space[2 * shift + k] +=
b *
c;
2311 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2312 working_space[3 * shift + k] +=
b *
c;
2316 b =
a * (yw -
f) / ywm;
2317 working_space[2 * shift + k] +=
b *
c;
2319 working_space[3 * shift + k] +=
b *
c;
2326 working_space, working_space[peak_vel],
2327 working_space[peak_vel + 1],
2328 working_space[peak_vel + 2],
2329 working_space[peak_vel + 6],
2330 working_space[peak_vel + 7],
2331 working_space[peak_vel + 9],
2332 working_space[peak_vel + 11],
2333 working_space[peak_vel + 12],
2334 working_space[peak_vel + 13]);
2338 working_space[peak_vel],
2339 working_space[peak_vel + 1],
2340 working_space[peak_vel + 2]);
2346 if (((
a +
d) <= 0 &&
a >= 0) || ((
a +
d) >= 0 &&
a <= 0))
2354 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2355 working_space[2 * shift + k] +=
b *
c;
2356 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2357 working_space[3 * shift + k] +=
b *
c;
2361 b =
a * (yw -
f) / ywm;
2362 working_space[2 * shift + k] +=
b *
c;
2364 working_space[3 * shift + k] +=
b *
c;
2371 working_space, working_space[peak_vel],
2372 working_space[peak_vel + 1],
2373 working_space[peak_vel + 2]);
2379 if (((
a +
d) <= 0 &&
a >= 0) || ((
a +
d) >= 0 &&
a <= 0))
2387 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2388 working_space[2 * shift + k] +=
b *
c;
2389 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2390 working_space[3 * shift + k] +=
b *
c;
2394 b =
a * (yw -
f) / ywm;
2395 working_space[2 * shift + k] +=
b *
c;
2397 working_space[3 * shift + k] +=
b *
c;
2407 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2408 working_space[2 * shift + k] +=
b *
c;
2409 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2410 working_space[3 * shift + k] +=
b *
c;
2414 b =
a * (yw -
f) / ywm;
2415 working_space[2 * shift + k] +=
b *
c;
2417 working_space[3 * shift + k] +=
b *
c;
2427 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2428 working_space[2 * shift + k] +=
b *
c;
2429 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2430 working_space[3 * shift + k] +=
b *
c;
2434 b =
a * (yw -
f) / ywm;
2435 working_space[2 * shift + k] +=
b *
c;
2437 working_space[3 * shift + k] +=
b *
c;
2447 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2448 working_space[2 * shift + k] +=
b *
c;
2449 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2450 working_space[3 * shift + k] +=
b *
c;
2454 b =
a * (yw -
f) / ywm;
2455 working_space[2 * shift + k] +=
b *
c;
2457 working_space[3 * shift + k] +=
b *
c;
2464 working_space, working_space[peak_vel],
2465 working_space[peak_vel + 1],
2466 working_space[peak_vel + 12],
2467 working_space[peak_vel + 13]);
2471 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2472 working_space[2 * shift + k] +=
b *
c;
2473 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2474 working_space[3 * shift + k] +=
b *
c;
2478 b =
a * (yw -
f) / ywm;
2479 working_space[2 * shift + k] +=
b *
c;
2481 working_space[3 * shift + k] +=
b *
c;
2488 working_space, working_space[peak_vel],
2489 working_space[peak_vel + 1]);
2493 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2494 working_space[2 * shift + k] +=
b *
c;
2495 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2496 working_space[3 * shift + k] +=
b *
c;
2500 b =
a * (yw -
f) / ywm;
2501 working_space[2 * shift + k] +=
b *
c;
2503 working_space[3 * shift + k] +=
b *
c;
2510 working_space[peak_vel],
2511 working_space[peak_vel + 12]);
2515 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2516 working_space[2 * shift + k] +=
b *
c;
2517 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2518 working_space[3 * shift + k] +=
b *
c;
2522 b =
a * (yw -
f) / ywm;
2523 working_space[2 * shift + k] +=
b *
c;
2525 working_space[3 * shift + k] +=
b *
c;
2532 working_space[peak_vel + 1],
2533 working_space[peak_vel + 13]);
2537 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2538 working_space[2 * shift + k] +=
b *
c;
2539 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2540 working_space[3 * shift + k] +=
b *
c;
2544 b =
a * (yw -
f) / ywm;
2545 working_space[2 * shift + k] +=
b *
c;
2547 working_space[3 * shift + k] +=
b *
c;
2554 working_space[peak_vel]);
2558 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2559 working_space[2 * shift + k] +=
b *
c;
2560 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2561 working_space[3 * shift + k] +=
b *
c;
2565 b =
a * (yw -
f) / ywm;
2566 working_space[2 * shift + k] +=
b *
c;
2568 working_space[3 * shift + k] +=
b *
c;
2575 working_space[peak_vel + 1]);
2579 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2580 working_space[2 * shift + k] +=
b *
c;
2581 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2582 working_space[3 * shift + k] +=
b *
c;
2586 b =
a * (yw -
f) / ywm;
2587 working_space[2 * shift + k] +=
b *
c;
2589 working_space[3 * shift + k] +=
b *
c;
2596 working_space, working_space[peak_vel],
2597 working_space[peak_vel + 1],
2598 working_space[peak_vel + 6],
2599 working_space[peak_vel + 8],
2600 working_space[peak_vel + 12],
2601 working_space[peak_vel + 13]);
2605 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2606 working_space[2 * shift + k] +=
b *
c;
2607 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2608 working_space[3 * shift + k] +=
b *
c;
2612 b =
a * (yw -
f) / ywm;
2613 working_space[2 * shift + k] +=
b *
c;
2615 working_space[3 * shift + k] +=
b *
c;
2622 working_space, working_space[peak_vel],
2623 working_space[peak_vel + 1],
2624 working_space[peak_vel + 6],
2625 working_space[peak_vel + 8],
2626 working_space[peak_vel + 12],
2627 working_space[peak_vel + 13]);
2631 b =
a * (yw * yw -
f *
f) / (ywm * ywm);
2632 working_space[2 * shift + k] +=
b *
c;
2633 b =
a *
a * (4 * yw - 2 *
f) / (ywm * ywm);
2634 working_space[3 * shift + k] +=
b *
c;
2638 b =
a * (yw -
f) / ywm;
2639 working_space[2 * shift + k] +=
b *
c;
2641 working_space[3 * shift + k] +=
b *
c;
2648 for (j = 0; j <
size; j++) {
2649 if (
TMath::Abs(working_space[3 * shift + j]) > 0.000001)
2650 working_space[2 * shift + j] = working_space[2 * shift + j] /
TMath::Abs(working_space[3 * shift + j]);
2652 working_space[2 * shift + j] = 0;
2661 for (j = 0; j <
size; j++) {
2662 working_space[4 * shift + j] = working_space[shift + j];
2668 chi_min = 10000 * chi2;
2671 chi_min = 0.1 * chi2;
2673 for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
2674 for (j = 0; j <
size; j++) {
2675 working_space[shift + j] = working_space[4 * shift + j] + pi * alpha * working_space[2 * shift + j];
2677 for (i = 0, j = 0; i <
fNPeaks; i++) {
2679 if (working_space[shift + j] < 0)
2680 working_space[shift + j] = 0;
2681 working_space[7 * i] = working_space[shift + j];
2685 if (working_space[shift + j] <
fXmin)
2686 working_space[shift + j] =
fXmin;
2687 if (working_space[shift + j] >
fXmax)
2688 working_space[shift + j] =
fXmax;
2689 working_space[7 * i + 1] = working_space[shift + j];
2693 if (working_space[shift + j] <
fYmin)
2694 working_space[shift + j] =
fYmin;
2695 if (working_space[shift + j] >
fYmax)
2696 working_space[shift + j] =
fYmax;
2697 working_space[7 * i + 2] = working_space[shift + j];
2701 if (working_space[shift + j] < 0)
2702 working_space[shift + j] = 0;
2703 working_space[7 * i + 3] = working_space[shift + j];
2707 if (working_space[shift + j] < 0)
2708 working_space[shift + j] = 0;
2709 working_space[7 * i + 4] = working_space[shift + j];
2713 if (working_space[shift + j] <
fXmin)
2714 working_space[shift + j] =
fXmin;
2715 if (working_space[shift + j] >
fXmax)
2716 working_space[shift + j] =
fXmax;
2717 working_space[7 * i + 5] = working_space[shift + j];
2721 if (working_space[shift + j] <
fYmin)
2722 working_space[shift + j] =
fYmin;
2723 if (working_space[shift + j] >
fYmax)
2724 working_space[shift + j] =
fYmax;
2725 working_space[7 * i + 6] = working_space[shift + j];
2730 if (working_space[shift + j] < 0.001) {
2731 working_space[shift + j] = 0.001;
2733 working_space[peak_vel] = working_space[shift + j];
2737 if (working_space[shift + j] < 0.001) {
2738 working_space[shift + j] = 0.001;
2740 working_space[peak_vel + 1] = working_space[shift + j];
2744 if (working_space[shift + j] < -1) {
2745 working_space[shift + j] = -1;
2747 if (working_space[shift + j] > 1) {
2748 working_space[shift + j] = 1;
2750 working_space[peak_vel + 2] = working_space[shift + j];
2754 working_space[peak_vel + 3] = working_space[shift + j];
2758 working_space[peak_vel + 4] = working_space[shift + j];
2762 working_space[peak_vel + 5] = working_space[shift + j];
2766 working_space[peak_vel + 6] = working_space[shift + j];
2770 working_space[peak_vel + 7] = working_space[shift + j];
2774 working_space[peak_vel + 8] = working_space[shift + j];
2778 working_space[peak_vel + 9] = working_space[shift + j];
2782 working_space[peak_vel + 10] = working_space[shift + j];
2786 working_space[peak_vel + 11] = working_space[shift + j];
2790 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
2791 if (working_space[shift + j] < 0)
2792 working_space[shift + j] = -0.001;
2794 working_space[shift + j] = 0.001;
2796 working_space[peak_vel + 12] = working_space[shift + j];
2800 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
2801 if (working_space[shift + j] < 0)
2802 working_space[shift + j] = -0.001;
2804 working_space[shift + j] = 0.001;
2806 working_space[peak_vel + 13] = working_space[shift + j];
2812 yw = source[i1][i2];
2816 working_space[peak_vel],
2817 working_space[peak_vel + 1],
2818 working_space[peak_vel + 2],
2819 working_space[peak_vel + 3],
2820 working_space[peak_vel + 4],
2821 working_space[peak_vel + 5],
2822 working_space[peak_vel + 6],
2823 working_space[peak_vel + 7],
2824 working_space[peak_vel + 8],
2825 working_space[peak_vel + 9],
2826 working_space[peak_vel + 10],
2827 working_space[peak_vel + 11],
2828 working_space[peak_vel + 12],
2829 working_space[peak_vel + 13]);
2842 chi2 += (yw -
f) * (yw -
f) / ywm;
2850 pmin = pi, chi_min = chi2;
2860 for (j = 0; j <
size; j++) {
2861 working_space[shift + j] = working_space[4 * shift + j] + pmin * alpha * working_space[2 * shift + j];
2863 for (i = 0, j = 0; i <
fNPeaks; i++) {
2865 if (working_space[shift + j] < 0)
2866 working_space[shift + j] = 0;
2867 working_space[7 * i] = working_space[shift + j];
2871 if (working_space[shift + j] <
fXmin)
2872 working_space[shift + j] =
fXmin;
2873 if (working_space[shift + j] >
fXmax)
2874 working_space[shift + j] =
fXmax;
2875 working_space[7 * i + 1] = working_space[shift + j];
2879 if (working_space[shift + j] <
fYmin)
2880 working_space[shift + j] =
fYmin;
2881 if (working_space[shift + j] >
fYmax)
2882 working_space[shift + j] =
fYmax;
2883 working_space[7 * i + 2] = working_space[shift + j];
2887 if (working_space[shift + j] < 0)
2888 working_space[shift + j] = 0;
2889 working_space[7 * i + 3] = working_space[shift + j];
2893 if (working_space[shift + j] < 0)
2894 working_space[shift + j] = 0;
2895 working_space[7 * i + 4] = working_space[shift + j];
2899 if (working_space[shift + j] <
fXmin)
2900 working_space[shift + j] =
fXmin;
2901 if (working_space[shift + j] >
fXmax)
2902 working_space[shift + j] =
fXmax;
2903 working_space[7 * i + 5] = working_space[shift + j];
2907 if (working_space[shift + j] <
fYmin)
2908 working_space[shift + j] =
fYmin;
2909 if (working_space[shift + j] >
fYmax)
2910 working_space[shift + j] =
fYmax;
2911 working_space[7 * i + 6] = working_space[shift + j];
2916 if (working_space[shift + j] < 0.001) {
2917 working_space[shift + j] = 0.001;
2919 working_space[peak_vel] = working_space[shift + j];
2923 if (working_space[shift + j] < 0.001) {
2924 working_space[shift + j] = 0.001;
2926 working_space[peak_vel + 1] = working_space[shift + j];
2930 if (working_space[shift + j] < -1) {
2931 working_space[shift + j] = -1;
2933 if (working_space[shift + j] > 1) {
2934 working_space[shift + j] = 1;
2936 working_space[peak_vel + 2] = working_space[shift + j];
2940 working_space[peak_vel + 3] = working_space[shift + j];
2944 working_space[peak_vel + 4] = working_space[shift + j];
2948 working_space[peak_vel + 5] = working_space[shift + j];
2952 working_space[peak_vel + 6] = working_space[shift + j];
2956 working_space[peak_vel + 7] = working_space[shift + j];
2960 working_space[peak_vel + 8] = working_space[shift + j];
2964 working_space[peak_vel + 9] = working_space[shift + j];
2968 working_space[peak_vel + 10] = working_space[shift + j];
2972 working_space[peak_vel + 11] = working_space[shift + j];
2976 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
2977 if (working_space[shift + j] < 0)
2978 working_space[shift + j] = -0.001;
2980 working_space[shift + j] = 0.001;
2982 working_space[peak_vel + 12] = working_space[shift + j];
2986 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
2987 if (working_space[shift + j] < 0)
2988 working_space[shift + j] = -0.001;
2990 working_space[shift + j] = 0.001;
2992 working_space[peak_vel + 13] = working_space[shift + j];
3000 for (j = 0; j <
size; j++) {
3001 working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j];
3003 for (i = 0, j = 0; i <
fNPeaks; i++) {
3005 if (working_space[shift + j] < 0)
3006 working_space[shift + j] = 0;
3007 working_space[7 * i] = working_space[shift + j];
3011 if (working_space[shift + j] <
fXmin)
3012 working_space[shift + j] =
fXmin;
3013 if (working_space[shift + j] >
fXmax)
3014 working_space[shift + j] =
fXmax;
3015 working_space[7 * i + 1] = working_space[shift + j];
3019 if (working_space[shift + j] <
fYmin)
3020 working_space[shift + j] =
fYmin;
3021 if (working_space[shift + j] >
fYmax)
3022 working_space[shift + j] =
fYmax;
3023 working_space[7 * i + 2] = working_space[shift + j];
3027 if (working_space[shift + j] < 0)
3028 working_space[shift + j] = 0;
3029 working_space[7 * i + 3] = working_space[shift + j];
3033 if (working_space[shift + j] < 0)
3034 working_space[shift + j] = 0;
3035 working_space[7 * i + 4] = working_space[shift + j];
3039 if (working_space[shift + j] <
fXmin)
3040 working_space[shift + j] =
fXmin;
3041 if (working_space[shift + j] >
fXmax)
3042 working_space[shift + j] =
fXmax;
3043 working_space[7 * i + 5] = working_space[shift + j];
3047 if (working_space[shift + j] <
fYmin)
3048 working_space[shift + j] =
fYmin;
3049 if (working_space[shift + j] >
fYmax)
3050 working_space[shift + j] =
fYmax;
3051 working_space[7 * i + 6] = working_space[shift + j];
3056 if (working_space[shift + j] < 0.001) {
3057 working_space[shift + j] = 0.001;
3059 working_space[peak_vel] = working_space[shift + j];
3063 if (working_space[shift + j] < 0.001) {
3064 working_space[shift + j] = 0.001;
3066 working_space[peak_vel + 1] = working_space[shift + j];
3070 if (working_space[shift + j] < -1) {
3071 working_space[shift + j] = -1;
3073 if (working_space[shift + j] > 1) {
3074 working_space[shift + j] = 1;
3076 working_space[peak_vel + 2] = working_space[shift + j];
3080 working_space[peak_vel + 3] = working_space[shift + j];
3084 working_space[peak_vel + 4] = working_space[shift + j];
3088 working_space[peak_vel + 5] = working_space[shift + j];
3092 working_space[peak_vel + 6] = working_space[shift + j];
3096 working_space[peak_vel + 7] = working_space[shift + j];
3100 working_space[peak_vel + 8] = working_space[shift + j];
3104 working_space[peak_vel + 9] = working_space[shift + j];
3108 working_space[peak_vel + 10] = working_space[shift + j];
3112 working_space[peak_vel + 11] = working_space[shift + j];
3116 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
3117 if (working_space[shift + j] < 0)
3118 working_space[shift + j] = -0.001;
3120 working_space[shift + j] = 0.001;
3122 working_space[peak_vel + 12] = working_space[shift + j];
3126 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
3127 if (working_space[shift + j] < 0)
3128 working_space[shift + j] = -0.001;
3130 working_space[shift + j] = 0.001;
3132 working_space[peak_vel + 13] = working_space[shift + j];
3138 yw = source[i1][i2];
3141 working_space, working_space[peak_vel],
3142 working_space[peak_vel + 1],
3143 working_space[peak_vel + 2],
3144 working_space[peak_vel + 3],
3145 working_space[peak_vel + 4],
3146 working_space[peak_vel + 5],
3147 working_space[peak_vel + 6],
3148 working_space[peak_vel + 7],
3149 working_space[peak_vel + 8],
3150 working_space[peak_vel + 9],
3151 working_space[peak_vel + 10],
3152 working_space[peak_vel + 11],
3153 working_space[peak_vel + 12],
3154 working_space[peak_vel + 13]);
3167 chi += (yw -
f) * (yw -
f) / ywm;
3175 alpha = alpha * chi_opt / (2 * chi);
3178 alpha = alpha / 10.0;
3181 }
while (((chi > chi_opt
3186 for (j = 0; j <
size; j++) {
3187 working_space[4 * shift + j] = 0;
3188 working_space[2 * shift + j] = 0;
3190 for (i1 =
fXmin, chi_cel = 0; i1 <=
fXmax; i1++) {
3192 yw = source[i1][i2];
3196 working_space, working_space[peak_vel],
3197 working_space[peak_vel + 1],
3198 working_space[peak_vel + 2],
3199 working_space[peak_vel + 3],
3200 working_space[peak_vel + 4],
3201 working_space[peak_vel + 5],
3202 working_space[peak_vel + 6],
3203 working_space[peak_vel + 7],
3204 working_space[peak_vel + 8],
3205 working_space[peak_vel + 9],
3206 working_space[peak_vel + 10],
3207 working_space[peak_vel + 11],
3208 working_space[peak_vel + 12],
3209 working_space[peak_vel + 13]);
3210 chi_opt = (yw -
f) * (yw -
f) / yw;
3211 chi_cel += (yw -
f) * (yw -
f) / yw;
3214 for (j = 0, k = 0; j <
fNPeaks; j++) {
3217 working_space[7 * j + 1],
3218 working_space[7 * j + 2],
3219 working_space[peak_vel],
3220 working_space[peak_vel + 1],
3221 working_space[peak_vel + 2],
3222 working_space[peak_vel + 6],
3223 working_space[peak_vel + 7],
3224 working_space[peak_vel + 12],
3225 working_space[peak_vel + 13]);
3228 working_space[2 * shift + k] += chi_opt *
c;
3230 working_space[4 * shift + k] +=
b *
c;
3236 working_space[7 * j],
3237 working_space[7 * j + 1],
3238 working_space[7 * j + 2],
3239 working_space[peak_vel],
3240 working_space[peak_vel + 1],
3241 working_space[peak_vel + 2],
3242 working_space[peak_vel + 6],
3243 working_space[peak_vel + 7],
3244 working_space[peak_vel + 12],
3245 working_space[peak_vel + 13]);
3248 working_space[2 * shift + k] += chi_opt *
c;
3250 working_space[4 * shift + k] +=
b *
c;
3256 working_space[7 * j],
3257 working_space[7 * j + 1],
3258 working_space[7 * j + 2],
3259 working_space[peak_vel],
3260 working_space[peak_vel + 1],
3261 working_space[peak_vel + 2],
3262 working_space[peak_vel + 6],
3263 working_space[peak_vel + 7],
3264 working_space[peak_vel + 12],
3265 working_space[peak_vel + 13]);
3268 working_space[2 * shift + k] += chi_opt *
c;
3270 working_space[4 * shift + k] +=
b *
c;
3275 a =
Derampx(i1, working_space[7 * j + 5],
3276 working_space[peak_vel],
3277 working_space[peak_vel + 8],
3278 working_space[peak_vel + 10],
3279 working_space[peak_vel + 12]);
3282 working_space[2 * shift + k] += chi_opt *
c;
3284 working_space[4 * shift + k] +=
b *
c;
3289 a =
Derampx(i2, working_space[7 * j + 6],
3290 working_space[peak_vel + 1],
3291 working_space[peak_vel + 9],
3292 working_space[peak_vel + 11],
3293 working_space[peak_vel + 13]);
3296 working_space[2 * shift + k] += chi_opt *
c;
3298 working_space[4 * shift + k] +=
b *
c;
3303 a =
Deri01(i1, working_space[7 * j + 3],
3304 working_space[7 * j + 5],
3305 working_space[peak_vel],
3306 working_space[peak_vel + 8],
3307 working_space[peak_vel + 10],
3308 working_space[peak_vel + 12]);
3311 working_space[2 * shift + k] += chi_opt *
c;
3313 working_space[4 * shift + k] +=
b *
c;
3318 a =
Deri01(i2, working_space[7 * j + 4],
3319 working_space[7 * j + 6],
3320 working_space[peak_vel + 1],
3321 working_space[peak_vel + 9],
3322 working_space[peak_vel + 11],
3323 working_space[peak_vel + 13]);
3326 working_space[2 * shift + k] += chi_opt *
c;
3328 working_space[4 * shift + k] +=
b *
c;
3335 working_space, working_space[peak_vel],
3336 working_space[peak_vel + 1],
3337 working_space[peak_vel + 2],
3338 working_space[peak_vel + 6],
3339 working_space[peak_vel + 7],
3340 working_space[peak_vel + 8],
3341 working_space[peak_vel + 10],
3342 working_space[peak_vel + 12],
3343 working_space[peak_vel + 13]);
3346 working_space[2 * shift + k] += chi_opt *
c;
3348 working_space[4 * shift + k] +=
b *
c;
3354 working_space, working_space[peak_vel],
3355 working_space[peak_vel + 1],
3356 working_space[peak_vel + 2],
3357 working_space[peak_vel + 6],
3358 working_space[peak_vel + 7],
3359 working_space[peak_vel + 9],
3360 working_space[peak_vel + 11],
3361 working_space[peak_vel + 12],
3362 working_space[peak_vel + 13]);
3365 working_space[2 * shift + k] += chi_opt *
c;
3367 working_space[4 * shift + k] +=
b *
c;
3373 working_space, working_space[peak_vel],
3374 working_space[peak_vel + 1],
3375 working_space[peak_vel + 2]);
3378 working_space[2 * shift + k] += chi_opt *
c;
3380 working_space[4 * shift + k] +=
b *
c;
3388 working_space[2 * shift + k] += chi_opt *
c;
3390 working_space[4 * shift + k] +=
b *
c;
3398 working_space[2 * shift + k] += chi_opt *
c;
3400 working_space[4 * shift + k] +=
b *
c;
3408 working_space[2 * shift + k] += chi_opt *
c;
3410 working_space[4 * shift + k] +=
b *
c;
3416 working_space, working_space[peak_vel],
3417 working_space[peak_vel + 1],
3418 working_space[peak_vel + 12],
3419 working_space[peak_vel + 13]);
3422 working_space[2 * shift + k] += chi_opt *
c;
3424 working_space[4 * shift + k] +=
b *
c;
3430 working_space, working_space[peak_vel],
3431 working_space[peak_vel + 1]);
3434 working_space[2 * shift + k] += chi_opt *
c;
3436 working_space[4 * shift + k] +=
b *
c;
3442 working_space[peak_vel],
3443 working_space[peak_vel + 12]);
3446 working_space[2 * shift + k] += chi_opt *
c;
3448 working_space[4 * shift + k] +=
b *
c;
3454 working_space[peak_vel + 1],
3455 working_space[peak_vel + 13]);
3458 working_space[2 * shift + k] += chi_opt *
c;
3460 working_space[4 * shift + k] +=
b *
c;
3466 working_space[peak_vel]);
3469 working_space[2 * shift + k] += chi_opt *
c;
3471 working_space[4 * shift + k] +=
b *
c;
3477 working_space[peak_vel + 1]);
3480 working_space[2 * shift + k] += chi_opt *
c;
3482 working_space[4 * shift + k] +=
b *
c;
3488 working_space, working_space[peak_vel],
3489 working_space[peak_vel + 1],
3490 working_space[peak_vel + 6],
3491 working_space[peak_vel + 8],
3492 working_space[peak_vel + 12],
3493 working_space[peak_vel + 13]);
3496 working_space[2 * shift + k] += chi_opt *
c;
3498 working_space[4 * shift + k] +=
b *
c;
3504 working_space, working_space[peak_vel],
3505 working_space[peak_vel + 1],
3506 working_space[peak_vel + 6],
3507 working_space[peak_vel + 8],
3508 working_space[peak_vel + 12],
3509 working_space[peak_vel + 13]);
3512 working_space[2 * shift + k] += chi_opt *
c;
3514 working_space[4 * shift + k] +=
b *
c;
3522 chi_er = chi_cel /
b;
3523 for (i = 0, j = 0; i <
fNPeaks; i++) {
3525 Volume(working_space[7 * i], working_space[peak_vel],
3526 working_space[peak_vel + 1], working_space[peak_vel + 2]);
3530 a =
Derpa2(working_space[peak_vel],
3531 working_space[peak_vel + 1],
3532 working_space[peak_vel + 2]);
3533 b = working_space[4 * shift + j];
3543 working_space[peak_vel + 1],
3544 working_space[peak_vel + 2]);
3545 b = working_space[4 * shift + peak_vel];
3555 working_space[peak_vel],
3556 working_space[peak_vel + 2]);
3557 b = working_space[4 * shift + peak_vel + 1];
3566 a =
Derpro(working_space[shift + j], working_space[peak_vel],
3567 working_space[peak_vel + 1],
3568 working_space[peak_vel + 2]);
3569 b = working_space[4 * shift + peak_vel + 2];
3584 fAmpCalc[i] = working_space[shift + j];
3585 if (working_space[3 * shift + j] != 0)
3596 if (working_space[3 * shift + j] != 0)
3607 if (working_space[3 * shift + j] != 0)
3618 if (working_space[3 * shift + j] != 0)
3629 if (working_space[3 * shift + j] != 0)
3640 if (working_space[3 * shift + j] != 0)
3651 if (working_space[3 * shift + j] != 0)
3663 if (working_space[3 * shift + j] != 0)
3674 if (working_space[3 * shift + j] != 0)
3684 fRoCalc = working_space[shift + j];
3685 if (working_space[3 * shift + j] != 0)
3695 fA0Calc = working_space[shift + j];
3696 if (working_space[3 * shift + j] != 0)
3706 fAxCalc = working_space[shift + j];
3707 if (working_space[3 * shift + j] != 0)
3717 fAyCalc = working_space[shift + j];
3718 if (working_space[3 * shift + j] != 0)
3728 fTxyCalc = working_space[shift + j];
3729 if (working_space[3 * shift + j] != 0)
3739 fSxyCalc = working_space[shift + j];
3740 if (working_space[3 * shift + j] != 0)
3750 fTxCalc = working_space[shift + j];
3751 if (working_space[3 * shift + j] != 0)
3761 fTyCalc = working_space[shift + j];
3762 if (working_space[3 * shift + j] != 0)
3772 fSxCalc = working_space[shift + j];
3773 if (working_space[3 * shift + j] != 0)
3783 fSyCalc = working_space[shift + j];
3784 if (working_space[3 * shift + j] != 0)
3794 fBxCalc = working_space[shift + j];
3795 if (working_space[3 * shift + j] != 0)
3805 fByCalc = working_space[shift + j];
3806 if (working_space[3 * shift + j] != 0)
3820 working_space, working_space[peak_vel],
3821 working_space[peak_vel + 1],
3822 working_space[peak_vel + 2],
3823 working_space[peak_vel + 3],
3824 working_space[peak_vel + 4],
3825 working_space[peak_vel + 5],
3826 working_space[peak_vel + 6],
3827 working_space[peak_vel + 7],
3828 working_space[peak_vel + 8],
3829 working_space[peak_vel + 9],
3830 working_space[peak_vel + 10],
3831 working_space[peak_vel + 11],
3832 working_space[peak_vel + 12],
3833 working_space[peak_vel + 13]);
3837 delete [] working_space;
3947 Int_t i, i1, i2, j, k, shift =
3948 7 *
fNPeaks + 14, peak_vel,
size, iter, regul_cycle,
3950 Double_t a,
b,
c, alpha, chi_opt, yw, ywm,
f, chi2, chi_min, chi = 0
3951 , pi, pmin = 0, chi_cel = 0, chi_er;
3953 for (i = 0, j = 0; i <
fNPeaks; i++) {
3954 working_space[7 * i] =
fAmpInit[i];
3956 working_space[shift + j] =
fAmpInit[i];
4001 working_space[7 * i + 2] =
fRoInit;
4003 working_space[shift + j] =
fRoInit;
4006 working_space[7 * i + 3] =
fA0Init;
4008 working_space[shift + j] =
fA0Init;
4011 working_space[7 * i + 4] =
fAxInit;
4013 working_space[shift + j] =
fAxInit;
4016 working_space[7 * i + 5] =
fAyInit;
4018 working_space[shift + j] =
fAyInit;
4021 working_space[7 * i + 6] =
fTxyInit;
4023 working_space[shift + j] =
fTxyInit;
4026 working_space[7 * i + 7] =
fSxyInit;
4028 working_space[shift + j] =
fSxyInit;
4031 working_space[7 * i + 8] =
fTxInit;
4033 working_space[shift + j] =
fTxInit;
4036 working_space[7 * i + 9] =
fTyInit;
4038 working_space[shift + j] =
fTyInit;
4041 working_space[7 * i + 10] =
fSxyInit;
4043 working_space[shift + j] =
fSxInit;
4046 working_space[7 * i + 11] =
fSyInit;
4048 working_space[shift + j] =
fSyInit;
4051 working_space[7 * i + 12] =
fBxInit;
4053 working_space[shift + j] =
fBxInit;
4056 working_space[7 * i + 13] =
fByInit;
4058 working_space[shift + j] =
fByInit;
4063 for (i = 0; i <
size; i++)
4066 for (j = 0; j <
size; j++) {
4067 working_space[3 * shift + j] = 0;
4068 for (k = 0; k < (
size + 4); k++) {
4069 working_matrix[j][k] = 0;
4079 for (j = 0, k = 0; j <
fNPeaks; j++) {
4081 working_space[2 * shift + k] =
4083 working_space[7 * j + 1],
4084 working_space[7 * j + 2],
4085 working_space[peak_vel],
4086 working_space[peak_vel + 1],
4087 working_space[peak_vel + 2],
4088 working_space[peak_vel + 6],
4089 working_space[peak_vel + 7],
4090 working_space[peak_vel + 12],
4091 working_space[peak_vel + 13]);
4095 working_space[2 * shift + k] =
4097 working_space[7 * j],
4098 working_space[7 * j + 1],
4099 working_space[7 * j + 2],
4100 working_space[peak_vel],
4101 working_space[peak_vel + 1],
4102 working_space[peak_vel + 2],
4103 working_space[peak_vel + 6],
4104 working_space[peak_vel + 7],
4105 working_space[peak_vel + 12],
4106 working_space[peak_vel + 13]);
4110 working_space[2 * shift + k] =
4112 working_space[7 * j],
4113 working_space[7 * j + 1],
4114 working_space[7 * j + 2],
4115 working_space[peak_vel],
4116 working_space[peak_vel + 1],
4117 working_space[peak_vel + 2],
4118 working_space[peak_vel + 6],
4119 working_space[peak_vel + 7],
4120 working_space[peak_vel + 12],
4121 working_space[peak_vel + 13]);
4125 working_space[2 * shift + k] =
4126 Derampx(i1, working_space[7 * j + 5],
4127 working_space[peak_vel],
4128 working_space[peak_vel + 8],
4129 working_space[peak_vel + 10],
4130 working_space[peak_vel + 12]);
4134 working_space[2 * shift + k] =
4135 Derampx(i2, working_space[7 * j + 6],
4136 working_space[peak_vel + 1],
4137 working_space[peak_vel + 9],
4138 working_space[peak_vel + 11],
4139 working_space[peak_vel + 13]);
4143 working_space[2 * shift + k] =
4144 Deri01(i1, working_space[7 * j + 3],
4145 working_space[7 * j + 5],
4146 working_space[peak_vel],
4147 working_space[peak_vel + 8],
4148 working_space[peak_vel + 10],
4149 working_space[peak_vel + 12]);
4153 working_space[2 * shift + k] =
4154 Deri01(i2, working_space[7 * j + 4],
4155 working_space[7 * j + 6],
4156 working_space[peak_vel + 1],
4157 working_space[peak_vel + 9],
4158 working_space[peak_vel + 11],
4159 working_space[peak_vel + 13]);
4163 working_space[2 * shift + k] =
4165 working_space, working_space[peak_vel],
4166 working_space[peak_vel + 1],
4167 working_space[peak_vel + 2],
4168 working_space[peak_vel + 6],
4169 working_space[peak_vel + 7],
4170 working_space[peak_vel + 8],
4171 working_space[peak_vel + 10],
4172 working_space[peak_vel + 12],
4173 working_space[peak_vel + 13]);
4177 working_space[2 * shift + k] =
4179 working_space, working_space[peak_vel],
4180 working_space[peak_vel + 1],
4181 working_space[peak_vel + 2],
4182 working_space[peak_vel + 6],
4183 working_space[peak_vel + 7],
4184 working_space[peak_vel + 9],
4185 working_space[peak_vel + 11],
4186 working_space[peak_vel + 12],
4187 working_space[peak_vel + 13]);
4191 working_space[2 * shift + k] =
4193 working_space, working_space[peak_vel],
4194 working_space[peak_vel + 1],
4195 working_space[peak_vel + 2]);
4199 working_space[2 * shift + k] = 1.;
4203 working_space[2 * shift + k] = i1;
4207 working_space[2 * shift + k] = i2;
4211 working_space[2 * shift + k] =
4213 working_space, working_space[peak_vel],
4214 working_space[peak_vel + 1],
4215 working_space[peak_vel + 12],
4216 working_space[peak_vel + 13]);
4220 working_space[2 * shift + k] =
4222 working_space, working_space[peak_vel],
4223 working_space[peak_vel + 1]);
4227 working_space[2 * shift + k] =
4229 working_space[peak_vel],
4230 working_space[peak_vel + 12]);
4234 working_space[2 * shift + k] =
4236 working_space[peak_vel + 1],
4237 working_space[peak_vel + 13]);
4241 working_space[2 * shift + k] =
4243 working_space[peak_vel]);
4247 working_space[2 * shift + k] =
4249 working_space[peak_vel + 1]);
4253 working_space[2 * shift + k] =
4255 working_space, working_space[peak_vel],
4256 working_space[peak_vel + 1],
4257 working_space[peak_vel + 6],
4258 working_space[peak_vel + 8],
4259 working_space[peak_vel + 12],
4260 working_space[peak_vel + 13]);
4264 working_space[2 * shift + k] =
4266 working_space, working_space[peak_vel],
4267 working_space[peak_vel + 1],
4268 working_space[peak_vel + 6],
4269 working_space[peak_vel + 8],
4270 working_space[peak_vel + 12],
4271 working_space[peak_vel + 13]);
4274 yw = source[i1][i2];
4277 working_space, working_space[peak_vel],
4278 working_space[peak_vel + 1],
4279 working_space[peak_vel + 2],
4280 working_space[peak_vel + 3],
4281 working_space[peak_vel + 4],
4282 working_space[peak_vel + 5],
4283 working_space[peak_vel + 6],
4284 working_space[peak_vel + 7],
4285 working_space[peak_vel + 8],
4286 working_space[peak_vel + 9],
4287 working_space[peak_vel + 10],
4288 working_space[peak_vel + 11],
4289 working_space[peak_vel + 12],
4290 working_space[peak_vel + 13]);
4298 chi_opt += (yw -
f) * (yw -
f) / ywm;
4316 for (j = 0; j <
size; j++) {
4317 for (k = 0; k <
size; k++) {
4318 b = working_space[2 * shift +
4319 j] * working_space[2 * shift +
4322 b =
b * (4 * yw - 2 *
f) / ywm;
4323 working_matrix[j][k] +=
b;
4325 working_space[3 * shift + j] +=
b;
4329 b = (
f *
f - yw * yw) / (ywm * ywm);
4333 for (j = 0; j <
size; j++) {
4334 working_matrix[j][
size] -=
4335 b * working_space[2 * shift + j];
4339 for (i = 0; i <
size; i++) {
4340 working_matrix[i][
size + 1] = 0;
4343 for (i = 0; i <
size; i++) {
4344 working_space[2 * shift + i] = working_matrix[i][
size + 1];
4353 for (j = 0; j <
size; j++) {
4354 working_space[4 * shift + j] = working_space[shift + j];
4360 chi_min = 10000 * chi2;
4363 chi_min = 0.1 * chi2;
4365 for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
4366 for (j = 0; j <
size; j++) {
4367 working_space[shift + j] = working_space[4 * shift + j] + pi * alpha * working_space[2 * shift + j];
4369 for (i = 0, j = 0; i <
fNPeaks; i++) {
4371 if (working_space[shift + j] < 0)
4372 working_space[shift + j] = 0;
4373 working_space[7 * i] = working_space[shift + j];
4377 if (working_space[shift + j] <
fXmin)
4378 working_space[shift + j] =
fXmin;
4379 if (working_space[shift + j] >
fXmax)
4380 working_space[shift + j] =
fXmax;
4381 working_space[7 * i + 1] = working_space[shift + j];
4385 if (working_space[shift + j] <
fYmin)
4386 working_space[shift + j] =
fYmin;
4387 if (working_space[shift + j] >
fYmax)
4388 working_space[shift + j] =
fYmax;
4389 working_space[7 * i + 2] = working_space[shift + j];
4393 if (working_space[shift + j] < 0)
4394 working_space[shift + j] = 0;
4395 working_space[7 * i + 3] = working_space[shift + j];
4399 if (working_space[shift + j] < 0)
4400 working_space[shift + j] = 0;
4401 working_space[7 * i + 4] = working_space[shift + j];
4405 if (working_space[shift + j] <
fXmin)
4406 working_space[shift + j] =
fXmin;
4407 if (working_space[shift + j] >
fXmax)
4408 working_space[shift + j] =
fXmax;
4409 working_space[7 * i + 5] = working_space[shift + j];
4413 if (working_space[shift + j] <
fYmin)
4414 working_space[shift + j] =
fYmin;
4415 if (working_space[shift + j] >
fYmax)
4416 working_space[shift + j] =
fYmax;
4417 working_space[7 * i + 6] = working_space[shift + j];
4422 if (working_space[shift + j] < 0.001) {
4423 working_space[shift + j] = 0.001;
4425 working_space[peak_vel] = working_space[shift + j];
4429 if (working_space[shift + j] < 0.001) {
4430 working_space[shift + j] = 0.001;
4432 working_space[peak_vel + 1] = working_space[shift + j];
4436 if (working_space[shift + j] < -1) {
4437 working_space[shift + j] = -1;
4439 if (working_space[shift + j] > 1) {
4440 working_space[shift + j] = 1;
4442 working_space[peak_vel + 2] = working_space[shift + j];
4446 working_space[peak_vel + 3] = working_space[shift + j];
4450 working_space[peak_vel + 4] = working_space[shift + j];
4454 working_space[peak_vel + 5] = working_space[shift + j];
4458 working_space[peak_vel + 6] = working_space[shift + j];
4462 working_space[peak_vel + 7] = working_space[shift + j];
4466 working_space[peak_vel + 8] = working_space[shift + j];
4470 working_space[peak_vel + 9] = working_space[shift + j];
4474 working_space[peak_vel + 10] = working_space[shift + j];
4478 working_space[peak_vel + 11] = working_space[shift + j];
4482 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
4483 if (working_space[shift + j] < 0)
4484 working_space[shift + j] = -0.001;
4486 working_space[shift + j] = 0.001;
4488 working_space[peak_vel + 12] = working_space[shift + j];
4492 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
4493 if (working_space[shift + j] < 0)
4494 working_space[shift + j] = -0.001;
4496 working_space[shift + j] = 0.001;
4498 working_space[peak_vel + 13] = working_space[shift + j];
4504 yw = source[i1][i2];
4508 working_space[peak_vel],
4509 working_space[peak_vel + 1],
4510 working_space[peak_vel + 2],
4511 working_space[peak_vel + 3],
4512 working_space[peak_vel + 4],
4513 working_space[peak_vel + 5],
4514 working_space[peak_vel + 6],
4515 working_space[peak_vel + 7],
4516 working_space[peak_vel + 8],
4517 working_space[peak_vel + 9],
4518 working_space[peak_vel + 10],
4519 working_space[peak_vel + 11],
4520 working_space[peak_vel + 12],
4521 working_space[peak_vel + 13]);
4534 chi2 += (yw -
f) * (yw -
f) / ywm;
4542 pmin = pi, chi_min = chi2;
4552 for (j = 0; j <
size; j++) {
4553 working_space[shift + j] = working_space[4 * shift + j] + pmin * alpha * working_space[2 * shift + j];
4555 for (i = 0, j = 0; i <
fNPeaks; i++) {
4557 if (working_space[shift + j] < 0)
4558 working_space[shift + j] = 0;
4559 working_space[7 * i] = working_space[shift + j];
4563 if (working_space[shift + j] <
fXmin)
4564 working_space[shift + j] =
fXmin;
4565 if (working_space[shift + j] >
fXmax)
4566 working_space[shift + j] =
fXmax;
4567 working_space[7 * i + 1] = working_space[shift + j];
4571 if (working_space[shift + j] <
fYmin)
4572 working_space[shift + j] =
fYmin;
4573 if (working_space[shift + j] >
fYmax)
4574 working_space[shift + j] =
fYmax;
4575 working_space[7 * i + 2] = working_space[shift + j];
4579 if (working_space[shift + j] < 0)
4580 working_space[shift + j] = 0;
4581 working_space[7 * i + 3] = working_space[shift + j];
4585 if (working_space[shift + j] < 0)
4586 working_space[shift + j] = 0;
4587 working_space[7 * i + 4] = working_space[shift + j];
4591 if (working_space[shift + j] <
fXmin)
4592 working_space[shift + j] =
fXmin;
4593 if (working_space[shift + j] >
fXmax)
4594 working_space[shift + j] =
fXmax;
4595 working_space[7 * i + 5] = working_space[shift + j];
4599 if (working_space[shift + j] <
fYmin)
4600 working_space[shift + j] =
fYmin;
4601 if (working_space[shift + j] >
fYmax)
4602 working_space[shift + j] =
fYmax;
4603 working_space[7 * i + 6] = working_space[shift + j];
4608 if (working_space[shift + j] < 0.001) {
4609 working_space[shift + j] = 0.001;
4611 working_space[peak_vel] = working_space[shift + j];
4615 if (working_space[shift + j] < 0.001) {
4616 working_space[shift + j] = 0.001;
4618 working_space[peak_vel + 1] = working_space[shift + j];
4622 if (working_space[shift + j] < -1) {
4623 working_space[shift + j] = -1;
4625 if (working_space[shift + j] > 1) {
4626 working_space[shift + j] = 1;
4628 working_space[peak_vel + 2] = working_space[shift + j];
4632 working_space[peak_vel + 3] = working_space[shift + j];
4636 working_space[peak_vel + 4] = working_space[shift + j];
4640 working_space[peak_vel + 5] = working_space[shift + j];
4644 working_space[peak_vel + 6] = working_space[shift + j];
4648 working_space[peak_vel + 7] = working_space[shift + j];
4652 working_space[peak_vel + 8] = working_space[shift + j];
4656 working_space[peak_vel + 9] = working_space[shift + j];
4660 working_space[peak_vel + 10] = working_space[shift + j];
4664 working_space[peak_vel + 11] = working_space[shift + j];
4668 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
4669 if (working_space[shift + j] < 0)
4670 working_space[shift + j] = -0.001;
4672 working_space[shift + j] = 0.001;
4674 working_space[peak_vel + 12] = working_space[shift + j];
4678 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
4679 if (working_space[shift + j] < 0)
4680 working_space[shift + j] = -0.001;
4682 working_space[shift + j] = 0.001;
4684 working_space[peak_vel + 13] = working_space[shift + j];
4692 for (j = 0; j <
size; j++) {
4693 working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j];
4695 for (i = 0, j = 0; i <
fNPeaks; i++) {
4697 if (working_space[shift + j] < 0)
4698 working_space[shift + j] = 0;
4699 working_space[7 * i] = working_space[shift + j];
4703 if (working_space[shift + j] <
fXmin)
4704 working_space[shift + j] =
fXmin;
4705 if (working_space[shift + j] >
fXmax)
4706 working_space[shift + j] =
fXmax;
4707 working_space[7 * i + 1] = working_space[shift + j];
4711 if (working_space[shift + j] <
fYmin)
4712 working_space[shift + j] =
fYmin;
4713 if (working_space[shift + j] >
fYmax)
4714 working_space[shift + j] =
fYmax;
4715 working_space[7 * i + 2] = working_space[shift + j];
4719 if (working_space[shift + j] < 0)
4720 working_space[shift + j] = 0;
4721 working_space[7 * i + 3] = working_space[shift + j];
4725 if (working_space[shift + j] < 0)
4726 working_space[shift + j] = 0;
4727 working_space[7 * i + 4] = working_space[shift + j];
4731 if (working_space[shift + j] <
fXmin)
4732 working_space[shift + j] =
fXmin;
4733 if (working_space[shift + j] >
fXmax)
4734 working_space[shift + j] =
fXmax;
4735 working_space[7 * i + 5] = working_space[shift + j];
4739 if (working_space[shift + j] <
fYmin)
4740 working_space[shift + j] =
fYmin;
4741 if (working_space[shift + j] >
fYmax)
4742 working_space[shift + j] =
fYmax;
4743 working_space[7 * i + 6] = working_space[shift + j];
4748 if (working_space[shift + j] < 0.001) {
4749 working_space[shift + j] = 0.001;
4751 working_space[peak_vel] = working_space[shift + j];
4755 if (working_space[shift + j] < 0.001) {
4756 working_space[shift + j] = 0.001;
4758 working_space[peak_vel + 1] = working_space[shift + j];
4762 if (working_space[shift + j] < -1) {
4763 working_space[shift + j] = -1;
4765 if (working_space[shift + j] > 1) {
4766 working_space[shift + j] = 1;
4768 working_space[peak_vel + 2] = working_space[shift + j];
4772 working_space[peak_vel + 3] = working_space[shift + j];
4776 working_space[peak_vel + 4] = working_space[shift + j];
4780 working_space[peak_vel + 5] = working_space[shift + j];
4784 working_space[peak_vel + 6] = working_space[shift + j];
4788 working_space[peak_vel + 7] = working_space[shift + j];
4792 working_space[peak_vel + 8] = working_space[shift + j];
4796 working_space[peak_vel + 9] = working_space[shift + j];
4800 working_space[peak_vel + 10] = working_space[shift + j];
4804 working_space[peak_vel + 11] = working_space[shift + j];
4808 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
4809 if (working_space[shift + j] < 0)
4810 working_space[shift + j] = -0.001;
4812 working_space[shift + j] = 0.001;
4814 working_space[peak_vel + 12] = working_space[shift + j];
4818 if (
TMath::Abs(working_space[shift + j]) < 0.001) {
4819 if (working_space[shift + j] < 0)
4820 working_space[shift + j] = -0.001;
4822 working_space[shift + j] = 0.001;
4824 working_space[peak_vel + 13] = working_space[shift + j];
4830 yw = source[i1][i2];
4833 working_space, working_space[peak_vel],
4834 working_space[peak_vel + 1],
4835 working_space[peak_vel + 2],
4836 working_space[peak_vel + 3],
4837 working_space[peak_vel + 4],
4838 working_space[peak_vel + 5],
4839 working_space[peak_vel + 6],
4840 working_space[peak_vel + 7],
4841 working_space[peak_vel + 8],
4842 working_space[peak_vel + 9],
4843 working_space[peak_vel + 10],
4844 working_space[peak_vel + 11],
4845 working_space[peak_vel + 12],
4846 working_space[peak_vel + 13]);
4859 chi += (yw -
f) * (yw -
f) / ywm;
4867 alpha = alpha * chi_opt / (2 * chi);
4870 alpha = alpha / 10.0;
4873 }
while (((chi > chi_opt
4878 for (j = 0; j <
size; j++) {
4879 working_space[4 * shift + j] = 0;
4880 working_space[2 * shift + j] = 0;
4882 for (i1 =
fXmin, chi_cel = 0; i1 <=
fXmax; i1++) {
4884 yw = source[i1][i2];
4888 working_space, working_space[peak_vel],
4889 working_space[peak_vel + 1],
4890 working_space[peak_vel + 2],
4891 working_space[peak_vel + 3],
4892 working_space[peak_vel + 4],
4893 working_space[peak_vel + 5],
4894 working_space[peak_vel + 6],
4895 working_space[peak_vel + 7],
4896 working_space[peak_vel + 8],
4897 working_space[peak_vel + 9],
4898 working_space[peak_vel + 10],
4899 working_space[peak_vel + 11],
4900 working_space[peak_vel + 12],
4901 working_space[peak_vel + 13]);
4902 chi_opt = (yw -
f) * (yw -
f) / yw;
4903 chi_cel += (yw -
f) * (yw -
f) / yw;
4906 for (j = 0, k = 0; j <
fNPeaks; j++) {
4909 working_space[7 * j + 1],
4910 working_space[7 * j + 2],
4911 working_space[peak_vel],
4912 working_space[peak_vel + 1],
4913 working_space[peak_vel + 2],
4914 working_space[peak_vel + 6],
4915 working_space[peak_vel + 7],
4916 working_space[peak_vel + 12],
4917 working_space[peak_vel + 13]);
4919 working_space[2 * shift + k] += chi_opt;
4921 working_space[4 * shift + k] +=
b;
4927 working_space[7 * j],
4928 working_space[7 * j + 1],
4929 working_space[7 * j + 2],
4930 working_space[peak_vel],
4931 working_space[peak_vel + 1],
4932 working_space[peak_vel + 2],
4933 working_space[peak_vel + 6],
4934 working_space[peak_vel + 7],
4935 working_space[peak_vel + 12],
4936 working_space[peak_vel + 13]);
4938 working_space[2 * shift + k] += chi_opt;
4940 working_space[4 * shift + k] +=
b;
4946 working_space[7 * j],
4947 working_space[7 * j + 1],
4948 working_space[7 * j + 2],
4949 working_space[peak_vel],
4950 working_space[peak_vel + 1],
4951 working_space[peak_vel + 2],
4952 working_space[peak_vel + 6],
4953 working_space[peak_vel + 7],
4954 working_space[peak_vel + 12],
4955 working_space[peak_vel + 13]);
4957 working_space[2 * shift + k] += chi_opt;
4959 working_space[4 * shift + k] +=
b;
4964 a =
Derampx(i1, working_space[7 * j + 5],
4965 working_space[peak_vel],
4966 working_space[peak_vel + 8],
4967 working_space[peak_vel + 10],
4968 working_space[peak_vel + 12]);
4970 working_space[2 * shift + k] += chi_opt;
4972 working_space[4 * shift + k] +=
b;
4977 a =
Derampx(i2, working_space[7 * j + 6],
4978 working_space[peak_vel + 1],
4979 working_space[peak_vel + 9],
4980 working_space[peak_vel + 11],
4981 working_space[peak_vel + 13]);
4983 working_space[2 * shift + k] += chi_opt;
4985 working_space[4 * shift + k] +=
b;
4990 a =
Deri01(i1, working_space[7 * j + 3],
4991 working_space[7 * j + 5],
4992 working_space[peak_vel],
4993 working_space[peak_vel + 8],
4994 working_space[peak_vel + 10],
4995 working_space[peak_vel + 12]);
4997 working_space[2 * shift + k] += chi_opt;
4999 working_space[4 * shift + k] +=
b;
5004 a =
Deri01(i2, working_space[7 * j + 4],
5005 working_space[7 * j + 6],
5006 working_space[peak_vel + 1],
5007 working_space[peak_vel + 9],
5008 working_space[peak_vel + 11],
5009 working_space[peak_vel + 13]);
5011 working_space[2 * shift + k] += chi_opt;
5013 working_space[4 * shift + k] +=
b;
5020 working_space, working_space[peak_vel],
5021 working_space[peak_vel + 1],
5022 working_space[peak_vel + 2],
5023 working_space[peak_vel + 6],
5024 working_space[peak_vel + 7],
5025 working_space[peak_vel + 8],
5026 working_space[peak_vel + 10],
5027 working_space[peak_vel + 12],
5028 working_space[peak_vel + 13]);
5030 working_space[2 * shift + k] += chi_opt;
5032 working_space[4 * shift + k] +=
b;
5038 working_space, working_space[peak_vel],
5039 working_space[peak_vel + 1],
5040 working_space[peak_vel + 2],
5041 working_space[peak_vel + 6],
5042 working_space[peak_vel + 7],
5043 working_space[peak_vel + 9],
5044 working_space[peak_vel + 11],
5045 working_space[peak_vel + 12],
5046 working_space[peak_vel + 13]);
5048 working_space[2 * shift + k] += chi_opt;
5050 working_space[4 * shift + k] +=
b;
5056 working_space, working_space[peak_vel],
5057 working_space[peak_vel + 1],
5058 working_space[peak_vel + 2]);
5060 working_space[2 * shift + k] += chi_opt;
5062 working_space[4 * shift + k] +=
b;
5069 working_space[2 * shift + k] += chi_opt;
5071 working_space[4 * shift + k] +=
b;
5078 working_space[2 * shift + k] += chi_opt;
5080 working_space[4 * shift + k] +=
b;
5087 working_space[2 * shift + k] += chi_opt;
5089 working_space[4 * shift + k] +=
b;
5095 working_space, working_space[peak_vel],
5096 working_space[peak_vel + 1],
5097 working_space[peak_vel + 12],
5098 working_space[peak_vel + 13]);
5100 working_space[2 * shift + k] += chi_opt;
5102 working_space[4 * shift + k] +=
b;
5108 working_space, working_space[peak_vel],
5109 working_space[peak_vel + 1]);
5111 working_space[2 * shift + k] += chi_opt;
5113 working_space[4 * shift + k] +=
b;
5119 working_space[peak_vel],
5120 working_space[peak_vel + 12]);
5122 working_space[2 * shift + k] += chi_opt;
5124 working_space[4 * shift + k] +=
b;
5130 working_space[peak_vel + 1],
5131 working_space[peak_vel + 13]);
5133 working_space[2 * shift + k] += chi_opt;
5135 working_space[4 * shift + k] +=
b;
5141 working_space[peak_vel]);
5143 working_space[2 * shift + k] += chi_opt;
5145 working_space[4 * shift + k] +=
b;
5151 working_space[peak_vel + 1]);
5153 working_space[2 * shift + k] += chi_opt;
5155 working_space[4 * shift + k] +=
b;
5161 working_space, working_space[peak_vel],
5162 working_space[peak_vel + 1],
5163 working_space[peak_vel + 6],
5164 working_space[peak_vel + 8],
5165 working_space[peak_vel + 12],
5166 working_space[peak_vel + 13]);
5168 working_space[2 * shift + k] += chi_opt;
5170 working_space[4 * shift + k] +=
b;
5176 working_space, working_space[peak_vel],
5177 working_space[peak_vel + 1],
5178 working_space[peak_vel + 6],
5179 working_space[peak_vel + 8],
5180 working_space[peak_vel + 12],
5181 working_space[peak_vel + 13]);
5183 working_space[2 * shift + k] += chi_opt;
5185 working_space[4 * shift + k] +=
b;
5193 chi_er = chi_cel /
b;
5194 for (i = 0, j = 0; i <
fNPeaks; i++) {
5196 Volume(working_space[7 * i], working_space[peak_vel],
5197 working_space[peak_vel + 1], working_space[peak_vel + 2]);
5201 a =
Derpa2(working_space[peak_vel],
5202 working_space[peak_vel + 1],
5203 working_space[peak_vel + 2]);
5204 b = working_space[4 * shift + j];
5214 working_space[peak_vel + 1],
5215 working_space[peak_vel + 2]);
5216 b = working_space[4 * shift + peak_vel];
5226 working_space[peak_vel],
5227 working_space[peak_vel + 2]);
5228 b = working_space[4 * shift + peak_vel + 1];
5237 a =
Derpro(working_space[shift + j], working_space[peak_vel],
5238 working_space[peak_vel + 1],
5239 working_space[peak_vel + 2]);
5240 b = working_space[4 * shift + peak_vel + 2];
5255 fAmpCalc[i] = working_space[shift + j];
5256 if (working_space[3 * shift + j] != 0)
5267 if (working_space[3 * shift + j] != 0)
5278 if (working_space[3 * shift + j] != 0)
5289 if (working_space[3 * shift + j] != 0)
5300 if (working_space[3 * shift + j] != 0)
5311 if (working_space[3 * shift + j] != 0)
5322 if (working_space[3 * shift + j] != 0)
5334 if (working_space[3 * shift + j] != 0)
5345 if (working_space[3 * shift + j] != 0)
5355 fRoCalc = working_space[shift + j];
5356 if (working_space[3 * shift + j] != 0)
5366 fA0Calc = working_space[shift + j];
5367 if (working_space[3 * shift + j] != 0)
5377 fAxCalc = working_space[shift + j];
5378 if (working_space[3 * shift + j] != 0)
5388 fAyCalc = working_space[shift + j];
5389 if (working_space[3 * shift + j] != 0)
5399 fTxyCalc = working_space[shift + j];
5400 if (working_space[3 * shift + j] != 0)
5410 fSxyCalc = working_space[shift + j];
5411 if (working_space[3 * shift + j] != 0)
5421 fTxCalc = working_space[shift + j];
5422 if (working_space[3 * shift + j] != 0)
5432 fTyCalc = working_space[shift + j];
5433 if (working_space[3 * shift + j] != 0)
5443 fSxCalc = working_space[shift + j];
5444 if (working_space[3 * shift + j] != 0)
5454 fSyCalc = working_space[shift + j];
5455 if (working_space[3 * shift + j] != 0)
5465 fBxCalc = working_space[shift + j];
5466 if (working_space[3 * shift + j] != 0)
5476 fByCalc = working_space[shift + j];
5477 if (working_space[3 * shift + j] != 0)
5491 working_space, working_space[peak_vel],
5492 working_space[peak_vel + 1],
5493 working_space[peak_vel + 2],
5494 working_space[peak_vel + 3],
5495 working_space[peak_vel + 4],
5496 working_space[peak_vel + 5],
5497 working_space[peak_vel + 6],
5498 working_space[peak_vel + 7],
5499 working_space[peak_vel + 8],
5500 working_space[peak_vel + 9],
5501 working_space[peak_vel + 10],
5502 working_space[peak_vel + 11],
5503 working_space[peak_vel + 12],
5504 working_space[peak_vel + 13]);
5509 for (i = 0; i <
size; i++)
delete [] working_matrix[i];
5510 delete [] working_matrix;
5511 delete [] working_space;
5528 Error(
"SetFitParameters",
"Wrong range");
5531 if (numberIterations <= 0){
5532 Error(
"SetFitParameters",
"Invalid number of iterations, must be positive");
5535 if (alpha <= 0 || alpha > 1){
5536 Error (
"SetFitParameters",
"Invalid step coefficient alpha, must be > than 0 and <=1");
5542 Error(
"SetFitParameters",
"Wrong type of statistic");
5547 Error(
"SetFitParameters",
"Wrong optimization algorithm");
5553 Error(
"SetFitParameters",
"Wrong power");
5558 Error(
"SetFitParameters",
"Wrong order of Taylor development");
5583void 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)
5585 if (sigmaX <= 0 || sigmaY <= 0){
5586 Error (
"SetPeakParameters",
"Invalid sigma, must be > than 0");
5589 if (ro < -1 || ro > 1){
5590 Error (
"SetPeakParameters",
"Invalid ro, must be from region <-1,1>");
5595 if(positionInitX[i] <
fXmin || positionInitX[i] >
fXmax){
5596 Error (
"SetPeakParameters",
"Invalid peak position, must be in the range fXmin, fXmax");
5599 if(positionInitY[i] <
fYmin || positionInitY[i] >
fYmax){
5600 Error (
"SetPeakParameters",
"Invalid peak position, must be in the range fYmin, fYmax");
5603 if(positionInitX1[i] <
fXmin || positionInitX1[i] >
fXmax){
5604 Error (
"SetPeakParameters",
"Invalid ridge position, must be in the range fXmin, fXmax");
5607 if(positionInitY1[i] <
fYmin || positionInitY1[i] >
fYmax){
5608 Error (
"SetPeakParameters",
"Invalid ridge position, must be in the range fYmin, fYmax");
5612 Error (
"SetPeakParameters",
"Invalid peak amplitude, must be > than 0");
5615 if(ampInitX1[i] < 0){
5616 Error (
"SetPeakParameters",
"Invalid x ridge amplitude, must be > than 0");
5619 if(ampInitY1[i] < 0){
5620 Error (
"SetPeakParameters",
"Invalid y ridge amplitude, must be > than 0");
5681void 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)
5759 amplitudeErrors[i] =
fAmpErr[i];
5858void 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)
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
The TNamed class is the base class for all named ROOT classes.
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Advanced 2-dimensional spectra fitting functions.
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 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...
Bool_t fFixBx
logical value of b parameter for 1D ridges in x direction, which allows to fix the parameter (not to ...
Bool_t * fFixPositionY
[fNPeaks] array of logical values which allow to fix appropriate y positions of 2D peaks (not fit)....
Double_t fSigmaErrY
error value of sigma y parameter
Bool_t * fFixAmpY1
[fNPeaks] array of logical values which allow to fix appropriate amplitudes of 1D ridges in y directi...
Double_t * fAmpCalcX1
[fNPeaks] array of calculated values of amplitudes of 1D ridges in x direction, output parameters
Double_t Volume(Double_t a, Double_t sx, Double_t sy, Double_t ro)
This function calculates volume of a peak.
Bool_t fFixSxy
logical value of s parameter for 2D peaks, which allows to fix the parameter (not to fit).
Double_t fTxErr
error value of t parameter for 1D ridges in x direction
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...
Double_t * fAmpErrY1
[fNPeaks] array of amplitudes errors of 1D ridges in y direction, output parameters
Double_t fBxInit
initial value of b parameter for 1D ridges in x direction (slope), for details see html manual and re...
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 fA0Err
error value of background a0 parameter
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 fSigmaCalcX
calculated value of sigma x parameter
Bool_t * fFixAmpX1
[fNPeaks] array of logical values which allow to fix appropriate amplitudes of 1D ridges in x directi...
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...
Bool_t fFixRo
logical value of correlation coefficient, which allows to fix the parameter (not to fit).
Double_t * fAmpErr
[fNPeaks] array of amplitudes errors of 2D peaks, output parameters
Bool_t fFixSx
logical value of s parameter for 1D ridges in x direction, which allows to fix the parameter (not to ...
Double_t * fAmpInit
[fNPeaks] array of initial values of amplitudes of 2D peaks, input parameters
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.
Double_t * fVolume
[fNPeaks] array of calculated volumes of 2D peaks, output parameters
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 * fPositionCalcY
[fNPeaks] array of calculated values of y positions of 2D peaks, output parameters
Int_t fFitTaylor
order of Taylor expansion, possible values kFitTaylorOrderFirst, kFitTaylorOrderSecond....
Double_t fSxyErr
error value of s parameter for 2D peaks
Double_t * fPositionErrY1
[fNPeaks] array of y positions errors of 1D ridges, output parameters
Double_t fBxCalc
calculated value of b parameter for 1D ridges in x direction
Double_t fRoErr
error value of correlation coefficient
Double_t fTxCalc
calculated value of t parameter for 1D ridges in x direction
Double_t fRoCalc
calculated value of correlation coefficient
Bool_t fFixBy
logical value of b parameter for 1D ridges in y direction, which allows to fix the parameter (not to ...
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 fTxyInit
initial value of t parameter for 2D peaks (relative amplitude of tail), for details see html manual a...
Bool_t fFixSigmaY
logical value of sigma y parameter, which allows to fix the parameter (not to fit).
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...
void GetRo(Double_t &ro, Double_t &roErr)
This function gets the ro parameter and its error.
Int_t fXmax
last fitted channel in x direction
Double_t fChi
here the fitting functions return resulting chi square
void GetVolumes(Double_t *volumes)
This function gets the volumes of fitted 2D peaks.
Double_t fByErr
error value of b parameter for 1D ridges in y direction
Double_t fA0Init
initial value of background a0 parameter(backgroud is estimated as a0+ax*x+ay*y)
Double_t fRoInit
initial value of correlation coefficient
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 fSigmaInitX
initial value of sigma x parameter
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 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...
Double_t * fPositionInitY1
[fNPeaks] array of initial y positions of 1D ridges, input parameters
Bool_t fFixSy
logical value of s parameter for 1D ridges in y direction, which allows to fix the parameter (not to ...
Bool_t * fFixAmp
[fNPeaks] array of logical values which allow to fix appropriate amplitudes of 2D peaks (not fit)....
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 fAyErr
error value of background ay parameter
Double_t fSxyCalc
calculated value of s parameter for 2D peaks
Double_t * fPositionInitX1
[fNPeaks] array of initial x positions of 1D ridges, input parameters
Double_t fTyCalc
calculated value of t parameter for 1D ridges in y direction
Double_t * fPositionErrX1
[fNPeaks] array of x positions errors of 1D ridges, output parameters
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 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 * fPositionInitY
[fNPeaks] array of initial values of y positions of 2D peaks, input parameters
Double_t fSxCalc
calculated value of s parameter for 1D ridges in x direction
void GetAmplitudeErrors(Double_t *amplitudeErrors, Double_t *amplitudeErrorsX1, Double_t *amplitudeErrorsY1)
This function gets the amplitudes of fitted 2D peaks and 1D ridges.
Double_t fAxCalc
calculated value of background ax parameter
Double_t fSyErr
error value of s parameter for 1D ridges in y direction
Int_t fStatisticType
type of statistics, possible values kFitOptimChiCounts (chi square statistics with counts as weightin...
Int_t fYmax
last fitted channel in y direction
Double_t * fPositionErrY
[fNPeaks] array of error values of y positions of 2D peaks, output parameters
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...
Double_t fByCalc
calculated value of b parameter for 1D ridges in y direction
Double_t fTyInit
initial value of t parameter for 1D ridges in y direction (relative amplitude of tail),...
Double_t fTxyCalc
calculated value of t parameter for 2D peaks
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 fByInit
initial value of b parameter for 1D ridges in y direction (slope), for details see html manual and re...
Double_t fSigmaCalcY
calculated value of sigma y parameter
Double_t * fAmpErrX1
[fNPeaks] array of amplitudes errors of 1D ridges in x direction, output parameters
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...
Int_t fNumberIterations
number of iterations in fitting procedure, input parameter, it should be > 0
Bool_t fFixTx
logical value of t parameter for 1D ridges in x direction, which allows to fix the parameter (not to ...
Double_t * fPositionCalcY1
[fNPeaks] array of calculated y positions of 1D ridges, output parameters
Bool_t fFixTxy
logical value of t parameter for 2D peaks, which allows to fix the parameter (not to fit).
void FitAwmi(Double_t **source)
This function fits the source spectrum.
Bool_t fFixTy
logical value of t parameter for 1D ridges in y direction, which allows to fix the parameter (not to ...
Double_t fAlpha
convergence coefficient, input parameter, it should be positive number and <=1, for details see refer...
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.
Bool_t fFixAx
logical value of ax parameter, which allows to fix the parameter (not to fit).
Double_t fTyErr
error value of t parameter for 1D ridges in y direction
Double_t fSyCalc
calculated value of s parameter for 1D ridges in y direction
Int_t fXmin
first fitted channel in x direction
Double_t fTxyErr
error value of t parameter for 2D peaks
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:
void GetSigmaX(Double_t &sigmaX, Double_t &sigmaErrX)
This function gets the sigma x parameter and its error.
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 Erfc(Double_t x)
This function calculates error function of x.
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 * fFixPositionX1
[fNPeaks] array of logical values which allow to fix appropriate x positions of 1D ridges (not fit)....
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 Ourpowl(Double_t a, Int_t pw)
power function
Double_t * fAmpInitY1
[fNPeaks] array of initial values of amplitudes of 1D ridges in y direction, input parameters
Double_t fAxErr
error value of background ax parameter
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 fAyCalc
calculated value of background ay parameter
Double_t fSigmaInitY
initial value of sigma y parameter
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.
Double_t * fPositionErrX
[fNPeaks] array of error values of x positions of 2D peaks, output parameters
Double_t * fAmpInitX1
[fNPeaks] array of initial values of amplitudes of 1D ridges in x direction, input parameters
Bool_t * fFixPositionY1
[fNPeaks] array of logical values which allow to fix appropriate y positions of 1D ridges (not fit)....
void GetVolumeErrors(Double_t *volumeErrors)
This function gets errors of the volumes of fitted 2D peaks.
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...
~TSpectrum2Fit() override
Destructor.
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 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...
Double_t fTxInit
initial value of t parameter for 1D ridges in x direction (relative amplitude of tail),...
Double_t fSyInit
initial value of s parameter for 1D ridges in y direction (relative amplitude of step),...
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 * fAmpCalc
[fNPeaks] array of calculated values of amplitudes of 2D peaks, output parameters
void GetSigmaY(Double_t &sigmaY, Double_t &sigmaErrY)
This function gets the sigma y parameter and its error.
TSpectrum2Fit(void)
Default constructor.
Double_t fA0Calc
calculated value of background a0 parameter
Double_t * fPositionInitX
[fNPeaks] array of initial values of x positions of 2D peaks, input parameters
Int_t fAlphaOptim
optimization of convergence algorithm, possible values kFitAlphaHalving, kFitAlphaOptimal
Bool_t fFixSigmaX
logical value of sigma x parameter, which allows to fix the parameter (not to fit).
Double_t fSxErr
error value of s parameter for 1D ridges in x direction
Int_t fNPeaks
number of peaks present in fit, input parameter, it should be > 0
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:
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 * fAmpCalcY1
[fNPeaks] array of calculated values of amplitudes of 1D ridges in y direction, output parameters
Bool_t fFixAy
logical value of ay parameter, which allows to fix the parameter (not to fit).
void FitStiefel(Double_t **source)
This function fits the source spectrum.
Double_t Derfc(Double_t x)
This function calculates derivative of error function of x.
Double_t fSxInit
initial value of s parameter for 1D ridges in x direction (relative amplitude of step),...
Double_t fSigmaErrX
error value of sigma x parameter
Int_t fYmin
first fitted channel in y direction
Double_t * fPositionCalcX1
[fNPeaks] array of calculated x positions of 1D ridges, output 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 fBxErr
error value of b parameter for 1D ridges in x direction
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.
Double_t * fPositionCalcX
[fNPeaks] array of calculated values of x positions of 2D peaks, output parameters
Double_t fAxInit
initial value of background ax parameter(backgroud is estimated as a0+ax*x+ay*y)
Bool_t * fFixPositionX
[fNPeaks] array of logical values which allow to fix appropriate x positions of 2D peaks (not fit)....
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.
void StiefelInversion(Double_t **a, Int_t size)
This function calculates solution of the system of linear equations.
Double_t fSxyInit
initial value of s parameter for 2D peaks (relative amplitude of step), for details see html manual a...
Int_t fPower
possible values kFitPower2,4,6,8,10,12, for details see references. It applies only for Awmi fitting ...
Bool_t fFixA0
logical value of a0 parameter, which allows to fix the parameter (not to fit).
Double_t fAyInit
initial value of background ay parameter(backgroud is estimated as a0+ax*x+ay*y)
Double_t * fVolumeErr
[fNPeaks] array of volumes errors of 2D peaks, output parameters
Double_t Log(Double_t x)
Returns the natural logarithm of x.
Double_t Sqrt(Double_t x)
Returns the square root of x.
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.