]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDtrack.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removal of the dependence on TPC (Haavard)
[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 * *
43bfc8af 7 * Permission to use , copy, modify and distribute this software and its *
46d29e70 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
43bfc8af 31// A. Bercuci - used for PID calculations
32#include "AliTRDcalibDB.h"
720a0a16 33#include "Cal/AliTRDCalPID.h"
43bfc8af 34
46d29e70 35ClassImp(AliTRDtrack)
36
53c17fbf 37///////////////////////////////////////////////////////////////////////////////
38// //
39// Represents a reconstructed TRD track //
40// Local TRD Kalman track //
41// //
42///////////////////////////////////////////////////////////////////////////////
7ad19338 43
fe33d239 44//_____________________________________________________________________________
45AliTRDtrack::AliTRDtrack()
46 :AliKalmanTrack()
47 ,fSeedLab(-1)
48 ,fdEdx(0)
49 ,fDE(0)
43bfc8af 50 ,fClusterOwner(kFALSE) // A.Bercuci
fe33d239 51 ,fStopped(kFALSE)
52 ,fLhElectron(0)
53 ,fNWrong(0)
54 ,fNRotate(0)
55 ,fNCross(0)
56 ,fNExpected(0)
57 ,fNLast(0)
58 ,fNExpectedLast(0)
59 ,fNdedx(0)
60 ,fChi2Last(1e10)
61 ,fBackupTrack(0x0)
f146da82 62{
fe33d239 63 //
64 // AliTRDtrack default constructor
65 //
66
67 for (Int_t i = 0; i < kNplane; i++) {
68 for (Int_t j = 0; j < kNslice; j++) {
69 fdEdxPlane[i][j] = 0.0;
6d45eaef 70 }
f146da82 71 fTimBinPlane[i] = -1;
43bfc8af 72 // A.Bercuci additions
73 fMom[i] = -1.;
74 fSnp[i] = 0.;
75 fTgl[i] = 0.;
76 }
fe33d239 77
78 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
79 fIndex[i] = 0;
6c94f330 80 fIndexBackup[i] = 0;
fe33d239 81 fdQdl[i] = 0;
43bfc8af 82 //A.Bercuci additions
83 fClusters[i] = 0x0;
84 }
fe33d239 85
86 for (Int_t i = 0; i < 3; i++) {
87 fBudget[i] = 0;
f146da82 88 }
fe33d239 89
f146da82 90}
6d45eaef 91
46d29e70 92//_____________________________________________________________________________
43bfc8af 93AliTRDtrack::AliTRDtrack(AliTRDcluster *c, Int_t index
fe33d239 94 , const Double_t p[5], const Double_t cov[15]
95 , Double_t x, Double_t alpha)
96 :AliKalmanTrack()
97 ,fSeedLab(-1)
98 ,fdEdx(0)
99 ,fDE(0)
43bfc8af 100 ,fClusterOwner(kFALSE) // A.Bercuci
fe33d239 101 ,fStopped(kFALSE)
102 ,fLhElectron(0)
103 ,fNWrong(0)
104 ,fNRotate(0)
105 ,fNCross(0)
106 ,fNExpected(0)
107 ,fNLast(0)
108 ,fNExpectedLast(0)
109 ,fNdedx(0)
110 ,fChi2Last(1e10)
43bfc8af 111 ,fBackupTrack(0x0)
15ed8ba1 112{
fe33d239 113 //
114 // The main AliTRDtrack constructor.
115 //
116
117 Double_t cnv = 1.0/(GetBz() * kB2C);
118
119 Double_t pp[5] = { p[0]
120 , p[1]
121 , x*p[4] - p[2]
122 , p[3]
123 , p[4]*cnv };
6c94f330 124
125 Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[5];
126 Double_t c32 = x*cov[13] - cov[8];
fe33d239 127 Double_t c20 = x*cov[10] - cov[3];
128 Double_t c21 = x*cov[11] - cov[4];
129 Double_t c42 = x*cov[14] - cov[12];
6c94f330 130
fe33d239 131 Double_t cc[15] = { cov[0 ]
132 , cov[1 ], cov[2 ]
133 , c20, c21, c22
134 , cov[6 ], cov[7 ], c32, cov[9 ]
135 , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
6c94f330 136
137 Set(x,alpha,pp,cc);
5443e65e 138 SetNumberOfClusters(1);
fe33d239 139 fIndex[0] = index;
43bfc8af 140 fClusters[0] = c; // A.Bercuci additions
5443e65e 141
fe33d239 142 for (Int_t i = 0; i < kNplane; i++) {
143 for (Int_t j = 0; j < kNslice; j++) {
144 fdEdxPlane[i][j] = 0.0;
6d45eaef 145 }
146 fTimBinPlane[i] = -1;
43bfc8af 147 // A.Bercuci additions
148 fMom[i] = -1.;
149 fSnp[i] = 0.;
150 fTgl[i] = 0.;
eab5961e 151 }
a2b90f83 152
5443e65e 153 Double_t q = TMath::Abs(c->GetQ());
fe33d239 154 Double_t s = GetSnp();
155 Double_t t = GetTgl();
156 if (s*s < 1) {
157 q *= TMath::Sqrt((1-s*s)/(1+t*t));
158 }
53c17fbf 159
43bfc8af 160 fdQdl[0] = q;
fe33d239 161 for (UInt_t i = 1; i < kMAXCLUSTERSPERTRACK; i++) {
162 fdQdl[i] = 0;
163 fIndex[i] = 0;
164 fIndexBackup[i] = 0;
43bfc8af 165 // A.Bercuci additions
166 fClusters[i] = 0x0;
15ed8ba1 167 }
fe33d239 168
169 for (Int_t i = 0; i < 3;i++) {
170 fBudget[i] = 0;
171 }
172
46d29e70 173}
174
175//_____________________________________________________________________________
43bfc8af 176AliTRDtrack::AliTRDtrack(const AliTRDtrack &t/*, const Bool_t owner*/)
fe33d239 177 :AliKalmanTrack(t)
178 ,fSeedLab(t.GetSeedLabel())
179 ,fdEdx(t.fdEdx)
180 ,fDE(t.fDE)
43bfc8af 181 ,fClusterOwner(kTRUE) // A.Bercuci
fe33d239 182 ,fStopped(t.fStopped)
183 ,fLhElectron(0)
184 ,fNWrong(t.fNWrong)
185 ,fNRotate(t.fNRotate)
186 ,fNCross(t.fNCross)
187 ,fNExpected(t.fNExpected)
188 ,fNLast(t.fNLast)
189 ,fNExpectedLast(t.fNExpectedLast)
190 ,fNdedx(t.fNdedx)
191 ,fChi2Last(t.fChi2Last)
43bfc8af 192 ,fBackupTrack(0x0)
53c17fbf 193{
46d29e70 194 //
195 // Copy constructor.
196 //
fe33d239 197
198 for (Int_t i = 0; i < kNplane; i++) {
199 for (Int_t j = 0; j < kNslice; j++) {
6d45eaef 200 fdEdxPlane[i][j] = t.fdEdxPlane[i][j];
201 }
202 fTimBinPlane[i] = t.fTimBinPlane[i];
203 fTracklets[i] = t.fTracklets[i];
43bfc8af 204 // A.Bercuci additions
205 fMom[i] = t.fMom[i];
206 fSnp[i] = t.fSnp[i];
207 fTgl[i] = t.fTgl[i];
eab5961e 208 }
46d29e70 209
fe33d239 210 Int_t n = t.GetNumberOfClusters();
6c94f330 211 SetNumberOfClusters(n);
fe33d239 212
213 for (Int_t i = 0; i < n; i++) {
214 fIndex[i] = t.fIndex[i];
215 fIndexBackup[i] = t.fIndex[i];
216 fdQdl[i] = t.fdQdl[i];
43bfc8af 217 // A.Bercuci additions
218 if(fClusterOwner && t.fClusters[i]) fClusters[i] = new AliTRDcluster(*(t.fClusters[i]));
219 else fClusters[i] = t.fClusters[i];
220 }
fe33d239 221
222 for (UInt_t i = n; i < kMAXCLUSTERSPERTRACK; i++) {
223 fdQdl[i] = 0;
224 fIndex[i] = 0;
225 fIndexBackup[i] = 0;
43bfc8af 226 // A.Bercuci additions
227 fClusters[i] = 0x0;
03b0452e 228 }
15ed8ba1 229
fe33d239 230 for (Int_t i = 0; i < 3;i++) {
231 fBudget[i] = t.fBudget[i];
15ed8ba1 232 }
fe33d239 233
5443e65e 234}
235
236//_____________________________________________________________________________
fe33d239 237AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t /*alpha*/)
238 :AliKalmanTrack(t)
239 ,fSeedLab(-1)
240 ,fdEdx(t.GetPIDsignal())
241 ,fDE(0)
43bfc8af 242 ,fClusterOwner(kFALSE) // A.Bercuci
fe33d239 243 ,fStopped(kFALSE)
244 ,fLhElectron(0.0)
245 ,fNWrong(0)
246 ,fNRotate(0)
247 ,fNCross(0)
248 ,fNExpected(0)
249 ,fNLast(0)
250 ,fNExpectedLast(0)
251 ,fNdedx(0)
252 ,fChi2Last(0.0)
253 ,fBackupTrack(0x0)
53c17fbf 254{
5443e65e 255 //
fe33d239 256 // Constructor from AliTPCtrack or AliITStrack
5443e65e 257 //
258
6c94f330 259 SetLabel(t.GetLabel());
fe33d239 260 SetChi2(0.0);
6c94f330 261 SetMass(t.GetMass());
5443e65e 262 SetNumberOfClusters(0);
263
fe33d239 264 for (Int_t i = 0; i < kNplane; i++) {
265 for (Int_t j = 0; j < kNslice; j++) {
6d45eaef 266 fdEdxPlane[i][j] = 0.0;
267 }
eab5961e 268 fTimBinPlane[i] = -1;
43bfc8af 269 // A.Bercuci additions
270 fMom[i] = -1.;
271 fSnp[i] = 0.;
272 fTgl[i] = 0.;
eab5961e 273 }
5443e65e 274
fe33d239 275 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
276 fdQdl[i] = 0;
277 fIndex[i] = 0;
278 fIndexBackup[i] = 0;
43bfc8af 279 // A.Bercuci additions
280 fClusters[i] = 0x0;
79e94bf8 281 }
4f1c04d3 282
fe33d239 283 for (Int_t i = 0; i < 3; i++) {
284 fBudget[i] = 0;
285 }
286
79e94bf8 287}
53c17fbf 288
79e94bf8 289//_____________________________________________________________________________
fe33d239 290AliTRDtrack::AliTRDtrack(const AliESDtrack &t)
291 :AliKalmanTrack()
292 ,fSeedLab(-1)
293 ,fdEdx(t.GetTRDsignal())
294 ,fDE(0)
43bfc8af 295 ,fClusterOwner(kFALSE) // A.Bercuci
fe33d239 296 ,fStopped(kFALSE)
297 ,fLhElectron(0)
298 ,fNWrong(0)
299 ,fNRotate(0)
300 ,fNCross(0)
301 ,fNExpected(0)
302 ,fNLast(0)
303 ,fNExpectedLast(0)
304 ,fNdedx(0)
305 ,fChi2Last(1e10)
306 ,fBackupTrack(0x0)
53c17fbf 307{
79e94bf8 308 //
309 // Constructor from AliESDtrack
310 //
fe33d239 311
79e94bf8 312 SetLabel(t.GetLabel());
fe33d239 313 SetChi2(0.0);
79e94bf8 314 SetMass(t.GetMass());
1e9bb598 315 SetNumberOfClusters(t.GetTRDclusters(fIndex));
15ed8ba1 316
46e2d86c 317 Int_t ncl = t.GetTRDclusters(fIndexBackup);
fe33d239 318 for (UInt_t i = ncl; i < kMAXCLUSTERSPERTRACK; i++) {
319 fIndexBackup[i] = 0;
320 fIndex[i] = 0;
15ed8ba1 321 }
fe33d239 322
323 for (Int_t i = 0; i < kNplane; i++) {
324 for (Int_t j = 0; j < kNslice; j++) {
6d45eaef 325 fdEdxPlane[i][j] = t.GetTRDsignals(i,j);
326 }
eab5961e 327 fTimBinPlane[i] = t.GetTRDTimBin(i);
43bfc8af 328 // A.Bercuci additions
329 fMom[i] = -1.;
330 fSnp[i] = 0.;
331 fTgl[i] = 0.;
eab5961e 332 }
79e94bf8 333
fe33d239 334 const AliExternalTrackParam *par = &t;
335 if (t.GetStatus() & AliESDtrack::kTRDbackup) {
336 par = t.GetOuterParam();
337 if (!par) {
338 AliError("No backup info!");
339 par = &t;
340 }
c5a8e3df 341 }
6c94f330 342 Set(par->GetX(),par->GetAlpha(),par->GetParameter(),par->GetCovariance());
79e94bf8 343
6c94f330 344
fe33d239 345 for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
43bfc8af 346 fdQdl[i] = 0;
347 // A.Bercuci additions
348 fClusters[i] = 0x0;
349 }
fe33d239 350
351 for (Int_t i = 0; i < 3; i++) {
352 fBudget[i] = 0;
353 }
5443e65e 354
fe33d239 355 if ((t.GetStatus() & AliESDtrack::kTIME) == 0) {
356 return;
357 }
c630aafd 358
c630aafd 359 StartTimeIntegral();
fe33d239 360 Double_t times[10];
361 t.GetIntegratedTimes(times);
362 SetIntegratedTimes(times);
c630aafd 363 SetIntegratedLength(t.GetIntegratedLength());
364
16d9fbba 365}
366
53c17fbf 367//____________________________________________________________________________
368AliTRDtrack::~AliTRDtrack()
3fad3d32 369{
370 //
53c17fbf 371 // Destructor
372 //
3fad3d32 373
fe33d239 374 if (fBackupTrack) {
375 delete fBackupTrack;
376 }
43bfc8af 377 fBackupTrack = 0x0;
378 if (fClusterOwner){
379 UInt_t icluster=0;
380 while(icluster<kMAXCLUSTERSPERTRACK && fClusters[icluster]){
381 delete fClusters[icluster];
382 fClusters[icluster] = 0x0;
383 icluster++;
384 }
385 }
3fad3d32 386
53c17fbf 387}
388
53c17fbf 389//____________________________________________________________________________
390Float_t AliTRDtrack::StatusForTOF()
7ad19338 391{
53c17fbf 392 //
393 // Defines the status of the TOF extrapolation
394 //
03b0452e 395
fe33d239 396 // Definition of res ????
397 Float_t res = (0.2 + 0.8 * (fN / (fNExpected + 5.0)))
398 * (0.4 + 0.6 * fTracklets[5].GetN() / 20.0);
399 res *= (0.25 + 0.8 * 40.0 / (40.0 + fBudget[2]));
03b0452e 400 return res;
401
fe33d239 402 // This part of the function is never reached ????
403 // What defines these parameters ????
404 Int_t status = 0;
405 if (GetNumberOfClusters() < 20) return 0;
406 if ((fN > 110) &&
407 (fChi2/(Float_t(fN)) < 3)) return 3; // Gold
408 if ((fNLast > 30) &&
409 (fChi2Last/(Float_t(fNLast)) < 3)) return 3; // Gold
410 if ((fNLast > 20) &&
411 (fChi2Last/(Float_t(fNLast)) < 2)) return 3; // Gold
412 if ((fNLast/(fNExpectedLast+3.0) > 0.8) &&
413 (fChi2Last/Float_t(fNLast) < 5) &&
414 (fNLast > 20)) return 2; // Silber
415 if ((fNLast > 5) &&
416 (((fNLast+1.0)/(fNExpectedLast+1.0)) > 0.8) &&
417 (fChi2Last/(fNLast-5.0) < 6)) return 1;
53c17fbf 418
4f1c04d3 419 return status;
53c17fbf 420
4f1c04d3 421}
46d29e70 422
423//_____________________________________________________________________________
53c17fbf 424Int_t AliTRDtrack::Compare(const TObject *o) const
425{
426 //
427 // Compares tracks according to their Y2 or curvature
428 //
46d29e70 429
d49e463b 430 AliTRDtrack *t = (AliTRDtrack *) o;
46d29e70 431
d49e463b 432 Double_t co = TMath::Abs(t->GetC());
433 Double_t c = TMath::Abs(GetC());
46d29e70 434
d49e463b 435 if (c > co) {
436 return 1;
437 }
438 else if (c < co) {
439 return -1;
440 }
46d29e70 441 return 0;
53c17fbf 442
46d29e70 443}
444
a819a5f7 445//_____________________________________________________________________________
43bfc8af 446void AliTRDtrack::CookdEdx(Double_t low, Double_t up)
15ed8ba1 447{
448 //
d49e463b 449 // Calculates the truncated dE/dx within the "low" and "up" cuts.
15ed8ba1 450 //
451
d49e463b 452 // Array to sort the dEdx values according to amplitude
43bfc8af 453
454 Float_t sorted[kMAXCLUSTERSPERTRACK];
455 fdEdx = 0.;
456
15ed8ba1 457 // Require at least 10 clusters for a dedx measurement
c6011b06 458 if (fNdedx < 10) return;
15ed8ba1 459
43bfc8af 460
d49e463b 461 // Can fdQdl be negative ????
43bfc8af 462 for (Int_t i = 0; i < fNdedx; i++) sorted[i] = TMath::Abs(fdQdl[i]);
15ed8ba1 463 // Sort the dedx values by amplitude
c6011b06 464 Int_t *index = new Int_t[fNdedx];
465 TMath::Sort(fNdedx, sorted, index, kFALSE);
466
43bfc8af 467 // Sum up the truncated charge between lower and upper bounds
c6011b06 468 Int_t nl = Int_t(low * fNdedx);
469 Int_t nu = Int_t( up * fNdedx);
43bfc8af 470 for (Int_t i = nl; i <= nu; i++) fdEdx += sorted[index[i]];
471 fdEdx /= (nu - nl + 1.0);
472
473 delete[] index;
474}
475
476//_____________________________________________________________________________
477void AliTRDtrack::CookdEdxTimBin()
478{
479 //
480 // Set fdEdxPlane and fTimBinPlane and also get the
481 // Time bin for Max. Cluster
482 //
483 // Authors:
484 // Prashant Shukla (shukla@physi.uni-heidelberg.de)
485 // Alexandru Bercuci (A.Bercuci@gsi.de)
486
487
488 Double_t maxcharge[AliESDtrack::kNPlane]; // max charge in chamber
489 // number of clusters attached to track per chamber and slice
490 Int_t nCluster[AliESDtrack::kNPlane][AliESDtrack::kNSlice];
491 //number of time bins in chamber
492 Int_t ntb = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
493 Int_t plane; // plane of current cluster
494 Int_t tb; // time bin of current cluster
495 Int_t slice; // curent slice
496 AliTRDcluster *cluster = 0x0; // pointer to current cluster
497
498 // Reset class and local contors/variables
499 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
500 fTimBinPlane[iPlane] = -1;
501 maxcharge[iPlane] = 0.;
502 for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
503 fdEdxPlane[iPlane][iSlice] = 0.;
504 nCluster[iPlane][iSlice] = 0;
505 }
15ed8ba1 506 }
c6011b06 507
43bfc8af 508 // start looping over clusters attached to this track
509 for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) {
510 cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]);
511 if(!cluster) continue;
512
513 // Read info from current cluster
514 plane = AliTRDgeometry::GetPlane(cluster->GetDetector());
515 if (plane < 0 || plane >= AliESDtrack::kNPlane) {
516 AliError(Form("Wrong plane %d", plane));
517 continue;
518 }
519
520 tb = cluster->GetLocalTimeBin();
521 if(tb == 0 || tb >= ntb){
522 AliWarning(Form("time bin 0 or > %d in cluster %d", ntb, iClus));
523 AliInfo(Form("dQ/dl %f", fdQdl[iClus]));
524 continue;
525 }
526
527 slice = tb * AliESDtrack::kNSlice / ntb;
528
529 fdEdxPlane[plane][slice] += fdQdl[iClus];
530 if(fdQdl[iClus] > maxcharge[plane]) {
531 maxcharge[plane] = fdQdl[iClus];
532 fTimBinPlane[plane] = tb;
533 }
534 nCluster[plane][slice]++;
535 } // End of loop over cluster
536
537
538
539 // Normalize fdEdxPlane to number of clusters and set track segments
540 for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
541 for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
542 if (nCluster[iPlane][iSlice]) fdEdxPlane[iPlane][iSlice] /= nCluster[iPlane][iSlice];
543 }
544 }
545}
546
547
548//_____________________________________________________________________________
549void AliTRDtrack::SetTrackSegmentDirMom(const Int_t plane)
550{
551 if(plane<0 || plane>= kNplane){
552 AliError(Form("Trying to access out of range plane (%d)", plane));
553 return;
554 }
555
556 fSnp[plane] = GetSnp();
557 fTgl[plane] = GetTgl();
558 Double_t p[3]; GetPxPyPz(p);
559 fMom[plane] = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
560}
561
562//_____________________________________________________________________________
563Float_t AliTRDtrack::GetTrackLengthPlane(Int_t plane) const
564{
565 if(plane < 0 || plane >= kNplane) return 0.;
566 return (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/TMath::Sqrt((1.
567- fSnp[plane]*fSnp[plane]) / (1. + fTgl[plane]*fTgl[plane]));
568
569}
570
571
572//_____________________________________________________________________________
573Int_t AliTRDtrack::CookPID(AliESDtrack *esd)
574{
575 //
576 // This function calculates the PID probabilities based on TRD signals
577 //
578 // The method produces probabilities based on the charge
579 // and the position of the maximum time bin in each layer.
580 // The dE/dx information can be used as global charge or 2 to 3
720a0a16 581 // slices. Check AliTRDCalPID and AliTRDCalPIDRefMaker for the actual
43bfc8af 582 // implementation.
583 //
584 // Author
585 // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
586
587
588 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
589 if (!calibration) {
590 AliError("No access to calibration data");
591 return -1;
592 }
593
594 // Retrieve the CDB container class with the probability distributions
720a0a16 595 const AliTRDCalPID *pd = calibration->GetPIDLQObject();
43bfc8af 596 if (!pd) {
720a0a16 597 AliError("No access to AliTRDCalPID");
43bfc8af 598 return -1;
599 }
600
601
602 Double_t p0 = 1./AliPID::kSPECIES;
603 if(AliPID::kSPECIES != 5){
604 AliError("Probabilities array defined only for 5 values. Please modify !!");
605 return -1;
606 }
607 Double_t p[] = {p0, p0, p0, p0, p0};
608 Float_t length; // track segment length in chamber
609 Int_t nPlanePID = 0;
610
611 // Skip tracks which have no TRD signal at all
612 if (fdEdx == 0.) return -1;
613
614 for (Int_t iPlane = 0; iPlane < AliTRDgeometry::kNplan; iPlane++) {
615
616 length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/TMath::Sqrt((1. - fSnp[iPlane]*fSnp[iPlane]) / (1. + fTgl[iPlane]*fTgl[iPlane]));
617
618 // check data
619 if((fdEdxPlane[iPlane][0] + fdEdxPlane[iPlane][1] + fdEdxPlane[iPlane][2]) <= 0.
620 || fTimBinPlane[iPlane] == -1.) continue;
621
622
623 // this track segment has fulfilled all requierments
624 nPlanePID++;
625
626 // Get the probabilities for the different particle species
627 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
628 p[iSpecies] *= pd->GetProbability(iSpecies, fMom[iPlane], fdEdxPlane[iPlane], length);
629 //p[iSpecies] *= pd->GetProbabilityT(iSpecies, fMom[iPlane], fTimBinPlane[iPlane]);
630 }
631 }
632 if(nPlanePID == 0) return 0;
633
634 // normalize probabilities
635 Double_t probTotal = 0.;
636 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += p[iSpecies];
637
638 if(probTotal <= 0.) {
639 AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference histograms.");
640 return -1;
641 }
642 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal;
643
644
645 // book PID to the ESD track
646 esd->SetTRDpid(p);
647 esd->SetTRDpidQuality(nPlanePID);
648
649 return 0;
650}
651
3551db50 652
a819a5f7 653
46d29e70 654//_____________________________________________________________________________
fe33d239 655Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t x0, Double_t rho)
46d29e70 656{
fe33d239 657 //
46d29e70 658 // Propagates a track of particle with mass=pm to a reference plane
659 // defined by x=xk through media of density=rho and radiationLength=x0
fe33d239 660 //
661
662 if (xk == GetX()) {
663 return kTRUE;
664 }
15ed8ba1 665
fe33d239 666 Double_t oldX = GetX();
667 Double_t oldY = GetY();
668 Double_t oldZ = GetZ();
15ed8ba1 669
fe33d239 670 Double_t bz = GetBz();
46d29e70 671
fe33d239 672 if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
673 return kFALSE;
674 }
9c9d2487 675
fe33d239 676 Double_t x = GetX();
677 Double_t y = GetY();
678 Double_t z = GetZ();
679 Double_t d = TMath::Sqrt((x-oldX)*(x-oldX)
680 + (y-oldY)*(y-oldY)
681 + (z-oldZ)*(z-oldZ));
682 if (oldX < xk) {
683 if (IsStartedTimeIntegral()) {
684 Double_t l2 = d;
685 Double_t crv = GetC();
686 if (TMath::Abs(l2*crv) > 0.0001) {
687 // Make correction for curvature if neccesary
688 l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY));
689 l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
690 l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ));
691 }
692 AddTimeStep(l2);
6c94f330 693 }
46d29e70 694 }
15ed8ba1 695
6c94f330 696 Double_t ll = (oldX < xk) ? -d : d;
fe33d239 697 if (!AliExternalTrackParam::CorrectForMaterial(ll*rho/x0,x0,GetMass())) {
698 return kFALSE;
699 }
15ed8ba1 700
fe33d239 701 {
702
703 // Energy losses************************
704 Double_t p2 = (1.0 + GetTgl()*GetTgl()) / (Get1Pt()*Get1Pt());
705 Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
d248a39d 706 if (beta2<1.0e-10 || (5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0) {
fe33d239 707 return kFALSE;
708 }
709
710 Double_t dE = 0.153e-3 / beta2
711 * (log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
712 * d * rho;
713 Float_t budget = d * rho;
714 fBudget[0] += budget;
715
716 /*
717 // Suspicious part - think about it ?
718 Double_t kinE = TMath::Sqrt(p2);
719 if (dE > 0.8*kinE) dE = 0.8 * kinE; //
720 if (dE < 0) dE = 0.0; // Not valid region for Bethe bloch
721 */
722
723 fDE += dE;
724
725 /*
726 // Suspicious ! I.B.
727 Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE)); // Energy loss fluctuation
728 Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
729 fCcc += sigmac2;
730 fCee += fX*fX * sigmac2;
731 */
6c94f330 732
0d5b5c27 733 }
734
6c94f330 735 return kTRUE;
6d45eaef 736
46d29e70 737}
738
46d29e70 739//_____________________________________________________________________________
fe33d239 740Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index
741 , Double_t h01)
46d29e70 742{
fe33d239 743 //
6c94f330 744 // Assignes found cluster to the track and updates track information
fe33d239 745 //
46d29e70 746
b8dc2353 747 Bool_t fNoTilt = kTRUE;
fe33d239 748 if (TMath::Abs(h01) > 0.003) {
749 fNoTilt = kFALSE;
750 }
751
752 // Add angular effect to the error contribution - MI
6c94f330 753 Float_t tangent2 = GetSnp()*GetSnp();
fe33d239 754 if (tangent2 < 0.90000) {
755 tangent2 = tangent2 / (1.0 - tangent2);
b8dc2353 756 }
fe33d239 757 //Float_t errang = tangent2 * 0.04;
b8dc2353 758
fe33d239 759 Double_t p[2] = {c->GetY(), c->GetZ() };
760 //Double_t cov[3] = {c->GetSigmaY2()+errang, 0.0, c->GetSigmaZ2()*100.0 };
761 Double_t sy2 = c->GetSigmaY2() * 4.0;
762 Double_t sz2 = c->GetSigmaZ2() * 4.0;
763 Double_t cov[3] = {sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
46e2d86c 764
fe33d239 765 if (!AliExternalTrackParam::Update(p,cov)) {
766 return kFALSE;
767 }
46d29e70 768
fe33d239 769 Int_t n = GetNumberOfClusters();
770 fIndex[n] = index;
b8dc2353 771 SetNumberOfClusters(n+1);
43bfc8af 772
773
b8dc2353 774 SetChi2(GetChi2()+chisq);
5443e65e 775
6c94f330 776 return kTRUE;
fe33d239 777
778}
53c17fbf 779
46e2d86c 780//_____________________________________________________________________________
43bfc8af 781Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index, Double_t h01, Int_t /*plane*/)
fe33d239 782{
783 //
46e2d86c 784 // Assignes found cluster to the track and updates track information
fe33d239 785 //
786
46e2d86c 787 Bool_t fNoTilt = kTRUE;
fe33d239 788 if (TMath::Abs(h01) > 0.003) {
789 fNoTilt = kFALSE;
790 }
791
792 // Add angular effect to the error contribution and make correction - MI
6c94f330 793 Double_t tangent2 = GetSnp()*GetSnp();
794 if (tangent2 < 0.90000){
fe33d239 795 tangent2 = tangent2 / (1.0-tangent2);
15ed8ba1 796 }
6c94f330 797 Double_t tangent = TMath::Sqrt(tangent2);
fe33d239 798 if (GetSnp() < 0) {
799 tangent *= -1;
800 }
801
802 //
803 // Is the following still needed ????
804 //
805
6c94f330 806 // Double_t correction = 0*plane;
807 /*
808 Double_t errang = tangent2*0.04; //
809 Double_t errsys =0.025*0.025*20; //systematic error part
15ed8ba1 810
6c94f330 811 Float_t extend =1;
812 if (c->GetNPads()==4) extend=2;
813 */
814 //if (c->GetNPads()==5) extend=3;
815 //if (c->GetNPads()==6) extend=3;
816 //if (c->GetQ()<15) return 1;
15ed8ba1 817
c5a8e3df 818 /*
6c94f330 819 if (corrector!=0){
3c625a9b 820 //if (0){
821 correction = corrector->GetCorrection(plane,c->GetLocalTimeBin(),tangent);
822 if (TMath::Abs(correction)>0){
823 //if we have info
824 errang = corrector->GetSigma(plane,c->GetLocalTimeBin(),tangent);
825 errang *= errang;
826 errang += tangent2*0.04;
827 }
828 }
c5a8e3df 829 */
6c94f330 830 //
831 //Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
832 /*
833 {
834 Double_t dy=c->GetY() - GetY(), dz=c->GetZ() - GetZ();
835 printf("%e %e %e %e\n",dy,dz,padlength/2,h01);
c84a5e9e 836 }
6c94f330 837 */
15ed8ba1 838
fe33d239 839 Double_t p[2] = { c->GetY(), c->GetZ() };
840 //Double_t cov[3]={ (c->GetSigmaY2()+errang+errsys)*extend, 0.0, c->GetSigmaZ2()*10000.0 };
841 Double_t sy2 = c->GetSigmaY2() * 4.0;
842 Double_t sz2 = c->GetSigmaZ2() * 4.0;
843 Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
46e2d86c 844
fe33d239 845 if (!AliExternalTrackParam::Update(p,cov)) {
846 return kFALSE;
847 }
848
43bfc8af 849 // Register cluster to track
fe33d239 850 Int_t n = GetNumberOfClusters();
851 fIndex[n] = index;
43bfc8af 852 fClusters[n] = c; // A.Bercuci 25.07.07
46e2d86c 853 SetNumberOfClusters(n+1);
fe33d239 854 SetChi2(GetChi2() + chisq);
855
856 return kTRUE;
46e2d86c 857
53c17fbf 858}
7ad19338 859
fe33d239 860// //_____________________________________________________________________________
861// Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet)
862// {
863// //
864// // Assignes found tracklet to the track and updates track information
865// //
866// // Can this be removed ????
867// //
868//
869// Double_t r00=(tracklet.GetTrackletSigma2()), r01=0., r11= 10000.;
870// r00+=fCyy; r01+=fCzy; r11+=fCzz;
871// //
872// Double_t det=r00*r11 - r01*r01;
873// Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
874// //
875
876// Double_t dy=tracklet.GetY() - fY, dz=tracklet.GetZ() - fZ;
7ad19338 877
15ed8ba1 878
fe33d239 879// Double_t s00 = tracklet.GetTrackletSigma2(); // error pad
880// Double_t s11 = 100000; // error pad-row
881// Float_t h01 = tracklet.GetTilt();
882// //
883// // r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00;
884// r00 = fCyy + fCzz*h01*h01+s00;
885// // r01 = fCzy + fCzz*h01;
886// r01 = fCzy ;
887// r11 = fCzz + s11;
888// det = r00*r11 - r01*r01;
889// // inverse matrix
890// tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
891
892// Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11;
893// Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11;
894// Double_t k20=fCey*r00+fCez*r01, k21=fCey*r01+fCez*r11;
895// Double_t k30=fCty*r00+fCtz*r01, k31=fCty*r01+fCtz*r11;
896// Double_t k40=fCcy*r00+fCcz*r01, k41=fCcy*r01+fCcz*r11;
15ed8ba1 897
fe33d239 898// // K matrix
899// // k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01);
900// // k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01);
901// // k20=fCey*r00+fCez*(r01+h01*r00),k21=fCey*r01+fCez*(r11+h01*r01);
902// // k30=fCty*r00+fCtz*(r01+h01*r00),k31=fCty*r01+fCtz*(r11+h01*r01);
903// // k40=fCcy*r00+fCcz*(r01+h01*r00),k41=fCcy*r01+fCcz*(r11+h01*r01);
904// //
905// //Update measurement
906// Double_t cur=fC + k40*dy + k41*dz, eta=fE + k20*dy + k21*dz;
907// // cur=fC + k40*dy + k41*dz; eta=fE + k20*dy + k21*dz;
908// if (TMath::Abs(cur*fX-eta) >= 0.90000) {
909// //Int_t n=GetNumberOfClusters();
910// // if (n>4) cerr<<n<<" AliTRDtrack warning: Filtering failed !\n";
911// return 0;
912// }
913// // k01+=h01*k00;
914// // k11+=h01*k10;
915// // k21+=h01*k20;
916// // k31+=h01*k30;
917// // k41+=h01*k40;
918
919
920// fY += k00*dy + k01*dz;
921// fZ += k10*dy + k11*dz;
922// fE = eta;
923// fT += k30*dy + k31*dz;
924// fC = cur;
7ad19338 925
6c94f330 926
fe33d239 927// //Update covariance
928// //
929// //
930// Double_t oldyy = fCyy, oldzz = fCzz; //, oldee=fCee, oldcc =fCcc;
931// Double_t oldzy = fCzy, oldey = fCey, oldty=fCty, oldcy =fCcy;
932// Double_t oldez = fCez, oldtz = fCtz, oldcz=fCcz;
933// //Double_t oldte = fCte, oldce = fCce;
934// //Double_t oldct = fCct;
935
936// fCyy-=k00*oldyy+k01*oldzy;
937// fCzy-=k10*oldyy+k11*oldzy;
938// fCey-=k20*oldyy+k21*oldzy;
939// fCty-=k30*oldyy+k31*oldzy;
940// fCcy-=k40*oldyy+k41*oldzy;
941// //
942// fCzz-=k10*oldzy+k11*oldzz;
943// fCez-=k20*oldzy+k21*oldzz;
944// fCtz-=k30*oldzy+k31*oldzz;
945// fCcz-=k40*oldzy+k41*oldzz;
946// //
947// fCee-=k20*oldey+k21*oldez;
948// fCte-=k30*oldey+k31*oldez;
949// fCce-=k40*oldey+k41*oldez;
950// //
951// fCtt-=k30*oldty+k31*oldtz;
952// fCct-=k40*oldty+k41*oldtz;
953// //
954// fCcc-=k40*oldcy+k41*oldcz;
955// //
15ed8ba1 956
fe33d239 957// //Int_t n=GetNumberOfClusters();
958// //fIndex[n]=index;
959// //SetNumberOfClusters(n+1);
6c94f330 960
fe33d239 961// //SetChi2(GetChi2()+chisq);
962// // cerr<<"in update: fIndex["<<fN<<"] = "<<index<<endl;
7ad19338 963
fe33d239 964// return 1;
7ad19338 965
fe33d239 966// }
c84a5e9e 967
46d29e70 968//_____________________________________________________________________________
6c94f330 969Bool_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
46d29e70 970{
fe33d239 971 //
6c94f330 972 // Rotates track parameters in R*phi plane
973 // if absolute rotation alpha is in global system
974 // otherwise alpha rotation is relative to the current rotation angle
fe33d239 975 //
976
3fad3d32 977 if (absolute) {
6c94f330 978 alpha -= GetAlpha();
3fad3d32 979 }
980 else{
981 fNRotate++;
982 }
46d29e70 983
6c94f330 984 return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
fe33d239 985
53c17fbf 986}
7ad19338 987
46d29e70 988//_____________________________________________________________________________
fd621f36 989Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
46d29e70 990{
53c17fbf 991 //
992 // Returns the track chi2
993 //
994
fe33d239 995 Double_t p[2] = { c->GetY(), c->GetZ() };
996 Double_t sy2 = c->GetSigmaY2() * 4.0;
997 Double_t sz2 = c->GetSigmaZ2() * 4.0;
998 Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
6c94f330 999
1000 return AliExternalTrackParam::GetPredictedChi2(p,cov);
15ed8ba1 1001
fe33d239 1002 //
1003 // Can the following be removed ????
1004 //
6c94f330 1005 /*
1006 Bool_t fNoTilt = kTRUE;
1007 if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE;
15ed8ba1 1008
6c94f330 1009 return (c->GetY() - GetY())*(c->GetY() - GetY())/c->GetSigmaY2();
1010 */
15ed8ba1 1011
6c94f330 1012 /*
1013 Double_t chi2, dy, r00, r01, r11;
b8dc2353 1014
6c94f330 1015 if(fNoTilt) {
1016 dy=c->GetY() - fY;
1017 r00=c->GetSigmaY2();
1018 chi2 = (dy*dy)/r00;
46d29e70 1019 }
b8dc2353 1020 else {
6c94f330 1021 Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12);
1022 //
1023 r00=c->GetSigmaY2(); r01=0.; r11=c->GetSigmaZ2();
1024 r00+=fCyy; r01+=fCzy; r11+=fCzz;
b8dc2353 1025
6c94f330 1026 Double_t det=r00*r11 - r01*r01;
b8dc2353 1027 if (TMath::Abs(det) < 1.e-10) {
6c94f330 1028 Int_t n=GetNumberOfClusters();
1029 if (n>4) cerr<<n<<" AliTRDtrack warning: Singular matrix !\n";
b8dc2353 1030 return 1e10;
1031 }
6c94f330 1032 Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
1033 Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
4f1c04d3 1034 Double_t tiltdz = dz;
6c94f330 1035 if (TMath::Abs(tiltdz)>padlength/2.) {
1036 tiltdz = TMath::Sign(padlength/2,dz);
4f1c04d3 1037 }
6c94f330 1038 // dy=dy+h01*dz;
1039 dy=dy+h01*tiltdz;
a819a5f7 1040
6c94f330 1041 chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
b8dc2353 1042 }
53c17fbf 1043
b8dc2353 1044 return chi2;
6c94f330 1045 */
fe33d239 1046
1047}
7ad19338 1048
53c17fbf 1049//_____________________________________________________________________________
16d9fbba 1050void AliTRDtrack::MakeBackupTrack()
1051{
1052 //
53c17fbf 1053 // Creates a backup track
16d9fbba 1054 //
53c17fbf 1055
fe33d239 1056 if (fBackupTrack) {
1057 delete fBackupTrack;
1058 }
16d9fbba 1059 fBackupTrack = new AliTRDtrack(*this);
4f1c04d3 1060
1061}
1062
53c17fbf 1063//_____________________________________________________________________________
1064Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
1065{
4f1c04d3 1066 //
fe33d239 1067 // Find a prolongation at given x
1068 // Return 0 if it does not exist
1069 //
1070
1071 Double_t bz = GetBz();
15ed8ba1 1072
fe33d239 1073 if (!AliExternalTrackParam::GetYAt(xk,bz,y)) {
1074 return 0;
1075 }
1076 if (!AliExternalTrackParam::GetZAt(xk,bz,z)) {
1077 return 0;
1078 }
4f1c04d3 1079
6c94f330 1080 return 1;
fe33d239 1081
16d9fbba 1082}
3fad3d32 1083
53c17fbf 1084//_____________________________________________________________________________
fe33d239 1085Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
3fad3d32 1086{
1087 //
6c94f330 1088 // Propagate track to given x position
fe33d239 1089 // Works inside of the 20 degree segmentation (local cooordinate frame for TRD , TPC, TOF)
3fad3d32 1090 //
fe33d239 1091 // Material budget from geo manager
6c94f330 1092 //
fe33d239 1093
1094 Double_t xyz0[3];
1095 Double_t xyz1[3];
1096 Double_t y;
1097 Double_t z;
1098
1099 const Double_t kAlphac = TMath::Pi()/9.0;
15ed8ba1 1100 const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
fe33d239 1101
1102 // Critical alpha - cross sector indication
1103 Double_t dir = (GetX()>xr) ? -1.0 : 1.0;
1104
1105 // Direction +-
1106 for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) {
1107
6c94f330 1108 GetXYZ(xyz0);
3fad3d32 1109 GetProlongation(x,y,z);
fe33d239 1110 xyz1[0] = x * TMath::Cos(GetAlpha()) + y * TMath::Sin(GetAlpha());
1111 xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha());
3fad3d32 1112 xyz1[2] = z;
1113 Double_t param[7];
1114 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
fe33d239 1115
1116 if ((param[0] > 0) &&
1117 (param[1] > 0)) {
1118 PropagateTo(x,param[1],param[0]);
1119 }
1120
1121 if (GetY() > GetX()*kTalphac) {
53c17fbf 1122 Rotate(-kAlphac);
3fad3d32 1123 }
fe33d239 1124 if (GetY() < -GetX()*kTalphac) {
1125 Rotate( kAlphac);
3fad3d32 1126 }
fe33d239 1127
3fad3d32 1128 }
fe33d239 1129
3fad3d32 1130 PropagateTo(xr);
53c17fbf 1131
3fad3d32 1132 return 0;
3fad3d32 1133
53c17fbf 1134}
3fad3d32 1135
53c17fbf 1136//_____________________________________________________________________________
6c94f330 1137Int_t AliTRDtrack::PropagateToR(Double_t r,Double_t step)
3fad3d32 1138{
1139 //
fe33d239 1140 // Propagate track to the radial position
1141 // Rotation always connected to the last track position
3fad3d32 1142 //
fe33d239 1143
1144 Double_t xyz0[3];
1145 Double_t xyz1[3];
1146 Double_t y;
1147 Double_t z;
1148
6c94f330 1149 Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
fe33d239 1150 // Direction +-
1151 Double_t dir = (radius>r) ? -1.0 : 1.0;
1152
1153 for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
1154
6c94f330 1155 GetXYZ(xyz0);
3fad3d32 1156 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1157 Rotate(alpha,kTRUE);
6c94f330 1158 GetXYZ(xyz0);
3fad3d32 1159 GetProlongation(x,y,z);
fe33d239 1160 xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1161 xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
3fad3d32 1162 xyz1[2] = z;
1163 Double_t param[7];
1164 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
fe33d239 1165 if (param[1] <= 0) {
1166 param[1] = 100000000;
1167 }
3fad3d32 1168 PropagateTo(x,param[1],param[0]);
fe33d239 1169
3fad3d32 1170 }
fe33d239 1171
6c94f330 1172 GetXYZ(xyz0);
3fad3d32 1173 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1174 Rotate(alpha,kTRUE);
6c94f330 1175 GetXYZ(xyz0);
3fad3d32 1176 GetProlongation(r,y,z);
fe33d239 1177 xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1178 xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
3fad3d32 1179 xyz1[2] = z;
1180 Double_t param[7];
1181 AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
fe33d239 1182
1183 if (param[1] <= 0) {
1184 param[1] = 100000000;
1185 }
3fad3d32 1186 PropagateTo(r,param[1],param[0]);
53c17fbf 1187
3fad3d32 1188 return 0;
53c17fbf 1189
1190}
1191
1192//_____________________________________________________________________________
4e28c495 1193Int_t AliTRDtrack::GetSector() const
53c17fbf 1194{
1195 //
1196 // Return the current sector
1197 //
1198
fe33d239 1199 return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha())
1200 % AliTRDgeometry::kNsect;
53c17fbf 1201
3fad3d32 1202}
1203
53c17fbf 1204//_____________________________________________________________________________
4e28c495 1205void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
53c17fbf 1206{
1207 //
1208 // The sampled energy loss
1209 //
1210
1211 Double_t s = GetSnp();
1212 Double_t t = GetTgl();
fe33d239 1213 q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
53c17fbf 1214 fdQdl[i] = q;
1215
1216}
1217
fe33d239 1218//_____________________________________________________________________________
4e28c495 1219void AliTRDtrack::SetSampledEdx(Float_t q)
53c17fbf 1220{
1221 //
1222 // The sampled energy loss
1223 //
1224
1225 Double_t s = GetSnp();
1226 Double_t t = GetTgl();
fe33d239 1227 q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
53c17fbf 1228 fdQdl[fNdedx] = q;
1229 fNdedx++;
1230
1231}
1232
fe33d239 1233//_____________________________________________________________________________
1234Double_t AliTRDtrack::GetBz() const
1235{
15ed8ba1 1236 //
fe33d239 1237 // Returns Bz component of the magnetic field (kG)
15ed8ba1 1238 //
53c17fbf 1239
fe33d239 1240 if (AliTracker::UniformField()) {
1241 return AliTracker::GetBz();
1242 }
1243 Double_t r[3];
1244 GetXYZ(r);
1245
1246 return AliTracker::GetBz(r);
53c17fbf 1247
fe33d239 1248}