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