CbmRoot
CbmRichRingFitterTAU.cxx
Go to the documentation of this file.
1 
7 #include "CbmRichRingFitterTAU.h"
8 
9 #include <cmath>
10 #include <iostream>
11 #include <vector>
12 
13 using std::cout;
14 using std::endl;
15 using std::fabs;
16 using std::pow;
17 using std::sqrt;
18 using std::vector;
19 
21 
23 
25  int nofHits = ring->GetNofHits();
26 
27  double radius = -1.;
28  double centerX = -1.;
29  double centerY = -1.;
30 
31  if (nofHits < 3) {
32  ring->SetRadius(0.);
33  ring->SetCenterX(0.);
34  ring->SetCenterY(0.);
35  return;
36  }
37 
38  int iter, iterMax = 20;
39 
40  int riterMax = 4;
41 
42  double Xi, Yi, Zi;
43  double M0, Mx, My, Mz, Mxy, Mxx, Myy, Mxz, Myz, Mzz, Mxz2, Myz2, Cov_xy;
44  double A0, A1, A2, A22, A3, A33, epsilon = 0.000000000001;
45  double Dy, xnew, xold, ynew, yold = 100000000000.;
46  double GAM, DET;
47  //double Xcenter = -1.,Ycenter = -1.,Radius = -1.;
48 
49  double sigma;
50  double ctsigma;
51  double dist, weight;
52  const double ct = 3.;
53  const double copt = 0.2;
54  const double zero = 0.0001;
55  const double SigmaMin = 0.18;
56  double amax = -100000;
57  double amin = 100000;
58  double sig1 = 0.;
59  double sig2 = 0.;
60 
61  vector<double> x; // x coordinates of hits;
62  vector<double> y; // y coordinates of hits;
63  vector<double> a; // amplitudes of hits;
64  vector<double> w; // weights of points;
65  vector<double> d; // distance between points and circle;
66 
67  if (fRobust < 1 || fRobust > 2) { riterMax = 1; }
68  if (fRobust == 1) { riterMax = 4; }
69  if (fRobust == 2) { riterMax = 4; }
70 
71  for (int i = 0; i < nofHits; i++) {
72  x.push_back(ring->GetHit(i).fX);
73  y.push_back(ring->GetHit(i).fY);
74  a.push_back(1.);
75  }
76 
77  for (int i = 0; i < nofHits; i++) {
78  if (a[i] > amax) amax = a[i];
79  if (a[i] < amin) amin = a[i];
80  }
81 
82  for (int i = 0; i < nofHits; i++)
83  w.push_back(1.);
84 
85  for (int riter = 0; riter < riterMax; riter++) {
86  M0 = 0;
87  Mx = My = 0.;
88 
89  for (int i = 0; i < nofHits; i++) {
90  Mx += x[i] * w[i];
91  My += y[i] * w[i];
92  M0 += w[i];
93  }
94 
95  Mx /= M0;
96  My /= M0;
97 
98  //computing moments (note: all moments are normed, i.e. divided by N)
99  Mxx = Myy = Mxy = Mxz = Myz = Mzz = 0.;
100 
101  for (int i = 0; i < nofHits; i++) {
102  Xi = x[i] - Mx;
103  Yi = y[i] - My;
104  Zi = Xi * Xi + Yi * Yi;
105 
106  Mxy += Xi * Yi * w[i];
107  Mxx += Xi * Xi * w[i];
108  Myy += Yi * Yi * w[i];
109  Mxz += Xi * Zi * w[i];
110  Myz += Yi * Zi * w[i];
111  Mzz += Zi * Zi * w[i];
112  }
113  Mxx /= M0;
114  Myy /= M0;
115  Mxy /= M0;
116  Mxz /= M0;
117  Myz /= M0;
118  Mzz /= M0;
119 
120  // computing the coefficients of the characteristic polynomial
121  Mz = Mxx + Myy;
122  Cov_xy = Mxx * Myy - Mxy * Mxy;
123  Mxz2 = Mxz * Mxz;
124  Myz2 = Myz * Myz;
125 
126  A3 = 4. * Mz;
127  A2 = -3. * Mz * Mz - Mzz;
128  A1 = Mzz * Mz + 4. * Cov_xy * Mz - Mxz2 - Myz2 - Mz * Mz * Mz;
129  A0 = Mxz2 * Myy + Myz2 * Mxx - Mzz * Cov_xy - 2. * Mxz * Myz * Mxy
130  + Mz * Mz * Cov_xy;
131 
132  A22 = A2 + A2;
133  A33 = A3 + A3 + A3;
134  iter = 0;
135  xnew = 0.;
136 
137  // Newton's method starting at x=0
138  for (iter = 0; iter < iterMax; iter++) {
139  ynew = A0 + xnew * (A1 + xnew * (A2 + xnew * A3));
140 
141  if (fabs(ynew) > fabs(yold)) {
142  xnew = 0.;
143  break;
144  }
145 
146  Dy = A1 + xnew * (A22 + xnew * A33);
147  xold = xnew;
148  xnew = xold - ynew / Dy;
149 
150  if (fabs((xnew - xold) / xnew) < epsilon) break;
151  }
152 
153  if (iter == iterMax - 1) {
154  //printf("Newton3 does not converge in %d iterations\n",iterMax);
155  xnew = 0.;
156  }
157  if (xnew < 0.) {
158  iter = 30;
159  // printf("Negative root: x=%f\n",xnew);
160  }
161 
162  // computing the circle parameters
163  GAM = -Mz;
164  DET = xnew * xnew - xnew * Mz + Cov_xy;
165  centerX = (Mxz * (Myy - xnew) - Myz * Mxy) / DET / 2.;
166  centerY = (Myz * (Mxx - xnew) - Mxz * Mxy) / DET / 2.;
167  radius = sqrt(centerX * centerX + centerY * centerY - GAM);
168  centerX = (Mxz * (Myy - xnew) - Myz * Mxy) / DET / 2. + Mx;
169  centerY = (Myz * (Mxx - xnew) - Mxz * Mxy) / DET / 2. + My;
170 
171  if (riter < riterMax - 1) {
172  for (int i = 0; i < nofHits; i++) {
173  dist =
174  sqrt(pow((centerX - x[i]), 2) + pow((centerY - y[i]), 2)) - radius;
175  dist = fabs(dist);
176  d.push_back(dist);
177  }
178 
179  for (unsigned int i = 0; i < d.size(); i++) {
180  sig1 += w[i] * d[i] * d[i];
181  sig2 += w[i];
182  }
183  sigma = sqrt(sig1 / sig2);
184  if (sigma < SigmaMin) sigma = SigmaMin;
185  ctsigma = ct * sigma;
186  sig1 = 0.;
187  sig2 = 0.;
188 
189  w.clear();
190  for (unsigned int i = 0; i < d.size(); i++) {
191  if (fRobust == 1) {
192  if (d[i] <= ctsigma) {
193  weight = pow((1 - pow((d[i] / ctsigma), 2)), 2);
194  if (weight < zero) weight = zero;
195  } else {
196  weight = zero;
197  }
198  w.push_back(weight);
199  }
200  if (fRobust == 2 || fRobust == 3) {
201  sigma *= 2;
202  weight = 1 + copt;
203  weight /= 1 + copt * exp(pow((d[i] / sigma), 2));
204  if (weight < zero) weight = zero;
205  w.push_back(weight);
206  }
207  }
208  d.clear();
209  }
210  }
211  x.clear();
212  y.clear();
213  a.clear();
214  w.clear();
215 
216  ring->SetRadius(radius);
217  ring->SetCenterX(centerX);
218  ring->SetCenterY(centerY);
219 
220  CalcChi2(ring);
221 }
exp
friend F32vec4 exp(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:134
CbmRichRingFitterTAU::fRobust
int fRobust
Definition: CbmRichRingFitterTAU.h:42
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmRichRingFitterTAU::~CbmRichRingFitterTAU
virtual ~CbmRichRingFitterTAU()
Destructor.
Definition: CbmRichRingFitterTAU.cxx:22
CbmRichHitLight::fX
float fX
Definition: CbmRichRingLight.h:34
CbmRichRingLight::SetRadius
void SetRadius(float r)
Definition: CbmRichRingLight.h:124
CbmRichRingFitterBase::CalcChi2
virtual void CalcChi2(CbmRichRingLight *ring)
Calculate chi2 for circle fitting algorithms.
Definition: CbmRichRingFitterBase.h:48
CbmRichRingLight::GetNofHits
int GetNofHits() const
Return number of hits in ring.
Definition: CbmRichRingLight.h:108
CbmRichRingLight::SetCenterY
void SetCenterY(float y)
Definition: CbmRichRingLight.h:123
CbmRichRingLight::SetCenterX
void SetCenterX(float x)
Definition: CbmRichRingLight.h:122
d
double d
Definition: P4_F64vec2.h:24
CbmRichRingFitterTAU::DoFit
virtual void DoFit(CbmRichRingLight *ring)
Inherited from CbmRichRingFitterBase.
Definition: CbmRichRingFitterTAU.cxx:24
CbmRichHitLight::fY
float fY
Definition: CbmRichRingLight.h:35
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmRichRingLight::GetHit
CbmRichHitLight GetHit(int ind)
Return hit by the index.
Definition: CbmRichRingLight.h:114
CbmRichRingFitterTAU::CbmRichRingFitterTAU
CbmRichRingFitterTAU()
Default constructor.
Definition: CbmRichRingFitterTAU.cxx:20
CbmRichRingLight
Definition: CbmRichRingLight.h:39
CbmRichRingFitterTAU.h
Here the ring is fitted with the TAU algorithm from A. Ayriyan/ G. Ososkov.