Update for Ds
[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()) {
c168e271 276 if (n>fgkWARN) AliDebug(1,"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){
c168e271 323 if (n>fgkWARN) AliDebug(1,Form("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)) {
c168e271 328 if (n>fgkWARN) AliDebug(1,Form("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)) {
c168e271 333 if (n>fgkWARN) AliDebug(1,Form("fC11=%f\n",sC11));
6c94f330 334 return kFALSE;
a9a2d814 335 }
6c94f330 336 Double_t sC22=GetCovariance()[5];
337 if (sC22<=0 || sC22>1.) {
c168e271 338 if (n>fgkWARN) AliDebug(1,Form("fC22=%f\n",sC22));
6c94f330 339 return kFALSE;
a9a2d814 340 }
6c94f330 341 Double_t sC33=GetCovariance()[9];
342 if (sC33<=0 || sC33>1.) {
c168e271 343 if (n>fgkWARN) AliDebug(1,Form("fC33=%f\n",sC33));
6c94f330 344 return kFALSE;
a9a2d814 345 }
6c94f330 346 Double_t sC44=GetCovariance()[14];
347 if (sC44<=0 /*|| sC44>6e-5*/) {
15fa26ff 348 if (n>fgkWARN) AliDebug(1,Form("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();
c168e271 367 if (n>fgkWARN) AliDebug(1,"Wrong invariant !");
8a84886c 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
367c7d1a 491void AliITStrackV2::CookdEdx(Double_t /*low*/, Double_t /*up*/) {
23efe5f1 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
367c7d1a 525 Double_t sumamp=0,sumweight=0;
526 Double_t weight[4]={1.,1.,0.,0.};
527 if(nc==3) weight[1]=0.5;
528 else if(nc<3) weight[1]=0.;
09cf9d66 529 for (Int_t i=0; i<nc; i++) {
367c7d1a 530 sumamp+= dedx[i]*weight[i];
531 sumweight+=weight[i];
09cf9d66 532 }
533 SetdEdx(sumamp/sumweight);
23efe5f1 534}
c7bafca9 535
8602c008 536//____________________________________________________________________________
537Bool_t AliITStrackV2::
e50912db 538GetPhiZat(Double_t r, Double_t &phi, Double_t &z) const {
539 //------------------------------------------------------------------
540 // This function returns the global cylindrical (phi,z) of the track
541 // position estimated at the radius r.
542 // The track curvature is neglected.
543 //------------------------------------------------------------------
544 Double_t d=GetD(0.,0.);
bc02d571 545 if (TMath::Abs(d) > r) {
546 if (r>1e-1) return kFALSE;
547 r = TMath::Abs(d);
548 }
e50912db 549
550 Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
494a930c 551 if (TMath::Abs(d) > rcurr) return kFALSE;
552 Double_t globXYZcurr[3]; GetXYZ(globXYZcurr);
553 Double_t phicurr=TMath::ATan2(globXYZcurr[1],globXYZcurr[0]);
e50912db 554
2089d338 555 if (GetX()>=0.) {
556 phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr);
557 } else {
558 phi=phicurr+TMath::ASin(d/r)+TMath::ASin(d/rcurr)-TMath::Pi();
559 }
560
494a930c 561 // return a phi in [0,2pi[
562 if (phi<0.) phi+=2.*TMath::Pi();
563 else if (phi>=2.*TMath::Pi()) phi-=2.*TMath::Pi();
c7de5c65 564 z=GetZ()+GetTgl()*(TMath::Sqrt((r-d)*(r+d))-TMath::Sqrt((rcurr-d)*(rcurr+d)));
e50912db 565 return kTRUE;
566}
567//____________________________________________________________________________
568Bool_t AliITStrackV2::
569GetLocalXat(Double_t r,Double_t &xloc) const {
8602c008 570 //------------------------------------------------------------------
e50912db 571 // This function returns the local x of the track
572 // position estimated at the radius r.
8602c008 573 // The track curvature is neglected.
574 //------------------------------------------------------------------
575 Double_t d=GetD(0.,0.);
6a14afad 576 if (TMath::Abs(d) > r) {
577 if (r>1e-1) return kFALSE;
578 r = TMath::Abs(d);
579 }
8602c008 580
e50912db 581 Double_t rcurr=TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
494a930c 582 Double_t globXYZcurr[3]; GetXYZ(globXYZcurr);
583 Double_t phicurr=TMath::ATan2(globXYZcurr[1],globXYZcurr[0]);
2089d338 584 Double_t phi;
585 if (GetX()>=0.) {
586 phi=phicurr+TMath::ASin(d/r)-TMath::ASin(d/rcurr);
587 } else {
588 phi=phicurr+TMath::ASin(d/r)+TMath::ASin(d/rcurr)-TMath::Pi();
589 }
8602c008 590
e50912db 591 xloc=r*(TMath::Cos(phi)*TMath::Cos(GetAlpha())
592 +TMath::Sin(phi)*TMath::Sin(GetAlpha()));
8602c008 593
594 return kTRUE;
595}
63126a46 596
a4354152 597//____________________________________________________________________________
598Bool_t AliITStrackV2::ImproveKalman(Double_t xyz[3],Double_t ers[3], const Double_t* xlMS, const Double_t* x2X0MS, Int_t nMS)
599{
600 // Substitute the state of the track (p_{k|k},C_{k|k}) at the k-th measumerent by its
601 // smoothed value from the k-th measurement + measurement at the vertex.
602 // Account for the MS on nMS layers at x-postions xlMS with x/x0 = x2X0MS
603 // p_{k|kv} = p_{k|k} + C_{k|k}*D^Tr_{k+1} B^{-1}_{k+1} ( vtx - D_{k+1}*p_{k|k})
604 // C_{k|kv} = C_{k|k}*( I - D^Tr_{k+1} B^{-1}_{k+1} D_{k+1} C^Tr_{k|k})
605 //
606 // where D_{k} = H_{k} F_{k} with H being the matrix converting the tracks parameters
607 // to measurements m_{k} = H_{k} p_{k} and F_{k} the matrix propagating the track between the
608 // the point k-1 and k: p_{k|k-1} = F_{k} p_{k-1|k-1}
609 //
610 // 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
611 // at point k+1 (i.e. vertex), and C_{k+1|k} - error matrix extrapolated from k-th measurement to
612 // k+1 (vtx) and accounting for the MS inbetween
613 //
614 // H = {{1,0,0,0,0},{0,1,0,0,0}}
615 //
616 double covc[15], *cori = (double*) GetCovariance(),par[5] = {GetY(),GetZ(),GetSnp(),GetTgl(),GetSigned1Pt()},
617 &c00=cori[0],
618 &c01=cori[1],&c11=cori[2],
619 &c02=cori[3],&c12=cori[4],&c22=cori[5],
620 &c03=cori[6],&c13=cori[7],&c23=cori[8],&c33=cori[9],
621 &c04=cori[10],&c14=cori[11],&c24=cori[12],&c34=cori[13],&c44=cori[14],
622 // for smoothed cov matrix
623 &cov00=covc[0],
624 &cov01=covc[1],&cov11=covc[2],
625 &cov02=covc[3],&cov12=covc[4],&cov22=covc[5],
626 &cov03=covc[6],&cov13=covc[7],&cov23=covc[8],&cov33=covc[9],
627 &cov04=covc[10],&cov14=covc[11],&cov24=covc[12],&cov34=covc[13],&cov44=covc[14];
628 //
629 double x = GetX(), alpha = GetAlpha();
630 // vertex in the track frame
631 double cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
632 double xv = xyz[0]*cs + xyz[1]*sn, yv =-xyz[0]*sn + xyz[1]*cs, zv = xyz[2];
633 double dx = xv - GetX();
634 if (TMath::Abs(dx)<=kAlmost0) return kTRUE;
635 //
636 double cnv=GetBz()*kB2C, x2r=cnv*par[4]*dx, f1=par[2], f2=f1+x2r;
637 if (TMath::Abs(f1) >= kAlmost1 || TMath::Abs(f2) >= kAlmost1) {
638 AliInfo(Form("Fail: %+e %+e",f1,f2));
639 return kFALSE;
640 }
641 double r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2)), dx2r=dx/(r1+r2);
642 // elements of matrix F_{k+1} (1s on diagonal)
643 double f02 = 2*dx2r, f04 = cnv*dx*dx2r, f13/*, f24 = cnv*dx*/;
644 if (TMath::Abs(x2r)<0.05) f13 = dx*r2+f2*(f1+f2)*dx2r; // see AliExternalTrackParam::PropagateTo
645 else {
646 double dy2dx = (f1+f2)/(r1+r2);
647 f13 = 2*TMath::ASin(0.5*TMath::Sqrt(1+dy2dx*dy2dx)*x2r)/(cnv*par[4]);
648 }
649 // elements of matrix D_{k+1} = H_{k+1} * F_{k+1}
650 // double d00 = 1., d11 = 1.;
651 double &d02 = f02, &d04 = f04, &d13 = f13;
652 //
653 // elements of matrix DC = D_{k+1}*C_{kk}^T
654 double dc00 = c00+c02*d02+c04*d04, dc10 = c01+c03*d13;
655 double dc01 = c01+c12*d02+c14*d04, dc11 = c11+c13*d13;
656 double dc02 = c02+c22*d02+c24*d04, dc12 = c12+c23*d13;
657 double dc03 = c03+c23*d02+c34*d04, dc13 = c13+c33*d13;
658 double dc04 = c04+c24*d02+c44*d04, dc14 = c14+c34*d13;
659 //
660 // difference between the vertex and the the track extrapolated to vertex
661 yv -= par[0] + par[2]*d02 + par[4]*d04;
662 zv -= par[1] + par[3]*d13;
663 //
664 // y,z part of the cov.matrix extrapolated to vtx (w/o MS contribution)
665 // C_{k+1,k} = H F_{k+1} C_{k,k} F^Tr_{k+1} H^Tr = D C D^Tr
666 double cv00 = dc00+dc02*d02+dc04*d04, cv01 = dc01+dc03*d13, cv11 = dc11+dc13*d13;
667 //
668 // add MS contribution layer by layer
669 double xCurr = x;
670 double p2Curr = par[2];
671 //
672 // precalculated factors of MS contribution matrix:
673 double ms22t = (1. + par[3]*par[3]);
674 double ms33t = ms22t*ms22t;
675 double p34 = par[3]*par[4];
676 double ms34t = p34*ms22t;
677 double ms44t = p34*p34;
678 //
679 double p2=(1.+ par[3]*par[3])/(par[4]*par[4]);
680 double beta2 = p2/(p2+GetMass()*GetMass());
681 double theta2t = 14.1*14.1/(beta2*p2*1e6) * (1. + par[3]*par[3]);
682 //
683 // account for the MS in the layers between the last measurement and the vertex
684 for (int il=0;il<nMS;il++) {
685 double dfx = xlMS[il] - xCurr;
686 xCurr = xlMS[il];
687 p2Curr += dfx*cnv*par[4]; // p2 at the scattering layer
688 double dxL=xv - xCurr; // distance from scatering layer to vtx
689 double x2rL=cnv*par[4]*dxL, f1L=p2Curr, f2L=f1L+x2rL;
690 if (TMath::Abs(f1L) >= kAlmost1 || TMath::Abs(f2L) >= kAlmost1) {
691 AliInfo(Form("FailMS at step %d of %d: dfx:%e dxL:%e %e %e",il,nMS,dfx,dxL,f1L,f2L));
692 return kFALSE;
693 }
694 double r1L=TMath::Sqrt((1.-f1L)*(1.+f1L)), r2L=TMath::Sqrt((1.-f2L)*(1.+f2L)), dx2rL=dxL/(r1L+r2L);
695 // elements of matrix for propagation from scatering layer to vertex
696 double f02L = 2*dx2rL, f04L = cnv*dxL*dx2rL, f13L/*, f24L = cnv*dxL*/;
697 if (TMath::Abs(x2rL)<0.05) f13L = dxL*r2L+f2L*(f1L+f2L)*dx2rL; // see AliExternalTrackParam::PropagateTo
698 else {
699 double dy2dxL = (f1L+f2L)/(r1L+r2L);
700 f13L = 2*TMath::ASin(0.5*TMath::Sqrt(1+dy2dxL*dy2dxL)*x2rL)/(cnv*par[4]);
701 }
702 // MS contribution matrix:
703 double theta2 = theta2t*TMath::Abs(x2X0MS[il]);
704 double ms22 = theta2*(1.-p2Curr)*(1.+p2Curr)*ms22t;
705 double ms33 = theta2*ms33t;
706 double ms34 = theta2*ms34t;
707 double ms44 = theta2*ms44t;
708 //
709 // add H F MS F^Tr H^Tr to cv
710 cv00 += f02L*f02L*ms22 + f04L*f04L*ms44;
711 cv01 += f04L*f13L*ms34;
712 cv11 += f13L*f13L*ms33;
713 }
714 //
715 // inverse of matrix B
716 double b11 = ers[1]*ers[1] + cv00;
717 double b00 = ers[2]*ers[2] + cv11;
718 double det = b11*b00 - cv01*cv01;
719 if (TMath::Abs(det)<kAlmost0) {
720 AliInfo(Form("Fail on det %e: %e %e %e",det,cv00,cv11,cv01));
721 return kFALSE;
722 }
723 det = 1./det;
724 b00 *= det; b11 *= det;
725 double b01 = -cv01*det;
726 //
727 // elements of matrix DC^Tr * B^-1
728 double dcb00 = b00*dc00+b01*dc10, dcb01 = b01*dc00+b11*dc10;
729 double dcb10 = b00*dc01+b01*dc11, dcb11 = b01*dc01+b11*dc11;
730 double dcb20 = b00*dc02+b01*dc12, dcb21 = b01*dc02+b11*dc12;
731 double dcb30 = b00*dc03+b01*dc13, dcb31 = b01*dc03+b11*dc13;
732 double dcb40 = b00*dc04+b01*dc14, dcb41 = b01*dc04+b11*dc14;
733 //
734 // 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})
735 par[0] += dcb00*yv + dcb01*zv;
736 par[1] += dcb10*yv + dcb11*zv;
737 par[2] += dcb20*yv + dcb21*zv;
738 par[3] += dcb30*yv + dcb31*zv;
739 par[4] += dcb40*yv + dcb41*zv;
740 //
741 // C_{k|kv} = C_{k|k} - C_{k|k} D^Tr_{k+1} B^{-1}_{k+1} D_{k+1} C^Tr_{k|k})
742 //
743 cov00 = c00 - (dc00*dcb00 + dc10*dcb01);
744 cov01 = c01 - (dc01*dcb00 + dc11*dcb01);
745 cov02 = c02 - (dc02*dcb00 + dc12*dcb01);
746 cov03 = c03 - (dc03*dcb00 + dc13*dcb01);
747 cov04 = c04 - (dc04*dcb00 + dc14*dcb01);
748 //
749 cov11 = c11 - (dc01*dcb10 + dc11*dcb11);
750 cov12 = c12 - (dc02*dcb10 + dc12*dcb11);
751 cov13 = c13 - (dc03*dcb10 + dc13*dcb11);
752 cov14 = c14 - (dc04*dcb10 + dc14*dcb11);
753 //
754 cov22 = c22 - (dc02*dcb20 + dc12*dcb21);
755 cov23 = c23 - (dc03*dcb20 + dc13*dcb21);
756 cov24 = c24 - (dc04*dcb20 + dc14*dcb21);
757 //
758 cov33 = c33 - (dc03*dcb30 + dc13*dcb31);
759 cov34 = c34 - (dc04*dcb30 + dc14*dcb31);
760 //
761 cov44 = c44 - (dc04*dcb40 + dc14*dcb41);
762 //
763 Set(x,alpha,par,covc);
764 if (!Invariant()) {
765 AliInfo(Form("Fail on Invariant, X=%e",GetX()));
766 return kFALSE;
767 }
768 return kTRUE;
769 //
770}