]> 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
Rearranged conditions for faster execution
[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"
e43c066c 28#include "AliTPCtrack.h"
83d73500 29#include "AliESDtrack.h"
006b5f7f 30#include "AliITStrackV2.h"
31
1aedd96e 32ClassImp(AliITStrackV2)
8676d691 33
e43c066c 34const Int_t kWARN=5;
35
22b13da0 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),
babd135a 61 fNUsed(0),
62 fNSkipped(0),
e43c066c 63 fNDeadZone(0),
64 fDeadZoneProbability(0),
65 fReconstructed(kFALSE),
66 fConstrain(kFALSE),
83d73500 67 fESDtrack(0)
e43c066c 68{
69 for(Int_t i=0; i<kMaxLayer; i++) {fIndex[i]=0;fClIndex[i]=-1;}
70 for(Int_t i=0; i<4; i++) fdEdxSample[i]=0;
71 for(Int_t i=0; i<6; i++) { fNy[i]=0; fNz[i]=0; fNormQ[i]=0; fNormChi2[i]=1000;}
72 for(Int_t i=0; i<12; i++) {fDy[i]=0; fDz[i]=0; fSigmaY[i]=0; fSigmaZ[i]=0; fChi2MIP[i]=0;}
73 fD[0]=0; fD[1]=0;
74 fExpQ=40;
75 fdEdxMismatch=0;
76 fChi22=0;
77}
78
79//____________________________________________________________________________
80AliITStrackV2::AliITStrackV2(const AliTPCtrack& t) throw (const Char_t *) :
81AliKalmanTrack(t) {
1aedd96e 82 //------------------------------------------------------------------
e43c066c 83 //Conversion TPC track -> ITS track
1aedd96e 84 //------------------------------------------------------------------
e43c066c 85 SetChi2(0.);
86 SetNumberOfClusters(0);
87
88 fdEdx = t.GetdEdx();
89 SetMass(t.GetMass());
90
91 fAlpha = t.GetAlpha();
92 if (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
93 else if (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
94
95 //Conversion of the track parameters
96 Double_t x,p[5]; t.GetExternalParameters(x,p);
97 fX=x; x=GetConvConst();
98 fP0=p[0];
99 fP1=p[1];
100 fP2=p[2];
101 fP3=p[3];
102 fP4=p[4]/x;
103
104 //Conversion of the covariance matrix
105 Double_t c[15]; t.GetExternalCovariance(c);
106
107 fC00=c[0 ];
108 fC10=c[1 ]; fC11=c[2 ];
109 fC20=c[3 ]; fC21=c[4 ]; fC22=c[5 ];
110 fC30=c[6 ]; fC31=c[7 ]; fC32=c[8 ]; fC33=c[9 ];
111 fC40=c[10]/x; fC41=c[11]/x; fC42=c[12]/x; fC43=c[13]/x; fC44=c[14]/x/x;
112
113 for(Int_t i=0; i<6; i++) { fNy[i]=0; fNz[i]=0; fNormQ[i]=0; fNormChi2[i]=1000;}
114 for(Int_t i=0; i<12; i++) {fDy[i]=0; fDz[i]=0; fSigmaY[i]=0; fSigmaZ[i]=0; }
115 fConstrain=kFALSE;
116 //
117 if (!Invariant()) throw "AliITStrackV2: conversion failed !\n";
118
22b13da0 119}
83d73500 120
83d73500 121//____________________________________________________________________________
67c3dcbe 122AliITStrackV2::AliITStrackV2(AliESDtrack& t,Bool_t c) throw (const Char_t *) :
83d73500 123AliKalmanTrack() {
124 //------------------------------------------------------------------
67c3dcbe 125 // Conversion ESD track -> ITS track.
126 // If c==kTRUE, create the ITS track out of the constrained params.
83d73500 127 //------------------------------------------------------------------
128 SetNumberOfClusters(t.GetITSclusters(fIndex));
129 SetLabel(t.GetLabel());
130 SetMass(t.GetMass());
e43c066c 131 //
132 //
83d73500 133
134 fdEdx=t.GetITSsignal();
135 fAlpha = t.GetAlpha();
136 if (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
137 else if (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
138
9b280d80 139 //Conversion of the track parameters
67c3dcbe 140 Double_t x,p[5];
141 if (c) t.GetConstrainedExternalParameters(x,p);
142 else t.GetExternalParameters(x,p);
9b280d80 143 fX=x; x=GetConvConst();
006b5f7f 144 fP0=p[0];
145 fP1=p[1];
146 fP2=p[2];
147 fP3=p[3];
148 fP4=p[4]/x;
149
9b280d80 150 //Conversion of the covariance matrix
67c3dcbe 151 Double_t cv[15];
152 if (c) t.GetConstrainedExternalCovariance(cv);
153 else t.GetExternalCovariance(cv);
154 fC00=cv[0 ];
155 fC10=cv[1 ]; fC11=cv[2 ];
156 fC20=cv[3 ]; fC21=cv[4 ]; fC22=cv[5 ];
157 fC30=cv[6 ]; fC31=cv[7 ]; fC32=cv[8 ]; fC33=cv[9 ];
158 fC40=cv[10]/x; fC41=cv[11]/x; fC42=cv[12]/x; fC43=cv[13]/x; fC44=cv[14]/x/x;
006b5f7f 159
83d73500 160 if (t.GetStatus()&AliESDtrack::kTIME) {
161 StartTimeIntegral();
162 Double_t times[10]; t.GetIntegratedTimes(times); SetIntegratedTimes(times);
163 SetIntegratedLength(t.GetIntegratedLength());
164 }
165 fESDtrack=&t;
babd135a 166 fNUsed = 0;
167 fReconstructed = kFALSE;
e43c066c 168 fNSkipped =0;
169 fNDeadZone = 0;
170 fDeadZoneProbability = 0;
171 for(Int_t i=0; i<6; i++) {fClIndex[i]=-1; fNy[i]=0; fNz[i]=0; fNormQ[i]=0; fNormChi2[i]=1000;}
172 for(Int_t i=0; i<12; i++) {fDy[i]=0; fDz[i]=0; fSigmaY[i]=0; fSigmaZ[i]=0;fChi2MIP[i]=0;}
173 fD[0]=0; fD[1]=0;
174 fExpQ=40;
175 fConstrain = kFALSE;
176 fdEdxMismatch=0;
177 fChi22 =0;
178
0aa211b1 179 //if (!Invariant()) throw "AliITStrackV2: conversion failed !\n";
e43c066c 180
006b5f7f 181}
182
83d73500 183void AliITStrackV2::UpdateESDtrack(ULong_t flags) {
184 fESDtrack->UpdateTrackParams(this,flags);
babd135a 185 if (flags == AliESDtrack::kITSin) fESDtrack->SetITSChi2MIP(fChi2MIP);
83d73500 186}
67c3dcbe 187void AliITStrackV2::SetConstrainedESDtrack(Double_t chi2) {
188 fESDtrack->SetConstrainedTrackParams(this,chi2);
189}
83d73500 190
006b5f7f 191//____________________________________________________________________________
192AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) : AliKalmanTrack(t) {
193 //------------------------------------------------------------------
194 //Copy constructor
195 //------------------------------------------------------------------
196 fX=t.fX;
197 fAlpha=t.fAlpha;
198 fdEdx=t.fdEdx;
199
200 fP0=t.fP0; fP1=t.fP1; fP2=t.fP2; fP3=t.fP3; fP4=t.fP4;
201
202 fC00=t.fC00;
203 fC10=t.fC10; fC11=t.fC11;
204 fC20=t.fC20; fC21=t.fC21; fC22=t.fC22;
205 fC30=t.fC30; fC31=t.fC31; fC32=t.fC32; fC33=t.fC33;
206 fC40=t.fC40; fC41=t.fC41; fC42=t.fC42; fC43=t.fC43; fC44=t.fC44;
207
208 Int_t n=GetNumberOfClusters();
23efe5f1 209 for (Int_t i=0; i<n; i++) {
e43c066c 210 fIndex[i]=t.fIndex[i];
211 if (i<4) fdEdxSample[i]=t.fdEdxSample[i];
23efe5f1 212 }
83d73500 213 fESDtrack=t.fESDtrack;
babd135a 214 fNUsed = t.fNUsed;
215 fReconstructed = t.fReconstructed;
216 fNSkipped = t.fNSkipped;
e43c066c 217 fNDeadZone = t.fNDeadZone;
218 fDeadZoneProbability = t.fDeadZoneProbability;
219 fLab = t.fLab;
220 fFakeRatio = t.fFakeRatio;
221 fdEdxMismatch = t.fdEdxMismatch;
222 fChi22 = t.fChi22;
223
224
225 fD[0]=t.fD[0]; fD[1]=t.fD[1];
226 fExpQ= t.fExpQ;
227 for(Int_t i=0; i<6; i++) {
228 fClIndex[i]= t.fClIndex[i]; fNy[i]=t.fNy[i]; fNz[i]=t.fNz[i]; fNormQ[i]=t.fNormQ[i]; fNormChi2[i] = t.fNormChi2[i];
229 }
230 for(Int_t i=0; i<12; i++) {fDy[i]=t.fDy[i]; fDz[i]=t.fDz[i];
231 fSigmaY[i]=t.fSigmaY[i]; fSigmaZ[i]=t.fSigmaZ[i];fChi2MIP[i]=t.fChi2MIP[i];}
232 fConstrain = t.fConstrain;
233 //memcpy(fDy,t.fDy,6*sizeof(Float_t));
234 //memcpy(fDz,t.fDz,6*sizeof(Float_t));
235 //memcpy(fSigmaY,t.fSigmaY,6*sizeof(Float_t));
236 //memcpy(fSigmaZ,t.fSigmaZ,6*sizeof(Float_t));
237 //memcpy(fChi2MIP,t.fChi2MIP,12*sizeof(Float_t));
006b5f7f 238}
a9a2d814 239
006b5f7f 240//_____________________________________________________________________________
241Int_t AliITStrackV2::Compare(const TObject *o) const {
242 //-----------------------------------------------------------------
243 // This function compares tracks according to the their curvature
244 //-----------------------------------------------------------------
7f6ddf58 245 AliITStrackV2 *t=(AliITStrackV2*)o;
a9a2d814 246 //Double_t co=TMath::Abs(t->Get1Pt());
247 //Double_t c =TMath::Abs(Get1Pt());
e43c066c 248 Double_t co=t->GetSigmaY2()*t->GetSigmaZ2()*(0.5+TMath::Sqrt(0.5*t->fD[0]*t->fD[0]+t->fD[1]*t->fD[1]));
249 Double_t c =GetSigmaY2()*GetSigmaZ2()*(0.5+TMath::Sqrt(0.5*fD[0]*fD[0]+fD[1]*fD[1]));
006b5f7f 250 if (c>co) return 1;
251 else if (c<co) return -1;
252 return 0;
253}
254
255//_____________________________________________________________________________
256void AliITStrackV2::GetExternalCovariance(Double_t cc[15]) const {
257 //-------------------------------------------------------------------------
258 // This function returns an external representation of the covriance matrix.
259 // (See comments in AliTPCtrack.h about external track representation)
260 //-------------------------------------------------------------------------
9b280d80 261 Double_t a=GetConvConst();
006b5f7f 262
263 cc[0 ]=fC00;
264 cc[1 ]=fC10; cc[2 ]=fC11;
265 cc[3 ]=fC20; cc[4 ]=fC21; cc[5 ]=fC22;
266 cc[6 ]=fC30; cc[7 ]=fC31; cc[8 ]=fC32; cc[9 ]=fC33;
267 cc[10]=fC40*a; cc[11]=fC41*a; cc[12]=fC42*a; cc[13]=fC43*a; cc[14]=fC44*a*a;
268}
269
270//____________________________________________________________________________
a9a2d814 271Int_t AliITStrackV2::PropagateToVertex(Double_t d,Double_t x0) {
006b5f7f 272 //------------------------------------------------------------------
273 //This function propagates a track to the minimal distance from the origin
274 //------------------------------------------------------------------
5b316485 275 //Double_t xv=fP2*(fX*fP2 - fP0*TMath::Sqrt(1.- fP2*fP2)); //linear approxim.
276 Double_t tgf=-(fP4*fX - fP2)/(fP4*fP0 + TMath::Sqrt(1 - fP2*fP2));
277 Double_t snf=tgf/TMath::Sqrt(1.+ tgf*tgf);
278 Double_t xv=(snf - fP2)/fP4 + fX;
279 return PropagateTo(xv,d,x0);
006b5f7f 280}
281
282//____________________________________________________________________________
283Int_t AliITStrackV2::
284GetGlobalXYZat(Double_t xk, Double_t &x, Double_t &y, Double_t &z) const {
285 //------------------------------------------------------------------
286 //This function returns a track position in the global system
287 //------------------------------------------------------------------
288 Double_t dx=xk-fX;
289 Double_t f1=fP2, f2=f1 + fP4*dx;
a9a2d814 290 if (TMath::Abs(f2) >= 0.9999) {
7f6ddf58 291 Int_t n=GetNumberOfClusters();
292 if (n>kWARN)
8676d691 293 Warning("GetGlobalXYZat","Propagation failed (%d) !\n",n);
006b5f7f 294 return 0;
295 }
296
297 Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
298
299 Double_t yk = fP0 + dx*(f1+f2)/(r1+r2);
300 Double_t zk = fP1 + dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
301
302 Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
303 x = xk*cs - yk*sn;
304 y = xk*sn + yk*cs;
305 z = zk;
306
307 return 1;
308}
309
310//_____________________________________________________________________________
311Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c) const
312{
313 //-----------------------------------------------------------------
314 // This function calculates a predicted chi2 increment.
315 //-----------------------------------------------------------------
316 Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
317 r00+=fC00; r01+=fC10; r11+=fC11;
babd135a 318 //
006b5f7f 319 Double_t det=r00*r11 - r01*r01;
7f6ddf58 320 if (TMath::Abs(det) < 1.e-30) {
006b5f7f 321 Int_t n=GetNumberOfClusters();
7f6ddf58 322 if (n>kWARN)
8676d691 323 Warning("GetPredictedChi2","Singular matrix (%d) !\n",n);
006b5f7f 324 return 1e10;
325 }
326 Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
33329428 327
006b5f7f 328 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
33329428 329
006b5f7f 330 return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
331}
332
e43c066c 333Double_t AliITStrackV2::GetPredictedChi2MI(Double_t cy, Double_t cz, Double_t cerry, Double_t cerrz) const
334{
335 //-----------------------------------------------------------------
336 // This function calculates a predicted chi2 increment.
337 //-----------------------------------------------------------------
338 Double_t r00=cerry*cerry, r01=0., r11=cerrz*cerrz;
339 r00+=fC00; r01+=fC10; r11+=fC11;
340 //
341 Double_t det=r00*r11 - r01*r01;
342 if (TMath::Abs(det) < 1.e-30) {
343 Int_t n=GetNumberOfClusters();
344 if (n>kWARN)
345 Warning("GetPredictedChi2","Singular matrix (%d) !\n",n);
346 return 1e10;
347 }
348 Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
349
350 Double_t dy=cy - fP0, dz=cz - fP1;
351
352 return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
353}
354
006b5f7f 355//____________________________________________________________________________
a9a2d814 356Int_t AliITStrackV2::CorrectForMaterial(Double_t d, Double_t x0) {
357 //------------------------------------------------------------------
358 //This function corrects the track parameters for crossed material
359 //------------------------------------------------------------------
e43c066c 360 // Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
361 Double_t p2=(1.+ fP3*fP3)/(Get1Pt()*Get1Pt());
362 Double_t et = p2 + GetMass()*GetMass();
363 Double_t beta2=p2/et;
364 et = sqrt(et);
a9a2d814 365 d*=TMath::Sqrt((1.+ fP3*fP3)/(1.- fP2*fP2));
e43c066c 366 //d*=TMath::Sqrt(1.+ fP3*fP3 +fP2*fP2/(1.- fP2*fP2));
a9a2d814 367
368 //Multiple scattering******************
369 if (d!=0) {
33329428 370 Double_t theta2=14.1*14.1/(beta2*p2*1e6)*TMath::Abs(d);
371 //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*TMath::Abs(d)*9.36*2.33;
a9a2d814 372 fC22 += theta2*(1.- fP2*fP2)*(1. + fP3*fP3);
373 fC33 += theta2*(1. + fP3*fP3)*(1. + fP3*fP3);
374 fC43 += theta2*fP3*fP4*(1. + fP3*fP3);
375 fC44 += theta2*fP3*fP4*fP3*fP4;
376 }
377
378 //Energy losses************************
379 if (x0!=0.) {
380 d*=x0;
e43c066c 381 // Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d;
382 //Double_t dE=0;
383 Double_t dE = 0.265*0.153e-3*(39.2-55.6*beta2+28.7*beta2*beta2+27.41/beta2)*d;
384 /*
385 if (beta2/(1-beta2)>3.5*3.5){
83d73500 386 dE=0.153e-3/beta2*(log(3.5*5940)+0.5*log(beta2/(1-beta2)) - beta2)*d;
e43c066c 387 }
388 else{
389 dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d;
390 dE+=0.06e-3/(beta2*beta2)*d;
391 }
392 */
393 fP4*=(1.- et/p2*dE);
394 Double_t delta44 = (dE*fP4*et/p2);
395 delta44*=delta44;
396 fC44+= delta44/400.;
a9a2d814 397 }
398
399 if (!Invariant()) return 0;
400
401 return 1;
402}
403
404//____________________________________________________________________________
405Int_t AliITStrackV2::PropagateTo(Double_t xk, Double_t d, Double_t x0) {
006b5f7f 406 //------------------------------------------------------------------
407 //This function propagates a track
408 //------------------------------------------------------------------
409 Double_t x1=fX, x2=xk, dx=x2-x1;
410 Double_t f1=fP2, f2=f1 + fP4*dx;
a9a2d814 411 if (TMath::Abs(f2) >= 0.9999) {
006b5f7f 412 Int_t n=GetNumberOfClusters();
7f6ddf58 413 if (n>kWARN)
8676d691 414 Warning("PropagateTo","Propagation failed !\n",n);
006b5f7f 415 return 0;
416 }
417
f29dbb4b 418 // old position [SR, GSI, 17.02.2003]
419 Double_t oldX = fX, oldY = fP0, oldZ = fP1;
420
006b5f7f 421 Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
422
423 fP0 += dx*(f1+f2)/(r1+r2);
424 fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
425 fP2 += dx*fP4;
426
427 //f = F - 1
428
429 Double_t f02= dx/(r1*r1*r1);
430 Double_t f04=0.5*dx*dx/(r1*r1*r1);
431 Double_t f12= dx*fP3*f1/(r1*r1*r1);
432 Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1);
433 Double_t f13= dx/r1;
434 Double_t f24= dx;
435
436 //b = C*ft
437 Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
438 Double_t b02=f24*fC40;
439 Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
440 Double_t b12=f24*fC41;
441 Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
442 Double_t b22=f24*fC42;
443 Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
444 Double_t b42=f24*fC44;
445 Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
446 Double_t b32=f24*fC43;
447
448 //a = f*b = f*C*ft
449 Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
450 Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
451 Double_t a22=f24*b42;
452
453 //F*C*Ft = C + (b + bt + a)
454 fC00 += b00 + b00 + a00;
455 fC10 += b10 + b01 + a01;
456 fC20 += b20 + b02 + a02;
457 fC30 += b30;
458 fC40 += b40;
459 fC11 += b11 + b11 + a11;
460 fC21 += b21 + b12 + a12;
461 fC31 += b31;
462 fC41 += b41;
463 fC22 += b22 + b22 + a22;
464 fC32 += b32;
465 fC42 += b42;
466
467 fX=x2;
468
a9a2d814 469 if (!CorrectForMaterial(d,x0)) return 0;
006b5f7f 470
f29dbb4b 471 // Integrated Time [SR, GSI, 17.02.2003]
880f41b9 472 if (IsStartedTimeIntegral() && fX>oldX) {
8676d691 473 Double_t l2 = (fX-oldX)*(fX-oldX)+(fP0-oldY)*(fP0-oldY)+
474 (fP1-oldZ)*(fP1-oldZ);
f29dbb4b 475 AddTimeStep(TMath::Sqrt(l2));
476 }
477 //
478
006b5f7f 479 return 1;
480}
481
482//____________________________________________________________________________
483Int_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, UInt_t index) {
484 //------------------------------------------------------------------
485 //This function updates track parameters
486 //------------------------------------------------------------------
487 Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
488 Double_t c00=fC00;
489 Double_t c10=fC10, c11=fC11;
490 Double_t c20=fC20, c21=fC21, c22=fC22;
491 Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
492 Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
493
494
495 Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
496 r00+=fC00; r01+=fC10; r11+=fC11;
497 Double_t det=r00*r11 - r01*r01;
498 Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
499
babd135a 500
006b5f7f 501 Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
502 Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
503 Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
504 Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
505 Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
506
507 Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
babd135a 508 Int_t layer = (index & 0xf0000000) >> 28;
509 fDy[layer] = dy;
510 fDz[layer] = dz;
511 fSigmaY[layer] = TMath::Sqrt(c->GetSigmaY2()+fC00);
512 fSigmaZ[layer] = TMath::Sqrt(c->GetSigmaZ2()+fC11);
513
006b5f7f 514 Double_t sf=fP2 + k20*dy + k21*dz;
7f6ddf58 515
006b5f7f 516 fP0 += k00*dy + k01*dz;
517 fP1 += k10*dy + k11*dz;
518 fP2 = sf;
519 fP3 += k30*dy + k31*dz;
520 fP4 += k40*dy + k41*dz;
521
522 Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
523 Double_t c12=fC21, c13=fC31, c14=fC41;
524
525 fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
526 fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
527 fC40-=k00*c04+k01*c14;
528
529 fC11-=k10*c01+k11*fC11;
530 fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
531 fC41-=k10*c04+k11*c14;
532
533 fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
534 fC42-=k20*c04+k21*c14;
535
536 fC33-=k30*c03+k31*c13;
537 fC43-=k30*c04+k31*c14;
538
539 fC44-=k40*c04+k41*c14;
540
541 if (!Invariant()) {
542 fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
543 fC00=c00;
544 fC10=c10; fC11=c11;
545 fC20=c20; fC21=c21; fC22=c22;
546 fC30=c30; fC31=c31; fC32=c32; fC33=c33;
547 fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
548 return 0;
549 }
550
67c3dcbe 551 if (chi2<0) return 1;
552
006b5f7f 553 Int_t n=GetNumberOfClusters();
554 fIndex[n]=index;
555 SetNumberOfClusters(n+1);
556 SetChi2(GetChi2()+chi2);
557
558 return 1;
559}
560
e43c066c 561//____________________________________________________________________________
562Int_t AliITStrackV2::UpdateMI(Double_t cy, Double_t cz, Double_t cerry, Double_t cerrz, Double_t chi2,UInt_t index) {
563 //------------------------------------------------------------------
564 //This function updates track parameters
565 //------------------------------------------------------------------
566 Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
567 Double_t c00=fC00;
568 Double_t c10=fC10, c11=fC11;
569 Double_t c20=fC20, c21=fC21, c22=fC22;
570 Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
571 Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
572
573
574 Double_t r00=cerry*cerry, r01=0., r11=cerrz*cerrz;
575 r00+=fC00; r01+=fC10; r11+=fC11;
576 Double_t det=r00*r11 - r01*r01;
577 Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
578
579
580 Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
581 Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
582 Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
583 Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
584 Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
585
586 Double_t dy=cy - fP0, dz=cz - fP1;
587 Int_t layer = (index & 0xf0000000) >> 28;
588 fDy[layer] = dy;
589 fDz[layer] = dz;
590 fSigmaY[layer] = TMath::Sqrt(cerry*cerry+fC00);
591 fSigmaZ[layer] = TMath::Sqrt(cerrz*cerrz+fC11);
592
593 Double_t sf=fP2 + k20*dy + k21*dz;
594
595 fP0 += k00*dy + k01*dz;
596 fP1 += k10*dy + k11*dz;
597 fP2 = sf;
598 fP3 += k30*dy + k31*dz;
599 fP4 += k40*dy + k41*dz;
600
601 Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
602 Double_t c12=fC21, c13=fC31, c14=fC41;
603
604 fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
605 fC20-=k00*c02+k01*c12; fC30-=k00*c03+k01*c13;
606 fC40-=k00*c04+k01*c14;
607
608 fC11-=k10*c01+k11*fC11;
609 fC21-=k10*c02+k11*c12; fC31-=k10*c03+k11*c13;
610 fC41-=k10*c04+k11*c14;
611
612 fC22-=k20*c02+k21*c12; fC32-=k20*c03+k21*c13;
613 fC42-=k20*c04+k21*c14;
614
615 fC33-=k30*c03+k31*c13;
616 fC43-=k30*c04+k31*c14;
617
618 fC44-=k40*c04+k41*c14;
619
620 if (!Invariant()) {
621 fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
622 fC00=c00;
623 fC10=c10; fC11=c11;
624 fC20=c20; fC21=c21; fC22=c22;
625 fC30=c30; fC31=c31; fC32=c32; fC33=c33;
626 fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
627 return 0;
628 }
629
630 if (chi2<0) return 1;
631 Int_t n=GetNumberOfClusters();
632 fIndex[n]=index;
633 SetNumberOfClusters(n+1);
634 SetChi2(GetChi2()+chi2);
635
636 return 1;
637}
638
006b5f7f 639Int_t AliITStrackV2::Invariant() const {
640 //------------------------------------------------------------------
641 // This function is for debugging purpose only
642 //------------------------------------------------------------------
7f6ddf58 643 Int_t n=GetNumberOfClusters();
644
a9a2d814 645 if (TMath::Abs(fP2)>=0.9999){
8676d691 646 if (n>kWARN) Warning("Invariant","fP2=%f\n",fP2);
a9a2d814 647 return 0;
648 }
649 if (fC00<=0 || fC00>9.) {
8676d691 650 if (n>kWARN) Warning("Invariant","fC00=%f\n",fC00);
a9a2d814 651 return 0;
652 }
653 if (fC11<=0 || fC11>9.) {
8676d691 654 if (n>kWARN) Warning("Invariant","fC11=%f\n",fC11);
a9a2d814 655 return 0;
656 }
657 if (fC22<=0 || fC22>1.) {
8676d691 658 if (n>kWARN) Warning("Invariant","fC22=%f\n",fC22);
a9a2d814 659 return 0;
660 }
661 if (fC33<=0 || fC33>1.) {
8676d691 662 if (n>kWARN) Warning("Invariant","fC33=%f\n",fC33);
a9a2d814 663 return 0;
664 }
665 if (fC44<=0 || fC44>6e-5) {
8676d691 666 if (n>kWARN) Warning("Invariant","fC44=%f\n",fC44);
a9a2d814 667 return 0;
14825d5a 668 }
14825d5a 669 return 1;
006b5f7f 670}
671
672//____________________________________________________________________________
a9a2d814 673Int_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) {
006b5f7f 674 //------------------------------------------------------------------
675 //This function propagates a track
676 //------------------------------------------------------------------
33329428 677 Double_t alpha=fAlpha, x=fX;
006b5f7f 678 Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
679 Double_t c00=fC00;
680 Double_t c10=fC10, c11=fC11;
681 Double_t c20=fC20, c21=fC21, c22=fC22;
682 Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
683 Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
684
33329428 685 if (alp < -TMath::Pi()) alp += 2*TMath::Pi();
686 else if (alp >= TMath::Pi()) alp -= 2*TMath::Pi();
687 Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha);
688 Double_t sf=fP2, cf=TMath::Sqrt(1.- fP2*fP2);
006b5f7f 689
e43c066c 690 TMatrixD *T=0;
33329428 691 // **** rotation **********************
692 {
006b5f7f 693 fAlpha = alp;
33329428 694 fX = x*ca + p0*sa;
695 fP0= -x*sa + p0*ca;
696 fP2= sf*ca - cf*sa;
697
e43c066c 698 static TMatrixD C(5,5);
699 C(0,0)=c00;
700 C(1,0)=c10; C(1,1)=c11;
701 C(2,0)=c20; C(2,1)=c21; C(2,2)=c22;
702 C(3,0)=c30; C(3,1)=c31; C(3,2)=c32; C(3,3)=c33;
703 C(4,0)=c40; C(4,1)=c41; C(4,2)=c42; C(4,3)=c43; C(4,4)=c44;
704 C(0,1)=C(1,0);
705 C(0,2)=C(2,0); C(1,2)=C(2,1);
706 C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
707 C(0,4)=C(4,0); C(1,4)=C(4,1); C(2,4)=C(4,2); C(3,4)=C(4,3);
708
709
710 static TMatrixD F(6,5);
711 F(0,0)=sa;
712 F(1,0)=ca;
713 F(2,1)=F(4,3)=F(5,4)=1;
714 F(3,2)=ca + sf/cf*sa;
715
716 //TMatrixD tmp(C,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F));
717
718 static TMatrixD Ft(5,6);
719 Ft(0,0)=sa;
720 Ft(0,1)=ca;
721 Ft(1,2)=Ft(3,4)=Ft(4,5)=1;
722 Ft(2,3)=ca + sf/cf*sa;
723
724 TMatrixD tmp(C,TMatrixD::kMult,Ft);
725 T=new TMatrixD(F,TMatrixD::kMult,tmp);
726
727
33329428 728 }
006b5f7f 729
33329428 730 // **** translation ******************
731 {
732 Double_t dx=xk-fX;
006b5f7f 733 Double_t f1=fP2, f2=f1 + fP4*dx;
a9a2d814 734 if (TMath::Abs(f2) >= 0.9999) {
006b5f7f 735 Int_t n=GetNumberOfClusters();
7f6ddf58 736 if (n>kWARN)
8676d691 737 Warning("Propagate","Propagation failed (%d) !\n",n);
006b5f7f 738 return 0;
739 }
33329428 740 Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
006b5f7f 741
33329428 742 fX=xk;
006b5f7f 743 fP0 += dx*(f1+f2)/(r1+r2);
744 fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
745 fP2 += dx*fP4;
746
e43c066c 747 static TMatrixD F(5,6);
748 F(0,1)=F(1,2)=F(2,3)=F(3,4)=F(4,5)=1;
749 F(0,3)=dx/(r1+r2)*(2+(f1+f2)*(f2/r2+f1/r1)/(r1+r2));
750 F(0,5)=dx*dx/(r1+r2)*(1+(f1+f2)*f2/(r1+r2));
751 F(1,3)=dx*fP3/(f1*r2 + f2*r1)*(2-(f1+f2)*(r2-f1*f2/r2+r1-f2*f1/r1)/(f1*r2 + f2*r1));
752 F(1,4)=dx*(f1+f2)/(f1*r2 + f2*r1);
753 F(1,5)=dx*dx*fP3/(f1*r2 + f2*r1)*(1-(f1+f2)*(-f1*f2/r2+r1)/(f1*r2 + f2*r1));
754 F(2,5)=dx;
755 F(0,0)=-1/(r1+r2)*((f1+f2)+dx*fP4*(1+(f1+f2)/(r1+r2)*f2/r2));
756 F(1,0)=-fP3/(f1*r2 + f2*r1)*((f1+f2)+dx*fP4*(1+(f1+f2)/(f1*r2 + f2*r1)*(f1*f2/r2-r1)));
757 F(2,0)=-fP4;
758
759 TMatrixD tmp(*T,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F));
760 delete T;
761 TMatrixD C(F,TMatrixD::kMult,tmp);
762
763 fC00=C(0,0);
764 fC10=C(1,0); fC11=C(1,1);
765 fC20=C(2,0); fC21=C(2,1); fC22=C(2,2);
766 fC30=C(3,0); fC31=C(3,1); fC32=C(3,2); fC33=C(3,3);
767 fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4);
006b5f7f 768
006b5f7f 769 if (!Invariant()) {
33329428 770 fAlpha=alpha;
771 fX=x;
006b5f7f 772 fP0=p0; fP1=p1; fP2=p2; fP3=p3; fP4=p4;
773 fC00=c00;
774 fC10=c10; fC11=c11;
775 fC20=c20; fC21=c21; fC22=c22;
776 fC30=c30; fC31=c31; fC32=c32; fC33=c33;
777 fC40=c40; fC41=c41; fC42=c42; fC43=c43; fC44=c44;
778 return 0;
779 }
33329428 780 }
006b5f7f 781
782 return 1;
783}
784
e43c066c 785
786
787Int_t AliITStrackV2::GetProlongationFast(Double_t alp, Double_t xk,Double_t &y, Double_t &z)
788{
789 //-----------------------------------------------------------------------------
790 //get fast prolongation
791 //-----------------------------------------------------------------------------
792 Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha);
793 Double_t cf=TMath::Sqrt(1.- fP2*fP2);
794 // **** rotation **********************
795 y= -fX*sa + fP0*ca;
796 // **** translation ******************
797 Double_t dx = xk- fX*ca - fP0*sa;
798 Double_t f1=fP2*ca - cf*sa, f2=f1 + fP4*dx;
799 if (TMath::Abs(f2) >= 0.9999) {
800 return 0;
801 }
802 Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
803 y += dx*(f1+f2)/(r1+r2);
804 z = fP1+dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
805 return 1;
806}
807
808
a6a6f799 809Double_t AliITStrackV2::GetD(Double_t x, Double_t y) const {
006b5f7f 810 //------------------------------------------------------------------
a6a6f799 811 // This function calculates the transverse impact parameter
812 // with respect to a point with global coordinates (x,y)
006b5f7f 813 //------------------------------------------------------------------
a6a6f799 814 Double_t xt=fX, yt=fP0;
815
816 Double_t sn=TMath::Sin(fAlpha), cs=TMath::Cos(fAlpha);
817 Double_t a = x*cs + y*sn;
818 y = -x*sn + y*cs; x=a;
819 xt-=x; yt-=y;
820
821 sn=fP4*xt - fP2; cs=fP4*yt + TMath::Sqrt(1.- fP2*fP2);
822 a=2*(xt*fP2 - yt*TMath::Sqrt(1.- fP2*fP2))-fP4*(xt*xt + yt*yt);
006b5f7f 823 if (fP4<0) a=-a;
824 return a/(1 + TMath::Sqrt(sn*sn + cs*cs));
a9a2d814 825}
006b5f7f 826
babd135a 827Double_t AliITStrackV2::GetZat(Double_t x) const {
828 //------------------------------------------------------------------
829 // This function calculates the z at given x point - in current coordinate system
830 //------------------------------------------------------------------
831 Double_t x1=fX, x2=x, dx=x2-x1;
832 //
833 Double_t f1=fP2, f2=f1 + fP4*dx;
834 if (TMath::Abs(f2) >= 0.9999) {
835 return 10000000;
836 }
babd135a 837 Double_t r1=sqrt(1.- f1*f1), r2=sqrt(1.- f2*f2);
838 Double_t z = fP1 + dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
839 return z;
840}
841
842
843
844
8676d691 845Int_t AliITStrackV2::Improve(Double_t x0,Double_t xyz[3],Double_t ers[3]) {
a9a2d814 846 //------------------------------------------------------------------
847 //This function improves angular track parameters
848 //------------------------------------------------------------------
8676d691 849 Double_t cs=TMath::Cos(fAlpha), sn=TMath::Sin(fAlpha);
850 //Double_t xv = xyz[0]*cs + xyz[1]*sn; // vertex
851 Double_t yv =-xyz[0]*sn + xyz[1]*cs; // in the
852 Double_t zv = xyz[2]; // local frame
a9a2d814 853 Double_t dy=fP0-yv, dz=fP1-zv;
854 Double_t r2=fX*fX+dy*dy;
855 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
856 Double_t beta2=p2/(p2 + GetMass()*GetMass());
857 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
8676d691 858 Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
859 //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
a9a2d814 860 {
861 Double_t parp=0.5*(fP4*fX + dy*TMath::Sqrt(4/r2-fP4*fP4));
862 Double_t sigma2p = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
863 sigma2p += fC00/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
8676d691 864 sigma2p += ers[1]*ers[1]/r2;
a9a2d814 865 sigma2p += 0.25*fC44*fX*fX;
866 Double_t eps2p=sigma2p/(fC22+sigma2p);
867 fP0 += fC20/(fC22+sigma2p)*(parp-fP2);
868 fP2 = eps2p*fP2 + (1-eps2p)*parp;
869 fC22 *= eps2p;
870 fC20 *= eps2p;
871 }
872 {
873 Double_t parl=0.5*fP4*dz/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
874 Double_t sigma2l=theta2;
875 sigma2l += fC11/r2+fC00*dy*dy*dz*dz/(r2*r2*r2);
8676d691 876 sigma2l += ers[2]*ers[2]/r2;
a9a2d814 877 Double_t eps2l=sigma2l/(fC33+sigma2l);
878 fP1 += fC31/(fC33+sigma2l)*(parl-fP3);
879 fP4 += fC43/(fC33+sigma2l)*(parl-fP3);
880 fP3 = eps2l*fP3 + (1-eps2l)*parl;
881 fC33 *= eps2l; fC43 *= eps2l;
882 fC31 *= eps2l;
883 }
884 if (!Invariant()) return 0;
885 return 1;
886}
006b5f7f 887
a9a2d814 888/*
006b5f7f 889Int_t AliITStrackV2::Improve(Double_t x0,Double_t yv,Double_t zv) {
890 //------------------------------------------------------------------
891 //This function improves angular track parameters
892 //------------------------------------------------------------------
893 Double_t dy=fP0-yv, dz=fP1-zv;
894 Double_t r2=fX*fX+dy*dy;
895 Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
7f6ddf58 896 Double_t beta2=p2/(p2 + GetMass()*GetMass());
006b5f7f 897 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
a9a2d814 898 //Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
899 Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
006b5f7f 900
901 Double_t par=0.5*(fP4*fX + dy*TMath::Sqrt(4/r2-fP4*fP4));
902 Double_t sigma2 = theta2*(1.- GetSnp()*GetSnp())*(1. + GetTgl()*GetTgl());
903 sigma2 += fC00/r2*(1.- dy*dy/r2)*(1.- dy*dy/r2);
904 sigma2 += kSigmaYV*kSigmaYV/r2;
905 sigma2 += 0.25*fC44*fX*fX;
906 Double_t eps2=sigma2/(fC22+sigma2), eps=TMath::Sqrt(eps2);
907 if (10*r2*fC44<fC22) {
908 fP2 = eps2*fP2 + (1-eps2)*par;
909 fC22*=eps2; fC21*=eps; fC20*=eps; fC32*=eps; fC42*=eps;
910 }
911
912 par=0.5*fP4*dz/TMath::ASin(0.5*fP4*TMath::Sqrt(r2));
913 sigma2=theta2;
914 sigma2 += fC11/r2+fC00*dy*dy*dz*dz/(r2*r2*r2);
915 sigma2 += kSigmaZV*kSigmaZV/r2;
916 eps2=sigma2/(fC33+sigma2); eps=TMath::Sqrt(eps2);
917 Double_t tgl=fP3;
918 fP3 = eps2*fP3 + (1-eps2)*par;
919 fC33*=eps2; fC32*=eps; fC31*=eps; fC30*=eps; fC43*=eps;
920
921 eps=TMath::Sqrt((1+fP3*fP3)/(1+tgl*tgl));
922 fP4*=eps;
923 fC44*=eps*eps; fC43*=eps;fC42*=eps; fC41*=eps; fC40*=eps;
924
925 if (!Invariant()) return 0;
926 return 1;
927}
006b5f7f 928*/
14825d5a 929void AliITStrackV2::ResetCovariance() {
930 //------------------------------------------------------------------
931 //This function makes a track forget its history :)
932 //------------------------------------------------------------------
933
934 fC00*=10.;
935 fC10=0.; fC11*=10.;
936 fC20=0.; fC21=0.; fC22*=10.;
937 fC30=0.; fC31=0.; fC32=0.; fC33*=10.;
938 fC40=0.; fC41=0.; fC42=0.; fC43=0.; fC44*=10.;
939
940}
23efe5f1 941
942void AliITStrackV2::CookdEdx(Double_t low, Double_t up) {
943 //-----------------------------------------------------------------
a9a2d814 944 // This function calculates dE/dX within the "low" and "up" cuts.
945 // Origin: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
23efe5f1 946 //-----------------------------------------------------------------
23efe5f1 947 // The clusters order is: SSD-2, SSD-1, SDD-2, SDD-1, SPD-2, SPD-1
c0dd5278 948
949 Int_t i;
950 Int_t nc=0;
951 for (i=0; i<GetNumberOfClusters(); i++) {
952 Int_t idx=GetClusterIndex(i);
953 idx=(idx&0xf0000000)>>28;
954 if (idx>1) nc++; // Take only SSD and SDD
955 }
23efe5f1 956
a9a2d814 957 Int_t swap;//stupid sorting
23efe5f1 958 do {
959 swap=0;
960 for (i=0; i<nc-1; i++) {
961 if (fdEdxSample[i]<=fdEdxSample[i+1]) continue;
962 Float_t tmp=fdEdxSample[i];
963 fdEdxSample[i]=fdEdxSample[i+1]; fdEdxSample[i+1]=tmp;
964 swap++;
965 }
966 } while (swap);
967
a9a2d814 968 Int_t nl=Int_t(low*nc), nu=Int_t(up*nc); //b.b. to take two lowest dEdX
969 // values from four ones choose
970 // nu=2
23efe5f1 971 Float_t dedx=0;
a9a2d814 972 for (i=nl; i<nu; i++) dedx += fdEdxSample[i];
c0dd5278 973 if (nu-nl>0) dedx /= (nu-nl);
23efe5f1 974
23efe5f1 975 SetdEdx(dedx);
976}