We noticed that we did not provide in AliESDtrackCuts the methods with the "standard...
[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
ae00569a 16/* $Id$ */
17
66e811e2 18///////////////////////////////////////////////////////////////////////////
006b5f7f 19// Implementation of the ITS track class
20//
21// Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
a9a2d814 22// dEdx analysis by: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
66e811e2 23///////////////////////////////////////////////////////////////////////////
006b5f7f 24#include <TMath.h>
006b5f7f 25
26#include "AliCluster.h"
58e536c5 27#include "AliESDVertex.h"
8a84886c 28#include "AliITSReconstructor.h"
006b5f7f 29#include "AliITStrackV2.h"
f7a1cc68 30#include "AliTracker.h"
a4354152 31#include "AliLog.h"
006b5f7f 32
8602c008 33const Int_t AliITStrackV2::fgkWARN = 5;
34
1aedd96e 35ClassImp(AliITStrackV2)
8676d691 36
e43c066c 37
006b5f7f 38//____________________________________________________________________________
6c94f330 39AliITStrackV2::AliITStrackV2() : AliKalmanTrack(),
545af9ef 40 fCheckInvariant(kTRUE),
22b13da0 41 fdEdx(0),
83d73500 42 fESDtrack(0)
e43c066c 43{
ae00569a 44 for(Int_t i=0; i<2*AliITSgeomTGeo::kNLayers; i++) {fIndex[i]=-1; fModule[i]=-1;}
25015f7a 45 for(Int_t i=0; i<AliITSgeomTGeo::kNLayers; i++) {fSharedWeight[i]=0;}
ae00569a 46 for(Int_t i=0; i<4; i++) fdEdxSample[i]=0;
e43c066c 47}
48
83d73500 49
22b13da0 50//____________________________________________________________________________
763bdf3e 51AliITStrackV2::AliITStrackV2(AliESDtrack& t,Bool_t c):
6c94f330 52 AliKalmanTrack(),
545af9ef 53 fCheckInvariant(kTRUE),
6c94f330 54 fdEdx(t.GetITSsignal()),
55 fESDtrack(&t)
56{
83d73500 57 //------------------------------------------------------------------
67c3dcbe 58 // Conversion ESD track -> ITS track.
59 // If c==kTRUE, create the ITS track out of the constrained params.
83d73500 60 //------------------------------------------------------------------
6c94f330 61 const AliExternalTrackParam *par=&t;
62 if (c) {
6c4ef2ed 63 par=t.GetConstrainedParam();
763bdf3e 64 if (!par) AliError("AliITStrackV2: conversion failed !\n");
6c94f330 65 }
66 Set(par->GetX(),par->GetAlpha(),par->GetParameter(),par->GetCovariance());
67
83d73500 68 SetLabel(t.GetLabel());
69 SetMass(t.GetMass());
6c94f330 70 SetNumberOfClusters(t.GetITSclusters(fIndex));
006b5f7f 71
83d73500 72 if (t.GetStatus()&AliESDtrack::kTIME) {
73 StartTimeIntegral();
74 Double_t times[10]; t.GetIntegratedTimes(times); SetIntegratedTimes(times);
75 SetIntegratedLength(t.GetIntegratedLength());
76 }
e43c066c 77
25015f7a 78 for(Int_t i=0; i<AliITSgeomTGeo::kNLayers; i++) {fSharedWeight[i]=0;}
68b8060b 79 for(Int_t i=0; i<4; i++) fdEdxSample[i]=0;
006b5f7f 80}
81
7ac23097 82//____________________________________________________________________________
83void AliITStrackV2::ResetClusters() {
84 //------------------------------------------------------------------
85 // Reset the array of attached clusters.
86 //------------------------------------------------------------------
87 for (Int_t i=0; i<2*AliITSgeomTGeo::kNLayers; i++) fIndex[i]=-1;
25015f7a 88 for (Int_t i=0; i<AliITSgeomTGeo::kNLayers; i++) {fSharedWeight[i]=0;}
7ac23097 89 SetChi2(0.);
90 SetNumberOfClusters(0);
91}
92
09cf9d66 93//____________________________________________________________________________
15dd636f 94void AliITStrackV2::UpdateESDtrack(ULong_t flags) const {
09cf9d66 95 // Update track params
83d73500 96 fESDtrack->UpdateTrackParams(this,flags);
ae00569a 97 // copy the module indices
25015f7a 98 Int_t i;
99 for(i=0;i<2*AliITSgeomTGeo::kNLayers;i++) {
ae00569a 100 // printf(" %d\n",GetModuleIndex(i));
101 fESDtrack->SetITSModuleIndex(i,GetModuleIndex(i));
102 }
25015f7a 103 // copy the map of shared clusters
104 if(flags==AliESDtrack::kITSin) {
105 UChar_t itsSharedMap=0;
106 for(i=0;i<AliITSgeomTGeo::kNLayers;i++) {
63126a46 107 if(fSharedWeight[i]>0) SETBIT(itsSharedMap,i);
25015f7a 108
109 }
110 fESDtrack->SetITSSharedMap(itsSharedMap);
111 }
25015f7a 112
09cf9d66 113 // copy the 4 dedx samples
114 Double_t sdedx[4]={0.,0.,0.,0.};
25015f7a 115 for(i=0; i<4; i++) sdedx[i]=fdEdxSample[i];
09cf9d66 116 fESDtrack->SetITSdEdxSamples(sdedx);
83d73500 117}
118
006b5f7f 119//____________________________________________________________________________
6c94f330 120AliITStrackV2::AliITStrackV2(const AliITStrackV2& t) :
121 AliKalmanTrack(t),
545af9ef 122 fCheckInvariant(t.fCheckInvariant),
6c94f330 123 fdEdx(t.fdEdx),
124 fESDtrack(t.fESDtrack)
125{
006b5f7f 126 //------------------------------------------------------------------
127 //Copy constructor
128 //------------------------------------------------------------------
ef7253ac 129 Int_t i;
ef7253ac 130 for (i=0; i<4; i++) fdEdxSample[i]=t.fdEdxSample[i];
ae00569a 131 for (i=0; i<2*AliITSgeomTGeo::GetNLayers(); i++) {
132 fIndex[i]=t.fIndex[i];
133 fModule[i]=t.fModule[i];
134 }
25015f7a 135 for (i=0; i<AliITSgeomTGeo::kNLayers; i++) {fSharedWeight[i]=t.fSharedWeight[i];}
006b5f7f 136}
a9a2d814 137
006b5f7f 138//_____________________________________________________________________________
139Int_t AliITStrackV2::Compare(const TObject *o) const {
140 //-----------------------------------------------------------------
141 // This function compares tracks according to the their curvature
142 //-----------------------------------------------------------------
7f6ddf58 143 AliITStrackV2 *t=(AliITStrackV2*)o;
6c23ffed 144 //Double_t co=OneOverPt();
145 //Double_t c =OneOverPt();
15dd636f 146 Double_t co=t->GetSigmaY2()*t->GetSigmaZ2();
147 Double_t c =GetSigmaY2()*GetSigmaZ2();
006b5f7f 148 if (c>co) return 1;
149 else if (c<co) return -1;
150 return 0;
151}
152
006b5f7f 153//____________________________________________________________________________
6c94f330 154Bool_t
155AliITStrackV2::PropagateToVertex(const AliESDVertex *v,Double_t d,Double_t x0)
156{
006b5f7f 157 //------------------------------------------------------------------
158 //This function propagates a track to the minimal distance from the origin
6c94f330 159 //------------------------------------------------------------------
160 Double_t bz=GetBz();
e50912db 161 if (PropagateToDCA(v,bz,kVeryBig)) {
162 Double_t xOverX0,xTimesRho;
163 xOverX0 = d; xTimesRho = d*x0;
164 if (CorrectForMeanMaterial(xOverX0,xTimesRho,kTRUE)) return kTRUE;
165 }
6c94f330 166 return kFALSE;
006b5f7f 167}
168
169//____________________________________________________________________________
6c94f330 170Bool_t AliITStrackV2::
e50912db 171GetGlobalXYZat(Double_t xloc, Double_t &x, Double_t &y, Double_t &z) const {
006b5f7f 172 //------------------------------------------------------------------
173 //This function returns a track position in the global system
174 //------------------------------------------------------------------
6c94f330 175 Double_t r[3];
f7a1cc68 176 Bool_t rc=GetXYZAt(xloc, GetBz(), r);
6c94f330 177 x=r[0]; y=r[1]; z=r[2];
178 return rc;
006b5f7f 179}
180
181//_____________________________________________________________________________
6c94f330 182Double_t AliITStrackV2::GetPredictedChi2(const AliCluster *c) const {
006b5f7f 183 //-----------------------------------------------------------------
184 // This function calculates a predicted chi2 increment.
185 //-----------------------------------------------------------------
6c94f330 186 Double_t p[2]={c->GetY(), c->GetZ()};
187 Double_t cov[3]={c->GetSigmaY2(), 0., c->GetSigmaZ2()};
188 return AliExternalTrackParam::GetPredictedChi2(p,cov);
a9a2d814 189}
190
191//____________________________________________________________________________
6c94f330 192Bool_t AliITStrackV2::PropagateTo(Double_t xk, Double_t d, Double_t x0) {
006b5f7f 193 //------------------------------------------------------------------
194 //This function propagates a track
195 //------------------------------------------------------------------
afd25725 196
6c94f330 197 Double_t oldX=GetX(), oldY=GetY(), oldZ=GetZ();
006b5f7f 198
1defa402 199 //Double_t bz=GetBz();
200 //if (!AliExternalTrackParam::PropagateTo(xk,bz)) return kFALSE;
201 Double_t b[3]; GetBxByBz(b);
202 if (!AliExternalTrackParam::PropagateToBxByBz(xk,b)) return kFALSE;
e50912db 203 Double_t xOverX0,xTimesRho;
204 xOverX0 = d; xTimesRho = d*x0;
205 if (!CorrectForMeanMaterial(xOverX0,xTimesRho,kTRUE)) return kFALSE;
006b5f7f 206
dd44614b 207 Double_t x=GetX(), y=GetY(), z=GetZ();
6c94f330 208 if (IsStartedTimeIntegral() && x>oldX) {
209 Double_t l2 = (x-oldX)*(x-oldX) + (y-oldY)*(y-oldY) + (z-oldZ)*(z-oldZ);
f29dbb4b 210 AddTimeStep(TMath::Sqrt(l2));
211 }
f29dbb4b 212
6c94f330 213 return kTRUE;
006b5f7f 214}
215
216//____________________________________________________________________________
db355ee7 217Bool_t AliITStrackV2::PropagateToTGeo(Double_t xToGo, Int_t nstep, Double_t &xOverX0, Double_t &xTimesRho, Bool_t addTime) {
afd25725 218 //-------------------------------------------------------------------
219 // Propagates the track to a reference plane x=xToGo in n steps.
220 // These n steps are only used to take into account the curvature.
221 // The material is calculated with TGeo. (L.Gaudichet)
222 //-------------------------------------------------------------------
223
224 Double_t startx = GetX(), starty = GetY(), startz = GetZ();
225 Double_t sign = (startx<xToGo) ? -1.:1.;
226 Double_t step = (xToGo-startx)/TMath::Abs(nstep);
227
1defa402 228 Double_t start[3], end[3], mparam[7];
229 //Double_t bz = GetBz();
230 Double_t b[3]; GetBxByBz(b);
231 Double_t bz = b[2];
232
afd25725 233 Double_t x = startx;
234
235 for (Int_t i=0; i<nstep; i++) {
236
237 GetXYZ(start); //starting global position
238 x += step;
239 if (!GetXYZAt(x, bz, end)) return kFALSE;
1defa402 240 //if (!AliExternalTrackParam::PropagateTo(x, bz)) return kFALSE;
241 if (!AliExternalTrackParam::PropagateToBxByBz(x, b)) return kFALSE;
afd25725 242 AliTracker::MeanMaterialBudget(start, end, mparam);
8fc9f3dc 243 xTimesRho = sign*mparam[4]*mparam[0];
244 xOverX0 = mparam[1];
afd25725 245 if (mparam[1]<900000) {
afd25725 246 if (!AliExternalTrackParam::CorrectForMeanMaterial(xOverX0,
8fc9f3dc 247 xTimesRho,GetMass())) return kFALSE;
248 } else { // this happens when MeanMaterialBudget cannot cross a boundary
249 return kFALSE;
afd25725 250 }
251 }
252
db355ee7 253 if (addTime && IsStartedTimeIntegral() && GetX()>startx) {
afd25725 254 Double_t l2 = ( (GetX()-startx)*(GetX()-startx) +
255 (GetY()-starty)*(GetY()-starty) +
256 (GetZ()-startz)*(GetZ()-startz) );
257 AddTimeStep(TMath::Sqrt(l2));
258 }
259
260 return kTRUE;
261}
262
263//____________________________________________________________________________
6c94f330 264Bool_t AliITStrackV2::Update(const AliCluster* c, Double_t chi2, Int_t index)
265{
006b5f7f 266 //------------------------------------------------------------------
267 //This function updates track parameters
268 //------------------------------------------------------------------
6c94f330 269 Double_t p[2]={c->GetY(), c->GetZ()};
d9ead1a0 270 Double_t cov[3]={c->GetSigmaY2(), c->GetSigmaYZ(), c->GetSigmaZ2()};
006b5f7f 271
6c94f330 272 if (!AliExternalTrackParam::Update(p,cov)) return kFALSE;
006b5f7f 273
28d5b13d 274 Int_t n=GetNumberOfClusters();
006b5f7f 275 if (!Invariant()) {
28d5b13d 276 if (n>fgkWARN) AliWarning("Wrong invariant !");
6c94f330 277 return kFALSE;
006b5f7f 278 }
279
6c94f330 280 if (chi2<0) return kTRUE;
67c3dcbe 281
b5205c90 282 // fill residuals for ITS+TPC tracks
badb951a 283 if (fESDtrack) {
284 if (fESDtrack->GetStatus()&AliESDtrack::kTPCin) {
285 AliTracker::FillResiduals(this,p,cov,c->GetVolumeId());
286 }
b5205c90 287 }
150f264c 288
006b5f7f 289 fIndex[n]=index;
290 SetNumberOfClusters(n+1);
291 SetChi2(GetChi2()+chi2);
292
6c94f330 293 return kTRUE;
006b5f7f 294}
295
6c94f330 296Bool_t AliITStrackV2::Invariant() const {
006b5f7f 297 //------------------------------------------------------------------
298 // This function is for debugging purpose only
299 //------------------------------------------------------------------
545af9ef 300 if(!fCheckInvariant) return kTRUE;
301
7f6ddf58 302 Int_t n=GetNumberOfClusters();
a4354152 303 static Float_t bz = GetBz();
8a84886c 304 // take into account the misalignment error
305 Float_t maxMisalErrY2=0,maxMisalErrZ2=0;
a4354152 306 //RS
307 const AliITSRecoParam* recopar = AliITSReconstructor::GetRecoParam();
308 if (!recopar) recopar = AliITSRecoParam::GetHighFluxParam();
309
8a84886c 310 for (Int_t lay=0; lay<AliITSgeomTGeo::kNLayers; lay++) {
a4354152 311 maxMisalErrY2 = TMath::Max(maxMisalErrY2,recopar->GetClusterMisalErrorY(lay,bz));
312 maxMisalErrZ2 = TMath::Max(maxMisalErrZ2,recopar->GetClusterMisalErrorZ(lay,bz));
8a84886c 313 }
314 maxMisalErrY2 *= maxMisalErrY2;
315 maxMisalErrZ2 *= maxMisalErrZ2;
316 // this is because when we reset before refitting, we multiply the
317 // matrix by 10
318 maxMisalErrY2 *= 10.;
319 maxMisalErrZ2 *= 10.;
320
6c94f330 321 Double_t sP2=GetParameter()[2];
322 if (TMath::Abs(sP2) >= kAlmost1){
8602c008 323 if (n>fgkWARN) Warning("Invariant","fP2=%f\n",sP2);
6c94f330 324 return kFALSE;
a9a2d814 325 }
6c94f330 326 Double_t sC00=GetCovariance()[0];
8a84886c 327 if (sC00<=0 || sC00>(9.+maxMisalErrY2)) {
8602c008 328 if (n>fgkWARN) Warning("Invariant","fC00=%f\n",sC00);
6c94f330 329 return kFALSE;
a9a2d814 330 }
6c94f330 331 Double_t sC11=GetCovariance()[2];
8a84886c 332 if (sC11<=0 || sC11>(9.+maxMisalErrZ2)) {
8602c008 333 if (n>fgkWARN) Warning("Invariant","fC11=%f\n",sC11);
6c94f330 334 return kFALSE;
a9a2d814 335 }
6c94f330 336 Double_t sC22=GetCovariance()[5];
337 if (sC22<=0 || sC22>1.) {
8602c008 338 if (n>fgkWARN) Warning("Invariant","fC22=%f\n",sC22);
6c94f330 339 return kFALSE;
a9a2d814 340 }
6c94f330 341 Double_t sC33=GetCovariance()[9];
342 if (sC33<=0 || sC33>1.) {
8602c008 343 if (n>fgkWARN) Warning("Invariant","fC33=%f\n",sC33);
6c94f330 344 return kFALSE;
a9a2d814 345 }
6c94f330 346 Double_t sC44=GetCovariance()[14];
347 if (sC44<=0 /*|| sC44>6e-5*/) {
8602c008 348 if (n>fgkWARN) Warning("Invariant","fC44=%f\n",sC44);
6c94f330 349 return kFALSE;
14825d5a 350 }
6c94f330 351
352 return kTRUE;
006b5f7f 353}
354
355//____________________________________________________________________________
6c94f330 356Bool_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) {
006b5f7f 357 //------------------------------------------------------------------
358 //This function propagates a track
359 //------------------------------------------------------------------
1defa402 360 //Double_t bz=GetBz();
361 //if (!AliExternalTrackParam::Propagate(alp,xk,bz)) return kFALSE;
362 Double_t b[3]; GetBxByBz(b);
363 if (!AliExternalTrackParam::PropagateBxByBz(alp,xk,b)) return kFALSE;
006b5f7f 364
006b5f7f 365 if (!Invariant()) {
8a84886c 366 Int_t n=GetNumberOfClusters();
367 if (n>fgkWARN) AliWarning("Wrong invariant !");
368 return kFALSE;
791f9a2a 369 }
370
6c94f330 371 return kTRUE;
a9a2d814 372}
006b5f7f 373
afd25725 374Bool_t AliITStrackV2::MeanBudgetToPrimVertex(Double_t xyz[3], Double_t step, Double_t &d) const {
375
376 //-------------------------------------------------------------------
377 // Get the mean material budget between the actual point and the
378 // primary vertex. (L.Gaudichet)
379 //-------------------------------------------------------------------
380
381 Double_t cs=TMath::Cos(GetAlpha()), sn=TMath::Sin(GetAlpha());
382 Double_t vertexX = xyz[0]*cs + xyz[1]*sn;
383
384 Int_t nstep = Int_t((GetX()-vertexX)/step);
385 if (nstep<1) nstep = 1;
386 step = (GetX()-vertexX)/nstep;
387
c7719399 388 // Double_t mparam[7], densMean=0, radLength=0, length=0;
389 Double_t mparam[7];
afd25725 390 Double_t p1[3], p2[3], x = GetX(), bz = GetBz();
391 GetXYZ(p1);
392
393 d=0.;
394
395 for (Int_t i=0; i<nstep; i++) {
396 x += step;
397 if (!GetXYZAt(x, bz, p2)) return kFALSE;
398 AliTracker::MeanMaterialBudget(p1, p2, mparam);
399 if (mparam[1]>900000) return kFALSE;
400 d += mparam[1];
401
402 p1[0] = p2[0];
403 p1[1] = p2[1];
404 p1[2] = p2[2];
405 }
406
407 return kTRUE;
408}
409
6c94f330 410Bool_t AliITStrackV2::Improve(Double_t x0,Double_t xyz[3],Double_t ers[3]) {
babd135a 411 //------------------------------------------------------------------
6c94f330 412 //This function improves angular track parameters
babd135a 413 //------------------------------------------------------------------
86be8934 414 //Store the initail track parameters
415
86be8934 416 Double_t x = GetX();
417 Double_t alpha = GetAlpha();
418 Double_t par[] = {GetY(),GetZ(),GetSnp(),GetTgl(),GetSigned1Pt()};
419 Double_t cov[] = {
420 GetSigmaY2(),
421 GetSigmaZY(),
422 GetSigmaZ2(),
423 GetSigmaSnpY(),
424 GetSigmaSnpZ(),
425 GetSigmaSnp2(),
426 GetSigmaTglY(),
427 GetSigmaTglZ(),
428 GetSigmaTglSnp(),
429 GetSigmaTgl2(),
430 GetSigma1PtY(),
431 GetSigma1PtZ(),
432 GetSigma1PtSnp(),
433 GetSigma1PtTgl(),
434 GetSigma1Pt2()
435 };
436
437
6c94f330 438 Double_t cs=TMath::Cos(GetAlpha()), sn=TMath::Sin(GetAlpha());
cb4706fe 439 Double_t xv = xyz[0]*cs + xyz[1]*sn; // vertex
6c94f330 440 Double_t yv =-xyz[0]*sn + xyz[1]*cs; // in the
441 Double_t zv = xyz[2]; // local frame
babd135a 442
cb4706fe 443 Double_t dx = x - xv, dy = par[0] - yv, dz = par[1] - zv;
444 Double_t r2=dx*dx + dy*dy;
6c23ffed 445 Double_t p2=(1.+ GetTgl()*GetTgl())/(GetSigned1Pt()*GetSigned1Pt());
a9a2d814 446 Double_t beta2=p2/(p2 + GetMass()*GetMass());
447 x0*=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(1.- GetSnp()*GetSnp()));
8676d691 448 Double_t theta2=14.1*14.1/(beta2*p2*1e6)*x0;
449 //Double_t theta2=1.0259e-6*14*14/28/(beta2*p2)*x0*9.36*2.33;
6c94f330 450
63126a46 451 Double_t bz=GetBz();
452 Double_t cnv=bz*kB2C;
453 Double_t curv=GetC(bz);
a9a2d814 454 {
63126a46 455 Double_t dummy = 4/r2 - curv*curv;
6c94f330 456 if (dummy < 0) return kFALSE;
63126a46 457 Double_t parp = 0.5*(curv*dx + dy*TMath::Sqrt(dummy));
bfd20868 458 Double_t sigma2p = theta2*(1.-GetSnp())*(1.+GetSnp())*(1. + GetTgl()*GetTgl());
459 Double_t ovSqr2 = 1./TMath::Sqrt(r2);
460 Double_t tfact = ovSqr2*(1.-dy*ovSqr2)*(1.+dy*ovSqr2);
86be8934 461 sigma2p += cov[0]*tfact*tfact;
6c94f330 462 sigma2p += ers[1]*ers[1]/r2;
cb4706fe 463 sigma2p += 0.25*cov[14]*cnv*cnv*dx*dx;
86be8934 464 Double_t eps2p=sigma2p/(cov[5] + sigma2p);
465 par[0] += cov[3]/(cov[5] + sigma2p)*(parp - GetSnp());
466 par[2] = eps2p*GetSnp() + (1 - eps2p)*parp;
467 cov[5] *= eps2p;
468 cov[3] *= eps2p;
a9a2d814 469 }
470 {
63126a46 471 Double_t parl=0.5*curv*dz/TMath::ASin(0.5*curv*TMath::Sqrt(r2));
6c94f330 472 Double_t sigma2l=theta2;
86be8934 473 sigma2l += cov[2]/r2 + cov[0]*dy*dy*dz*dz/(r2*r2*r2);
6c94f330 474 sigma2l += ers[2]*ers[2]/r2;
86be8934 475 Double_t eps2l = sigma2l/(cov[9] + sigma2l);
476 par[1] += cov[7 ]/(cov[9] + sigma2l)*(parl - par[3]);
477 par[4] += cov[13]/(cov[9] + sigma2l)*(parl - par[3]);
478 par[3] = eps2l*par[3] + (1-eps2l)*parl;
479 cov[9] *= eps2l;
480 cov[13]*= eps2l;
481 cov[7] *= eps2l;
a9a2d814 482 }
86be8934 483
484 Set(x,alpha,par,cov);
485
6c94f330 486 if (!Invariant()) return kFALSE;
14825d5a 487
6c94f330 488 return kTRUE;
14825d5a 489}
23efe5f1 490
491void AliITStrackV2::CookdEdx(Double_t low, Double_t up) {
492 //-----------------------------------------------------------------
a9a2d814 493 // This function calculates dE/dX within the "low" and "up" cuts.
494 // Origin: Boris Batyunya, JINR, Boris.Batiounia@cern.ch
09cf9d66 495 // Updated: F. Prino 8-June-2009
23efe5f1 496 //-----------------------------------------------------------------
09cf9d66 497 // The cluster order is: SDD-1, SDD-2, SSD-1, SSD-2
c0dd5278 498
c0dd5278 499 Int_t nc=0;
09cf9d66 500 Float_t dedx[4];
501 for (Int_t il=0; il<4; il++) { // count good (>0) dE/dx values
502 if(fdEdxSample[il]>0.){
503 dedx[nc]= fdEdxSample[il];
504 nc++;
505 }
506 }
507 if(nc<1){
508 SetdEdx(0.);
509 return;
c0dd5278 510 }
23efe5f1 511
09cf9d66 512 Int_t swap; // sort in ascending order
23efe5f1 513 do {
514 swap=0;
09cf9d66 515 for (Int_t i=0; i<nc-1; i++) {
516 if (dedx[i]<=dedx[i+1]) continue;
517 Float_t tmp=dedx[i];
518 dedx[i]=dedx[i+1];
519 dedx[i+1]=tmp;
23efe5f1 520 swap++;
521 }
522 } while (swap);
523
23efe5f1 524
09cf9d66 525 Double_t nl=low*nc, nu =up*nc;
526 Float_t sumamp = 0;
527 Float_t sumweight =0;
528 for (Int_t i=0; i<nc; i++) {
529 Float_t weight =1;
530 if (i<nl+0.1) weight = TMath::Max(1.-(nl-i),0.);
531 if (i>nu-1) weight = TMath::Max(nu-i,0.);
532 sumamp+= dedx[i]*weight;
533 sumweight+=weight;
534 }
535 SetdEdx(sumamp/sumweight);
23efe5f1 536}
c7bafca9 537
8602c008 538//____________________________________________________________________________
539Bool_t AliITStrackV2::
e50912db 540GetPhiZat(Double_t r, Double_t &phi, Double_t &z) const {
541 //------------------------------------------------------------------
542 // This function returns the global cylindrical (phi,z) of the track
543 // position estimated at the radius r.
544 // The track curvature is neglected.
545 //------------------------------------------------------------------
546 Double_t d=GetD(0.,0.);
bc02d571 547 if (TMath::Abs(d) > r) {
548 if (r>1e-1) return kFALSE;
549 r = TMath::Abs(d);
550 }
e50912db 551
552 Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
494a930c 553 if (TMath::Abs(d) > rcurr) return kFALSE;
554 Double_t globXYZcurr[3]; GetXYZ(globXYZcurr);
555 Double_t phicurr=TMath::ATan2(globXYZcurr[1],globXYZcurr[0]);
e50912db 556
2089d338 557 if (GetX()>=0.) {
558 phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr);
559 } else {
560 phi=phicurr+TMath::ASin(d/r)+TMath::ASin(d/rcurr)-TMath::Pi();
561 }
562
494a930c 563 // return a phi in [0,2pi[
564 if (phi<0.) phi+=2.*TMath::Pi();
565 else if (phi>=2.*TMath::Pi()) phi-=2.*TMath::Pi();
c7de5c65 566 z=GetZ()+GetTgl()*(TMath::Sqrt((r-d)*(r+d))-TMath::Sqrt((rcurr-d)*(rcurr+d)));
e50912db 567 return kTRUE;
568}
569//____________________________________________________________________________
570Bool_t AliITStrackV2::
571GetLocalXat(Double_t r,Double_t &xloc) const {
8602c008 572 //------------------------------------------------------------------
e50912db 573 // This function returns the local x of the track
574 // position estimated at the radius r.
8602c008 575 // The track curvature is neglected.
576 //------------------------------------------------------------------
577 Double_t d=GetD(0.,0.);
6a14afad 578 if (TMath::Abs(d) > r) {
579 if (r>1e-1) return kFALSE;
580 r = TMath::Abs(d);
581 }
8602c008 582
e50912db 583 Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
494a930c 584 Double_t globXYZcurr[3]; GetXYZ(globXYZcurr);
585 Double_t phicurr=TMath::ATan2(globXYZcurr[1],globXYZcurr[0]);
2089d338 586 Double_t phi;
587 if (GetX()>=0.) {
588 phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr);
589 } else {
590 phi=phicurr+TMath::ASin(d/r)+TMath::ASin(d/rcurr)-TMath::Pi();
591 }
8602c008 592
e50912db 593 xloc=r*(TMath::Cos(phi)*TMath::Cos(GetAlpha())
594 +TMath::Sin(phi)*TMath::Sin(GetAlpha()));
8602c008 595
596 return kTRUE;
597}
63126a46 598
a4354152 599//____________________________________________________________________________
600Bool_t AliITStrackV2::ImproveKalman(Double_t xyz[3],Double_t ers[3], const Double_t* xlMS, const Double_t* x2X0MS, Int_t nMS)
601{
602 // Substitute the state of the track (p_{k|k},C_{k|k}) at the k-th measumerent by its
603 // smoothed value from the k-th measurement + measurement at the vertex.
604 // Account for the MS on nMS layers at x-postions xlMS with x/x0 = x2X0MS
605 // p_{k|kv} = p_{k|k} + C_{k|k}*D^Tr_{k+1} B^{-1}_{k+1} ( vtx - D_{k+1}*p_{k|k})
606 // C_{k|kv} = C_{k|k}*( I - D^Tr_{k+1} B^{-1}_{k+1} D_{k+1} C^Tr_{k|k})
607 //
608 // where D_{k} = H_{k} F_{k} with H being the matrix converting the tracks parameters
609 // to measurements m_{k} = H_{k} p_{k} and F_{k} the matrix propagating the track between the
610 // the point k-1 and k: p_{k|k-1} = F_{k} p_{k-1|k-1}
611 //
612 // B_{k+1} = V_{k+1} + H_{k+1} C_{k+1|k} H^Tr_{k+1} with V_{k+1} being the error of the measurment
613 // at point k+1 (i.e. vertex), and C_{k+1|k} - error matrix extrapolated from k-th measurement to
614 // k+1 (vtx) and accounting for the MS inbetween
615 //
616 // H = {{1,0,0,0,0},{0,1,0,0,0}}
617 //
618 double covc[15], *cori = (double*) GetCovariance(),par[5] = {GetY(),GetZ(),GetSnp(),GetTgl(),GetSigned1Pt()},
619 &c00=cori[0],
620 &c01=cori[1],&c11=cori[2],
621 &c02=cori[3],&c12=cori[4],&c22=cori[5],
622 &c03=cori[6],&c13=cori[7],&c23=cori[8],&c33=cori[9],
623 &c04=cori[10],&c14=cori[11],&c24=cori[12],&c34=cori[13],&c44=cori[14],
624 // for smoothed cov matrix
625 &cov00=covc[0],
626 &cov01=covc[1],&cov11=covc[2],
627 &cov02=covc[3],&cov12=covc[4],&cov22=covc[5],
628 &cov03=covc[6],&cov13=covc[7],&cov23=covc[8],&cov33=covc[9],
629 &cov04=covc[10],&cov14=covc[11],&cov24=covc[12],&cov34=covc[13],&cov44=covc[14];
630 //
631 double x = GetX(), alpha = GetAlpha();
632 // vertex in the track frame
633 double cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
634 double xv = xyz[0]*cs + xyz[1]*sn, yv =-xyz[0]*sn + xyz[1]*cs, zv = xyz[2];
635 double dx = xv - GetX();
636 if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
637 //
638 double cnv=GetBz()*kB2C, x2r=cnv*par[4]*dx, f1=par[2], f2=f1+x2r;
639 if (TMath::Abs(f1) >= kAlmost1 || TMath::Abs(f2) >= kAlmost1) {
640 AliInfo(Form("Fail: %+e %+e",f1,f2));
641 return kFALSE;
642 }
643 double r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2)), dx2r=dx/(r1+r2);
644 // elements of matrix F_{k+1} (1s on diagonal)
645 double f02 = 2*dx2r, f04 = cnv*dx*dx2r, f13/*, f24 = cnv*dx*/;
646 if (TMath::Abs(x2r)<0.05) f13 = dx*r2+f2*(f1+f2)*dx2r; // see AliExternalTrackParam::PropagateTo
647 else {
648 double dy2dx = (f1+f2)/(r1+r2);
649 f13 = 2*TMath::ASin(0.5*TMath::Sqrt(1+dy2dx*dy2dx)*x2r)/(cnv*par[4]);
650 }
651 // elements of matrix D_{k+1} = H_{k+1} * F_{k+1}
652 // double d00 = 1., d11 = 1.;
653 double &d02 = f02, &d04 = f04, &d13 = f13;
654 //
655 // elements of matrix DC = D_{k+1}*C_{kk}^T
656 double dc00 = c00+c02*d02+c04*d04, dc10 = c01+c03*d13;
657 double dc01 = c01+c12*d02+c14*d04, dc11 = c11+c13*d13;
658 double dc02 = c02+c22*d02+c24*d04, dc12 = c12+c23*d13;
659 double dc03 = c03+c23*d02+c34*d04, dc13 = c13+c33*d13;
660 double dc04 = c04+c24*d02+c44*d04, dc14 = c14+c34*d13;
661 //
662 // difference between the vertex and the the track extrapolated to vertex
663 yv -= par[0] + par[2]*d02 + par[4]*d04;
664 zv -= par[1] + par[3]*d13;
665 //
666 // y,z part of the cov.matrix extrapolated to vtx (w/o MS contribution)
667 // C_{k+1,k} = H F_{k+1} C_{k,k} F^Tr_{k+1} H^Tr = D C D^Tr
668 double cv00 = dc00+dc02*d02+dc04*d04, cv01 = dc01+dc03*d13, cv11 = dc11+dc13*d13;
669 //
670 // add MS contribution layer by layer
671 double xCurr = x;
672 double p2Curr = par[2];
673 //
674 // precalculated factors of MS contribution matrix:
675 double ms22t = (1. + par[3]*par[3]);
676 double ms33t = ms22t*ms22t;
677 double p34 = par[3]*par[4];
678 double ms34t = p34*ms22t;
679 double ms44t = p34*p34;
680 //
681 double p2=(1.+ par[3]*par[3])/(par[4]*par[4]);
682 double beta2 = p2/(p2+GetMass()*GetMass());
683 double theta2t = 14.1*14.1/(beta2*p2*1e6) * (1. + par[3]*par[3]);
684 //
685 // account for the MS in the layers between the last measurement and the vertex
686 for (int il=0;il<nMS;il++) {
687 double dfx = xlMS[il] - xCurr;
688 xCurr = xlMS[il];
689 p2Curr += dfx*cnv*par[4]; // p2 at the scattering layer
690 double dxL=xv - xCurr; // distance from scatering layer to vtx
691 double x2rL=cnv*par[4]*dxL, f1L=p2Curr, f2L=f1L+x2rL;
692 if (TMath::Abs(f1L) >= kAlmost1 || TMath::Abs(f2L) >= kAlmost1) {
693 AliInfo(Form("FailMS at step %d of %d: dfx:%e dxL:%e %e %e",il,nMS,dfx,dxL,f1L,f2L));
694 return kFALSE;
695 }
696 double r1L=TMath::Sqrt((1.-f1L)*(1.+f1L)), r2L=TMath::Sqrt((1.-f2L)*(1.+f2L)), dx2rL=dxL/(r1L+r2L);
697 // elements of matrix for propagation from scatering layer to vertex
698 double f02L = 2*dx2rL, f04L = cnv*dxL*dx2rL, f13L/*, f24L = cnv*dxL*/;
699 if (TMath::Abs(x2rL)<0.05) f13L = dxL*r2L+f2L*(f1L+f2L)*dx2rL; // see AliExternalTrackParam::PropagateTo
700 else {
701 double dy2dxL = (f1L+f2L)/(r1L+r2L);
702 f13L = 2*TMath::ASin(0.5*TMath::Sqrt(1+dy2dxL*dy2dxL)*x2rL)/(cnv*par[4]);
703 }
704 // MS contribution matrix:
705 double theta2 = theta2t*TMath::Abs(x2X0MS[il]);
706 double ms22 = theta2*(1.-p2Curr)*(1.+p2Curr)*ms22t;
707 double ms33 = theta2*ms33t;
708 double ms34 = theta2*ms34t;
709 double ms44 = theta2*ms44t;
710 //
711 // add H F MS F^Tr H^Tr to cv
712 cv00 += f02L*f02L*ms22 + f04L*f04L*ms44;
713 cv01 += f04L*f13L*ms34;
714 cv11 += f13L*f13L*ms33;
715 }
716 //
717 // inverse of matrix B
718 double b11 = ers[1]*ers[1] + cv00;
719 double b00 = ers[2]*ers[2] + cv11;
720 double det = b11*b00 - cv01*cv01;
721 if (TMath::Abs(det)<kAlmost0) {
722 AliInfo(Form("Fail on det %e: %e %e %e",det,cv00,cv11,cv01));
723 return kFALSE;
724 }
725 det = 1./det;
726 b00 *= det; b11 *= det;
727 double b01 = -cv01*det;
728 //
729 // elements of matrix DC^Tr * B^-1
730 double dcb00 = b00*dc00+b01*dc10, dcb01 = b01*dc00+b11*dc10;
731 double dcb10 = b00*dc01+b01*dc11, dcb11 = b01*dc01+b11*dc11;
732 double dcb20 = b00*dc02+b01*dc12, dcb21 = b01*dc02+b11*dc12;
733 double dcb30 = b00*dc03+b01*dc13, dcb31 = b01*dc03+b11*dc13;
734 double dcb40 = b00*dc04+b01*dc14, dcb41 = b01*dc04+b11*dc14;
735 //
736 // p_{k|k+1} = p_{k|k} + C_{k|k}*D^Tr_{k+1} B^{-1}_{k+1} ( vtx - D_{k+1}*p_{k|k})
737 par[0] += dcb00*yv + dcb01*zv;
738 par[1] += dcb10*yv + dcb11*zv;
739 par[2] += dcb20*yv + dcb21*zv;
740 par[3] += dcb30*yv + dcb31*zv;
741 par[4] += dcb40*yv + dcb41*zv;
742 //
743 // C_{k|kv} = C_{k|k} - C_{k|k} D^Tr_{k+1} B^{-1}_{k+1} D_{k+1} C^Tr_{k|k})
744 //
745 cov00 = c00 - (dc00*dcb00 + dc10*dcb01);
746 cov01 = c01 - (dc01*dcb00 + dc11*dcb01);
747 cov02 = c02 - (dc02*dcb00 + dc12*dcb01);
748 cov03 = c03 - (dc03*dcb00 + dc13*dcb01);
749 cov04 = c04 - (dc04*dcb00 + dc14*dcb01);
750 //
751 cov11 = c11 - (dc01*dcb10 + dc11*dcb11);
752 cov12 = c12 - (dc02*dcb10 + dc12*dcb11);
753 cov13 = c13 - (dc03*dcb10 + dc13*dcb11);
754 cov14 = c14 - (dc04*dcb10 + dc14*dcb11);
755 //
756 cov22 = c22 - (dc02*dcb20 + dc12*dcb21);
757 cov23 = c23 - (dc03*dcb20 + dc13*dcb21);
758 cov24 = c24 - (dc04*dcb20 + dc14*dcb21);
759 //
760 cov33 = c33 - (dc03*dcb30 + dc13*dcb31);
761 cov34 = c34 - (dc04*dcb30 + dc14*dcb31);
762 //
763 cov44 = c44 - (dc04*dcb40 + dc14*dcb41);
764 //
765 Set(x,alpha,par,covc);
766 if (!Invariant()) {
767 AliInfo(Form("Fail on Invariant, X=%e",GetX()));
768 return kFALSE;
769 }
770 return kTRUE;
771 //
772}