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