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