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