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