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