00001
00002
00003 #include "MdcTrkRecon/MdcLine.h"
00004 #include <stdlib.h>
00005 #include <math.h>
00006
00007 #include <iostream>
00008 using namespace std;
00009
00010 int MdcLine::fit(int nUse) {
00011
00012 double s1, sx, sy, sxx, sxy;
00013 double delinv, denom;
00014 double weight = 0.013;
00015 int ihit;
00016 if (nUse > nPoint) nUse = nPoint;
00017 s1 = sx = sy = sxx = sxy = 0.0;
00018 chisq = 0.0;
00019 for (ihit = 0; ihit < nUse; ihit++) {
00020
00021 if (sigma[ihit] < 0.0) continue;
00022 weight = 1. / (sigma[ihit] * sigma[ihit]);
00023 s1 += weight;
00024 sx += x[ihit] * weight;
00025 sy += y[ihit] * weight;
00026 sxx += x[ihit] * (x[ihit] * weight);
00027 sxy += y[ihit] * (x[ihit] * weight);
00028
00029
00030
00031
00032
00033
00034 }
00035
00036
00037 denom = s1 * sxx - sx * sx;
00038 delinv = (denom == 0.0) ? 1.e20 : 1. / denom;
00039 intercept = (sy * sxx - sx * sxy) * delinv;
00040 slope = (s1 * sxy - sx * sy) * delinv;
00041 errmat[0] = sxx * delinv;
00042 errmat[1] = -sx * delinv;
00043 errmat[2] = s1 * delinv;
00044
00045
00046 for (ihit = 0; ihit < nUse; ihit++) {
00047 if (sigma[ihit] < 0.0) continue;
00048 resid[ihit] = ( y[ihit] - intercept - slope * x[ihit] );
00049 chisq += resid[ihit] * resid[ihit] / (sigma[ihit] * sigma[ihit]);
00050
00051 }
00052
00053
00054
00055
00056 return 0;
00057 }
00058
00059
00060
00061
00062