Track time measurement (S.Radomski)
[u/mrichter/AliRoot.git] / ITS / AliITStrackV2.cxx
CommitLineData
006b5f7f 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//-------------------------------------------------------------------------
17// Implementation of the ITS track class
18//
19// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
a9a2d814 20// dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
006b5f7f 21//-------------------------------------------------------------------------
22
23#include <TMatrixD.h>
24
25#include <TMath.h>
4ae5bbc4 26#include <Riostream.h>
006b5f7f 27
28#include "AliCluster.h"
517b130f 29#include "AliTPCtrack.h"
006b5f7f 30#include "AliITStrackV2.h"
31
32ClassImp(AliITStrackV2)
33
7f6ddf58 34const Int_t kWARN=5;
006b5f7f 35
36//____________________________________________________________________________
e6ba4672 37AliITStrackV2::AliITStrackV2():AliKalmanTrack(),
22b13da0 38 fX(0),
39 fAlpha(0),
40 fdEdx(0),
41 fP0(0),
42 fP1(0),
43 fP2(0),
44 fP3(0),
45 fP4(0),
46 fC00(0),
47 fC10(0),
48 fC11(0),
49 fC20(0),
50 fC21(0),
51 fC22(0),
52 fC30(0),
53 fC31(0),
54 fC32(0),
55 fC33(0),
56 fC40(0),
57 fC41(0),
58 fC42(0),
59 fC43(0),
e6ba4672 60 fC44(0)
61 {
22b13da0 62 for(Int_t i=0; i<kMaxLayer; i++) fIndex[i]=0;
63 for(Int_t i=0; i<4; i++) fdEdxSample[i]=0;
64}
65//____________________________________________________________________________
e1a65e92 66AliITStrackV2::AliITStrackV2(const AliTPCtrack& t) throw (const Char_t *) :
67AliKalmanTrack(t) {
006b5f7f 68 //------------------------------------------------------------------
9b280d80 69 //Conversion TPC track -> ITS track
006b5f7f 70 //------------------------------------------------------------------
006b5f7f 71 SetChi2(0.);
72 SetNumberOfClusters(0);
14825d5a 73
7f6ddf58 74 fdEdx = t.GetdEdx();
75 SetMass(t.GetMass());
76
006b5f7f 77 fAlpha = t.GetAlpha();
78 if (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
79 else if (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
80
9b280d80 81 //Conversion of the track parameters
006b5f7f 82 Double_t x,p[5]; t.GetExternalParameters(x,p);
9b280d80 83 fX=x; x=GetConvConst();
006b5f7f 84 fP0=p[0];
85 fP1=p[1];
86 fP2=p[2];
87 fP3=p[3];
88 fP4=p[4]/x;
89
9b280d80 90 //Conversion of the covariance matrix
006b5f7f 91 Double_t c[15]; t.GetExternalCovariance(c);
92
93 fC00=c[0 ];
94 fC10=c[1 ]; fC11=c[2 ];
95 fC20=c[3 ]; fC21=c[4 ]; fC22=c[5 ];
96 fC30=c[6 ]; fC31=c[7 ]; fC32=c[8 ]; fC33=c[9 ];
97 fC40=c[10]/x; fC41=c[11]/x; fC42=c[12]/x; fC43=c[13]/x; fC44=c[14]/x/x;
98
99 if (!Invariant()) throw "AliITStrackV2: conversion failed !\n";
649ab0de 100
006b5f7f 101}
102
103//____________________________________________________________________________
104AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) : AliKalmanTrack(t) {
105 //------------------------------------------------------------------
106 //Copy constructor
107 //------------------------------------------------------------------
108 fX=t.fX;
109 fAlpha=t.fAlpha;
110 fdEdx=t.fdEdx;
111
112 fP0=t.fP0; fP1=t.fP1; fP2=t.fP2; fP3=t.fP3; fP4=t.fP4;
113
114 fC00=t.fC00;
115 fC10=t.fC10; fC11=t.fC11;
116 fC20=t.fC20; fC21=t.fC21; fC22=t.fC22;
117 fC30=t.fC30; fC31=t.fC31; fC32=t.fC32; fC33=t.fC33;
118 fC40=t.fC40; fC41=t.fC41; fC42=t.fC42; fC43=t.fC43; fC44=t.fC44;
119
120 Int_t n=GetNumberOfClusters();
23efe5f1 121 for (Int_t i=0; i<n; i++) {
a9a2d814 122 fIndex[i]=t.fIndex[i];
123 if (i<4) fdEdxSample[i]=t.fdEdxSample[i];
23efe5f1 124 }
006b5f7f 125}
a9a2d814 126
006b5f7f 127//_____________________________________________________________________________
128Int_t AliITStrackV2::Compare(const TObject *o) const {
129 //-----------------------------------------------------------------
130 // This function compares tracks according to the their curvature
131 //-----------------------------------------------------------------
7f6ddf58 132 AliITStrackV2 *t=(AliITStrackV2*)o;
a9a2d814 133 //Double_t co=TMath::Abs(t->Get1Pt());
134 //Double_t c =TMath::Abs(Get1Pt());
135 Double_t co=t->GetSigmaY2()*t->GetSigmaZ2();
136 Double_t c =GetSigmaY2()*GetSigmaZ2();
006b5f7f 137 if (c>co) return 1;
138 else if (c<co) return -1;
139 return 0;
140}
141
142//_____________________________________________________________________________
143void AliITStrackV2::GetExternalCovariance(Double_t cc[15]) const {
144 //-------------------------------------------------------------------------
145 // This function returns an external representation of the covriance matrix.
146 // (See comments in AliTPCtrack.h about external track representation)
147 //-------------------------------------------------------------------------
9b280d80 148 Double_t a=GetConvConst();
006b5f7f 149
150 cc[0 ]=fC00;
151 cc[1 ]=fC10; cc[2 ]=fC11;
152 cc[3 ]=fC20; cc[4 ]=fC21; cc[5 ]=fC22;
153 cc[6 ]=fC30; cc[7 ]=fC31; cc[8 ]=fC32; cc[9 ]=fC33;
154 cc[10]=fC40*a; cc[11]=fC41*a; cc[12]=fC42*a; cc[13]=fC43*a; cc[14]=fC44*a*a;
155}
156
157//____________________________________________________________________________
a9a2d814 158Int_t AliITStrackV2::PropagateToVertex(Double_t d,Double_t x0) {
006b5f7f 159 //------------------------------------------------------------------
160 //This function propagates a track to the minimal distance from the origin
161 //------------------------------------------------------------------
162 Double_t xv=fP2*(fX*fP2 - fP0*TMath::Sqrt(1.- fP2*fP2)); //linear approxim.
a9a2d814 163 PropagateTo(xv,d,x0);
006b5f7f 164 return 0;
165}
166
167//____________________________________________________________________________
168Int_t AliITStrackV2::
169GetGlobalXYZat(Double_t xk, Double_t &x, Double_t &y, Double_t &z) const {
170 //------------------------------------------------------------------
171 //This function returns a track position in the global system
172 //------------------------------------------------------------------
173 Double_t dx=xk-fX;
174 Double_t f1=fP2, f2=f1 + fP4*dx;
a9a2d814 175 if (TMath::Abs(f2) >= 0.9999) {
7f6ddf58 176 Int_t n=GetNumberOfClusters();
177 if (n>kWARN)
178 cerr<<n<<" AliITStrackV2::GetGlobalXYZat: Propagation failed !\n";
006b5f7f 179 return 0;
180 }
181
182 Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
183
184 Double_t yk = fP0 + dx*(f1+f2)/(r1+r2);
185 Double_t zk = fP1 + dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
186
187 Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
188 x = xk*cs - yk*sn;
189 y = xk*sn + yk*cs;
190 z = zk;
191
192 return 1;
193}
194
195//_____________________________________________________________________________
196Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c) const
197{
198 //-----------------------------------------------------------------
199 // This function calculates a predicted chi2 increment.
200 //-----------------------------------------------------------------
201 Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
202 r00+=fC00; r01+=fC10; r11+=fC11;
203
204 Double_t det=r00*r11 - r01*r01;
7f6ddf58 205 if (TMath::Abs(det) < 1.e-30) {
006b5f7f 206 Int_t n=GetNumberOfClusters();
7f6ddf58 207 if (n>kWARN)
208 cerr<<n<<" AliKalmanTrack::GetPredictedChi2: Singular matrix !\n";
006b5f7f 209 return 1e10;
210 }
211 Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
33329428 212
006b5f7f 213 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
33329428 214
006b5f7f 215 return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
216}
217
218//_____________________________________________________________________________
219Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c,Double_t *m,
7f6ddf58 220Double_t x0) const {
006b5f7f 221 //-----------------------------------------------------------------
222 // This function calculates a chi2 increment with a vertex contraint
223 //-----------------------------------------------------------------
224 TVectorD x(5); x(0)=fP0; x(1)=fP1; x(2)=fP2; x(3)=fP3; x(4)=fP4;
225 TMatrixD C(5,5);
226 C(0,0)=fC00;
227 C(1,0)=fC10; C(1,1)=fC11;
228 C(2,0)=fC20; C(2,1)=fC21; C(2,2)=fC22;
229 C(3,0)=fC30; C(3,1)=fC31; C(3,2)=fC32; C(3,3)=fC33;
230 C(4,0)=fC40; C(4,1)=fC41; C(4,2)=fC42; C(4,3)=fC43; C(4,4)=fC44;
231
232 C(0,1)=C(1,0);
233 C(0,2)=C(2,0); C(1,2)=C(2,1);
234 C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
235 C(0,4)=C(4,0); C(1,4)=C(4,1); C(2,4)=C(4,2); C(3,4)=C(4,3);
236
237 TMatrixD H(4,5); H.UnitMatrix();
238 Double_t dy=(c->GetY() - m[0]), dz=(c->GetZ() - m[1]);
239
240 Double_t dr=TMath::Sqrt(fX*fX + dy*dy);
241 Double_t r =TMath::Sqrt(4/dr/dr - fP4*fP4);
242 Double_t sn=0.5*(fP4*fX + dy*r);
243 Double_t tg=0.5*fP4*dz/TMath::ASin(0.5*fP4*dr);
244 TVectorD mm(4);
245 mm(0)=m[0]=c->GetY(); mm(1)=m[1]=c->GetZ(); mm(2)=m[2]=sn; mm(3)=m[3]=tg;
246
247 Double_t v22=0.,v33=0.;
248 //x0=0.;
249 if (x0!=0.) {
250 Double_t pp2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
7f6ddf58 251 Double_t beta2=pp2/(pp2 + GetMass()*GetMass());
006b5f7f 252 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
253 Double_t theta2=14.1*14.1/(beta2*pp2*1e6)*x0;
254 v22 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
255 v33 = theta2*(1.+ GetTgl()*GetTgl())*(1. + GetTgl()*GetTgl());
256 }
257 Double_t sy2=c->GetSigmaY2(), sz2=c->GetSigmaZ2();
258 v22+=kSigmaYV*kSigmaYV/dr/dr;
259 v22+=sy2/dr/dr;
260 Double_t v20=sy2/dr;
261
262 v33+=kSigmaZV*kSigmaZV/dr/dr;
263 v33+=sz2/dr/dr;
264 Double_t v31=sz2/dr;
265
266 TMatrixD V(4,4);
267 V(0,0)=m[4 ]=sy2; V(0,1)=m[5 ]=0.; V(0,2)=m[6 ]=v20; V(0,3)=m[7 ]=0.;
268 V(1,0)=m[8 ]=0.; V(1,1)=m[9 ]=sz2; V(1,2)=m[10]=0.; V(1,3)=m[11]=v31;
269 V(2,0)=m[12]=v20; V(2,1)=m[13]=0.; V(2,2)=m[14]=v22; V(2,3)=m[15]=0.;
270 V(3,0)=m[16]=0.; V(3,1)=m[17]=v31; V(3,2)=m[18]=0.; V(3,3)=m[19]=v33;
271
272 TVectorD res=x; res*=H; res-=mm; //res*=-1;
273 TMatrixD tmp(H,TMatrixD::kMult,C);
274 TMatrixD R(tmp,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed,H)); R+=V;
275
276 Double_t det=R.Determinant();
7f6ddf58 277 if (TMath::Abs(det) < 1.e-30) {
006b5f7f 278 Int_t n=GetNumberOfClusters();
7f6ddf58 279 if (n>kWARN)
280 cerr<<n<<" AliITStrackV2::GetPredictedChi2: Singular matrix !\n";
006b5f7f 281 return 1e10;
282 }
283
284 R.Invert();
285
286 TVectorD rs=res;
287 res*=R;
288 return rs*res;
289}
290
291//____________________________________________________________________________
a9a2d814 292Int_t AliITStrackV2::CorrectForMaterial(Double_t d, Double_t x0) {
293 //------------------------------------------------------------------
294 //This function corrects the track parameters for crossed material
295 //------------------------------------------------------------------
296 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
297 Double_t beta2=p2/(p2 + GetMass()*GetMass());
298 d*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
299
300 //Multiple scattering******************
301 if (d!=0) {
33329428 302 Double_t theta2=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(d);
303 //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*TMath::Abs(d)*9.36*2.33;
a9a2d814 304 fC22 += theta2*(1.- fP2*fP2)*(1. + fP3*fP3);
305 fC33 += theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
306 fC43 += theta2*fP3*fP4*(1. + fP3*fP3);
307 fC44 += theta2*fP3*fP4*fP3*fP4;
308 }
309
310 //Energy losses************************
311 if (x0!=0.) {
312 d*=x0;
313 Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d;
314 fP4*=(1.- sqrt(p2+GetMass()*GetMass())/p2*dE);
315 }
316
317 if (!Invariant()) return 0;
318
319 return 1;
320}
321
322//____________________________________________________________________________
323Int_t AliITStrackV2::PropagateTo(Double_t xk, Double_t d, Double_t x0) {
006b5f7f 324 //------------------------------------------------------------------
325 //This function propagates a track
326 //------------------------------------------------------------------
327 Double_t x1=fX, x2=xk, dx=x2-x1;
328 Double_t f1=fP2, f2=f1 + fP4*dx;
a9a2d814 329 if (TMath::Abs(f2) >= 0.9999) {
006b5f7f 330 Int_t n=GetNumberOfClusters();
7f6ddf58 331 if (n>kWARN)
332 cerr<<n<<" AliITStrackV2::PropagateTo: Propagation failed !\n";
006b5f7f 333 return 0;
334 }
335
336 Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
337
338 fP0 += dx*(f1+f2)/(r1+r2);
339 fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
340 fP2 += dx*fP4;
341
342 //f = F - 1
343
344 Double_t f02= dx/(r1*r1*r1);
345 Double_t f04=0.5*dx*dx/(r1*r1*r1);
346 Double_t f12= dx*fP3*f1/(r1*r1*r1);
347 Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1);
348 Double_t f13= dx/r1;
349 Double_t f24= dx;
350
351 //b = C*ft
352 Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
353 Double_t b02=f24*fC40;
354 Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
355 Double_t b12=f24*fC41;
356 Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
357 Double_t b22=f24*fC42;
358 Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
359 Double_t b42=f24*fC44;
360 Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
361 Double_t b32=f24*fC43;
362
363 //a = f*b = f*C*ft
364 Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
365 Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
366 Double_t a22=f24*b42;
367
368 //F*C*Ft = C + (b + bt + a)
369 fC00 += b00 + b00 + a00;
370 fC10 += b10 + b01 + a01;
371 fC20 += b20 + b02 + a02;
372 fC30 += b30;
373 fC40 += b40;
374 fC11 += b11 + b11 + a11;
375 fC21 += b21 + b12 + a12;
376 fC31 += b31;
377 fC41 += b41;
378 fC22 += b22 + b22 + a22;
379 fC32 += b32;
380 fC42 += b42;
381
382 fX=x2;
383
a9a2d814 384 if (!CorrectForMaterial(d,x0)) return 0;
006b5f7f 385
386 return 1;
387}
388
389//____________________________________________________________________________
390Int_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, UInt_t index) {
391 //------------------------------------------------------------------
392 //This function updates track parameters
393 //------------------------------------------------------------------
394 Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
395 Double_t c00=fC00;
396 Double_t c10=fC10, c11=fC11;
397 Double_t c20=fC20, c21=fC21, c22=fC22;
398 Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
399 Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
400
401
402 Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
403 r00+=fC00; r01+=fC10; r11+=fC11;
404 Double_t det=r00*r11 - r01*r01;
405 Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
406
407 Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
408 Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
409 Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
410 Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
411 Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
412
413 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
414 Double_t sf=fP2 + k20*dy + k21*dz;
7f6ddf58 415
006b5f7f 416 fP0 += k00*dy + k01*dz;
417 fP1 += k10*dy + k11*dz;
418 fP2 = sf;
419 fP3 += k30*dy + k31*dz;
420 fP4 += k40*dy + k41*dz;
421
422 Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
423 Double_t c12=fC21, c13=fC31, c14=fC41;
424
425 fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
426 fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
427 fC40-=k00*c04+k01*c14;
428
429 fC11-=k10*c01+k11*fC11;
430 fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
431 fC41-=k10*c04+k11*c14;
432
433 fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
434 fC42-=k20*c04+k21*c14;
435
436 fC33-=k30*c03+k31*c13;
437 fC43-=k30*c04+k31*c14;
438
439 fC44-=k40*c04+k41*c14;
440
441 if (!Invariant()) {
442 fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
443 fC00=c00;
444 fC10=c10; fC11=c11;
445 fC20=c20; fC21=c21; fC22=c22;
446 fC30=c30; fC31=c31; fC32=c32; fC33=c33;
447 fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
448 return 0;
449 }
450
451 Int_t n=GetNumberOfClusters();
452 fIndex[n]=index;
453 SetNumberOfClusters(n+1);
454 SetChi2(GetChi2()+chi2);
455
456 return 1;
457}
458
006b5f7f 459Int_t AliITStrackV2::Invariant() const {
460 //------------------------------------------------------------------
461 // This function is for debugging purpose only
462 //------------------------------------------------------------------
7f6ddf58 463 Int_t n=GetNumberOfClusters();
464
a9a2d814 465 if (TMath::Abs(fP2)>=0.9999){
466 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fP2="<<fP2<<endl;
467 return 0;
468 }
469 if (fC00<=0 || fC00>9.) {
470 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC00="<<fC00<<endl;
471 return 0;
472 }
473 if (fC11<=0 || fC11>9.) {
474 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC11="<<fC11<<endl;
475 return 0;
476 }
477 if (fC22<=0 || fC22>1.) {
478 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC22="<<fC22<<endl;
479 return 0;
480 }
481 if (fC33<=0 || fC33>1.) {
482 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC33="<<fC33<<endl;
483 return 0;
484 }
485 if (fC44<=0 || fC44>6e-5) {
486 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC44="<<fC44<<endl;
487 return 0;
14825d5a 488 }
14825d5a 489 return 1;
006b5f7f 490}
491
492//____________________________________________________________________________
a9a2d814 493Int_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) {
006b5f7f 494 //------------------------------------------------------------------
495 //This function propagates a track
496 //------------------------------------------------------------------
33329428 497 Double_t alpha=fAlpha, x=fX;
006b5f7f 498 Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
499 Double_t c00=fC00;
500 Double_t c10=fC10, c11=fC11;
501 Double_t c20=fC20, c21=fC21, c22=fC22;
502 Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
503 Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
504
33329428 505 if (alp < -TMath::Pi()) alp += 2*TMath::Pi();
506 else if (alp >= TMath::Pi()) alp -= 2*TMath::Pi();
507 Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha);
508 Double_t sf=fP2, cf=TMath::Sqrt(1.- fP2*fP2);
006b5f7f 509
33329428 510 TMatrixD *T=0;
511 // **** rotation **********************
512 {
006b5f7f 513 fAlpha = alp;
33329428 514 fX = x*ca + p0*sa;
515 fP0= -x*sa + p0*ca;
516 fP2= sf*ca - cf*sa;
517
518 TMatrixD C(5,5);
519 C(0,0)=c00;
520 C(1,0)=c10; C(1,1)=c11;
521 C(2,0)=c20; C(2,1)=c21; C(2,2)=c22;
522 C(3,0)=c30; C(3,1)=c31; C(3,2)=c32; C(3,3)=c33;
523 C(4,0)=c40; C(4,1)=c41; C(4,2)=c42; C(4,3)=c43; C(4,4)=c44;
524 C(0,1)=C(1,0);
525 C(0,2)=C(2,0); C(1,2)=C(2,1);
526 C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
527 C(0,4)=C(4,0); C(1,4)=C(4,1); C(2,4)=C(4,2); C(3,4)=C(4,3);
006b5f7f 528
33329428 529 TMatrixD F(6,5);
530 F(0,0)=sa;
531 F(1,0)=ca;
532 F(2,1)=F(4,3)=F(5,4)=1;
533 F(3,2)=ca + sf/cf*sa;
006b5f7f 534
33329428 535 TMatrixD tmp(C,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F));
536 T=new TMatrixD(F,TMatrixD::kMult,tmp);
537 }
006b5f7f 538
33329428 539 // **** translation ******************
540 {
541 Double_t dx=xk-fX;
006b5f7f 542 Double_t f1=fP2, f2=f1 + fP4*dx;
a9a2d814 543 if (TMath::Abs(f2) >= 0.9999) {
006b5f7f 544 Int_t n=GetNumberOfClusters();
7f6ddf58 545 if (n>kWARN)
546 cerr<<n<<" AliITStrackV2::Propagate: Propagation failed !\n";
006b5f7f 547 return 0;
548 }
33329428 549 Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
006b5f7f 550
33329428 551 fX=xk;
006b5f7f 552 fP0 += dx*(f1+f2)/(r1+r2);
553 fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
554 fP2 += dx*fP4;
555
33329428 556 TMatrixD F(5,6);
557 F(0,1)=F(1,2)=F(2,3)=F(3,4)=F(4,5)=1;
558 F(0,3)=dx/(r1+r2)*(2+(f1+f2)*(f2/r2+f1/r1)/(r1+r2));
559 F(0,5)=dx*dx/(r1+r2)*(1+(f1+f2)*f2/(r1+r2));
560 F(1,3)=dx*fP3/(f1*r2 + f2*r1)*(2-(f1+f2)*(r2-f1*f2/r2+r1-f2*f1/r1)/(f1*r2 + f2*r1));
561 F(1,4)=dx*(f1+f2)/(f1*r2 + f2*r1);
562 F(1,5)=dx*dx*fP3/(f1*r2 + f2*r1)*(1-(f1+f2)*(-f1*f2/r2+r1)/(f1*r2 + f2*r1));
563 F(2,5)=dx;
564 F(0,0)=-1/(r1+r2)*((f1+f2)+dx*fP4*(1+(f1+f2)/(r1+r2)*f2/r2));
565 F(1,0)=-fP3/(f1*r2 + f2*r1)*((f1+f2)+dx*fP4*(1+(f1+f2)/(f1*r2 + f2*r1)*(f1*f2/r2-r1)));
566 F(2,0)=-fP4;
567
568 TMatrixD tmp(*T,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F));
569 delete T;
570 TMatrixD C(F,TMatrixD::kMult,tmp);
006b5f7f 571
572 fC00=C(0,0);
573 fC10=C(1,0); fC11=C(1,1);
574 fC20=C(2,0); fC21=C(2,1); fC22=C(2,2);
575 fC30=C(3,0); fC31=C(3,1); fC32=C(3,2); fC33=C(3,3);
576 fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4);
577
006b5f7f 578 if (!Invariant()) {
33329428 579 fAlpha=alpha;
580 fX=x;
006b5f7f 581 fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
582 fC00=c00;
583 fC10=c10; fC11=c11;
584 fC20=c20; fC21=c21; fC22=c22;
585 fC30=c30; fC31=c31; fC32=c32; fC33=c33;
586 fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
587 return 0;
588 }
33329428 589 }
006b5f7f 590
591 return 1;
592}
593
a6a6f799 594Double_t AliITStrackV2::GetD(Double_t x, Double_t y) const {
006b5f7f 595 //------------------------------------------------------------------
a6a6f799 596 // This function calculates the transverse impact parameter
597 // with respect to a point with global coordinates (x,y)
006b5f7f 598 //------------------------------------------------------------------
a6a6f799 599 Double_t xt=fX, yt=fP0;
600
601 Double_t sn=TMath::Sin(fAlpha), cs=TMath::Cos(fAlpha);
602 Double_t a = x*cs + y*sn;
603 y = -x*sn + y*cs; x=a;
604 xt-=x; yt-=y;
605
606 sn=fP4*xt - fP2; cs=fP4*yt + TMath::Sqrt(1.- fP2*fP2);
607 a=2*(xt*fP2 - yt*TMath::Sqrt(1.- fP2*fP2))-fP4*(xt*xt + yt*yt);
006b5f7f 608 if (fP4<0) a=-a;
609 return a/(1 + TMath::Sqrt(sn*sn + cs*cs));
a9a2d814 610}
006b5f7f 611
a9a2d814 612Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
613 //------------------------------------------------------------------
614 //This function improves angular track parameters
615 //------------------------------------------------------------------
616 Double_t dy=fP0-yv, dz=fP1-zv;
617 Double_t r2=fX*fX+dy*dy;
618 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
619 Double_t beta2=p2/(p2 + GetMass()*GetMass());
620 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
621 //Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
622 Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
623 {
624 Double_t parp=0.5*(fP4*fX + dy*TMath::Sqrt(4/r2-fP4*fP4));
625 Double_t sigma2p = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
626 sigma2p += fC00/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
627 sigma2p += kSigmaYV*kSigmaYV/r2;
628 sigma2p += 0.25*fC44*fX*fX;
629 Double_t eps2p=sigma2p/(fC22+sigma2p);
630 fP0 += fC20/(fC22+sigma2p)*(parp-fP2);
631 fP2 = eps2p*fP2 + (1-eps2p)*parp;
632 fC22 *= eps2p;
633 fC20 *= eps2p;
634 }
635 {
636 Double_t parl=0.5*fP4*dz/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
637 Double_t sigma2l=theta2;
638 sigma2l += fC11/r2+fC00*dy*dy*dz*dz/(r2*r2*r2);
639 sigma2l += kSigmaZV*kSigmaZV/r2;
640 Double_t eps2l=sigma2l/(fC33+sigma2l);
641 fP1 += fC31/(fC33+sigma2l)*(parl-fP3);
642 fP4 += fC43/(fC33+sigma2l)*(parl-fP3);
643 fP3 = eps2l*fP3 + (1-eps2l)*parl;
644 fC33 *= eps2l; fC43 *= eps2l;
645 fC31 *= eps2l;
646 }
647 if (!Invariant()) return 0;
648 return 1;
649}
006b5f7f 650
a9a2d814 651/*
006b5f7f 652Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
653 //------------------------------------------------------------------
654 //This function improves angular track parameters
655 //------------------------------------------------------------------
656 Double_t dy=fP0-yv, dz=fP1-zv;
657 Double_t r2=fX*fX+dy*dy;
658 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
7f6ddf58 659 Double_t beta2=p2/(p2 + GetMass()*GetMass());
006b5f7f 660 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
a9a2d814 661 //Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
662 Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
006b5f7f 663
664 Double_t par=0.5*(fP4*fX + dy*TMath::Sqrt(4/r2-fP4*fP4));
665 Double_t sigma2 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
666 sigma2 += fC00/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
667 sigma2 += kSigmaYV*kSigmaYV/r2;
668 sigma2 += 0.25*fC44*fX*fX;
669 Double_t eps2=sigma2/(fC22+sigma2), eps=TMath::Sqrt(eps2);
670 if (10*r2*fC44<fC22) {
671 fP2 = eps2*fP2 + (1-eps2)*par;
672 fC22*=eps2; fC21*=eps; fC20*=eps; fC32*=eps; fC42*=eps;
673 }
674
675 par=0.5*fP4*dz/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
676 sigma2=theta2;
677 sigma2 += fC11/r2+fC00*dy*dy*dz*dz/(r2*r2*r2);
678 sigma2 += kSigmaZV*kSigmaZV/r2;
679 eps2=sigma2/(fC33+sigma2); eps=TMath::Sqrt(eps2);
680 Double_t tgl=fP3;
681 fP3 = eps2*fP3 + (1-eps2)*par;
682 fC33*=eps2; fC32*=eps; fC31*=eps; fC30*=eps; fC43*=eps;
683
684 eps=TMath::Sqrt((1+fP3*fP3)/(1+tgl*tgl));
685 fP4*=eps;
686 fC44*=eps*eps; fC43*=eps;fC42*=eps; fC41*=eps; fC40*=eps;
687
688 if (!Invariant()) return 0;
689 return 1;
690}
006b5f7f 691*/
14825d5a 692void AliITStrackV2::ResetCovariance() {
693 //------------------------------------------------------------------
694 //This function makes a track forget its history :)
695 //------------------------------------------------------------------
696
697 fC00*=10.;
698 fC10=0.; fC11*=10.;
699 fC20=0.; fC21=0.; fC22*=10.;
700 fC30=0.; fC31=0.; fC32=0.; fC33*=10.;
701 fC40=0.; fC41=0.; fC42=0.; fC43=0.; fC44*=10.;
702
703}
23efe5f1 704
705void AliITStrackV2::CookdEdx(Double_t low, Double_t up) {
706 //-----------------------------------------------------------------
a9a2d814 707 // This function calculates dE/dX within the "low" and "up" cuts.
708 // Origin: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
23efe5f1 709 //-----------------------------------------------------------------
710 Int_t i;
a9a2d814 711 Int_t nc=4;
23efe5f1 712 // The clusters order is: SSD-2, SSD-1, SDD-2, SDD-1, SPD-2, SPD-1
713 // Take only SSD and SDD
23efe5f1 714
a9a2d814 715 Int_t swap;//stupid sorting
23efe5f1 716 do {
717 swap=0;
718 for (i=0; i<nc-1; i++) {
719 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
720 Float_t tmp=fdEdxSample[i];
721 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
722 swap++;
723 }
724 } while (swap);
725
a9a2d814 726 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc); //b.b. to take two lowest dEdX
727 // values from four ones choose
728 // nu=2
23efe5f1 729 Float_t dedx=0;
a9a2d814 730 for (i=nl; i<nu; i++) dedx += fdEdxSample[i];
731 dedx /= (nu-nl);
23efe5f1 732
23efe5f1 733 SetdEdx(dedx);
734}