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