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 x0, Double_t rho)
 654 {
 655   //
 656   // Propagates a track of particle with mass=pm to a reference plane
 657   // defined by x=xk through media of density=rho and radiationLength=x0
 658   //
 659
 660   if (xk == GetX()) {
 661     return kTRUE;
 662   }
 663
 664   Double_t oldX = GetX();
 665   Double_t oldY = GetY();
 666   Double_t oldZ = GetZ();
 667
 668   Double_t bz   = GetBz();
 669
 670   if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
 671     return kFALSE;
 672   }
 673
 674   Double_t x = GetX();
 675   Double_t y = GetY();
 676   Double_t z = GetZ();
 677   Double_t d = TMath::Sqrt((x-oldX)*(x-oldX)
 678                          + (y-oldY)*(y-oldY)
 679                          + (z-oldZ)*(z-oldZ));
 680   if (oldX < xk) {
 681     if (IsStartedTimeIntegral()) {
 682       Double_t l2  = d;
 683       Double_t crv = GetC();
 684       if (TMath::Abs(l2*crv) > 0.0001) {
 685         // Make correction for curvature if neccesary
 686         l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY));
 687         l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
 688         l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ));
 689       }
 690       AddTimeStep(l2);
 691     }
 692   }
 693
 694   Double_t ll = (oldX < xk) ? -d : d;
 695   if (!AliExternalTrackParam::CorrectForMaterial(ll*rho/x0,x0,GetMass())) {
 696     return kFALSE;
 697   }
 698
 699   {
 700
 701     // Energy losses************************
 702     Double_t p2    = (1.0 + GetTgl()*GetTgl()) / (Get1Pt()*Get1Pt());
 703     Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
 704     if (beta2<1.0e-10 || (5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0) {
 705       return kFALSE;
 706     }
 707
 708     Double_t dE    = 0.153e-3 / beta2
 709                    * (log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
 710                    * d * rho;
 711     Float_t budget = d * rho;
 712     fBudget[0] += budget;
 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 // Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet)
 860 // {
 861 //   //
 862 //   // Assignes found tracklet to the track and updates track information
 863 //   //
 864 //   // Can this be removed ????
 865 //   //
 866 //
 867 //   Double_t r00=(tracklet.GetTrackletSigma2()), r01=0., r11= 10000.;
 868 //   r00+=fCyy; r01+=fCzy; r11+=fCzz;
 869 //   //
 870 //   Double_t det=r00*r11 - r01*r01;
 871 //   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 872 //   //
 873
 874 //   Double_t dy=tracklet.GetY() - fY, dz=tracklet.GetZ() - fZ;
 875
 876
 877 //   Double_t s00 = tracklet.GetTrackletSigma2();  // error pad
 878 //   Double_t s11 = 100000;   // error pad-row
 879 //   Float_t  h01 = tracklet.GetTilt();
 880 //   //
 881 //   //  r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00;
 882 //   r00 = fCyy + fCzz*h01*h01+s00;
 883 //   //  r01 = fCzy + fCzz*h01;
 884 //   r01 = fCzy ;
 885 //   r11 = fCzz + s11;
 886 //   det = r00*r11 - r01*r01;
 887 //   // inverse matrix
 888 //   tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 889
 890 //   Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11;
 891 //   Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11;
 892 //   Double_t k20=fCey*r00+fCez*r01, k21=fCey*r01+fCez*r11;
 893 //   Double_t k30=fCty*r00+fCtz*r01, k31=fCty*r01+fCtz*r11;
 894 //   Double_t k40=fCcy*r00+fCcz*r01, k41=fCcy*r01+fCcz*r11;
 895
 896 //   // K matrix
 897 // //   k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01);
 898 // //   k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01);
 899 // //   k20=fCey*r00+fCez*(r01+h01*r00),k21=fCey*r01+fCez*(r11+h01*r01);
 900 // //   k30=fCty*r00+fCtz*(r01+h01*r00),k31=fCty*r01+fCtz*(r11+h01*r01);
 901 // //   k40=fCcy*r00+fCcz*(r01+h01*r00),k41=fCcy*r01+fCcz*(r11+h01*r01);
 902 //   //
 903 //   //Update measurement
 904 //   Double_t cur=fC + k40*dy + k41*dz, eta=fE + k20*dy + k21*dz;
 905 //   //  cur=fC + k40*dy + k41*dz; eta=fE + k20*dy + k21*dz;
 906 //   if (TMath::Abs(cur*fX-eta) >= 0.90000) {
 907 //     //Int_t n=GetNumberOfClusters();
 908 //     //      if (n>4) cerr<<n<<" AliTRDtrack warning: Filtering failed !\n";
 909 //     return 0;
 910 //   }
 911 // //   k01+=h01*k00;
 912 // //   k11+=h01*k10;
 913 // //   k21+=h01*k20;
 914 // //   k31+=h01*k30;
 915 // //   k41+=h01*k40;
 916
 917
 918 //   fY += k00*dy + k01*dz;
 919 //   fZ += k10*dy + k11*dz;
 920 //   fE  = eta;
 921 //   fT += k30*dy + k31*dz;
 922 //   fC  = cur;
 923
 924
 925 //   //Update covariance
 926 //   //
 927 //   //
 928 //   Double_t oldyy = fCyy, oldzz = fCzz; //, oldee=fCee, oldcc =fCcc;
 929 //   Double_t oldzy = fCzy, oldey = fCey, oldty=fCty, oldcy =fCcy;
 930 //   Double_t oldez = fCez, oldtz = fCtz, oldcz=fCcz;
 931 //   //Double_t oldte = fCte, oldce = fCce;
 932 //   //Double_t oldct = fCct;
 933
 934 //   fCyy-=k00*oldyy+k01*oldzy;
 935 //   fCzy-=k10*oldyy+k11*oldzy;
 936 //   fCey-=k20*oldyy+k21*oldzy;
 937 //   fCty-=k30*oldyy+k31*oldzy;
 938 //   fCcy-=k40*oldyy+k41*oldzy;
 939 //   //
 940 //   fCzz-=k10*oldzy+k11*oldzz;
 941 //   fCez-=k20*oldzy+k21*oldzz;
 942 //   fCtz-=k30*oldzy+k31*oldzz;
 943 //   fCcz-=k40*oldzy+k41*oldzz;
 944 //   //
 945 //   fCee-=k20*oldey+k21*oldez;
 946 //   fCte-=k30*oldey+k31*oldez;
 947 //   fCce-=k40*oldey+k41*oldez;
 948 //   //
 949 //   fCtt-=k30*oldty+k31*oldtz;
 950 //   fCct-=k40*oldty+k41*oldtz;
 951 //   //
 952 //   fCcc-=k40*oldcy+k41*oldcz;
 953 //   //
 954
 955 //   //Int_t n=GetNumberOfClusters();
 956 //   //fIndex[n]=index;
 957 //   //SetNumberOfClusters(n+1);
 958
 959 //   //SetChi2(GetChi2()+chisq);
 960 //   //  cerr<<"in update: fIndex["<<fN<<"] = "<<index<<endl;
 961
 962 //   return 1;
 963
 964 // }
 965
 966 //_____________________________________________________________________________
 967 Bool_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
 968 {
 969   //
 970   // Rotates track parameters in R*phi plane
 971   // if absolute rotation alpha is in global system
 972   // otherwise alpha rotation is relative to the current rotation angle
 973   //
 974
 975   if (absolute) {
 976     alpha -= GetAlpha();
 977   }
 978   else{
 979     fNRotate++;
 980   }
 981
 982   return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
 983
 984 }
 985
 986 //_____________________________________________________________________________
 987 Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
 988 {
 989   //
 990   // Returns the track chi2
 991   //
 992
 993   Double_t p[2]   = { c->GetY(), c->GetZ() };
 994   Double_t sy2    = c->GetSigmaY2() * 4.0;
 995   Double_t sz2    = c->GetSigmaZ2() * 4.0;
 996   Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
 997
 998   return AliExternalTrackParam::GetPredictedChi2(p,cov);
 999
1000   //
1001   // Can the following be removed ????
1002   //
1003   /*
1004   Bool_t fNoTilt = kTRUE;
1005   if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE;
1006
1007   return (c->GetY() - GetY())*(c->GetY() - GetY())/c->GetSigmaY2();
1008   */
1009
1010   /*
1011   Double_t chi2, dy, r00, r01, r11;
1012
1013   if(fNoTilt) {
1014     dy=c->GetY() - fY;
1015     r00=c->GetSigmaY2();
1016     chi2 = (dy*dy)/r00;
1017   }
1018   else {
1019     Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12);
1020     //
1021     r00=c->GetSigmaY2(); r01=0.; r11=c->GetSigmaZ2();
1022     r00+=fCyy; r01+=fCzy; r11+=fCzz;
1023
1024     Double_t det=r00*r11 - r01*r01;
1025     if (TMath::Abs(det) < 1.e-10) {
1026       Int_t n=GetNumberOfClusters();
1027       if (n>4) cerr<<n<<" AliTRDtrack warning: Singular matrix !\n";
1028       return 1e10;
1029     }
1030     Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
1031     Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
1032     Double_t tiltdz = dz;
1033     if (TMath::Abs(tiltdz)>padlength/2.) {
1034       tiltdz = TMath::Sign(padlength/2,dz);
1035     }
1036     //    dy=dy+h01*dz;
1037     dy=dy+h01*tiltdz;
1038
1039     chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
1040   }
1041
1042   return chi2;
1043   */
1044
1045 }
1046
1047 //_____________________________________________________________________________
1048 void AliTRDtrack::MakeBackupTrack()
1049 {
1050   //
1051   // Creates a backup track
1052   //
1053
1054   if (fBackupTrack) {
1055     delete fBackupTrack;
1056   }
1057   fBackupTrack = new AliTRDtrack(*this);
1058
1059 }
1060
1061 //_____________________________________________________________________________
1062 Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
1063 {
1064   //
1065   // Find a prolongation at given x
1066   // Return 0 if it does not exist
1067   //
1068
1069   Double_t bz = GetBz();
1070
1071   if (!AliExternalTrackParam::GetYAt(xk,bz,y)) {
1072     return 0;
1073   }
1074   if (!AliExternalTrackParam::GetZAt(xk,bz,z)) {
1075     return 0;
1076   }
1077
1078   return 1;
1079
1080 }
1081
1082 //_____________________________________________________________________________
1083 Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
1084 {
1085   //
1086   // Propagate track to given x  position
1087   // Works inside of the 20 degree segmentation (local cooordinate frame for TRD , TPC, TOF)
1088   //
1089   // Material budget from geo manager
1090   //
1091
1092   Double_t xyz0[3];
1093   Double_t xyz1[3];
1094   Double_t y;
1095   Double_t z;
1096
1097   const Double_t kAlphac  = TMath::Pi()/9.0;
1098   const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
1099
1100   // Critical alpha  - cross sector indication
1101   Double_t dir = (GetX()>xr) ? -1.0 : 1.0;
1102
1103   // Direction +-
1104   for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) {
1105
1106     GetXYZ(xyz0);
1107     GetProlongation(x,y,z);
1108     xyz1[0] = x * TMath::Cos(GetAlpha()) + y * TMath::Sin(GetAlpha());
1109     xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha());
1110     xyz1[2] = z;
1111     Double_t param[7];
1112     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1113
1114     if ((param[0] > 0) &&
1115         (param[1] > 0)) {
1116       PropagateTo(x,param[1],param[0]);
1117     }
1118
1119     if (GetY() >  GetX()*kTalphac) {
1120       Rotate(-kAlphac);
1121     }
1122     if (GetY() < -GetX()*kTalphac) {
1123       Rotate( kAlphac);
1124     }
1125
1126   }
1127
1128   PropagateTo(xr);
1129
1130   return 0;
1131
1132 }
1133
1134 //_____________________________________________________________________________
1135 Int_t   AliTRDtrack::PropagateToR(Double_t r,Double_t step)
1136 {
1137   //
1138   // Propagate track to the radial position
1139   // Rotation always connected to the last track position
1140   //
1141
1142   Double_t xyz0[3];
1143   Double_t xyz1[3];
1144   Double_t y;
1145   Double_t z;
1146
1147   Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
1148   // Direction +-
1149   Double_t dir    = (radius>r) ? -1.0 : 1.0;
1150
1151   for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
1152
1153     GetXYZ(xyz0);
1154     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1155     Rotate(alpha,kTRUE);
1156     GetXYZ(xyz0);
1157     GetProlongation(x,y,z);
1158     xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1159     xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1160     xyz1[2] = z;
1161     Double_t param[7];
1162     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1163     if (param[1] <= 0) {
1164       param[1] = 100000000;
1165     }
1166     PropagateTo(x,param[1],param[0]);
1167
1168   }
1169
1170   GetXYZ(xyz0);
1171   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1172   Rotate(alpha,kTRUE);
1173   GetXYZ(xyz0);
1174   GetProlongation(r,y,z);
1175   xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1176   xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1177   xyz1[2] = z;
1178   Double_t param[7];
1179   AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1180
1181   if (param[1] <= 0) {
1182     param[1] = 100000000;
1183   }
1184   PropagateTo(r,param[1],param[0]);
1185
1186   return 0;
1187
1188 }
1189
1190 //_____________________________________________________________________________
1191 Int_t AliTRDtrack::GetSector() const
1192 {
1193   //
1194   // Return the current sector
1195   //
1196
1197   return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha())
1198          % AliTRDgeometry::kNsect;
1199
1200 }
1201
1202 //_____________________________________________________________________________
1203 void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
1204 {
1205   //
1206   // The sampled energy loss
1207   //
1208
1209   Double_t s = GetSnp();
1210   Double_t t = GetTgl();
1211   q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1212   fdQdl[i] = q;
1213
1214 }
1215
1216 //_____________________________________________________________________________
1217 void AliTRDtrack::SetSampledEdx(Float_t q)
1218 {
1219   //
1220   // The sampled energy loss
1221   //
1222
1223   Double_t s = GetSnp();
1224   Double_t t = GetTgl();
1225   q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1226   fdQdl[fNdedx] = q;
1227   fNdedx++;
1228
1229 }
1230
1231 //_____________________________________________________________________________
1232 Double_t AliTRDtrack::GetBz() const
1233 {
1234   //
1235   // Returns Bz component of the magnetic field (kG)
1236   //
1237
1238   if (AliTracker::UniformField()) {
1239     return AliTracker::GetBz();
1240   }
1241   Double_t r[3];
1242   GetXYZ(r);
1243
1244   return AliTracker::GetBz(r);
1245
1246 }