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