MuonCalib::MuCCaFitterImplementation Class Reference

#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 More...
double	m_bfpar [4]
double	m_xout0
	track points from 2 circles More...
double	m_yout0
double	m_xout1
double	m_yout1
double	m_xpoint [100]
	track points More...
double	m_ypoint [100]
double	m_aoutn
	results More...
double	m_sig2a
	track slope's variance More...
double	m_bout
	track constant term More...
double	m_sig2b
	error on track constant term More...
double	m_corrab
	correlation term More...
double	m_chi2f
	chisquared More...

Detailed Description

The straight line fitter for drift circles used by Calib

Author

Fabri.nosp@m.zio..nosp@m.Petru.nosp@m.cci@.nosp@m.cern..nosp@m.ch, Paolo.nosp@m..Bra.nosp@m.nchin.nosp@m.i@ce.nosp@m.rn.ch

Definition at line 67 of file MuCCaFitter.h.

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.

                                                                                                         {
        /************************/
        /* Compute TRACK POINTS */
        /************************/
        int i, nhitc;
        //    double xpoint[nhit],ypoint[nhit];
        /* loop over hit couples */
        nhitc = nhit / 2;
        for (i = 0; i < nhitc; i++) {
            /* Compute TRACK POINTS for 2 circles */
            computehitsfromcircles(xcirc[2 * i], ycirc[2 * i], rcirc[2 * i], xcirc[2 * i + 1], ycirc[2 * i + 1], rcirc[2 * i + 1], a, b);
            m_xpoint[2 * i] = m_xout0;
            m_ypoint[2 * i] = m_yout0;
            m_xpoint[2 * i + 1] = m_xout1;
            m_ypoint[2 * i + 1] = m_yout1;
        }
        if (2 * nhitc != nhit) {
            /* Compute TRACK POINTS for 2 circles */
            computehitsfromcircles(xcirc[nhit - 2], ycirc[nhit - 2], rcirc[nhit - 2], xcirc[nhit - 1], ycirc[nhit - 1], rcirc[nhit - 1], a,
                                   b);
            m_xpoint[nhit - 2] = m_xout0;
            m_ypoint[nhit - 2] = m_yout0;
            m_xpoint[nhit - 1] = m_xout1;
            m_ypoint[nhit - 1] = m_yout1;
        }
    }

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.

                                                                     {
        /**************************************/
        /* Compute TRACK POINTS for 2 circles */
        /**************************************/
        double ang[4], bpar[4];
        int i;
        double xout0 = 0, yout0 = 0;
        double xout1 = 0, yout1 = 0;
        double dist, dist1, xint, yint, xref, yref;
 
        /* Compute the 4 lines tangent to 2 circles */
        Computelin(x0, y0, r0, x1, y1, r1);
        for (i = 0; i < 4; i++) {
            bpar[i] = m_bfpar[i];
            ang[i] = m_angularcoefficient[i];
        }
 
        dist = 9999999;
        /* for each of the 4 tangents : */
        for (i = 0; i < 4; i++) {
            /* compute tangent point */
            xint = (x0 + ang[i] * y0 - ang[i] * bpar[i]) / (1 + ang[i] * ang[i]);
            yint = ang[i] * (xint) + bpar[i];
            /* compute point from reference line */
            double bprime;
            if (a != 0) {
                bprime = y0 + x0 / a;
                xref = (bprime - b) * a / (a * a + 1);
                yref = (b + bprime * a * a) / (a * a + 1);
            } else {
                xref = x0;
                yref = b;
            }
            /* choose tangent point closest to that of the reference lines (TRACK POINT)*/
            dist1 = std::sqrt((xint - xref) * (xint - xref) + (yint - yref) * (yint - yref));
            if (dist1 < dist) {
                dist = dist1;
                xout0 = xint;
                yout0 = yint;
            }
        }
        dist = 9999999;
        for (i = 0; i < 4; i++) {
            /* compute tangent point */
            xint = (x1 + ang[i] * y1 - ang[i] * bpar[i]) / (1 + ang[i] * ang[i]);
            yint = ang[i] * (xint) + bpar[i];
            /* compute point from reference line */
            double bprime;
            if (a != 0) {
                bprime = y1 + x1 / a;
                xref = (bprime - b) * a / (a * a + 1);
                yref = (b + bprime * a * a) / (a * a + 1);
            } else {
                xref = x1;
                yref = b;
            }
            /* choose tangent point closest to that of the reference line (TRACK POINT)*/
            dist1 = std::sqrt((xint - xref) * (xint - xref) + (yint - yref) * (yint - yref));
            if (dist1 < dist) {
                dist = dist1;
                xout1 = xint;
                yout1 = yint;
            }
        }
        /* set global variables (method output)*/
        m_xout0 = xout0;
        m_yout0 = yout0;
        m_xout1 = xout1;
        m_yout1 = yout1;
    }

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.

                                                                                                               {
        /********************************************/
        /* Compute the 4 lines tangent to 2 circles */
        /********************************************/
        double delta{0}, averagephi{0}, dist{0}, dphiin{0}, dphiex{0};
        double bpar[4], phi{0};
        double bfparn[2];
        bfparn[0] = 0;
        bfparn[1] = 0;
        double bcand[2][4];
        int segnob[] = {1, -1, 1, -1};
        int ncandid[4], ncand;
        int i{0}, firsttime{0};
        double angularcoefficient[4];
        double bfpar[4];
        double dy{0}, dx{0};
 
        /* compute a parameters */
        dx = x2 - x1;
        dy = y2 - y1;
        averagephi = std::atan2(dy, dx);
        dist = std::hypot(dx, dy);
 
        delta = r2 + r1;
        dphiin = std::asin(delta / dist);
 
        delta = r2 - r1;
        dphiex = std::asin(delta / dist);
 
        int f = 1;
        phi = averagephi + f * dphiex;
        if (phi < 0) { phi = 2 * M_PI + (phi); }
        angularcoefficient[0] = std::tan(phi);
 
        f = -1;
        phi = averagephi + f * dphiex;
        if (phi < 0) { phi = 2 * M_PI + (phi); }
        angularcoefficient[1] = std::tan(phi);
 
        f = 1;
        phi = averagephi + f * dphiin;
        if (phi < 0) { phi = 2 * M_PI + (phi); }
        angularcoefficient[2] = std::tan(phi);
 
        f = -1;
        phi = averagephi + f * dphiin;
        if (phi < 0) { phi = 2 * M_PI + (phi); }
        angularcoefficient[3] = std::tan(phi);
 
        /* Compute b parameters */
        for (i = 0; i < 4; i++) {
            bpar[0] = y1 - angularcoefficient[i] * x1 + segnob[0] * r1 * std::sqrt(1 + angularcoefficient[i] * angularcoefficient[i]);
            bpar[1] = y1 - angularcoefficient[i] * x1 + segnob[1] * r1 * std::sqrt(1 + angularcoefficient[i] * angularcoefficient[i]);
            bpar[2] = y2 - angularcoefficient[i] * x2 + segnob[2] * r2 * std::sqrt(1 + angularcoefficient[i] * angularcoefficient[i]);
            bpar[3] = y2 - angularcoefficient[i] * x2 + segnob[3] * r2 * std::sqrt(1 + angularcoefficient[i] * angularcoefficient[i]);
            double delta = 0.00001;
            ncand = 0;
            if (std::abs(bpar[0] - bpar[2]) < delta) {
                bfparn[ncand] = bpar[0];
                ncand = ncand + 1;
            }
            if (std::abs(bpar[0] - bpar[3]) < delta) {
                bfparn[ncand] = bpar[0];
                ncand = ncand + 1;
            }
            if (std::abs(bpar[1] - bpar[2]) < delta) {
                bfparn[ncand] = bpar[1];
                ncand = ncand + 1;
            }
            if (std::abs(bpar[1] - bpar[3]) < delta) {
                bfparn[ncand] = bpar[1];
                ncand = ncand + 1;
            }
            bcand[0][i] = bfparn[0];
            bcand[1][i] = bfparn[1];
            ncandid[i] = ncand;
        }
        firsttime = 0;
        for (i = 0; i < 4; i++) {
            if (ncandid[i] == 1) {
                bfpar[i] = bcand[0][i];
            } else {
                bfpar[i] = bcand[firsttime][i];
                firsttime++;
            }
        }
        /* set global variables (method output)*/
        for (int i = 0; i < 4; i++) {
            m_bfpar[i] = bfpar[i];
            m_angularcoefficient[i] = angularcoefficient[i];
        }
    }

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.

                                                                                                                     {
        /********************************************/
        /* Compute the 4 lines tangent to 2 circles */
        /********************************************/
        double delta;
        double averagephi, dist, dphiin, dphiex;
        double bpar[4], phi;
        double bfparn[2];
        bfparn[0] = 0;
        bfparn[1] = 0;
        // int segnobfn[2];
        // int segnoc[2][4],ncandid[4],ncand;
        int ncandid[4], ncand;
        double bcand[2][4];
        int segnob[] = {1, -1, 1, -1};
        int firsttime, i;
        double angularcoefficient[4];
        double bfpar[4];
        // int segnobf[4];
        double dy, dx;
 
        /* compute a parameters */
        dx = x2 - x1;
        dy = y2 - y1;
        averagephi = std::atan2(dy, dx);
        dist = std::sqrt((x2 - x1) * (x2 - x1) + (y2 - y1) * (y2 - y1));
        delta = r2 + r1;
 
        dphiin = std::asin(delta / dist);
 
        delta = r2 - r1;
        dphiex = std::asin(delta / dist);
 
        int f = 1;
        phi = averagephi + f * dphiex;
        if (phi < 0) { phi = 2 * M_PI + (phi); }
        angularcoefficient[0] = std::tan(phi);
 
        f = -1;
        phi = averagephi + f * dphiex;
        if (phi < 0) { phi = 2 * M_PI + (phi); }
        angularcoefficient[1] = std::tan(phi);
 
        f = 1;
        phi = averagephi + f * dphiin;
        if (phi < 0) { phi = 2 * M_PI + (phi); }
        angularcoefficient[2] = std::tan(phi);
 
        f = -1;
        phi = averagephi + f * dphiin;
        if (phi < 0) { phi = 2 * M_PI + (phi); }
        angularcoefficient[3] = std::tan(phi);
 
        /* Compute b parameters */
        for (i = 0; i < 4; i++) {
            bpar[0] = y1 - angularcoefficient[i] * x1 + segnob[0] * r1 * std::sqrt(1 + angularcoefficient[i] * angularcoefficient[i]);
            bpar[1] = y1 - angularcoefficient[i] * x1 + segnob[1] * r1 * std::sqrt(1 + angularcoefficient[i] * angularcoefficient[i]);
            bpar[2] = y2 - angularcoefficient[i] * x2 + segnob[2] * r2 * std::sqrt(1 + angularcoefficient[i] * angularcoefficient[i]);
            bpar[3] = y2 - angularcoefficient[i] * x2 + segnob[3] * r2 * std::sqrt(1 + angularcoefficient[i] * angularcoefficient[i]);
            double delta = 0.00001;
            ncand = 0;
            if (std::abs(bpar[0] - bpar[2]) < delta) {
                bfparn[ncand] = bpar[0];
                ncand = ncand + 1;
            }
            if (std::abs(bpar[0] - bpar[3]) < delta) {
                bfparn[ncand] = bpar[0];
                ncand = ncand + 1;
            }
            if (std::abs(bpar[1] - bpar[2]) < delta) {
                bfparn[ncand] = bpar[1];
                ncand = ncand + 1;
            }
            if (std::abs(bpar[1] - bpar[3]) < delta) {
                bfparn[ncand] = bpar[1];
                ncand = ncand + 1;
            }
            bcand[0][i] = bfparn[0];
            bcand[1][i] = bfparn[1];
            ncandid[i] = ncand;
        }
        firsttime = 0;
        for (i = 0; i < 4; i++) {
            if (ncandid[i] == 1) {
                bfpar[i] = bcand[0][i];
            } else {
                bfpar[i] = bcand[firsttime][i];
                firsttime++;
            }
        }
        /* set global variables (method output)*/
        for (i = 0; i < 4; i++) {
            m_bfpar[i] = bfpar[i];
            m_angularcoefficient[i] = angularcoefficient[i];
        }
    }

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.

                                                                               {
        /***************************************/
        /* Fit a line to n=number_of_hits hits */
        /***************************************/
        double xout[100], yout[100];
        double dist[100], rsub[100];
        double chi2outn[4], chi2min;
        double aref[4], bref[4];
        int fhit, i, j, icouple = 0, lhit;
        /* compute 4 tanget lines from first and last hit */
        fhit = 0;
        lhit = number_of_hits - 1;
        Computelinparnew(x[fhit], y[fhit], r[fhit], x[lhit], y[lhit], r[lhit]);
        for (int ii = 0; ii < 4; ii++) {
            bref[ii] = m_bfpar[ii];
            aref[ii] = m_angularcoefficient[ii];
        }
        /* choose the REFERENCE LINE (aref[icouple],bref[icouple])*/
        for (i = 0; i < 4; i++) {
            for (j = 0; j < number_of_hits; j++) {
                double d;
                d = std::abs(aref[i] * x[j] + bref[i] - y[j]);
                dist[j] = d / std::sqrt(1. + aref[i] * aref[i]);
                rsub[j] = r[j] - dist[j];
            }
            double chi2out = 0;
            for (int ii = 1; ii < number_of_hits - 1; ii++) { chi2out += rsub[ii] * rsub[ii] / (sr[ii] * sr[ii]); }
            chi2outn[i] = chi2out;
        }
        chi2min = 999999999.;
        for (i = 0; i < 4; i++) {
            if (chi2outn[i] < chi2min) {
                chi2min = chi2outn[i];
                icouple = i;
            }
        }
        /* Compute TRACK POINTS */
        Computehitpsemes(number_of_hits, x, y, r, aref[icouple], bref[icouple]);
        for (i = 0; i < number_of_hits; i++) {
            xout[i] = m_xpoint[i];
            yout[i] = m_ypoint[i];
        }
        /* Compute a & b parameters of the track from TRACK POINTS */
        // int idim;
        //  ifail;
        double aoutn, bout, sig2a, sig2b, corrab;
        double temp, det;
        double hesse[2][2];
        double W, WX, WX2, WY, WXY;
        //  double errormatrix[2][2];
        W = 0.;
        WX = 0.;
        WX2 = 0.;
        WY = 0.;
        WXY = 0;
        for (int i = 0; i < number_of_hits; i++) {
            temp = 1. / (sr[i] * sr[i]);
            W += temp;
            WX += xout[i] * temp;
            WX2 += xout[i] * xout[i] * temp;
            WY += yout[i] * temp;
            WXY += xout[i] * yout[i] * temp;
        }
        det = W * WX2 - WX * WX;
        aoutn = (W * WXY - WY * WX) / det;
        bout = (WY * WX2 - WX * WXY) / det;
        hesse[1][1] = W;
        hesse[0][0] = WX2;
        hesse[1][0] = WX;
        hesse[0][1] = WX;
        // idim        = 2;
        /* invert hessian matrix */
        hesse[1][1] = hesse[0][0] / det;
        hesse[0][0] = hesse[1][1] / det;
        hesse[1][0] = -1. * (hesse[0][1] / det);
        hesse[0][1] = -1. * (hesse[0][1] / det);
        sig2a = hesse[0][0];
        sig2b = hesse[1][1];
        corrab = hesse[1][0];
        double deno, btest, bs1, bs2, db1, db2;
        btest = bout;
        bout = 0;
        deno = 0;
        for (int i = 0; i < number_of_hits; i++) {
            bs1 = (y[i] - aoutn * x[i] - r[i] * std::sqrt(1 + aoutn * aoutn));
            bs2 = (y[i] - aoutn * x[i] + r[i] * std::sqrt(1 + aoutn * aoutn));
            db1 = std::abs(bs1 - btest);
            db2 = std::abs(bs2 - btest);
            if (db1 < db2) {
                bout += bs1 / (sr[i] * sr[i]);
            } else {
                bout += bs2 / (sr[i] * sr[i]);
            }
            deno += 1 / (sr[i] * sr[i]);
        }
        bout /= deno;
        /* compute chi2 from residuals*/
        double resd;
        double chi2f = 0;
        for (int i = 0; i < number_of_hits; i++) {
            double xi, yi;
            xi = (aoutn * y[i] - aoutn * bout + x[i]) / (1. + aoutn * aoutn);
            yi = aoutn * xi + bout;
            resd = (std::sqrt((xi - x[i]) * (xi - x[i]) + (yi - y[i]) * (yi - y[i])) - r[i]) / sr[i];
            chi2f += resd * resd;
        }
        /* set global variables (method output)*/
        m_aoutn = aoutn;
        m_sig2a = sig2a;
        m_bout = bout;
        m_sig2b = sig2b;
        m_corrab = corrab;
        m_chi2f = chi2f;
    }

get_a()

double MuonCalib::MuCCaFitterImplementation::get_a

(

)

Definition at line 271 of file MuCCaFitter.cxx.

{ return m_aoutn; }

get_b()

double MuonCalib::MuCCaFitterImplementation::get_b

(

)

Definition at line 272 of file MuCCaFitter.cxx.

{ return m_bout; }

get_chi2f()

double MuonCalib::MuCCaFitterImplementation::get_chi2f

(

)

Definition at line 276 of file MuCCaFitter.cxx.

{ return m_chi2f; }

get_corrab()

double MuonCalib::MuCCaFitterImplementation::get_corrab

(

)

Definition at line 275 of file MuCCaFitter.cxx.

{ return m_corrab; }

get_da()

double MuonCalib::MuCCaFitterImplementation::get_da

(

)

Definition at line 273 of file MuCCaFitter.cxx.

{ return m_sig2a; }

get_db()

double MuonCalib::MuCCaFitterImplementation::get_db

(

)

Definition at line 274 of file MuCCaFitter.cxx.

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

• MuCCaFitter.h
• MuCCaFitter.cxx