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 <Riostream.h>
  19
  20 #include <TMath.h>
  21 #include <TVector2.h>
  22
  23 #include "AliTracker.h"
  24 #include "AliESDtrack.h"
  25
  26 #include "AliTRDgeometry.h"
  27 #include "AliTRDcluster.h"
  28 #include "AliTRDtrack.h"
  29 #include "AliTRDtracklet.h"
  30
  31 ClassImp(AliTRDtrack)
  32
  33 ///////////////////////////////////////////////////////////////////////////////
  34 //                                                                           //
  35 //  Represents a reconstructed TRD track                                     //
  36 //  Local TRD Kalman track                                                   //
  37 //                                                                           //
  38 ///////////////////////////////////////////////////////////////////////////////
  39
  40 //_____________________________________________________________________________
  41 AliTRDtrack::AliTRDtrack()
  42   :AliKalmanTrack()
  43   ,fSeedLab(-1)
  44   ,fdEdx(0)
  45   ,fDE(0)
  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   }
  68
  69   for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
  70     fIndex[i]       = 0;
  71     fIndexBackup[i] = 0;
  72     fdQdl[i]        = 0;
  73   }
  74
  75   for (Int_t i = 0; i < 3; i++) {
  76     fBudget[i] = 0;
  77   }
  78
  79 }
  80
  81 //_____________________________________________________________________________
  82 AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, Int_t index
  83                        , const Double_t p[5], const Double_t cov[15]
  84                        , Double_t x, Double_t alpha)
  85   :AliKalmanTrack()
  86   ,fSeedLab(-1)
  87   ,fdEdx(0)
  88   ,fDE(0)
  89   ,fStopped(kFALSE)
  90   ,fLhElectron(0)
  91   ,fNWrong(0)
  92   ,fNRotate(0)
  93   ,fNCross(0)
  94   ,fNExpected(0)
  95   ,fNLast(0)
  96   ,fNExpectedLast(0)
  97   ,fNdedx(0)
  98   ,fChi2Last(1e10)
  99   ,fBackupTrack(0x0)
 100 {
 101   //
 102   // The main AliTRDtrack constructor.
 103   //
 104
 105   Double_t cnv   = 1.0/(GetBz() * kB2C);
 106
 107   Double_t pp[5] = { p[0]
 108                    , p[1]
 109                    , x*p[4] - p[2]
 110                    , p[3]
 111                    , p[4]*cnv      };
 112
 113   Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[5];
 114   Double_t c32 = x*cov[13] - cov[8];
 115   Double_t c20 = x*cov[10] - cov[3];
 116   Double_t c21 = x*cov[11] - cov[4];
 117   Double_t c42 = x*cov[14] - cov[12];
 118
 119   Double_t cc[15] = { cov[0 ]
 120                     , cov[1 ],     cov[2 ]
 121                     , c20,         c21,         c22
 122                     , cov[6 ],     cov[7 ],     c32,     cov[9 ]
 123                     , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
 124
 125   Set(x,alpha,pp,cc);
 126   SetNumberOfClusters(1);
 127   fIndex[0] = index;
 128
 129   for (Int_t i = 0; i < kNplane; i++) {
 130     for (Int_t j = 0; j < kNslice; j++) {
 131       fdEdxPlane[i][j] = 0.0;
 132     }
 133     fTimBinPlane[i] = -1;
 134   }
 135
 136   Double_t q = TMath::Abs(c->GetQ());
 137   Double_t s = GetSnp();
 138   Double_t t = GetTgl();
 139   if (s*s < 1) {
 140     q *= TMath::Sqrt((1-s*s)/(1+t*t));
 141   }
 142
 143   fdQdl[0] = q;
 144
 145   for (UInt_t i = 1; i < kMAXCLUSTERSPERTRACK; i++) {
 146     fdQdl[i]        = 0;
 147     fIndex[i]       = 0;
 148     fIndexBackup[i] = 0;
 149   }
 150
 151   for (Int_t i = 0; i < 3;i++) {
 152     fBudget[i]      = 0;
 153   }
 154
 155 }
 156
 157 //_____________________________________________________________________________
 158 AliTRDtrack::AliTRDtrack(const AliTRDtrack &t)
 159   :AliKalmanTrack(t)
 160   ,fSeedLab(t.GetSeedLabel())
 161   ,fdEdx(t.fdEdx)
 162   ,fDE(t.fDE)
 163   ,fStopped(t.fStopped)
 164   ,fLhElectron(0)
 165   ,fNWrong(t.fNWrong)
 166   ,fNRotate(t.fNRotate)
 167   ,fNCross(t.fNCross)
 168   ,fNExpected(t.fNExpected)
 169   ,fNLast(t.fNLast)
 170   ,fNExpectedLast(t.fNExpectedLast)
 171   ,fNdedx(t.fNdedx)
 172   ,fChi2Last(t.fChi2Last)
 173   ,fBackupTrack(0x0)
 174 {
 175   //
 176   // Copy constructor.
 177   //
 178
 179   for (Int_t i = 0; i < kNplane; i++) {
 180     for (Int_t j = 0; j < kNslice; j++) {
 181       fdEdxPlane[i][j] = t.fdEdxPlane[i][j];
 182     }
 183     fTimBinPlane[i] = t.fTimBinPlane[i];
 184     fTracklets[i]   = t.fTracklets[i];
 185   }
 186
 187   Int_t n = t.GetNumberOfClusters();
 188   SetNumberOfClusters(n);
 189
 190   for (Int_t i = 0; i < n; i++) {
 191     fIndex[i]       = t.fIndex[i];
 192     fIndexBackup[i] = t.fIndex[i];
 193     fdQdl[i]        = t.fdQdl[i];
 194   }
 195
 196   for (UInt_t i = n; i < kMAXCLUSTERSPERTRACK; i++) {
 197     fdQdl[i]        = 0;
 198     fIndex[i]       = 0;
 199     fIndexBackup[i] = 0;
 200   }
 201
 202   for (Int_t i = 0; i < 3;i++) {
 203     fBudget[i]      = t.fBudget[i];
 204   }
 205
 206 }
 207
 208 //_____________________________________________________________________________
 209 AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t /*alpha*/)
 210   :AliKalmanTrack(t)
 211   ,fSeedLab(-1)
 212   ,fdEdx(t.GetPIDsignal())
 213   ,fDE(0)
 214   ,fStopped(kFALSE)
 215   ,fLhElectron(0.0)
 216   ,fNWrong(0)
 217   ,fNRotate(0)
 218   ,fNCross(0)
 219   ,fNExpected(0)
 220   ,fNLast(0)
 221   ,fNExpectedLast(0)
 222   ,fNdedx(0)
 223   ,fChi2Last(0.0)
 224   ,fBackupTrack(0x0)
 225 {
 226   //
 227   // Constructor from AliTPCtrack or AliITStrack
 228   //
 229
 230   SetLabel(t.GetLabel());
 231   SetChi2(0.0);
 232   SetMass(t.GetMass());
 233   SetNumberOfClusters(0);
 234
 235   for (Int_t i = 0; i < kNplane; i++) {
 236     for (Int_t j = 0; j < kNslice; j++) {
 237       fdEdxPlane[i][j] = 0.0;
 238     }
 239     fTimBinPlane[i] = -1;
 240   }
 241
 242   for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
 243     fdQdl[i]        = 0;
 244     fIndex[i]       = 0;
 245     fIndexBackup[i] = 0;
 246   }
 247
 248   for (Int_t i = 0; i < 3; i++) {
 249     fBudget[i]      = 0;
 250   }
 251
 252 }
 253
 254 //_____________________________________________________________________________
 255 AliTRDtrack::AliTRDtrack(const AliESDtrack &t)
 256   :AliKalmanTrack()
 257   ,fSeedLab(-1)
 258   ,fdEdx(t.GetTRDsignal())
 259   ,fDE(0)
 260   ,fStopped(kFALSE)
 261   ,fLhElectron(0)
 262   ,fNWrong(0)
 263   ,fNRotate(0)
 264   ,fNCross(0)
 265   ,fNExpected(0)
 266   ,fNLast(0)
 267   ,fNExpectedLast(0)
 268   ,fNdedx(0)
 269   ,fChi2Last(1e10)
 270   ,fBackupTrack(0x0)
 271 {
 272   //
 273   // Constructor from AliESDtrack
 274   //
 275
 276   SetLabel(t.GetLabel());
 277   SetChi2(0.0);
 278   SetMass(t.GetMass());
 279   SetNumberOfClusters(t.GetTRDclusters(fIndex));
 280
 281   Int_t ncl = t.GetTRDclusters(fIndexBackup);
 282   for (UInt_t i = ncl; i < kMAXCLUSTERSPERTRACK; i++) {
 283     fIndexBackup[i] = 0;
 284     fIndex[i]       = 0;
 285   }
 286
 287   for (Int_t i = 0; i < kNplane; i++) {
 288     for (Int_t j = 0; j < kNslice; j++) {
 289       fdEdxPlane[i][j] = t.GetTRDsignals(i,j);
 290     }
 291     fTimBinPlane[i] = t.GetTRDTimBin(i);
 292   }
 293
 294   const AliExternalTrackParam *par = &t;
 295   if (t.GetStatus() & AliESDtrack::kTRDbackup) {
 296     par = t.GetOuterParam();
 297     if (!par) {
 298       AliError("No backup info!");
 299       par = &t;
 300     }
 301   }
 302   Set(par->GetX(),par->GetAlpha(),par->GetParameter(),par->GetCovariance());
 303
 304
 305   for (UInt_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
 306     fdQdl[i]   = 0;
 307   }
 308
 309   for (Int_t i = 0; i < 3; i++) {
 310     fBudget[i] = 0;
 311   }
 312
 313   if ((t.GetStatus() & AliESDtrack::kTIME) == 0) {
 314     return;
 315   }
 316
 317   StartTimeIntegral();
 318   Double_t times[10];
 319   t.GetIntegratedTimes(times);
 320   SetIntegratedTimes(times);
 321   SetIntegratedLength(t.GetIntegratedLength());
 322
 323 }
 324
 325 //____________________________________________________________________________
 326 AliTRDtrack::~AliTRDtrack()
 327 {
 328   //
 329   // Destructor
 330   //
 331
 332   if (fBackupTrack) {
 333     delete fBackupTrack;
 334   }
 335   fBackupTrack = 0;
 336
 337 }
 338
 339 //____________________________________________________________________________
 340 Float_t AliTRDtrack::StatusForTOF()
 341 {
 342   //
 343   // Defines the status of the TOF extrapolation
 344   //
 345
 346   // Definition of res ????
 347   Float_t res = (0.2 + 0.8 * (fN / (fNExpected + 5.0)))
 348                      * (0.4 + 0.6 * fTracklets[5].GetN() / 20.0);
 349   res *= (0.25 + 0.8 * 40.0 / (40.0 + fBudget[2]));
 350   return res;
 351
 352   // This part of the function is never reached ????
 353   // What defines these parameters ????
 354   Int_t status = 0;
 355   if (GetNumberOfClusters()                <  20)  return 0;
 356   if ((fN                                  > 110) &&
 357       (fChi2/(Float_t(fN))                 <   3)) return 3; // Gold
 358   if ((fNLast                              >  30) &&
 359       (fChi2Last/(Float_t(fNLast))         <   3)) return 3; // Gold
 360   if ((fNLast                              >  20) &&
 361       (fChi2Last/(Float_t(fNLast))         <   2)) return 3; // Gold
 362   if ((fNLast/(fNExpectedLast+3.0)         > 0.8) &&
 363       (fChi2Last/Float_t(fNLast)           <   5) &&
 364       (fNLast                              >  20)) return 2; // Silber
 365   if ((fNLast                              >   5) &&
 366       (((fNLast+1.0)/(fNExpectedLast+1.0)) > 0.8) &&
 367       (fChi2Last/(fNLast-5.0)              <   6)) return 1;
 368
 369   return status;
 370
 371 }
 372
 373 //_____________________________________________________________________________
 374 Int_t AliTRDtrack::Compare(const TObject *o) const
 375 {
 376   //
 377   // Compares tracks according to their Y2 or curvature
 378   //
 379
 380   AliTRDtrack *t = (AliTRDtrack *) o;
 381
 382   Double_t co = TMath::Abs(t->GetC());
 383   Double_t c  = TMath::Abs(GetC());
 384
 385   if      (c > co) {
 386     return 1;
 387   }
 388   else if (c < co) {
 389     return -1;
 390   }
 391   return 0;
 392
 393 }
 394
 395 //_____________________________________________________________________________
 396 void AliTRDtrack::CookdEdx(Double_t low, Double_t up)
 397 {
 398   //
 399   // Calculates the truncated dE/dx within the "low" and "up" cuts.
 400   //
 401
 402   Int_t   i  = 0;
 403
 404   // Array to sort the dEdx values according to amplitude
 405   Float_t sorted[kMAXCLUSTERSPERTRACK];
 406
 407   // Number of clusters used for dedx
 408   Int_t   nc = fNdedx;
 409
 410   // Require at least 10 clusters for a dedx measurement
 411   if (nc < 10) {
 412     SetdEdx(0);
 413     return;
 414   }
 415
 416   // Lower and upper bound
 417   Int_t nl = Int_t(low * nc);
 418   Int_t nu = Int_t( up * nc);
 419
 420   // Can fdQdl be negative ????
 421   for (i = 0; i < nc; i++) {
 422     sorted[i] = TMath::Abs(fdQdl[i]);
 423   }
 424
 425   // Sort the dedx values by amplitude
 426   Int_t *index = new Int_t[nc];
 427   TMath::Sort(nc,sorted,index,kFALSE);
 428
 429   // Sum up the truncated charge between nl and nu
 430   Float_t dedx = 0.0;
 431   for (i = nl; i <= nu; i++) {
 432     dedx += sorted[index[i]];
 433   }
 434   dedx /= (nu - nl + 1.0);
 435   SetdEdx(dedx);
 436
 437 }
 438
 439 //_____________________________________________________________________________
 440 Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t x0, Double_t rho)
 441 {
 442   //
 443   // Propagates a track of particle with mass=pm to a reference plane
 444   // defined by x=xk through media of density=rho and radiationLength=x0
 445   //
 446
 447   if (xk == GetX()) {
 448     return kTRUE;
 449   }
 450
 451   Double_t oldX = GetX();
 452   Double_t oldY = GetY();
 453   Double_t oldZ = GetZ();
 454
 455   Double_t bz   = GetBz();
 456
 457   if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
 458     return kFALSE;
 459   }
 460
 461   Double_t x = GetX();
 462   Double_t y = GetY();
 463   Double_t z = GetZ();
 464   Double_t d = TMath::Sqrt((x-oldX)*(x-oldX)
 465                          + (y-oldY)*(y-oldY)
 466                          + (z-oldZ)*(z-oldZ));
 467   if (oldX < xk) {
 468     if (IsStartedTimeIntegral()) {
 469       Double_t l2  = d;
 470       Double_t crv = GetC();
 471       if (TMath::Abs(l2*crv) > 0.0001) {
 472         // Make correction for curvature if neccesary
 473         l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY));
 474         l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
 475         l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ));
 476       }
 477       AddTimeStep(l2);
 478     }
 479   }
 480
 481   Double_t ll = (oldX < xk) ? -d : d;
 482   if (!AliExternalTrackParam::CorrectForMaterial(ll*rho/x0,x0,GetMass())) {
 483     return kFALSE;
 484   }
 485
 486   {
 487
 488     // Energy losses************************
 489     Double_t p2    = (1.0 + GetTgl()*GetTgl()) / (Get1Pt()*Get1Pt());
 490     Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
 491     if (beta2<1.0e-10 || (5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0) {
 492       return kFALSE;
 493     }
 494
 495     Double_t dE    = 0.153e-3 / beta2
 496                    * (log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
 497                    * d * rho;
 498     Float_t budget = d * rho;
 499     fBudget[0] += budget;
 500
 501     /*
 502     // Suspicious part - think about it ?
 503     Double_t kinE =  TMath::Sqrt(p2);
 504     if (dE > 0.8*kinE) dE = 0.8 * kinE;  //
 505     if (dE < 0)        dE = 0.0;         // Not valid region for Bethe bloch
 506     */
 507
 508     fDE += dE;
 509
 510     /*
 511     // Suspicious ! I.B.
 512     Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE));   // Energy loss fluctuation
 513     Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
 514     fCcc += sigmac2;
 515     fCee += fX*fX * sigmac2;
 516     */
 517
 518   }
 519
 520   return kTRUE;
 521
 522 }
 523
 524 //_____________________________________________________________________________
 525 Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index
 526                          , Double_t h01)
 527 {
 528   //
 529   // Assignes found cluster to the track and updates track information
 530   //
 531
 532   Bool_t fNoTilt = kTRUE;
 533   if (TMath::Abs(h01) > 0.003) {
 534     fNoTilt = kFALSE;
 535   }
 536
 537   // Add angular effect to the error contribution -  MI
 538   Float_t tangent2 = GetSnp()*GetSnp();
 539   if (tangent2 < 0.90000) {
 540     tangent2 = tangent2 / (1.0 - tangent2);
 541   }
 542   //Float_t errang = tangent2 * 0.04;
 543
 544   Double_t p[2]   = {c->GetY(), c->GetZ() };
 545   //Double_t cov[3] = {c->GetSigmaY2()+errang, 0.0, c->GetSigmaZ2()*100.0 };
 546   Double_t sy2    = c->GetSigmaY2() * 4.0;
 547   Double_t sz2    = c->GetSigmaZ2() * 4.0;
 548   Double_t cov[3] = {sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
 549
 550   if (!AliExternalTrackParam::Update(p,cov)) {
 551     return kFALSE;
 552   }
 553
 554   Int_t n = GetNumberOfClusters();
 555   fIndex[n] = index;
 556   SetNumberOfClusters(n+1);
 557
 558   SetChi2(GetChi2()+chisq);
 559
 560   return kTRUE;
 561
 562 }
 563
 564 //_____________________________________________________________________________
 565 Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, Int_t index
 566                           , Double_t h01, Int_t /*plane*/)
 567 {
 568   //
 569   // Assignes found cluster to the track and updates track information
 570   //
 571
 572   Bool_t fNoTilt = kTRUE;
 573   if (TMath::Abs(h01) > 0.003) {
 574     fNoTilt = kFALSE;
 575   }
 576
 577   // Add angular effect to the error contribution and make correction  -  MI
 578   Double_t tangent2 = GetSnp()*GetSnp();
 579   if (tangent2 < 0.90000){
 580     tangent2 = tangent2 / (1.0-tangent2);
 581   }
 582   Double_t tangent = TMath::Sqrt(tangent2);
 583   if (GetSnp() < 0) {
 584     tangent *= -1;
 585   }
 586
 587   //
 588   // Is the following still needed ????
 589   //
 590
 591   //  Double_t correction = 0*plane;
 592   /*
 593   Double_t errang = tangent2*0.04;  //
 594   Double_t errsys =0.025*0.025*20;  //systematic error part
 595
 596   Float_t extend =1;
 597   if (c->GetNPads()==4) extend=2;
 598   */
 599   //if (c->GetNPads()==5)  extend=3;
 600   //if (c->GetNPads()==6)  extend=3;
 601   //if (c->GetQ()<15) return 1;
 602
 603   /*
 604   if (corrector!=0){
 605   //if (0){
 606     correction = corrector->GetCorrection(plane,c->GetLocalTimeBin(),tangent);
 607     if (TMath::Abs(correction)>0){
 608       //if we have info
 609       errang     = corrector->GetSigma(plane,c->GetLocalTimeBin(),tangent);
 610       errang    *= errang;
 611       errang    += tangent2*0.04;
 612     }
 613   }
 614   */
 615   //
 616   //Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
 617   /*
 618     {
 619       Double_t dy=c->GetY() - GetY(), dz=c->GetZ() - GetZ();
 620       printf("%e %e %e %e\n",dy,dz,padlength/2,h01);
 621     }
 622   */
 623
 624   Double_t p[2]   = { c->GetY(), c->GetZ() };
 625   //Double_t cov[3]={ (c->GetSigmaY2()+errang+errsys)*extend, 0.0, c->GetSigmaZ2()*10000.0 };
 626   Double_t sy2    = c->GetSigmaY2() * 4.0;
 627   Double_t sz2    = c->GetSigmaZ2() * 4.0;
 628   Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
 629
 630   if (!AliExternalTrackParam::Update(p,cov)) {
 631     return kFALSE;
 632   }
 633
 634   Int_t n = GetNumberOfClusters();
 635   fIndex[n] = index;
 636   SetNumberOfClusters(n+1);
 637   SetChi2(GetChi2() + chisq);
 638
 639   return kTRUE;
 640
 641 }
 642
 643 // //_____________________________________________________________________________
 644 // Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet)
 645 // {
 646 //   //
 647 //   // Assignes found tracklet to the track and updates track information
 648 //   //
 649 //   // Can this be removed ????
 650 //   //
 651 //
 652 //   Double_t r00=(tracklet.GetTrackletSigma2()), r01=0., r11= 10000.;
 653 //   r00+=fCyy; r01+=fCzy; r11+=fCzz;
 654 //   //
 655 //   Double_t det=r00*r11 - r01*r01;
 656 //   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 657 //   //
 658
 659 //   Double_t dy=tracklet.GetY() - fY, dz=tracklet.GetZ() - fZ;
 660
 661
 662 //   Double_t s00 = tracklet.GetTrackletSigma2();  // error pad
 663 //   Double_t s11 = 100000;   // error pad-row
 664 //   Float_t  h01 = tracklet.GetTilt();
 665 //   //
 666 //   //  r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00;
 667 //   r00 = fCyy + fCzz*h01*h01+s00;
 668 //   //  r01 = fCzy + fCzz*h01;
 669 //   r01 = fCzy ;
 670 //   r11 = fCzz + s11;
 671 //   det = r00*r11 - r01*r01;
 672 //   // inverse matrix
 673 //   tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 674
 675 //   Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11;
 676 //   Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11;
 677 //   Double_t k20=fCey*r00+fCez*r01, k21=fCey*r01+fCez*r11;
 678 //   Double_t k30=fCty*r00+fCtz*r01, k31=fCty*r01+fCtz*r11;
 679 //   Double_t k40=fCcy*r00+fCcz*r01, k41=fCcy*r01+fCcz*r11;
 680
 681 //   // K matrix
 682 // //   k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01);
 683 // //   k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01);
 684 // //   k20=fCey*r00+fCez*(r01+h01*r00),k21=fCey*r01+fCez*(r11+h01*r01);
 685 // //   k30=fCty*r00+fCtz*(r01+h01*r00),k31=fCty*r01+fCtz*(r11+h01*r01);
 686 // //   k40=fCcy*r00+fCcz*(r01+h01*r00),k41=fCcy*r01+fCcz*(r11+h01*r01);
 687 //   //
 688 //   //Update measurement
 689 //   Double_t cur=fC + k40*dy + k41*dz, eta=fE + k20*dy + k21*dz;
 690 //   //  cur=fC + k40*dy + k41*dz; eta=fE + k20*dy + k21*dz;
 691 //   if (TMath::Abs(cur*fX-eta) >= 0.90000) {
 692 //     //Int_t n=GetNumberOfClusters();
 693 //     //      if (n>4) cerr<<n<<" AliTRDtrack warning: Filtering failed !\n";
 694 //     return 0;
 695 //   }
 696 // //   k01+=h01*k00;
 697 // //   k11+=h01*k10;
 698 // //   k21+=h01*k20;
 699 // //   k31+=h01*k30;
 700 // //   k41+=h01*k40;
 701
 702
 703 //   fY += k00*dy + k01*dz;
 704 //   fZ += k10*dy + k11*dz;
 705 //   fE  = eta;
 706 //   fT += k30*dy + k31*dz;
 707 //   fC  = cur;
 708
 709
 710 //   //Update covariance
 711 //   //
 712 //   //
 713 //   Double_t oldyy = fCyy, oldzz = fCzz; //, oldee=fCee, oldcc =fCcc;
 714 //   Double_t oldzy = fCzy, oldey = fCey, oldty=fCty, oldcy =fCcy;
 715 //   Double_t oldez = fCez, oldtz = fCtz, oldcz=fCcz;
 716 //   //Double_t oldte = fCte, oldce = fCce;
 717 //   //Double_t oldct = fCct;
 718
 719 //   fCyy-=k00*oldyy+k01*oldzy;
 720 //   fCzy-=k10*oldyy+k11*oldzy;
 721 //   fCey-=k20*oldyy+k21*oldzy;
 722 //   fCty-=k30*oldyy+k31*oldzy;
 723 //   fCcy-=k40*oldyy+k41*oldzy;
 724 //   //
 725 //   fCzz-=k10*oldzy+k11*oldzz;
 726 //   fCez-=k20*oldzy+k21*oldzz;
 727 //   fCtz-=k30*oldzy+k31*oldzz;
 728 //   fCcz-=k40*oldzy+k41*oldzz;
 729 //   //
 730 //   fCee-=k20*oldey+k21*oldez;
 731 //   fCte-=k30*oldey+k31*oldez;
 732 //   fCce-=k40*oldey+k41*oldez;
 733 //   //
 734 //   fCtt-=k30*oldty+k31*oldtz;
 735 //   fCct-=k40*oldty+k41*oldtz;
 736 //   //
 737 //   fCcc-=k40*oldcy+k41*oldcz;
 738 //   //
 739
 740 //   //Int_t n=GetNumberOfClusters();
 741 //   //fIndex[n]=index;
 742 //   //SetNumberOfClusters(n+1);
 743
 744 //   //SetChi2(GetChi2()+chisq);
 745 //   //  cerr<<"in update: fIndex["<<fN<<"] = "<<index<<endl;
 746
 747 //   return 1;
 748
 749 // }
 750
 751 //_____________________________________________________________________________
 752 Bool_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
 753 {
 754   //
 755   // Rotates track parameters in R*phi plane
 756   // if absolute rotation alpha is in global system
 757   // otherwise alpha rotation is relative to the current rotation angle
 758   //
 759
 760   if (absolute) {
 761     alpha -= GetAlpha();
 762   }
 763   else{
 764     fNRotate++;
 765   }
 766
 767   return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
 768
 769 }
 770
 771 //_____________________________________________________________________________
 772 Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
 773 {
 774   //
 775   // Returns the track chi2
 776   //
 777
 778   Double_t p[2]   = { c->GetY(), c->GetZ() };
 779   Double_t sy2    = c->GetSigmaY2() * 4.0;
 780   Double_t sz2    = c->GetSigmaZ2() * 4.0;
 781   Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
 782
 783   return AliExternalTrackParam::GetPredictedChi2(p,cov);
 784
 785   //
 786   // Can the following be removed ????
 787   //
 788   /*
 789   Bool_t fNoTilt = kTRUE;
 790   if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE;
 791
 792   return (c->GetY() - GetY())*(c->GetY() - GetY())/c->GetSigmaY2();
 793   */
 794
 795   /*
 796   Double_t chi2, dy, r00, r01, r11;
 797
 798   if(fNoTilt) {
 799     dy=c->GetY() - fY;
 800     r00=c->GetSigmaY2();
 801     chi2 = (dy*dy)/r00;
 802   }
 803   else {
 804     Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12);
 805     //
 806     r00=c->GetSigmaY2(); r01=0.; r11=c->GetSigmaZ2();
 807     r00+=fCyy; r01+=fCzy; r11+=fCzz;
 808
 809     Double_t det=r00*r11 - r01*r01;
 810     if (TMath::Abs(det) < 1.e-10) {
 811       Int_t n=GetNumberOfClusters();
 812       if (n>4) cerr<<n<<" AliTRDtrack warning: Singular matrix !\n";
 813       return 1e10;
 814     }
 815     Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
 816     Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
 817     Double_t tiltdz = dz;
 818     if (TMath::Abs(tiltdz)>padlength/2.) {
 819       tiltdz = TMath::Sign(padlength/2,dz);
 820     }
 821     //    dy=dy+h01*dz;
 822     dy=dy+h01*tiltdz;
 823
 824     chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
 825   }
 826
 827   return chi2;
 828   */
 829
 830 }
 831
 832 //_____________________________________________________________________________
 833 void AliTRDtrack::MakeBackupTrack()
 834 {
 835   //
 836   // Creates a backup track
 837   //
 838
 839   if (fBackupTrack) {
 840     delete fBackupTrack;
 841   }
 842   fBackupTrack = new AliTRDtrack(*this);
 843
 844 }
 845
 846 //_____________________________________________________________________________
 847 Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
 848 {
 849   //
 850   // Find a prolongation at given x
 851   // Return 0 if it does not exist
 852   //
 853
 854   Double_t bz = GetBz();
 855
 856   if (!AliExternalTrackParam::GetYAt(xk,bz,y)) {
 857     return 0;
 858   }
 859   if (!AliExternalTrackParam::GetZAt(xk,bz,z)) {
 860     return 0;
 861   }
 862
 863   return 1;
 864
 865 }
 866
 867 //_____________________________________________________________________________
 868 Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
 869 {
 870   //
 871   // Propagate track to given x  position
 872   // Works inside of the 20 degree segmentation (local cooordinate frame for TRD , TPC, TOF)
 873   //
 874   // Material budget from geo manager
 875   //
 876
 877   Double_t xyz0[3];
 878   Double_t xyz1[3];
 879   Double_t y;
 880   Double_t z;
 881
 882   const Double_t kAlphac  = TMath::Pi()/9.0;
 883   const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
 884
 885   // Critical alpha  - cross sector indication
 886   Double_t dir = (GetX()>xr) ? -1.0 : 1.0;
 887
 888   // Direction +-
 889   for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) {
 890
 891     GetXYZ(xyz0);
 892     GetProlongation(x,y,z);
 893     xyz1[0] = x * TMath::Cos(GetAlpha()) + y * TMath::Sin(GetAlpha());
 894     xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha());
 895     xyz1[2] = z;
 896     Double_t param[7];
 897     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
 898
 899     if ((param[0] > 0) &&
 900         (param[1] > 0)) {
 901       PropagateTo(x,param[1],param[0]);
 902     }
 903
 904     if (GetY() >  GetX()*kTalphac) {
 905       Rotate(-kAlphac);
 906     }
 907     if (GetY() < -GetX()*kTalphac) {
 908       Rotate( kAlphac);
 909     }
 910
 911   }
 912
 913   PropagateTo(xr);
 914
 915   return 0;
 916
 917 }
 918
 919 //_____________________________________________________________________________
 920 Int_t   AliTRDtrack::PropagateToR(Double_t r,Double_t step)
 921 {
 922   //
 923   // Propagate track to the radial position
 924   // Rotation always connected to the last track position
 925   //
 926
 927   Double_t xyz0[3];
 928   Double_t xyz1[3];
 929   Double_t y;
 930   Double_t z;
 931
 932   Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
 933   // Direction +-
 934   Double_t dir    = (radius>r) ? -1.0 : 1.0;
 935
 936   for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
 937
 938     GetXYZ(xyz0);
 939     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 940     Rotate(alpha,kTRUE);
 941     GetXYZ(xyz0);
 942     GetProlongation(x,y,z);
 943     xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
 944     xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
 945     xyz1[2] = z;
 946     Double_t param[7];
 947     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
 948     if (param[1] <= 0) {
 949       param[1] = 100000000;
 950     }
 951     PropagateTo(x,param[1],param[0]);
 952
 953   }
 954
 955   GetXYZ(xyz0);
 956   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 957   Rotate(alpha,kTRUE);
 958   GetXYZ(xyz0);
 959   GetProlongation(r,y,z);
 960   xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
 961   xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
 962   xyz1[2] = z;
 963   Double_t param[7];
 964   AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
 965
 966   if (param[1] <= 0) {
 967     param[1] = 100000000;
 968   }
 969   PropagateTo(r,param[1],param[0]);
 970
 971   return 0;
 972
 973 }
 974
 975 //_____________________________________________________________________________
 976 Int_t AliTRDtrack::GetSector() const
 977 {
 978   //
 979   // Return the current sector
 980   //
 981
 982   return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha())
 983          % AliTRDgeometry::kNsect;
 984
 985 }
 986
 987 //_____________________________________________________________________________
 988 void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
 989 {
 990   //
 991   // The sampled energy loss
 992   //
 993
 994   Double_t s = GetSnp();
 995   Double_t t = GetTgl();
 996   q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
 997   fdQdl[i] = q;
 998
 999 }
1000
1001 //_____________________________________________________________________________
1002 void AliTRDtrack::SetSampledEdx(Float_t q)
1003 {
1004   //
1005   // The sampled energy loss
1006   //
1007
1008   Double_t s = GetSnp();
1009   Double_t t = GetTgl();
1010   q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1011   fdQdl[fNdedx] = q;
1012   fNdedx++;
1013
1014 }
1015
1016 //_____________________________________________________________________________
1017 Double_t AliTRDtrack::GetBz() const
1018 {
1019   //
1020   // Returns Bz component of the magnetic field (kG)
1021   //
1022
1023   if (AliTracker::UniformField()) {
1024     return AliTracker::GetBz();
1025   }
1026   Double_t r[3];
1027   GetXYZ(r);
1028
1029   return AliTracker::GetBz(r);
1030
1031 }