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 474 of file MuCCaFitter.cxx.

475 {
476 /************************/
477 /* Compute TRACK POINTS */
478 /************************/
479 int i, nhitc;
480 // double xpoint[nhit],ypoint[nhit];
481 /* loop over hit couples */
482 nhitc = nhit / 2;
483 for (i = 0; i < nhitc; i++) {
484 /* Compute TRACK POINTS for 2 circles */
485 computehitsfromcircles(xcirc[2 * i], ycirc[2 * i], rcirc[2 * i], xcirc[2 * i + 1], ycirc[2 * i + 1], rcirc[2 * i + 1], a, b);
486 m_xpoint[2 * i] = m_xout0;
487 m_ypoint[2 * i] = m_yout0;
488 m_xpoint[2 * i + 1] = m_xout1;
489 m_ypoint[2 * i + 1] = m_yout1;
490 }
491 if (2 * nhitc != nhit) {
492 /* Compute TRACK POINTS for 2 circles */
493 computehitsfromcircles(xcirc[nhit - 2], ycirc[nhit - 2], rcirc[nhit - 2], xcirc[nhit - 1], ycirc[nhit - 1], rcirc[nhit - 1], a,
494 b);
495 m_xpoint[nhit - 2] = m_xout0;
496 m_ypoint[nhit - 2] = m_yout0;
497 m_xpoint[nhit - 1] = m_xout1;
498 m_ypoint[nhit - 1] = m_yout1;
499 }
500 }
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 502 of file MuCCaFitter.cxx.

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

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

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

◆ 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 165 of file MuCCaFitter.cxx.

166 {
167 /***************************************/
168 /* Fit a line to n=number_of_hits hits */
169 /***************************************/
170 double xout[100], yout[100];
171 double dist[100], rsub[100];
172 double chi2outn[4], chi2min;
173 double aref[4], bref[4];
174 int fhit, i, j, icouple = 0, lhit;
175 /* compute 4 tanget lines from first and last hit */
176 fhit = 0;
177 lhit = number_of_hits - 1;
178 Computelinparnew(x[fhit], y[fhit], r[fhit], x[lhit], y[lhit], r[lhit]);
179 for (int ii = 0; ii < 4; ii++) {
180 bref[ii] = m_bfpar[ii];
181 aref[ii] = m_angularcoefficient[ii];
182 }
183 /* choose the REFERENCE LINE (aref[icouple],bref[icouple])*/
184 for (i = 0; i < 4; i++) {
185 for (j = 0; j < number_of_hits; j++) {
186 double d;
187 d = std::abs(aref[i] * x[j] + bref[i] - y[j]);
188 dist[j] = d / std::sqrt(1. + aref[i] * aref[i]);
189 rsub[j] = r[j] - dist[j];
190 }
191 double chi2out = 0;
192 for (int ii = 1; ii < number_of_hits - 1; ii++) { chi2out += rsub[ii] * rsub[ii] / (sr[ii] * sr[ii]); }
193 chi2outn[i] = chi2out;
194 }
195 chi2min = 999999999.;
196 for (i = 0; i < 4; i++) {
197 if (chi2outn[i] < chi2min) {
198 chi2min = chi2outn[i];
199 icouple = i;
200 }
201 }
202 /* Compute TRACK POINTS */
203 Computehitpsemes(number_of_hits, x, y, r, aref[icouple], bref[icouple]);
204 for (i = 0; i < number_of_hits; i++) {
205 xout[i] = m_xpoint[i];
206 yout[i] = m_ypoint[i];
207 }
208 /* Compute a & b parameters of the track from TRACK POINTS */
209 // int idim;
210 // ifail;
211 double aoutn, bout, sig2a, sig2b, corrab;
212 double temp, det;
213 double hesse[2][2];
214 double W, WX, WX2, WY, WXY;
215 // double errormatrix[2][2];
216 W = 0.;
217 WX = 0.;
218 WX2 = 0.;
219 WY = 0.;
220 WXY = 0;
221 for (int i = 0; i < number_of_hits; i++) {
222 temp = 1. / (sr[i] * sr[i]);
223 W += temp;
224 WX += xout[i] * temp;
225 WX2 += xout[i] * xout[i] * temp;
226 WY += yout[i] * temp;
227 WXY += xout[i] * yout[i] * temp;
228 }
229 det = W * WX2 - WX * WX;
230 aoutn = (W * WXY - WY * WX) / det;
231 bout = (WY * WX2 - WX * WXY) / det;
232 hesse[1][1] = W;
233 hesse[0][0] = WX2;
234 hesse[1][0] = WX;
235 hesse[0][1] = WX;
236 // idim = 2;
237 /* invert hessian matrix */
238 hesse[1][1] = hesse[0][0] / det;
239 hesse[0][0] = hesse[1][1] / det;
240 hesse[1][0] = -1. * (hesse[0][1] / det);
241 hesse[0][1] = -1. * (hesse[0][1] / det);
242 sig2a = hesse[0][0];
243 sig2b = hesse[1][1];
244 corrab = hesse[1][0];
245 double deno, btest, bs1, bs2, db1, db2;
246 btest = bout;
247 bout = 0;
248 deno = 0;
249 for (int i = 0; i < number_of_hits; i++) {
250 bs1 = (y[i] - aoutn * x[i] - r[i] * std::sqrt(1 + aoutn * aoutn));
251 bs2 = (y[i] - aoutn * x[i] + r[i] * std::sqrt(1 + aoutn * aoutn));
252 db1 = std::abs(bs1 - btest);
253 db2 = std::abs(bs2 - btest);
254 if (db1 < db2) {
255 bout += bs1 / (sr[i] * sr[i]);
256 } else {
257 bout += bs2 / (sr[i] * sr[i]);
258 }
259 deno += 1 / (sr[i] * sr[i]);
260 }
261 bout /= deno;
262 /* compute chi2 from residuals*/
263 double resd;
264 double chi2f = 0;
265 for (int i = 0; i < number_of_hits; i++) {
266 double xi, yi;
267 xi = (aoutn * y[i] - aoutn * bout + x[i]) / (1. + aoutn * aoutn);
268 yi = aoutn * xi + bout;
269 resd = (std::sqrt((xi - x[i]) * (xi - x[i]) + (yi - y[i]) * (yi - y[i])) - r[i]) / sr[i];
270 chi2f += resd * resd;
271 }
272 /* set global variables (method output)*/
273 m_aoutn = aoutn;
274 m_sig2a = sig2a;
275 m_bout = bout;
276 m_sig2b = sig2b;
277 m_corrab = corrab;
278 m_chi2f = chi2f;
279 }
#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
float j(const xAOD::IParticle &, const xAOD::TrackMeasurementValidation &hit, const Eigen::Matrix3d &jab_inv)

◆ get_a()

double MuonCalib::MuCCaFitterImplementation::get_a ( )

Definition at line 281 of file MuCCaFitter.cxx.

281{ return m_aoutn; }

◆ get_b()

double MuonCalib::MuCCaFitterImplementation::get_b ( )

Definition at line 282 of file MuCCaFitter.cxx.

282{ return m_bout; }

◆ get_chi2f()

double MuonCalib::MuCCaFitterImplementation::get_chi2f ( )

Definition at line 286 of file MuCCaFitter.cxx.

286{ return m_chi2f; }

◆ get_corrab()

double MuonCalib::MuCCaFitterImplementation::get_corrab ( )

Definition at line 285 of file MuCCaFitter.cxx.

285{ return m_corrab; }

◆ get_da()

double MuonCalib::MuCCaFitterImplementation::get_da ( )

Definition at line 283 of file MuCCaFitter.cxx.

283{ return m_sig2a; }

◆ get_db()

double MuonCalib::MuCCaFitterImplementation::get_db ( )

Definition at line 284 of file MuCCaFitter.cxx.

284{ 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: