ATLAS Offline Software
Loading...
Searching...
No Matches
MuonCalib::MuCCaFitterImplementation Class Reference

The straight line fitter for drift circles used by Calib. More...

#include <MuCCaFitter.h>

Collaboration diagram for MuonCalib::MuCCaFitterImplementation:

Public Member Functions

void Computelin (double x1, double y1, double r1, double x2, double y2, double r2)
void Computelinparnew (double x1, double y1, double r1, double x2, double y2, double r2)
void computehitsfromcircles (double x0, double y0, double r0, double x1, double y1, double r1, double a, double b)
void Computehitpsemes (int nhit, const std::vector< double > &xcirc, const std::vector< double > &ycirc, const std::vector< double > &rcirc, double a, double b)
void Computeparam3 (int number_of_hits, const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &r, const std::vector< double > &sr)
double get_a ()
double get_b ()
double get_da ()
double get_db ()
double get_corrab ()
double get_chi2f ()

Private Attributes

double m_angularcoefficient [4]
 parameters of the 4 tangent lines from first and last hit
double m_bfpar [4]
double m_xout0
 track points from 2 circles
double m_yout0
double m_xout1
double m_yout1
double m_xpoint [100]
 track points
double m_ypoint [100]
double m_aoutn
 results
double m_sig2a
 track slope's variance
double m_bout
 track constant term
double m_sig2b
 error on track constant term
double m_corrab
 correlation term
double m_chi2f
 chisquared

Detailed Description

Member Function Documentation

◆ Computehitpsemes()

void MuonCalib::MuCCaFitterImplementation::Computehitpsemes ( int nhit,
const std::vector< double > & xcirc,
const std::vector< double > & ycirc,
const std::vector< double > & rcirc,
double a,
double b )

Definition at line 467 of file MuCCaFitter.cxx.

468 {
469 /************************/
470 /* Compute TRACK POINTS */
471 /************************/
472 int i, nhitc;
473 // double xpoint[nhit],ypoint[nhit];
474 /* loop over hit couples */
475 nhitc = nhit / 2;
476 for (i = 0; i < nhitc; i++) {
477 /* Compute TRACK POINTS for 2 circles */
478 computehitsfromcircles(xcirc[2 * i], ycirc[2 * i], rcirc[2 * i], xcirc[2 * i + 1], ycirc[2 * i + 1], rcirc[2 * i + 1], a, b);
479 m_xpoint[2 * i] = m_xout0;
480 m_ypoint[2 * i] = m_yout0;
481 m_xpoint[2 * i + 1] = m_xout1;
482 m_ypoint[2 * i + 1] = m_yout1;
483 }
484 if (2 * nhitc != nhit) {
485 /* Compute TRACK POINTS for 2 circles */
486 computehitsfromcircles(xcirc[nhit - 2], ycirc[nhit - 2], rcirc[nhit - 2], xcirc[nhit - 1], ycirc[nhit - 1], rcirc[nhit - 1], a,
487 b);
488 m_xpoint[nhit - 2] = m_xout0;
489 m_ypoint[nhit - 2] = m_yout0;
490 m_xpoint[nhit - 1] = m_xout1;
491 m_ypoint[nhit - 1] = m_yout1;
492 }
493 }
static Double_t a
void computehitsfromcircles(double x0, double y0, double r0, double x1, double y1, double r1, double a, double b)
double m_xpoint[100]
track points
Definition MuCCaFitter.h:96
double m_xout0
track points from 2 circles
Definition MuCCaFitter.h:91

◆ computehitsfromcircles()

void MuonCalib::MuCCaFitterImplementation::computehitsfromcircles ( double x0,
double y0,
double r0,
double x1,
double y1,
double r1,
double a,
double b )

Definition at line 495 of file MuCCaFitter.cxx.

496 {
497 /**************************************/
498 /* Compute TRACK POINTS for 2 circles */
499 /**************************************/
500 double ang[4], bpar[4];
501 int i;
502 double xout0 = 0, yout0 = 0;
503 double xout1 = 0, yout1 = 0;
504 double dist, dist1, xint, yint, xref, yref;
505
506 /* Compute the 4 lines tangent to 2 circles */
507 Computelin(x0, y0, r0, x1, y1, r1);
508 for (i = 0; i < 4; i++) {
509 bpar[i] = m_bfpar[i];
510 ang[i] = m_angularcoefficient[i];
511 }
512
513 dist = 9999999;
514 /* for each of the 4 tangents : */
515 for (i = 0; i < 4; i++) {
516 /* compute tangent point */
517 xint = (x0 + ang[i] * y0 - ang[i] * bpar[i]) / (1 + ang[i] * ang[i]);
518 yint = ang[i] * (xint) + bpar[i];
519 /* compute point from reference line */
520 double bprime;
521 if (a != 0) {
522 bprime = y0 + x0 / a;
523 xref = (bprime - b) * a / (a * a + 1);
524 yref = (b + bprime * a * a) / (a * a + 1);
525 } else {
526 xref = x0;
527 yref = b;
528 }
529 /* choose tangent point closest to that of the reference lines (TRACK POINT)*/
530 dist1 = std::sqrt((xint - xref) * (xint - xref) + (yint - yref) * (yint - yref));
531 if (dist1 < dist) {
532 dist = dist1;
533 xout0 = xint;
534 yout0 = yint;
535 }
536 }
537 dist = 9999999;
538 for (i = 0; i < 4; i++) {
539 /* compute tangent point */
540 xint = (x1 + ang[i] * y1 - ang[i] * bpar[i]) / (1 + ang[i] * ang[i]);
541 yint = ang[i] * (xint) + bpar[i];
542 /* compute point from reference line */
543 double bprime;
544 if (a != 0) {
545 bprime = y1 + x1 / a;
546 xref = (bprime - b) * a / (a * a + 1);
547 yref = (b + bprime * a * a) / (a * a + 1);
548 } else {
549 xref = x1;
550 yref = b;
551 }
552 /* choose tangent point closest to that of the reference line (TRACK POINT)*/
553 dist1 = std::sqrt((xint - xref) * (xint - xref) + (yint - yref) * (yint - yref));
554 if (dist1 < dist) {
555 dist = dist1;
556 xout1 = xint;
557 yout1 = yint;
558 }
559 }
560 /* set global variables (method output)*/
561 m_xout0 = xout0;
562 m_yout0 = yout0;
563 m_xout1 = xout1;
564 m_yout1 = yout1;
565 }
double m_angularcoefficient[4]
parameters of the 4 tangent lines from first and last hit
Definition MuCCaFitter.h:88
void Computelin(double x1, double y1, double r1, double x2, double y2, double r2)

◆ Computelin()

void MuonCalib::MuCCaFitterImplementation::Computelin ( double x1,
double y1,
double r1,
double x2,
double y2,
double r2 )

Definition at line 375 of file MuCCaFitter.cxx.

375 {
376 /********************************************/
377 /* Compute the 4 lines tangent to 2 circles */
378 /********************************************/
379 double delta{0}, averagephi{0}, dist{0}, dphiin{0}, dphiex{0};
380 double bpar[4], phi{0};
381 double bfparn[2];
382 bfparn[0] = 0;
383 bfparn[1] = 0;
384 double bcand[2][4];
385 int segnob[] = {1, -1, 1, -1};
386 int ncandid[4], ncand;
387 int i{0}, firsttime{0};
388 double angularcoefficient[4];
389 double bfpar[4];
390 double dy{0}, dx{0};
391
392 /* compute a parameters */
393 dx = x2 - x1;
394 dy = y2 - y1;
395 averagephi = std::atan2(dy, dx);
396 dist = std::hypot(dx, dy);
397
398 delta = r2 + r1;
399 dphiin = std::asin(delta / dist);
400
401 delta = r2 - r1;
402 dphiex = std::asin(delta / dist);
403
404 int f = 1;
405 phi = averagephi + f * dphiex;
406 if (phi < 0) { phi = 2 * M_PI + (phi); }
407 angularcoefficient[0] = std::tan(phi);
408
409 f = -1;
410 phi = averagephi + f * dphiex;
411 if (phi < 0) { phi = 2 * M_PI + (phi); }
412 angularcoefficient[1] = std::tan(phi);
413
414 f = 1;
415 phi = averagephi + f * dphiin;
416 if (phi < 0) { phi = 2 * M_PI + (phi); }
417 angularcoefficient[2] = std::tan(phi);
418
419 f = -1;
420 phi = averagephi + f * dphiin;
421 if (phi < 0) { phi = 2 * M_PI + (phi); }
422 angularcoefficient[3] = std::tan(phi);
423
424 /* Compute b parameters */
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;
431 ncand = 0;
432 if (std::abs(bpar[0] - bpar[2]) < delta) {
433 bfparn[ncand] = bpar[0];
434 ncand = ncand + 1;
435 }
436 if (std::abs(bpar[0] - bpar[3]) < delta) {
437 bfparn[ncand] = bpar[0];
438 ncand = ncand + 1;
439 }
440 if (std::abs(bpar[1] - bpar[2]) < delta) {
441 bfparn[ncand] = bpar[1];
442 ncand = ncand + 1;
443 }
444 if (std::abs(bpar[1] - bpar[3]) < delta) {
445 bfparn[ncand] = bpar[1];
446 ncand = ncand + 1;
447 }
448 bcand[0][i] = bfparn[0];
449 bcand[1][i] = bfparn[1];
450 ncandid[i] = ncand;
451 }
452 firsttime = 0;
453 for (i = 0; i < 4; i++) {
454 if (ncandid[i] == 1) {
455 bfpar[i] = bcand[0][i];
456 } else {
457 bfpar[i] = bcand[firsttime][i];
458 firsttime++;
459 }
460 }
461 /* set global variables (method output)*/
462 for (int i = 0; i < 4; i++) {
463 m_bfpar[i] = bfpar[i];
464 m_angularcoefficient[i] = angularcoefficient[i];
465 }
466 }
#define M_PI
Scalar phi() const
phi method

◆ Computelinparnew()

void MuonCalib::MuCCaFitterImplementation::Computelinparnew ( double x1,
double y1,
double r1,
double x2,
double y2,
double r2 )

Definition at line 278 of file MuCCaFitter.cxx.

278 {
279 /********************************************/
280 /* Compute the 4 lines tangent to 2 circles */
281 /********************************************/
282 double delta;
283 double averagephi, dist, dphiin, dphiex;
284 double bpar[4], phi;
285 double bfparn[2];
286 bfparn[0] = 0;
287 bfparn[1] = 0;
288 // int segnobfn[2];
289 // int segnoc[2][4],ncandid[4],ncand;
290 int ncandid[4], ncand;
291 double bcand[2][4];
292 int segnob[] = {1, -1, 1, -1};
293 int firsttime, i;
294 double angularcoefficient[4];
295 double bfpar[4];
296 // int segnobf[4];
297 double dy, dx;
298
299 /* compute a parameters */
300 dx = x2 - x1;
301 dy = y2 - y1;
302 averagephi = std::atan2(dy, dx);
303 dist = std::sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));
304 delta = r2 + r1;
305
306 dphiin = std::asin(delta / dist);
307
308 delta = r2 - r1;
309 dphiex = std::asin(delta / dist);
310
311 int f = 1;
312 phi = averagephi + f * dphiex;
313 if (phi < 0) { phi = 2 * M_PI + (phi); }
314 angularcoefficient[0] = std::tan(phi);
315
316 f = -1;
317 phi = averagephi + f * dphiex;
318 if (phi < 0) { phi = 2 * M_PI + (phi); }
319 angularcoefficient[1] = std::tan(phi);
320
321 f = 1;
322 phi = averagephi + f * dphiin;
323 if (phi < 0) { phi = 2 * M_PI + (phi); }
324 angularcoefficient[2] = std::tan(phi);
325
326 f = -1;
327 phi = averagephi + f * dphiin;
328 if (phi < 0) { phi = 2 * M_PI + (phi); }
329 angularcoefficient[3] = std::tan(phi);
330
331 /* Compute b parameters */
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;
338 ncand = 0;
339 if (std::abs(bpar[0] - bpar[2]) < delta) {
340 bfparn[ncand] = bpar[0];
341 ncand = ncand + 1;
342 }
343 if (std::abs(bpar[0] - bpar[3]) < delta) {
344 bfparn[ncand] = bpar[0];
345 ncand = ncand + 1;
346 }
347 if (std::abs(bpar[1] - bpar[2]) < delta) {
348 bfparn[ncand] = bpar[1];
349 ncand = ncand + 1;
350 }
351 if (std::abs(bpar[1] - bpar[3]) < delta) {
352 bfparn[ncand] = bpar[1];
353 ncand = ncand + 1;
354 }
355 bcand[0][i] = bfparn[0];
356 bcand[1][i] = bfparn[1];
357 ncandid[i] = ncand;
358 }
359 firsttime = 0;
360 for (i = 0; i < 4; i++) {
361 if (ncandid[i] == 1) {
362 bfpar[i] = bcand[0][i];
363 } else {
364 bfpar[i] = bcand[firsttime][i];
365 firsttime++;
366 }
367 }
368 /* set global variables (method output)*/
369 for (i = 0; i < 4; i++) {
370 m_bfpar[i] = bfpar[i];
371 m_angularcoefficient[i] = angularcoefficient[i];
372 }
373 }

◆ Computeparam3()

void MuonCalib::MuCCaFitterImplementation::Computeparam3 ( int number_of_hits,
const std::vector< double > & x,
const std::vector< double > & y,
const std::vector< double > & r,
const std::vector< double > & sr )

Definition at line 155 of file MuCCaFitter.cxx.

156 {
157 /***************************************/
158 /* Fit a line to n=number_of_hits hits */
159 /***************************************/
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;
165 /* compute 4 tanget lines from first and last hit */
166 fhit = 0;
167 lhit = number_of_hits - 1;
168 Computelinparnew(x[fhit], y[fhit], r[fhit], x[lhit], y[lhit], r[lhit]);
169 for (int ii = 0; ii < 4; ii++) {
170 bref[ii] = m_bfpar[ii];
171 aref[ii] = m_angularcoefficient[ii];
172 }
173 /* choose the REFERENCE LINE (aref[icouple],bref[icouple])*/
174 for (i = 0; i < 4; i++) {
175 for (j = 0; j < number_of_hits; j++) {
176 double d;
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];
180 }
181 double chi2out = 0;
182 for (int ii = 1; ii < number_of_hits - 1; ii++) { chi2out += rsub[ii] * rsub[ii] / (sr[ii] * sr[ii]); }
183 chi2outn[i] = chi2out;
184 }
185 chi2min = 999999999.;
186 for (i = 0; i < 4; i++) {
187 if (chi2outn[i] < chi2min) {
188 chi2min = chi2outn[i];
189 icouple = i;
190 }
191 }
192 /* Compute TRACK POINTS */
193 Computehitpsemes(number_of_hits, x, y, r, aref[icouple], bref[icouple]);
194 for (i = 0; i < number_of_hits; i++) {
195 xout[i] = m_xpoint[i];
196 yout[i] = m_ypoint[i];
197 }
198 /* Compute a & b parameters of the track from TRACK POINTS */
199 // int idim;
200 // ifail;
201 double aoutn, bout, sig2a, sig2b, corrab;
202 double temp, det;
203 double hesse[2][2];
204 double W, WX, WX2, WY, WXY;
205 // double errormatrix[2][2];
206 W = 0.;
207 WX = 0.;
208 WX2 = 0.;
209 WY = 0.;
210 WXY = 0;
211 for (int i = 0; i < number_of_hits; i++) {
212 temp = 1. / (sr[i] * sr[i]);
213 W += temp;
214 WX += xout[i] * temp;
215 WX2 += xout[i] * xout[i] * temp;
216 WY += yout[i] * temp;
217 WXY += xout[i] * yout[i] * temp;
218 }
219 det = W * WX2 - WX * WX;
220 aoutn = (W * WXY - WY * WX) / det;
221 bout = (WY * WX2 - WX * WXY) / det;
222 hesse[1][1] = W;
223 hesse[0][0] = WX2;
224 hesse[1][0] = WX;
225 hesse[0][1] = WX;
226 // idim = 2;
227 /* invert hessian matrix */
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);
232 sig2a = hesse[0][0];
233 sig2b = hesse[1][1];
234 corrab = hesse[1][0];
235 double deno, btest, bs1, bs2, db1, db2;
236 btest = bout;
237 bout = 0;
238 deno = 0;
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);
244 if (db1 < db2) {
245 bout += bs1 / (sr[i] * sr[i]);
246 } else {
247 bout += bs2 / (sr[i] * sr[i]);
248 }
249 deno += 1 / (sr[i] * sr[i]);
250 }
251 bout /= deno;
252 /* compute chi2 from residuals*/
253 double resd;
254 double chi2f = 0;
255 for (int i = 0; i < number_of_hits; i++) {
256 double xi, yi;
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;
261 }
262 /* set global variables (method output)*/
263 m_aoutn = aoutn;
264 m_sig2a = sig2a;
265 m_bout = bout;
266 m_sig2b = sig2b;
267 m_corrab = corrab;
268 m_chi2f = chi2f;
269 }
#define y
#define x
void Computelinparnew(double x1, double y1, double r1, double x2, double y2, double r2)
double m_bout
track constant term
double m_sig2b
error on track constant term
double m_sig2a
track slope's variance
double m_corrab
correlation term
void Computehitpsemes(int nhit, const std::vector< double > &xcirc, const std::vector< double > &ycirc, const std::vector< double > &rcirc, double a, double b)
int r
Definition globals.cxx:22

◆ get_a()

double MuonCalib::MuCCaFitterImplementation::get_a ( )

Definition at line 271 of file MuCCaFitter.cxx.

271{ return m_aoutn; }

◆ get_b()

double MuonCalib::MuCCaFitterImplementation::get_b ( )

Definition at line 272 of file MuCCaFitter.cxx.

272{ return m_bout; }

◆ get_chi2f()

double MuonCalib::MuCCaFitterImplementation::get_chi2f ( )

Definition at line 276 of file MuCCaFitter.cxx.

276{ return m_chi2f; }

◆ get_corrab()

double MuonCalib::MuCCaFitterImplementation::get_corrab ( )

Definition at line 275 of file MuCCaFitter.cxx.

275{ return m_corrab; }

◆ get_da()

double MuonCalib::MuCCaFitterImplementation::get_da ( )

Definition at line 273 of file MuCCaFitter.cxx.

273{ return m_sig2a; }

◆ get_db()

double MuonCalib::MuCCaFitterImplementation::get_db ( )

Definition at line 274 of file MuCCaFitter.cxx.

274{ return m_sig2b; }

Member Data Documentation

◆ m_angularcoefficient

double MuonCalib::MuCCaFitterImplementation::m_angularcoefficient[4]
private

parameters of the 4 tangent lines from first and last hit

Definition at line 88 of file MuCCaFitter.h.

◆ m_aoutn

double MuonCalib::MuCCaFitterImplementation::m_aoutn
private

results

track slope

Definition at line 99 of file MuCCaFitter.h.

◆ m_bfpar

double MuonCalib::MuCCaFitterImplementation::m_bfpar[4]
private

Definition at line 89 of file MuCCaFitter.h.

◆ m_bout

double MuonCalib::MuCCaFitterImplementation::m_bout
private

track constant term

Definition at line 101 of file MuCCaFitter.h.

◆ m_chi2f

double MuonCalib::MuCCaFitterImplementation::m_chi2f
private

chisquared

Definition at line 104 of file MuCCaFitter.h.

◆ m_corrab

double MuonCalib::MuCCaFitterImplementation::m_corrab
private

correlation term

Definition at line 103 of file MuCCaFitter.h.

◆ m_sig2a

double MuonCalib::MuCCaFitterImplementation::m_sig2a
private

track slope's variance

Definition at line 100 of file MuCCaFitter.h.

◆ m_sig2b

double MuonCalib::MuCCaFitterImplementation::m_sig2b
private

error on track constant term

Definition at line 102 of file MuCCaFitter.h.

◆ m_xout0

double MuonCalib::MuCCaFitterImplementation::m_xout0
private

track points from 2 circles

Definition at line 91 of file MuCCaFitter.h.

◆ m_xout1

double MuonCalib::MuCCaFitterImplementation::m_xout1
private

Definition at line 93 of file MuCCaFitter.h.

◆ m_xpoint

double MuonCalib::MuCCaFitterImplementation::m_xpoint[100]
private

track points

Definition at line 96 of file MuCCaFitter.h.

◆ m_yout0

double MuonCalib::MuCCaFitterImplementation::m_yout0
private

Definition at line 92 of file MuCCaFitter.h.

◆ m_yout1

double MuonCalib::MuCCaFitterImplementation::m_yout1
private

Definition at line 94 of file MuCCaFitter.h.

◆ m_ypoint

double MuonCalib::MuCCaFitterImplementation::m_ypoint[100]
private

Definition at line 97 of file MuCCaFitter.h.


The documentation for this class was generated from the following files: