Add the number of local boards
[u/mrichter/AliRoot.git] / TRD / AliTRDtrack.cxx
CommitLineData
46d29e70 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
0fa7dfa7 16/* $Id$ */
46d29e70 17
a2cb5b3d 18#include <Riostream.h>
fe33d239 19
53c17fbf 20#include <TMath.h>
e1e6896f 21#include <TVector2.h>
46d29e70 22
6c94f330 23#include "AliTracker.h"
53c17fbf 24#include "AliESDtrack.h"
fe33d239 25
46d29e70 26#include "AliTRDgeometry.h"
27#include "AliTRDcluster.h"
28#include "AliTRDtrack.h"
53c17fbf 29#include "AliTRDtracklet.h"
b3a5a838 30
46d29e70 31ClassImp(AliTRDtrack)
32
53c17fbf 33///////////////////////////////////////////////////////////////////////////////
34// //
35// Represents a reconstructed TRD track //
36// Local TRD Kalman track //
37// //
38///////////////////////////////////////////////////////////////////////////////
7ad19338 39
fe33d239 40//_____________________________________________________________________________
41AliTRDtrack::AliTRDtrack()
42 :AliKalmanTrack()
43 ,fSeedLab(-1)
44 ,fdEdx(0)
45 ,fDE(0)
46 ,fStopped(kFALSE)
47 ,fLhElectron(0)
48 ,fNWrong(0)
49 ,fNRotate(0)
50 ,fNCross(0)
51 ,fNExpected(0)
52 ,fNLast(0)
53 ,fNExpectedLast(0)
54 ,fNdedx(0)
55 ,fChi2Last(1e10)
56 ,fBackupTrack(0x0)
f146da82 57{
fe33d239 58 //
59 // AliTRDtrack default constructor
60 //
61
62 for (Int_t i = 0; i < kNplane; i++) {
63 for (Int_t j = 0; j < kNslice; j++) {
64 fdEdxPlane[i][j] = 0.0;
6d45eaef 65 }
f146da82 66 fTimBinPlane[i] = -1;
67 }
fe33d239 68
69 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
70 fIndex[i] = 0;
6c94f330 71 fIndexBackup[i] = 0;
fe33d239 72 fdQdl[i] = 0;
73 }
74
75 for (Int_t i = 0; i < 3; i++) {
76 fBudget[i] = 0;
f146da82 77 }
fe33d239 78
f146da82 79}
6d45eaef 80
46d29e70 81//_____________________________________________________________________________
fe33d239 82AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, Int_t index
83 , const Double_t p[5], const Double_t cov[15]
84 , Double_t x, Double_t alpha)
85 :AliKalmanTrack()
86 ,fSeedLab(-1)
87 ,fdEdx(0)
88 ,fDE(0)
89 ,fStopped(kFALSE)
90 ,fLhElectron(0)
91 ,fNWrong(0)
92 ,fNRotate(0)
93 ,fNCross(0)
94 ,fNExpected(0)
95 ,fNLast(0)
96 ,fNExpectedLast(0)
97 ,fNdedx(0)
98 ,fChi2Last(1e10)
99 ,fBackupTrack(0x0)
15ed8ba1 100{
fe33d239 101 //
102 // The main AliTRDtrack constructor.
103 //
104
105 Double_t cnv = 1.0/(GetBz() * kB2C);
106
107 Double_t pp[5] = { p[0]
108 , p[1]
109 , x*p[4] - p[2]
110 , p[3]
111 , p[4]*cnv };
6c94f330 112
113 Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[5];
114 Double_t c32 = x*cov[13] - cov[8];
fe33d239 115 Double_t c20 = x*cov[10] - cov[3];
116 Double_t c21 = x*cov[11] - cov[4];
117 Double_t c42 = x*cov[14] - cov[12];
6c94f330 118
fe33d239 119 Double_t cc[15] = { cov[0 ]
120 , cov[1 ], cov[2 ]
121 , c20, c21, c22
122 , cov[6 ], cov[7 ], c32, cov[9 ]
123 , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
6c94f330 124
125 Set(x,alpha,pp,cc);
5443e65e 126 SetNumberOfClusters(1);
fe33d239 127 fIndex[0] = index;
5443e65e 128
fe33d239 129 for (Int_t i = 0; i < kNplane; i++) {
130 for (Int_t j = 0; j < kNslice; j++) {
131 fdEdxPlane[i][j] = 0.0;
6d45eaef 132 }
133 fTimBinPlane[i] = -1;
eab5961e 134 }
a2b90f83 135
5443e65e 136 Double_t q = TMath::Abs(c->GetQ());
fe33d239 137 Double_t s = GetSnp();
138 Double_t t = GetTgl();
139 if (s*s < 1) {
140 q *= TMath::Sqrt((1-s*s)/(1+t*t));
141 }
53c17fbf 142
6c94f330 143 fdQdl[0] = q;
fe33d239 144
145 for (UInt_t i = 1; i < kMAXCLUSTERSPERTRACK; i++) {
146 fdQdl[i] = 0;
147 fIndex[i] = 0;
148 fIndexBackup[i] = 0;
15ed8ba1 149 }
fe33d239 150
151 for (Int_t i = 0; i < 3;i++) {
152 fBudget[i] = 0;
153 }
154
46d29e70 155}
156
157//_____________________________________________________________________________
fe33d239 158AliTRDtrack::AliTRDtrack(const AliTRDtrack &t)
159 :AliKalmanTrack(t)
160 ,fSeedLab(t.GetSeedLabel())
161 ,fdEdx(t.fdEdx)
162 ,fDE(t.fDE)
163 ,fStopped(t.fStopped)
164 ,fLhElectron(0)
165 ,fNWrong(t.fNWrong)
166 ,fNRotate(t.fNRotate)
167 ,fNCross(t.fNCross)
168 ,fNExpected(t.fNExpected)
169 ,fNLast(t.fNLast)
170 ,fNExpectedLast(t.fNExpectedLast)
171 ,fNdedx(t.fNdedx)
172 ,fChi2Last(t.fChi2Last)
173 ,fBackupTrack(0x0)
53c17fbf 174{
46d29e70 175 //
176 // Copy constructor.
177 //
fe33d239 178
179 for (Int_t i = 0; i < kNplane; i++) {
180 for (Int_t j = 0; j < kNslice; j++) {
6d45eaef 181 fdEdxPlane[i][j] = t.fdEdxPlane[i][j];
182 }
183 fTimBinPlane[i] = t.fTimBinPlane[i];
184 fTracklets[i] = t.fTracklets[i];
eab5961e 185 }
46d29e70 186
fe33d239 187 Int_t n = t.GetNumberOfClusters();
6c94f330 188 SetNumberOfClusters(n);
fe33d239 189
190 for (Int_t i = 0; i < n; i++) {
191 fIndex[i] = t.fIndex[i];
192 fIndexBackup[i] = t.fIndex[i];
193 fdQdl[i] = t.fdQdl[i];
194 }
195
196 for (UInt_t i = n; i < kMAXCLUSTERSPERTRACK; i++) {
197 fdQdl[i] = 0;
198 fIndex[i] = 0;
199 fIndexBackup[i] = 0;
03b0452e 200 }
15ed8ba1 201
fe33d239 202 for (Int_t i = 0; i < 3;i++) {
203 fBudget[i] = t.fBudget[i];
15ed8ba1 204 }
fe33d239 205
5443e65e 206}
207
208//_____________________________________________________________________________
fe33d239 209AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t /*alpha*/)
210 :AliKalmanTrack(t)
211 ,fSeedLab(-1)
212 ,fdEdx(t.GetPIDsignal())
213 ,fDE(0)
214 ,fStopped(kFALSE)
215 ,fLhElectron(0.0)
216 ,fNWrong(0)
217 ,fNRotate(0)
218 ,fNCross(0)
219 ,fNExpected(0)
220 ,fNLast(0)
221 ,fNExpectedLast(0)
222 ,fNdedx(0)
223 ,fChi2Last(0.0)
224 ,fBackupTrack(0x0)
53c17fbf 225{
5443e65e 226 //
fe33d239 227 // Constructor from AliTPCtrack or AliITStrack
5443e65e 228 //
229
6c94f330 230 SetLabel(t.GetLabel());
fe33d239 231 SetChi2(0.0);
6c94f330 232 SetMass(t.GetMass());
5443e65e 233 SetNumberOfClusters(0);
234
fe33d239 235 for (Int_t i = 0; i < kNplane; i++) {
236 for (Int_t j = 0; j < kNslice; j++) {
6d45eaef 237 fdEdxPlane[i][j] = 0.0;
238 }
eab5961e 239 fTimBinPlane[i] = -1;
240 }
5443e65e 241
fe33d239 242 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
243 fdQdl[i] = 0;
244 fIndex[i] = 0;
245 fIndexBackup[i] = 0;
79e94bf8 246 }
4f1c04d3 247
fe33d239 248 for (Int_t i = 0; i < 3; i++) {
249 fBudget[i] = 0;
250 }
251
79e94bf8 252}
53c17fbf 253
79e94bf8 254//_____________________________________________________________________________
fe33d239 255AliTRDtrack::AliTRDtrack(const AliESDtrack &t)
256 :AliKalmanTrack()
257 ,fSeedLab(-1)
258 ,fdEdx(t.GetTRDsignal())
259 ,fDE(0)
260 ,fStopped(kFALSE)
261 ,fLhElectron(0)
262 ,fNWrong(0)
263 ,fNRotate(0)
264 ,fNCross(0)
265 ,fNExpected(0)
266 ,fNLast(0)
267 ,fNExpectedLast(0)
268 ,fNdedx(0)
269 ,fChi2Last(1e10)
270 ,fBackupTrack(0x0)
53c17fbf 271{
79e94bf8 272 //
273 // Constructor from AliESDtrack
274 //
fe33d239 275
79e94bf8 276 SetLabel(t.GetLabel());
fe33d239 277 SetChi2(0.0);
79e94bf8 278 SetMass(t.GetMass());
1e9bb598 279 SetNumberOfClusters(t.GetTRDclusters(fIndex));
15ed8ba1 280
46e2d86c 281 Int_t ncl = t.GetTRDclusters(fIndexBackup);
fe33d239 282 for (UInt_t i = ncl; i < kMAXCLUSTERSPERTRACK; i++) {
283 fIndexBackup[i] = 0;
284 fIndex[i] = 0;
15ed8ba1 285 }
fe33d239 286
287 for (Int_t i = 0; i < kNplane; i++) {
288 for (Int_t j = 0; j < kNslice; j++) {
6d45eaef 289 fdEdxPlane[i][j] = t.GetTRDsignals(i,j);
290 }
eab5961e 291 fTimBinPlane[i] = t.GetTRDTimBin(i);
292 }
79e94bf8 293
fe33d239 294 const AliExternalTrackParam *par = &t;
295 if (t.GetStatus() & AliESDtrack::kTRDbackup) {
296 par = t.GetOuterParam();
297 if (!par) {
298 AliError("No backup info!");
299 par = &t;
300 }
c5a8e3df 301 }
6c94f330 302 Set(par->GetX(),par->GetAlpha(),par->GetParameter(),par->GetCovariance());
79e94bf8 303
6c94f330 304
fe33d239 305 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
306 fdQdl[i] = 0;
307 }
308
309 for (Int_t i = 0; i < 3; i++) {
310 fBudget[i] = 0;
311 }
5443e65e 312
fe33d239 313 if ((t.GetStatus() & AliESDtrack::kTIME) == 0) {
314 return;
315 }
c630aafd 316
c630aafd 317 StartTimeIntegral();
fe33d239 318 Double_t times[10];
319 t.GetIntegratedTimes(times);
320 SetIntegratedTimes(times);
c630aafd 321 SetIntegratedLength(t.GetIntegratedLength());
322
16d9fbba 323}
324
53c17fbf 325//____________________________________________________________________________
326AliTRDtrack::~AliTRDtrack()
3fad3d32 327{
328 //
53c17fbf 329 // Destructor
330 //
3fad3d32 331
fe33d239 332 if (fBackupTrack) {
333 delete fBackupTrack;
334 }
53c17fbf 335 fBackupTrack = 0;
3fad3d32 336
53c17fbf 337}
338
339//____________________________________________________________________________
53c17fbf 340Float_t AliTRDtrack::StatusForTOF()
7ad19338 341{
53c17fbf 342 //
343 // Defines the status of the TOF extrapolation
344 //
03b0452e 345
fe33d239 346 // Definition of res ????
347 Float_t res = (0.2 + 0.8 * (fN / (fNExpected + 5.0)))
348 * (0.4 + 0.6 * fTracklets[5].GetN() / 20.0);
349 res *= (0.25 + 0.8 * 40.0 / (40.0 + fBudget[2]));
03b0452e 350 return res;
351
fe33d239 352 // This part of the function is never reached ????
353 // What defines these parameters ????
354 Int_t status = 0;
355 if (GetNumberOfClusters() < 20) return 0;
356 if ((fN > 110) &&
357 (fChi2/(Float_t(fN)) < 3)) return 3; // Gold
358 if ((fNLast > 30) &&
359 (fChi2Last/(Float_t(fNLast)) < 3)) return 3; // Gold
360 if ((fNLast > 20) &&
361 (fChi2Last/(Float_t(fNLast)) < 2)) return 3; // Gold
362 if ((fNLast/(fNExpectedLast+3.0) > 0.8) &&
363 (fChi2Last/Float_t(fNLast) < 5) &&
364 (fNLast > 20)) return 2; // Silber
365 if ((fNLast > 5) &&
366 (((fNLast+1.0)/(fNExpectedLast+1.0)) > 0.8) &&
367 (fChi2Last/(fNLast-5.0) < 6)) return 1;
53c17fbf 368
4f1c04d3 369 return status;
53c17fbf 370
4f1c04d3 371}
46d29e70 372
373//_____________________________________________________________________________
53c17fbf 374Int_t AliTRDtrack::Compare(const TObject *o) const
375{
376 //
377 // Compares tracks according to their Y2 or curvature
378 //
46d29e70 379
d49e463b 380 AliTRDtrack *t = (AliTRDtrack *) o;
46d29e70 381
d49e463b 382 Double_t co = TMath::Abs(t->GetC());
383 Double_t c = TMath::Abs(GetC());
46d29e70 384
d49e463b 385 if (c > co) {
386 return 1;
387 }
388 else if (c < co) {
389 return -1;
390 }
46d29e70 391 return 0;
53c17fbf 392
46d29e70 393}
394
395//_____________________________________________________________________________
15ed8ba1 396void AliTRDtrack::CookdEdx(Double_t low, Double_t up)
397{
398 //
d49e463b 399 // Calculates the truncated dE/dx within the "low" and "up" cuts.
15ed8ba1 400 //
401
d49e463b 402 Int_t i = 0;
403
404 // Array to sort the dEdx values according to amplitude
405 Float_t sorted[kMAXCLUSTERSPERTRACK];
15ed8ba1 406
407 // Number of clusters used for dedx
d49e463b 408 Int_t nc = fNdedx;
409
15ed8ba1 410 // Require at least 10 clusters for a dedx measurement
15ed8ba1 411 if (nc < 10) {
4f1c04d3 412 SetdEdx(0);
413 return;
414 }
a819a5f7 415
15ed8ba1 416 // Lower and upper bound
417 Int_t nl = Int_t(low * nc);
418 Int_t nu = Int_t( up * nc);
419
d49e463b 420 // Can fdQdl be negative ????
15ed8ba1 421 for (i = 0; i < nc; i++) {
422 sorted[i] = TMath::Abs(fdQdl[i]);
423 }
424
425 // Sort the dedx values by amplitude
426 Int_t *index = new Int_t[nc];
427 TMath::Sort(nc,sorted,index,kFALSE);
d49e463b 428
15ed8ba1 429 // Sum up the truncated charge between nl and nu
d49e463b 430 Float_t dedx = 0.0;
15ed8ba1 431 for (i = nl; i <= nu; i++) {
432 dedx += sorted[index[i]];
433 }
d49e463b 434 dedx /= (nu - nl + 1.0);
3551db50 435 SetdEdx(dedx);
436
a819a5f7 437}
438
a819a5f7 439//_____________________________________________________________________________
fe33d239 440Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t x0, Double_t rho)
46d29e70 441{
fe33d239 442 //
46d29e70 443 // Propagates a track of particle with mass=pm to a reference plane
444 // defined by x=xk through media of density=rho and radiationLength=x0
fe33d239 445 //
446
447 if (xk == GetX()) {
448 return kTRUE;
449 }
15ed8ba1 450
fe33d239 451 Double_t oldX = GetX();
452 Double_t oldY = GetY();
453 Double_t oldZ = GetZ();
15ed8ba1 454
fe33d239 455 Double_t bz = GetBz();
46d29e70 456
fe33d239 457 if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
458 return kFALSE;
459 }
9c9d2487 460
fe33d239 461 Double_t x = GetX();
462 Double_t y = GetY();
463 Double_t z = GetZ();
464 Double_t d = TMath::Sqrt((x-oldX)*(x-oldX)
465 + (y-oldY)*(y-oldY)
466 + (z-oldZ)*(z-oldZ));
467 if (oldX < xk) {
468 if (IsStartedTimeIntegral()) {
469 Double_t l2 = d;
470 Double_t crv = GetC();
471 if (TMath::Abs(l2*crv) > 0.0001) {
472 // Make correction for curvature if neccesary
473 l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY));
474 l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
475 l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ));
476 }
477 AddTimeStep(l2);
6c94f330 478 }
46d29e70 479 }
15ed8ba1 480
6c94f330 481 Double_t ll = (oldX < xk) ? -d : d;
fe33d239 482 if (!AliExternalTrackParam::CorrectForMaterial(ll*rho/x0,x0,GetMass())) {
483 return kFALSE;
484 }
15ed8ba1 485
fe33d239 486 {
487
488 // Energy losses************************
489 Double_t p2 = (1.0 + GetTgl()*GetTgl()) / (Get1Pt()*Get1Pt());
490 Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
d248a39d 491 if (beta2<1.0e-10 || (5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0) {
fe33d239 492 return kFALSE;
493 }
494
495 Double_t dE = 0.153e-3 / beta2
496 * (log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
497 * d * rho;
498 Float_t budget = d * rho;
499 fBudget[0] += budget;
500
501 /*
502 // Suspicious part - think about it ?
503 Double_t kinE = TMath::Sqrt(p2);
504 if (dE > 0.8*kinE) dE = 0.8 * kinE; //
505 if (dE < 0) dE = 0.0; // Not valid region for Bethe bloch
506 */
507
508 fDE += dE;
509
510 /*
511 // Suspicious ! I.B.
512 Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE)); // Energy loss fluctuation
513 Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
514 fCcc += sigmac2;
515 fCee += fX*fX * sigmac2;
516 */
6c94f330 517
0d5b5c27 518 }
519
6c94f330 520 return kTRUE;
6d45eaef 521
46d29e70 522}
523
46d29e70 524//_____________________________________________________________________________
fe33d239 525Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index
526 , Double_t h01)
46d29e70 527{
fe33d239 528 //
6c94f330 529 // Assignes found cluster to the track and updates track information
fe33d239 530 //
46d29e70 531
b8dc2353 532 Bool_t fNoTilt = kTRUE;
fe33d239 533 if (TMath::Abs(h01) > 0.003) {
534 fNoTilt = kFALSE;
535 }
536
537 // Add angular effect to the error contribution - MI
6c94f330 538 Float_t tangent2 = GetSnp()*GetSnp();
fe33d239 539 if (tangent2 < 0.90000) {
540 tangent2 = tangent2 / (1.0 - tangent2);
b8dc2353 541 }
fe33d239 542 //Float_t errang = tangent2 * 0.04;
b8dc2353 543
fe33d239 544 Double_t p[2] = {c->GetY(), c->GetZ() };
545 //Double_t cov[3] = {c->GetSigmaY2()+errang, 0.0, c->GetSigmaZ2()*100.0 };
546 Double_t sy2 = c->GetSigmaY2() * 4.0;
547 Double_t sz2 = c->GetSigmaZ2() * 4.0;
548 Double_t cov[3] = {sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
46e2d86c 549
fe33d239 550 if (!AliExternalTrackParam::Update(p,cov)) {
551 return kFALSE;
552 }
46d29e70 553
fe33d239 554 Int_t n = GetNumberOfClusters();
555 fIndex[n] = index;
b8dc2353 556 SetNumberOfClusters(n+1);
fd621f36 557
b8dc2353 558 SetChi2(GetChi2()+chisq);
5443e65e 559
6c94f330 560 return kTRUE;
fe33d239 561
562}
53c17fbf 563
46e2d86c 564//_____________________________________________________________________________
fe33d239 565Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, Int_t index
566 , Double_t h01, Int_t /*plane*/)
567{
568 //
46e2d86c 569 // Assignes found cluster to the track and updates track information
fe33d239 570 //
571
46e2d86c 572 Bool_t fNoTilt = kTRUE;
fe33d239 573 if (TMath::Abs(h01) > 0.003) {
574 fNoTilt = kFALSE;
575 }
576
577 // Add angular effect to the error contribution and make correction - MI
6c94f330 578 Double_t tangent2 = GetSnp()*GetSnp();
579 if (tangent2 < 0.90000){
fe33d239 580 tangent2 = tangent2 / (1.0-tangent2);
15ed8ba1 581 }
6c94f330 582 Double_t tangent = TMath::Sqrt(tangent2);
fe33d239 583 if (GetSnp() < 0) {
584 tangent *= -1;
585 }
586
587 //
588 // Is the following still needed ????
589 //
590
6c94f330 591 // Double_t correction = 0*plane;
592 /*
593 Double_t errang = tangent2*0.04; //
594 Double_t errsys =0.025*0.025*20; //systematic error part
15ed8ba1 595
6c94f330 596 Float_t extend =1;
597 if (c->GetNPads()==4) extend=2;
598 */
599 //if (c->GetNPads()==5) extend=3;
600 //if (c->GetNPads()==6) extend=3;
601 //if (c->GetQ()<15) return 1;
15ed8ba1 602
c5a8e3df 603 /*
6c94f330 604 if (corrector!=0){
3c625a9b 605 //if (0){
606 correction = corrector->GetCorrection(plane,c->GetLocalTimeBin(),tangent);
607 if (TMath::Abs(correction)>0){
608 //if we have info
609 errang = corrector->GetSigma(plane,c->GetLocalTimeBin(),tangent);
610 errang *= errang;
611 errang += tangent2*0.04;
612 }
613 }
c5a8e3df 614 */
6c94f330 615 //
616 //Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
617 /*
618 {
619 Double_t dy=c->GetY() - GetY(), dz=c->GetZ() - GetZ();
620 printf("%e %e %e %e\n",dy,dz,padlength/2,h01);
c84a5e9e 621 }
6c94f330 622 */
15ed8ba1 623
fe33d239 624 Double_t p[2] = { c->GetY(), c->GetZ() };
625 //Double_t cov[3]={ (c->GetSigmaY2()+errang+errsys)*extend, 0.0, c->GetSigmaZ2()*10000.0 };
626 Double_t sy2 = c->GetSigmaY2() * 4.0;
627 Double_t sz2 = c->GetSigmaZ2() * 4.0;
628 Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
46e2d86c 629
fe33d239 630 if (!AliExternalTrackParam::Update(p,cov)) {
631 return kFALSE;
632 }
633
634 Int_t n = GetNumberOfClusters();
635 fIndex[n] = index;
46e2d86c 636 SetNumberOfClusters(n+1);
fe33d239 637 SetChi2(GetChi2() + chisq);
638
639 return kTRUE;
46e2d86c 640
53c17fbf 641}
7ad19338 642
fe33d239 643// //_____________________________________________________________________________
644// Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet)
645// {
646// //
647// // Assignes found tracklet to the track and updates track information
648// //
649// // Can this be removed ????
650// //
651//
652// Double_t r00=(tracklet.GetTrackletSigma2()), r01=0., r11= 10000.;
653// r00+=fCyy; r01+=fCzy; r11+=fCzz;
654// //
655// Double_t det=r00*r11 - r01*r01;
656// Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
657// //
658
659// Double_t dy=tracklet.GetY() - fY, dz=tracklet.GetZ() - fZ;
7ad19338 660
15ed8ba1 661
fe33d239 662// Double_t s00 = tracklet.GetTrackletSigma2(); // error pad
663// Double_t s11 = 100000; // error pad-row
664// Float_t h01 = tracklet.GetTilt();
665// //
666// // r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00;
667// r00 = fCyy + fCzz*h01*h01+s00;
668// // r01 = fCzy + fCzz*h01;
669// r01 = fCzy ;
670// r11 = fCzz + s11;
671// det = r00*r11 - r01*r01;
672// // inverse matrix
673// tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
674
675// Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11;
676// Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11;
677// Double_t k20=fCey*r00+fCez*r01, k21=fCey*r01+fCez*r11;
678// Double_t k30=fCty*r00+fCtz*r01, k31=fCty*r01+fCtz*r11;
679// Double_t k40=fCcy*r00+fCcz*r01, k41=fCcy*r01+fCcz*r11;
15ed8ba1 680
fe33d239 681// // K matrix
682// // k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01);
683// // k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01);
684// // k20=fCey*r00+fCez*(r01+h01*r00),k21=fCey*r01+fCez*(r11+h01*r01);
685// // k30=fCty*r00+fCtz*(r01+h01*r00),k31=fCty*r01+fCtz*(r11+h01*r01);
686// // k40=fCcy*r00+fCcz*(r01+h01*r00),k41=fCcy*r01+fCcz*(r11+h01*r01);
687// //
688// //Update measurement
689// Double_t cur=fC + k40*dy + k41*dz, eta=fE + k20*dy + k21*dz;
690// // cur=fC + k40*dy + k41*dz; eta=fE + k20*dy + k21*dz;
691// if (TMath::Abs(cur*fX-eta) >= 0.90000) {
692// //Int_t n=GetNumberOfClusters();
693// // if (n>4) cerr<<n<<" AliTRDtrack warning: Filtering failed !\n";
694// return 0;
695// }
696// // k01+=h01*k00;
697// // k11+=h01*k10;
698// // k21+=h01*k20;
699// // k31+=h01*k30;
700// // k41+=h01*k40;
701
702
703// fY += k00*dy + k01*dz;
704// fZ += k10*dy + k11*dz;
705// fE = eta;
706// fT += k30*dy + k31*dz;
707// fC = cur;
7ad19338 708
6c94f330 709
fe33d239 710// //Update covariance
711// //
712// //
713// Double_t oldyy = fCyy, oldzz = fCzz; //, oldee=fCee, oldcc =fCcc;
714// Double_t oldzy = fCzy, oldey = fCey, oldty=fCty, oldcy =fCcy;
715// Double_t oldez = fCez, oldtz = fCtz, oldcz=fCcz;
716// //Double_t oldte = fCte, oldce = fCce;
717// //Double_t oldct = fCct;
718
719// fCyy-=k00*oldyy+k01*oldzy;
720// fCzy-=k10*oldyy+k11*oldzy;
721// fCey-=k20*oldyy+k21*oldzy;
722// fCty-=k30*oldyy+k31*oldzy;
723// fCcy-=k40*oldyy+k41*oldzy;
724// //
725// fCzz-=k10*oldzy+k11*oldzz;
726// fCez-=k20*oldzy+k21*oldzz;
727// fCtz-=k30*oldzy+k31*oldzz;
728// fCcz-=k40*oldzy+k41*oldzz;
729// //
730// fCee-=k20*oldey+k21*oldez;
731// fCte-=k30*oldey+k31*oldez;
732// fCce-=k40*oldey+k41*oldez;
733// //
734// fCtt-=k30*oldty+k31*oldtz;
735// fCct-=k40*oldty+k41*oldtz;
736// //
737// fCcc-=k40*oldcy+k41*oldcz;
738// //
15ed8ba1 739
fe33d239 740// //Int_t n=GetNumberOfClusters();
741// //fIndex[n]=index;
742// //SetNumberOfClusters(n+1);
6c94f330 743
fe33d239 744// //SetChi2(GetChi2()+chisq);
745// // cerr<<"in update: fIndex["<<fN<<"] = "<<index<<endl;
7ad19338 746
fe33d239 747// return 1;
7ad19338 748
fe33d239 749// }
c84a5e9e 750
46d29e70 751//_____________________________________________________________________________
6c94f330 752Bool_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
46d29e70 753{
fe33d239 754 //
6c94f330 755 // Rotates track parameters in R*phi plane
756 // if absolute rotation alpha is in global system
757 // otherwise alpha rotation is relative to the current rotation angle
fe33d239 758 //
759
3fad3d32 760 if (absolute) {
6c94f330 761 alpha -= GetAlpha();
3fad3d32 762 }
763 else{
764 fNRotate++;
765 }
46d29e70 766
6c94f330 767 return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
fe33d239 768
53c17fbf 769}
7ad19338 770
46d29e70 771//_____________________________________________________________________________
fd621f36 772Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
46d29e70 773{
53c17fbf 774 //
775 // Returns the track chi2
776 //
777
fe33d239 778 Double_t p[2] = { c->GetY(), c->GetZ() };
779 Double_t sy2 = c->GetSigmaY2() * 4.0;
780 Double_t sz2 = c->GetSigmaZ2() * 4.0;
781 Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
6c94f330 782
783 return AliExternalTrackParam::GetPredictedChi2(p,cov);
15ed8ba1 784
fe33d239 785 //
786 // Can the following be removed ????
787 //
6c94f330 788 /*
789 Bool_t fNoTilt = kTRUE;
790 if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE;
15ed8ba1 791
6c94f330 792 return (c->GetY() - GetY())*(c->GetY() - GetY())/c->GetSigmaY2();
793 */
15ed8ba1 794
6c94f330 795 /*
796 Double_t chi2, dy, r00, r01, r11;
b8dc2353 797
6c94f330 798 if(fNoTilt) {
799 dy=c->GetY() - fY;
800 r00=c->GetSigmaY2();
801 chi2 = (dy*dy)/r00;
46d29e70 802 }
b8dc2353 803 else {
6c94f330 804 Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12);
805 //
806 r00=c->GetSigmaY2(); r01=0.; r11=c->GetSigmaZ2();
807 r00+=fCyy; r01+=fCzy; r11+=fCzz;
b8dc2353 808
6c94f330 809 Double_t det=r00*r11 - r01*r01;
b8dc2353 810 if (TMath::Abs(det) < 1.e-10) {
6c94f330 811 Int_t n=GetNumberOfClusters();
812 if (n>4) cerr<<n<<" AliTRDtrack warning: Singular matrix !\n";
b8dc2353 813 return 1e10;
814 }
6c94f330 815 Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
816 Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
4f1c04d3 817 Double_t tiltdz = dz;
6c94f330 818 if (TMath::Abs(tiltdz)>padlength/2.) {
819 tiltdz = TMath::Sign(padlength/2,dz);
4f1c04d3 820 }
6c94f330 821 // dy=dy+h01*dz;
822 dy=dy+h01*tiltdz;
a819a5f7 823
6c94f330 824 chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
b8dc2353 825 }
53c17fbf 826
b8dc2353 827 return chi2;
6c94f330 828 */
fe33d239 829
830}
7ad19338 831
53c17fbf 832//_____________________________________________________________________________
16d9fbba 833void AliTRDtrack::MakeBackupTrack()
834{
835 //
53c17fbf 836 // Creates a backup track
16d9fbba 837 //
53c17fbf 838
fe33d239 839 if (fBackupTrack) {
840 delete fBackupTrack;
841 }
16d9fbba 842 fBackupTrack = new AliTRDtrack(*this);
4f1c04d3 843
844}
845
53c17fbf 846//_____________________________________________________________________________
847Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
848{
4f1c04d3 849 //
fe33d239 850 // Find a prolongation at given x
851 // Return 0 if it does not exist
852 //
853
854 Double_t bz = GetBz();
15ed8ba1 855
fe33d239 856 if (!AliExternalTrackParam::GetYAt(xk,bz,y)) {
857 return 0;
858 }
859 if (!AliExternalTrackParam::GetZAt(xk,bz,z)) {
860 return 0;
861 }
4f1c04d3 862
6c94f330 863 return 1;
fe33d239 864
16d9fbba 865}
3fad3d32 866
53c17fbf 867//_____________________________________________________________________________
fe33d239 868Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
3fad3d32 869{
870 //
6c94f330 871 // Propagate track to given x position
fe33d239 872 // Works inside of the 20 degree segmentation (local cooordinate frame for TRD , TPC, TOF)
3fad3d32 873 //
fe33d239 874 // Material budget from geo manager
6c94f330 875 //
fe33d239 876
877 Double_t xyz0[3];
878 Double_t xyz1[3];
879 Double_t y;
880 Double_t z;
881
882 const Double_t kAlphac = TMath::Pi()/9.0;
15ed8ba1 883 const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
fe33d239 884
885 // Critical alpha - cross sector indication
886 Double_t dir = (GetX()>xr) ? -1.0 : 1.0;
887
888 // Direction +-
889 for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) {
890
6c94f330 891 GetXYZ(xyz0);
3fad3d32 892 GetProlongation(x,y,z);
fe33d239 893 xyz1[0] = x * TMath::Cos(GetAlpha()) + y * TMath::Sin(GetAlpha());
894 xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha());
3fad3d32 895 xyz1[2] = z;
896 Double_t param[7];
897 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
fe33d239 898
899 if ((param[0] > 0) &&
900 (param[1] > 0)) {
901 PropagateTo(x,param[1],param[0]);
902 }
903
904 if (GetY() > GetX()*kTalphac) {
53c17fbf 905 Rotate(-kAlphac);
3fad3d32 906 }
fe33d239 907 if (GetY() < -GetX()*kTalphac) {
908 Rotate( kAlphac);
3fad3d32 909 }
fe33d239 910
3fad3d32 911 }
fe33d239 912
3fad3d32 913 PropagateTo(xr);
53c17fbf 914
3fad3d32 915 return 0;
3fad3d32 916
53c17fbf 917}
3fad3d32 918
53c17fbf 919//_____________________________________________________________________________
6c94f330 920Int_t AliTRDtrack::PropagateToR(Double_t r,Double_t step)
3fad3d32 921{
922 //
fe33d239 923 // Propagate track to the radial position
924 // Rotation always connected to the last track position
3fad3d32 925 //
fe33d239 926
927 Double_t xyz0[3];
928 Double_t xyz1[3];
929 Double_t y;
930 Double_t z;
931
6c94f330 932 Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
fe33d239 933 // Direction +-
934 Double_t dir = (radius>r) ? -1.0 : 1.0;
935
936 for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
937
6c94f330 938 GetXYZ(xyz0);
3fad3d32 939 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
940 Rotate(alpha,kTRUE);
6c94f330 941 GetXYZ(xyz0);
3fad3d32 942 GetProlongation(x,y,z);
fe33d239 943 xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
944 xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
3fad3d32 945 xyz1[2] = z;
946 Double_t param[7];
947 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
fe33d239 948 if (param[1] <= 0) {
949 param[1] = 100000000;
950 }
3fad3d32 951 PropagateTo(x,param[1],param[0]);
fe33d239 952
3fad3d32 953 }
fe33d239 954
6c94f330 955 GetXYZ(xyz0);
3fad3d32 956 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
957 Rotate(alpha,kTRUE);
6c94f330 958 GetXYZ(xyz0);
3fad3d32 959 GetProlongation(r,y,z);
fe33d239 960 xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
961 xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
3fad3d32 962 xyz1[2] = z;
963 Double_t param[7];
964 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
fe33d239 965
966 if (param[1] <= 0) {
967 param[1] = 100000000;
968 }
3fad3d32 969 PropagateTo(r,param[1],param[0]);
53c17fbf 970
3fad3d32 971 return 0;
53c17fbf 972
973}
974
975//_____________________________________________________________________________
4e28c495 976Int_t AliTRDtrack::GetSector() const
53c17fbf 977{
978 //
979 // Return the current sector
980 //
981
fe33d239 982 return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha())
983 % AliTRDgeometry::kNsect;
53c17fbf 984
3fad3d32 985}
986
53c17fbf 987//_____________________________________________________________________________
4e28c495 988void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
53c17fbf 989{
990 //
991 // The sampled energy loss
992 //
993
994 Double_t s = GetSnp();
995 Double_t t = GetTgl();
fe33d239 996 q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
53c17fbf 997 fdQdl[i] = q;
998
999}
1000
fe33d239 1001//_____________________________________________________________________________
4e28c495 1002void AliTRDtrack::SetSampledEdx(Float_t q)
53c17fbf 1003{
1004 //
1005 // The sampled energy loss
1006 //
1007
1008 Double_t s = GetSnp();
1009 Double_t t = GetTgl();
fe33d239 1010 q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
53c17fbf 1011 fdQdl[fNdedx] = q;
1012 fNdedx++;
1013
1014}
1015
fe33d239 1016//_____________________________________________________________________________
1017Double_t AliTRDtrack::GetBz() const
1018{
15ed8ba1 1019 //
fe33d239 1020 // Returns Bz component of the magnetic field (kG)
15ed8ba1 1021 //
53c17fbf 1022
fe33d239 1023 if (AliTracker::UniformField()) {
1024 return AliTracker::GetBz();
1025 }
1026 Double_t r[3];
1027 GetXYZ(r);
1028
1029 return AliTracker::GetBz(r);
53c17fbf 1030
fe33d239 1031}