38 int iter, iterMax = 20;
43 double M0, Mx, My, Mz, Mxy, Mxx, Myy, Mxz, Myz, Mzz, Mxz2, Myz2, Cov_xy;
44 double A0, A1, A2, A22, A3, A33, epsilon = 0.000000000001;
45 double Dy, xnew, xold, ynew, yold = 100000000000.;
53 const double copt = 0.2;
54 const double zero = 0.0001;
55 const double SigmaMin = 0.18;
56 double amax = -100000;
67 if (fRobust < 1 || fRobust > 2) { riterMax = 1; }
68 if (
fRobust == 1) { riterMax = 4; }
69 if (
fRobust == 2) { riterMax = 4; }
71 for (
int i = 0;
i < nofHits;
i++) {
77 for (
int i = 0;
i < nofHits;
i++) {
78 if (a[
i] > amax) amax = a[
i];
79 if (a[
i] < amin) amin = a[
i];
82 for (
int i = 0;
i < nofHits;
i++)
85 for (
int riter = 0; riter < riterMax; riter++) {
89 for (
int i = 0;
i < nofHits;
i++) {
99 Mxx = Myy = Mxy = Mxz = Myz = Mzz = 0.;
101 for (
int i = 0;
i < nofHits;
i++) {
104 Zi = Xi * Xi + Yi * Yi;
106 Mxy += Xi * Yi * w[
i];
107 Mxx += Xi * Xi * w[
i];
108 Myy += Yi * Yi * w[
i];
109 Mxz += Xi * Zi * w[
i];
110 Myz += Yi * Zi * w[
i];
111 Mzz += Zi * Zi * w[
i];
122 Cov_xy = Mxx * Myy - Mxy * Mxy;
127 A2 = -3. * Mz * Mz - Mzz;
128 A1 = Mzz * Mz + 4. * Cov_xy * Mz - Mxz2 - Myz2 - Mz * Mz * Mz;
129 A0 = Mxz2 * Myy + Myz2 * Mxx - Mzz * Cov_xy - 2. * Mxz * Myz * Mxy
138 for (iter = 0; iter < iterMax; iter++) {
139 ynew = A0 + xnew * (A1 + xnew * (A2 + xnew * A3));
146 Dy = A1 + xnew * (A22 + xnew * A33);
148 xnew = xold - ynew / Dy;
150 if (
fabs((xnew - xold) / xnew) < epsilon)
break;
153 if (iter == iterMax - 1) {
164 DET = xnew * xnew - xnew * Mz + Cov_xy;
165 centerX = (Mxz * (Myy - xnew) - Myz * Mxy) / DET / 2.;
166 centerY = (Myz * (Mxx - xnew) - Mxz * Mxy) / DET / 2.;
167 radius =
sqrt(centerX * centerX + centerY * centerY - GAM);
168 centerX = (Mxz * (Myy - xnew) - Myz * Mxy) / DET / 2. + Mx;
169 centerY = (Myz * (Mxx - xnew) - Mxz * Mxy) / DET / 2. + My;
171 if (riter < riterMax - 1) {
172 for (
int i = 0;
i < nofHits;
i++) {
174 sqrt(pow((centerX -
x[
i]), 2) + pow((centerY -
y[
i]), 2)) - radius;
179 for (
unsigned int i = 0;
i <
d.size();
i++) {
180 sig1 += w[
i] *
d[
i] *
d[
i];
183 sigma =
sqrt(sig1 / sig2);
184 if (sigma < SigmaMin) sigma = SigmaMin;
185 ctsigma = ct * sigma;
190 for (
unsigned int i = 0;
i <
d.size();
i++) {
192 if (
d[
i] <= ctsigma) {
193 weight = pow((1 - pow((
d[
i] / ctsigma), 2)), 2);
194 if (weight < zero) weight = zero;
203 weight /= 1 + copt *
exp(pow((
d[
i] / sigma), 2));
204 if (weight < zero) weight = zero;