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