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