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