CbmRoot
CbmRichRingFitterEllipseTau.cxx
Go to the documentation of this file.
1 
8 // GMij are indices for a 5x5 matrix.
9 #define GM00 0
10 #define GM01 1
11 #define GM02 2
12 #define GM03 3
13 #define GM04 4
14 
15 #define GM10 5
16 #define GM11 6
17 #define GM12 7
18 #define GM13 8
19 #define GM14 9
20 
21 #define GM20 10
22 #define GM21 11
23 #define GM22 12
24 #define GM23 13
25 #define GM24 14
26 
27 #define GM30 15
28 #define GM31 16
29 #define GM32 17
30 #define GM33 18
31 #define GM34 19
32 
33 #define GM40 20
34 #define GM41 21
35 #define GM42 22
36 #define GM43 23
37 #define GM44 24
38 
39 
40 // Aij are indices for a 6x6 matrix.
41 
42 #define GA00 0
43 #define GA01 1
44 #define GA02 2
45 #define GA03 3
46 #define GA04 4
47 #define GA05 5
48 
49 #define GA10 6
50 #define GA11 7
51 #define GA12 8
52 #define GA13 9
53 #define GA14 10
54 #define GA15 11
55 
56 #define GA20 12
57 #define GA21 13
58 #define GA22 14
59 #define GA23 15
60 #define GA24 16
61 #define GA25 17
62 
63 #define GA30 18
64 #define GA31 19
65 #define GA32 20
66 #define GA33 21
67 #define GA34 22
68 #define GA35 23
69 
70 #define GA40 24
71 #define GA41 25
72 #define GA42 26
73 #define GA43 27
74 #define GA44 28
75 #define GA45 29
76 
77 #define GA50 30
78 #define GA51 31
79 #define GA52 32
80 #define GA53 33
81 #define GA54 34
82 #define GA55 35
83 
85 
86 
87 using std::cout;
88 using std::endl;
89 
91 
93 
95  int nofHits = ring->GetNofHits();
96 
97  if (nofHits <= 5) {
98  ring->SetXYABP(-1., -1., -1., -1., -1.);
99  ring->SetRadius(-1.);
100  return;
101  }
102 
103  if (nofHits >= MAX_NOF_HITS_IN_RING) {
104  cout
105  << "-E- CbmRichRingFitterEllipseTau::DoFit(), too many hits in the ring:"
106  << nofHits << endl;
107  ring->SetXYABP(-1., -1., -1., -1., -1.);
108  ring->SetRadius(-1.);
109  return;
110  }
111 
112  //for (int i = 0; i < nofHits; i++) {
113  // CbmRichHit* hit = (CbmRichHit*) fHitsArray->At(ring->GetHit(i));
114  // fX.push_back(hit->GetX());
115  // fY.push_back(hit->GetY());
116  //}
117 
118  InitMatrices(ring);
119  Taubin();
120  TransformEllipse(ring);
121  CalcChi2(ring);
122 }
123 
125  // TMatrixD PQ(5,5); // fPQ = P^(-1) * Q
126 
127  Inv5x5();
128  Double_t PQa[25];
129  AMultB(fP, 25, 5, fQ, 25, 5, PQa); //PQ = fP*fQ;
130  TMatrixD PQ(5, 5, PQa);
131 
132  // Double_t PQa2[5][5];
133  // Double_t d[5];
134  // Double_t v[5][5];
135  // for(Int_t i = 0; i<5; i++){
136  // for (Int_t j = 0; j<5;j++){
137  // PQa2[i][j] = PQa[5*i+j];
138  // cout << PQa2[i][j] << " ";
139  // }
140  // cout << endl;
141  // }
142  // cout << endl << endl;
143  //
144  // Jacobi(PQa2, d, v);
145  // Eigsrt(d, v);
146 
147  TMatrixDEigen eig(PQ);
148  TMatrixD evc = eig.GetEigenVectors();
149 
150  Double_t AlgParF = 0.;
151  AlgParF -= evc(0, 0) * fM[GA05];
152  fAlgPar[0] = evc(0, 0);
153 
154  AlgParF -= evc(1, 0) * fM[GA15];
155  fAlgPar[1] = evc(1, 0);
156 
157  AlgParF -= evc(2, 0) * fM[GA25];
158  fAlgPar[2] = evc(2, 0);
159 
160  AlgParF -= evc(3, 0) * fM[GA35];
161  fAlgPar[3] = evc(3, 0);
162 
163  AlgParF -= evc(4, 0) * fM[GA45];
164  fAlgPar[4] = evc(4, 0);
165 
166  fAlgPar[5] = AlgParF;
167 }
168 
170  const unsigned int numHits = ring->GetNofHits();
171  const unsigned int numHits2 = 2 * numHits;
172  const unsigned int numHits3 = 3 * numHits;
173  const unsigned int numHits4 = 4 * numHits;
174  const unsigned int numHits5 = 5 * numHits;
175  const unsigned int numHits6 = 6 * numHits;
176  const double oneOverNumHits = 1. / numHits;
177  unsigned int i6;
178  for (unsigned int i = 0; i < numHits; i++) {
179  double x = ring->GetHit(i).fX;
180  double y = ring->GetHit(i).fY;
181  i6 = i * 6;
182  fZ[i6] = fZT[i] = x * x;
183  fZ[i6 + 1] = fZT[i + numHits] = x * y;
184  fZ[i6 + 2] = fZT[i + numHits2] = y * y;
185  fZ[i6 + 3] = fZT[i + numHits3] = x;
186  fZ[i6 + 4] = fZT[i + numHits4] = y;
187  fZ[i6 + 5] = fZT[i + numHits5] = 1.;
188  }
189  AMultB(fZT, numHits6, numHits, fZ, numHits6, 6, fM);
190  for (unsigned int i = 0; i < numHits6; i++)
191  fM[i] = oneOverNumHits * fM[i];
192 
193  for (int i = 0; i < 25; i++)
194  fP[i] = 0.;
195 
196  fP[GM00] = fM[GA00] - fM[GA05] * fM[GA05];
197  fP[GM01] = fP[GM10] = fM[GA01] - fM[GA05] * fM[GA15];
198  fP[GM02] = fP[GM20] = fM[GA02] - fM[GA05] * fM[GA25];
199  fP[GM03] = fP[GM30] = fM[GA03] - fM[GA05] * fM[GA35];
200  fP[GM04] = fP[GM40] = fM[GA04] - fM[GA05] * fM[GA45];
201 
202  fP[GM11] = fM[GA11] - fM[GA15] * fM[GA15];
203  fP[GM12] = fP[GM21] = fM[GA12] - fM[GA15] * fM[GA25];
204  fP[GM13] = fP[GM31] = fM[GA13] - fM[GA15] * fM[GA35];
205  fP[GM14] = fP[GM41] = fM[GA14] - fM[GA15] * fM[GA45];
206 
207  fP[GM22] = fM[GA22] - fM[GA25] * fM[GA25];
208  fP[GM23] = fP[GM32] = fM[GA23] - fM[GA25] * fM[GA35];
209  fP[GM24] = fP[GM42] = fM[GA24] - fM[GA25] * fM[GA45];
210 
211  fP[GM33] = fM[GA33] - fM[GA35] * fM[GA35];
212  fP[GM34] = fP[GM43] = fM[GA34] - fM[GA35] * fM[GA45];
213 
214  fP[GM44] = fM[GA44] - fM[GA45] * fM[GA45];
215 
216 
217  for (int i = 0; i < 25; i++)
218  fQ[i] = 0.;
219  fQ[GM00] = 4. * fM[GA50];
220  fQ[GM01] = fQ[GM10] = 2. * fM[GA51];
221  fQ[GM03] = fQ[GM30] = 2. * fM[GA53];
222  fQ[GM11] = fM[GA50] + fM[GA52];
223  fQ[GM12] = fQ[GM21] = 2. * fM[GA51];
224  fQ[GM13] = fQ[GM31] = fM[GA54];
225  fQ[GM14] = fQ[GM41] = fM[GA53];
226  fQ[GM22] = 4. * fM[GA52];
227  fQ[GM24] = fQ[GM42] = 2. * fM[GA54];
228  fQ[GM33] = fQ[GM44] = 1.;
229 }
230 
232  double Pxx = fAlgPar[0];
233  double Pxy = fAlgPar[1];
234  double Pyy = fAlgPar[2];
235  double Px = fAlgPar[3];
236  double Py = fAlgPar[4];
237  double P = fAlgPar[5];
238  CalcChi2(Pxx, Pxy, Pyy, Px, Py, P, ring);
239  ring->SetABCDEF(Pxx, Pxy, Pyy, Px, Py, P);
240 
241  double alpha;
242  double QQx, QQy, Qxx, Qyy, Qx, Qy, Q;
243  double cosa, sina, cca, ssa, sin2a;
244  double xc, yc;
245  if (fabs(Pxx - Pyy) > 0.1e-10) {
246  alpha = atan(Pxy / (Pxx - Pyy));
247  alpha = alpha / 2.0;
248  } else
249  alpha = 1.57079633;
250 
251  cosa = cos(alpha);
252  sina = sin(alpha);
253  cca = cosa * cosa;
254  ssa = sina * sina;
255  sin2a = sin(2. * alpha);
256  Pxy = Pxy * sin2a / 2.;
257  Qx = Px * cosa + Py * sina;
258  Qy = -Px * sina + Py * cosa;
259  QQx = Qx * Qx;
260  QQy = Qy * Qy;
261  Qxx = 1. / (Pxx * cca + Pxy + Pyy * ssa);
262  Qyy = 1. / (Pxx * ssa - Pxy + Pyy * cca);
263  Q = -P + Qxx * QQx / 4. + Qyy * QQy / 4.;
264  xc = Qx * Qxx;
265  yc = Qy * Qyy;
266 
267  double axisA = TMath::Sqrt(Q * Qxx);
268  double axisB = TMath::Sqrt(Q * Qyy);
269  double centerX = -xc * cosa / 2. + yc * sina / 2.;
270  double centerY = -xc * sina / 2. - yc * cosa / 2.;
271 
272  ring->SetXYABP(centerX, centerY, axisA, axisB, alpha);
273  ring->SetRadius((axisA + axisB) / 2.);
274 
275  if (ring->GetAaxis() < ring->GetBaxis()) {
276  double tmp = ring->GetAaxis();
277  ring->SetAaxis(ring->GetBaxis());
278  ring->SetBaxis(tmp);
279 
280  tmp = ring->GetPhi();
281  if (ring->GetPhi() <= 0)
282  ring->SetPhi(ring->GetPhi() + 1.57079633);
283  else
284  ring->SetPhi(ring->GetPhi() - 1.57079633);
285  }
286 }
287 
289  // Find all NECESSARY 2x2 dets: (30 of them)
290  const double det2_23_01 = fP[GM20] * fP[GM31] - fP[GM21] * fP[GM30];
291  const double det2_23_02 = fP[GM20] * fP[GM32] - fP[GM22] * fP[GM30];
292  const double det2_23_03 = fP[GM20] * fP[GM33] - fP[GM23] * fP[GM30];
293  const double det2_23_04 = fP[GM20] * fP[GM34] - fP[GM24] * fP[GM30];
294  const double det2_23_12 = fP[GM21] * fP[GM32] - fP[GM22] * fP[GM31];
295  const double det2_23_13 = fP[GM21] * fP[GM33] - fP[GM23] * fP[GM31];
296  const double det2_23_14 = fP[GM21] * fP[GM34] - fP[GM24] * fP[GM31];
297  const double det2_23_23 = fP[GM22] * fP[GM33] - fP[GM23] * fP[GM32];
298  const double det2_23_24 = fP[GM22] * fP[GM34] - fP[GM24] * fP[GM32];
299  const double det2_23_34 = fP[GM23] * fP[GM34] - fP[GM24] * fP[GM33];
300  const double det2_24_01 = fP[GM20] * fP[GM41] - fP[GM21] * fP[GM40];
301  const double det2_24_02 = fP[GM20] * fP[GM42] - fP[GM22] * fP[GM40];
302  const double det2_24_03 = fP[GM20] * fP[GM43] - fP[GM23] * fP[GM40];
303  const double det2_24_04 = fP[GM20] * fP[GM44] - fP[GM24] * fP[GM40];
304  const double det2_24_12 = fP[GM21] * fP[GM42] - fP[GM22] * fP[GM41];
305  const double det2_24_13 = fP[GM21] * fP[GM43] - fP[GM23] * fP[GM41];
306  const double det2_24_14 = fP[GM21] * fP[GM44] - fP[GM24] * fP[GM41];
307  const double det2_24_23 = fP[GM22] * fP[GM43] - fP[GM23] * fP[GM42];
308  const double det2_24_24 = fP[GM22] * fP[GM44] - fP[GM24] * fP[GM42];
309  const double det2_24_34 = fP[GM23] * fP[GM44] - fP[GM24] * fP[GM43];
310  const double det2_34_01 = fP[GM30] * fP[GM41] - fP[GM31] * fP[GM40];
311  const double det2_34_02 = fP[GM30] * fP[GM42] - fP[GM32] * fP[GM40];
312  const double det2_34_03 = fP[GM30] * fP[GM43] - fP[GM33] * fP[GM40];
313  const double det2_34_04 = fP[GM30] * fP[GM44] - fP[GM34] * fP[GM40];
314  const double det2_34_12 = fP[GM31] * fP[GM42] - fP[GM32] * fP[GM41];
315  const double det2_34_13 = fP[GM31] * fP[GM43] - fP[GM33] * fP[GM41];
316  const double det2_34_14 = fP[GM31] * fP[GM44] - fP[GM34] * fP[GM41];
317  const double det2_34_23 = fP[GM32] * fP[GM43] - fP[GM33] * fP[GM42];
318  const double det2_34_24 = fP[GM32] * fP[GM44] - fP[GM34] * fP[GM42];
319  const double det2_34_34 = fP[GM33] * fP[GM44] - fP[GM34] * fP[GM43];
320 
321  // Find all NECESSARY 3x3 dets: (40 of them)
322  const double det3_123_012 =
323  fP[GM10] * det2_23_12 - fP[GM11] * det2_23_02 + fP[GM12] * det2_23_01;
324  const double det3_123_013 =
325  fP[GM10] * det2_23_13 - fP[GM11] * det2_23_03 + fP[GM13] * det2_23_01;
326  const double det3_123_014 =
327  fP[GM10] * det2_23_14 - fP[GM11] * det2_23_04 + fP[GM14] * det2_23_01;
328  const double det3_123_023 =
329  fP[GM10] * det2_23_23 - fP[GM12] * det2_23_03 + fP[GM13] * det2_23_02;
330  const double det3_123_024 =
331  fP[GM10] * det2_23_24 - fP[GM12] * det2_23_04 + fP[GM14] * det2_23_02;
332  const double det3_123_034 =
333  fP[GM10] * det2_23_34 - fP[GM13] * det2_23_04 + fP[GM14] * det2_23_03;
334  const double det3_123_123 =
335  fP[GM11] * det2_23_23 - fP[GM12] * det2_23_13 + fP[GM13] * det2_23_12;
336  const double det3_123_124 =
337  fP[GM11] * det2_23_24 - fP[GM12] * det2_23_14 + fP[GM14] * det2_23_12;
338  const double det3_123_134 =
339  fP[GM11] * det2_23_34 - fP[GM13] * det2_23_14 + fP[GM14] * det2_23_13;
340  const double det3_123_234 =
341  fP[GM12] * det2_23_34 - fP[GM13] * det2_23_24 + fP[GM14] * det2_23_23;
342  const double det3_124_012 =
343  fP[GM10] * det2_24_12 - fP[GM11] * det2_24_02 + fP[GM12] * det2_24_01;
344  const double det3_124_013 =
345  fP[GM10] * det2_24_13 - fP[GM11] * det2_24_03 + fP[GM13] * det2_24_01;
346  const double det3_124_014 =
347  fP[GM10] * det2_24_14 - fP[GM11] * det2_24_04 + fP[GM14] * det2_24_01;
348  const double det3_124_023 =
349  fP[GM10] * det2_24_23 - fP[GM12] * det2_24_03 + fP[GM13] * det2_24_02;
350  const double det3_124_024 =
351  fP[GM10] * det2_24_24 - fP[GM12] * det2_24_04 + fP[GM14] * det2_24_02;
352  const double det3_124_034 =
353  fP[GM10] * det2_24_34 - fP[GM13] * det2_24_04 + fP[GM14] * det2_24_03;
354  const double det3_124_123 =
355  fP[GM11] * det2_24_23 - fP[GM12] * det2_24_13 + fP[GM13] * det2_24_12;
356  const double det3_124_124 =
357  fP[GM11] * det2_24_24 - fP[GM12] * det2_24_14 + fP[GM14] * det2_24_12;
358  const double det3_124_134 =
359  fP[GM11] * det2_24_34 - fP[GM13] * det2_24_14 + fP[GM14] * det2_24_13;
360  const double det3_124_234 =
361  fP[GM12] * det2_24_34 - fP[GM13] * det2_24_24 + fP[GM14] * det2_24_23;
362  const double det3_134_012 =
363  fP[GM10] * det2_34_12 - fP[GM11] * det2_34_02 + fP[GM12] * det2_34_01;
364  const double det3_134_013 =
365  fP[GM10] * det2_34_13 - fP[GM11] * det2_34_03 + fP[GM13] * det2_34_01;
366  const double det3_134_014 =
367  fP[GM10] * det2_34_14 - fP[GM11] * det2_34_04 + fP[GM14] * det2_34_01;
368  const double det3_134_023 =
369  fP[GM10] * det2_34_23 - fP[GM12] * det2_34_03 + fP[GM13] * det2_34_02;
370  const double det3_134_024 =
371  fP[GM10] * det2_34_24 - fP[GM12] * det2_34_04 + fP[GM14] * det2_34_02;
372  const double det3_134_034 =
373  fP[GM10] * det2_34_34 - fP[GM13] * det2_34_04 + fP[GM14] * det2_34_03;
374  const double det3_134_123 =
375  fP[GM11] * det2_34_23 - fP[GM12] * det2_34_13 + fP[GM13] * det2_34_12;
376  const double det3_134_124 =
377  fP[GM11] * det2_34_24 - fP[GM12] * det2_34_14 + fP[GM14] * det2_34_12;
378  const double det3_134_134 =
379  fP[GM11] * det2_34_34 - fP[GM13] * det2_34_14 + fP[GM14] * det2_34_13;
380  const double det3_134_234 =
381  fP[GM12] * det2_34_34 - fP[GM13] * det2_34_24 + fP[GM14] * det2_34_23;
382  const double det3_234_012 =
383  fP[GM20] * det2_34_12 - fP[GM21] * det2_34_02 + fP[GM22] * det2_34_01;
384  const double det3_234_013 =
385  fP[GM20] * det2_34_13 - fP[GM21] * det2_34_03 + fP[GM23] * det2_34_01;
386  const double det3_234_014 =
387  fP[GM20] * det2_34_14 - fP[GM21] * det2_34_04 + fP[GM24] * det2_34_01;
388  const double det3_234_023 =
389  fP[GM20] * det2_34_23 - fP[GM22] * det2_34_03 + fP[GM23] * det2_34_02;
390  const double det3_234_024 =
391  fP[GM20] * det2_34_24 - fP[GM22] * det2_34_04 + fP[GM24] * det2_34_02;
392  const double det3_234_034 =
393  fP[GM20] * det2_34_34 - fP[GM23] * det2_34_04 + fP[GM24] * det2_34_03;
394  const double det3_234_123 =
395  fP[GM21] * det2_34_23 - fP[GM22] * det2_34_13 + fP[GM23] * det2_34_12;
396  const double det3_234_124 =
397  fP[GM21] * det2_34_24 - fP[GM22] * det2_34_14 + fP[GM24] * det2_34_12;
398  const double det3_234_134 =
399  fP[GM21] * det2_34_34 - fP[GM23] * det2_34_14 + fP[GM24] * det2_34_13;
400  const double det3_234_234 =
401  fP[GM22] * det2_34_34 - fP[GM23] * det2_34_24 + fP[GM24] * det2_34_23;
402 
403  // Find all NECESSARY 4x4 dets: (25 of them)
404  const double det4_0123_0123 =
405  fP[GM00] * det3_123_123 - fP[GM01] * det3_123_023 + fP[GM02] * det3_123_013
406  - fP[GM03] * det3_123_012;
407  const double det4_0123_0124 =
408  fP[GM00] * det3_123_124 - fP[GM01] * det3_123_024 + fP[GM02] * det3_123_014
409  - fP[GM04] * det3_123_012;
410  const double det4_0123_0134 =
411  fP[GM00] * det3_123_134 - fP[GM01] * det3_123_034 + fP[GM03] * det3_123_014
412  - fP[GM04] * det3_123_013;
413  const double det4_0123_0234 =
414  fP[GM00] * det3_123_234 - fP[GM02] * det3_123_034 + fP[GM03] * det3_123_024
415  - fP[GM04] * det3_123_023;
416  const double det4_0123_1234 =
417  fP[GM01] * det3_123_234 - fP[GM02] * det3_123_134 + fP[GM03] * det3_123_124
418  - fP[GM04] * det3_123_123;
419  const double det4_0124_0123 =
420  fP[GM00] * det3_124_123 - fP[GM01] * det3_124_023 + fP[GM02] * det3_124_013
421  - fP[GM03] * det3_124_012;
422  const double det4_0124_0124 =
423  fP[GM00] * det3_124_124 - fP[GM01] * det3_124_024 + fP[GM02] * det3_124_014
424  - fP[GM04] * det3_124_012;
425  const double det4_0124_0134 =
426  fP[GM00] * det3_124_134 - fP[GM01] * det3_124_034 + fP[GM03] * det3_124_014
427  - fP[GM04] * det3_124_013;
428  const double det4_0124_0234 =
429  fP[GM00] * det3_124_234 - fP[GM02] * det3_124_034 + fP[GM03] * det3_124_024
430  - fP[GM04] * det3_124_023;
431  const double det4_0124_1234 =
432  fP[GM01] * det3_124_234 - fP[GM02] * det3_124_134 + fP[GM03] * det3_124_124
433  - fP[GM04] * det3_124_123;
434  const double det4_0134_0123 =
435  fP[GM00] * det3_134_123 - fP[GM01] * det3_134_023 + fP[GM02] * det3_134_013
436  - fP[GM03] * det3_134_012;
437  const double det4_0134_0124 =
438  fP[GM00] * det3_134_124 - fP[GM01] * det3_134_024 + fP[GM02] * det3_134_014
439  - fP[GM04] * det3_134_012;
440  const double det4_0134_0134 =
441  fP[GM00] * det3_134_134 - fP[GM01] * det3_134_034 + fP[GM03] * det3_134_014
442  - fP[GM04] * det3_134_013;
443  const double det4_0134_0234 =
444  fP[GM00] * det3_134_234 - fP[GM02] * det3_134_034 + fP[GM03] * det3_134_024
445  - fP[GM04] * det3_134_023;
446  const double det4_0134_1234 =
447  fP[GM01] * det3_134_234 - fP[GM02] * det3_134_134 + fP[GM03] * det3_134_124
448  - fP[GM04] * det3_134_123;
449  const double det4_0234_0123 =
450  fP[GM00] * det3_234_123 - fP[GM01] * det3_234_023 + fP[GM02] * det3_234_013
451  - fP[GM03] * det3_234_012;
452  const double det4_0234_0124 =
453  fP[GM00] * det3_234_124 - fP[GM01] * det3_234_024 + fP[GM02] * det3_234_014
454  - fP[GM04] * det3_234_012;
455  const double det4_0234_0134 =
456  fP[GM00] * det3_234_134 - fP[GM01] * det3_234_034 + fP[GM03] * det3_234_014
457  - fP[GM04] * det3_234_013;
458  const double det4_0234_0234 =
459  fP[GM00] * det3_234_234 - fP[GM02] * det3_234_034 + fP[GM03] * det3_234_024
460  - fP[GM04] * det3_234_023;
461  const double det4_0234_1234 =
462  fP[GM01] * det3_234_234 - fP[GM02] * det3_234_134 + fP[GM03] * det3_234_124
463  - fP[GM04] * det3_234_123;
464  const double det4_1234_0123 =
465  fP[GM10] * det3_234_123 - fP[GM11] * det3_234_023 + fP[GM12] * det3_234_013
466  - fP[GM13] * det3_234_012;
467  const double det4_1234_0124 =
468  fP[GM10] * det3_234_124 - fP[GM11] * det3_234_024 + fP[GM12] * det3_234_014
469  - fP[GM14] * det3_234_012;
470  const double det4_1234_0134 =
471  fP[GM10] * det3_234_134 - fP[GM11] * det3_234_034 + fP[GM13] * det3_234_014
472  - fP[GM14] * det3_234_013;
473  const double det4_1234_0234 =
474  fP[GM10] * det3_234_234 - fP[GM12] * det3_234_034 + fP[GM13] * det3_234_024
475  - fP[GM14] * det3_234_023;
476  const double det4_1234_1234 =
477  fP[GM11] * det3_234_234 - fP[GM12] * det3_234_134 + fP[GM13] * det3_234_124
478  - fP[GM14] * det3_234_123;
479 
480  // Find the 5x5 det:
481  const double det = fP[GM00] * det4_1234_1234 - fP[GM01] * det4_1234_0234
482  + fP[GM02] * det4_1234_0134 - fP[GM03] * det4_1234_0124
483  + fP[GM04] * det4_1234_0123;
484  // if (determ)
485  // *determ = det;
486  //
487  // if (det == 0) {
488  // Error("Inv5x5", "matrix is singular");
489  // return kFALSE;
490  // }
491  const double oneOverDet = 1.0 / det;
492  const double mn1OverDet = -oneOverDet;
493 
494  fP[GM00] = det4_1234_1234 * oneOverDet;
495  fP[GM01] = det4_0234_1234 * mn1OverDet;
496  fP[GM02] = det4_0134_1234 * oneOverDet;
497  fP[GM03] = det4_0124_1234 * mn1OverDet;
498  fP[GM04] = det4_0123_1234 * oneOverDet;
499 
500  fP[GM10] = det4_1234_0234 * mn1OverDet;
501  fP[GM11] = det4_0234_0234 * oneOverDet;
502  fP[GM12] = det4_0134_0234 * mn1OverDet;
503  fP[GM13] = det4_0124_0234 * oneOverDet;
504  fP[GM14] = det4_0123_0234 * mn1OverDet;
505 
506  fP[GM20] = det4_1234_0134 * oneOverDet;
507  fP[GM21] = det4_0234_0134 * mn1OverDet;
508  fP[GM22] = det4_0134_0134 * oneOverDet;
509  fP[GM23] = det4_0124_0134 * mn1OverDet;
510  fP[GM24] = det4_0123_0134 * oneOverDet;
511 
512  fP[GM30] = det4_1234_0124 * mn1OverDet;
513  fP[GM31] = det4_0234_0124 * oneOverDet;
514  fP[GM32] = det4_0134_0124 * mn1OverDet;
515  fP[GM33] = det4_0124_0124 * oneOverDet;
516  fP[GM34] = det4_0123_0124 * mn1OverDet;
517 
518  fP[GM40] = det4_1234_0123 * oneOverDet;
519  fP[GM41] = det4_0234_0123 * mn1OverDet;
520  fP[GM42] = det4_0134_0123 * oneOverDet;
521  fP[GM43] = det4_0124_0123 * mn1OverDet;
522  fP[GM44] = det4_0123_0123 * oneOverDet;
523 }
524 
525 void CbmRichRingFitterEllipseTau::AMultB(const double* const ap,
526  int na,
527  int ncolsa,
528  const double* const bp,
529  int nb,
530  int ncolsb,
531  double* cp) {
532  // Elementary routine to calculate matrix multiplication A*B
533 
534  const double* arp0 = ap; // Pointer to A[i,0];
535  while (arp0 < ap + na) {
536  for (
537  const double* bcp = bp;
538  bcp
539  < bp + ncolsb;) { // Pointer to the j-th column of B, Start bcp = B[0,0]
540  const double* arp =
541  arp0; // Pointer to the i-th row of A, reset to A[i,0]
542  double cij = 0;
543  while (bcp < bp + nb) { // Scan the i-th row of A and
544  cij += *arp++ * *bcp; // the j-th col of B
545  bcp += ncolsb;
546  }
547  *cp++ = cij;
548  bcp -= nb - 1; // Set bcp to the (j+1)-th col
549  }
550  arp0 += ncolsa; // Set ap to the (i+1)-th row
551  }
552 }
553 
554 #define ROTATE(a, i, j, k, l) \
555  g = a[i][j]; \
556  h = a[k][l]; \
557  a[i][j] = g - s * (h + g * tau); \
558  a[k][l] = h + s * (g - h * tau)
559 #define MAXSWEEP 50
561  double d[5],
562  double v[5][5]) {
563  double tresh, theta, tau, t, sm, s, h, g, c;
564 
565  double b[5], z[5];
566  unsigned int ip, iq, i, j;
567  for (ip = 0; ip < 5; ip++) {
568  for (iq = 0; iq < 5; iq++)
569  v[ip][iq] = 0.;
570  v[ip][ip] = 1.;
571  }
572 
573  for (ip = 0; ip < 5; ip++) {
574  b[ip] = a[ip][ip];
575  z[ip] = 0.;
576  }
577 
578  int nrot = 0;
579 
580  for (i = 1; i <= MAXSWEEP; i++) {
581 
582  for (sm = 0., ip = 0; ip < 5; ip++)
583  for (iq = ip + 1; iq < 5; iq++)
584  sm += fabs(a[ip][iq]);
585  if (sm == 0.) { return; }
586 
587  tresh = (i < 4 ? 0.2 * sm / (5 * 5) : 0.);
588 
589  for (ip = 0; ip < 4; ip++)
590  for (iq = ip + 1; iq < 5; iq++) {
591 
592  g = 100. * fabs(a[ip][iq]);
593  if (i > 4 && (float) fabs(d[ip] + g) == (float) fabs(d[ip])
594  && (float) fabs(d[iq] + g) == (float) fabs(d[iq]))
595  a[ip][iq] = 0.;
596 
597  else if (fabs(a[ip][iq]) > tresh) {
598  h = d[ip] - d[iq];
599 
600  if ((float) (fabs(h) + g) == (float) fabs(h))
601  t = a[ip][iq] / h;
602  else {
603  theta = 0.5 * h / a[ip][iq];
604  t = 1. / (fabs(theta) + sqrt(1. + theta * theta));
605  if (theta < 0.) t = -t;
606  }
607  c = 1. / sqrt(1 + t * t);
608  s = t * c;
609  tau = s / (1. + c);
610  h = t * a[ip][iq];
611  z[ip] -= h;
612  z[iq] += h;
613  d[ip] -= h;
614  d[iq] += h;
615  a[ip][iq] = 0.;
616  for (j = 0; j < ip; j++) {
617  ROTATE(a, j, ip, j, iq);
618  }
619  for (j = ip + 1; j < iq; j++) {
620  ROTATE(a, ip, j, iq, j);
621  }
622  for (j = iq + 1; j < 5; j++) {
623  ROTATE(a, ip, j, j, iq);
624  }
625  for (j = 0; j < 5; j++) {
626  ROTATE(v, j, ip, j, iq);
627  }
628  ++nrot;
629  }
630  }
631  for (ip = 0; ip < 5; ip++) {
632  b[ip] += z[ip];
633  d[ip] = b[ip];
634  z[ip] = 0.;
635  }
636  } //i rot
637 }
638 
639 void CbmRichRingFitterEllipseTau::Eigsrt(double d[5], double v[5][5]) {
640  double p;
641  int i, k, j;
642  for (i = 0; i < 5; i++) {
643  p = d[k = i];
644  for (j = i + 1; j < 5; j++)
645  if (d[j] >= p) p = d[k = j];
646  if (k != i) {
647  d[k] = d[i];
648  d[i] = p;
649  for (j = 0; j < 5; j++) {
650  p = v[j][i];
651  v[j][i] = v[j][k];
652  v[j][k] = p;
653  }
654  }
655  }
656 }
CbmRichRingFitterEllipseBase::CalcChi2
virtual void CalcChi2(CbmRichRingLight *ring)
Calculate chi2 of the ellipse fit.
Definition: CbmRichRingFitterEllipseBase.h:42
CbmRichRingFitterEllipseTau::TransformEllipse
void TransformEllipse(CbmRichRingLight *ring)
Transform fitted curve to ellipse parameters.
Definition: CbmRichRingFitterEllipseTau.cxx:231
GA13
#define GA13
Definition: CbmRichRingFitterEllipseTau.cxx:52
CbmRichRingFitterEllipseTau::fM
double fM[36]
Definition: CbmRichRingFitterEllipseTau.h:53
sin
friend F32vec4 sin(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:136
GA52
#define GA52
Definition: CbmRichRingFitterEllipseTau.cxx:79
GM24
#define GM24
Definition: CbmRichRingFitterEllipseTau.cxx:25
GA23
#define GA23
Definition: CbmRichRingFitterEllipseTau.cxx:59
GA22
#define GA22
Definition: CbmRichRingFitterEllipseTau.cxx:58
GM00
#define GM00
Definition: CbmRichRingFitterEllipseTau.cxx:9
GM42
#define GM42
Definition: CbmRichRingFitterEllipseTau.cxx:35
CbmRichRingFitterEllipseTau::InitMatrices
void InitMatrices(CbmRichRingLight *ring)
Initialize all matrices.
Definition: CbmRichRingFitterEllipseTau.cxx:169
GA35
#define GA35
Definition: CbmRichRingFitterEllipseTau.cxx:68
CbmRichRingFitterEllipseTau.h
Here the ring is fitted with Taubin algorithm from A. Ayriyan, G. Ososkov, N. Chernov.
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
GM14
#define GM14
Definition: CbmRichRingFitterEllipseTau.cxx:19
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
GM02
#define GM02
Definition: CbmRichRingFitterEllipseTau.cxx:11
GM04
#define GM04
Definition: CbmRichRingFitterEllipseTau.cxx:13
CbmRichRingLight::SetABCDEF
void SetABCDEF(float a, float b, float c, float d, float e, float f)
Set all 6 parameters of curve equation Axx+Bxy+Cyy+Dx+Ey+F.
Definition: CbmRichRingLight.h:225
GM44
#define GM44
Definition: CbmRichRingFitterEllipseTau.cxx:37
CbmRichRingLight::SetAaxis
void SetAaxis(double a)
Definition: CbmRichRingLight.h:166
CbmRichRingFitterEllipseTau::Inv5x5
void Inv5x5()
Invert 5x5 matrix.
Definition: CbmRichRingFitterEllipseTau.cxx:288
GM01
#define GM01
Definition: CbmRichRingFitterEllipseTau.cxx:10
CbmRichRingFitterEllipseTau::Eigsrt
void Eigsrt(double d[5], double v[5][5])
Find eigenvalues.
Definition: CbmRichRingFitterEllipseTau.cxx:639
MAXSWEEP
#define MAXSWEEP
Definition: CbmRichRingFitterEllipseTau.cxx:559
GA04
#define GA04
Definition: CbmRichRingFitterEllipseTau.cxx:46
CbmRichHitLight::fX
float fX
Definition: CbmRichRingLight.h:34
GA50
#define GA50
Definition: CbmRichRingFitterEllipseTau.cxx:77
CbmRichRingLight::SetRadius
void SetRadius(float r)
Definition: CbmRichRingLight.h:124
CbmRichRingFitterEllipseTau::fP
double fP[25]
Definition: CbmRichRingFitterEllipseTau.h:54
GA44
#define GA44
Definition: CbmRichRingFitterEllipseTau.cxx:74
GA15
#define GA15
Definition: CbmRichRingFitterEllipseTau.cxx:54
GM23
#define GM23
Definition: CbmRichRingFitterEllipseTau.cxx:24
CbmRichRingFitterEllipseTau::fAlgPar
double fAlgPar[6]
Definition: CbmRichRingFitterEllipseTau.h:58
GA33
#define GA33
Definition: CbmRichRingFitterEllipseTau.cxx:66
CbmRichRingLight::GetNofHits
int GetNofHits() const
Return number of hits in ring.
Definition: CbmRichRingLight.h:108
CbmRichRingLight::SetXYABP
void SetXYABP(float x, float y, float a, float b, float p)
Set all 5 ellipse parameters.
Definition: CbmRichRingLight.h:146
GA45
#define GA45
Definition: CbmRichRingFitterEllipseTau.cxx:75
CbmRichRingLight::GetAaxis
float GetAaxis() const
Definition: CbmRichRingLight.h:163
GA01
#define GA01
Definition: CbmRichRingFitterEllipseTau.cxx:43
h
Data class with information on a STS local track.
CbmRichRingFitterEllipseTau::fQ
double fQ[25]
Definition: CbmRichRingFitterEllipseTau.h:55
d
double d
Definition: P4_F64vec2.h:24
GM12
#define GM12
Definition: CbmRichRingFitterEllipseTau.cxx:17
GA24
#define GA24
Definition: CbmRichRingFitterEllipseTau.cxx:60
CbmRichRingFitterEllipseTau::Taubin
void Taubin()
Perform Taubin method.
Definition: CbmRichRingFitterEllipseTau.cxx:124
GM41
#define GM41
Definition: CbmRichRingFitterEllipseTau.cxx:34
GM43
#define GM43
Definition: CbmRichRingFitterEllipseTau.cxx:36
ROTATE
#define ROTATE(a, i, j, k, l)
Definition: CbmRichRingFitterEllipseTau.cxx:554
GM13
#define GM13
Definition: CbmRichRingFitterEllipseTau.cxx:18
CbmRichRingFitterEllipseTau::AMultB
void AMultB(const double *const ap, int na, int ncolsa, const double *const bp, int nb, int ncolsb, double *cp)
Matrices multiplication.
Definition: CbmRichRingFitterEllipseTau.cxx:525
GM11
#define GM11
Definition: CbmRichRingFitterEllipseTau.cxx:16
GM21
#define GM21
Definition: CbmRichRingFitterEllipseTau.cxx:22
GA12
#define GA12
Definition: CbmRichRingFitterEllipseTau.cxx:51
CbmRichRingLight::GetBaxis
float GetBaxis() const
Definition: CbmRichRingLight.h:164
GM10
#define GM10
Definition: CbmRichRingFitterEllipseTau.cxx:15
CbmRichRingFitterEllipseTau::fZ
double fZ[MAX_NOF_HITS_IN_RING *6]
Definition: CbmRichRingFitterEllipseTau.h:56
GA53
#define GA53
Definition: CbmRichRingFitterEllipseTau.cxx:80
CbmRichRingFitterEllipseTau::CbmRichRingFitterEllipseTau
CbmRichRingFitterEllipseTau()
Default constructor.
Definition: CbmRichRingFitterEllipseTau.cxx:90
CbmRichRingFitterBase::MAX_NOF_HITS_IN_RING
static const int MAX_NOF_HITS_IN_RING
Definition: CbmRichRingFitterBase.h:70
GM31
#define GM31
Definition: CbmRichRingFitterEllipseTau.cxx:28
GA11
#define GA11
Definition: CbmRichRingFitterEllipseTau.cxx:50
GA00
#define GA00
Definition: CbmRichRingFitterEllipseTau.cxx:42
v
__m128 v
Definition: L1/vectors/P4_F32vec4.h:1
GA14
#define GA14
Definition: CbmRichRingFitterEllipseTau.cxx:53
GA02
#define GA02
Definition: CbmRichRingFitterEllipseTau.cxx:44
CbmRichHitLight::fY
float fY
Definition: CbmRichRingLight.h:35
GM33
#define GM33
Definition: CbmRichRingFitterEllipseTau.cxx:30
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
cos
friend F32vec4 cos(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:137
CbmRichRingLight::SetPhi
void SetPhi(double phi)
Definition: CbmRichRingLight.h:168
GM32
#define GM32
Definition: CbmRichRingFitterEllipseTau.cxx:29
GA25
#define GA25
Definition: CbmRichRingFitterEllipseTau.cxx:61
GM22
#define GM22
Definition: CbmRichRingFitterEllipseTau.cxx:23
GA51
#define GA51
Definition: CbmRichRingFitterEllipseTau.cxx:78
CbmRichRingFitterEllipseTau::fZT
double fZT[MAX_NOF_HITS_IN_RING *6]
Definition: CbmRichRingFitterEllipseTau.h:57
GM20
#define GM20
Definition: CbmRichRingFitterEllipseTau.cxx:21
CbmRichRingLight::SetBaxis
void SetBaxis(double b)
Definition: CbmRichRingLight.h:167
CbmRichRingFitterEllipseTau::~CbmRichRingFitterEllipseTau
virtual ~CbmRichRingFitterEllipseTau()
Destructor.
Definition: CbmRichRingFitterEllipseTau.cxx:92
GM03
#define GM03
Definition: CbmRichRingFitterEllipseTau.cxx:12
GA05
#define GA05
Definition: CbmRichRingFitterEllipseTau.cxx:47
CbmRichRingFitterEllipseTau::DoFit
virtual void DoFit(CbmRichRingLight *ring)
Inherited from CbmRichRingFitterBase.
Definition: CbmRichRingFitterEllipseTau.cxx:94
GA34
#define GA34
Definition: CbmRichRingFitterEllipseTau.cxx:67
GA03
#define GA03
Definition: CbmRichRingFitterEllipseTau.cxx:45
CbmRichRingLight
Definition: CbmRichRingLight.h:39
GA54
#define GA54
Definition: CbmRichRingFitterEllipseTau.cxx:81
GM40
#define GM40
Definition: CbmRichRingFitterEllipseTau.cxx:33
GM30
#define GM30
Definition: CbmRichRingFitterEllipseTau.cxx:27
GM34
#define GM34
Definition: CbmRichRingFitterEllipseTau.cxx:31
CbmRichRingLight::GetPhi
float GetPhi() const
Definition: CbmRichRingLight.h:165
CbmRichRingFitterEllipseTau::Jacobi
void Jacobi(double a[5][5], double d[5], double v[5][5])
Jacobi method.
Definition: CbmRichRingFitterEllipseTau.cxx:560