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