Using TMath functions
[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>
4ae5bbc4 26#include <Riostream.h>
006b5f7f 27
28#include "AliCluster.h"
517b130f 29#include "AliTPCtrack.h"
006b5f7f 30#include "AliITStrackV2.h"
31
32ClassImp(AliITStrackV2)
33
7f6ddf58 34const Int_t kWARN=5;
1a947abf 35Double_t AliITStrackV2::fSigmaYV = 0.005;
36Double_t AliITStrackV2::fSigmaZV = 0.01;
006b5f7f 37//____________________________________________________________________________
e6ba4672 38AliITStrackV2::AliITStrackV2():AliKalmanTrack(),
22b13da0 39 fX(0),
40 fAlpha(0),
41 fdEdx(0),
42 fP0(0),
43 fP1(0),
44 fP2(0),
45 fP3(0),
46 fP4(0),
47 fC00(0),
48 fC10(0),
49 fC11(0),
50 fC20(0),
51 fC21(0),
52 fC22(0),
53 fC30(0),
54 fC31(0),
55 fC32(0),
56 fC33(0),
57 fC40(0),
58 fC41(0),
59 fC42(0),
60 fC43(0),
e6ba4672 61 fC44(0)
62 {
22b13da0 63 for(Int_t i=0; i<kMaxLayer; i++) fIndex[i]=0;
64 for(Int_t i=0; i<4; i++) fdEdxSample[i]=0;
65}
66//____________________________________________________________________________
e1a65e92 67AliITStrackV2::AliITStrackV2(const AliTPCtrack& t) throw (const Char_t *) :
68AliKalmanTrack(t) {
006b5f7f 69 //------------------------------------------------------------------
9b280d80 70 //Conversion TPC track -> ITS track
006b5f7f 71 //------------------------------------------------------------------
006b5f7f 72 SetChi2(0.);
73 SetNumberOfClusters(0);
14825d5a 74
7f6ddf58 75 fdEdx = t.GetdEdx();
76 SetMass(t.GetMass());
77
006b5f7f 78 fAlpha = t.GetAlpha();
79 if (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
80 else if (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
81
9b280d80 82 //Conversion of the track parameters
006b5f7f 83 Double_t x,p[5]; t.GetExternalParameters(x,p);
9b280d80 84 fX=x; x=GetConvConst();
006b5f7f 85 fP0=p[0];
86 fP1=p[1];
87 fP2=p[2];
88 fP3=p[3];
89 fP4=p[4]/x;
90
9b280d80 91 //Conversion of the covariance matrix
006b5f7f 92 Double_t c[15]; t.GetExternalCovariance(c);
93
94 fC00=c[0 ];
95 fC10=c[1 ]; fC11=c[2 ];
96 fC20=c[3 ]; fC21=c[4 ]; fC22=c[5 ];
97 fC30=c[6 ]; fC31=c[7 ]; fC32=c[8 ]; fC33=c[9 ];
98 fC40=c[10]/x; fC41=c[11]/x; fC42=c[12]/x; fC43=c[13]/x; fC44=c[14]/x/x;
99
100 if (!Invariant()) throw "AliITStrackV2: conversion failed !\n";
649ab0de 101
006b5f7f 102}
103
104//____________________________________________________________________________
105AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) : AliKalmanTrack(t) {
106 //------------------------------------------------------------------
107 //Copy constructor
108 //------------------------------------------------------------------
109 fX=t.fX;
110 fAlpha=t.fAlpha;
111 fdEdx=t.fdEdx;
112
113 fP0=t.fP0; fP1=t.fP1; fP2=t.fP2; fP3=t.fP3; fP4=t.fP4;
114
115 fC00=t.fC00;
116 fC10=t.fC10; fC11=t.fC11;
117 fC20=t.fC20; fC21=t.fC21; fC22=t.fC22;
118 fC30=t.fC30; fC31=t.fC31; fC32=t.fC32; fC33=t.fC33;
119 fC40=t.fC40; fC41=t.fC41; fC42=t.fC42; fC43=t.fC43; fC44=t.fC44;
120
121 Int_t n=GetNumberOfClusters();
23efe5f1 122 for (Int_t i=0; i<n; i++) {
a9a2d814 123 fIndex[i]=t.fIndex[i];
124 if (i<4) fdEdxSample[i]=t.fdEdxSample[i];
23efe5f1 125 }
006b5f7f 126}
a9a2d814 127
006b5f7f 128//_____________________________________________________________________________
129Int_t AliITStrackV2::Compare(const TObject *o) const {
130 //-----------------------------------------------------------------
131 // This function compares tracks according to the their curvature
132 //-----------------------------------------------------------------
7f6ddf58 133 AliITStrackV2 *t=(AliITStrackV2*)o;
a9a2d814 134 //Double_t co=TMath::Abs(t->Get1Pt());
135 //Double_t c =TMath::Abs(Get1Pt());
136 Double_t co=t->GetSigmaY2()*t->GetSigmaZ2();
137 Double_t c =GetSigmaY2()*GetSigmaZ2();
006b5f7f 138 if (c>co) return 1;
139 else if (c<co) return -1;
140 return 0;
141}
142
143//_____________________________________________________________________________
144void AliITStrackV2::GetExternalCovariance(Double_t cc[15]) const {
145 //-------------------------------------------------------------------------
146 // This function returns an external representation of the covriance matrix.
147 // (See comments in AliTPCtrack.h about external track representation)
148 //-------------------------------------------------------------------------
9b280d80 149 Double_t a=GetConvConst();
006b5f7f 150
151 cc[0 ]=fC00;
152 cc[1 ]=fC10; cc[2 ]=fC11;
153 cc[3 ]=fC20; cc[4 ]=fC21; cc[5 ]=fC22;
154 cc[6 ]=fC30; cc[7 ]=fC31; cc[8 ]=fC32; cc[9 ]=fC33;
155 cc[10]=fC40*a; cc[11]=fC41*a; cc[12]=fC42*a; cc[13]=fC43*a; cc[14]=fC44*a*a;
156}
157
158//____________________________________________________________________________
a9a2d814 159Int_t AliITStrackV2::PropagateToVertex(Double_t d,Double_t x0) {
006b5f7f 160 //------------------------------------------------------------------
161 //This function propagates a track to the minimal distance from the origin
162 //------------------------------------------------------------------
5b316485 163 //Double_t xv=fP2*(fX*fP2 - fP0*TMath::Sqrt(1.- fP2*fP2)); //linear approxim.
164 Double_t tgf=-(fP4*fX - fP2)/(fP4*fP0 + TMath::Sqrt(1 - fP2*fP2));
165 Double_t snf=tgf/TMath::Sqrt(1.+ tgf*tgf);
166 Double_t xv=(snf - fP2)/fP4 + fX;
167 return PropagateTo(xv,d,x0);
006b5f7f 168}
169
170//____________________________________________________________________________
171Int_t AliITStrackV2::
172GetGlobalXYZat(Double_t xk, Double_t &x, Double_t &y, Double_t &z) const {
173 //------------------------------------------------------------------
174 //This function returns a track position in the global system
175 //------------------------------------------------------------------
176 Double_t dx=xk-fX;
177 Double_t f1=fP2, f2=f1 + fP4*dx;
a9a2d814 178 if (TMath::Abs(f2) >= 0.9999) {
7f6ddf58 179 Int_t n=GetNumberOfClusters();
180 if (n>kWARN)
181 cerr<<n<<" AliITStrackV2::GetGlobalXYZat: Propagation failed !\n";
006b5f7f 182 return 0;
183 }
184
185 Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
186
187 Double_t yk = fP0 + dx*(f1+f2)/(r1+r2);
188 Double_t zk = fP1 + dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
189
190 Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
191 x = xk*cs - yk*sn;
192 y = xk*sn + yk*cs;
193 z = zk;
194
195 return 1;
196}
197
198//_____________________________________________________________________________
199Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c) const
200{
201 //-----------------------------------------------------------------
202 // This function calculates a predicted chi2 increment.
203 //-----------------------------------------------------------------
204 Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
205 r00+=fC00; r01+=fC10; r11+=fC11;
206
207 Double_t det=r00*r11 - r01*r01;
7f6ddf58 208 if (TMath::Abs(det) < 1.e-30) {
006b5f7f 209 Int_t n=GetNumberOfClusters();
7f6ddf58 210 if (n>kWARN)
211 cerr<<n<<" AliKalmanTrack::GetPredictedChi2: Singular matrix !\n";
006b5f7f 212 return 1e10;
213 }
214 Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
33329428 215
006b5f7f 216 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
33329428 217
006b5f7f 218 return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
219}
220
221//_____________________________________________________________________________
222Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c,Double_t *m,
7f6ddf58 223Double_t x0) const {
006b5f7f 224 //-----------------------------------------------------------------
225 // This function calculates a chi2 increment with a vertex contraint
226 //-----------------------------------------------------------------
227 TVectorD x(5); x(0)=fP0; x(1)=fP1; x(2)=fP2; x(3)=fP3; x(4)=fP4;
228 TMatrixD C(5,5);
229 C(0,0)=fC00;
230 C(1,0)=fC10; C(1,1)=fC11;
231 C(2,0)=fC20; C(2,1)=fC21; C(2,2)=fC22;
232 C(3,0)=fC30; C(3,1)=fC31; C(3,2)=fC32; C(3,3)=fC33;
233 C(4,0)=fC40; C(4,1)=fC41; C(4,2)=fC42; C(4,3)=fC43; C(4,4)=fC44;
234
235 C(0,1)=C(1,0);
236 C(0,2)=C(2,0); C(1,2)=C(2,1);
237 C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
238 C(0,4)=C(4,0); C(1,4)=C(4,1); C(2,4)=C(4,2); C(3,4)=C(4,3);
239
240 TMatrixD H(4,5); H.UnitMatrix();
241 Double_t dy=(c->GetY() - m[0]), dz=(c->GetZ() - m[1]);
242
243 Double_t dr=TMath::Sqrt(fX*fX + dy*dy);
244 Double_t r =TMath::Sqrt(4/dr/dr - fP4*fP4);
245 Double_t sn=0.5*(fP4*fX + dy*r);
246 Double_t tg=0.5*fP4*dz/TMath::ASin(0.5*fP4*dr);
247 TVectorD mm(4);
248 mm(0)=m[0]=c->GetY(); mm(1)=m[1]=c->GetZ(); mm(2)=m[2]=sn; mm(3)=m[3]=tg;
249
250 Double_t v22=0.,v33=0.;
251 //x0=0.;
252 if (x0!=0.) {
253 Double_t pp2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
7f6ddf58 254 Double_t beta2=pp2/(pp2 + GetMass()*GetMass());
006b5f7f 255 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
256 Double_t theta2=14.1*14.1/(beta2*pp2*1e6)*x0;
257 v22 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
258 v33 = theta2*(1.+ GetTgl()*GetTgl())*(1. + GetTgl()*GetTgl());
259 }
260 Double_t sy2=c->GetSigmaY2(), sz2=c->GetSigmaZ2();
1a947abf 261 v22+=fSigmaYV*fSigmaYV/dr/dr;
006b5f7f 262 v22+=sy2/dr/dr;
263 Double_t v20=sy2/dr;
264
1a947abf 265 v33+=fSigmaZV*fSigmaZV/dr/dr;
006b5f7f 266 v33+=sz2/dr/dr;
267 Double_t v31=sz2/dr;
268
269 TMatrixD V(4,4);
270 V(0,0)=m[4 ]=sy2; V(0,1)=m[5 ]=0.; V(0,2)=m[6 ]=v20; V(0,3)=m[7 ]=0.;
271 V(1,0)=m[8 ]=0.; V(1,1)=m[9 ]=sz2; V(1,2)=m[10]=0.; V(1,3)=m[11]=v31;
272 V(2,0)=m[12]=v20; V(2,1)=m[13]=0.; V(2,2)=m[14]=v22; V(2,3)=m[15]=0.;
273 V(3,0)=m[16]=0.; V(3,1)=m[17]=v31; V(3,2)=m[18]=0.; V(3,3)=m[19]=v33;
274
275 TVectorD res=x; res*=H; res-=mm; //res*=-1;
276 TMatrixD tmp(H,TMatrixD::kMult,C);
277 TMatrixD R(tmp,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed,H)); R+=V;
278
279 Double_t det=R.Determinant();
7f6ddf58 280 if (TMath::Abs(det) < 1.e-30) {
006b5f7f 281 Int_t n=GetNumberOfClusters();
7f6ddf58 282 if (n>kWARN)
283 cerr<<n<<" AliITStrackV2::GetPredictedChi2: Singular matrix !\n";
006b5f7f 284 return 1e10;
285 }
286
287 R.Invert();
288
289 TVectorD rs=res;
290 res*=R;
291 return rs*res;
292}
293
294//____________________________________________________________________________
a9a2d814 295Int_t AliITStrackV2::CorrectForMaterial(Double_t d, Double_t x0) {
296 //------------------------------------------------------------------
297 //This function corrects the track parameters for crossed material
298 //------------------------------------------------------------------
299 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
300 Double_t beta2=p2/(p2 + GetMass()*GetMass());
301 d*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
302
303 //Multiple scattering******************
304 if (d!=0) {
33329428 305 Double_t theta2=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(d);
306 //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*TMath::Abs(d)*9.36*2.33;
a9a2d814 307 fC22 += theta2*(1.- fP2*fP2)*(1. + fP3*fP3);
308 fC33 += theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
309 fC43 += theta2*fP3*fP4*(1. + fP3*fP3);
310 fC44 += theta2*fP3*fP4*fP3*fP4;
311 }
312
313 //Energy losses************************
314 if (x0!=0.) {
315 d*=x0;
316 Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d;
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)
335 cerr<<n<<" AliITStrackV2::PropagateTo: Propagation failed !\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]
393 if (IsStartedTimeIntegral()) {
394 Double_t l2 = (fX-oldX)*(fX-oldX)+(fP0-oldY)*(fP0-oldY)+(fP1-oldZ)*(fP1-oldZ);
395 AddTimeStep(TMath::Sqrt(l2));
396 }
397 //
398
006b5f7f 399 return 1;
400}
401
402//____________________________________________________________________________
403Int_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, UInt_t index) {
404 //------------------------------------------------------------------
405 //This function updates track parameters
406 //------------------------------------------------------------------
407 Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
408 Double_t c00=fC00;
409 Double_t c10=fC10, c11=fC11;
410 Double_t c20=fC20, c21=fC21, c22=fC22;
411 Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
412 Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
413
414
415 Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
416 r00+=fC00; r01+=fC10; r11+=fC11;
417 Double_t det=r00*r11 - r01*r01;
418 Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
419
420 Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
421 Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
422 Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
423 Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
424 Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
425
426 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
427 Double_t sf=fP2 + k20*dy + k21*dz;
7f6ddf58 428
006b5f7f 429 fP0 += k00*dy + k01*dz;
430 fP1 += k10*dy + k11*dz;
431 fP2 = sf;
432 fP3 += k30*dy + k31*dz;
433 fP4 += k40*dy + k41*dz;
434
435 Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
436 Double_t c12=fC21, c13=fC31, c14=fC41;
437
438 fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
439 fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
440 fC40-=k00*c04+k01*c14;
441
442 fC11-=k10*c01+k11*fC11;
443 fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
444 fC41-=k10*c04+k11*c14;
445
446 fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
447 fC42-=k20*c04+k21*c14;
448
449 fC33-=k30*c03+k31*c13;
450 fC43-=k30*c04+k31*c14;
451
452 fC44-=k40*c04+k41*c14;
453
454 if (!Invariant()) {
455 fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
456 fC00=c00;
457 fC10=c10; fC11=c11;
458 fC20=c20; fC21=c21; fC22=c22;
459 fC30=c30; fC31=c31; fC32=c32; fC33=c33;
460 fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
461 return 0;
462 }
463
464 Int_t n=GetNumberOfClusters();
465 fIndex[n]=index;
466 SetNumberOfClusters(n+1);
467 SetChi2(GetChi2()+chi2);
468
469 return 1;
470}
471
006b5f7f 472Int_t AliITStrackV2::Invariant() const {
473 //------------------------------------------------------------------
474 // This function is for debugging purpose only
475 //------------------------------------------------------------------
7f6ddf58 476 Int_t n=GetNumberOfClusters();
477
a9a2d814 478 if (TMath::Abs(fP2)>=0.9999){
479 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fP2="<<fP2<<endl;
480 return 0;
481 }
482 if (fC00<=0 || fC00>9.) {
483 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC00="<<fC00<<endl;
484 return 0;
485 }
486 if (fC11<=0 || fC11>9.) {
487 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC11="<<fC11<<endl;
488 return 0;
489 }
490 if (fC22<=0 || fC22>1.) {
491 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC22="<<fC22<<endl;
492 return 0;
493 }
494 if (fC33<=0 || fC33>1.) {
495 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC33="<<fC33<<endl;
496 return 0;
497 }
498 if (fC44<=0 || fC44>6e-5) {
499 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC44="<<fC44<<endl;
500 return 0;
14825d5a 501 }
14825d5a 502 return 1;
006b5f7f 503}
504
505//____________________________________________________________________________
a9a2d814 506Int_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) {
006b5f7f 507 //------------------------------------------------------------------
508 //This function propagates a track
509 //------------------------------------------------------------------
33329428 510 Double_t alpha=fAlpha, x=fX;
006b5f7f 511 Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
512 Double_t c00=fC00;
513 Double_t c10=fC10, c11=fC11;
514 Double_t c20=fC20, c21=fC21, c22=fC22;
515 Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
516 Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
517
33329428 518 if (alp < -TMath::Pi()) alp += 2*TMath::Pi();
519 else if (alp >= TMath::Pi()) alp -= 2*TMath::Pi();
520 Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha);
521 Double_t sf=fP2, cf=TMath::Sqrt(1.- fP2*fP2);
006b5f7f 522
33329428 523 TMatrixD *T=0;
524 // **** rotation **********************
525 {
006b5f7f 526 fAlpha = alp;
33329428 527 fX = x*ca + p0*sa;
528 fP0= -x*sa + p0*ca;
529 fP2= sf*ca - cf*sa;
530
531 TMatrixD C(5,5);
532 C(0,0)=c00;
533 C(1,0)=c10; C(1,1)=c11;
534 C(2,0)=c20; C(2,1)=c21; C(2,2)=c22;
535 C(3,0)=c30; C(3,1)=c31; C(3,2)=c32; C(3,3)=c33;
536 C(4,0)=c40; C(4,1)=c41; C(4,2)=c42; C(4,3)=c43; C(4,4)=c44;
537 C(0,1)=C(1,0);
538 C(0,2)=C(2,0); C(1,2)=C(2,1);
539 C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
540 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 541
33329428 542 TMatrixD F(6,5);
543 F(0,0)=sa;
544 F(1,0)=ca;
545 F(2,1)=F(4,3)=F(5,4)=1;
546 F(3,2)=ca + sf/cf*sa;
006b5f7f 547
33329428 548 TMatrixD tmp(C,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F));
549 T=new TMatrixD(F,TMatrixD::kMult,tmp);
550 }
006b5f7f 551
33329428 552 // **** translation ******************
553 {
554 Double_t dx=xk-fX;
006b5f7f 555 Double_t f1=fP2, f2=f1 + fP4*dx;
a9a2d814 556 if (TMath::Abs(f2) >= 0.9999) {
006b5f7f 557 Int_t n=GetNumberOfClusters();
7f6ddf58 558 if (n>kWARN)
559 cerr<<n<<" AliITStrackV2::Propagate: Propagation failed !\n";
006b5f7f 560 return 0;
561 }
33329428 562 Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
006b5f7f 563
33329428 564 fX=xk;
006b5f7f 565 fP0 += dx*(f1+f2)/(r1+r2);
566 fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
567 fP2 += dx*fP4;
568
33329428 569 TMatrixD F(5,6);
570 F(0,1)=F(1,2)=F(2,3)=F(3,4)=F(4,5)=1;
571 F(0,3)=dx/(r1+r2)*(2+(f1+f2)*(f2/r2+f1/r1)/(r1+r2));
572 F(0,5)=dx*dx/(r1+r2)*(1+(f1+f2)*f2/(r1+r2));
573 F(1,3)=dx*fP3/(f1*r2 + f2*r1)*(2-(f1+f2)*(r2-f1*f2/r2+r1-f2*f1/r1)/(f1*r2 + f2*r1));
574 F(1,4)=dx*(f1+f2)/(f1*r2 + f2*r1);
575 F(1,5)=dx*dx*fP3/(f1*r2 + f2*r1)*(1-(f1+f2)*(-f1*f2/r2+r1)/(f1*r2 + f2*r1));
576 F(2,5)=dx;
577 F(0,0)=-1/(r1+r2)*((f1+f2)+dx*fP4*(1+(f1+f2)/(r1+r2)*f2/r2));
578 F(1,0)=-fP3/(f1*r2 + f2*r1)*((f1+f2)+dx*fP4*(1+(f1+f2)/(f1*r2 + f2*r1)*(f1*f2/r2-r1)));
579 F(2,0)=-fP4;
580
581 TMatrixD tmp(*T,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F));
582 delete T;
583 TMatrixD C(F,TMatrixD::kMult,tmp);
006b5f7f 584
585 fC00=C(0,0);
586 fC10=C(1,0); fC11=C(1,1);
587 fC20=C(2,0); fC21=C(2,1); fC22=C(2,2);
588 fC30=C(3,0); fC31=C(3,1); fC32=C(3,2); fC33=C(3,3);
589 fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4);
590
006b5f7f 591 if (!Invariant()) {
33329428 592 fAlpha=alpha;
593 fX=x;
006b5f7f 594 fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
595 fC00=c00;
596 fC10=c10; fC11=c11;
597 fC20=c20; fC21=c21; fC22=c22;
598 fC30=c30; fC31=c31; fC32=c32; fC33=c33;
599 fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
600 return 0;
601 }
33329428 602 }
006b5f7f 603
604 return 1;
605}
606
a6a6f799 607Double_t AliITStrackV2::GetD(Double_t x, Double_t y) const {
006b5f7f 608 //------------------------------------------------------------------
a6a6f799 609 // This function calculates the transverse impact parameter
610 // with respect to a point with global coordinates (x,y)
006b5f7f 611 //------------------------------------------------------------------
a6a6f799 612 Double_t xt=fX, yt=fP0;
613
614 Double_t sn=TMath::Sin(fAlpha), cs=TMath::Cos(fAlpha);
615 Double_t a = x*cs + y*sn;
616 y = -x*sn + y*cs; x=a;
617 xt-=x; yt-=y;
618
619 sn=fP4*xt - fP2; cs=fP4*yt + TMath::Sqrt(1.- fP2*fP2);
620 a=2*(xt*fP2 - yt*TMath::Sqrt(1.- fP2*fP2))-fP4*(xt*xt + yt*yt);
006b5f7f 621 if (fP4<0) a=-a;
622 return a/(1 + TMath::Sqrt(sn*sn + cs*cs));
a9a2d814 623}
006b5f7f 624
a9a2d814 625Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
626 //------------------------------------------------------------------
627 //This function improves angular track parameters
628 //------------------------------------------------------------------
629 Double_t dy=fP0-yv, dz=fP1-zv;
630 Double_t r2=fX*fX+dy*dy;
631 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
632 Double_t beta2=p2/(p2 + GetMass()*GetMass());
633 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
634 //Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
635 Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
636 {
637 Double_t parp=0.5*(fP4*fX + dy*TMath::Sqrt(4/r2-fP4*fP4));
638 Double_t sigma2p = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
639 sigma2p += fC00/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
1a947abf 640 sigma2p += fSigmaYV*fSigmaYV/r2;
a9a2d814 641 sigma2p += 0.25*fC44*fX*fX;
642 Double_t eps2p=sigma2p/(fC22+sigma2p);
643 fP0 += fC20/(fC22+sigma2p)*(parp-fP2);
644 fP2 = eps2p*fP2 + (1-eps2p)*parp;
645 fC22 *= eps2p;
646 fC20 *= eps2p;
647 }
648 {
649 Double_t parl=0.5*fP4*dz/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
650 Double_t sigma2l=theta2;
651 sigma2l += fC11/r2+fC00*dy*dy*dz*dz/(r2*r2*r2);
1a947abf 652 sigma2l += fSigmaZV*fSigmaZV/r2;
a9a2d814 653 Double_t eps2l=sigma2l/(fC33+sigma2l);
654 fP1 += fC31/(fC33+sigma2l)*(parl-fP3);
655 fP4 += fC43/(fC33+sigma2l)*(parl-fP3);
656 fP3 = eps2l*fP3 + (1-eps2l)*parl;
657 fC33 *= eps2l; fC43 *= eps2l;
658 fC31 *= eps2l;
659 }
660 if (!Invariant()) return 0;
661 return 1;
662}
006b5f7f 663
a9a2d814 664/*
006b5f7f 665Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
666 //------------------------------------------------------------------
667 //This function improves angular track parameters
668 //------------------------------------------------------------------
669 Double_t dy=fP0-yv, dz=fP1-zv;
670 Double_t r2=fX*fX+dy*dy;
671 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
7f6ddf58 672 Double_t beta2=p2/(p2 + GetMass()*GetMass());
006b5f7f 673 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
a9a2d814 674 //Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
675 Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
006b5f7f 676
677 Double_t par=0.5*(fP4*fX + dy*TMath::Sqrt(4/r2-fP4*fP4));
678 Double_t sigma2 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
679 sigma2 += fC00/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
680 sigma2 += kSigmaYV*kSigmaYV/r2;
681 sigma2 += 0.25*fC44*fX*fX;
682 Double_t eps2=sigma2/(fC22+sigma2), eps=TMath::Sqrt(eps2);
683 if (10*r2*fC44<fC22) {
684 fP2 = eps2*fP2 + (1-eps2)*par;
685 fC22*=eps2; fC21*=eps; fC20*=eps; fC32*=eps; fC42*=eps;
686 }
687
688 par=0.5*fP4*dz/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
689 sigma2=theta2;
690 sigma2 += fC11/r2+fC00*dy*dy*dz*dz/(r2*r2*r2);
691 sigma2 += kSigmaZV*kSigmaZV/r2;
692 eps2=sigma2/(fC33+sigma2); eps=TMath::Sqrt(eps2);
693 Double_t tgl=fP3;
694 fP3 = eps2*fP3 + (1-eps2)*par;
695 fC33*=eps2; fC32*=eps; fC31*=eps; fC30*=eps; fC43*=eps;
696
697 eps=TMath::Sqrt((1+fP3*fP3)/(1+tgl*tgl));
698 fP4*=eps;
699 fC44*=eps*eps; fC43*=eps;fC42*=eps; fC41*=eps; fC40*=eps;
700
701 if (!Invariant()) return 0;
702 return 1;
703}
006b5f7f 704*/
14825d5a 705void AliITStrackV2::ResetCovariance() {
706 //------------------------------------------------------------------
707 //This function makes a track forget its history :)
708 //------------------------------------------------------------------
709
710 fC00*=10.;
711 fC10=0.; fC11*=10.;
712 fC20=0.; fC21=0.; fC22*=10.;
713 fC30=0.; fC31=0.; fC32=0.; fC33*=10.;
714 fC40=0.; fC41=0.; fC42=0.; fC43=0.; fC44*=10.;
715
716}
23efe5f1 717
718void AliITStrackV2::CookdEdx(Double_t low, Double_t up) {
719 //-----------------------------------------------------------------
a9a2d814 720 // This function calculates dE/dX within the "low" and "up" cuts.
721 // Origin: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
23efe5f1 722 //-----------------------------------------------------------------
723 Int_t i;
a9a2d814 724 Int_t nc=4;
23efe5f1 725 // The clusters order is: SSD-2, SSD-1, SDD-2, SDD-1, SPD-2, SPD-1
726 // Take only SSD and SDD
23efe5f1 727
a9a2d814 728 Int_t swap;//stupid sorting
23efe5f1 729 do {
730 swap=0;
731 for (i=0; i<nc-1; i++) {
732 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
733 Float_t tmp=fdEdxSample[i];
734 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
735 swap++;
736 }
737 } while (swap);
738
a9a2d814 739 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc); //b.b. to take two lowest dEdX
740 // values from four ones choose
741 // nu=2
23efe5f1 742 Float_t dedx=0;
a9a2d814 743 for (i=nl; i<nu; i++) dedx += fdEdxSample[i];
744 dedx /= (nu-nl);
23efe5f1 745
23efe5f1 746 SetdEdx(dedx);
747}