]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSNeuralTrack.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Now the full chain includes raw data.
[u/mrichter/AliRoot.git] / ITS / AliITSNeuralTrack.cxx
CommitLineData
f77f13c8 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17
853a0f19 18// AliITSNeuralTrack
19//
20// The format of output data from Neural Tracker
21// It can export data in the format of AliITSIOTrack
22// (V1) tracking.
23// Compatibility adaptation to V2 tracking is on the way.
24// Author: A. Pulvirenti
25
b9d722bc 26#include <Riostream.h>
853a0f19 27//#include <cstdlib>
28//#include <cstring>
b9d722bc 29
853a0f19 30//#include <TObject.h>
31//#include <TROOT.h>
b9d722bc 32#include <TMath.h>
33#include <TString.h>
853a0f19 34//#include <TObjArray.h>
35//#include <TH1.h>
b9d722bc 36#include <TMatrixD.h>
49370114 37#if ROOT_VERSION_CODE >= 262146
38#include <TMatrixDEigen.h>
39#endif
b9d722bc 40
41//#include "AliITSVertex.h"
42#include "AliITSIOTrack.h"
43#include "AliITSNeuralPoint.h"
44
45#include "AliITSNeuralTrack.h"
46
47
48
49ClassImp(AliITSNeuralTrack)
50//
51//
52//
53AliITSNeuralTrack::AliITSNeuralTrack() : fMatrix(5,5), fVertex()
54{
853a0f19 55 // Default constructor
56
b9d722bc 57 Int_t i;
58
59 fMass = 0.1396; // default assumption: pion
60 fField = 2.0; // default assumption: B = 0.4 Tesla
61
62 fXC = fYC = fR = fC = 0.0;
63 fTanL = fG0 = fDt = fDz = 0.0;
64 fStateR = fStatePhi = fStateZ = fChi2 = fNSteps = 0.0;
65
66 fLabel = 0;
67 fCount = 0;
68 for (i = 0; i < 6; i++) fPoint[i] = 0;
69
70 fVertex.X() = 0.0;
71 fVertex.Y() = 0.0;
72 fVertex.Z() = 0.0;
73 fVertex.ErrX() = 0.0;
74 fVertex.ErrY() = 0.0;
75 fVertex.ErrZ() = 0.0;
76}
77//
78//
79//
853a0f19 80AliITSNeuralTrack::AliITSNeuralTrack(AliITSNeuralTrack &track)
81: TObject((TObject&)track), fMatrix(5,5), fVertex()
82{
83// copy constructor
84
85 Int_t i;
86
87 fMass = 0.1396; // default assumption: pion
88 fField = 2.0; // default assumption: B = 0.4 Tesla
89
90 fXC = fYC = fR = fC = 0.0;
91 fTanL = fG0 = fDt = fDz = 0.0;
92 fStateR = fStatePhi = fStateZ = fChi2 = fNSteps = 0.0;
93
94 fLabel = 0;
95 fCount = 0;
96 for (i = 0; i < 6; i++) fPoint[i] = track.fPoint[i];
97
98 fVertex.X() = 0.0;
99 fVertex.Y() = 0.0;
100 fVertex.Z() = 0.0;
101 fVertex.ErrX() = 0.0;
102 fVertex.ErrY() = 0.0;
103 fVertex.ErrZ() = 0.0;
104}
105//
106//
107//
b9d722bc 108AliITSNeuralTrack::~AliITSNeuralTrack()
109{
110 Int_t l;
111 for (l = 0; l < 6; l++) fPoint[l] = 0;
112}
113//
114//
115//
116void AliITSNeuralTrack::AssignLabel()
853a0f19 117{
b9d722bc 118// Assigns a GEANT label to the found track.
119// Every cluster has up to three labels (it can have less). Then each label is
120// recorded for each point. Then, counts are made to check if some of the labels
121// appear more than once. Finally, the label which appears most times is assigned
122// to the track in the field fLabel.
123// The number of points containing that label is counted in the fCount data-member.
853a0f19 124
b9d722bc 125 Bool_t found;
126 Int_t i, j, l, lab, max = 0;
127
128 // We have up to 6 points for 3 labels each => up to 18 possible different values
129 Int_t idx[18] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
130 Int_t count[18] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
131
132 for (l = 0; l < 6; l++) {
133 if (!fPoint[l]) continue;
134 // Sometimes the same label appears two times in the same recpoint.
135 // With these if statements, such problem is solved by turning
136 // one of them to -1.
137 if (fPoint[l]->GetLabel(1) >= 0 && fPoint[l]->GetLabel(1) == fPoint[l]->GetLabel(0))
138 fPoint[l]->SetLabel(1, -1);
139 if (fPoint[l]->GetLabel(2) >= 0 && fPoint[l]->GetLabel(2) == fPoint[l]->GetLabel(0))
140 fPoint[l]->SetLabel(2, -1);
141 if (fPoint[l]->GetLabel(2) >= 0 && fPoint[l]->GetLabel(2) == fPoint[l]->GetLabel(1))
142 fPoint[l]->SetLabel(2, -1);
143 for (i = 0; i < 3; i++) {
144 lab = fPoint[l]->GetLabel(i);
145 if (lab < 0) continue;
146 found = kFALSE;
147 for (j = 0; j < max; j++) {
148 if (idx[j] == lab) {
149 count[j]++;
150 found = kTRUE;
151 }
152 }
153 if(!found) {
154 max++;
155 idx[max - 1] = lab;
156 count[max - 1] = 1;
157 }
158 }
159 }
160
161 j = 0, max = count[0];
162 for (i = 0; i < 18; i++) {
163 if (count[i] > max) {
164 j = i;
165 max = count[i];
166 }
167 }
168 fLabel = idx[j];
169 fCount = count[j];
170}
171//
172//
173//
174void AliITSNeuralTrack::CleanSlot(Int_t i, Bool_t del)
853a0f19 175{
b9d722bc 176// Removes a point from the corresponding layer slot in the found track.
177// If the argument is TRUE, the point object is also deleted from heap.
853a0f19 178
b9d722bc 179 if (i >= 0 && i < 6) {
180 if (del) delete fPoint[i];
181 fPoint[i] = 0;
182 }
183}
184//
185//
186//
187void AliITSNeuralTrack::GetModuleData(Int_t layer, Int_t &mod, Int_t &pos)
853a0f19 188{
b9d722bc 189// Returns the point coordinates according to the TreeR philosophy in galice.root files
190// that consist in the module number (mod) and the position in the TClonesArray of
191// the points reconstructed in that module for the run being examined.
853a0f19 192
b9d722bc 193 if (layer < 0 || layer > 5) {
194 Error("GetModuleData", "Layer out of range: %d", layer);
195 return;
196 }
197 mod = fPoint[layer]->GetModule();
198 pos = fPoint[layer]->GetIndex();
199}
200//
201//
202//
203void AliITSNeuralTrack::Insert(AliITSNeuralPoint *point)
853a0f19 204{
b9d722bc 205// A trivial method to insert a point in the tracks;
206// the point is inserted to the slot corresponding to its ITS layer.
853a0f19 207
b9d722bc 208 if (!point) return;
209
210 Int_t layer = point->GetLayer();
211 if (layer < 0 || layer > 6) {
212 Error("Insert", "Layer index %d out of range", layer);
213 return;
214 }
215
216 fPoint[layer] = point;
217}
218//
219//
220//
221Int_t AliITSNeuralTrack::OccupationMask()
853a0f19 222{
b9d722bc 223// Returns a byte which maps the occupied slots.
224// Each bit represents a layer going from the less significant on.
853a0f19 225
b9d722bc 226 Int_t i, check, mask = 0;
227 for (i = 0; i < 6; i++) {
228 check = 1 << i;
229 if (fPoint[i]) mask |= check;
230 }
231 return mask;
232}
233//
234//
235//
236void AliITSNeuralTrack::PrintLabels()
853a0f19 237{
b9d722bc 238// Prints the results of the AssignLabel() method, together with
239// the GEANT labels assigned to each point, in order to evaluate
240// how the assigned label is distributed among points.
853a0f19 241
b9d722bc 242 cout << "Assigned label = " << fLabel << " -- counted " << fCount << " times: " << endl << endl;
243 for (Int_t i = 0; i < 6; i++) {
244 cout << "Point #" << i + 1 << " --> ";
245 if (fPoint[i]) {
246 cout << "labels = " << fPoint[i]->GetLabel(0) << ", ";
247 cout << fPoint[i]->GetLabel(1) << ", ";
248 cout << fPoint[i]->GetLabel(2) << endl;
249 }
250 else {
251 cout << "not assigned" << endl;
252 }
253 }
254 cout << endl;
255}
256//
257//
258//
259Bool_t AliITSNeuralTrack::AddEL(Int_t layer, Double_t sign)
260{
853a0f19 261// Calculates the correction for energy loss
262
b9d722bc 263 Double_t width = 0.0;
264 switch (layer) {
265 case 0: width = 0.00260 + 0.00283; break;
266 case 1: width = 0.0180; break;
267 case 2: width = 0.0094; break;
268 case 3: width = 0.0095; break;
269 case 4: width = 0.0091; break;
270 case 5: width = 0.0087; break;
271 default:
272 Error("AddEL", "Layer value %d out of range!", layer);
273 return kFALSE;
274 }
275 width *= 1.7;
276
277 if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) {
278 width += 0.012;
279 }
280
281 Double_t invSqCosL = 1. + fTanL * fTanL; // = 1 / (cos(lambda)^2) = 1 + tan(lambda)^2
282 Double_t invCosL = TMath::Sqrt(invSqCosL); // = 1 / cos(lambda)
283 Double_t pt = GetPt(); // = transverse momentum
284 Double_t p2 = pt *pt * invSqCosL; // = square modulus of momentum
285 Double_t energy = TMath::Sqrt(p2 + fMass * fMass); // = energy
286 Double_t beta2 = p2 / (p2 + fMass * fMass); // = (v / c) ^ 2
287 if (beta2 == 0.0) {
288 printf("Anomaly in AddEL: pt=%8.6f invSqCosL=%8.6f fMass=%8.7f --> beta2 = %8.7f\n", pt, invSqCosL, fMass, beta2);
289 return kFALSE;
290 }
291
292 Double_t dE = 0.153 / beta2 * (log(5940. * beta2 / (1. - beta2)) - beta2) * width * 21.82 * invCosL;
293 dE = sign * dE * 0.001;
294
295 energy += dE;
296 p2 = energy * energy - fMass * fMass;
297 pt = TMath::Sqrt(p2) / invCosL;
298 if (fC < 0.) pt = -pt;
299 fC = (0.299792458 * 0.2 * fField) / (pt * 100.);
300
301 return kTRUE;
302}
303//
304//
305//
306Bool_t AliITSNeuralTrack::AddMS(Int_t layer)
307{
853a0f19 308// Calculates the noise perturbation due to multiple scattering
309
b9d722bc 310 Double_t width = 0.0;
311 switch (layer) {
312 case 0: width = 0.00260 + 0.00283; break;
313 case 1: width = 0.0180; break;
314 case 2: width = 0.0094; break;
315 case 3: width = 0.0095; break;
316 case 4: width = 0.0091; break;
317 case 5: width = 0.0087; break;
318 default:
319 Error("AddEL", "Layer value %d out of range!", layer);
320 return kFALSE;
321 }
322 width *= 1.7;
323
324 if((layer == 5) && (fStatePhi < 0.174 || fStatePhi > 6.100 || (fStatePhi > 2.960 && fStatePhi < 3.31))) {
325 width += 0.012;
326 }
327
328 Double_t cosL = TMath::Cos(TMath::ATan(fTanL));
329 Double_t halfC = fC / 2.;
330 Double_t q20 = 1. / (cosL * cosL);
331 Double_t q30 = fC * fTanL;
332
333 Double_t q40 = halfC * (fStateR * fStateR - fDt * fDt) / (1. + 2. * halfC * fDt);
334 Double_t dd = fDt + halfC * fDt * fDt - halfC * fStateR * fStateR;
335 Double_t dprova = fStateR * fStateR - dd * dd;
336 Double_t q41 = 0.;
337 if(dprova > 0.) q41 = -1. / cosL * TMath::Sqrt(dprova) / (1. + 2. * halfC *fDt);
338
339 Double_t p2 = (GetPt()*GetPt()) / (cosL * cosL);
340 Double_t beta2 = p2 / (p2 + fMass * fMass);
341 Double_t theta2 = 14.1 * 14.1 / (beta2 * p2 * 1.e6) * (width / TMath::Abs(cosL));
342
343 fMatrix(2,2) += theta2 * (q40 * q40 + q41 * q41);
344 fMatrix(3,2) += theta2 * q20 * q40;
345 fMatrix(2,3) += theta2 * q20 * q40;
346 fMatrix(3,3) += theta2 * q20 * q20;
347 fMatrix(4,2) += theta2 * q30 * q40;
348 fMatrix(2,4) += theta2 * q30 * q40;
349 fMatrix(4,3) += theta2 * q30 * q20;
350 fMatrix(3,4) += theta2 * q30 * q20;
351 fMatrix(4,4) += theta2 * q30 * q30;
352
353 return kTRUE;
354}
355//
356//
357//
358Int_t AliITSNeuralTrack::PropagateTo(Double_t rk)
359{
360 // Propagation method.
361 // Changes the state vector according to a new radial position
362 // which is specified by the passed 'r' value (in cylindircal coordinates).
363 // The covariance matrix is also propagated (and enlarged) according to
364 // the FCFt technique, where F is the jacobian of the new parameters
365 // w.r.t. their old values.
366 // The option argument forces the method to add also the energy loss
367 // and the multiple scattering effects, which respectively have the effect
368 // of changing the curvature and widening the covariance matrix.
369
b9f05b32 370 if (rk < TMath::Abs(fDt)) {
b9d722bc 371 Error("PropagateTo", Form("Impossible propagation to r (=%17.15g) < Dt (=%17.15g)", rk, fDt));
372 return 0;
373 }
374
375 Double_t duepi = 2. * TMath::Pi();
376 Double_t rkm1 = fStateR;
377 Double_t aAk = ArgPhi(rk), aAkm1 = ArgPhi(rkm1);
378 Double_t ak = ArgZ(rk), akm1 = ArgZ(rkm1);
379
380 fStatePhi += TMath::ASin(aAk) - TMath::ASin(aAkm1);
381 if(fStatePhi > duepi) fStatePhi -= duepi;
382 if(fStatePhi < 0.) fStatePhi += duepi;
383
384 Double_t halfC = 0.5 * fC;
385 fStateZ += fTanL / halfC * (TMath::ASin(ak)-TMath::ASin(akm1));
386
387 Double_t bk = ArgB(rk), bkm1 = ArgB(rkm1);
388 Double_t ck = ArgC(rk), ckm1 = ArgC(rkm1);
389
390 Double_t f02 = ck / TMath::Sqrt(1. - aAk * aAk) - ckm1 / TMath::Sqrt(1. - aAkm1 * aAkm1);
391 Double_t f04 = bk / TMath::Sqrt(1. - aAk * aAk) - bkm1 / TMath::Sqrt(1. - aAkm1 * aAkm1);
392 Double_t f12 = fTanL * fDt * (1. / rk - 1. / rkm1);
393 Double_t f13 = rk - rkm1;
394
395 Double_t c00 = fMatrix(0,0);
396 Double_t c10 = fMatrix(1,0);
397 Double_t c11 = fMatrix(1,1);
398 Double_t c20 = fMatrix(2,0);
399 Double_t c21 = fMatrix(2,1);
400 Double_t c22 = fMatrix(2,2);
401 Double_t c30 = fMatrix(3,0);
402 Double_t c31 = fMatrix(3,1);
403 Double_t c32 = fMatrix(3,2);
404 Double_t c33 = fMatrix(3,3);
405 Double_t c40 = fMatrix(4,0);
406 Double_t c41 = fMatrix(4,1);
407 Double_t c42 = fMatrix(4,2);
408 Double_t c43 = fMatrix(4,3);
409 Double_t c44 = fMatrix(4,4);
410
411 Double_t r10 = c10 + c21*f02 + c41*f04;
412 Double_t r20 = c20 + c22*f02 + c42*f04;
413 Double_t r30 = c30 + c32*f02 + c43*f04;
414 Double_t r40 = c40 + c42*f02 + c44*f04;
415 Double_t r21 = c21 + c22*f12 + c32*f13;
416 Double_t r31 = c31 + c32*f12 + c33*f13;
417 Double_t r41 = c41 + c42*f12 + c43*f13;
418
419 fMatrix(0,0) = c00 + c20*f02 + c40*f04 + f02*r20 + f04*r40;
420 fMatrix(1,0) = fMatrix(0,1) = r10 + f12*r20 + f13*r30;
421 fMatrix(1,1) = c11 + c21*f12 + c31*f13 + f12*r21 + f13*r31;
422 fMatrix(2,0) = fMatrix(0,2) = r20;
423 fMatrix(2,1) = fMatrix(1,2) = r21;
424 fMatrix(3,0) = fMatrix(0,3) = r30;
425 fMatrix(3,1) = fMatrix(1,3) = r31;
426 fMatrix(4,0) = fMatrix(0,4) = r40;
427 fMatrix(4,1) = fMatrix(1,4) = r41;
428
429 fStateR = rk;
430
431 if (rkm1 < fStateR) // going to greater R --> energy LOSS
432 return -1;
433 else // going to smaller R --> energy GAIN
434 return 1;
435}
436//
437//
438//
439Bool_t AliITSNeuralTrack::SeedCovariance()
440{
441 // generate a covariance matrix depending on the results obtained from
442 // the preliminary seeding fit procedure.
443 // It calculates the variances for C, D ans TanL, according to the
444 // differences of the fitted values from the requested ones necessary
445 // to make the curve exactly pass through each point.
446
447 /*
448 Int_t i, j;
449 AliITSNeuralPoint *p = 0;
450 Double_t r, argPhi, phiC, phiD, argZ, zL;
451 Double_t sumC = 0.0, sumD = 0.0, sumphi = 0., sumz = 0., sumL = 0.;
452 for (i = 0; i < fNum; i++) {
453 p = At(i);
454 if (!p) continue;
455 r = p->GetR2();
456 // weight and derivatives of phi and zeta w.r.t. various params
457 sumphi += 1./ p->ErrorGetPhi();
458 argPhi = ArgPhi(r);
459 argZ = ArgZ(r);
460 if (argPhi > 100.0 || argZ > 100.0) {
461 Error("InitCovariance", "Argument error");
462 return kFALSE;
463 }
464 phiC = DerArgPhiC(r) / TMath::Sqrt(1.0 - argPhi * argPhi);
465 phiD = DerArgPhiD(r) / TMath::Sqrt(1.0 - argPhi * argPhi);
466 if (phiC > 100.0 || phiD > 100.0) {
467 Error("InitCovariance", "Argument error");
468 return kFALSE;
469 }
470 zL = asin(argZ) / fC;
471 sumL += zL * zL;
472 sumC += phiC * phiC;
473 sumD += phiD * phiD;
474 sumz += 1.0 / (p->fError[2] * p->fError[2]);
475 }
476
477 for (i = 0; i < 5; i++) for (j = 0; j < 5; j++) fMatrix(i,j) = 0.;
478 fMatrix(0,0) = 1. / sumphi;
479 fMatrix(1,1) = 1. / sumz;
480 fMatrix(2,2) = 1. / sumD;
481 fMatrix(3,3) = 1. / sumL;
482 fMatrix(4,4) = 1. / sumC;
483 fMatrix.Print();
484 */
485
486 AliITSNeuralPoint *p = 0;
487 Double_t delta, cs, sn, r, argz;
488 Double_t diffC, diffD, diffL, calcC, calcD, calcL;
489
490 Int_t l;
491 for (l = 0; l < 6; l++) {
492 p = fPoint[l];
493 if (!p) break;
494 sn = TMath::Sin(p->GetPhi() - fG0);
495 cs = TMath::Cos(p->GetPhi() - fG0);
496 r = p->GetR2();
497 calcC = (fDt/r - sn) / (2.*fDt*sn - r - fDt*fDt/r);
498 argz = ArgZ(r);
499 if (argz > 1000.0) {
500 Error("Covariance", "Value too high");
501 return kFALSE;
502 }
503 calcL = (p->Z() - fDz) * fC / asin(argz);
504 delta = fR*fR + r*r + 2.0*fR*r*sin(p->GetPhi() - fG0);
505 if (delta < 0.E0) {
506 if (delta >= -0.5)
507 delta = 0.;
508 else {
509 Error("Covariance", Form("Discriminant = %g --- Dt = %g", delta, fDt));
510 return kFALSE;
511 }
512 }
513 delta = sqrt(delta);
514 if (fC >= 0)
515 calcD = delta - fR;
516 else
517 calcD = fR - delta;
518 diffD = calcD - fDt;
519 diffL = calcL - fTanL;
520 diffC = fC - calcC;
521 fMatrix(0,0) += 100000000.0 * p->GetError("phi") * p->GetError("phi");
522 fMatrix(1,1) += 10000.0 * p->ErrZ() * p->ErrZ();
523 fMatrix(2,2) += 100000.0 * diffD * diffD;
524 fMatrix(3,3) += diffL * diffL;
525 fMatrix(4,4) += 100000000.0 * diffC * diffC;
526 }
853a0f19 527 Double_t n = 0.;
528 for (l = 0; l < 6; l++) if (fPoint[l]) n++;
f77f13c8 529 fMatrix *= 1./(n * (n+1));
b9d722bc 530 return kTRUE;
531}
532//
533//
534//
535Bool_t AliITSNeuralTrack::Filter(AliITSNeuralPoint *test)
536{
537 // Makes all calculations which apply the Kalman filter to the
538 // stored guess of the state vector, after propagation to a new layer
539
540 if (!test) {
541 Error("Filter", "Null pointer passed");
542 return kFALSE;
543 }
544
545 Double_t m[2];
546 Double_t rk, phik, zk;
547 rk = test->GetR2();
548 phik = test->GetPhi();
549 zk = test->Z();
550 m[0]=phik;
551 m[1]=zk;
552
553 //////////////////////// Evaluation of the error matrix V /////////
554 Double_t v00 = test->GetError("phi") * rk;
555 Double_t v11 = test->ErrZ();
556 ////////////////////////////////////////////////////////////////////
557
558 // Get the covariance matrix
559 Double_t cin00, cin10, cin20, cin30, cin40;
560 Double_t cin11, cin21, cin31, cin41, cin22;
561 Double_t cin32, cin42, cin33, cin43, cin44;
562 cin00 = fMatrix(0,0);
563 cin10 = fMatrix(1,0);
564 cin20 = fMatrix(2,0);
565 cin30 = fMatrix(3,0);
566 cin40 = fMatrix(4,0);
567 cin11 = fMatrix(1,1);
568 cin21 = fMatrix(2,1);
569 cin31 = fMatrix(3,1);
570 cin41 = fMatrix(4,1);
571 cin22 = fMatrix(2,2);
572 cin32 = fMatrix(3,2);
573 cin42 = fMatrix(4,2);
574 cin33 = fMatrix(3,3);
575 cin43 = fMatrix(4,3);
576 cin44 = fMatrix(4,4);
577
578 // Calculate R matrix
579 Double_t rold00 = cin00 + v00;
580 Double_t rold10 = cin10;
581 Double_t rold11 = cin11 + v11;
582
583 ////////////////////// R matrix inversion /////////////////////////
584 Double_t det = rold00*rold11 - rold10*rold10;
585 Double_t r00 = rold11/det;
586 Double_t r10 = -rold10/det;
587 Double_t r11 = rold00/det;
588 ////////////////////////////////////////////////////////////////////
589
590 // Calculate Kalman matrix
591 Double_t k00 = cin00*r00 + cin10*r10;
592 Double_t k01 = cin00*r10 + cin10*r11;
593 Double_t k10 = cin10*r00 + cin11*r10;
594 Double_t k11 = cin10*r10 + cin11*r11;
595 Double_t k20 = cin20*r00 + cin21*r10;
596 Double_t k21 = cin20*r10 + cin21*r11;
597 Double_t k30 = cin30*r00 + cin31*r10;
598 Double_t k31 = cin30*r10 + cin31*r11;
599 Double_t k40 = cin40*r00 + cin41*r10;
600 Double_t k41 = cin40*r10 + cin41*r11;
601
602 // Get state vector (will keep the old values for phi and z)
603 Double_t x0, x1, x2, x3, x4, savex0, savex1;
604 x0 = savex0 = fStatePhi;
605 x1 = savex1 = fStateZ;
606 x2 = fDt;
607 x3 = fTanL;
608 x4 = fC;
609
610 // Update the state vector
611 x0 += k00*(m[0]-savex0) + k01*(m[1]-savex1);
612 x1 += k10*(m[0]-savex0) + k11*(m[1]-savex1);
613 x2 += k20*(m[0]-savex0) + k21*(m[1]-savex1);
614 x3 += k30*(m[0]-savex0) + k31*(m[1]-savex1);
615 x4 += k40*(m[0]-savex0) + k41*(m[1]-savex1);
616
617 // Update the covariance matrix
618 Double_t cout00, cout10, cout20, cout30, cout40;
619 Double_t cout11, cout21, cout31, cout41, cout22;
620 Double_t cout32, cout42, cout33, cout43, cout44;
621
622 cout00 = cin00 - k00*cin00 - k01*cin10;
623 cout10 = cin10 - k00*cin10 - k01*cin11;
624 cout20 = cin20 - k00*cin20 - k01*cin21;
625 cout30 = cin30 - k00*cin30 - k01*cin31;
626 cout40 = cin40 - k00*cin40 - k01*cin41;
627 cout11 = cin11 - k10*cin10 - k11*cin11;
628 cout21 = cin21 - k10*cin20 - k11*cin21;
629 cout31 = cin31 - k10*cin30 - k11*cin31;
630 cout41 = cin41 - k10*cin40 - k11*cin41;
631 cout22 = cin22 - k20*cin20 - k21*cin21;
632 cout32 = cin32 - k20*cin30 - k21*cin31;
633 cout42 = cin42 - k20*cin40 - k21*cin41;
634 cout33 = cin33 - k30*cin30 - k31*cin31;
635 cout43 = cin43 - k30*cin40 - k31*cin41;
636 cout44 = cin44 - k40*cin40 - k41*cin41;
637
638 // Store the new covariance matrix
639 fMatrix(0,0) = cout00;
640 fMatrix(1,0) = fMatrix(0,1) = cout10;
641 fMatrix(2,0) = fMatrix(0,2) = cout20;
642 fMatrix(3,0) = fMatrix(0,3) = cout30;
643 fMatrix(4,0) = fMatrix(0,4) = cout40;
644 fMatrix(1,1) = cout11;
645 fMatrix(2,1) = fMatrix(1,2) = cout21;
646 fMatrix(3,1) = fMatrix(1,3) = cout31;
647 fMatrix(4,1) = fMatrix(1,4) = cout41;
648 fMatrix(2,2) = cout22;
649 fMatrix(3,2) = fMatrix(2,3) = cout32;
650 fMatrix(4,2) = fMatrix(2,4) = cout42;
651 fMatrix(3,3) = cout33;
652 fMatrix(4,3) = fMatrix(3,4) = cout43;
653 fMatrix(4,4) = cout44;
654
655 // Calculation of the chi2 increment
656 Double_t vmcold00 = v00 - cout00;
657 Double_t vmcold10 = -cout10;
658 Double_t vmcold11 = v11 - cout11;
659 ////////////////////// Matrix vmc inversion ///////////////////////
660 det = vmcold00*vmcold11 - vmcold10*vmcold10;
661 Double_t vmc00=vmcold11/det;
662 Double_t vmc10 = -vmcold10/det;
663 Double_t vmc11 = vmcold00/det;
664 ////////////////////////////////////////////////////////////////////
665 Double_t chi2 = (m[0] - x0)*( vmc00*(m[0] - x0) + 2.*vmc10*(m[1] - x1) ) + (m[1] - x1)*vmc11*(m[1] - x1);
666 fChi2 += chi2;
667 fNSteps++;
668
669 return kTRUE;
670}
671//
672//
673//
674Bool_t AliITSNeuralTrack::KalmanFit()
853a0f19 675{
b9d722bc 676// Applies the Kalman Filter to improve the track parameters resolution.
677// First, thre point which lies closer to the estimated helix is chosen.
678// Then, a fit is performed towards the 6th layer
679// Finally, the track is refitted to the 1st layer
853a0f19 680
b9d722bc 681 Double_t rho;
682 Int_t l, layer, sign;
683
684 fStateR = fPoint[0]->GetR2();
685 fStatePhi = fPoint[0]->GetPhi();
686 fStateZ = fPoint[0]->Z();
687
688 if (!PropagateTo(3.0)) {
689 Error("KalmanFit", "Unsuccessful initialization");
690 return kFALSE;
691 }
692 l=0;
693
694 // Performs a Kalman filter going from the actual state position
695 // towards layer 6 position
696 // Now, the propagation + filtering operations can be performed
697 Double_t argPhi = 0.0, argZ = 0.0;
698 while (l <= 5) {
699 if (!fPoint[l]) {
700 Error("KalmanFit", "Not six points!");
701 return kFALSE;
702 }
703 layer = fPoint[l]->GetLayer();
704 rho = fPoint[l]->GetR2();
705 sign = PropagateTo(rho);
706 if (!sign) return kFALSE;
707 AddEL(layer, -1.0);
708 AddMS(layer);
709 if (!Filter(fPoint[l])) return kFALSE;
710 // these two parameters are update according to the filtered values
711 argPhi = ArgPhi(fStateR);
712 argZ = ArgZ(fStateR);
713 if (argPhi > 1.0 || argPhi < -1.0 || argZ > 1.0 || argZ < -1.0) {
714 Error("Filter", Form("Filtering returns too large values: %g, %g", argPhi, argZ));
715 return kFALSE;
716 }
717 fG0 = fStatePhi - asin(argPhi);
718 fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
719 l++;
720 }
721
722 // Now a Kalman filter i performed going from the actual state position
723 // towards layer 1 position and then propagates to vertex
724 if (l >= 5) l = 5;
725 while (l >= 1) {
726 layer = fPoint[l]->GetLayer();
727 rho = fPoint[l]->GetR2();
728 AddEL(layer, 1.0);
729 sign = PropagateTo(rho);
730 if (!sign) return kFALSE;
731 AddMS(layer);
732 if (!Filter(fPoint[l])) return kFALSE;
733 // these two parameters are update according to the filtered values
734 argPhi = ArgPhi(fStateR);
735 argZ = ArgZ(fStateR);
736 if (argPhi > 1.0 || argPhi < -1.0 || argZ > 1.0 || argZ < -1.0) {
737 Error("Filter", Form("Filtering returns too large values: %g, %g", argPhi, argZ));
738 return kFALSE;
739 }
740 fG0 = fStatePhi - asin(argPhi);
741 fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
742 l--;
743 }
744 return kTRUE;
745}
746//
747//
748//
749Bool_t AliITSNeuralTrack::RiemannFit()
750{
751 // Method which executes the circle fit via a Riemann Sphere projection
752 // with the only improvement of a weighted mean, due to different errors
753 // over different point measurements.
754 // As an output, it returns kTRUE or kFALSE respectively if the fit succeeded or not
755 // in fact, if some variables assume strange values, the fit is aborted,
756 // in order to prevent the class from raising a floating point error;
757
758 Int_t i, j;
759
760 // M1 - matrix of ones
761 TMatrixD m1(6,1);
762 for (i = 0; i < 6; i++) m1(i,0) = 1.0;
763
764 // X - matrix of Rieman projection coordinates
853a0f19 765 TMatrixD mX(6,3);
b9d722bc 766 for (i = 0; i < 6; i++) {
853a0f19 767 mX(i,0) = fPoint[i]->X();
768 mX(i,1) = fPoint[i]->Y();
769 mX(i,2) = fPoint[i]->GetR2sq();
b9d722bc 770 }
771
772 // W - matrix of weights
773 Double_t xterm, yterm, ex, ey;
853a0f19 774 TMatrixD mW(6,6);
b9d722bc 775 for (i = 0; i < 6; i++) {
776 xterm = fPoint[i]->X() * fPoint[i]->GetPhi() - fPoint[i]->Y() / fPoint[i]->GetR2();
777 ex = fPoint[i]->ErrX();
778 yterm = fPoint[i]->Y() * fPoint[i]->GetPhi() + fPoint[i]->X() / fPoint[i]->GetR2();
779 ey = fPoint[i]->ErrY();
853a0f19 780 mW(i,i) = fPoint[i]->GetR2sq() / (xterm * xterm * ex * ex + yterm * yterm * ey * ey);
b9d722bc 781 }
782
783 // Xm - weighted sample mean
853a0f19 784 Double_t meanX = 0.0, meanY = 0.0, meanW = 0.0, sw = 0.0;
b9d722bc 785 for (i = 0; i < 6; i++) {
853a0f19 786 meanX += mW(i,i) * mX(i,0);
787 meanY += mW(i,i) * mX(i,1);
788 meanW += mW(i,i) * mX(i,2);
789 sw += mW(i,i);
b9d722bc 790 }
853a0f19 791 meanX /= sw;
792 meanY /= sw;
793 meanW /= sw;
b9d722bc 794
795 // V - sample covariance matrix
796 for (i = 0; i < 6; i++) {
853a0f19 797 mX(i,0) -= meanX;
798 mX(i,1) -= meanY;
799 mX(i,2) -= meanW;
b9d722bc 800 }
853a0f19 801 TMatrixD mXt(TMatrixD::kTransposed, mX);
802 TMatrixD mWX(mW, TMatrixD::kMult, mX);
803 TMatrixD mV(mXt, TMatrixD::kMult, mWX);
b9d722bc 804 for (i = 0; i < 3; i++) {
805 for (j = i + 1; j < 3; j++) {
853a0f19 806 mV(i,j) = mV(j,i) = (mV(i,j) + mV(j,i)) * 0.5;
b9d722bc 807 }
808 }
809
810 // Eigenvalue problem solving for V matrix
811 Int_t ileast = 0;
853a0f19 812 TVectorD eval(3), n(3);
49370114 813 // TMatrixD evec = mV.EigenVectors(eval);
814#if ROOT_VERSION_CODE >= 262146
815 TMatrixDEigen ei(mV);
816 TMatrixD evec = ei.GetEigenVectors();
817 eval = ei.GetEigenValues();
818#else
853a0f19 819 TMatrixD evec = mV.EigenVectors(eval);
49370114 820#endif
821
853a0f19 822 if (eval(1) < eval(ileast)) ileast = 1;
823 if (eval(2) < eval(ileast)) ileast = 2;
824 n(0) = evec(0, ileast);
825 n(1) = evec(1, ileast);
826 n(2) = evec(2, ileast);
b9d722bc 827
828 // c - known term in the plane intersection with Riemann axes
853a0f19 829 Double_t c = -(meanX * n(0) + meanY * n(1) + meanW * n(2));
b9d722bc 830
831 fXC = -n(0) / (2. * n(2));
832 fYC = -n(1) / (2. * n(2));
833 fR = (1. - n(2)*n(2) - 4.*c*n(2)) / (4. * n(2) * n(2));
834
835 if (fR <= 0.E0) {
836 Error("RiemannFit", "Radius comed less than zero!!!");
837 return kFALSE;
838 }
839 fR = TMath::Sqrt(fR);
840 fC = 1.0 / fR;
841
842 // evaluating signs for curvature and others
843 Double_t phi1 = 0.0, phi2, temp1, temp2, sumdphi = 0.0, ref = TMath::Pi();
844 AliITSNeuralPoint *p = fPoint[0];
845 phi1 = p->GetPhi();
846 for (i = 1; i < 6; i++) {
847 p = (AliITSNeuralPoint*)fPoint[i];
848 if (!p) break;
849 phi2 = p->GetPhi();
850 temp1 = phi1;
851 temp2 = phi2;
852 if (temp1 > ref && temp2 < ref)
853 temp2 += 2.0 * ref;
854 else if (temp1 < ref && temp2 > ref)
855 temp1 += 2.0 * ref;
856 sumdphi += temp2 - temp1;
857 phi1 = phi2;
858 }
859 if (sumdphi < 0.E0) fC = -fC;
860 Double_t diff, angle = TMath::ATan2(fYC, fXC);
861 if (fC < 0.E0)
862 fG0 = angle + 0.5 * TMath::Pi();
863 else
864 fG0 = angle - 0.5 * TMath::Pi();
865 diff = angle - fG0;
866
853a0f19 867 Double_t d = TMath::Sqrt(fXC*fXC + fYC*fYC) - fR;
b9d722bc 868 if (fC >= 0.E0)
853a0f19 869 fDt = d;
b9d722bc 870 else
853a0f19 871 fDt = -d;
b9d722bc 872
853a0f19 873 Int_t nn = 6;
b9d722bc 874 Double_t halfC = 0.5 * fC;
853a0f19 875 Double_t *s = new Double_t[nn], *z = new Double_t[nn], *ws = new Double_t[nn];
b9d722bc 876 for (j = 0; j < 6; j++) {
877 p = fPoint[j];
878 if (!p) break;
879 s[j] = ArgZ(p->GetR2());
880 if (s[j] > 100.0) return kFALSE;
881 z[j] = p->Z();
882 s[j] = asin(s[j]) / halfC;
883 ws[j] = 1.0 / (p->ErrZ()*p->ErrZ());
884 }
885
886 // second tep final fit
853a0f19 887 Double_t sums2 = 0.0, sumz = 0.0, sumsz = 0.0, sums = 0.0, sumw = 0.0;
888 for (i = 0; i < nn; i++) {
889 sums2 += ws[i] * s[i] * s[i];
890 sumz += ws[i] * z[i];
891 sums += ws[i] * s[i];
892 sumsz += ws[i] * s[i] * z[i];
b9d722bc 893 sumw += ws[i];
894 }
853a0f19 895 sums2 /= sumw;
896 sumz /= sumw;
897 sums /= sumw;
898 sumsz /= sumw;
899 d = sums2 - sums*sums;
b9d722bc 900
853a0f19 901 fDz = (sums2*sumz - sums*sumsz) / d;
902 fTanL = (sumsz - sums*sumz) / d;
b9d722bc 903
904 delete [] s;
905 delete [] z;
906 delete [] ws;
907
908 return kTRUE;
909}
910//
911//
912//
913void AliITSNeuralTrack::PrintState(Bool_t matrix)
853a0f19 914{
b9d722bc 915// Prints the state vector values.
916// The argument switches on or off the printing of the covariance matrix.
853a0f19 917
b9d722bc 918 cout << "\nState vector: " << endl;
919 cout << " Rho = " << fStateR << "\n";
920 cout << " Phi = " << fStatePhi << "\n";
921 cout << " Z = " << fStateZ << "\n";
922 cout << " Dt = " << fDt << "\n";
923 cout << " Dz = " << fDz << "\n";
924 cout << "TanL = " << fTanL << "\n";
925 cout << " C = " << fC << "\n";
926 cout << " G0 = " << fG0 << "\n";
927 cout << " XC = " << fXC << "\n";
928 cout << " YC = " << fYC << "\n";
929 if (matrix) {
930 cout << "\nCovariance Matrix: " << endl;
931 fMatrix.Print();
932 }
933 cout << "Actual square chi = " << fChi2;
934}
935//
936//
937//
853a0f19 938Double_t AliITSNeuralTrack::GetDz() const
b9d722bc 939{
940// Double_t argZ = ArgZ(fStateR);
941// if (argZ > 9.9) {
942// Error("GetDz", "Too large value: %g", argZ);
943// return 0.0;
944// }
945// fDz = fStateZ - (2.0 * fTanL / fC) * asin(argZ);
946 return fDz;
947}
948//
949//
950//
853a0f19 951Double_t AliITSNeuralTrack::GetGamma() const
b9d722bc 952{
953// these two parameters are update according to the filtered values
954// Double_t argPhi = ArgPhi(fStateR);
955// if (argPhi > 9.9) {
956// Error("Filter", "Too large value: %g", argPhi);
957// return kFALSE;
958// }
959// fG0 = fStatePhi - asin(argPhi);
960 return fG0;
961}
962//
963//
964//
853a0f19 965Double_t AliITSNeuralTrack::GetPhi(Double_t r) const
966{
b9d722bc 967// Gives the value of azymuthal coordinate in the helix
968// as a function of cylindric radius
853a0f19 969
b9d722bc 970 Double_t arg = ArgPhi(r);
971 if (arg > 0.9) return 0.0;
972 arg = fG0 + asin(arg);
973 while (arg >= 2.0 * TMath::Pi()) { arg -= 2.0 * TMath::Pi(); }
974 while (arg < 0.0) { arg += 2.0 * TMath::Pi(); }
975 return arg;
976}
977//
978//
979//
853a0f19 980Double_t AliITSNeuralTrack::GetZ(Double_t r) const
981{
b9d722bc 982// gives the value of Z in the helix
983// as a function of cylindric radius
853a0f19 984
b9d722bc 985 Double_t arg = ArgZ(r);
986 if (arg > 0.9) return 0.0;
987 return fDz + fTanL * asin(arg) / fC;
988}
989//
990//
991//
992Double_t AliITSNeuralTrack::GetdEdX()
993{
853a0f19 994// total energy loss of the track
995
b9d722bc 996 Double_t q[4] = {0., 0., 0., 0.}, dedx = 0.0;
997 Int_t i = 0, swap = 0;
998 for (i = 2; i < 6; i++) {
999 if (!fPoint[i]) continue;
1000 q[i - 2] = (Double_t)fPoint[i]->GetCharge();
1001 q[i - 2] /= (1 + fTanL*fTanL);
1002 }
1003 q[0] /= 280.;
1004 q[1] /= 280.;
1005 q[2] /= 38.;
1006 q[3] /= 38.;
1007 do {
1008 swap = 0;
1009 for (i = 0; i < 3; i++) {
1010 if (q[i] <= q[i + 1]) continue;
1011 Double_t tmp = q[i];
1012 q[i] = q[i + 1];
1013 q[i+1] = tmp;
1014 swap++;
1015 }
1016 } while(swap);
1017 if(q[0] < 0.) {
1018 q[0] = q[1];
1019 q[1] = q[2];
1020 q[2] = q[3];
1021 q[3] = -1.;
1022 }
1023 dedx = (q[0] + q[1]) / 2.;
1024 return dedx;
1025}
1026//
1027//
1028//
1029void AliITSNeuralTrack::SetVertex(Double_t *pos, Double_t *err)
1030{
1031 // Stores vertex data
1032
1033 if (!pos || !err) return;
1034 fVertex.ErrX() = err[0];
1035 fVertex.ErrY() = err[1];
1036 fVertex.ErrZ() = err[2];
1037 fVertex.SetLayer(0);
1038 fVertex.SetModule(0);
1039 fVertex.SetIndex(0);
1040 fVertex.SetLabel(0, -1);
1041 fVertex.SetLabel(1, -1);
1042 fVertex.SetLabel(2, -1);
1043 fVertex.SetUser(1);
1044}
1045//
1046//
1047//
1048AliITSIOTrack* AliITSNeuralTrack::ExportIOtrack(Int_t min)
b9d722bc 1049{
853a0f19 1050// Exports an object in the standard format for reconstructed tracks
1051
b9d722bc 1052 Int_t layer = 0;
1053 AliITSIOTrack *track = new AliITSIOTrack;
1054
1055 // covariance matrix
1056 track->SetCovMatrix(fMatrix(0,0), fMatrix(1,0), fMatrix(1,1),
1057 fMatrix(2,0), fMatrix(2,1), fMatrix(2,2),
1058 fMatrix(3,0), fMatrix(3,1), fMatrix(3,2),
1059 fMatrix(3,3), fMatrix(4,0), fMatrix(4,1),
1060 fMatrix(4,2), fMatrix(4,3), fMatrix(4,4));
1061
1062 // labels
1063 track->SetLabel(IsGood(min) ? fLabel : -fLabel);
1064 track->SetTPCLabel(-1);
1065
1066 // points characteristics
1067 for (layer = 0; layer < 6; layer++) {
1068 if (fPoint[layer]) {
1069 track->SetIdModule(layer, fPoint[layer]->GetModule());
1070 track->SetIdPoint(layer, fPoint[layer]->GetIndex());
1071 }
1072 }
1073
1074 // state vector
1075 track->SetStatePhi(fStatePhi);
1076 track->SetStateZ(fStateZ);
1077 track->SetStateD(fDt);
1078 track->SetStateTgl(fTanL);
1079 track->SetStateC(fC);
1080 track->SetRadius(fStateR);
1081 track->SetCharge((fC > 0.0) ? -1 : 1);
1082 track->SetDz(fDz);
1083
1084 // track parameters in the closest point
1085 track->SetX(fStateR * cos(fStatePhi));
1086 track->SetY(fStateR * cos(fStatePhi));
1087 track->SetZ(fStateZ);
1088 track->SetPx(GetPt() * cos(fG0));
1089 track->SetPy(GetPt() * sin(fG0));
1090 track->SetPz(GetPt() * fTanL);
1091
1092 // PID
1093 track->SetPid(fPDG);
1094 track->SetMass(fMass);
1095
1096 return track;
1097}
1098//
1099//
1100//====================================================================================
1101//============================ PRIVATE METHODS ============================
1102//====================================================================================
1103//
1104//
1105Double_t AliITSNeuralTrack::ArgPhi(Double_t r) const
1106{
1107 // calculates the expression ((1/2)Cr + (1 + (1/2)CD) D/r) / (1 + CD)
1108
1109 Double_t arg, num, den;
1110 num = (0.5 * fC * r) + (1. + (0.5 * fC * fDt)) * (fDt / r);
1111 den = 1. + fC * fDt;
1112 if (den == 0.) {
1113 Error("ArgPhi", "Denominator = 0!");
1114 return 10.0;
1115 }
1116 arg = num / den;
1117 if (TMath::Abs(arg) < 1.) return arg;
1118 if (TMath::Abs(arg) <= 1.00001) return (arg > 0.) ? 0.99999999999 : -0.9999999999;
1119 Error("ArgPhi", "Value too large: %17.15g", arg);
1120 return 10.0;
1121}
1122//
1123//
1124//
1125Double_t AliITSNeuralTrack::ArgZ(Double_t r) const
1126{
1127 // calculates the expression (1/2)C * sqrt( (r^2 - Dt^2) / (1 + CD) )
1128
1129 Double_t arg;
1130 arg = (r * r - fDt * fDt) / (1. + fC * fDt);
1131 if (arg < 0.) {
b9f05b32 1132 if (TMath::Abs(arg) < 1.E-6) arg = 0.;
b9d722bc 1133 else {
1134 Error("ArgZ", "Square root argument error: %17.15g < 0", arg);
1135 return 10.;
1136 }
1137 }
1138 arg = 0.5 * fC * TMath::Sqrt(arg);
1139 if (TMath::Abs(arg) < 1.) return arg;
1140 if (TMath::Abs(arg) <= 1.00001) return (arg > 0.) ? 0.99999999999 : -0.9999999999;
1141 Error("ArgZ", "Value too large: %17.15g", arg);
1142 return 10.0;
1143}
1144//
1145//
1146//
1147Double_t AliITSNeuralTrack::ArgB(Double_t r) const
1148{
853a0f19 1149// UTILITY FUNCTION
1150
b9d722bc 1151 Double_t arg;
1152 arg = (r*r - fDt*fDt);
1153 arg /= (r*(1.+ fC*fDt)*(1.+ fC*fDt));
1154 return arg;
1155}
1156//
1157//
1158//
1159Double_t AliITSNeuralTrack::ArgC(Double_t r) const
1160{
853a0f19 1161// UTILITY FUNCTION
1162
b9d722bc 1163 Double_t arg;
1164 arg = (1./r - fC * ArgPhi(r));
1165 arg /= 1.+ fC*fDt;
1166 return arg;
1167}