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