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