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