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   AliTracker::FillResiduals(this,p,cov,c->GetVolumeId());
 880
 881   // Register cluster to track
 882   Int_t n      = GetNumberOfClusters();
 883   fIndex[n]    = index;
 884   fClusters[n] = c;
 885   SetNumberOfClusters(n+1);
 886   SetChi2(GetChi2() + chisq);
 887
 888   return kTRUE;
 889
 890 }
 891
 892 //_____________________________________________________________________________
 893 Bool_t AliTRDtrack::Update(const AliTRDtracklet &t, Double_t chisq, Int_t index)
 894 {
 895   //
 896   // Assignes the found tracklet <t> to the track
 897   // and updates track information
 898   // <chisq> : predicted chi2
 899   // <index> : ??
 900   //
 901
 902   Double_t h01    = t.GetTilt(); // Is this correct????
 903   Double_t p[2]   = { t.GetY(), t.GetZ() };
 904   Double_t sy2    = t.GetTrackletSigma2();  // Error pad-column
 905   Double_t sz2    = 100000.0;               // Error pad-row (????)
 906   Double_t cov[3] = { sy2 + h01*h01*sz2  // Does this have any sense????
 907                     , h01 * (sy2 - sz2)
 908                     , sz2 + h01*h01*sy2 };
 909   if (!AliExternalTrackParam::Update(p,cov)) {
 910     return kFALSE;
 911   }
 912
 913   Int_t n   = GetNumberOfClusters();
 914   fIndex[n] = index;
 915   SetNumberOfClusters(n+1);
 916   SetChi2(GetChi2()+chisq);
 917
 918   return kTRUE;
 919
 920 }
 921
 922 //_____________________________________________________________________________
 923 Bool_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
 924 {
 925   //
 926   // Rotates track parameters in R*phi plane
 927   // if absolute rotation alpha is in global system
 928   // otherwise alpha rotation is relative to the current rotation angle
 929   //
 930
 931   if (absolute) {
 932     alpha -= GetAlpha();
 933   }
 934   else{
 935     fNRotate++;
 936   }
 937
 938   return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
 939
 940 }
 941
 942 //_____________________________________________________________________________
 943 Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
 944 {
 945   //
 946   // Returns the track chi2
 947   //
 948
 949   Double_t p[2]   = { c->GetY()
 950                     , c->GetZ() };
 951   Double_t sy2    = c->GetSigmaY2() * 4.0;
 952   Double_t sz2    = c->GetSigmaZ2() * 4.0;
 953   Double_t cov[3] = { sy2 + h01*h01*sz2
 954                     , h01*(sy2-sz2)
 955                     , sz2 + h01*h01*sy2 };
 956
 957   return AliExternalTrackParam::GetPredictedChi2(p,cov);
 958
 959 }
 960
 961 //_____________________________________________________________________________
 962 void AliTRDtrack::MakeBackupTrack()
 963 {
 964   //
 965   // Creates a backup track
 966   //
 967
 968   if (fBackupTrack) {
 969     delete fBackupTrack;
 970   }
 971   fBackupTrack = new AliTRDtrack(*this);
 972
 973 }
 974
 975 //_____________________________________________________________________________
 976 Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
 977 {
 978   //
 979   // Find a prolongation at given x
 980   // Return 0 if it does not exist
 981   //
 982
 983   Double_t bz = GetBz();
 984
 985   if (!AliExternalTrackParam::GetYAt(xk,bz,y)) {
 986     return 0;
 987   }
 988   if (!AliExternalTrackParam::GetZAt(xk,bz,z)) {
 989     return 0;
 990   }
 991
 992   return 1;
 993
 994 }
 995
 996 //_____________________________________________________________________________
 997 Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
 998 {
 999   //
1000   // Propagate track to given x  position
1001   // Works inside of the 20 degree segmentation
1002   // (local cooordinate frame for TRD , TPC, TOF)
1003   //
1004   // Material budget from geo manager
1005   //
1006
1007   Double_t xyz0[3];
1008   Double_t xyz1[3];
1009   Double_t y;
1010   Double_t z;
1011
1012   const Double_t kAlphac  = TMath::Pi()/9.0;
1013   const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
1014
1015   // Critical alpha  - cross sector indication
1016   Double_t dir = (GetX() > xr) ? -1.0 : 1.0;
1017
1018   // Direction +-
1019   for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) {
1020
1021     GetXYZ(xyz0);
1022     GetProlongation(x,y,z);
1023     xyz1[0] = x * TMath::Cos(GetAlpha()) + y * TMath::Sin(GetAlpha());
1024     xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha());
1025     xyz1[2] = z;
1026     Double_t param[7];
1027     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1028
1029     if ((param[0] > 0) &&
1030         (param[1] > 0)) {
1031       PropagateTo(x,param[1],param[0]*param[4]);
1032     }
1033
1034     if (GetY() >  GetX()*kTalphac) {
1035       Rotate(-kAlphac);
1036     }
1037     if (GetY() < -GetX()*kTalphac) {
1038       Rotate( kAlphac);
1039     }
1040
1041   }
1042
1043   PropagateTo(xr);
1044
1045   return 0;
1046
1047 }
1048
1049 //_____________________________________________________________________________
1050 Int_t   AliTRDtrack::PropagateToR(Double_t r,Double_t step)
1051 {
1052   //
1053   // Propagate track to the radial position
1054   // Rotation always connected to the last track position
1055   //
1056
1057   Double_t xyz0[3];
1058   Double_t xyz1[3];
1059   Double_t y;
1060   Double_t z;
1061
1062   Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
1063   // Direction +-
1064   Double_t dir    = (radius > r) ? -1.0 : 1.0;
1065
1066   for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
1067
1068     GetXYZ(xyz0);
1069     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1070     Rotate(alpha,kTRUE);
1071     GetXYZ(xyz0);
1072     GetProlongation(x,y,z);
1073     xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1074     xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1075     xyz1[2] = z;
1076     Double_t param[7];
1077     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1078     if (param[1] <= 0) {
1079       param[1] = 100000000;
1080     }
1081     PropagateTo(x,param[1],param[0]*param[4]);
1082
1083   }
1084
1085   GetXYZ(xyz0);
1086   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1087   Rotate(alpha,kTRUE);
1088   GetXYZ(xyz0);
1089   GetProlongation(r,y,z);
1090   xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1091   xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1092   xyz1[2] = z;
1093   Double_t param[7];
1094   AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
1095
1096   if (param[1] <= 0) {
1097     param[1] = 100000000;
1098   }
1099   PropagateTo(r,param[1],param[0]*param[4]);
1100
1101   return 0;
1102
1103 }
1104
1105 //_____________________________________________________________________________
1106 Int_t AliTRDtrack::GetSector() const
1107 {
1108   //
1109   // Return the current sector
1110   //
1111
1112   return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha())
1113              % AliTRDgeometry::kNsect;
1114
1115 }
1116
1117 //_____________________________________________________________________________
1118 void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
1119 {
1120   //
1121   // The sampled energy loss
1122   //
1123
1124   Double_t s = GetSnp();
1125   Double_t t = GetTgl();
1126   q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1127   fdQdl[i] = q;
1128
1129 }
1130
1131 //_____________________________________________________________________________
1132 void AliTRDtrack::SetSampledEdx(Float_t q)
1133 {
1134   //
1135   // The sampled energy loss
1136   //
1137
1138   Double_t s = GetSnp();
1139   Double_t t = GetTgl();
1140   q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1141   fdQdl[fNdedx] = q;
1142   fNdedx++;
1143
1144 }
1145
1146 //_____________________________________________________________________________
1147 Double_t AliTRDtrack::GetBz() const
1148 {
1149   //
1150   // Returns Bz component of the magnetic field (kG)
1151   //
1152
1153   if (AliTracker::UniformField()) {
1154     return AliTracker::GetBz();
1155   }
1156   Double_t r[3];
1157   GetXYZ(r);
1158
1159   return AliTracker::GetBz(r);
1160
1161 }