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