CbmRoot
CbmRichRingFitterCircle.cxx
Go to the documentation of this file.
1 
9 #include "CbmRichRingLight.h"
10 
12 
14 
16  int nofHits = ring->GetNofHits();
17  if (nofHits < 3) return;
18 
19  float c[3], a[3][3];
20 
21  float b1 = 0.f;
22  float b2 = 0.f;
23  float b3 = 0.f;
24 
25  float b12 = 0.f;
26  float b22 = 0.f;
27  float b32 = nofHits;
28 
29  float a11 = 0.f;
30  float a12 = 0.f;
31  float a13 = 0.f;
32  float a21 = 0.f;
33  float a22 = 0.f;
34  float a23 = 0.f;
35  float a31 = 0.f;
36  float a32 = 0.f;
37  float a33 = 0.f;
38 
39  float meanX = 0.f;
40  float meanY = 0.f;
41 
42  for (int iHit = 0; iHit < nofHits; iHit++) {
43  float hx = ring->GetHit(iHit).fX;
44  float hy = ring->GetHit(iHit).fY;
45 
46  b1 += (hx * hx + hy * hy) * hx;
47  b2 += (hx * hx + hy * hy) * hy;
48  b3 += (hx * hx + hy * hy);
49 
50  b12 += hx;
51  b22 += hy;
52 
53  a11 += 2 * hx * hx;
54  a12 += 2 * hx * hy;
55  a22 += 2 * hy * hy;
56 
57  meanX += hx;
58  meanY += hy;
59  }
60 
61  if (nofHits != 0) {
62  meanX = meanX / (float) (nofHits);
63  meanY = meanY / (float) (nofHits);
64  }
65 
66  a21 = a12;
67 
68  a13 = b12;
69  a23 = b22;
70 
71  a31 = 2 * b12;
72  a32 = 2 * b22;
73  a33 = b32;
74 
75  c[0] = b1 * b22 - b2 * b12;
76  c[1] = b1 * b32 - b3 * b12;
77  c[2] = b2 * b32 - b3 * b22;
78 
79  a[0][0] = a11 * b22 - a21 * b12;
80  a[1][0] = a12 * b22 - a22 * b12;
81  a[2][0] = a13 * b22 - a23 * b12;
82 
83  a[0][1] = a11 * b32 - a31 * b12;
84  a[1][1] = a12 * b32 - a32 * b12;
85  a[2][1] = a13 * b32 - a33 * b12;
86 
87  a[0][2] = a21 * b32 - a31 * b22;
88  a[1][2] = a22 * b32 - a32 * b22;
89  a[2][2] = a23 * b32 - a33 * b22;
90 
91  float det1 = a[0][0] * a[1][1] - a[0][1] * a[1][0];
92 
93  float x11 = (c[0] * a[1][1] - c[1] * a[1][0]) / det1;
94  float x21 = (a[0][0] * c[1] - a[0][1] * c[0]) / det1;
95 
96  // Float_t det2 = a[0][1]*a[1][2] - a[0][2]*a[1][1];
97  // Float_t det3 = a[0][0]*a[1][2] - a[0][2]*a[1][0];
98  // Float_t x12 = (c[1]*a[1][2] - c[2]*a[1][1])/det2;
99  // Float_t x22 = (a[0][1]*c[2] - a[0][2]*c[1])/det2;
100 
101  float radius =
102  sqrt((b3 + b32 * (x11 * x11 + x21 * x21) - a31 * x11 - a32 * x21) / a33);
103  float centerX = x11;
104  float centerY = x21;
105 
106  ring->SetRadius(radius);
107  ring->SetCenterX(centerX);
108  ring->SetCenterY(centerY);
109 
110  //if (TMath::IsNaN(radius) == 1) ring->SetRadius(0.);
111  //if (TMath::IsNaN(centerX) == 1) ring->SetCenterX(0.);
112  //if (TMath::IsNaN(centerY) == 1) ring->SetCenterY(0.);
113 
114  CalcChi2(ring);
115 }
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmRichHitLight::fX
float fX
Definition: CbmRichRingLight.h:34
CbmRichRingFitterCircle::CbmRichRingFitterCircle
CbmRichRingFitterCircle()
Default constructor.
Definition: CbmRichRingFitterCircle.cxx:11
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
CbmRichRingLight.h
CbmRichRingFitterCircle.h
Implementation of a ring fitting algorithm with equation of a circle. Algorithm from F77 subroutine o...
CbmRichHitLight::fY
float fY
Definition: CbmRichRingLight.h:35
CbmRichRingLight::GetHit
CbmRichHitLight GetHit(int ind)
Return hit by the index.
Definition: CbmRichRingLight.h:114
CbmRichRingFitterCircle::DoFit
void DoFit(CbmRichRingLight *ring)
Inherited from CbmRichRingFitterBase.
Definition: CbmRichRingFitterCircle.cxx:15
CbmRichRingFitterCircle::~CbmRichRingFitterCircle
virtual ~CbmRichRingFitterCircle()
Destructor.
Definition: CbmRichRingFitterCircle.cxx:13
CbmRichRingLight
Definition: CbmRichRingLight.h:39