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()
 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::SetTrackSegmentDirMom(const Int_t plane)
 548 {
 549   //
 550   // Set the momenta for a track segment in a given plane
 551   //
 552
 553   if ((plane <       0) ||
 554       (plane>= kNplane)) {
 555     AliError(Form("Trying to access out of range plane (%d)", plane));
 556     return;
 557   }
 558
 559   fSnp[plane] = GetSnp();
 560   fTgl[plane] = GetTgl();
 561   Double_t p[3];
 562   GetPxPyPz(p);
 563   fMom[plane] = TMath::Sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
 564
 565 }
 566
 567 //_____________________________________________________________________________
 568 Float_t AliTRDtrack::GetTrackLengthPlane(Int_t plane) const
 569 {
 570   //
 571   // Calculate the track length of a track segment in a given plane
 572   //
 573
 574   if ((plane < 0) || (plane >= kNplane)) return 0.;
 575
 576   return (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())
 577         / TMath::Sqrt((1.0 - fSnp[plane]*fSnp[plane])
 578                     / (1.0 + fTgl[plane]*fTgl[plane]));
 579
 580 }
 581
 582 //_____________________________________________________________________________
 583 Int_t AliTRDtrack::CookPID(AliESDtrack *esd)
 584 {
 585   //
 586   // This function calculates the PID probabilities based on TRD signals
 587   //
 588   // The method produces probabilities based on the charge
 589   // and the position of the maximum time bin in each layer.
 590   // The dE/dx information can be used as global charge or 2 to 3
 591   // slices. Check AliTRDCalPID and AliTRDCalPIDRefMaker for the actual
 592   // implementation.
 593   //
 594   // Author
 595   // Alex Bercuci (A.Bercuci@gsi.de) 2nd May 2007
 596   //
 597
 598   AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 599   if (!calibration) {
 600     AliError("No access to calibration data");
 601     return -1;
 602   }
 603
 604         // Retrieve the CDB container class with the probability distributions
 605         const AliTRDCalPID *pd = calibration->GetPIDLQObject();
 606         if (!pd) {
 607                 AliError("No access to AliTRDCalPID");
 608                 return -1;
 609         }
 610
 611
 612         Double_t p0 = 1./AliPID::kSPECIES;
 613         if(AliPID::kSPECIES != 5){
 614                 AliError("Probabilities array defined only for 5 values. Please modify !!");
 615                 return -1;
 616         }
 617         Double_t p[] = {p0, p0, p0, p0, p0};
 618         Float_t length;  // track segment length in chamber
 619         Int_t nPlanePID = 0;
 620         // Skip tracks which have no TRD signal at all
 621         if (fdEdx == 0.) return -1;
 622
 623         for (Int_t iPlane = 0; iPlane < AliTRDgeometry::kNplan; iPlane++) {
 624
 625                 length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/TMath::Sqrt((1. - fSnp[iPlane]*fSnp[iPlane]) / (1. + fTgl[iPlane]*fTgl[iPlane]));
 626
 627                 // check data
 628                 if((fdEdxPlane[iPlane][0] + fdEdxPlane[iPlane][1] + fdEdxPlane[iPlane][2]) <=  0.
 629                 || fTimBinPlane[iPlane] == -1.) continue;
 630
 631                 // this track segment has fulfilled all requierments
 632                 nPlanePID++;
 633
 634                 // Get the probabilities for the different particle species
 635                 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
 636                         p[iSpecies] *= pd->GetProbability(iSpecies, fMom[iPlane], fdEdxPlane[iPlane], length);
 637                         //p[iSpecies] *= pd->GetProbabilityT(iSpecies, fMom[iPlane], fTimBinPlane[iPlane]);
 638                 }
 639         }
 640         if(nPlanePID == 0) return 0;
 641
 642         // normalize probabilities
 643         Double_t probTotal = 0.;
 644         for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += p[iSpecies];
 645
 646         if(probTotal <= 0.) {
 647                 AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference histograms.");
 648                 return -1;
 649         }
 650         for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal;
 651
 652
 653         // book PID to the ESD track
 654         esd->SetTRDpid(p);
 655         esd->SetTRDpidQuality(nPlanePID);
 656
 657         return 0;
 658 }
 659
 660 //_____________________________________________________________________________
 661 Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
 662 {
 663   //
 664   // Propagates this track to a reference plane defined by "xk" [cm]
 665   // correcting for the mean crossed material.
 666   //
 667   // "xx0"  - thickness/rad.length [units of the radiation length]
 668   // "xrho" - thickness*density    [g/cm^2]
 669   //
 670
 671   if (xk == GetX()) {
 672     return kTRUE;
 673   }
 674
 675   Double_t oldX = GetX();
 676   Double_t oldY = GetY();
 677   Double_t oldZ = GetZ();
 678
 679   Double_t bz   = GetBz();
 680
 681   if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
 682     return kFALSE;
 683   }
 684
 685   Double_t x = GetX();
 686   Double_t y = GetY();
 687   Double_t z = GetZ();
 688
 689   if (oldX < xk) {
 690     xrho = -xrho;
 691     if (IsStartedTimeIntegral()) {
 692       Double_t l2  = TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY) + (z-oldZ)*(z-oldZ));
 693       Double_t crv = GetC();
 694       if (TMath::Abs(l2*crv) > 0.0001) {
 695         // Make correction for curvature if neccesary
 696         l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY));
 697         l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
 698         l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ));
 699       }
 700       AddTimeStep(l2);
 701     }
 702   }
 703
 704   if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0,xrho,GetMass())) {
 705     return kFALSE;
 706   }
 707
 708   {
 709
 710     // Energy losses************************
 711     Double_t p2    = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt());
 712     Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
 713     if (beta2<1.0e-10 || (5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0) {
 714       return kFALSE;
 715     }
 716
 717     Double_t dE    = 0.153e-3 / beta2
 718                    * (log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
 719                    * xrho;
 720     fBudget[0] += xrho;
 721
 722     /*
 723     // Suspicious part - think about it ?
 724     Double_t kinE =  TMath::Sqrt(p2);
 725     if (dE > 0.8*kinE) dE = 0.8 * kinE;  //
 726     if (dE < 0)        dE = 0.0;         // Not valid region for Bethe bloch
 727     */
 728
 729     fDE += dE;
 730
 731     /*
 732     // Suspicious ! I.B.
 733     Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE));   // Energy loss fluctuation
 734     Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
 735     fCcc += sigmac2;
 736     fCee += fX*fX * sigmac2;
 737     */
 738
 739   }
 740
 741   return kTRUE;
 742
 743 }
 744
 745 //_____________________________________________________________________________
 746 Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index
 747                          , Double_t h01)
 748 {
 749   //
 750   // Assignes found cluster to the track and updates track information
 751   //
 752
 753   Bool_t fNoTilt = kTRUE;
 754   if (TMath::Abs(h01) > 0.003) {
 755     fNoTilt = kFALSE;
 756   }
 757
 758   // Add angular effect to the error contribution -  MI
 759   Float_t tangent2 = GetSnp()*GetSnp();
 760   if (tangent2 < 0.90000) {
 761     tangent2 = tangent2 / (1.0 - tangent2);
 762   }
 763   //Float_t errang = tangent2 * 0.04;
 764
 765   Double_t p[2]   = {c->GetY(), c->GetZ() };
 766   //Double_t cov[3] = {c->GetSigmaY2()+errang, 0.0, c->GetSigmaZ2()*100.0 };
 767   Double_t sy2    = c->GetSigmaY2() * 4.0;
 768   Double_t sz2    = c->GetSigmaZ2() * 4.0;
 769   Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
 770
 771   if (!AliExternalTrackParam::Update(p,cov)) {
 772     return kFALSE;
 773   }
 774
 775   Int_t n   = GetNumberOfClusters();
 776   fIndex[n] = index;
 777   SetNumberOfClusters(n+1);
 778
 779   SetChi2(GetChi2()+chisq);
 780
 781   return kTRUE;
 782
 783 }
 784
 785 //_____________________________________________________________________________
 786 Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index
 787                           , Double_t h01, Int_t /*plane*/)
 788 {
 789   //
 790   // Assignes found cluster to the track and updates track information
 791   //
 792
 793   Bool_t fNoTilt = kTRUE;
 794   if (TMath::Abs(h01) > 0.003) {
 795     fNoTilt = kFALSE;
 796   }
 797
 798   // Add angular effect to the error contribution and make correction  -  MI
 799   Double_t tangent2 = GetSnp()*GetSnp();
 800   if (tangent2 < 0.90000){
 801     tangent2 = tangent2 / (1.0-tangent2);
 802   }
 803   Double_t tangent = TMath::Sqrt(tangent2);
 804   if (GetSnp() < 0) {
 805     tangent *= -1;
 806   }
 807
 808   //
 809   // Is the following still needed ????
 810   //
 811
 812   //  Double_t correction = 0*plane;
 813   /*
 814   Double_t errang = tangent2*0.04;  //
 815   Double_t errsys =0.025*0.025*20;  //systematic error part
 816
 817   Float_t extend =1;
 818   if (c->GetNPads()==4) extend=2;
 819   */
 820   //if (c->GetNPads()==5)  extend=3;
 821   //if (c->GetNPads()==6)  extend=3;
 822   //if (c->GetQ()<15) return 1;
 823
 824   /*
 825   if (corrector!=0){
 826   //if (0){
 827     correction = corrector->GetCorrection(plane,c->GetLocalTimeBin(),tangent);
 828     if (TMath::Abs(correction)>0){
 829       //if we have info
 830       errang     = corrector->GetSigma(plane,c->GetLocalTimeBin(),tangent);
 831       errang    *= errang;
 832       errang    += tangent2*0.04;
 833     }
 834   }
 835   */
 836   //
 837   //Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
 838   /*
 839     {
 840       Double_t dy=c->GetY() - GetY(), dz=c->GetZ() - GetZ();
 841       printf("%e %e %e %e\n",dy,dz,padlength/2,h01);
 842     }
 843   */
 844
 845   Double_t p[2]   = { c->GetY(), c->GetZ() };
 846   //Double_t cov[3]={ (c->GetSigmaY2()+errang+errsys)*extend, 0.0, c->GetSigmaZ2()*10000.0 };
 847   Double_t sy2    = c->GetSigmaY2() * 4.0;
 848   Double_t sz2    = c->GetSigmaZ2() * 4.0;
 849   Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
 850
 851   if (!AliExternalTrackParam::Update(p,cov)) {
 852     return kFALSE;
 853   }
 854
 855   // Register cluster to track
 856   Int_t n      = GetNumberOfClusters();
 857   fIndex[n]    = index;
 858   fClusters[n] = c;
 859   SetNumberOfClusters(n+1);
 860   SetChi2(GetChi2() + chisq);
 861
 862   return kTRUE;
 863
 864 }
 865
 866 // //_____________________________________________________________________________
 867 // Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet)
 868 // {
 869 //   //
 870 //   // Assignes found tracklet to the track and updates track information
 871 //   //
 872 //   // Can this be removed ????
 873 //   //
 874 //
 875 //   Double_t r00=(tracklet.GetTrackletSigma2()), r01=0., r11= 10000.;
 876 //   r00+=fCyy; r01+=fCzy; r11+=fCzz;
 877 //   //
 878 //   Double_t det=r00*r11 - r01*r01;
 879 //   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 880 //   //
 881
 882 //   Double_t dy=tracklet.GetY() - fY, dz=tracklet.GetZ() - fZ;
 883
 884
 885 //   Double_t s00 = tracklet.GetTrackletSigma2();  // error pad
 886 //   Double_t s11 = 100000;   // error pad-row
 887 //   Float_t  h01 = tracklet.GetTilt();
 888 //   //
 889 //   //  r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00;
 890 //   r00 = fCyy + fCzz*h01*h01+s00;
 891 //   //  r01 = fCzy + fCzz*h01;
 892 //   r01 = fCzy ;
 893 //   r11 = fCzz + s11;
 894 //   det = r00*r11 - r01*r01;
 895 //   // inverse matrix
 896 //   tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 897
 898 //   Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11;
 899 //   Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11;
 900 //   Double_t k20=fCey*r00+fCez*r01, k21=fCey*r01+fCez*r11;
 901 //   Double_t k30=fCty*r00+fCtz*r01, k31=fCty*r01+fCtz*r11;
 902 //   Double_t k40=fCcy*r00+fCcz*r01, k41=fCcy*r01+fCcz*r11;
 903
 904 //   // K matrix
 905 // //   k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01);
 906 // //   k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01);
 907 // //   k20=fCey*r00+fCez*(r01+h01*r00),k21=fCey*r01+fCez*(r11+h01*r01);
 908 // //   k30=fCty*r00+fCtz*(r01+h01*r00),k31=fCty*r01+fCtz*(r11+h01*r01);
 909 // //   k40=fCcy*r00+fCcz*(r01+h01*r00),k41=fCcy*r01+fCcz*(r11+h01*r01);
 910 //   //
 911 //   //Update measurement
 912 //   Double_t cur=fC + k40*dy + k41*dz, eta=fE + k20*dy + k21*dz;
 913 //   //  cur=fC + k40*dy + k41*dz; eta=fE + k20*dy + k21*dz;
 914 //   if (TMath::Abs(cur*fX-eta) >= 0.90000) {
 915 //     //Int_t n=GetNumberOfClusters();
 916 //     //      if (n>4) cerr<<n<<" AliTRDtrack warning: Filtering failed !\n";
 917 //     return 0;
 918 //   }
 919 // //   k01+=h01*k00;
 920 // //   k11+=h01*k10;
 921 // //   k21+=h01*k20;
 922 // //   k31+=h01*k30;
 923 // //   k41+=h01*k40;
 924
 925
 926 //   fY += k00*dy + k01*dz;
 927 //   fZ += k10*dy + k11*dz;
 928 //   fE  = eta;
 929 //   fT += k30*dy + k31*dz;
 930 //   fC  = cur;
 931
 932
 933 //   //Update covariance
 934 //   //
 935 //   //
 936 //   Double_t oldyy = fCyy, oldzz = fCzz; //, oldee=fCee, oldcc =fCcc;
 937 //   Double_t oldzy = fCzy, oldey = fCey, oldty=fCty, oldcy =fCcy;
 938 //   Double_t oldez = fCez, oldtz = fCtz, oldcz=fCcz;
 939 //   //Double_t oldte = fCte, oldce = fCce;
 940 //   //Double_t oldct = fCct;
 941
 942 //   fCyy-=k00*oldyy+k01*oldzy;
 943 //   fCzy-=k10*oldyy+k11*oldzy;
 944 //   fCey-=k20*oldyy+k21*oldzy;
 945 //   fCty-=k30*oldyy+k31*oldzy;
 946 //   fCcy-=k40*oldyy+k41*oldzy;
 947 //   //
 948 //   fCzz-=k10*oldzy+k11*oldzz;
 949 //   fCez-=k20*oldzy+k21*oldzz;
 950 //   fCtz-=k30*oldzy+k31*oldzz;
 951 //   fCcz-=k40*oldzy+k41*oldzz;
 952 //   //
 953 //   fCee-=k20*oldey+k21*oldez;
 954 //   fCte-=k30*oldey+k31*oldez;
 955 //   fCce-=k40*oldey+k41*oldez;
 956 //   //
 957 //   fCtt-=k30*oldty+k31*oldtz;
 958 //   fCct-=k40*oldty+k41*oldtz;
 959 //   //
 960
 961 //   fCcc-=k40*oldcy+k41*oldcz;
 962 //   //
 963
 964 //   //Int_t n=GetNumberOfClusters();
 965 //   //fIndex[n]=index;
 966 //   //SetNumberOfClusters(n+1);
 967
 968 //   //SetChi2(GetChi2()+chisq);
 969 //   //  cerr<<"in update: fIndex["<<fN<<"] = "<<index<<endl;
 970
 971 //   return 1;
 972
 973 // }
 974
 975 //_____________________________________________________________________________
 976 Bool_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
 977 {
 978   //
 979   // Rotates track parameters in R*phi plane
 980   // if absolute rotation alpha is in global system
 981   // otherwise alpha rotation is relative to the current rotation angle
 982   //
 983
 984   if (absolute) {
 985     alpha -= GetAlpha();
 986   }
 987   else{
 988     fNRotate++;
 989   }
 990
 991   return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
 992
 993 }
 994
 995 //_____________________________________________________________________________
 996 Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
 997 {
 998   //
 999   // Returns the track chi2
1000   //
1001
1002   Double_t p[2]   = { c->GetY(), c->GetZ() };
1003   Double_t sy2    = c->GetSigmaY2() * 4.0;
1004   Double_t sz2    = c->GetSigmaZ2() * 4.0;
1005   Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
1006
1007   return AliExternalTrackParam::GetPredictedChi2(p,cov);
1008
1009   //
1010   // Can the following be removed ????
1011   //
1012   /*
1013   Bool_t fNoTilt = kTRUE;
1014   if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE;
1015
1016   return (c->GetY() - GetY())*(c->GetY() - GetY())/c->GetSigmaY2();
1017   */
1018
1019   /*
1020   Double_t chi2, dy, r00, r01, r11;
1021
1022   if(fNoTilt) {
1023     dy=c->GetY() - fY;
1024     r00=c->GetSigmaY2();
1025     chi2 = (dy*dy)/r00;
1026   }
1027   else {
1028     Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12);
1029     //
1030     r00=c->GetSigmaY2(); r01=0.; r11=c->GetSigmaZ2();
1031     r00+=fCyy; r01+=fCzy; r11+=fCzz;
1032
1033     Double_t det=r00*r11 - r01*r01;
1034     if (TMath::Abs(det) < 1.e-10) {
1035       Int_t n=GetNumberOfClusters();
1036       if (n>4) cerr<<n<<" AliTRDtrack warning: Singular matrix !\n";
1037       return 1e10;
1038     }
1039     Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
1040     Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
1041     Double_t tiltdz = dz;
1042     if (TMath::Abs(tiltdz)>padlength/2.) {
1043       tiltdz = TMath::Sign(padlength/2,dz);
1044     }
1045     //    dy=dy+h01*dz;
1046     dy=dy+h01*tiltdz;
1047
1048     chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
1049   }
1050
1051   return chi2;
1052   */
1053
1054 }
1055
1056 //_____________________________________________________________________________
1057 void AliTRDtrack::MakeBackupTrack()
1058 {
1059   //
1060   // Creates a backup track
1061   //
1062
1063   if (fBackupTrack) {
1064     delete fBackupTrack;
1065   }
1066   fBackupTrack = new AliTRDtrack(*this);
1067
1068 }
1069
1070 //_____________________________________________________________________________
1071 Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
1072 {
1073   //
1074   // Find a prolongation at given x
1075   // Return 0 if it does not exist
1076   //
1077
1078   Double_t bz = GetBz();
1079
1080   if (!AliExternalTrackParam::GetYAt(xk,bz,y)) {
1081     return 0;
1082   }
1083   if (!AliExternalTrackParam::GetZAt(xk,bz,z)) {
1084     return 0;
1085   }
1086
1087   return 1;
1088
1089 }
1090
1091 //_____________________________________________________________________________
1092 Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
1093 {
1094   //
1095   // Propagate track to given x  position
1096   // Works inside of the 20 degree segmentation
1097   // (local cooordinate frame for TRD , TPC, TOF)
1098   //
1099   // Material budget from geo manager
1100   //
1101
1102   Double_t xyz0[3];
1103   Double_t xyz1[3];
1104   Double_t y;
1105   Double_t z;
1106
1107   const Double_t kAlphac  = TMath::Pi()/9.0;
1108   const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
1109
1110   // Critical alpha  - cross sector indication
1111   Double_t dir = (GetX() > xr) ? -1.0 : 1.0;
1112
1113   // Direction +-
1114   for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) {
1115
1116     GetXYZ(xyz0);
1117     GetProlongation(x,y,z);
1118     xyz1[0] = x * TMath::Cos(GetAlpha()) + y * TMath::Sin(GetAlpha());
1119     xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha());
1120     xyz1[2] = z;
1121     Double_t param[7];
1122     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1123
1124     if ((param[0] > 0) &&
1125         (param[1] > 0)) {
1126       PropagateTo(x,param[1],param[0]*param[4]);
1127     }
1128
1129     if (GetY() >  GetX()*kTalphac) {
1130       Rotate(-kAlphac);
1131     }
1132     if (GetY() < -GetX()*kTalphac) {
1133       Rotate( kAlphac);
1134     }
1135
1136   }
1137
1138   PropagateTo(xr);
1139
1140   return 0;
1141
1142 }
1143
1144 //_____________________________________________________________________________
1145 Int_t   AliTRDtrack::PropagateToR(Double_t r,Double_t step)
1146 {
1147   //
1148   // Propagate track to the radial position
1149   // Rotation always connected to the last track position
1150   //
1151
1152   Double_t xyz0[3];
1153   Double_t xyz1[3];
1154   Double_t y;
1155   Double_t z;
1156
1157   Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
1158   // Direction +-
1159   Double_t dir    = (radius > r) ? -1.0 : 1.0;
1160
1161   for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
1162
1163     GetXYZ(xyz0);
1164     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1165     Rotate(alpha,kTRUE);
1166     GetXYZ(xyz0);
1167     GetProlongation(x,y,z);
1168     xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1169     xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1170     xyz1[2] = z;
1171     Double_t param[7];
1172     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1173     if (param[1] <= 0) {
1174       param[1] = 100000000;
1175     }
1176     PropagateTo(x,param[1],param[0]*param[4]);
1177
1178   }
1179
1180   GetXYZ(xyz0);
1181   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1182   Rotate(alpha,kTRUE);
1183   GetXYZ(xyz0);
1184   GetProlongation(r,y,z);
1185   xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1186   xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1187   xyz1[2] = z;
1188   Double_t param[7];
1189   AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1190
1191   if (param[1] <= 0) {
1192     param[1] = 100000000;
1193   }
1194   PropagateTo(r,param[1],param[0]*param[4]);
1195
1196   return 0;
1197
1198 }
1199
1200 //_____________________________________________________________________________
1201 Int_t AliTRDtrack::GetSector() const
1202 {
1203   //
1204   // Return the current sector
1205   //
1206
1207   return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha())
1208              % AliTRDgeometry::kNsect;
1209
1210 }
1211
1212 //_____________________________________________________________________________
1213 void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
1214 {
1215   //
1216   // The sampled energy loss
1217   //
1218
1219   Double_t s = GetSnp();
1220   Double_t t = GetTgl();
1221   q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1222   fdQdl[i] = q;
1223
1224 }
1225
1226 //_____________________________________________________________________________
1227 void AliTRDtrack::SetSampledEdx(Float_t q)
1228 {
1229   //
1230   // The sampled energy loss
1231   //
1232
1233   Double_t s = GetSnp();
1234   Double_t t = GetTgl();
1235   q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1236   fdQdl[fNdedx] = q;
1237   fNdedx++;
1238
1239 }
1240
1241 //_____________________________________________________________________________
1242 Double_t AliTRDtrack::GetBz() const
1243 {
1244   //
1245   // Returns Bz component of the magnetic field (kG)
1246   //
1247
1248   if (AliTracker::UniformField()) {
1249     return AliTracker::GetBz();
1250   }
1251   Double_t r[3];
1252   GetXYZ(r);
1253
1254   return AliTracker::GetBz(r);
1255
1256 }