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/AliTRDCalPIDLQ.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) // A.Bercuci
  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) // A.Bercuci
 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) // A.Bercuci
 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) // A.Bercuci
 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) // A.Bercuci
 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 AliTRDCalPIDLQ and AliTRDCalPIDLQRef 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 AliTRDCalPIDLQ *pd = calibration->GetPIDLQObject();
 596         if (!pd) {
 597                 AliError("No access to AliTRDCalPIDLQ");
 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
 611         // Skip tracks which have no TRD signal at all
 612         if (fdEdx == 0.) return -1;
 613
 614         for (Int_t iPlane = 0; iPlane < AliTRDgeometry::kNplan; iPlane++) {
 615
 616                 length = (AliTRDgeometry::AmThick() + AliTRDgeometry::DrThick())/TMath::Sqrt((1. - fSnp[iPlane]*fSnp[iPlane]) / (1. + fTgl[iPlane]*fTgl[iPlane]));
 617
 618                 // check data
 619                 if((fdEdxPlane[iPlane][0] + fdEdxPlane[iPlane][1] + fdEdxPlane[iPlane][2]) <=  0.
 620                 || fTimBinPlane[iPlane] == -1.) continue;
 621
 622
 623                 // this track segment has fulfilled all requierments
 624                 nPlanePID++;
 625
 626                 // Get the probabilities for the different particle species
 627                 for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) {
 628                         p[iSpecies] *= pd->GetProbability(iSpecies, fMom[iPlane], fdEdxPlane[iPlane], length);
 629                         //p[iSpecies] *= pd->GetProbabilityT(iSpecies, fMom[iPlane], fTimBinPlane[iPlane]);
 630                 }
 631         }
 632         if(nPlanePID == 0) return 0;
 633
 634         // normalize probabilities
 635         Double_t probTotal = 0.;
 636         for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) probTotal += p[iSpecies];
 637
 638         if(probTotal <= 0.) {
 639                 AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference histograms.");
 640                 return -1;
 641         }
 642         for(Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) p[iSpecies] /= probTotal;
 643
 644
 645         // book PID to the ESD track
 646         esd->SetTRDpid(p);
 647         esd->SetTRDpidQuality(nPlanePID);
 648
 649         return 0;
 650 }
 651
 652
 653
 654 //_____________________________________________________________________________
 655 Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t x0, Double_t rho)
 656 {
 657   //
 658   // Propagates a track of particle with mass=pm to a reference plane
 659   // defined by x=xk through media of density=rho and radiationLength=x0
 660   //
 661
 662   if (xk == GetX()) {
 663     return kTRUE;
 664   }
 665
 666   Double_t oldX = GetX();
 667   Double_t oldY = GetY();
 668   Double_t oldZ = GetZ();
 669
 670   Double_t bz   = GetBz();
 671
 672   if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
 673     return kFALSE;
 674   }
 675
 676   Double_t x = GetX();
 677   Double_t y = GetY();
 678   Double_t z = GetZ();
 679   Double_t d = TMath::Sqrt((x-oldX)*(x-oldX)
 680                          + (y-oldY)*(y-oldY)
 681                          + (z-oldZ)*(z-oldZ));
 682   if (oldX < xk) {
 683     if (IsStartedTimeIntegral()) {
 684       Double_t l2  = d;
 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   Double_t ll = (oldX < xk) ? -d : d;
 697   if (!AliExternalTrackParam::CorrectForMaterial(ll*rho/x0,x0,GetMass())) {
 698     return kFALSE;
 699   }
 700
 701   {
 702
 703     // Energy losses************************
 704     Double_t p2    = (1.0 + GetTgl()*GetTgl()) / (Get1Pt()*Get1Pt());
 705     Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
 706     if (beta2<1.0e-10 || (5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0) {
 707       return kFALSE;
 708     }
 709
 710     Double_t dE    = 0.153e-3 / beta2
 711                    * (log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
 712                    * d * rho;
 713     Float_t budget = d * rho;
 714     fBudget[0] += budget;
 715
 716     /*
 717     // Suspicious part - think about it ?
 718     Double_t kinE =  TMath::Sqrt(p2);
 719     if (dE > 0.8*kinE) dE = 0.8 * kinE;  //
 720     if (dE < 0)        dE = 0.0;         // Not valid region for Bethe bloch
 721     */
 722
 723     fDE += dE;
 724
 725     /*
 726     // Suspicious ! I.B.
 727     Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE));   // Energy loss fluctuation
 728     Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
 729     fCcc += sigmac2;
 730     fCee += fX*fX * sigmac2;
 731     */
 732
 733   }
 734
 735   return kTRUE;
 736
 737 }
 738
 739 //_____________________________________________________________________________
 740 Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index
 741                          , Double_t h01)
 742 {
 743   //
 744   // Assignes found cluster to the track and updates track information
 745   //
 746
 747   Bool_t fNoTilt = kTRUE;
 748   if (TMath::Abs(h01) > 0.003) {
 749     fNoTilt = kFALSE;
 750   }
 751
 752   // Add angular effect to the error contribution -  MI
 753   Float_t tangent2 = GetSnp()*GetSnp();
 754   if (tangent2 < 0.90000) {
 755     tangent2 = tangent2 / (1.0 - tangent2);
 756   }
 757   //Float_t errang = tangent2 * 0.04;
 758
 759   Double_t p[2]   = {c->GetY(), c->GetZ() };
 760   //Double_t cov[3] = {c->GetSigmaY2()+errang, 0.0, c->GetSigmaZ2()*100.0 };
 761   Double_t sy2    = c->GetSigmaY2() * 4.0;
 762   Double_t sz2    = c->GetSigmaZ2() * 4.0;
 763   Double_t cov[3] = {sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
 764
 765   if (!AliExternalTrackParam::Update(p,cov)) {
 766     return kFALSE;
 767   }
 768
 769   Int_t n = GetNumberOfClusters();
 770   fIndex[n] = index;
 771   SetNumberOfClusters(n+1);
 772
 773
 774   SetChi2(GetChi2()+chisq);
 775
 776   return kTRUE;
 777
 778 }
 779
 780 //_____________________________________________________________________________
 781 Int_t AliTRDtrack::UpdateMI(AliTRDcluster *c, Double_t chisq, Int_t index, Double_t h01, Int_t /*plane*/)
 782 {
 783   //
 784   // Assignes found cluster to the track and updates track information
 785   //
 786
 787   Bool_t fNoTilt = kTRUE;
 788   if (TMath::Abs(h01) > 0.003) {
 789     fNoTilt = kFALSE;
 790   }
 791
 792   // Add angular effect to the error contribution and make correction  -  MI
 793   Double_t tangent2 = GetSnp()*GetSnp();
 794   if (tangent2 < 0.90000){
 795     tangent2 = tangent2 / (1.0-tangent2);
 796   }
 797   Double_t tangent = TMath::Sqrt(tangent2);
 798   if (GetSnp() < 0) {
 799     tangent *= -1;
 800   }
 801
 802   //
 803   // Is the following still needed ????
 804   //
 805
 806   //  Double_t correction = 0*plane;
 807   /*
 808   Double_t errang = tangent2*0.04;  //
 809   Double_t errsys =0.025*0.025*20;  //systematic error part
 810
 811   Float_t extend =1;
 812   if (c->GetNPads()==4) extend=2;
 813   */
 814   //if (c->GetNPads()==5)  extend=3;
 815   //if (c->GetNPads()==6)  extend=3;
 816   //if (c->GetQ()<15) return 1;
 817
 818   /*
 819   if (corrector!=0){
 820   //if (0){
 821     correction = corrector->GetCorrection(plane,c->GetLocalTimeBin(),tangent);
 822     if (TMath::Abs(correction)>0){
 823       //if we have info
 824       errang     = corrector->GetSigma(plane,c->GetLocalTimeBin(),tangent);
 825       errang    *= errang;
 826       errang    += tangent2*0.04;
 827     }
 828   }
 829   */
 830   //
 831   //Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
 832   /*
 833     {
 834       Double_t dy=c->GetY() - GetY(), dz=c->GetZ() - GetZ();
 835       printf("%e %e %e %e\n",dy,dz,padlength/2,h01);
 836     }
 837   */
 838
 839   Double_t p[2]   = { c->GetY(), c->GetZ() };
 840   //Double_t cov[3]={ (c->GetSigmaY2()+errang+errsys)*extend, 0.0, c->GetSigmaZ2()*10000.0 };
 841   Double_t sy2    = c->GetSigmaY2() * 4.0;
 842   Double_t sz2    = c->GetSigmaZ2() * 4.0;
 843   Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
 844
 845   if (!AliExternalTrackParam::Update(p,cov)) {
 846     return kFALSE;
 847   }
 848
 849         // Register cluster to track
 850   Int_t n = GetNumberOfClusters();
 851   fIndex[n] = index;
 852         fClusters[n] = c; // A.Bercuci 25.07.07
 853   SetNumberOfClusters(n+1);
 854   SetChi2(GetChi2() + chisq);
 855
 856   return kTRUE;
 857
 858 }
 859
 860 // //_____________________________________________________________________________
 861 // Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet)
 862 // {
 863 //   //
 864 //   // Assignes found tracklet to the track and updates track information
 865 //   //
 866 //   // Can this be removed ????
 867 //   //
 868 //
 869 //   Double_t r00=(tracklet.GetTrackletSigma2()), r01=0., r11= 10000.;
 870 //   r00+=fCyy; r01+=fCzy; r11+=fCzz;
 871 //   //
 872 //   Double_t det=r00*r11 - r01*r01;
 873 //   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 874 //   //
 875
 876 //   Double_t dy=tracklet.GetY() - fY, dz=tracklet.GetZ() - fZ;
 877
 878
 879 //   Double_t s00 = tracklet.GetTrackletSigma2();  // error pad
 880 //   Double_t s11 = 100000;   // error pad-row
 881 //   Float_t  h01 = tracklet.GetTilt();
 882 //   //
 883 //   //  r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00;
 884 //   r00 = fCyy + fCzz*h01*h01+s00;
 885 //   //  r01 = fCzy + fCzz*h01;
 886 //   r01 = fCzy ;
 887 //   r11 = fCzz + s11;
 888 //   det = r00*r11 - r01*r01;
 889 //   // inverse matrix
 890 //   tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 891
 892 //   Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11;
 893 //   Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11;
 894 //   Double_t k20=fCey*r00+fCez*r01, k21=fCey*r01+fCez*r11;
 895 //   Double_t k30=fCty*r00+fCtz*r01, k31=fCty*r01+fCtz*r11;
 896 //   Double_t k40=fCcy*r00+fCcz*r01, k41=fCcy*r01+fCcz*r11;
 897
 898 //   // K matrix
 899 // //   k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01);
 900 // //   k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01);
 901 // //   k20=fCey*r00+fCez*(r01+h01*r00),k21=fCey*r01+fCez*(r11+h01*r01);
 902 // //   k30=fCty*r00+fCtz*(r01+h01*r00),k31=fCty*r01+fCtz*(r11+h01*r01);
 903 // //   k40=fCcy*r00+fCcz*(r01+h01*r00),k41=fCcy*r01+fCcz*(r11+h01*r01);
 904 //   //
 905 //   //Update measurement
 906 //   Double_t cur=fC + k40*dy + k41*dz, eta=fE + k20*dy + k21*dz;
 907 //   //  cur=fC + k40*dy + k41*dz; eta=fE + k20*dy + k21*dz;
 908 //   if (TMath::Abs(cur*fX-eta) >= 0.90000) {
 909 //     //Int_t n=GetNumberOfClusters();
 910 //     //      if (n>4) cerr<<n<<" AliTRDtrack warning: Filtering failed !\n";
 911 //     return 0;
 912 //   }
 913 // //   k01+=h01*k00;
 914 // //   k11+=h01*k10;
 915 // //   k21+=h01*k20;
 916 // //   k31+=h01*k30;
 917 // //   k41+=h01*k40;
 918
 919
 920 //   fY += k00*dy + k01*dz;
 921 //   fZ += k10*dy + k11*dz;
 922 //   fE  = eta;
 923 //   fT += k30*dy + k31*dz;
 924 //   fC  = cur;
 925
 926
 927 //   //Update covariance
 928 //   //
 929 //   //
 930 //   Double_t oldyy = fCyy, oldzz = fCzz; //, oldee=fCee, oldcc =fCcc;
 931 //   Double_t oldzy = fCzy, oldey = fCey, oldty=fCty, oldcy =fCcy;
 932 //   Double_t oldez = fCez, oldtz = fCtz, oldcz=fCcz;
 933 //   //Double_t oldte = fCte, oldce = fCce;
 934 //   //Double_t oldct = fCct;
 935
 936 //   fCyy-=k00*oldyy+k01*oldzy;
 937 //   fCzy-=k10*oldyy+k11*oldzy;
 938 //   fCey-=k20*oldyy+k21*oldzy;
 939 //   fCty-=k30*oldyy+k31*oldzy;
 940 //   fCcy-=k40*oldyy+k41*oldzy;
 941 //   //
 942 //   fCzz-=k10*oldzy+k11*oldzz;
 943 //   fCez-=k20*oldzy+k21*oldzz;
 944 //   fCtz-=k30*oldzy+k31*oldzz;
 945 //   fCcz-=k40*oldzy+k41*oldzz;
 946 //   //
 947 //   fCee-=k20*oldey+k21*oldez;
 948 //   fCte-=k30*oldey+k31*oldez;
 949 //   fCce-=k40*oldey+k41*oldez;
 950 //   //
 951 //   fCtt-=k30*oldty+k31*oldtz;
 952 //   fCct-=k40*oldty+k41*oldtz;
 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     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1115
1116     if ((param[0] > 0) &&
1117         (param[1] > 0)) {
1118       PropagateTo(x,param[1],param[0]);
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     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1165     if (param[1] <= 0) {
1166       param[1] = 100000000;
1167     }
1168     PropagateTo(x,param[1],param[0]);
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   AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1182
1183   if (param[1] <= 0) {
1184     param[1] = 100000000;
1185   }
1186   PropagateTo(r,param[1],param[0]);
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 }