BmnRoot
Loading...
Searching...
No Matches
BmnMath.cxx
Go to the documentation of this file.
1#include "BmnMath.h"
2
3#include "BmnFieldMap.h"
4#include "CbmGlobalTrack.h"
5#include "CbmHit.h"
6#include "CbmStripHit.h"
7#include "FairRunAna.h"
8#include "FairTrackParam.h"
9#include "Math/BrentRootFinder.h"
10#include "Math/WrappedTF1.h"
11#include "TArc.h"
12#include "TCanvas.h"
13#include "TClonesArray.h"
14#include "TF1.h"
15#include "TFitResult.h"
16#include "TGraph.h"
17#include "TH2F.h"
18#include "TLine.h"
19#include "TMath.h"
20
21#include <cmath>
22#include <iostream>
23
24using namespace TMath;
25
26namespace lit
27{
28
29Float_t ChiSq(const FairTrackParam* par, const BmnHit* hit)
30{
31 Float_t dxx = hit->GetDx() * hit->GetDx();
32 Float_t dxy = 0.0; // hit->GetDxy();
33 Float_t dyy = hit->GetDy() * hit->GetDy();
34 Float_t xmx = hit->GetX() - par->GetX();
35 Float_t ymy = hit->GetY() - par->GetY();
36 Float_t C0 = par->GetCovariance(0, 0);
37 Float_t C1 = par->GetCovariance(0, 1);
38 Float_t C5 = par->GetCovariance(1, 1);
39
40 Float_t norm = dxx * dyy - dxx * C5 - dyy * C0 + C0 * C5 - dxy * dxy + 2 * dxy * C1 - C1 * C1;
41 if (norm == 0.) {
42 norm = 1e-10;
43 }
44 return ((xmx * (dyy - C5) - ymy * (dxy - C1)) * xmx + (-xmx * (dxy - C1) + ymy * (dxx - C0)) * ymy) / norm;
45}
46
47Int_t NDF(const BmnTrack* track)
48{
49 Int_t ndf = 0;
50 for (Int_t i = 0; i < track->GetNHits(); i++)
51 ndf += 2;
52 ndf -= 5;
53 return ndf;
54}
55
56} // namespace lit
57
58// checking with flag isField is needed for backward compatibility
59Bool_t IsParCorrect(const FairTrackParam* par, const Bool_t isField)
60{
61 return IsParCorrect(par);
62}
63
64Bool_t IsParCorrect(const FairTrackParam* par)
65{
66 const Float_t maxSlopeX = 5.;
67 const Float_t maxSlopeY = 1.;
68 const Float_t maxX = 300.0;
69 const Float_t maxY = 200.0;
70 // const Float_t minSlope = 1e-10;
71 const Float_t maxQp = 100.; // p = 10 MeV
72 const Float_t minQp = 0.01; // p = 100 GeV
73
74 if (abs(par->GetTx()) > maxSlopeX || abs(par->GetTy()) > maxSlopeY)
75 return kFALSE;
76 if (abs(par->GetX()) > maxX || abs(par->GetY()) > maxY)
77 return kFALSE;
78 if (IsNaN(par->GetX()) || IsNaN(par->GetY()) || IsNaN(par->GetTx()) || IsNaN(par->GetTy()))
79 return kFALSE;
80
81 BmnFieldMap* field = (BmnFieldMap*)FairRunAna::Instance()->GetField();
82 Bool_t isField = !(field->IsFieldOff());
83
84 if (isField) {
85 if (abs(par->GetQp()) > maxQp)
86 return kFALSE;
87 if (abs(par->GetQp()) < minQp)
88 return kFALSE;
89 if (IsNaN(par->GetQp()))
90 return kFALSE;
91 }
92
93 return kTRUE;
94}
95
96Float_t NumericalRootFinder(TF1 f, Float_t left, Float_t right)
97{
98 // Create the wrapper for function
99 ROOT::Math::WrappedTF1 wf1(f);
100
101 // Create the Integrator
102 ROOT::Math::BrentRootFinder brf;
103
104 // Set parameters of the method
105 brf.SetFunction(wf1, left, right);
106 brf.Solve();
107
108 // cout << brf.Root() << endl;
109
110 return brf.Root();
111}
112
113TVector3 LineFit(BmnTrack* track, const TClonesArray* arr, TString type)
114{
115 // Weighted Least Square Method//
116 Float_t Xi = 0.0, Yi = 0.0; // coordinates of current track point
117 Float_t a = 0.0, b = 0.0; // parameters of line: y = a * x + b
118
119 Float_t Si = 0.0; // sigma
120 Float_t Wi = 0.0; // weight = 1 / sigma^2
121 Float_t SumW = 0.0; // sum of weights
122 Float_t SumWX = 0.0; // sum of (weight * x)
123 Float_t SumWY = 0.0; // sum of (weight * y)
124 Float_t SumWXY = 0.0; // sum of (weight * x * y)
125 Float_t SumWX2 = 0.0; // sum of (weight * x * x)
126
127 const Float_t nHits = track->GetNHits();
128
129 for (Int_t i = 0; i < nHits; ++i) {
130 BmnHit* hit = (BmnHit*)arr->At(track->GetHitIndex(i));
131 if (type.Contains("XY")) {
132 Xi = hit->GetX();
133 Yi = hit->GetY();
134 Si = hit->GetDy();
135 } else if (type.Contains("ZX")) {
136 Xi = hit->GetZ();
137 Yi = hit->GetX();
138 Si = hit->GetDx();
139 // h_Hits->Fill(hit->GetZ(), hit->GetX());
140 } else if (type.Contains("ZY")) {
141 Xi = hit->GetZ();
142 Yi = hit->GetY();
143 Si = hit->GetDy();
144 }
145
146 if (Si == 0.0)
147 return TVector3(0.0, 0.0, 0.0);
148
149 Wi = 1.0 / Si / Si;
150 SumW += Wi;
151 SumWXY += Wi * Xi * Yi;
152 SumWX += Wi * Xi;
153 SumWX2 += Wi * Xi * Xi;
154 SumWY += Wi * Yi;
155 }
156
157 a = (SumW * SumWXY - SumWX * SumWY) / (SumW * SumWX2 - SumWX * SumWX);
158 b = (SumWX2 * SumWY - SumWX * SumWXY) / (SumW * SumWX2 - SumWX * SumWX);
159
160 // z1, x1, z2, x2
161 // TLine* line = new TLine();
162 // line->SetFillStyle(0);
163 // line->SetLineWidth(2);
164 // line->SetLineColor(kBlue);
165 // line->SetNoEdges(1);
166
167 // TCanvas* c_New = new TCanvas("c", "c", 1000, 500);
168 // c_New->cd();
169 // h_Hits->SetMarkerStyle(20);
170 // h_Hits->SetMarkerSize(1.5);
171 // h_Hits->SetMarkerColor(kRed);
172 // h_Hits->Draw("P");
173 // h_HitsAll->SetMarkerStyle(20);
174 // h_HitsAll->SetMarkerSize(0.75);
175 // h_HitsAll->SetMarkerColor(kBlue);
176 // h_HitsAll->Draw("Psame");
177 // line->Draw("same");
178 // c_New->SaveAs("hits.png");
179 // getchar();
180 // delete h_Hits;
181 // delete c_New;
182 // delete line;
183
184 Float_t chi2 = 0.0;
185
186 for (Int_t i = 0; i < nHits; ++i) {
187 BmnHit* hit = (BmnHit*)arr->At(track->GetHitIndex(i));
188 if (type.Contains("XY")) {
189 Xi = hit->GetX();
190 Yi = hit->GetY();
191 Si = hit->GetDy();
192 } else if (type.Contains("ZX")) {
193 Xi = hit->GetZ();
194 Yi = hit->GetX();
195 Si = hit->GetDx();
196 } else if (type.Contains("ZY")) {
197 Xi = hit->GetZ();
198 Yi = hit->GetY();
199 Si = hit->GetDy();
200 }
201
202 chi2 += Sq((Yi - a * Xi - b) / Si);
203 // chi2 += Sq((Yi - a * Xi - b) / (a * Xi + b));
204 }
205
206 return TVector3(a, b, chi2);
207}
208
209TVector3 LineFitBy3Hits(const BmnHit* h0, const BmnHit* h1, const BmnHit* h2)
210{
211 // Weighted Least Square Method//
212
213 // sigma
214 Double_t S0 = h0->GetDy();
215 Double_t S1 = h1->GetDy();
216 Double_t S2 = h2->GetDy();
217 // weight = 1 / sigma^2
218 Double_t W0 = 1.0 / S0 / S0;
219 Double_t W1 = 1.0 / S1 / S1;
220 Double_t W2 = 1.0 / S2 / S2;
221
222 Double_t X0 = h0->GetZ();
223 Double_t X1 = h1->GetZ();
224 Double_t X2 = h2->GetZ();
225
226 Double_t Y0 = h0->GetY();
227 Double_t Y1 = h1->GetY();
228 Double_t Y2 = h2->GetY();
229
230 Double_t SumW = W0 + W1 + W2; // sum of weights
231 Double_t SumWX = W0 * X0 + W1 * X1 + W2 * X2; // sum of (weight * x)
232 Double_t SumWY = W0 * Y0 + W1 * Y1 + W2 * Y2; // sum of (weight * y)
233 Double_t SumWXY = W0 * X0 * Y0 + W1 * X1 * Y1 + W2 * X2 * Y2; // sum of (weight * x * y)
234 Double_t SumWX2 = W0 * X0 * X0 + W1 * X1 * X1 + W2 * X2 * X2; // sum of (weight * x * x)
235
236 // parameters of line: y = a * x + b
237 Double_t a = (SumW * SumWXY - SumWX * SumWY) / (SumW * SumWX2 - SumWX * SumWX);
238 Double_t b = (SumWX2 * SumWY - SumWX * SumWXY) / (SumW * SumWX2 - SumWX * SumWX);
239
240 Double_t chi2 = Sq((Y0 - a * X0 - b) / S0) + Sq((Y1 - a * X1 - b) / S1) + Sq((Y2 - a * X2 - b) / S2);
241
242 return TVector3(a, b, chi2);
243}
244
245void LineFit(Double_t& par1, Double_t& par2, BmnTrack* track, TClonesArray* arr, Int_t type, Int_t idSkip)
246{
247 // Weighted Least Square Method//
248 Float_t Xi = 0.0, Yi = 0.0; // coordinates of current track point
249 Float_t a = 0.0, b = 0.0; // parameters of line: y = a * x + b
250
251 Float_t Si = 0.0; // sigma
252 Float_t Wi = 0.0; // weight = 1 / sigma^2
253 Float_t SumW = 0.0; // sum of weights
254 Float_t SumWX = 0.0; // sum of (weight * x)
255 Float_t SumWY = 0.0; // sum of (weight * y)
256 Float_t SumWXY = 0.0; // sum of (weight * x * y)
257 Float_t SumWX2 = 0.0; // sum of (weight * x * x)
258
259 const Float_t nHits = track->GetNHits();
260 for (Int_t i = 0; i < nHits; ++i) {
261 BmnHit* hit = (BmnHit*)arr->At(track->GetHitIndex(i));
262
263 if (i == idSkip)
264 continue;
265
266 else if (type == 1) {
267 Xi = hit->GetZ();
268 Yi = hit->GetX();
269 Si = hit->GetDx();
270 } else {
271 Xi = hit->GetZ();
272 Yi = hit->GetY();
273 Si = hit->GetDy();
274 }
275
276 Wi = 1.0 / Si / Si;
277 SumW += Wi;
278 SumWXY += Wi * Xi * Yi;
279 SumWX += Wi * Xi;
280 SumWX2 += Wi * Xi * Xi;
281 SumWY += Wi * Yi;
282 }
283
284 a = (SumW * SumWXY - SumWX * SumWY) / (SumW * SumWX2 - SumWX * SumWX);
285 b = (SumWX2 * SumWY - SumWX * SumWXY) / (SumW * SumWX2 - SumWX * SumWX);
286
287 Float_t chi2 = 0.0;
288
289 for (Int_t i = 0; i < nHits; ++i) {
290 BmnHit* hit = (BmnHit*)arr->At(track->GetHitIndex(i));
291
292 if (i == idSkip)
293 continue;
294
295 if (type == 1) {
296 Xi = hit->GetZ();
297 Yi = hit->GetX();
298 Si = hit->GetDx();
299 } else {
300 Xi = hit->GetZ();
301 Yi = hit->GetY();
302 Si = hit->GetDy();
303 }
304
305 chi2 += ((Yi - a * Xi - b) / Si) * ((Yi - a * Xi - b) / Si);
306 }
307
308 par1 = a;
309 par2 = b;
310}
311
312TVector3 CircleFit(BmnTrack* track, const TClonesArray* arr, Double_t& chi2)
313{
314 // Weighted Least Square Method//
315 Double_t Xi = 0.0, Yi = 0.0, Zi = 0.0; // coordinates of current track point
316 Double_t Xc = 0.0, Zc = 0.0, R = 0.0;
317 chi2 = 0.0;
318
319 Double_t Wi = 0.0; // weight = 1 / sigma^2
320
321 Double_t Sx = 0.0; // sum of weights
322 Double_t Sxx = 0.0;
323 Double_t Syy = 0.0;
324 Double_t Sxy = 0.0;
325 Double_t Sy = 0.0;
326 Double_t Sz = 0.0;
327 Double_t Szx = 0.0;
328 Double_t Szy = 0.0;
329
330 // TH2F* h_Hits = new TH2F("h_Hits", "h_Hits", 400, 0.0, 250.0, 400, -10.0, 10.0);
331
332 const Float_t nHits = track->GetNHits();
333 for (Int_t i = 0; i < nHits; ++i) {
334 BmnHit* hit = (BmnHit*)arr->At(track->GetHitIndex(i));
335 // Use Z and X coordinates of hits to fit in ZX plane
336 Yi = hit->GetZ();
337 Xi = hit->GetX();
338 // h_Hits->Fill(hit->GetZ(), hit->GetX());
339 Zi = Xi * Xi + Yi * Yi;
340 Wi = 1.0 / hit->GetDx() / hit->GetDx();
341
342 Sx += Wi * Xi;
343 Sy += Wi * Yi;
344 Sz += Wi * Zi;
345 Sxx += Wi * Xi * Xi;
346 Sxy += Wi * Xi * Yi;
347 Syy += Wi * Yi * Yi;
348 Szx += Wi * Zi * Xi;
349 Szy += Wi * Zi * Yi;
350 }
351
352 Double_t C = ((Sz * Sx - Szx) / (Sxx - Sx * Sx) - (Sz * Sy - Szy) / (Sxy - Sx * Sy))
353 / ((Sxy - Sx * Sy) / (Sxx - Sx * Sx) - (Syy - Sy * Sy) / (Sxy - Sx * Sy));
354 Double_t B = ((Sz * Sx - Szx) - C * (Sxy - Sx * Sy)) / (Sxx - Sx * Sx);
355 Double_t D = -Sz - B * Sx - C * Sy;
356
357 Xc = -0.5 * B;
358 Zc = -0.5 * C;
359 R = Sqrt(0.25 * B * B + 0.25 * C * C - D);
360
361 for (Int_t i = 0; i < nHits; ++i) {
362 BmnHit* hit = (BmnHit*)arr->At(track->GetHitIndex(i));
363 chi2 += Sq((R - Sqrt(Sq(hit->GetX() - Xc) + Sq(hit->GetZ() - Zc))) / hit->GetDx());
364 }
365
366 return TVector3(Zc, Xc, R);
367}
368
369TVector3 CircleFit(vector<BmnHit*> hits, Int_t idSkip)
370{
371 // Weighted Least Square Method//
372 Double_t Xi = 0.0, Yi = 0.0, Zi = 0.0; // coordinates of current track point
373 Double_t Xc = 0.0, Zc = 0.0, R = 0.0;
374
375 Double_t Wi = 0.0; // weight = 1 / sigma^2
376
377 Double_t Sx = 0.0; // sum of weights
378 Double_t Sxx = 0.0;
379 Double_t Syy = 0.0;
380 Double_t Sxy = 0.0;
381 Double_t Sy = 0.0;
382 Double_t Sz = 0.0;
383 Double_t Szx = 0.0;
384 Double_t Szy = 0.0;
385
386 // TH2F* h_Hits = new TH2F("h_Hits", "h_Hits", 400, 0.0, 250.0, 400, -10.0, 10.0);
387
388 const Float_t nHits = hits.size();
389 for (Int_t i = 0; i < nHits; ++i) {
390 if (i == idSkip)
391 continue;
392 BmnHit* hit = hits[i];
393 // Use Z and X coordinates of hits to fit in ZX plane
394 Yi = hit->GetZ();
395 Xi = hit->GetX();
396 // h_Hits->Fill(hit->GetZ(), hit->GetX());
397 Zi = Xi * Xi + Yi * Yi;
398 Wi = 1.0 / hit->GetDx() / hit->GetDx();
399
400 Sx += Wi * Xi;
401 Sy += Wi * Yi;
402 Sz += Wi * Zi;
403 Sxx += Wi * Xi * Xi;
404 Sxy += Wi * Xi * Yi;
405 Syy += Wi * Yi * Yi;
406 Szx += Wi * Zi * Xi;
407 Szy += Wi * Zi * Yi;
408 }
409
410 Double_t C = ((Sz * Sx - Szx) / (Sxx - Sx * Sx) - (Sz * Sy - Szy) / (Sxy - Sx * Sy))
411 / ((Sxy - Sx * Sy) / (Sxx - Sx * Sx) - (Syy - Sy * Sy) / (Sxy - Sx * Sy));
412 Double_t B = ((Sz * Sx - Szx) - C * (Sxy - Sx * Sy)) / (Sxx - Sx * Sx);
413 Double_t D = -Sz - B * Sx - C * Sy;
414
415 Xc = -0.5 * B;
416 Zc = -0.5 * C;
417 R = Sqrt(0.25 * B * B + 0.25 * C * C - D);
418
419 return TVector3(Zc, Xc, R);
420}
421
422Double_t CalcTx(const BmnHit* h0, const BmnHit* h1, const BmnHit* h2)
423{
424 // function calculates Tx in point h0, so for TX_last use reverse order: CalcTx(h2, h1, h0)
425 TVector3 CircParZX = CircleBy3Hit(h0, h1, h2);
426 Double_t Xc = CircParZX.Y(); // x-coordinate of fit-circle center
427 Double_t Zc = CircParZX.X(); // z-coordinate of fit-circle center
428
429 return (-1.0 * (h0->GetZ() - Zc) / (h0->GetX() - Xc));
430}
431
432Float_t Dist(Float_t x1, Float_t y1, Float_t x2, Float_t y2)
433{
434 if (Sq(x1 - x2) + Sq(y1 - y2) <= 0.0) {
435 return 0.0;
436 } else {
437 return Sqrt(Sq(x1 - x2) + Sq(y1 - y2));
438 }
439}
440
441Float_t NewtonSolver(Float_t A0, Float_t A1, Float_t A2, Float_t A22)
442{
443 Double_t Dy = 0.0;
444 Double_t xnew = 0.0;
445 Double_t ynew = 0.0;
446 Double_t yold = 1e+11;
447 Double_t xold = 0.0;
448 const Double_t eps = 1e-12;
449 Int_t iter = 0;
450 const Int_t iterMax = 20;
451 do {
452 ynew = A0 + xnew * (A1 + xnew * (A2 + 4.0 * xnew * xnew));
453 if (fabs(ynew) > fabs(yold)) {
454 xnew = 0.0;
455 break;
456 }
457 Dy = A1 + xnew * (A22 + 16.0 * xnew * xnew);
458 xold = xnew;
459 xnew = xold - ynew / Dy;
460 iter++;
461 } while (Abs((xnew - xold) / xnew) > eps && iter < iterMax);
462
463 if (iter == iterMax - 1) {
464 xnew = 0.0;
465 }
466
467 return xnew;
468}
469
470TVector3 CircleBy3Hit(BmnTrack* track, const TClonesArray* arr)
471{
472 const Float_t nHits = track->GetNHits();
473 if (nHits < 3)
474 return TVector3(0.0, 0.0, 0.0);
475 BmnHit* hit0 = (BmnHit*)arr->At(track->GetHitIndex(0));
476 BmnHit* hit1 = (BmnHit*)arr->At(track->GetHitIndex(1));
477 BmnHit* hit2 = (BmnHit*)arr->At(track->GetHitIndex(2));
478
479 Float_t x1 = hit0->GetX();
480 Float_t z1 = hit0->GetZ();
481 Float_t x2 = hit1->GetX();
482 Float_t z2 = hit1->GetZ();
483 Float_t x3 = hit2->GetX();
484 Float_t z3 = hit2->GetZ();
485
486 Float_t x1_2 = x1 * x1;
487 Float_t z1_2 = z1 * z1;
488 Float_t x2_2 = x2 * x2;
489 Float_t z2_2 = z2 * z2;
490 Float_t x3_2 = x3 * x3;
491 Float_t z3_2 = z3 * z3;
492
493 Float_t B = ((x1 - x3) * (x2_2 + z2_2) + (x2 - x1) * (x3_2 + z3_2) + (x3 - x2) * (x1_2 + z1_2))
494 / (x1 * (z3 - z2) + x2 * (z1 - z3) + x3 * (z2 - z1));
495 Float_t A = ((x2_2 + z2_2) - (x1_2 + z1_2) - B * (z1 - z2)) / (x1 - x2);
496 Float_t C = -x1_2 - z1_2 - A * x1 - B * z1;
497
498 Float_t Xc = -A / 2;
499 Float_t Zc = -B / 2;
500 Float_t R = Sqrt(A * A + B * B - 4 * C) / 2;
501
502 return TVector3(Zc, Xc, R);
503}
504
505TVector3 Pol2By3Hit(BmnTrack* track, const TClonesArray* arr)
506{
507 const Int_t nHits = track->GetNHits();
508 if (nHits < 3)
509 return TVector3(0.0, 0.0, 0.0);
510 BmnHit* hit0 = (BmnHit*)arr->At(track->GetHitIndex(0));
511 BmnHit* hit1 = (BmnHit*)arr->At(track->GetHitIndex(1));
512 BmnHit* hit2 = (BmnHit*)arr->At(track->GetHitIndex(2));
513
514 Float_t x0 = hit0->GetX();
515 Float_t z0 = hit0->GetZ();
516 Float_t x1 = hit1->GetX();
517 Float_t z1 = hit1->GetZ();
518 Float_t x2 = hit2->GetX();
519 Float_t z2 = hit2->GetZ();
520
521 Float_t z0_2 = z0 * z0;
522 Float_t z1_2 = z1 * z1;
523 Float_t z2_2 = z2 * z2;
524
525 Float_t B = (x2 * (z1_2 - z0_2) + x1 * (z0_2 - z2_2) + x0 * (z2_2 - z1_2)) / ((z2 - z1) * (z1 - z0) * (z0 - z2));
526 Float_t A = ((x2 - x1) - B * (z2 - z1)) / (z2_2 - z1_2);
527 Float_t C = x0 - B * z0 - A * z0_2;
528
529 return TVector3(A, B, C);
530}
531
532void DrawHits(BmnTrack* track, const TClonesArray* arr)
533{
534 const Int_t nHits = track->GetNHits();
535 Float_t z[nHits];
536 Float_t x[nHits];
537 for (Int_t i = 0; i < nHits; ++i) {
538 z[i] = ((BmnHit*)arr->At(track->GetHitIndex(i)))->GetZ();
539 x[i] = ((BmnHit*)arr->At(track->GetHitIndex(i)))->GetX();
540 }
541 TCanvas* c = new TCanvas("c", "c", 1000, 1000);
542 TGraph* gr = new TGraph(nHits, z, x);
543 gr->Fit("pol2");
544 gr->Draw("AL*");
545 c->Update();
546 getchar();
547 delete gr;
548 delete c;
549}
550
551TVector3 CircleBy3Hit(const BmnHit* h0, const BmnHit* h1, const BmnHit* h2)
552{
553 Float_t x0 = h0->GetX();
554 Float_t z0 = h0->GetZ();
555 Float_t x1 = h1->GetX();
556 Float_t z1 = h1->GetZ();
557 Float_t x2 = h2->GetX();
558 Float_t z2 = h2->GetZ();
559
560 // cout << x1 << " " << z1 << " " << x2 << " " << z2 << " " << x3 << " " << z3 << " " << endl;
561
562 Float_t x0_2 = x0 * x0;
563 Float_t z0_2 = z0 * z0;
564 Float_t x1_2 = x1 * x1;
565 Float_t z1_2 = z1 * z1;
566 Float_t x2_2 = x2 * x2;
567 Float_t z2_2 = z2 * z2;
568 Float_t dx10 = x1 - x0;
569 Float_t dx21 = x2 - x1;
570 Float_t dz10 = z1 - z0;
571 Float_t dz21 = z2 - z1;
572 Float_t dx10_2 = x1_2 - x0_2;
573 Float_t dx21_2 = x2_2 - x1_2;
574 Float_t dz10_2 = z1_2 - z0_2;
575 Float_t dz21_2 = z2_2 - z1_2;
576
577 Float_t B = ((dx10_2 + dz10_2) / dx10 - (dx21_2 + dz21_2) / dx21) / (dz21 / dx21 - dz10 / dx10);
578 Float_t A = -(dx21_2 + dz21_2 + B * dz21) / dx21;
579 Float_t C = -x0_2 - z0_2 - A * x0 - B * z0;
580
581 Float_t Xc = -A / 2;
582 Float_t Zc = -B / 2;
583 Float_t R = Sqrt(A * A + B * B - 4 * C) / 2;
584
585 return TVector3(Zc, Xc, R);
586}
587// Пусть пока этот крокодил поживет здесь:)
588
589void fit_seg(Double_t* z_loc,
590 Double_t* rh_seg,
591 Double_t* rh_sigm_seg,
592 Double_t* par_ab,
593 Int_t skip_first,
594 Int_t skip_second)
595{
596 Double_t sqrt_2 = sqrt(2.);
597 // linear fit
598 Float_t A[4][4] = {{0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0, 0}}; // coef matrix
599 Float_t f[4] = {0}; // free coef
600 // Float_t sigm_sq[8] = {1,1,1,1,1,1,1,1};
601 Int_t h[8] = {0, 0, 0, 0, 0, 0, 0, 0};
602
603 for (Int_t i = 0; i < 4; i++)
604 par_ab[i] = 999;
605
606 for (Int_t i = 0; i < 8; i++) {
607 h[i] = 1;
608 // out1<<"setting h[i]"<<endl;
609 if (i == skip_first || i == skip_second || Abs(rh_seg[i] + 999.) < FLT_EPSILON) {
610 h[i] = 0;
611 }
612 } // i
613
614 A[0][0] = 2 * z_loc[0] * z_loc[0] * h[0] / rh_sigm_seg[0] + z_loc[4] * z_loc[4] * h[4] / rh_sigm_seg[4]
615 + z_loc[6] * z_loc[6] * h[6] / rh_sigm_seg[6] + 2 * z_loc[1] * z_loc[1] * h[1] / rh_sigm_seg[1]
616 + z_loc[5] * z_loc[5] * h[5] / rh_sigm_seg[5] + z_loc[7] * z_loc[7] * h[7] / rh_sigm_seg[7]; // aX_a
617
618 A[0][1] = 2 * z_loc[0] * h[0] / rh_sigm_seg[0] + z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[6] * h[6] / rh_sigm_seg[6]
619 + 2 * z_loc[1] * h[1] / rh_sigm_seg[1] + z_loc[5] * h[5] / rh_sigm_seg[5]
620 + z_loc[7] * h[7] / rh_sigm_seg[7]; // bX_a
621
622 A[0][2] = z_loc[6] * z_loc[6] * h[6] / rh_sigm_seg[6] - z_loc[4] * z_loc[4] * h[4] / rh_sigm_seg[4]
623 + z_loc[7] * z_loc[7] * h[7] / rh_sigm_seg[7] - z_loc[5] * z_loc[5] * h[5] / rh_sigm_seg[5]; // aY
624
625 A[0][3] = z_loc[6] * h[6] / rh_sigm_seg[6] - z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[7] * h[7] / rh_sigm_seg[7]
626 - z_loc[5] * h[5] / rh_sigm_seg[5]; // bY_a
627
628 // dChi2/d_b_x
629
630 A[1][0] = 2 * z_loc[0] * h[0] / rh_sigm_seg[0] + z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[6] * h[6] / rh_sigm_seg[6]
631 + 2 * z_loc[1] * h[1] / rh_sigm_seg[1] + z_loc[5] * h[5] / rh_sigm_seg[5]
632 + z_loc[7] * h[7] / rh_sigm_seg[7];
633
634 A[1][1] = 2 * h[0] / rh_sigm_seg[0] + 1 * h[4] / rh_sigm_seg[4] + 1 * h[6] / rh_sigm_seg[6]
635 + 2 * h[1] / rh_sigm_seg[1] + 1 * h[5] / rh_sigm_seg[5] + 1 * h[7] / rh_sigm_seg[7]; // bX_a
636
637 A[1][2] = z_loc[6] * h[6] / rh_sigm_seg[6] - z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[7] * h[7] / rh_sigm_seg[7]
638 - z_loc[5] * h[5] / rh_sigm_seg[5]; // aY_a
639
640 A[1][3] = 1 * h[7] / rh_sigm_seg[7] - 1 * h[5] / rh_sigm_seg[5] + 1 * h[6] / rh_sigm_seg[6]
641 - 1 * h[4] / rh_sigm_seg[4]; // bY_a
642
643 // dChi2/da_y
644
645 A[2][0] = z_loc[6] * z_loc[6] * h[6] / rh_sigm_seg[6] - z_loc[4] * z_loc[4] * h[4] / rh_sigm_seg[4]
646 + z_loc[7] * z_loc[7] * h[7] / rh_sigm_seg[7] - z_loc[5] * z_loc[5] * h[5] / rh_sigm_seg[5]; // aX_a
647
648 A[2][1] = z_loc[6] * h[6] / rh_sigm_seg[6] - z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[7] * h[7] / rh_sigm_seg[7]
649 - z_loc[5] * h[5] / rh_sigm_seg[5]; // bX_a
650
651 A[2][2] = 2 * z_loc[2] * z_loc[2] * h[2] / rh_sigm_seg[2] + z_loc[4] * z_loc[4] * h[4] / rh_sigm_seg[4]
652 + z_loc[6] * z_loc[6] * h[6] / rh_sigm_seg[6] + 2 * z_loc[3] * z_loc[3] * h[3] / rh_sigm_seg[3]
653 + z_loc[5] * z_loc[5] * h[5] / rh_sigm_seg[5] + z_loc[7] * z_loc[7] * h[7] / rh_sigm_seg[7]; // aY_a
654
655 A[2][3] = 2 * z_loc[2] * h[2] / rh_sigm_seg[2] + z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[6] * h[6] / rh_sigm_seg[6]
656 + 2 * z_loc[3] * h[3] / rh_sigm_seg[3] + z_loc[5] * h[5] / rh_sigm_seg[5]
657 + z_loc[7] * h[7] / rh_sigm_seg[7];
658
660
661 A[3][0] = z_loc[6] * h[6] / rh_sigm_seg[6] - z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[7] * h[7] / rh_sigm_seg[7]
662 - z_loc[5] * h[5] / rh_sigm_seg[5]; // aX_a
663
664 A[3][1] = 1 * h[6] / rh_sigm_seg[6] - 1 * h[4] / rh_sigm_seg[4] + 1 * h[7] / rh_sigm_seg[7]
665 - 1 * h[5] / rh_sigm_seg[5]; // bX_a
666
667 A[3][2] = 2 * z_loc[2] * h[2] / rh_sigm_seg[2] + z_loc[4] * h[4] / rh_sigm_seg[4] + z_loc[6] * h[6] / rh_sigm_seg[6]
668 + 2 * z_loc[3] * h[3] / rh_sigm_seg[3] + z_loc[5] * h[5] / rh_sigm_seg[5]
669 + z_loc[7] * h[7] / rh_sigm_seg[7]; // aY_a
670
671 A[3][3] = 2 * h[2] / rh_sigm_seg[2] + 1 * h[4] / rh_sigm_seg[4] + 1 * h[6] / rh_sigm_seg[6]
672 + 2 * h[3] / rh_sigm_seg[3] + 1 * h[5] / rh_sigm_seg[5] + 1 * h[7] / rh_sigm_seg[7]; // bY_a
673
674 // free coef
675
676 // dChi2/da_x
677
678 f[0] = 2 * z_loc[0] * rh_seg[0] * h[0] / rh_sigm_seg[0] + sqrt_2 * z_loc[6] * rh_seg[6] * h[6] / rh_sigm_seg[6]
679 - sqrt_2 * z_loc[4] * rh_seg[4] * h[4] / rh_sigm_seg[4] + 2 * z_loc[1] * rh_seg[1] * h[1] / rh_sigm_seg[1]
680 + sqrt_2 * z_loc[7] * rh_seg[7] * h[7] / rh_sigm_seg[7]
681 - sqrt_2 * z_loc[5] * rh_seg[5] * h[5] / rh_sigm_seg[5]; // j = nr of seg
682 // dChi2/db_x
683 f[1] = 2 * rh_seg[0] * h[0] / rh_sigm_seg[0] + sqrt_2 * rh_seg[6] * h[6] / rh_sigm_seg[6]
684 - sqrt_2 * rh_seg[4] * h[4] / rh_sigm_seg[4] + 2 * rh_seg[1] * h[1] / rh_sigm_seg[1]
685 + sqrt_2 * rh_seg[7] * h[7] / rh_sigm_seg[7] - sqrt_2 * rh_seg[5] * h[5] / rh_sigm_seg[5]; // j = nr of seg
686 // dChi2/da_y
687 f[2] = 2 * z_loc[2] * rh_seg[2] * h[2] / rh_sigm_seg[2] + sqrt_2 * z_loc[6] * rh_seg[6] * h[6] / rh_sigm_seg[6]
688 + sqrt_2 * z_loc[4] * rh_seg[4] * h[4] / rh_sigm_seg[4] + 2 * z_loc[3] * rh_seg[3] * h[3] / rh_sigm_seg[3]
689 + sqrt_2 * z_loc[7] * rh_seg[7] * h[7] / rh_sigm_seg[7]
690 + sqrt_2 * z_loc[5] * rh_seg[5] * h[5] / rh_sigm_seg[5]; // j = nr of seg
692 f[3] = 2 * rh_seg[2] * h[2] / rh_sigm_seg[2] + sqrt_2 * rh_seg[6] * h[6] / rh_sigm_seg[6]
693 + sqrt_2 * rh_seg[4] * h[4] / rh_sigm_seg[4] + 2 * rh_seg[3] * h[3] / rh_sigm_seg[3]
694 + sqrt_2 * rh_seg[7] * h[7] / rh_sigm_seg[7] + sqrt_2 * rh_seg[5] * h[5] / rh_sigm_seg[5]; // j = nr of seg
695
696 // inverse the matrix
697
698 /**** Gaussian algorithm for 4x4 matrix inversion ****/
699 Int_t i1, j1, k1, l1;
700 Double_t factor;
701 Double_t temp[4];
702 Double_t b[4][4];
703 Double_t A0[4][4];
704
705 for (i1 = 0; i1 < 4; i1++)
706 for (j1 = 0; j1 < 4; j1++)
707 A0[i1][j1] = A[i1][j1];
708
709 // Set b to I
710 for (i1 = 0; i1 < 4; i1++)
711 for (j1 = 0; j1 < 4; j1++)
712 if (i1 == j1)
713 b[i1][j1] = 1.0;
714 else
715 b[i1][j1] = 0.0;
716
717 for (i1 = 0; i1 < 4; i1++) {
718 for (j1 = i1 + 1; j1 < 4; j1++)
719 if (Abs(A[i1][i1]) < Abs(A[j1][i1])) {
720 for (l1 = 0; l1 < 4; l1++)
721 temp[l1] = A[i1][l1];
722 for (l1 = 0; l1 < 4; l1++)
723 A[i1][l1] = A[j1][l1];
724 for (l1 = 0; l1 < 4; l1++)
725 A[j1][l1] = temp[l1];
726 for (l1 = 0; l1 < 4; l1++)
727 temp[l1] = b[i1][l1];
728 for (l1 = 0; l1 < 4; l1++)
729 b[i1][l1] = b[j1][l1];
730 for (l1 = 0; l1 < 4; l1++)
731 b[j1][l1] = temp[l1];
732 }
733 factor = A[i1][i1];
734 for (j1 = 4 - 1; j1 > -1; j1--) {
735 b[i1][j1] /= factor;
736 A[i1][j1] /= factor;
737 }
738 for (j1 = i1 + 1; j1 < 4; j1++) {
739 factor = -A[j1][i1];
740 for (k1 = 0; k1 < 4; k1++) {
741 A[j1][k1] += A[i1][k1] * factor;
742 b[j1][k1] += b[i1][k1] * factor;
743 }
744 }
745 }
746 for (i1 = 3; i1 > 0; i1--) {
747 for (j1 = i1 - 1; j1 > -1; j1--) {
748 factor = -A[j1][i1];
749 for (k1 = 0; k1 < 4; k1++) {
750 A[j1][k1] += A[i1][k1] * factor;
751 b[j1][k1] += b[i1][k1] * factor;
752 }
753 }
754 }
755
756 Float_t sum;
757
758 /*Float_t A1[4][4] = {
759 {0, 0, 0, 0},
760 {0, 0, 0, 0},
761 {0, 0, 0, 0},
762 {0, 0, 0, 0}};
763*/
764 for (i1 = 0; i1 < 4; ++i1)
765 for (j1 = 0; j1 < 4; ++j1) {
766 sum = 0;
767
768 for (k1 = 0; k1 < 4; ++k1)
769 sum += A0[i1][k1] * b[k1][j1];
770 // A1[i1][j1] = sum;
771 }
772
773 for (i1 = 0; i1 < 4; i1++) {
774 par_ab[i1] = 0;
775 for (j1 = 0; j1 < 4; j1++) {
776 par_ab[i1] += b[i1][j1] * f[j1];
777 }
778 }
779}
780
781void Pol2Fit(BmnTrack* track, const TClonesArray* arr, Double_t& A, Double_t& B, Double_t& C, Int_t idSkip)
782{
783 const Float_t nHits = track->GetNHits();
784 TGraph* gr = new TGraph(nHits - 1);
785 Int_t iPoint = 0;
786 for (Int_t i = 0; i < nHits; ++i) {
787 if (i == idSkip)
788 continue;
789 BmnHit* hit = (BmnHit*)arr->At(track->GetHitIndex(i));
790 gr->SetPoint(iPoint++, hit->GetZ(), hit->GetX());
791 }
792 TFitResultPtr ptr = gr->Fit("pol2", "SQ");
793 delete gr;
794 A = ptr->Parameter(2);
795 B = ptr->Parameter(1);
796 C = ptr->Parameter(0);
797}
798
799TVector3 Pol2Fit(vector<BmnHit*> hits, Int_t idSkip)
800{
801 const Int_t nHits = hits.size();
802 TGraph gr(nHits - 1);
803 Int_t iPoint = 0;
804 for (Int_t i = 0; i < nHits; ++i) {
805 if (i == idSkip)
806 continue;
807 BmnHit* hit = hits[i];
808 gr.SetPoint(iPoint++, hit->GetZ(), hit->GetX());
809 }
810 TFitResultPtr ptr = gr.Fit("pol2", "SQ");
811 Double_t A = ptr->Parameter(2);
812 Double_t B = ptr->Parameter(1);
813 Double_t C = ptr->Parameter(0);
814
815 return TVector3(A, B, C);
816}
817
818TVector2 LineFit(vector<BmnHit*> hits, Int_t idSkip, TString type)
819{
820 // Weighted Least Square Method//
821 Float_t Xi = 0.0, Yi = 0.0; // coordinates of current track point
822 Float_t a = 0.0, b = 0.0; // parameters of line: y = a * x + b
823
824 Float_t Si = 0.0; // sigma
825 Float_t Wi = 0.0; // weight = 1 / sigma^2
826 Float_t SumW = 0.0; // sum of weights
827 Float_t SumWX = 0.0; // sum of (weight * x)
828 Float_t SumWY = 0.0; // sum of (weight * y)
829 Float_t SumWXY = 0.0; // sum of (weight * x * y)
830 Float_t SumWX2 = 0.0; // sum of (weight * x * x)
831
832 const Int_t nHits = hits.size();
833
834 for (Int_t i = 0; i < nHits; ++i) {
835 BmnHit* hit = hits[i];
836 if (i == idSkip)
837 continue;
838 if (type.Contains("XY")) {
839 Xi = hit->GetX();
840 Yi = hit->GetY();
841 Si = hit->GetDy();
842 } else if (type.Contains("ZX")) {
843 Xi = hit->GetZ();
844 Yi = hit->GetX();
845 Si = hit->GetDx();
846 } else if (type.Contains("ZY")) {
847 Xi = hit->GetZ();
848 Yi = hit->GetY();
849 Si = hit->GetDy();
850 }
851
852 if (Si == 0.0)
853 return TVector2(0.0, 0.0);
854
855 Wi = 1.0 / Si / Si;
856 SumW += Wi;
857 SumWXY += Wi * Xi * Yi;
858 SumWX += Wi * Xi;
859 SumWX2 += Wi * Xi * Xi;
860 SumWY += Wi * Yi;
861 }
862
863 a = (SumW * SumWXY - SumWX * SumWY) / (SumW * SumWX2 - SumWX * SumWX);
864 b = (SumWX2 * SumWY - SumWX * SumWXY) / (SumW * SumWX2 - SumWX * SumWX);
865
866 return TVector2(a, b);
867}
868
869vector<Double_t> dist(vector<Double_t> qp, Double_t mu)
870{
871 vector<Double_t> res;
872 for (size_t iEle = 0; iEle < qp.size(); iEle++)
873 res.push_back(Abs(qp[iEle] - mu));
874
875 return res;
876}
877
878vector<Double_t> W(vector<Double_t> dist, Double_t sig)
879{
880 vector<Double_t> res;
881 const Int_t C = 10;
882 for (size_t iEle = 0; iEle < dist.size(); iEle++) {
883 Double_t w = (dist[iEle] > C * sig) ? 0. : Power(1 - Power(dist[iEle] / C / sig, 2), 2);
884 res.push_back(w);
885 }
886
887 return res;
888}
889
890Double_t Sigma(vector<Double_t> dist, vector<Double_t> w)
891{
892 if (dist.size() != w.size())
893 throw;
894
895 Double_t WiDi2Sum = 0.;
896 Double_t WiSum = 0.;
897
898 for (size_t iEle = 0; iEle < dist.size(); iEle++) {
899 WiSum += w[iEle];
900 WiDi2Sum += w[iEle] * dist[iEle] * dist[iEle];
901 }
902
903 return Sqrt(WiDi2Sum / WiSum);
904}
905
906Double_t Mu(vector<Double_t> qp, vector<Double_t> w)
907{
908 if (qp.size() != w.size())
909 throw;
910
911 Double_t WiQpSum = 0.;
912 Double_t WiSum = 0.;
913
914 for (size_t iEle = 0; iEle < qp.size(); iEle++) {
915 WiSum += w[iEle];
916 WiQpSum += w[iEle] * qp[iEle];
917 }
918
919 return WiQpSum / WiSum;
920}
921
922// void DrawBar(Int_t iEv, Int_t nEv) {
923// cout.flush();
924// Float_t progress = iEv * 1.0 / nEv;
925// Int_t barWidth = 70;
926// printf(ANSI_COLOR_BLUE "[");
927//
928// Int_t pos = barWidth * progress;
929// for (Int_t i = 0; i < barWidth; ++i) {
930// if (i < pos) printf("=");
931// else if (i == pos) printf(">");
932// else printf(" ");
933// }
934//
935// printf("] " ANSI_COLOR_RESET);
936// printf(ANSI_COLOR_RED "%d%%\r" ANSI_COLOR_RESET, Int_t(progress * 100.0 + 0.5));
937// cout.flush();
938// }
939
940void DrawBar(UInt_t iEv, UInt_t nEv)
941{
942 if (nEv < 100)
943 return;
944 if ((iEv % (nEv / 100)) > 100.0 / nEv)
945 return;
946 Float_t progress = iEv * 1.0 / nEv;
947 Int_t barWidth = 70;
948 Int_t pos = barWidth * progress;
949 cout.flush();
950 for (Int_t i = 0; i < barWidth; ++i) {
951 if (i <= pos)
952 printf(ANSI_COLOR_BLUE_BG " " ANSI_COLOR_RESET);
953 else
954 printf(ANSI_COLOR_YELLOW_BG " " ANSI_COLOR_RESET);
955 }
956
957 printf(ANSI_COLOR_RED "[%d%%]\r" ANSI_COLOR_RESET, Int_t(progress * 100.0 + 0.5));
958 cout.flush();
959}
960
961void DrawBar(Long64_t iEv, Long64_t nEv)
962{
963 cout.flush();
964 Float_t progress = iEv * 1.0 / nEv;
965 Int_t barWidth = 70;
966
967 Int_t pos = barWidth * progress;
968 for (Int_t i = 0; i < barWidth; ++i) {
969 if (i <= pos)
970 printf(ANSI_COLOR_BLUE_BG " " ANSI_COLOR_RESET);
971 else
972 printf(ANSI_COLOR_YELLOW_BG " " ANSI_COLOR_RESET);
973 }
974
975 printf(ANSI_COLOR_RED "[%d%%]\r" ANSI_COLOR_RESET, Int_t(progress * 100.0 + 0.5));
976 cout.flush();
977}
978
979Double_t GetVZByTwoStraightTracks(BmnTrack* tr0, BmnTrack* tr1, Double_t& dist)
980{
981 // method returns z coordinate of two tracks crossing point in case of absent magnetic field
982 // dist is filling by distance between this two tracks
983 FairTrackParam* p0 = tr0->GetParamFirst();
984 FairTrackParam* p1 = tr1->GetParamFirst();
985 Double_t axl = p0->GetTx();
986 Double_t ayl = p0->GetTy();
987 Double_t axr = p1->GetTx();
988 Double_t ayr = p1->GetTy();
989 Double_t bxl = p0->GetX() - axl * p0->GetZ();
990 Double_t byl = p0->GetY() - ayl * p0->GetZ();
991 Double_t bxr = p1->GetX() - axr * p1->GetZ();
992 Double_t byr = p1->GetY() - ayr * p1->GetZ();
993 Double_t vz = -((bxl - bxr) * (axl - axr) + (byl - byr) * (ayl - ayr))
994 / ((axl - axr) * (axl - axr) + (ayl - ayr) * (ayl - ayr));
995 dist = Sqrt(Sq(axl * vz + bxl - axr * vz - bxr) + Sq(ayl * vz + byl - ayr * vz - byr));
996 return vz;
997}
998
999Double_t GetVzByVectorStraightTracks(vector<BmnTrack> tr, Double_t& dist)
1000{
1001 if (tr.size() < 2) {
1002 dist = 1000;
1003 return 1000;
1004 }
1005 Double_t A = 0;
1006 Double_t B = 0;
1007 for (size_t i = 0; i < tr.size() - 1; ++i) {
1008 FairTrackParam* pi = tr[i].GetParamFirst();
1009 Double_t axi = pi->GetTx();
1010 Double_t ayi = pi->GetTy();
1011 Double_t bxi = pi->GetX() - axi * pi->GetZ();
1012 Double_t byi = pi->GetY() - ayi * pi->GetZ();
1013 for (size_t j = i + 1; j < tr.size(); ++j) {
1014 FairTrackParam* pj = tr[j].GetParamFirst();
1015 Double_t axj = pj->GetTx();
1016 Double_t ayj = pj->GetTy();
1017 Double_t bxj = pj->GetX() - axj * pj->GetZ();
1018 Double_t byj = pj->GetY() - ayj * pj->GetZ();
1019 A += (Sq(axi - axj) + Sq(ayi - ayj));
1020 B += ((bxi - bxj) * (axi - axj) + (byi - byj) * (ayi - ayj));
1021 }
1022 }
1023 Double_t vz = -B / A;
1024
1025 Double_t SumDist = 0;
1026 Int_t cntr = 0;
1027 for (size_t i = 0; i < tr.size() - 1; ++i) {
1028 FairTrackParam* pi = tr[i].GetParamFirst();
1029 Double_t axi = pi->GetTx();
1030 Double_t ayi = pi->GetTy();
1031 Double_t bxi = pi->GetX() - axi * pi->GetZ();
1032 Double_t byi = pi->GetY() - ayi * pi->GetZ();
1033 for (size_t j = i + 1; j < tr.size(); ++j) {
1034 FairTrackParam* pj = tr[j].GetParamFirst();
1035 Double_t axj = pj->GetTx();
1036 Double_t ayj = pj->GetTy();
1037 Double_t bxj = pj->GetX() - axj * pj->GetZ();
1038 Double_t byj = pj->GetY() - ayj * pj->GetZ();
1039 SumDist += Sqrt(Sq(axi * vz + bxi - axj * vz - bxj) + Sq(ayi * vz + byi - ayj * vz - byj));
1040 cntr++;
1041 }
1042 }
1043 dist = SumDist / cntr;
1044
1045 return vz;
1046}
1047
1048void UpdateTrackParam(FairTrackParam* initPar, const FairTrackParam* detPar, Double_t& chiSq)
1049{
1050
1051 // pardet(5) - det track parameters
1052 // covdet(5) - det track parameters covariation matrix
1053 // parinit(5) - init track parameters
1054 // covinit(5) - init track parameters covariation matrix
1055 // parinitdet(5) - init+det track parameters
1056 // covinitdet(5) - init+det track parameters covariation matrix
1057 // detPar - det track param
1058 // initTr - Global track param
1059
1060 TMatrixDSym covdet(5), icovdet(5);
1061 covdet *= 0.;
1062 for (int ir = 0; ir < 5; ir++) {
1063 covdet(ir, ir) = detPar->GetCovariance(ir, ir);
1064 if (detPar->GetCovariance(ir, ir) == 0)
1065 covdet(ir, ir) = 0.00001;
1066 }
1067 covdet(4, 4) = 10000.;
1068
1069 icovdet = covdet;
1070 icovdet.Invert();
1071
1072 TMatrixDSym covinit(5), icovinit(5);
1073 covinit *= 0.;
1074 for (int ir = 0; ir < 5; ir++)
1075 for (int ic = 0; ic < 5; ic++)
1076 covinit(ir, ic) = initPar->GetCovariance(ir, ic);
1077 icovinit = covinit;
1078 icovinit.Invert();
1079
1080 TMatrixDSym covinitdet(5), icovinitdet(5);
1081 icovinitdet *= 0.;
1082 icovinitdet += icovinit;
1083 icovinitdet += icovdet;
1084 covinitdet = icovinitdet;
1085 covinitdet.Invert();
1086
1087 TVectorD suminitdet(5), parinitdet(5);
1088 TVectorD parinit(5), pardet(5);
1089 pardet(0) = detPar->GetX();
1090 pardet(1) = detPar->GetY();
1091 pardet(2) = detPar->GetTx();
1092 pardet(3) = detPar->GetTy();
1093 pardet(4) = 1.;
1094 parinit(0) = initPar->GetX();
1095 parinit(1) = initPar->GetY();
1096 parinit(2) = initPar->GetTx();
1097 parinit(3) = initPar->GetTy();
1098 parinit(4) = initPar->GetQp();
1099 suminitdet *= 0.;
1100 suminitdet += icovinit * parinit;
1101 suminitdet += icovdet * pardet;
1102 parinitdet = covinitdet * suminitdet;
1103
1104 // Set Updated parameters
1105
1106 initPar->SetX(parinitdet(0));
1107 initPar->SetY(parinitdet(1));
1108 initPar->SetTx(parinitdet(2));
1109 initPar->SetTy(parinitdet(3));
1110 initPar->SetQp(parinitdet(4));
1111 initPar->SetCovMatrix(covinitdet);
1112
1113 // Calculate chiSq
1114
1115 TMatrixD dP2(5, 1);
1116
1117 dP2(0, 0) = parinitdet(0) - pardet(0);
1118 dP2(1, 0) = parinitdet(1) - pardet(1);
1119 dP2(2, 0) = parinitdet(2) - pardet(2);
1120 dP2(3, 0) = parinitdet(3) - pardet(3);
1121 dP2(4, 0) = parinitdet(4) - pardet(4);
1122
1123 TMatrixD dP2_T(1, 5);
1124
1125 dP2_T.Transpose(dP2);
1126
1127 //-----------------------
1128
1129 TMatrixD dP1(5, 1);
1130
1131 dP1(0, 0) = parinitdet(0) - parinit(0);
1132 dP1(1, 0) = parinitdet(1) - parinit(1);
1133 dP1(2, 0) = parinitdet(2) - parinit(2);
1134 dP1(3, 0) = parinitdet(3) - parinit(3);
1135 dP1(4, 0) = parinitdet(4) - parinit(4);
1136
1137 TMatrixD dP1_T(1, 5);
1138
1139 dP1_T.Transpose(dP1);
1140
1141 TMatrixD chi(1, 1);
1142 chi = dP2_T * icovdet * dP2 + dP1_T * icovinit * dP1;
1143 // chi = dP2_T*icovdet*dP2 + dP1_T*icovinitdet*dP1;
1144
1145 chiSq = chi(0, 0);
1146}
LineFit
TVector3 LineFit(BmnTrack *track, const TClonesArray *arr, TString type)
Definition BmnMath.cxx:113
DrawHits
void DrawHits(BmnTrack *track, const TClonesArray *arr)
Definition BmnMath.cxx:532
DrawBar
void DrawBar(UInt_t iEv, UInt_t nEv)
Definition BmnMath.cxx:940
CircleBy3Hit
TVector3 CircleBy3Hit(BmnTrack *track, const TClonesArray *arr)
Definition BmnMath.cxx:470
Mu
Double_t Mu(vector< Double_t > qp, vector< Double_t > w)
Definition BmnMath.cxx:906
GetVZByTwoStraightTracks
Double_t GetVZByTwoStraightTracks(BmnTrack *tr0, BmnTrack *tr1, Double_t &dist)
Definition BmnMath.cxx:979
CircleFit
TVector3 CircleFit(BmnTrack *track, const TClonesArray *arr, Double_t &chi2)
Definition BmnMath.cxx:312
CalcTx
Double_t CalcTx(const BmnHit *h0, const BmnHit *h1, const BmnHit *h2)
Definition BmnMath.cxx:422
UpdateTrackParam
void UpdateTrackParam(FairTrackParam *initPar, const FairTrackParam *detPar, Double_t &chiSq)
Definition BmnMath.cxx:1048
dist
vector< Double_t > dist(vector< Double_t > qp, Double_t mu)
Definition BmnMath.cxx:869
NewtonSolver
Float_t NewtonSolver(Float_t A0, Float_t A1, Float_t A2, Float_t A22)
Definition BmnMath.cxx:441
IsParCorrect
Bool_t IsParCorrect(const FairTrackParam *par, const Bool_t isField)
Definition BmnMath.cxx:59
Pol2Fit
void Pol2Fit(BmnTrack *track, const TClonesArray *arr, Double_t &A, Double_t &B, Double_t &C, Int_t idSkip)
Definition BmnMath.cxx:781
LineFitBy3Hits
TVector3 LineFitBy3Hits(const BmnHit *h0, const BmnHit *h1, const BmnHit *h2)
Definition BmnMath.cxx:209
Sigma
Double_t Sigma(vector< Double_t > dist, vector< Double_t > w)
Definition BmnMath.cxx:890
Dist
Float_t Dist(Float_t x1, Float_t y1, Float_t x2, Float_t y2)
Definition BmnMath.cxx:432
fit_seg
void fit_seg(Double_t *z_loc, Double_t *rh_seg, Double_t *rh_sigm_seg, Double_t *par_ab, Int_t skip_first, Int_t skip_second)
Definition BmnMath.cxx:589
W
vector< Double_t > W(vector< Double_t > dist, Double_t sig)
Definition BmnMath.cxx:878
NumericalRootFinder
Float_t NumericalRootFinder(TF1 f, Float_t left, Float_t right)
Definition BmnMath.cxx:96
GetVzByVectorStraightTracks
Double_t GetVzByVectorStraightTracks(vector< BmnTrack > tr, Double_t &dist)
Definition BmnMath.cxx:999
Pol2By3Hit
TVector3 Pol2By3Hit(BmnTrack *track, const TClonesArray *arr)
Definition BmnMath.cxx:505
z0
const Float_t z0
Definition BmnMwpcHit.cxx:6
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition P4_F32vec4.h:34
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
i
int i
Definition P4_F32vec4.h:22
f
float f
Definition P4_F32vec4.h:21
ANSI_COLOR_RED
#define ANSI_COLOR_RED
Definition BmnMath.h:16
ANSI_COLOR_YELLOW_BG
#define ANSI_COLOR_YELLOW_BG
Definition BmnMath.h:19
ANSI_COLOR_BLUE_BG
#define ANSI_COLOR_BLUE_BG
Definition BmnMath.h:20
ANSI_COLOR_RESET
#define ANSI_COLOR_RESET
Definition BmnMath.h:18
BmnFieldMap::IsFieldOff
Bool_t IsFieldOff() const
Definition BmnFieldMap.h:93
BmnTrack::GetParamFirst
FairTrackParam * GetParamFirst()
Definition BmnTrack.h:72
BmnTrack::GetNHits
Int_t GetNHits() const
Definition BmnTrack.h:44
BmnTrack::GetHitIndex
Int_t GetHitIndex(Int_t iHit) const
Definition BmnTrack.h:48
lit::NDF
Int_t NDF(const BmnTrack *track)
Definition BmnMath.cxx:47
lit::ChiSq
Float_t ChiSq(const FairTrackParam *par, const BmnHit *hit)
Definition BmnMath.cxx:29