StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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();
44  mRadius = 0;
45  mXCenter = 0;
46  mYCenter = 0;
47  mVariance = 0;
48  mRC = 0;
49 }
50 
51 unsigned int StFastCircleFitter::numberOfPoints() const {return mX.size();}
52 
53 double StFastCircleFitter::radius() const {return mRadius;}
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;
79  double radius2 = 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 
144  radius2 = xd*xd+yd*yd+g1;
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);
151  dradius2 = dx*dx+dy*dy;
152  mVariance += dradius2+radius2-2.*::sqrt(dradius2*radius2);
153  }
154  mVariance /= npoints-3.0;
155 
156  mRadius = ::sqrt(radius2);
157  mRC = 0;
158 
159  return true;
160 }
Generated by   doxygen 1.8.5