Muon::Fit2D Class Reference

A 2D linear regression calculator. More...

#include <Fit2D.h>

Collaboration diagram for Muon::Fit2D:

Classes
struct	LinStats
	A structure to hold linear fit statistics. More...
struct	Point
	A 2D point used in statistics and fit. More...
struct	SimpleStats

Public Types
typedef std::vector< Point * >	PointArray
	A vector of points.

Public Member Functions
	Fit2D ()
	Constructor.

Static Public Member Functions
static void	SimpleStatistics (const PointArray &points, SimpleStats &stats)
	Calculate simple statistics for the Y values of a set of points.
static void	fitLine (PointArray &points, double fExclChi2, bool bDump, LinStats &stats)
	Fit a straight line through the given points.
static void	fitPoint (PointArray &points, double fExclChi2, bool bDump, SimpleStats &stats)
	Estimate a new point from the given points.

Detailed Description

A 2D linear regression calculator.

The class accepts a list of X/ points and computes the least-squares straight line fit. It reports the intercept, slope and goodness of fit.

Points can be excluded from the calculation, and can also have a weight associated with them.

The class also computes simple statistics (meand, standard deviation) of the Y values.

Definition at line 28 of file Fit2D.h.

Member Typedef Documentation

PointArray

typedef std::vector<Point*> Muon::Fit2D::PointArray

A vector of points.

Definition at line 59 of file Fit2D.h.

Constructor & Destructor Documentation

Fit2D()

Muon::Fit2D::Fit2D ( )

inline

Constructor.

Definition at line 179 of file Fit2D.h.

    {
    }

Member Function Documentation

fitLine()

void Muon::Fit2D::fitLine ( Fit2D::PointArray & points, double fExclChi2, bool bDump, LinStats & stats )

static

Fit a straight line through the given points.

Parameters

points	The list of data points
fExclChi2	CHi2 value for excluding outliers
bDump	Write details to log
stats	[output] The fit results

Definition at line 132 of file Fit2D.cxx.

{
  MsgStream log(Athena::getMessageSvc(),"Fit2D::fitLine(");
    double *x = new double[points.size()];
    double *y = new double[points.size()];
    double *w = new double[points.size()];
    for (;;)
    {
        stats.clear();
        for (unsigned iPt = 0; iPt < points.size(); iPt++)
        {
            if (!points[iPt]->bExclude)
            {
                x[stats.n] = points[iPt]->fX;
                y[stats.n] = points[iPt]->fY;
                w[stats.n] = points[iPt]->fW;
                stats.n++;
            }
        }
        if (stats.n < 2) goto done;
 
        gsl_fit_wlinear(x, 1,
                        w, 1,
                        y, 1,
                        stats.n,
                        &stats.fIntercept, &stats.fSlope,
                        &stats.fCov[0][0], &stats.fCov[0][1], &stats.fCov[1][1],
                        &stats.fChi2);
        stats.fCov[1][0] = stats.fCov[0][1];
        if (bDump && log.level()<=MSG::DEBUG)
      log << MSG::DEBUG << stats.toString() << endmsg;
 
        double fY, fYerr, fMaxChi2 = 0;
        int i = 0, iPtMax = 0;
        for (unsigned iPt = 0; iPt < points.size(); iPt++)
        {
            if (!points[iPt]->bExclude)
            {
                stats.eval(x[i], fY, fYerr);
                double fDiff = y[i] - fY;
                points[iPt]->fChi2 = w[i] * fDiff * fDiff;
                if (fMaxChi2 < points[iPt]->fChi2)
                {
                    fMaxChi2 = points[iPt]->fChi2;
                    iPtMax = iPt;
                }
                if (bDump && log.level()<=MSG::DEBUG)
                {
                    log << MSG::DEBUG << "idx=" << points[iPt]->nIdx
                    << ", x=" << x[i]
                    << ", y=" << y[i]
                    << ", w=" << w[i]
                    << ", Y=" << fY
                    << ", chi2=" << points[iPt]->fChi2
                    << endmsg;
                }
                i++;
            }
        }
        if (fMaxChi2 < fExclChi2 || stats.n <= 2)
            break;
        points[iPtMax]->bExclude = true;
    }
done:
    delete [] x;
    delete [] y;
    delete [] w;
}

fitPoint()

void Muon::Fit2D::fitPoint ( PointArray & points, double fExclChi2, bool bDump, SimpleStats & stats )

static

Estimate a new point from the given points.

Parameters

points	The list of data points
fExclChi2	CHi2 value for excluding outliers
bDump	Write details to log
stats	[output] The fit results

Definition at line 73 of file Fit2D.cxx.

{
    MsgStream log(Athena::getMessageSvc(),"Fit2D::fitPoint(");
    double *y = new double[points.size()];
    double *w = new double[points.size()];
    for (;;)
    {
        stats.clear();
        for (unsigned iPt = 0; iPt < points.size(); iPt++)
        {
            // Exclude outliers.
            if (!points[iPt]->bExclude)
            {
                y[stats.n] = points[iPt]->fY;
                w[stats.n] = points[iPt]->fW;
                stats.n++;
            }
        }
        // Calculate mean and standard deviation.
        if (stats.n == 0) goto done;
        stats.fMean = gsl_stats_wmean(w, 1, y, 1, stats.n);
        if (stats.n == 1) goto done;
        stats.fStd = gsl_stats_wsd(w, 1, y, 1, stats.n);
        double fMaxChi2 = 0;
        int i = 0, iPtMax = 0;
        for (unsigned iPt = 0; iPt < points.size(); iPt++)
        {
            if (!points[iPt]->bExclude)
            {
                double fDiff = y[i] - stats.fMean;
                points[iPt]->fChi2 = w[i] * fDiff * fDiff;
                stats.fChi2 += points[iPt]->fChi2;
                if (fMaxChi2 < points[iPt]->fChi2)
                {
                    fMaxChi2 = points[iPt]->fChi2;
                    iPtMax = iPt;
                }
                if (bDump && log.level()<=MSG::DEBUG)
                {
                    log << MSG::DEBUG  << "idx=" << points[iPt]->nIdx
                    << ", y=" << y[i]
                    << ", w=" << w[i]
                    << ", chi2=" << points[iPt]->fChi2
                    << endmsg;
                }
                i++;
            }
        }
        if (bDump && log.level()<=MSG::DEBUG)
      log << MSG::DEBUG << stats.toString() << endmsg;
        if (fMaxChi2 < fExclChi2 || stats.n <= 2)
            break;
        points[iPtMax]->bExclude = true;
    }
done:
    delete [] y;
    delete [] w;
}

SimpleStatistics()

void Muon::Fit2D::SimpleStatistics ( const PointArray & points, SimpleStats & stats )

static

Calculate simple statistics for the Y values of a set of points.

Parameters

points	The list of data points
stats	[output] The statistics results

Definition at line 43 of file Fit2D.cxx.

{
    stats.clear();
    double *y = new double[points.size()];
    double *w = new double[points.size()];
    for (unsigned iPt = 0; iPt < points.size(); iPt++)
    {
        // Exclude outliers.
        if (!points[iPt]->bExclude)
        {
            y[stats.n] = points[iPt]->fY;
            w[stats.n] = points[iPt]->fW;
            stats.n++;
        }
    }
    if (stats.n == 0) goto done;
    // Calculate mean and standard deviation.
    stats.fMean = gsl_stats_wmean(w, 1, y, 1, stats.n);
    if (stats.n == 1) goto done;
    stats.fStd  = gsl_stats_wsd  (w, 1, y, 1, stats.n);
    for (int i = 0; i < stats.n; i++)
    {
        double fDiff = y[i] - stats.fMean;
        stats.fChi2 += w[i] * fDiff * fDiff;
    }
done:
    delete [] y;
    delete [] w;
}

---

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

• Fit2D.h
• Fit2D.cxx