10 #include "GaudiKernel/MsgStream.h"
31 if (
N < 2) {
return false; }
38 for (
int i = 0;
i <
N; ++
i) {
42 log << MSG::WARNING <<
"fit() TO FEW HITS SELECTED" <<
endmsg;
49 double S(0), Sy(0), Sz(0), Zc{0}, Yc{0};
50 std::vector<double>
y(
N),
x(
N),
z(
N);
56 y[ii] =
getY(
h.localPosition());
57 z[ii] =
getZ(
h.localPosition());
58 r[ii] = std::abs(
h.driftRadius());
60 if (
h.sigma2DriftRadius() > 0) {
61 w[ii] = 1. / (
h.sigma2DriftRadius());
62 sr[ii] = std::sqrt(
h.sigma2DriftRadius());
68 log <<
MSG::DEBUG <<
"fit() MuCCaFitter: (" <<
x[ii] <<
"," <<
y[ii] <<
") R = " <<
r[ii] <<
" W " <<
w[ii] <<
endmsg;
69 rw[ii] =
r[ii] *
w[ii];
84 std::unique_ptr<MuCCaFitterImplementation> Fitter = std::make_unique<MuCCaFitterImplementation>();
86 for (
int i = 0;
i != ii;
i++) {
87 log <<
MSG::DEBUG <<
"fit() MuCCaFitter hits passed to computepam Z=" <<
x[
i] <<
" phi=" <<
y[
i] <<
" zzz=" <<
z[
i]
88 <<
" drift=" <<
r[
i] <<
" error=" <<
w[
i] <<
endmsg;
91 Fitter->Computeparam3(ii,
z,
y,
r,
sr);
93 log <<
MSG::DEBUG <<
"fit() MuCCaFitter computed track a=" << Fitter->get_a() <<
" da= " << Fitter->get_da()
94 <<
" b= " << Fitter->get_b() <<
" db= " << Fitter->get_db() <<
" corrab= " << Fitter->get_corrab()
95 <<
" chi2f= " << Fitter->get_chi2f() <<
endmsg;
98 double afit = Fitter->get_a();
99 double chi2f = Fitter->get_chi2f();
101 std::vector<double> dist(
N), ddist(
N);
104 if (theta < 0.) theta =
M_PI + theta;
109 log <<
MSG::DEBUG <<
"fit() MuCCaFitter>>> theta= " << theta <<
" sinus= " << sinus <<
" cosin= " << cosin <<
endmsg;
111 <<
" Yc= " << Yc <<
endmsg;
113 double Syy(0), Szy(0), Syyzz(0), Att(0);
115 for (
int i = 0;
i <
N; ++
i) {
117 Syy += (
y[
i] - Yc) * (
y[
i] - Yc) *
w[
i];
118 Szy += (
y[
i] - Yc) * (
z[
i] - Zc) *
w[
i];
119 Syyzz += ((
y[
i] - Yc) - (
z[
i] - Zc)) * ((
y[
i] - Yc) + (
z[
i] - Zc)) *
w[
i];
120 dis = (
y[
i] - Yc) * cosin - (
z[
i] - Zc) * sinus;
127 Att = Syy + cosin * (2 * sinus * Szy - cosin * Syyzz);
130 log <<
MSG::DEBUG <<
"fit() MuCCaFitter>>> d= " <<
d <<
" R= " << R <<
" S= " <<
S <<
" Att= " << Att <<
endmsg;
131 for (
int i = 0;
i <
N; ++
i) {
133 dist[
i] = cosin * (
y[
i] - Yc) - sinus * (
z[
i] - Zc) -
d;
134 double dth = -(sinus * (
y[
i] - Yc) + cosin * (
z[
i] - Zc)) * (std::sqrt(1. / Att));
135 ddist[
i] = std::sqrt(dth * dth + (1. /
S));
142 log <<
MSG::DEBUG <<
"fit() New line: position " << npos <<
" direction " << ndir <<
" chi2f " << chi2f <<
endmsg;
144 seg.
set(chi2f / (
N - 2), npos, ndir);
148 hit_ptr->setDistanceToTrack(dist[
i], ddist[
i]);
156 const std::vector<double>&
sr) {
160 double xout[100], yout[100];
161 double dist[100], rsub[100];
162 double chi2outn[4], chi2min;
163 double aref[4], bref[4];
164 int fhit,
i, j, icouple = 0, lhit;
167 lhit = number_of_hits - 1;
169 for (
int ii = 0; ii < 4; ii++) {
174 for (
i = 0;
i < 4;
i++) {
175 for (j = 0; j < number_of_hits; j++) {
177 d = std::abs(aref[
i] *
x[j] + bref[
i] -
y[j]);
178 dist[j] =
d / std::sqrt(1. + aref[
i] * aref[
i]);
179 rsub[j] =
r[j] - dist[j];
182 for (
int ii = 1; ii < number_of_hits - 1; ii++) { chi2out += rsub[ii] * rsub[ii] / (
sr[ii] *
sr[ii]); }
183 chi2outn[
i] = chi2out;
185 chi2min = 999999999.;
186 for (
i = 0;
i < 4;
i++) {
187 if (chi2outn[
i] < chi2min) {
188 chi2min = chi2outn[
i];
194 for (
i = 0;
i < number_of_hits;
i++) {
201 double aoutn, bout, sig2a, sig2b, corrab;
204 double W, WX, WX2, WY, WXY;
211 for (
int i = 0;
i < number_of_hits;
i++) {
212 temp = 1. / (
sr[
i] *
sr[
i]);
214 WX += xout[
i] * temp;
215 WX2 += xout[
i] * xout[
i] * temp;
216 WY += yout[
i] * temp;
217 WXY += xout[
i] * yout[
i] * temp;
219 det =
W * WX2 - WX * WX;
220 aoutn = (
W * WXY - WY * WX) /
det;
221 bout = (WY * WX2 - WX * WXY) /
det;
228 hesse[1][1] = hesse[0][0] /
det;
229 hesse[0][0] = hesse[1][1] /
det;
230 hesse[1][0] = -1. * (hesse[0][1] /
det);
231 hesse[0][1] = -1. * (hesse[0][1] /
det);
234 corrab = hesse[1][0];
235 double deno, btest, bs1, bs2,
db1,
db2;
239 for (
int i = 0;
i < number_of_hits;
i++) {
240 bs1 = (
y[
i] - aoutn *
x[
i] -
r[
i] * std::sqrt(1 + aoutn * aoutn));
241 bs2 = (
y[
i] - aoutn *
x[
i] +
r[
i] * std::sqrt(1 + aoutn * aoutn));
242 db1 = std::abs(bs1 - btest);
243 db2 = std::abs(bs2 - btest);
245 bout += bs1 / (
sr[
i] *
sr[
i]);
247 bout += bs2 / (
sr[
i] *
sr[
i]);
249 deno += 1 / (
sr[
i] *
sr[
i]);
255 for (
int i = 0;
i < number_of_hits;
i++) {
257 xi = (aoutn *
y[
i] - aoutn * bout +
x[
i]) / (1. + aoutn * aoutn);
258 yi = aoutn * xi + bout;
259 resd = (std::sqrt((xi -
x[
i]) * (xi -
x[
i]) + (yi -
y[
i]) * (yi -
y[
i])) -
r[
i]) /
sr[
i];
260 chi2f += resd * resd;
283 double averagephi, dist, dphiin, dphiex;
290 int ncandid[4], ncand;
292 int segnob[] = {1, -1, 1, -1};
294 double angularcoefficient[4];
302 averagephi = std::atan2(
dy,
dx);
306 dphiin = std::asin(delta / dist);
309 dphiex = std::asin(delta / dist);
312 phi = averagephi +
f * dphiex;
313 if (phi < 0) { phi = 2 *
M_PI + (phi); }
314 angularcoefficient[0] =
std::tan(phi);
317 phi = averagephi +
f * dphiex;
318 if (phi < 0) { phi = 2 *
M_PI + (phi); }
319 angularcoefficient[1] =
std::tan(phi);
322 phi = averagephi +
f * dphiin;
323 if (phi < 0) { phi = 2 *
M_PI + (phi); }
324 angularcoefficient[2] =
std::tan(phi);
327 phi = averagephi +
f * dphiin;
328 if (phi < 0) { phi = 2 *
M_PI + (phi); }
329 angularcoefficient[3] =
std::tan(phi);
332 for (
i = 0;
i < 4;
i++) {
333 bpar[0] =
y1 - angularcoefficient[
i] *
x1 + segnob[0] * r1 * std::sqrt(1 + angularcoefficient[
i] * angularcoefficient[
i]);
334 bpar[1] =
y1 - angularcoefficient[
i] *
x1 + segnob[1] * r1 * std::sqrt(1 + angularcoefficient[
i] * angularcoefficient[
i]);
335 bpar[2] =
y2 - angularcoefficient[
i] *
x2 + segnob[2] * r2 * std::sqrt(1 + angularcoefficient[
i] * angularcoefficient[
i]);
336 bpar[3] =
y2 - angularcoefficient[
i] *
x2 + segnob[3] * r2 * std::sqrt(1 + angularcoefficient[
i] * angularcoefficient[
i]);
337 double delta = 0.00001;
339 if (std::abs(bpar[0] - bpar[2]) < delta) {
340 bfparn[ncand] = bpar[0];
343 if (std::abs(bpar[0] - bpar[3]) < delta) {
344 bfparn[ncand] = bpar[0];
347 if (std::abs(bpar[1] - bpar[2]) < delta) {
348 bfparn[ncand] = bpar[1];
351 if (std::abs(bpar[1] - bpar[3]) < delta) {
352 bfparn[ncand] = bpar[1];
355 bcand[0][
i] = bfparn[0];
356 bcand[1][
i] = bfparn[1];
360 for (
i = 0;
i < 4;
i++) {
361 if (ncandid[
i] == 1) {
362 bfpar[
i] = bcand[0][
i];
364 bfpar[
i] = bcand[firsttime][
i];
369 for (
i = 0;
i < 4;
i++) {
379 double delta{0}, averagephi{0}, dist{0}, dphiin{0}, dphiex{0};
380 double bpar[4], phi{0};
385 int segnob[] = {1, -1, 1, -1};
386 int ncandid[4], ncand;
387 int i{0}, firsttime{0};
388 double angularcoefficient[4];
395 averagephi = std::atan2(
dy,
dx);
396 dist = std::hypot(
dx,
dy);
399 dphiin = std::asin(delta / dist);
402 dphiex = std::asin(delta / dist);
405 phi = averagephi +
f * dphiex;
406 if (phi < 0) { phi = 2 *
M_PI + (phi); }
407 angularcoefficient[0] =
std::tan(phi);
410 phi = averagephi +
f * dphiex;
411 if (phi < 0) { phi = 2 *
M_PI + (phi); }
412 angularcoefficient[1] =
std::tan(phi);
415 phi = averagephi +
f * dphiin;
416 if (phi < 0) { phi = 2 *
M_PI + (phi); }
417 angularcoefficient[2] =
std::tan(phi);
420 phi = averagephi +
f * dphiin;
421 if (phi < 0) { phi = 2 *
M_PI + (phi); }
422 angularcoefficient[3] =
std::tan(phi);
425 for (
i = 0;
i < 4;
i++) {
426 bpar[0] =
y1 - angularcoefficient[
i] *
x1 + segnob[0] * r1 * std::sqrt(1 + angularcoefficient[
i] * angularcoefficient[
i]);
427 bpar[1] =
y1 - angularcoefficient[
i] *
x1 + segnob[1] * r1 * std::sqrt(1 + angularcoefficient[
i] * angularcoefficient[
i]);
428 bpar[2] =
y2 - angularcoefficient[
i] *
x2 + segnob[2] * r2 * std::sqrt(1 + angularcoefficient[
i] * angularcoefficient[
i]);
429 bpar[3] =
y2 - angularcoefficient[
i] *
x2 + segnob[3] * r2 * std::sqrt(1 + angularcoefficient[
i] * angularcoefficient[
i]);
430 double delta = 0.00001;
432 if (std::abs(bpar[0] - bpar[2]) < delta) {
433 bfparn[ncand] = bpar[0];
436 if (std::abs(bpar[0] - bpar[3]) < delta) {
437 bfparn[ncand] = bpar[0];
440 if (std::abs(bpar[1] - bpar[2]) < delta) {
441 bfparn[ncand] = bpar[1];
444 if (std::abs(bpar[1] - bpar[3]) < delta) {
445 bfparn[ncand] = bpar[1];
448 bcand[0][
i] = bfparn[0];
449 bcand[1][
i] = bfparn[1];
453 for (
i = 0;
i < 4;
i++) {
454 if (ncandid[
i] == 1) {
455 bfpar[
i] = bcand[0][
i];
457 bfpar[
i] = bcand[firsttime][
i];
462 for (
int i = 0;
i < 4;
i++) {
468 const std::vector<double>& rcirc,
double a,
double b) {
476 for (
i = 0;
i < nhitc;
i++) {
484 if (2 * nhitc != nhit) {
486 computehitsfromcircles(xcirc[nhit - 2], ycirc[nhit - 2], rcirc[nhit - 2], xcirc[nhit - 1], ycirc[nhit - 1], rcirc[nhit - 1],
a,
500 double ang[4], bpar[4];
502 double xout0 = 0, yout0 = 0;
503 double xout1 = 0, yout1 = 0;
504 double dist, dist1, xint, yint, xref, yref;
508 for (
i = 0;
i < 4;
i++) {
515 for (
i = 0;
i < 4;
i++) {
517 xint = (x0 + ang[
i] * y0 - ang[
i] * bpar[
i]) / (1 + ang[
i] * ang[
i]);
518 yint = ang[
i] * (xint) + bpar[
i];
522 bprime = y0 + x0 /
a;
523 xref = (bprime -
b) *
a / (
a *
a + 1);
524 yref = (
b + bprime *
a *
a) / (
a *
a + 1);
530 dist1 = std::sqrt((xint - xref) * (xint - xref) + (yint - yref) * (yint - yref));
538 for (
i = 0;
i < 4;
i++) {
540 xint = (
x1 + ang[
i] *
y1 - ang[
i] * bpar[
i]) / (1 + ang[
i] * ang[
i]);
541 yint = ang[
i] * (xint) + bpar[
i];
545 bprime =
y1 +
x1 /
a;
546 xref = (bprime -
b) *
a / (
a *
a + 1);
547 yref = (
b + bprime *
a *
a) / (
a *
a + 1);
553 dist1 = std::sqrt((xint - xref) * (xint - xref) + (yint - yref) * (yint - yref));