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