StRoot  1
StFastCircleFitter.cc
1 /***************************************************************************
2  *
3  * \$Id: StFastCircleFitter.cc,v 1.2 2003/09/02 17:59:34 perev Exp \$
4  *
5  * Author: Thomas Ullrich, Dec 1999
6  ***************************************************************************
7  *
8  * Description:
9  *
10  * Fast fitting routine using a iterational linear regression
11  * method (ILRM). Reference: N.Chernov, G.A.Ososkov, Computer
12  * Physics Communication 33 (1984) 329-333.
13  *
14  ***************************************************************************
15  *
16  * \$Log: StFastCircleFitter.cc,v \$
17  * Revision 1.2 2003/09/02 17:59:34 perev
18  * gcc 3.2 updates + WarnOff
19  *
20  * Revision 1.1 1999/12/21 16:28:48 ullrich
21  * Initial Revision
22  *
23  **************************************************************************/
24 #include "StFastCircleFitter.hh"
25 #include <math.h>
26
27 StFastCircleFitter::StFastCircleFitter()
28 {
29  clear();
30 }
31
32 StFastCircleFitter::~StFastCircleFitter() {/* nop */}
33
34 void StFastCircleFitter::addPoint(double x, double y)
35 {
36  mX.push_back(x);
37  mY.push_back(y);
38 }
39
40 void StFastCircleFitter::clear()
41 {
42  mX.clear();
43  mY.clear();
45  mXCenter = 0;
46  mYCenter = 0;
47  mVariance = 0;
48  mRC = 0;
49 }
50
51 unsigned int StFastCircleFitter::numberOfPoints() const {return mX.size();}
52
54
55 double StFastCircleFitter::xcenter() const {return mXCenter;}
56
57 double StFastCircleFitter::ycenter() const {return mYCenter;}
58
59 double StFastCircleFitter::variance() const {return mVariance;}
60
61 int StFastCircleFitter::rc() const {return mRC;}
62
63 bool StFastCircleFitter::fit()
64 {
65  int i;
66  double xx, yy, xx2, yy2;
67  double f, g, h, p, q, t, g0, g02, a, b, c, d;
68  double xroot, ff, fp, xd, yd, g1;
69  double dx, dy, dradius2, xnom;
70
71  double xgravity = 0.0;
72  double ygravity = 0.0;
73  double x2 = 0.0;
74  double y2 = 0.0;
75  double xy = 0.0;
76  double xx2y2 = 0.0;
77  double yx2y2 = 0.0;
78  double x2y22 = 0.0;
80
81  mRC = 0;
82
83  int npoints = mX.size();
84
85  if (npoints <= 3) {
86  mRC = 1;
87  return false;
88  }
89  for (i=0; i<npoints; i++) {
90  xgravity += mX[i];
91  ygravity += mY[i];
92  }
93  xgravity /= npoints;
94  ygravity /= npoints;
95
96  for (i=0; i<npoints; i++) {
97  xx = mX[i]-xgravity;
98  yy = mY[i]-ygravity;
99  xx2 = xx*xx;
100  yy2 = yy*yy;
101  x2 += xx2;
102  y2 += yy2;
103  xy += xx*yy;
104  xx2y2 += xx*(xx2+yy2);
105  yx2y2 += yy*(xx2+yy2);
106  x2y22 += (xx2+yy2)*(xx2+yy2);
107  }
108  if (xy == 0.) {
109  mRC = 2;
110  return false;
111  }
112  f = (3.*x2+y2)/npoints;
113  g = (x2+3.*y2)/npoints;
114  h = 2*xy/npoints;
115  p = xx2y2/npoints;
116  q = yx2y2/npoints;
117  t = x2y22/npoints;
118  g0 = (x2+y2)/npoints;
119  g02 = g0*g0;
120  a = -4.0;
121  b = (f*g-t-h*h)/g02;
122  c = (t*(f+g)-2.*(p*p+q*q))/(g02*g0);
123  d = (t*(h*h-f*g)+2.*(p*p*g+q*q*f)-4.*p*q*h)/(g02*g02);
124  xroot = 1.0;
125  for (i=0; i<5; i++) {
126  ff = (((xroot+a)*xroot+b)*xroot+c)*xroot+d;
127  fp = ((4.*xroot+3.*a)*xroot+2.*b)*xroot+c;
128  xroot -= ff/fp;
129  }
130  g1 = xroot*g0;
131  xnom = (g-g1)*(f-g1)-h*h;
132  if (xnom == 0.) {
133  mRC = 3;
134  return false;
135  }
136  yd = (q*(f-g1)-h*p)/xnom;
137  xnom = f-g1;
138  if (xnom == 0.) {
139  mRC = 4;
140  return false;
141  }
142  xd = (p-h*yd )/xnom;
143
145  mXCenter = xd+xgravity;
146  mYCenter = yd+ygravity;
147
148  for (i=0; i<npoints; i++) {
149  dx = mX[i]-(mXCenter);
150  dy = mY[i]-(mYCenter);