/
CbmRichRingFitterCircle.cxx
122 lines (91 loc) · 2.45 KB
/
CbmRichRingFitterCircle.cxx
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
/**
* \file CbmRichRingFitterCircle.cxx
*
* \author Supriya Das
* \date 2006
**/
#include "CbmRichRingFitterCircle.h"
#include "CbmRichRingLight.h"
CbmRichRingFitterCircle::CbmRichRingFitterCircle()
{
}
CbmRichRingFitterCircle::~CbmRichRingFitterCircle()
{
}
void CbmRichRingFitterCircle::DoFit(
CbmRichRingLight* ring)
{
int nofHits = ring->GetNofHits();
if (nofHits < 3) return;
float c[3], a[3][3];
float b1 = 0.f;
float b2 = 0.f;
float b3 = 0.f;
float b12 = 0.f;
float b22 = 0.f;
float b32 = nofHits;
float a11 = 0.f;
float a12 = 0.f;
float a13 = 0.f;
float a21 = 0.f;
float a22 = 0.f;
float a23 = 0.f;
float a31 = 0.f;
float a32 = 0.f;
float a33 = 0.f;
float meanX = 0.f;
float meanY = 0.f;
for (int iHit = 0; iHit < nofHits; iHit++) {
float hx = ring->GetHit(iHit).fX;
float hy = ring->GetHit(iHit).fY;
b1 += (hx*hx + hy*hy) * hx;
b2 += (hx*hx + hy*hy) * hy;
b3 += (hx*hx + hy*hy);
b12 += hx;
b22 += hy;
a11 += 2*hx*hx;
a12 += 2*hx*hy;
a22 += 2*hy*hy;
meanX += hx;
meanY += hy;
}
if (nofHits != 0) {
meanX = meanX/(float)(nofHits);
meanY = meanY/(float)(nofHits);
}
a21 = a12;
a13 = b12;
a23 = b22;
a31 = 2*b12;
a32 = 2*b22;
a33 = b32;
c[0] = b1*b22 - b2*b12;
c[1] = b1*b32 - b3*b12;
c[2] = b2*b32 - b3*b22;
a[0][0] = a11*b22 - a21*b12;
a[1][0] = a12*b22 - a22*b12;
a[2][0] = a13*b22 - a23*b12;
a[0][1] = a11*b32 - a31*b12;
a[1][1] = a12*b32 - a32*b12;
a[2][1] = a13*b32 - a33*b12;
a[0][2] = a21*b32-a31*b22;
a[1][2] = a22*b32-a32*b22;
a[2][2] = a23*b32-a33*b22;
float det1 = a[0][0]*a[1][1] - a[0][1]*a[1][0];
float x11 = (c[0]*a[1][1] - c[1]*a[1][0]) / det1;
float x21 = (a[0][0]*c[1] - a[0][1]*c[0]) / det1;
// Float_t det2 = a[0][1]*a[1][2] - a[0][2]*a[1][1];
// Float_t det3 = a[0][0]*a[1][2] - a[0][2]*a[1][0];
// Float_t x12 = (c[1]*a[1][2] - c[2]*a[1][1])/det2;
// Float_t x22 = (a[0][1]*c[2] - a[0][2]*c[1])/det2;
float radius = sqrt((b3 + b32*(x11*x11 + x21*x21) - a31*x11 - a32*x21)/a33);
float centerX = x11;
float centerY = x21;
ring->SetRadius(radius);
ring->SetCenterX(centerX);
ring->SetCenterY(centerY);
//if (TMath::IsNaN(radius) == 1) ring->SetRadius(0.);
//if (TMath::IsNaN(centerX) == 1) ring->SetCenterX(0.);
//if (TMath::IsNaN(centerY) == 1) ring->SetCenterY(0.);
CalcChi2(ring);
}