]> 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
AliAlignObjAngles becomes AliAlignObjParams (Raffaele)
[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)
79143cbf 50 ,fClusterOwner(kFALSE)
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)
79143cbf 100 ,fClusterOwner(kFALSE)
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)
79143cbf 181 ,fClusterOwner(kTRUE)
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)
79143cbf 242 ,fClusterOwner(kFALSE)
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)
79143cbf 295 ,fClusterOwner(kFALSE)
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;
43bfc8af 610 // Skip tracks which have no TRD signal at all
611 if (fdEdx == 0.) return -1;
612
613 for (Int_t iPlane = 0; iPlane < AliTRDgeometry::kNplan; iPlane++) {
614
615 length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/TMath::Sqrt((1. - fSnp[iPlane]*fSnp[iPlane]) / (1. + fTgl[iPlane]*fTgl[iPlane]));
616
617 // check data
618 if((fdEdxPlane[iPlane][0] + fdEdxPlane[iPlane][1] + fdEdxPlane[iPlane][2]) <= 0.
619 || fTimBinPlane[iPlane] == -1.) continue;
620
43bfc8af 621 // this track segment has fulfilled all requierments
622 nPlanePID++;
623
624 // Get the probabilities for the different particle species
625 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
626 p[iSpecies] *= pd->GetProbability(iSpecies, fMom[iPlane], fdEdxPlane[iPlane], length);
627 //p[iSpecies] *= pd->GetProbabilityT(iSpecies, fMom[iPlane], fTimBinPlane[iPlane]);
628 }
629 }
630 if(nPlanePID == 0) return 0;
631
632 // normalize probabilities
633 Double_t probTotal = 0.;
634 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += p[iSpecies];
635
636 if(probTotal <= 0.) {
637 AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference histograms.");
638 return -1;
639 }
640 for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal;
641
642
643 // book PID to the ESD track
644 esd->SetTRDpid(p);
645 esd->SetTRDpidQuality(nPlanePID);
646
647 return 0;
648}
649
3551db50 650
a819a5f7 651
46d29e70 652//_____________________________________________________________________________
826fe01c 653Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
46d29e70 654{
fe33d239 655 //
826fe01c 656 // Propagates this track to a reference plane defined by "xk" [cm]
657 // correcting for the mean crossed material.
fe33d239 658 //
826fe01c 659 // "xx0" - thickness/rad.length [units of the radiation length]
660 // "xrho" - thickness*density [g/cm^2]
661 //
fe33d239 662
663 if (xk == GetX()) {
664 return kTRUE;
665 }
15ed8ba1 666
fe33d239 667 Double_t oldX = GetX();
668 Double_t oldY = GetY();
669 Double_t oldZ = GetZ();
15ed8ba1 670
fe33d239 671 Double_t bz = GetBz();
46d29e70 672
fe33d239 673 if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
674 return kFALSE;
675 }
9c9d2487 676
fe33d239 677 Double_t x = GetX();
678 Double_t y = GetY();
679 Double_t z = GetZ();
826fe01c 680
fe33d239 681 if (oldX < xk) {
826fe01c 682 xrho = -xrho;
fe33d239 683 if (IsStartedTimeIntegral()) {
8234f064 684 Double_t l2 = TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY) + (z-oldZ)*(z-oldZ));
fe33d239 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
826fe01c 696 if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0,xrho,GetMass())) {
fe33d239 697 return kFALSE;
698 }
15ed8ba1 699
fe33d239 700 {
701
702 // Energy losses************************
6c23ffed 703 Double_t p2 = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt());
fe33d239 704 Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
d248a39d 705 if (beta2<1.0e-10 || (5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0) {
fe33d239 706 return kFALSE;
707 }
708
709 Double_t dE = 0.153e-3 / beta2
710 * (log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
826fe01c 711 * xrho;
712 fBudget[0] += xrho;
fe33d239 713
714 /*
715 // Suspicious part - think about it ?
716 Double_t kinE = TMath::Sqrt(p2);
717 if (dE > 0.8*kinE) dE = 0.8 * kinE; //
718 if (dE < 0) dE = 0.0; // Not valid region for Bethe bloch
719 */
720
721 fDE += dE;
722
723 /*
724 // Suspicious ! I.B.
725 Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE)); // Energy loss fluctuation
726 Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
727 fCcc += sigmac2;
728 fCee += fX*fX * sigmac2;
729 */
6c94f330 730
0d5b5c27 731 }
732
6c94f330 733 return kTRUE;
6d45eaef 734
46d29e70 735}
736
46d29e70 737//_____________________________________________________________________________
fe33d239 738Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index
739 , Double_t h01)
46d29e70 740{
fe33d239 741 //
6c94f330 742 // Assignes found cluster to the track and updates track information
fe33d239 743 //
46d29e70 744
b8dc2353 745 Bool_t fNoTilt = kTRUE;
fe33d239 746 if (TMath::Abs(h01) > 0.003) {
747 fNoTilt = kFALSE;
748 }
749
750 // Add angular effect to the error contribution - MI
6c94f330 751 Float_t tangent2 = GetSnp()*GetSnp();
fe33d239 752 if (tangent2 < 0.90000) {
753 tangent2 = tangent2 / (1.0 - tangent2);
b8dc2353 754 }
fe33d239 755 //Float_t errang = tangent2 * 0.04;
b8dc2353 756
fe33d239 757 Double_t p[2] = {c->GetY(), c->GetZ() };
758 //Double_t cov[3] = {c->GetSigmaY2()+errang, 0.0, c->GetSigmaZ2()*100.0 };
759 Double_t sy2 = c->GetSigmaY2() * 4.0;
760 Double_t sz2 = c->GetSigmaZ2() * 4.0;
761 Double_t cov[3] = {sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
46e2d86c 762
fe33d239 763 if (!AliExternalTrackParam::Update(p,cov)) {
764 return kFALSE;
765 }
46d29e70 766
fe33d239 767 Int_t n = GetNumberOfClusters();
768 fIndex[n] = index;
b8dc2353 769 SetNumberOfClusters(n+1);
43bfc8af 770
771
b8dc2353 772 SetChi2(GetChi2()+chisq);
5443e65e 773
6c94f330 774 return kTRUE;
fe33d239 775
776}
53c17fbf 777
46e2d86c 778//_____________________________________________________________________________
43bfc8af 779Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index, Double_t h01, Int_t /*plane*/)
fe33d239 780{
781 //
46e2d86c 782 // Assignes found cluster to the track and updates track information
fe33d239 783 //
784
46e2d86c 785 Bool_t fNoTilt = kTRUE;
fe33d239 786 if (TMath::Abs(h01) > 0.003) {
787 fNoTilt = kFALSE;
788 }
789
790 // Add angular effect to the error contribution and make correction - MI
6c94f330 791 Double_t tangent2 = GetSnp()*GetSnp();
792 if (tangent2 < 0.90000){
fe33d239 793 tangent2 = tangent2 / (1.0-tangent2);
15ed8ba1 794 }
6c94f330 795 Double_t tangent = TMath::Sqrt(tangent2);
fe33d239 796 if (GetSnp() < 0) {
797 tangent *= -1;
798 }
799
800 //
801 // Is the following still needed ????
802 //
803
6c94f330 804 // Double_t correction = 0*plane;
805 /*
806 Double_t errang = tangent2*0.04; //
807 Double_t errsys =0.025*0.025*20; //systematic error part
15ed8ba1 808
6c94f330 809 Float_t extend =1;
810 if (c->GetNPads()==4) extend=2;
811 */
812 //if (c->GetNPads()==5) extend=3;
813 //if (c->GetNPads()==6) extend=3;
814 //if (c->GetQ()<15) return 1;
15ed8ba1 815
c5a8e3df 816 /*
6c94f330 817 if (corrector!=0){
3c625a9b 818 //if (0){
819 correction = corrector->GetCorrection(plane,c->GetLocalTimeBin(),tangent);
820 if (TMath::Abs(correction)>0){
821 //if we have info
822 errang = corrector->GetSigma(plane,c->GetLocalTimeBin(),tangent);
823 errang *= errang;
824 errang += tangent2*0.04;
825 }
826 }
c5a8e3df 827 */
6c94f330 828 //
829 //Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
830 /*
831 {
832 Double_t dy=c->GetY() - GetY(), dz=c->GetZ() - GetZ();
833 printf("%e %e %e %e\n",dy,dz,padlength/2,h01);
c84a5e9e 834 }
6c94f330 835 */
15ed8ba1 836
fe33d239 837 Double_t p[2] = { c->GetY(), c->GetZ() };
838 //Double_t cov[3]={ (c->GetSigmaY2()+errang+errsys)*extend, 0.0, c->GetSigmaZ2()*10000.0 };
839 Double_t sy2 = c->GetSigmaY2() * 4.0;
840 Double_t sz2 = c->GetSigmaZ2() * 4.0;
841 Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
46e2d86c 842
fe33d239 843 if (!AliExternalTrackParam::Update(p,cov)) {
844 return kFALSE;
845 }
846
43bfc8af 847 // Register cluster to track
fe33d239 848 Int_t n = GetNumberOfClusters();
849 fIndex[n] = index;
43bfc8af 850 fClusters[n] = c; // A.Bercuci 25.07.07
46e2d86c 851 SetNumberOfClusters(n+1);
fe33d239 852 SetChi2(GetChi2() + chisq);
853
854 return kTRUE;
46e2d86c 855
6c23ffed 856
53c17fbf 857}
7ad19338 858
fe33d239 859// //_____________________________________________________________________________
860// Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet)
861// {
862// //
863// // Assignes found tracklet to the track and updates track information
864// //
865// // Can this be removed ????
866// //
867//
868// Double_t r00=(tracklet.GetTrackletSigma2()), r01=0., r11= 10000.;
869// r00+=fCyy; r01+=fCzy; r11+=fCzz;
870// //
871// Double_t det=r00*r11 - r01*r01;
872// Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
873// //
874
875// Double_t dy=tracklet.GetY() - fY, dz=tracklet.GetZ() - fZ;
7ad19338 876
15ed8ba1 877
fe33d239 878// Double_t s00 = tracklet.GetTrackletSigma2(); // error pad
879// Double_t s11 = 100000; // error pad-row
880// Float_t h01 = tracklet.GetTilt();
881// //
882// // r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00;
883// r00 = fCyy + fCzz*h01*h01+s00;
884// // r01 = fCzy + fCzz*h01;
885// r01 = fCzy ;
886// r11 = fCzz + s11;
887// det = r00*r11 - r01*r01;
888// // inverse matrix
889// tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
890
891// Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11;
892// Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11;
893// Double_t k20=fCey*r00+fCez*r01, k21=fCey*r01+fCez*r11;
894// Double_t k30=fCty*r00+fCtz*r01, k31=fCty*r01+fCtz*r11;
895// Double_t k40=fCcy*r00+fCcz*r01, k41=fCcy*r01+fCcz*r11;
15ed8ba1 896
fe33d239 897// // K matrix
898// // k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01);
899// // k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01);
900// // k20=fCey*r00+fCez*(r01+h01*r00),k21=fCey*r01+fCez*(r11+h01*r01);
901// // k30=fCty*r00+fCtz*(r01+h01*r00),k31=fCty*r01+fCtz*(r11+h01*r01);
902// // k40=fCcy*r00+fCcz*(r01+h01*r00),k41=fCcy*r01+fCcz*(r11+h01*r01);
903// //
904// //Update measurement
905// Double_t cur=fC + k40*dy + k41*dz, eta=fE + k20*dy + k21*dz;
906// // cur=fC + k40*dy + k41*dz; eta=fE + k20*dy + k21*dz;
907// if (TMath::Abs(cur*fX-eta) >= 0.90000) {
908// //Int_t n=GetNumberOfClusters();
909// // if (n>4) cerr<<n<<" AliTRDtrack warning: Filtering failed !\n";
910// return 0;
911// }
912// // k01+=h01*k00;
913// // k11+=h01*k10;
914// // k21+=h01*k20;
915// // k31+=h01*k30;
916// // k41+=h01*k40;
917
918
919// fY += k00*dy + k01*dz;
920// fZ += k10*dy + k11*dz;
921// fE = eta;
922// fT += k30*dy + k31*dz;
923// fC = cur;
7ad19338 924
6c94f330 925
fe33d239 926// //Update covariance
927// //
928// //
929// Double_t oldyy = fCyy, oldzz = fCzz; //, oldee=fCee, oldcc =fCcc;
930// Double_t oldzy = fCzy, oldey = fCey, oldty=fCty, oldcy =fCcy;
931// Double_t oldez = fCez, oldtz = fCtz, oldcz=fCcz;
932// //Double_t oldte = fCte, oldce = fCce;
933// //Double_t oldct = fCct;
934
935// fCyy-=k00*oldyy+k01*oldzy;
936// fCzy-=k10*oldyy+k11*oldzy;
937// fCey-=k20*oldyy+k21*oldzy;
938// fCty-=k30*oldyy+k31*oldzy;
939// fCcy-=k40*oldyy+k41*oldzy;
940// //
941// fCzz-=k10*oldzy+k11*oldzz;
942// fCez-=k20*oldzy+k21*oldzz;
943// fCtz-=k30*oldzy+k31*oldzz;
944// fCcz-=k40*oldzy+k41*oldzz;
945// //
946// fCee-=k20*oldey+k21*oldez;
947// fCte-=k30*oldey+k31*oldez;
948// fCce-=k40*oldey+k41*oldez;
949// //
950// fCtt-=k30*oldty+k31*oldtz;
951// fCct-=k40*oldty+k41*oldtz;
952// //
6c23ffed 953
fe33d239 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];
826fe01c 1114 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
fe33d239 1115
1116 if ((param[0] > 0) &&
1117 (param[1] > 0)) {
826fe01c 1118 PropagateTo(x,param[1],param[0]*param[4]);
fe33d239 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];
826fe01c 1164 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
fe33d239 1165 if (param[1] <= 0) {
1166 param[1] = 100000000;
1167 }
826fe01c 1168 PropagateTo(x,param[1],param[0]*param[4]);
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];
826fe01c 1181 AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
fe33d239 1182
1183 if (param[1] <= 0) {
1184 param[1] = 100000000;
1185 }
826fe01c 1186 PropagateTo(r,param[1],param[0]*param[4]);
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}