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