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);
 
   51         std::vector<double> 
r(
N), 
sr(
N), 
w(
N), rw(
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);
 
  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;
 
  317         phi = averagephi + 
f * dphiex;
 
  322         phi = averagephi + 
f * dphiin;
 
  327         phi = averagephi + 
f * dphiin;
 
  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;
 
  410         phi = averagephi + 
f * dphiex;
 
  415         phi = averagephi + 
f * dphiin;
 
  420         phi = averagephi + 
f * dphiin;
 
  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));