TRD/AliTRDtrack.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use    , copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 #include <Riostream.h>
  19
  20 #include <TMath.h>
  21 #include <TVector2.h>
  22
  23 #include "AliTracker.h"
  24 #include "AliESDtrack.h"
  25
  26 #include "AliTRDgeometry.h"
  27 #include "AliTRDcluster.h"
  28 #include "AliTRDtrack.h"
  29 #include "AliTRDtracklet.h"
  30
  31 // A. Bercuci - used for PID calculations
  32 #include "AliTRDcalibDB.h"
  33 #include "Cal/AliTRDCalPID.h"
  34
  35 ClassImp(AliTRDtrack)
  36
  37 ///////////////////////////////////////////////////////////////////////////////
  38 //                                                                           //
  39 //  Represents a reconstructed TRD track                                     //
  40 //  Local TRD Kalman track                                                   //
  41 //                                                                           //
  42 ///////////////////////////////////////////////////////////////////////////////
  43
  44 //_____________________________________________________________________________
  45 AliTRDtrack::AliTRDtrack()
  46   :AliKalmanTrack()
  47   ,fSeedLab(-1)
  48   ,fdEdx(0)
  49   ,fDE(0)
  50   ,fClusterOwner(kFALSE)
  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)
  62 {
  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;
  70     }
  71     fTimBinPlane[i] = -1;
  72     fMom[i]         = -1.;
  73     fSnp[i]         = 0.;
  74     fTgl[i]         = 0.;
  75   }
  76
  77   for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
  78     fIndex[i]       = 0;
  79     fIndexBackup[i] = 0;
  80     fdQdl[i]        = 0;
  81     fClusters[i]    = 0x0;
  82   }
  83
  84   for (Int_t i = 0; i < 3; i++) {
  85     fBudget[i] = 0;
  86   }
  87
  88 }
  89
  90 //_____________________________________________________________________________
  91 AliTRDtrack::AliTRDtrack(AliTRDcluster *c, Int_t index
  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)
  98   ,fClusterOwner(kFALSE)
  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)
 109   ,fBackupTrack(0x0)
 110 {
 111   //
 112   // The main AliTRDtrack constructor.
 113   //
 114
 115   Double_t cnv   = 1.0 / (GetBz() * kB2C);
 116
 117   Double_t pp[5] = { p[0]
 118                    , p[1]
 119                    , x*p[4] - p[2]
 120                    , p[3]
 121                    , p[4]*cnv      };
 122
 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];
 128
 129   Double_t cc[15] = { cov[ 0]
 130                     , cov[ 1],     cov[ 2]
 131                     , c20,         c21,         c22
 132                     , cov[ 6],     cov[ 7],     c32,     cov[ 9]
 133                     , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
 134
 135   Set(x,alpha,pp,cc);
 136   SetNumberOfClusters(1);
 137   fIndex[0]    = index;
 138   fClusters[0] = c;
 139
 140   for (Int_t i = 0; i < kNplane; i++) {
 141     for (Int_t j = 0; j < kNslice; j++) {
 142       fdEdxPlane[i][j] = 0.0;
 143     }
 144     fTimBinPlane[i] = -1;
 145     fMom[i]         = -1.;
 146     fSnp[i]         =  0.;
 147     fTgl[i]         =  0.;
 148   }
 149
 150   Double_t q = TMath::Abs(c->GetQ());
 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   }
 156
 157   fdQdl[0] = q;
 158   for (UInt_t i = 1; i < kMAXCLUSTERSPERTRACK; i++) {
 159     fdQdl[i]        = 0;
 160     fIndex[i]       = 0;
 161     fIndexBackup[i] = 0;
 162     fClusters[i]    = 0x0;
 163   }
 164
 165   for (Int_t i = 0; i < 3;i++) {
 166     fBudget[i]      = 0;
 167   }
 168
 169 }
 170
 171 //_____________________________________________________________________________
 172 AliTRDtrack::AliTRDtrack(const AliTRDtrack &t/*, const Bool_t owner*/)
 173   :AliKalmanTrack(t)
 174   ,fSeedLab(t.GetSeedLabel())
 175   ,fdEdx(t.fdEdx)
 176   ,fDE(t.fDE)
 177   ,fClusterOwner(kTRUE)
 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)
 188   ,fBackupTrack(0x0)
 189 {
 190   //
 191   // Copy constructor.
 192   //
 193
 194   for (Int_t i = 0; i < kNplane; i++) {
 195     for (Int_t j = 0; j < kNslice; j++) {
 196       fdEdxPlane[i][j] = t.fdEdxPlane[i][j];
 197     }
 198     fTimBinPlane[i] = t.fTimBinPlane[i];
 199     fTracklets[i]   = t.fTracklets[i];
 200     fMom[i]         = t.fMom[i];
 201     fSnp[i]         = t.fSnp[i];
 202     fTgl[i]         = t.fTgl[i];
 203   }
 204
 205   Int_t n = t.GetNumberOfClusters();
 206   SetNumberOfClusters(n);
 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];
 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   }
 219
 220   for (UInt_t i = n; i < kMAXCLUSTERSPERTRACK; i++) {
 221     fdQdl[i]        = 0;
 222     fIndex[i]       = 0;
 223     fIndexBackup[i] = 0;
 224     fClusters[i]    = 0x0;
 225   }
 226
 227   for (Int_t i = 0; i < 3;i++) {
 228     fBudget[i]      = t.fBudget[i];
 229   }
 230
 231 }
 232
 233 //_____________________________________________________________________________
 234 AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t /*alpha*/)
 235   :AliKalmanTrack(t)
 236   ,fSeedLab(-1)
 237   ,fdEdx(t.GetPIDsignal())
 238   ,fDE(0)
 239   ,fClusterOwner(kFALSE)
 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)
 251 {
 252   //
 253   // Constructor from AliTPCtrack or AliITStrack
 254   //
 255
 256   SetLabel(t.GetLabel());
 257   SetChi2(0.0);
 258   SetMass(t.GetMass());
 259   SetNumberOfClusters(0);
 260
 261   for (Int_t i = 0; i < kNplane; i++) {
 262     for (Int_t j = 0; j < kNslice; j++) {
 263       fdEdxPlane[i][j] = 0.0;
 264     }
 265     fTimBinPlane[i] = -1;
 266     fMom[i]         = -1.;
 267     fSnp[i]         =  0.;
 268     fTgl[i]         =  0.;
 269   }
 270
 271   for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
 272     fdQdl[i]        = 0;
 273     fIndex[i]       = 0;
 274     fIndexBackup[i] = 0;
 275     fClusters[i]    = 0x0;
 276   }
 277
 278   for (Int_t i = 0; i < 3; i++) {
 279     fBudget[i]      = 0;
 280   }
 281
 282 }
 283
 284 //_____________________________________________________________________________
 285 AliTRDtrack::AliTRDtrack(const AliESDtrack &t)
 286   :AliKalmanTrack()
 287   ,fSeedLab(-1)
 288   ,fdEdx(t.GetTRDsignal())
 289   ,fDE(0)
 290   ,fClusterOwner(kFALSE)
 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)
 302 {
 303   //
 304   // Constructor from AliESDtrack
 305   //
 306
 307   SetLabel(t.GetLabel());
 308   SetChi2(0.0);
 309   SetMass(t.GetMass());
 310   SetNumberOfClusters(t.GetTRDclusters(fIndex));
 311
 312   Int_t ncl = t.GetTRDclusters(fIndexBackup);
 313   for (UInt_t i = ncl; i < kMAXCLUSTERSPERTRACK; i++) {
 314     fIndexBackup[i] = 0;
 315     fIndex[i]       = 0;
 316   }
 317
 318   for (Int_t i = 0; i < kNplane; i++) {
 319     for (Int_t j = 0; j < kNslice; j++) {
 320       fdEdxPlane[i][j] = t.GetTRDsignals(i,j);
 321     }
 322     fTimBinPlane[i] = t.GetTRDTimBin(i);
 323     fMom[i]         = -1.;
 324     fSnp[i]         =  0.;
 325     fTgl[i]         =  0.;
 326   }
 327
 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     }
 335   }
 336   Set(par->GetX(),par->GetAlpha(),par->GetParameter(),par->GetCovariance());
 337
 338   for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
 339     fdQdl[i]     = 0;
 340     fClusters[i] = 0x0;
 341   }
 342
 343   for (Int_t i = 0; i < 3; i++) {
 344     fBudget[i] = 0;
 345   }
 346
 347   if ((t.GetStatus() & AliESDtrack::kTIME) == 0) {
 348     return;
 349   }
 350
 351   StartTimeIntegral();
 352   Double_t times[10];
 353   t.GetIntegratedTimes(times);
 354   SetIntegratedTimes(times);
 355   SetIntegratedLength(t.GetIntegratedLength());
 356
 357 }
 358
 359 //____________________________________________________________________________
 360 AliTRDtrack::~AliTRDtrack()
 361 {
 362   //
 363   // Destructor
 364   //
 365
 366   if (fBackupTrack) {
 367     delete fBackupTrack;
 368   }
 369   fBackupTrack = 0x0;
 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     }
 377   }
 378
 379 }
 380
 381 //____________________________________________________________________________
 382 Float_t AliTRDtrack::StatusForTOF()
 383 {
 384   //
 385   // Defines the status of the TOF extrapolation
 386   //
 387
 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]));
 392   return res;
 393
 394   // This part of the function is never reached ????
 395   // What defines these parameters ????
 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;
 412
 413 }
 414
 415 //_____________________________________________________________________________
 416 Int_t AliTRDtrack::Compare(const TObject *o) const
 417 {
 418   //
 419   // Compares tracks according to their Y2 or curvature
 420   //
 421
 422   AliTRDtrack *t = (AliTRDtrack *) o;
 423
 424   Double_t co = TMath::Abs(t->GetC());
 425   Double_t c  = TMath::Abs(GetC());
 426
 427   if      (c > co) {
 428     return 1;
 429   }
 430   else if (c < co) {
 431     return -1;
 432   }
 433   return 0;
 434
 435 }
 436
 437 //_____________________________________________________________________________
 438 void AliTRDtrack::CookdEdx(Double_t low, Double_t up)
 439 {
 440   //
 441   // Calculates the truncated dE/dx within the "low" and "up" cuts.
 442   //
 443
 444   // Array to sort the dEdx values according to amplitude
 445   Float_t sorted[kMAXCLUSTERSPERTRACK];
 446   fdEdx = 0.;
 447
 448   // Require at least 10 clusters for a dedx measurement
 449   if (fNdedx < 10) return;
 450
 451   // Can fdQdl be negative ????
 452   for (Int_t i = 0; i < fNdedx; i++) {
 453     sorted[i] = TMath::Abs(fdQdl[i]);
 454   }
 455   // Sort the dedx values by amplitude
 456   Int_t *index = new Int_t[fNdedx];
 457   TMath::Sort(fNdedx, sorted, index, kFALSE);
 458
 459   // Sum up the truncated charge between lower and upper bounds
 460   Int_t nl = Int_t(low * fNdedx);
 461   Int_t nu = Int_t( up * fNdedx);
 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;
 468
 469 }
 470
 471 //_____________________________________________________________________________
 472 void AliTRDtrack::CookdEdxTimBin(const Int_t/* tid*/)
 473 {
 474   //
 475   // Set fdEdxPlane and fTimBinPlane and also get the
 476   // Time bin for Max. Cluster
 477   //
 478   // Authors:
 479   // Prashant Shukla (shukla@physi.uni-heidelberg.de)
 480   // Alexandru Bercuci (A.Bercuci@gsi.de)
 481   //
 482
 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
 494   for (Int_t iPlane = 0; iPlane < AliESDtrack::kNPlane; iPlane++) {
 495     fTimBinPlane[iPlane] = -1;
 496     maxcharge[iPlane]    = 0.;
 497     for (Int_t iSlice = 0; iSlice < AliESDtrack::kNSlice; iSlice++) {
 498       fdEdxPlane[iPlane][iSlice] = 0.;
 499       nCluster[iPlane][iSlice]   = 0;
 500     }
 501   }
 502
 503   // Start looping over clusters attached to this track
 504   for (Int_t iClus = 0; iClus < GetNumberOfClusters(); iClus++) {
 505
 506     cluster = fClusters[iClus]; //(AliTRDcluster*)tracker->GetCluster(fIndex[iClus]);
 507     if(!cluster) continue;
 508
 509     // Read info from current cluster
 510     plane  = AliTRDgeometry::GetPlane(cluster->GetDetector());
 511     if (plane < 0 || plane >= AliESDtrack::kNPlane) {
 512       AliError(Form("Wrong plane %d", plane));
 513       continue;
 514     }
 515
 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     }
 522
 523     slice = tb * AliESDtrack::kNSlice / ntb;
 524
 525     fdEdxPlane[plane][slice] += fdQdl[iClus];
 526     if (fdQdl[iClus] > maxcharge[plane]) {
 527       maxcharge[plane]    = fdQdl[iClus];
 528       fTimBinPlane[plane] = tb;
 529     }
 530
 531     nCluster[plane][slice]++;
 532
 533   } // End of loop over cluster
 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++) {
 538       if (nCluster[iPlane][iSlice]) {
 539         fdEdxPlane[iPlane][iSlice] /= nCluster[iPlane][iSlice];
 540       }
 541     }
 542   }
 543
 544 }
 545
 546 //_____________________________________________________________________________
 547 void    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 //_____________________________________________________________________________
 594 void    AliTRDtrack::SetTrackSegmentDirMom(const Int_t plane)
 595 {
 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   }
 605
 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
 612 }
 613
 614 //_____________________________________________________________________________
 615 Float_t AliTRDtrack::GetTrackLengthPlane(Int_t plane) const
 616 {
 617   //
 618   // Calculate the track length of a track segment in a given plane
 619   //
 620
 621   if ((plane < 0) || (plane >= kNplane)) return 0.;
 622
 623   return (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())
 624         / TMath::Sqrt((1.0 - fSnp[plane]*fSnp[plane])
 625                     / (1.0 + fTgl[plane]*fTgl[plane]));
 626
 627 }
 628
 629 //_____________________________________________________________________________
 630 Bool_t AliTRDtrack::CookPID(Int_t &pidQuality)
 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
 638   // slices. Check AliTRDCalPID and AliTRDCalPIDRefMaker for the actual
 639   // implementation.
 640   //
 641   // Author
 642   // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
 643   //
 644
 645   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 646   if (!calibration) {
 647     AliError("No access to calibration data");
 648     return kFALSE;
 649   }
 650
 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   }
 657
 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   }
 663
 664   // Sets the a priori probabilities
 665   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
 666     fPID[ispec] = 1./AliPID::kSPECIES;
 667   }
 668
 669
 670   if(AliPID::kSPECIES != 5){
 671     AliError("Probabilities array defined only for 5 values. Please modify !!");
 672     return kFALSE;
 673   }
 674
 675
 676   pidQuality = 0;
 677   Float_t length;  // track segment length in chamber
 678
 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;
 707
 708   // normalize probabilities
 709   Double_t probTotal = 0.;
 710   for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += fPID[iSpecies];
 711
 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;
 717
 718   return kTRUE;
 719 }
 720
 721
 722 //_____________________________________________________________________________
 723 Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
 724 {
 725   //
 726   // Propagates this track to a reference plane defined by "xk" [cm]
 727   // correcting for the mean crossed material.
 728   //
 729   // "xx0"  - thickness/rad.length [units of the radiation length]
 730   // "xrho" - thickness*density    [g/cm^2]
 731   //
 732
 733   if (xk == GetX()) {
 734     return kTRUE;
 735   }
 736
 737   Double_t oldX = GetX();
 738   Double_t oldY = GetY();
 739   Double_t oldZ = GetZ();
 740
 741   Double_t bz   = GetBz();
 742
 743   if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
 744     return kFALSE;
 745   }
 746
 747   Double_t x = GetX();
 748   Double_t y = GetY();
 749   Double_t z = GetZ();
 750
 751   if (oldX < xk) {
 752     xrho = -xrho;
 753     if (IsStartedTimeIntegral()) {
 754       Double_t l2  = TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY) + (z-oldZ)*(z-oldZ));
 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);
 763     }
 764   }
 765
 766   if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0,xrho,GetMass())) {
 767     return kFALSE;
 768   }
 769
 770   {
 771
 772     // Energy losses************************
 773     Double_t p2    = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt());
 774     Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
 775     if (beta2<1.0e-10 || (5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0) {
 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)
 781                    * xrho;
 782     fBudget[0] += xrho;
 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     */
 800
 801   }
 802
 803   return kTRUE;
 804
 805 }
 806
 807 //_____________________________________________________________________________
 808 Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index
 809                          , Double_t h01)
 810 {
 811   //
 812   // Assignes found cluster to the track and updates track information
 813   //
 814
 815   Bool_t fNoTilt = kTRUE;
 816   if (TMath::Abs(h01) > 0.003) {
 817     fNoTilt = kFALSE;
 818   }
 819
 820   // Add angular effect to the error contribution -  MI
 821   Float_t tangent2 = GetSnp()*GetSnp();
 822   if (tangent2 < 0.90000) {
 823     tangent2 = tangent2 / (1.0 - tangent2);
 824   }
 825   //Float_t errang = tangent2 * 0.04;
 826
 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;
 831   Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
 832
 833   if (!AliExternalTrackParam::Update(p,cov)) {
 834     return kFALSE;
 835   }
 836
 837   Int_t n   = GetNumberOfClusters();
 838   fIndex[n] = index;
 839   SetNumberOfClusters(n+1);
 840
 841   SetChi2(GetChi2()+chisq);
 842
 843   return kTRUE;
 844
 845 }
 846
 847 //_____________________________________________________________________________
 848 Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index
 849                           , Double_t h01, Int_t /*plane*/, Int_t tid)
 850 {
 851   //
 852   // Assignes found cluster to the track and updates track information
 853   //
 854
 855   Bool_t fNoTilt = kTRUE;
 856   if (TMath::Abs(h01) > 0.003) {
 857     fNoTilt = kFALSE;
 858   }
 859
 860   // Add angular effect to the error contribution and make correction  -  MI
 861   Double_t tangent2 = GetSnp()*GetSnp();
 862   if (tangent2 < 0.90000){
 863     tangent2 = tangent2 / (1.0-tangent2);
 864   }
 865   Double_t tangent = TMath::Sqrt(tangent2);
 866   if (GetSnp() < 0) {
 867     tangent *= -1;
 868   }
 869
 870   //
 871   // Is the following still needed ????
 872   //
 873
 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
 878
 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;
 885
 886   /*
 887   if (corrector!=0){
 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   }
 897   */
 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);
 904     }
 905   */
 906
 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 };
 912
 913   if (!AliExternalTrackParam::Update(p,cov)) {
 914     return kFALSE;
 915   }
 916
 917   // Register cluster to track
 918   Int_t n      = GetNumberOfClusters();
 919   fIndex[n]    = index;
 920   fClusters[n] = c;
 921   SetNumberOfClusters(n+1);
 922   SetChi2(GetChi2() + chisq);
 923
 924   return kTRUE;
 925
 926 }
 927
 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;
 945
 946
 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;
 965
 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;
 993
 994
 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 //   //
1022
1023 //   fCcc-=k40*oldcy+k41*oldcz;
1024 //   //
1025
1026 //   //Int_t n=GetNumberOfClusters();
1027 //   //fIndex[n]=index;
1028 //   //SetNumberOfClusters(n+1);
1029
1030 //   //SetChi2(GetChi2()+chisq);
1031 //   //  cerr<<"in update: fIndex["<<fN<<"] = "<<index<<endl;
1032
1033 //   return 1;
1034
1035 // }
1036
1037 //_____________________________________________________________________________
1038 Bool_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
1039 {
1040   //
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
1044   //
1045
1046   if (absolute) {
1047     alpha -= GetAlpha();
1048   }
1049   else{
1050     fNRotate++;
1051   }
1052
1053   return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
1054
1055 }
1056
1057 //_____________________________________________________________________________
1058 Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
1059 {
1060   //
1061   // Returns the track chi2
1062   //
1063
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 };
1068
1069   return AliExternalTrackParam::GetPredictedChi2(p,cov);
1070
1071   //
1072   // Can the following be removed ????
1073   //
1074   /*
1075   Bool_t fNoTilt = kTRUE;
1076   if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE;
1077
1078   return (c->GetY() - GetY())*(c->GetY() - GetY())/c->GetSigmaY2();
1079   */
1080
1081   /*
1082   Double_t chi2, dy, r00, r01, r11;
1083
1084   if(fNoTilt) {
1085     dy=c->GetY() - fY;
1086     r00=c->GetSigmaY2();
1087     chi2 = (dy*dy)/r00;
1088   }
1089   else {
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;
1094
1095     Double_t det=r00*r11 - r01*r01;
1096     if (TMath::Abs(det) < 1.e-10) {
1097       Int_t n=GetNumberOfClusters();
1098       if (n>4) cerr<<n<<" AliTRDtrack warning: Singular matrix !\n";
1099       return 1e10;
1100     }
1101     Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
1102     Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
1103     Double_t tiltdz = dz;
1104     if (TMath::Abs(tiltdz)>padlength/2.) {
1105       tiltdz = TMath::Sign(padlength/2,dz);
1106     }
1107     //    dy=dy+h01*dz;
1108     dy=dy+h01*tiltdz;
1109
1110     chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
1111   }
1112
1113   return chi2;
1114   */
1115
1116 }
1117
1118 //_____________________________________________________________________________
1119 void AliTRDtrack::MakeBackupTrack()
1120 {
1121   //
1122   // Creates a backup track
1123   //
1124
1125   if (fBackupTrack) {
1126     delete fBackupTrack;
1127   }
1128   fBackupTrack = new AliTRDtrack(*this);
1129
1130 }
1131
1132 //_____________________________________________________________________________
1133 Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
1134 {
1135   //
1136   // Find a prolongation at given x
1137   // Return 0 if it does not exist
1138   //
1139
1140   Double_t bz = GetBz();
1141
1142   if (!AliExternalTrackParam::GetYAt(xk,bz,y)) {
1143     return 0;
1144   }
1145   if (!AliExternalTrackParam::GetZAt(xk,bz,z)) {
1146     return 0;
1147   }
1148
1149   return 1;
1150
1151 }
1152
1153 //_____________________________________________________________________________
1154 Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
1155 {
1156   //
1157   // Propagate track to given x  position
1158   // Works inside of the 20 degree segmentation
1159   // (local cooordinate frame for TRD , TPC, TOF)
1160   //
1161   // Material budget from geo manager
1162   //
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;
1170   const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
1171
1172   // Critical alpha  - cross sector indication
1173   Double_t dir = (GetX() > xr) ? -1.0 : 1.0;
1174
1175   // Direction +-
1176   for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) {
1177
1178     GetXYZ(xyz0);
1179     GetProlongation(x,y,z);
1180     xyz1[0] = x * TMath::Cos(GetAlpha()) + y * TMath::Sin(GetAlpha());
1181     xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha());
1182     xyz1[2] = z;
1183     Double_t param[7];
1184     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1185
1186     if ((param[0] > 0) &&
1187         (param[1] > 0)) {
1188       PropagateTo(x,param[1],param[0]*param[4]);
1189     }
1190
1191     if (GetY() >  GetX()*kTalphac) {
1192       Rotate(-kAlphac);
1193     }
1194     if (GetY() < -GetX()*kTalphac) {
1195       Rotate( kAlphac);
1196     }
1197
1198   }
1199
1200   PropagateTo(xr);
1201
1202   return 0;
1203
1204 }
1205
1206 //_____________________________________________________________________________
1207 Int_t   AliTRDtrack::PropagateToR(Double_t r,Double_t step)
1208 {
1209   //
1210   // Propagate track to the radial position
1211   // Rotation always connected to the last track position
1212   //
1213
1214   Double_t xyz0[3];
1215   Double_t xyz1[3];
1216   Double_t y;
1217   Double_t z;
1218
1219   Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
1220   // Direction +-
1221   Double_t dir    = (radius > r) ? -1.0 : 1.0;
1222
1223   for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
1224
1225     GetXYZ(xyz0);
1226     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1227     Rotate(alpha,kTRUE);
1228     GetXYZ(xyz0);
1229     GetProlongation(x,y,z);
1230     xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1231     xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1232     xyz1[2] = z;
1233     Double_t param[7];
1234     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1235     if (param[1] <= 0) {
1236       param[1] = 100000000;
1237     }
1238     PropagateTo(x,param[1],param[0]*param[4]);
1239
1240   }
1241
1242   GetXYZ(xyz0);
1243   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1244   Rotate(alpha,kTRUE);
1245   GetXYZ(xyz0);
1246   GetProlongation(r,y,z);
1247   xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1248   xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1249   xyz1[2] = z;
1250   Double_t param[7];
1251   AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1252
1253   if (param[1] <= 0) {
1254     param[1] = 100000000;
1255   }
1256   PropagateTo(r,param[1],param[0]*param[4]);
1257
1258   return 0;
1259
1260 }
1261
1262 //_____________________________________________________________________________
1263 Int_t AliTRDtrack::GetSector() const
1264 {
1265   //
1266   // Return the current sector
1267   //
1268
1269   return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha())
1270              % AliTRDgeometry::kNsect;
1271
1272 }
1273
1274 //_____________________________________________________________________________
1275 void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
1276 {
1277   //
1278   // The sampled energy loss
1279   //
1280
1281   Double_t s = GetSnp();
1282   Double_t t = GetTgl();
1283   q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1284   fdQdl[i] = q;
1285
1286 }
1287
1288 //_____________________________________________________________________________
1289 void AliTRDtrack::SetSampledEdx(Float_t q)
1290 {
1291   //
1292   // The sampled energy loss
1293   //
1294
1295   Double_t s = GetSnp();
1296   Double_t t = GetTgl();
1297   q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1298   fdQdl[fNdedx] = q;
1299   fNdedx++;
1300
1301 }
1302
1303 //_____________________________________________________________________________
1304 Double_t AliTRDtrack::GetBz() const
1305 {
1306   //
1307   // Returns Bz component of the magnetic field (kG)
1308   //
1309
1310   if (AliTracker::UniformField()) {
1311     return AliTracker::GetBz();
1312   }
1313   Double_t r[3];
1314   GetXYZ(r);
1315
1316   return AliTracker::GetBz(r);
1317
1318 }