]> 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
A simple macro to display a selected chamber a la STAR
[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;
22b13da0 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}
006b5f7f 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 //------------------------------------------------------------------
163 Double_t xv=fP2*(fX*fP2 - fP0*TMath::Sqrt(1.- fP2*fP2)); //linear approxim.
a9a2d814 164 PropagateTo(xv,d,x0);
006b5f7f 165 return 0;
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)
179 cerr<<n<<" AliITStrackV2::GetGlobalXYZat: Propagation failed !\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)
209 cerr<<n<<" AliKalmanTrack::GetPredictedChi2: Singular matrix !\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
219//_____________________________________________________________________________
220Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c,Double_t *m,
7f6ddf58 221Double_t x0) const {
006b5f7f 222 //-----------------------------------------------------------------
223 // This function calculates a chi2 increment with a vertex contraint
224 //-----------------------------------------------------------------
225 TVectorD x(5); x(0)=fP0; x(1)=fP1; x(2)=fP2; x(3)=fP3; x(4)=fP4;
226 TMatrixD C(5,5);
227 C(0,0)=fC00;
228 C(1,0)=fC10; C(1,1)=fC11;
229 C(2,0)=fC20; C(2,1)=fC21; C(2,2)=fC22;
230 C(3,0)=fC30; C(3,1)=fC31; C(3,2)=fC32; C(3,3)=fC33;
231 C(4,0)=fC40; C(4,1)=fC41; C(4,2)=fC42; C(4,3)=fC43; C(4,4)=fC44;
232
233 C(0,1)=C(1,0);
234 C(0,2)=C(2,0); C(1,2)=C(2,1);
235 C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
236 C(0,4)=C(4,0); C(1,4)=C(4,1); C(2,4)=C(4,2); C(3,4)=C(4,3);
237
238 TMatrixD H(4,5); H.UnitMatrix();
239 Double_t dy=(c->GetY() - m[0]), dz=(c->GetZ() - m[1]);
240
241 Double_t dr=TMath::Sqrt(fX*fX + dy*dy);
242 Double_t r =TMath::Sqrt(4/dr/dr - fP4*fP4);
243 Double_t sn=0.5*(fP4*fX + dy*r);
244 Double_t tg=0.5*fP4*dz/TMath::ASin(0.5*fP4*dr);
245 TVectorD mm(4);
246 mm(0)=m[0]=c->GetY(); mm(1)=m[1]=c->GetZ(); mm(2)=m[2]=sn; mm(3)=m[3]=tg;
247
248 Double_t v22=0.,v33=0.;
249 //x0=0.;
250 if (x0!=0.) {
251 Double_t pp2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
7f6ddf58 252 Double_t beta2=pp2/(pp2 + GetMass()*GetMass());
006b5f7f 253 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
254 Double_t theta2=14.1*14.1/(beta2*pp2*1e6)*x0;
255 v22 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
256 v33 = theta2*(1.+ GetTgl()*GetTgl())*(1. + GetTgl()*GetTgl());
257 }
258 Double_t sy2=c->GetSigmaY2(), sz2=c->GetSigmaZ2();
1a947abf 259 v22+=fSigmaYV*fSigmaYV/dr/dr;
006b5f7f 260 v22+=sy2/dr/dr;
261 Double_t v20=sy2/dr;
262
1a947abf 263 v33+=fSigmaZV*fSigmaZV/dr/dr;
006b5f7f 264 v33+=sz2/dr/dr;
265 Double_t v31=sz2/dr;
266
267 TMatrixD V(4,4);
268 V(0,0)=m[4 ]=sy2; V(0,1)=m[5 ]=0.; V(0,2)=m[6 ]=v20; V(0,3)=m[7 ]=0.;
269 V(1,0)=m[8 ]=0.; V(1,1)=m[9 ]=sz2; V(1,2)=m[10]=0.; V(1,3)=m[11]=v31;
270 V(2,0)=m[12]=v20; V(2,1)=m[13]=0.; V(2,2)=m[14]=v22; V(2,3)=m[15]=0.;
271 V(3,0)=m[16]=0.; V(3,1)=m[17]=v31; V(3,2)=m[18]=0.; V(3,3)=m[19]=v33;
272
273 TVectorD res=x; res*=H; res-=mm; //res*=-1;
274 TMatrixD tmp(H,TMatrixD::kMult,C);
275 TMatrixD R(tmp,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed,H)); R+=V;
276
277 Double_t det=R.Determinant();
7f6ddf58 278 if (TMath::Abs(det) < 1.e-30) {
006b5f7f 279 Int_t n=GetNumberOfClusters();
7f6ddf58 280 if (n>kWARN)
281 cerr<<n<<" AliITStrackV2::GetPredictedChi2: Singular matrix !\n";
006b5f7f 282 return 1e10;
283 }
284
285 R.Invert();
286
287 TVectorD rs=res;
288 res*=R;
289 return rs*res;
290}
291
292//____________________________________________________________________________
a9a2d814 293Int_t AliITStrackV2::CorrectForMaterial(Double_t d, Double_t x0) {
294 //------------------------------------------------------------------
295 //This function corrects the track parameters for crossed material
296 //------------------------------------------------------------------
297 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
298 Double_t beta2=p2/(p2 + GetMass()*GetMass());
299 d*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
300
301 //Multiple scattering******************
302 if (d!=0) {
33329428 303 Double_t theta2=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(d);
304 //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*TMath::Abs(d)*9.36*2.33;
a9a2d814 305 fC22 += theta2*(1.- fP2*fP2)*(1. + fP3*fP3);
306 fC33 += theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
307 fC43 += theta2*fP3*fP4*(1. + fP3*fP3);
308 fC44 += theta2*fP3*fP4*fP3*fP4;
309 }
310
311 //Energy losses************************
312 if (x0!=0.) {
313 d*=x0;
314 Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d;
315 fP4*=(1.- sqrt(p2+GetMass()*GetMass())/p2*dE);
316 }
317
318 if (!Invariant()) return 0;
319
320 return 1;
321}
322
323//____________________________________________________________________________
324Int_t AliITStrackV2::PropagateTo(Double_t xk, Double_t d, Double_t x0) {
006b5f7f 325 //------------------------------------------------------------------
326 //This function propagates a track
327 //------------------------------------------------------------------
328 Double_t x1=fX, x2=xk, dx=x2-x1;
329 Double_t f1=fP2, f2=f1 + fP4*dx;
a9a2d814 330 if (TMath::Abs(f2) >= 0.9999) {
006b5f7f 331 Int_t n=GetNumberOfClusters();
7f6ddf58 332 if (n>kWARN)
333 cerr<<n<<" AliITStrackV2::PropagateTo: Propagation failed !\n";
006b5f7f 334 return 0;
335 }
336
f29dbb4b 337 // old position [SR, GSI, 17.02.2003]
338 Double_t oldX = fX, oldY = fP0, oldZ = fP1;
339
006b5f7f 340 Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
341
342 fP0 += dx*(f1+f2)/(r1+r2);
343 fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
344 fP2 += dx*fP4;
345
346 //f = F - 1
347
348 Double_t f02= dx/(r1*r1*r1);
349 Double_t f04=0.5*dx*dx/(r1*r1*r1);
350 Double_t f12= dx*fP3*f1/(r1*r1*r1);
351 Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1);
352 Double_t f13= dx/r1;
353 Double_t f24= dx;
354
355 //b = C*ft
356 Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
357 Double_t b02=f24*fC40;
358 Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
359 Double_t b12=f24*fC41;
360 Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
361 Double_t b22=f24*fC42;
362 Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
363 Double_t b42=f24*fC44;
364 Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
365 Double_t b32=f24*fC43;
366
367 //a = f*b = f*C*ft
368 Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
369 Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
370 Double_t a22=f24*b42;
371
372 //F*C*Ft = C + (b + bt + a)
373 fC00 += b00 + b00 + a00;
374 fC10 += b10 + b01 + a01;
375 fC20 += b20 + b02 + a02;
376 fC30 += b30;
377 fC40 += b40;
378 fC11 += b11 + b11 + a11;
379 fC21 += b21 + b12 + a12;
380 fC31 += b31;
381 fC41 += b41;
382 fC22 += b22 + b22 + a22;
383 fC32 += b32;
384 fC42 += b42;
385
386 fX=x2;
387
a9a2d814 388 if (!CorrectForMaterial(d,x0)) return 0;
006b5f7f 389
f29dbb4b 390 // Integrated Time [SR, GSI, 17.02.2003]
391 if (IsStartedTimeIntegral()) {
392 Double_t l2 = (fX-oldX)*(fX-oldX)+(fP0-oldY)*(fP0-oldY)+(fP1-oldZ)*(fP1-oldZ);
393 AddTimeStep(TMath::Sqrt(l2));
394 }
395 //
396
006b5f7f 397 return 1;
398}
399
400//____________________________________________________________________________
401Int_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, UInt_t index) {
402 //------------------------------------------------------------------
403 //This function updates track parameters
404 //------------------------------------------------------------------
405 Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
406 Double_t c00=fC00;
407 Double_t c10=fC10, c11=fC11;
408 Double_t c20=fC20, c21=fC21, c22=fC22;
409 Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
410 Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
411
412
413 Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
414 r00+=fC00; r01+=fC10; r11+=fC11;
415 Double_t det=r00*r11 - r01*r01;
416 Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
417
418 Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
419 Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
420 Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
421 Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
422 Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
423
424 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
425 Double_t sf=fP2 + k20*dy + k21*dz;
7f6ddf58 426
006b5f7f 427 fP0 += k00*dy + k01*dz;
428 fP1 += k10*dy + k11*dz;
429 fP2 = sf;
430 fP3 += k30*dy + k31*dz;
431 fP4 += k40*dy + k41*dz;
432
433 Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
434 Double_t c12=fC21, c13=fC31, c14=fC41;
435
436 fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
437 fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
438 fC40-=k00*c04+k01*c14;
439
440 fC11-=k10*c01+k11*fC11;
441 fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
442 fC41-=k10*c04+k11*c14;
443
444 fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
445 fC42-=k20*c04+k21*c14;
446
447 fC33-=k30*c03+k31*c13;
448 fC43-=k30*c04+k31*c14;
449
450 fC44-=k40*c04+k41*c14;
451
452 if (!Invariant()) {
453 fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
454 fC00=c00;
455 fC10=c10; fC11=c11;
456 fC20=c20; fC21=c21; fC22=c22;
457 fC30=c30; fC31=c31; fC32=c32; fC33=c33;
458 fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
459 return 0;
460 }
461
462 Int_t n=GetNumberOfClusters();
463 fIndex[n]=index;
464 SetNumberOfClusters(n+1);
465 SetChi2(GetChi2()+chi2);
466
467 return 1;
468}
469
006b5f7f 470Int_t AliITStrackV2::Invariant() const {
471 //------------------------------------------------------------------
472 // This function is for debugging purpose only
473 //------------------------------------------------------------------
7f6ddf58 474 Int_t n=GetNumberOfClusters();
475
a9a2d814 476 if (TMath::Abs(fP2)>=0.9999){
477 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fP2="<<fP2<<endl;
478 return 0;
479 }
480 if (fC00<=0 || fC00>9.) {
481 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC00="<<fC00<<endl;
482 return 0;
483 }
484 if (fC11<=0 || fC11>9.) {
485 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC11="<<fC11<<endl;
486 return 0;
487 }
488 if (fC22<=0 || fC22>1.) {
489 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC22="<<fC22<<endl;
490 return 0;
491 }
492 if (fC33<=0 || fC33>1.) {
493 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC33="<<fC33<<endl;
494 return 0;
495 }
496 if (fC44<=0 || fC44>6e-5) {
497 if (n>kWARN) cout<<"AliITStrackV2::Invariant : fC44="<<fC44<<endl;
498 return 0;
14825d5a 499 }
14825d5a 500 return 1;
006b5f7f 501}
502
503//____________________________________________________________________________
a9a2d814 504Int_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) {
006b5f7f 505 //------------------------------------------------------------------
506 //This function propagates a track
507 //------------------------------------------------------------------
33329428 508 Double_t alpha=fAlpha, x=fX;
006b5f7f 509 Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
510 Double_t c00=fC00;
511 Double_t c10=fC10, c11=fC11;
512 Double_t c20=fC20, c21=fC21, c22=fC22;
513 Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
514 Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
515
33329428 516 if (alp < -TMath::Pi()) alp += 2*TMath::Pi();
517 else if (alp >= TMath::Pi()) alp -= 2*TMath::Pi();
518 Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha);
519 Double_t sf=fP2, cf=TMath::Sqrt(1.- fP2*fP2);
006b5f7f 520
33329428 521 TMatrixD *T=0;
522 // **** rotation **********************
523 {
006b5f7f 524 fAlpha = alp;
33329428 525 fX = x*ca + p0*sa;
526 fP0= -x*sa + p0*ca;
527 fP2= sf*ca - cf*sa;
528
529 TMatrixD C(5,5);
530 C(0,0)=c00;
531 C(1,0)=c10; C(1,1)=c11;
532 C(2,0)=c20; C(2,1)=c21; C(2,2)=c22;
533 C(3,0)=c30; C(3,1)=c31; C(3,2)=c32; C(3,3)=c33;
534 C(4,0)=c40; C(4,1)=c41; C(4,2)=c42; C(4,3)=c43; C(4,4)=c44;
535 C(0,1)=C(1,0);
536 C(0,2)=C(2,0); C(1,2)=C(2,1);
537 C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
538 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 539
33329428 540 TMatrixD F(6,5);
541 F(0,0)=sa;
542 F(1,0)=ca;
543 F(2,1)=F(4,3)=F(5,4)=1;
544 F(3,2)=ca + sf/cf*sa;
006b5f7f 545
33329428 546 TMatrixD tmp(C,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F));
547 T=new TMatrixD(F,TMatrixD::kMult,tmp);
548 }
006b5f7f 549
33329428 550 // **** translation ******************
551 {
552 Double_t dx=xk-fX;
006b5f7f 553 Double_t f1=fP2, f2=f1 + fP4*dx;
a9a2d814 554 if (TMath::Abs(f2) >= 0.9999) {
006b5f7f 555 Int_t n=GetNumberOfClusters();
7f6ddf58 556 if (n>kWARN)
557 cerr<<n<<" AliITStrackV2::Propagate: Propagation failed !\n";
006b5f7f 558 return 0;
559 }
33329428 560 Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
006b5f7f 561
33329428 562 fX=xk;
006b5f7f 563 fP0 += dx*(f1+f2)/(r1+r2);
564 fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
565 fP2 += dx*fP4;
566
33329428 567 TMatrixD F(5,6);
568 F(0,1)=F(1,2)=F(2,3)=F(3,4)=F(4,5)=1;
569 F(0,3)=dx/(r1+r2)*(2+(f1+f2)*(f2/r2+f1/r1)/(r1+r2));
570 F(0,5)=dx*dx/(r1+r2)*(1+(f1+f2)*f2/(r1+r2));
571 F(1,3)=dx*fP3/(f1*r2 + f2*r1)*(2-(f1+f2)*(r2-f1*f2/r2+r1-f2*f1/r1)/(f1*r2 + f2*r1));
572 F(1,4)=dx*(f1+f2)/(f1*r2 + f2*r1);
573 F(1,5)=dx*dx*fP3/(f1*r2 + f2*r1)*(1-(f1+f2)*(-f1*f2/r2+r1)/(f1*r2 + f2*r1));
574 F(2,5)=dx;
575 F(0,0)=-1/(r1+r2)*((f1+f2)+dx*fP4*(1+(f1+f2)/(r1+r2)*f2/r2));
576 F(1,0)=-fP3/(f1*r2 + f2*r1)*((f1+f2)+dx*fP4*(1+(f1+f2)/(f1*r2 + f2*r1)*(f1*f2/r2-r1)));
577 F(2,0)=-fP4;
578
579 TMatrixD tmp(*T,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F));
580 delete T;
581 TMatrixD C(F,TMatrixD::kMult,tmp);
006b5f7f 582
583 fC00=C(0,0);
584 fC10=C(1,0); fC11=C(1,1);
585 fC20=C(2,0); fC21=C(2,1); fC22=C(2,2);
586 fC30=C(3,0); fC31=C(3,1); fC32=C(3,2); fC33=C(3,3);
587 fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4);
588
006b5f7f 589 if (!Invariant()) {
33329428 590 fAlpha=alpha;
591 fX=x;
006b5f7f 592 fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
593 fC00=c00;
594 fC10=c10; fC11=c11;
595 fC20=c20; fC21=c21; fC22=c22;
596 fC30=c30; fC31=c31; fC32=c32; fC33=c33;
597 fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
598 return 0;
599 }
33329428 600 }
006b5f7f 601
602 return 1;
603}
604
a6a6f799 605Double_t AliITStrackV2::GetD(Double_t x, Double_t y) const {
006b5f7f 606 //------------------------------------------------------------------
a6a6f799 607 // This function calculates the transverse impact parameter
608 // with respect to a point with global coordinates (x,y)
006b5f7f 609 //------------------------------------------------------------------
a6a6f799 610 Double_t xt=fX, yt=fP0;
611
612 Double_t sn=TMath::Sin(fAlpha), cs=TMath::Cos(fAlpha);
613 Double_t a = x*cs + y*sn;
614 y = -x*sn + y*cs; x=a;
615 xt-=x; yt-=y;
616
617 sn=fP4*xt - fP2; cs=fP4*yt + TMath::Sqrt(1.- fP2*fP2);
618 a=2*(xt*fP2 - yt*TMath::Sqrt(1.- fP2*fP2))-fP4*(xt*xt + yt*yt);
006b5f7f 619 if (fP4<0) a=-a;
620 return a/(1 + TMath::Sqrt(sn*sn + cs*cs));
a9a2d814 621}
006b5f7f 622
a9a2d814 623Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
624 //------------------------------------------------------------------
625 //This function improves angular track parameters
626 //------------------------------------------------------------------
627 Double_t dy=fP0-yv, dz=fP1-zv;
628 Double_t r2=fX*fX+dy*dy;
629 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
630 Double_t beta2=p2/(p2 + GetMass()*GetMass());
631 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
632 //Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
633 Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
634 {
635 Double_t parp=0.5*(fP4*fX + dy*TMath::Sqrt(4/r2-fP4*fP4));
636 Double_t sigma2p = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
637 sigma2p += fC00/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
1a947abf 638 sigma2p += fSigmaYV*fSigmaYV/r2;
a9a2d814 639 sigma2p += 0.25*fC44*fX*fX;
640 Double_t eps2p=sigma2p/(fC22+sigma2p);
641 fP0 += fC20/(fC22+sigma2p)*(parp-fP2);
642 fP2 = eps2p*fP2 + (1-eps2p)*parp;
643 fC22 *= eps2p;
644 fC20 *= eps2p;
645 }
646 {
647 Double_t parl=0.5*fP4*dz/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
648 Double_t sigma2l=theta2;
649 sigma2l += fC11/r2+fC00*dy*dy*dz*dz/(r2*r2*r2);
1a947abf 650 sigma2l += fSigmaZV*fSigmaZV/r2;
a9a2d814 651 Double_t eps2l=sigma2l/(fC33+sigma2l);
652 fP1 += fC31/(fC33+sigma2l)*(parl-fP3);
653 fP4 += fC43/(fC33+sigma2l)*(parl-fP3);
654 fP3 = eps2l*fP3 + (1-eps2l)*parl;
655 fC33 *= eps2l; fC43 *= eps2l;
656 fC31 *= eps2l;
657 }
658 if (!Invariant()) return 0;
659 return 1;
660}
006b5f7f 661
a9a2d814 662/*
006b5f7f 663Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
664 //------------------------------------------------------------------
665 //This function improves angular track parameters
666 //------------------------------------------------------------------
667 Double_t dy=fP0-yv, dz=fP1-zv;
668 Double_t r2=fX*fX+dy*dy;
669 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
7f6ddf58 670 Double_t beta2=p2/(p2 + GetMass()*GetMass());
006b5f7f 671 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
a9a2d814 672 //Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
673 Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
006b5f7f 674
675 Double_t par=0.5*(fP4*fX + dy*TMath::Sqrt(4/r2-fP4*fP4));
676 Double_t sigma2 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
677 sigma2 += fC00/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
678 sigma2 += kSigmaYV*kSigmaYV/r2;
679 sigma2 += 0.25*fC44*fX*fX;
680 Double_t eps2=sigma2/(fC22+sigma2), eps=TMath::Sqrt(eps2);
681 if (10*r2*fC44<fC22) {
682 fP2 = eps2*fP2 + (1-eps2)*par;
683 fC22*=eps2; fC21*=eps; fC20*=eps; fC32*=eps; fC42*=eps;
684 }
685
686 par=0.5*fP4*dz/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
687 sigma2=theta2;
688 sigma2 += fC11/r2+fC00*dy*dy*dz*dz/(r2*r2*r2);
689 sigma2 += kSigmaZV*kSigmaZV/r2;
690 eps2=sigma2/(fC33+sigma2); eps=TMath::Sqrt(eps2);
691 Double_t tgl=fP3;
692 fP3 = eps2*fP3 + (1-eps2)*par;
693 fC33*=eps2; fC32*=eps; fC31*=eps; fC30*=eps; fC43*=eps;
694
695 eps=TMath::Sqrt((1+fP3*fP3)/(1+tgl*tgl));
696 fP4*=eps;
697 fC44*=eps*eps; fC43*=eps;fC42*=eps; fC41*=eps; fC40*=eps;
698
699 if (!Invariant()) return 0;
700 return 1;
701}
006b5f7f 702*/
14825d5a 703void AliITStrackV2::ResetCovariance() {
704 //------------------------------------------------------------------
705 //This function makes a track forget its history :)
706 //------------------------------------------------------------------
707
708 fC00*=10.;
709 fC10=0.; fC11*=10.;
710 fC20=0.; fC21=0.; fC22*=10.;
711 fC30=0.; fC31=0.; fC32=0.; fC33*=10.;
712 fC40=0.; fC41=0.; fC42=0.; fC43=0.; fC44*=10.;
713
714}
23efe5f1 715
716void AliITStrackV2::CookdEdx(Double_t low, Double_t up) {
717 //-----------------------------------------------------------------
a9a2d814 718 // This function calculates dE/dX within the "low" and "up" cuts.
719 // Origin: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
23efe5f1 720 //-----------------------------------------------------------------
721 Int_t i;
a9a2d814 722 Int_t nc=4;
23efe5f1 723 // The clusters order is: SSD-2, SSD-1, SDD-2, SDD-1, SPD-2, SPD-1
724 // Take only SSD and SDD
23efe5f1 725
a9a2d814 726 Int_t swap;//stupid sorting
23efe5f1 727 do {
728 swap=0;
729 for (i=0; i<nc-1; i++) {
730 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
731 Float_t tmp=fdEdxSample[i];
732 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
733 swap++;
734 }
735 } while (swap);
736
a9a2d814 737 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc); //b.b. to take two lowest dEdX
738 // values from four ones choose
739 // nu=2
23efe5f1 740 Float_t dedx=0;
a9a2d814 741 for (i=nl; i<nu; i++) dedx += fdEdxSample[i];
742 dedx /= (nu-nl);
23efe5f1 743
23efe5f1 744 SetdEdx(dedx);
745}