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