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
22b13da0 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}
006b5f7f 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
f29dbb4b 336 // old position [SR, GSI, 17.02.2003]
337 Double_t oldX = fX, oldY = fP0, oldZ = fP1;
338
006b5f7f 339 Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
340
341 fP0 += dx*(f1+f2)/(r1+r2);
342 fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
343 fP2 += dx*fP4;
344
345 //f = F - 1
346
347 Double_t f02= dx/(r1*r1*r1);
348 Double_t f04=0.5*dx*dx/(r1*r1*r1);
349 Double_t f12= dx*fP3*f1/(r1*r1*r1);
350 Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1);
351 Double_t f13= dx/r1;
352 Double_t f24= dx;
353
354 //b = C*ft
355 Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
356 Double_t b02=f24*fC40;
357 Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
358 Double_t b12=f24*fC41;
359 Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
360 Double_t b22=f24*fC42;
361 Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
362 Double_t b42=f24*fC44;
363 Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
364 Double_t b32=f24*fC43;
365
366 //a = f*b = f*C*ft
367 Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
368 Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
369 Double_t a22=f24*b42;
370
371 //F*C*Ft = C + (b + bt + a)
372 fC00 += b00 + b00 + a00;
373 fC10 += b10 + b01 + a01;
374 fC20 += b20 + b02 + a02;
375 fC30 += b30;
376 fC40 += b40;
377 fC11 += b11 + b11 + a11;
378 fC21 += b21 + b12 + a12;
379 fC31 += b31;
380 fC41 += b41;
381 fC22 += b22 + b22 + a22;
382 fC32 += b32;
383 fC42 += b42;
384
385 fX=x2;
386
a9a2d814 387 if (!CorrectForMaterial(d,x0)) return 0;
006b5f7f 388
f29dbb4b 389 // Integrated Time [SR, GSI, 17.02.2003]
390 if (IsStartedTimeIntegral()) {
391 Double_t l2 = (fX-oldX)*(fX-oldX)+(fP0-oldY)*(fP0-oldY)+(fP1-oldZ)*(fP1-oldZ);
392 AddTimeStep(TMath::Sqrt(l2));
393 }
394 //
395
006b5f7f 396 return 1;
397}
398
399//____________________________________________________________________________
400Int_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, UInt_t index) {
401 //------------------------------------------------------------------
402 //This function updates track parameters
403 //------------------------------------------------------------------
404 Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
405 Double_t c00=fC00;
406 Double_t c10=fC10, c11=fC11;
407 Double_t c20=fC20, c21=fC21, c22=fC22;
408 Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
409 Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
410
411
412 Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
413 r00+=fC00; r01+=fC10; r11+=fC11;
414 Double_t det=r00*r11 - r01*r01;
415 Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
416
417 Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
418 Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
419 Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
420 Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
421 Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
422
423 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
424 Double_t sf=fP2 + k20*dy + k21*dz;
7f6ddf58 425
006b5f7f 426 fP0 += k00*dy + k01*dz;
427 fP1 += k10*dy + k11*dz;
428 fP2 = sf;
429 fP3 += k30*dy + k31*dz;
430 fP4 += k40*dy + k41*dz;
431
432 Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
433 Double_t c12=fC21, c13=fC31, c14=fC41;
434
435 fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
436 fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
437 fC40-=k00*c04+k01*c14;
438
439 fC11-=k10*c01+k11*fC11;
440 fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
441 fC41-=k10*c04+k11*c14;
442
443 fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
444 fC42-=k20*c04+k21*c14;
445
446 fC33-=k30*c03+k31*c13;
447 fC43-=k30*c04+k31*c14;
448
449 fC44-=k40*c04+k41*c14;
450
451 if (!Invariant()) {
452 fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
453 fC00=c00;
454 fC10=c10; fC11=c11;
455 fC20=c20; fC21=c21; fC22=c22;
456 fC30=c30; fC31=c31; fC32=c32; fC33=c33;
457 fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
458 return 0;
459 }
460
461 Int_t n=GetNumberOfClusters();
462 fIndex[n]=index;
463 SetNumberOfClusters(n+1);
464 SetChi2(GetChi2()+chi2);
465
466 return 1;
467}
468
006b5f7f 469Int_t AliITStrackV2::Invariant() const {
470 //------------------------------------------------------------------
471 // This function is for debugging purpose only
472 //------------------------------------------------------------------
7f6ddf58 473 Int_t n=GetNumberOfClusters();
474
a9a2d814 475 if (TMath::Abs(fP2)>=0.9999){
476 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fP2="<<fP2<<endl;
477 return 0;
478 }
479 if (fC00<=0 || fC00>9.) {
480 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC00="<<fC00<<endl;
481 return 0;
482 }
483 if (fC11<=0 || fC11>9.) {
484 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC11="<<fC11<<endl;
485 return 0;
486 }
487 if (fC22<=0 || fC22>1.) {
488 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC22="<<fC22<<endl;
489 return 0;
490 }
491 if (fC33<=0 || fC33>1.) {
492 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC33="<<fC33<<endl;
493 return 0;
494 }
495 if (fC44<=0 || fC44>6e-5) {
496 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC44="<<fC44<<endl;
497 return 0;
14825d5a 498 }
14825d5a 499 return 1;
006b5f7f 500}
501
502//____________________________________________________________________________
a9a2d814 503Int_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) {
006b5f7f 504 //------------------------------------------------------------------
505 //This function propagates a track
506 //------------------------------------------------------------------
33329428 507 Double_t alpha=fAlpha, x=fX;
006b5f7f 508 Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
509 Double_t c00=fC00;
510 Double_t c10=fC10, c11=fC11;
511 Double_t c20=fC20, c21=fC21, c22=fC22;
512 Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
513 Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
514
33329428 515 if (alp < -TMath::Pi()) alp += 2*TMath::Pi();
516 else if (alp >= TMath::Pi()) alp -= 2*TMath::Pi();
517 Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha);
518 Double_t sf=fP2, cf=TMath::Sqrt(1.- fP2*fP2);
006b5f7f 519
33329428 520 TMatrixD *T=0;
521 // **** rotation **********************
522 {
006b5f7f 523 fAlpha = alp;
33329428 524 fX = x*ca + p0*sa;
525 fP0= -x*sa + p0*ca;
526 fP2= sf*ca - cf*sa;
527
528 TMatrixD C(5,5);
529 C(0,0)=c00;
530 C(1,0)=c10; C(1,1)=c11;
531 C(2,0)=c20; C(2,1)=c21; C(2,2)=c22;
532 C(3,0)=c30; C(3,1)=c31; C(3,2)=c32; C(3,3)=c33;
533 C(4,0)=c40; C(4,1)=c41; C(4,2)=c42; C(4,3)=c43; C(4,4)=c44;
534 C(0,1)=C(1,0);
535 C(0,2)=C(2,0); C(1,2)=C(2,1);
536 C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
537 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 538
33329428 539 TMatrixD F(6,5);
540 F(0,0)=sa;
541 F(1,0)=ca;
542 F(2,1)=F(4,3)=F(5,4)=1;
543 F(3,2)=ca + sf/cf*sa;
006b5f7f 544
33329428 545 TMatrixD tmp(C,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F));
546 T=new TMatrixD(F,TMatrixD::kMult,tmp);
547 }
006b5f7f 548
33329428 549 // **** translation ******************
550 {
551 Double_t dx=xk-fX;
006b5f7f 552 Double_t f1=fP2, f2=f1 + fP4*dx;
a9a2d814 553 if (TMath::Abs(f2) >= 0.9999) {
006b5f7f 554 Int_t n=GetNumberOfClusters();
7f6ddf58 555 if (n>kWARN)
556 cerr<<n<<" AliITStrackV2::Propagate: Propagation failed !\n";
006b5f7f 557 return 0;
558 }
33329428 559 Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
006b5f7f 560
33329428 561 fX=xk;
006b5f7f 562 fP0 += dx*(f1+f2)/(r1+r2);
563 fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
564 fP2 += dx*fP4;
565
33329428 566 TMatrixD F(5,6);
567 F(0,1)=F(1,2)=F(2,3)=F(3,4)=F(4,5)=1;
568 F(0,3)=dx/(r1+r2)*(2+(f1+f2)*(f2/r2+f1/r1)/(r1+r2));
569 F(0,5)=dx*dx/(r1+r2)*(1+(f1+f2)*f2/(r1+r2));
570 F(1,3)=dx*fP3/(f1*r2 + f2*r1)*(2-(f1+f2)*(r2-f1*f2/r2+r1-f2*f1/r1)/(f1*r2 + f2*r1));
571 F(1,4)=dx*(f1+f2)/(f1*r2 + f2*r1);
572 F(1,5)=dx*dx*fP3/(f1*r2 + f2*r1)*(1-(f1+f2)*(-f1*f2/r2+r1)/(f1*r2 + f2*r1));
573 F(2,5)=dx;
574 F(0,0)=-1/(r1+r2)*((f1+f2)+dx*fP4*(1+(f1+f2)/(r1+r2)*f2/r2));
575 F(1,0)=-fP3/(f1*r2 + f2*r1)*((f1+f2)+dx*fP4*(1+(f1+f2)/(f1*r2 + f2*r1)*(f1*f2/r2-r1)));
576 F(2,0)=-fP4;
577
578 TMatrixD tmp(*T,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F));
579 delete T;
580 TMatrixD C(F,TMatrixD::kMult,tmp);
006b5f7f 581
582 fC00=C(0,0);
583 fC10=C(1,0); fC11=C(1,1);
584 fC20=C(2,0); fC21=C(2,1); fC22=C(2,2);
585 fC30=C(3,0); fC31=C(3,1); fC32=C(3,2); fC33=C(3,3);
586 fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4);
587
006b5f7f 588 if (!Invariant()) {
33329428 589 fAlpha=alpha;
590 fX=x;
006b5f7f 591 fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
592 fC00=c00;
593 fC10=c10; fC11=c11;
594 fC20=c20; fC21=c21; fC22=c22;
595 fC30=c30; fC31=c31; fC32=c32; fC33=c33;
596 fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
597 return 0;
598 }
33329428 599 }
006b5f7f 600
601 return 1;
602}
603
a6a6f799 604Double_t AliITStrackV2::GetD(Double_t x, Double_t y) const {
006b5f7f 605 //------------------------------------------------------------------
a6a6f799 606 // This function calculates the transverse impact parameter
607 // with respect to a point with global coordinates (x,y)
006b5f7f 608 //------------------------------------------------------------------
a6a6f799 609 Double_t xt=fX, yt=fP0;
610
611 Double_t sn=TMath::Sin(fAlpha), cs=TMath::Cos(fAlpha);
612 Double_t a = x*cs + y*sn;
613 y = -x*sn + y*cs; x=a;
614 xt-=x; yt-=y;
615
616 sn=fP4*xt - fP2; cs=fP4*yt + TMath::Sqrt(1.- fP2*fP2);
617 a=2*(xt*fP2 - yt*TMath::Sqrt(1.- fP2*fP2))-fP4*(xt*xt + yt*yt);
006b5f7f 618 if (fP4<0) a=-a;
619 return a/(1 + TMath::Sqrt(sn*sn + cs*cs));
a9a2d814 620}
006b5f7f 621
a9a2d814 622Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
623 //------------------------------------------------------------------
624 //This function improves angular track parameters
625 //------------------------------------------------------------------
626 Double_t dy=fP0-yv, dz=fP1-zv;
627 Double_t r2=fX*fX+dy*dy;
628 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
629 Double_t beta2=p2/(p2 + GetMass()*GetMass());
630 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
631 //Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
632 Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
633 {
634 Double_t parp=0.5*(fP4*fX + dy*TMath::Sqrt(4/r2-fP4*fP4));
635 Double_t sigma2p = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
636 sigma2p += fC00/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
637 sigma2p += kSigmaYV*kSigmaYV/r2;
638 sigma2p += 0.25*fC44*fX*fX;
639 Double_t eps2p=sigma2p/(fC22+sigma2p);
640 fP0 += fC20/(fC22+sigma2p)*(parp-fP2);
641 fP2 = eps2p*fP2 + (1-eps2p)*parp;
642 fC22 *= eps2p;
643 fC20 *= eps2p;
644 }
645 {
646 Double_t parl=0.5*fP4*dz/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
647 Double_t sigma2l=theta2;
648 sigma2l += fC11/r2+fC00*dy*dy*dz*dz/(r2*r2*r2);
649 sigma2l += kSigmaZV*kSigmaZV/r2;
650 Double_t eps2l=sigma2l/(fC33+sigma2l);
651 fP1 += fC31/(fC33+sigma2l)*(parl-fP3);
652 fP4 += fC43/(fC33+sigma2l)*(parl-fP3);
653 fP3 = eps2l*fP3 + (1-eps2l)*parl;
654 fC33 *= eps2l; fC43 *= eps2l;
655 fC31 *= eps2l;
656 }
657 if (!Invariant()) return 0;
658 return 1;
659}
006b5f7f 660
a9a2d814 661/*
006b5f7f 662Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
663 //------------------------------------------------------------------
664 //This function improves angular track parameters
665 //------------------------------------------------------------------
666 Double_t dy=fP0-yv, dz=fP1-zv;
667 Double_t r2=fX*fX+dy*dy;
668 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
7f6ddf58 669 Double_t beta2=p2/(p2 + GetMass()*GetMass());
006b5f7f 670 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
a9a2d814 671 //Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
672 Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
006b5f7f 673
674 Double_t par=0.5*(fP4*fX + dy*TMath::Sqrt(4/r2-fP4*fP4));
675 Double_t sigma2 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
676 sigma2 += fC00/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
677 sigma2 += kSigmaYV*kSigmaYV/r2;
678 sigma2 += 0.25*fC44*fX*fX;
679 Double_t eps2=sigma2/(fC22+sigma2), eps=TMath::Sqrt(eps2);
680 if (10*r2*fC44<fC22) {
681 fP2 = eps2*fP2 + (1-eps2)*par;
682 fC22*=eps2; fC21*=eps; fC20*=eps; fC32*=eps; fC42*=eps;
683 }
684
685 par=0.5*fP4*dz/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
686 sigma2=theta2;
687 sigma2 += fC11/r2+fC00*dy*dy*dz*dz/(r2*r2*r2);
688 sigma2 += kSigmaZV*kSigmaZV/r2;
689 eps2=sigma2/(fC33+sigma2); eps=TMath::Sqrt(eps2);
690 Double_t tgl=fP3;
691 fP3 = eps2*fP3 + (1-eps2)*par;
692 fC33*=eps2; fC32*=eps; fC31*=eps; fC30*=eps; fC43*=eps;
693
694 eps=TMath::Sqrt((1+fP3*fP3)/(1+tgl*tgl));
695 fP4*=eps;
696 fC44*=eps*eps; fC43*=eps;fC42*=eps; fC41*=eps; fC40*=eps;
697
698 if (!Invariant()) return 0;
699 return 1;
700}
006b5f7f 701*/
14825d5a 702void AliITStrackV2::ResetCovariance() {
703 //------------------------------------------------------------------
704 //This function makes a track forget its history :)
705 //------------------------------------------------------------------
706
707 fC00*=10.;
708 fC10=0.; fC11*=10.;
709 fC20=0.; fC21=0.; fC22*=10.;
710 fC30=0.; fC31=0.; fC32=0.; fC33*=10.;
711 fC40=0.; fC41=0.; fC42=0.; fC43=0.; fC44*=10.;
712
713}
23efe5f1 714
715void AliITStrackV2::CookdEdx(Double_t low, Double_t up) {
716 //-----------------------------------------------------------------
a9a2d814 717 // This function calculates dE/dX within the "low" and "up" cuts.
718 // Origin: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
23efe5f1 719 //-----------------------------------------------------------------
720 Int_t i;
a9a2d814 721 Int_t nc=4;
23efe5f1 722 // The clusters order is: SSD-2, SSD-1, SDD-2, SDD-1, SPD-2, SPD-1
723 // Take only SSD and SDD
23efe5f1 724
a9a2d814 725 Int_t swap;//stupid sorting
23efe5f1 726 do {
727 swap=0;
728 for (i=0; i<nc-1; i++) {
729 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
730 Float_t tmp=fdEdxSample[i];
731 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
732 swap++;
733 }
734 } while (swap);
735
a9a2d814 736 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc); //b.b. to take two lowest dEdX
737 // values from four ones choose
738 // nu=2
23efe5f1 739 Float_t dedx=0;
a9a2d814 740 for (i=nl; i<nu; i++) dedx += fdEdxSample[i];
741 dedx /= (nu-nl);
23efe5f1 742
23efe5f1 743 SetdEdx(dedx);
744}