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