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