Print removed
[u/mrichter/AliRoot.git] / ITS / AliITSNeuralTrack.cxx
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
9 #include <Riostream.h>
10 //#include <cstdlib>
11 //#include <cstring>
12
13 //#include <TObject.h>
14 //#include <TROOT.h>
15 #include <TMath.h>
16 #include <TString.h>
17 //#include <TObjArray.h>
18 //#include <TH1.h>
19 #include <TMatrixD.h>
20 #if ROOT_VERSION_CODE >= 262146
21 #include <TMatrixDEigen.h>
22 #endif
23
24 //#include "AliITSVertex.h"
25 #include "AliITSIOTrack.h"
26 #include "AliITSNeuralPoint.h"
27
28 #include "AliITSNeuralTrack.h"
29
30
31
32 ClassImp(AliITSNeuralTrack)
33 //
34 //
35 //
36 AliITSNeuralTrack::AliITSNeuralTrack() : fMatrix(5,5), fVertex()
37 {
38         // Default constructor
39         
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 //
63 AliITSNeuralTrack::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 //
91 AliITSNeuralTrack::~AliITSNeuralTrack()
92 {
93         Int_t l;
94         for (l = 0; l < 6; l++) fPoint[l] = 0;
95 }
96 //
97 //
98 //
99 void AliITSNeuralTrack::AssignLabel()
100 {       
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. 
107
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 //
157 void AliITSNeuralTrack::CleanSlot(Int_t i, Bool_t del)
158 {
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.
161
162         if (i >= 0 && i < 6) {
163                 if (del) delete fPoint[i];
164                 fPoint[i] = 0;
165         }
166 }
167 //
168 //
169 //
170 void AliITSNeuralTrack::GetModuleData(Int_t layer, Int_t &mod, Int_t &pos)
171 {
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.
175
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 //
186 void AliITSNeuralTrack::Insert(AliITSNeuralPoint *point)
187 {
188 // A trivial method to insert a point in the tracks;
189 // the point is inserted to the slot corresponding to its ITS layer.
190
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 //
204 Int_t AliITSNeuralTrack::OccupationMask()
205 {
206 // Returns a byte which maps the occupied slots. 
207 // Each bit represents a layer going from the less significant on.
208
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 //
219 void AliITSNeuralTrack::PrintLabels()
220 {
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.
224
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 //
242 Bool_t AliITSNeuralTrack::AddEL(Int_t layer, Double_t sign)
243 {
244 // Calculates the correction for energy loss
245
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 //
289 Bool_t AliITSNeuralTrack::AddMS(Int_t layer)
290 {
291 // Calculates the noise perturbation due to multiple scattering
292
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 //
341 Int_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
353         if (rk < TMath::Abs(fDt)) {
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 //
422 Bool_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         }
510         Double_t n = 0.;
511         for (l = 0; l < 6; l++) if (fPoint[l]) n++;
512         fMatrix *= 1./(n++ * n);
513         return kTRUE;
514 }
515 //
516 //
517 //
518 Bool_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 //
657 Bool_t AliITSNeuralTrack::KalmanFit()
658 {
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
663
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 //
732 Bool_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
748         TMatrixD mX(6,3);
749         for (i = 0; i < 6; i++) {
750                 mX(i,0) = fPoint[i]->X();
751                 mX(i,1) = fPoint[i]->Y();
752                 mX(i,2) = fPoint[i]->GetR2sq();
753         }
754
755         // W - matrix of weights
756         Double_t xterm, yterm, ex, ey;
757         TMatrixD mW(6,6);
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();
763                 mW(i,i) = fPoint[i]->GetR2sq() / (xterm * xterm * ex * ex + yterm * yterm * ey * ey);
764         }
765
766         // Xm - weighted sample mean
767         Double_t meanX = 0.0, meanY = 0.0, meanW = 0.0, sw = 0.0;
768         for (i = 0; i < 6; i++) {
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);
773         }
774         meanX /= sw;
775         meanY /= sw;
776         meanW /= sw;
777
778         // V - sample covariance matrix
779         for (i = 0; i < 6; i++) {
780                 mX(i,0) -= meanX;
781                 mX(i,1) -= meanY;
782                 mX(i,2) -= meanW;
783         }
784         TMatrixD mXt(TMatrixD::kTransposed, mX);
785         TMatrixD mWX(mW, TMatrixD::kMult, mX);
786         TMatrixD mV(mXt, TMatrixD::kMult, mWX);
787         for (i = 0; i < 3; i++) {
788                 for (j = i + 1; j < 3; j++) {
789                         mV(i,j)  = mV(j,i)  = (mV(i,j) + mV(j,i)) * 0.5;
790                 }
791         }
792
793         // Eigenvalue problem solving for V matrix
794         Int_t ileast = 0;
795         TVectorD eval(3), n(3);
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
802         TMatrixD evec = mV.EigenVectors(eval);
803 #endif
804
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);
810
811         // c - known term in the plane intersection with Riemann axes
812         Double_t c = -(meanX * n(0) + meanY * n(1) + meanW * n(2));
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
850         Double_t d = TMath::Sqrt(fXC*fXC + fYC*fYC) - fR;
851         if (fC >= 0.E0)
852                 fDt = d;
853         else
854                 fDt = -d;
855
856         Int_t nn = 6;
857         Double_t halfC = 0.5 * fC;
858         Double_t *s = new Double_t[nn], *z = new Double_t[nn], *ws = new Double_t[nn];
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
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];
876                 sumw += ws[i];
877         }
878         sums2 /= sumw;
879         sumz /= sumw;
880         sums /= sumw;
881         sumsz /= sumw;
882         d = sums2 - sums*sums;
883
884         fDz = (sums2*sumz - sums*sumsz) / d;
885         fTanL = (sumsz - sums*sumz) / d;
886
887         delete [] s;
888         delete [] z;
889         delete [] ws;
890
891         return kTRUE;
892 }
893 //
894 //
895 //
896 void AliITSNeuralTrack::PrintState(Bool_t matrix)
897 {
898 // Prints the state vector values.
899 // The argument switches on or off the printing of the covariance matrix.
900
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 //
921 Double_t AliITSNeuralTrack::GetDz() const
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 //
934 Double_t AliITSNeuralTrack::GetGamma() const
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 //
948 Double_t AliITSNeuralTrack::GetPhi(Double_t r) const
949 {
950 // Gives the value of azymuthal coordinate in the helix
951 // as a function of cylindric radius
952
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 //
963 Double_t AliITSNeuralTrack::GetZ(Double_t r) const
964 {
965 // gives the value of Z in the helix
966 // as a function of cylindric radius
967
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 //
975 Double_t AliITSNeuralTrack::GetdEdX()
976 {
977 // total energy loss of the track
978
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 //
1012 void 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 //
1031 AliITSIOTrack* AliITSNeuralTrack::ExportIOtrack(Int_t min)
1032 {
1033 // Exports an object in the standard format for reconstructed tracks
1034
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 //
1088 Double_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 //
1108 Double_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.) {
1115                 if (TMath::Abs(arg) < 1.E-6) arg = 0.;
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 //
1130 Double_t AliITSNeuralTrack::ArgB(Double_t r) const 
1131 {
1132 // UTILITY FUNCTION 
1133
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 //
1142 Double_t AliITSNeuralTrack::ArgC(Double_t r) const 
1143 {
1144 // UTILITY FUNCTION
1145
1146         Double_t arg;
1147         arg = (1./r - fC * ArgPhi(r));
1148         arg /= 1.+ fC*fDt;
1149         return arg;
1150 }