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],