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