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