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 #include <TMath.h>
  20 #include <TVector2.h>
  21
  22 #include "AliESDtrack.h"
  23 #include "AliTRDgeometry.h"
  24 #include "AliTRDcluster.h"
  25 #include "AliTRDtrack.h"
  26 #include "AliTRDtracklet.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   ,fAlpha(0)
  44   ,fX(0)
  45   ,fStopped(kFALSE)
  46   ,fY(0)
  47   ,fZ(0)
  48   ,fE(0)
  49   ,fT(0)
  50   ,fC(0)
  51   ,fCyy(1e10)
  52   ,fCzy(0)
  53   ,fCzz(1e10)
  54   ,fCey(0)
  55   ,fCez(0)
  56   ,fCee(1e10)
  57   ,fCty(0)
  58   ,fCtz(0)
  59   ,fCte(0)
  60   ,fCtt(1e10)
  61   ,fCcy(0)
  62   ,fCcz(0)
  63   ,fCce(0)
  64   ,fCct(0)
  65   ,fCcc(1e10)
  66   ,fLhElectron(0)
  67   ,fNWrong(0)
  68   ,fNRotate(0)
  69   ,fNCross(0)
  70   ,fNExpected(0)
  71   ,fNLast(0)
  72   ,fNExpectedLast(0)
  73   ,fNdedx(0)
  74   ,fChi2Last(1e10)
  75   ,fBackupTrack(0x0)
  76 {
  77   //
  78   // AliTRDtrack default constructor
  79   //
  80
  81   Int_t  i = 0;
  82   Int_t  j = 0;
  83   UInt_t k = 0;
  84
  85   for (i = 0; i < kNplane; i++) {
  86     for (j = 0; j < kNslice; j++) {
  87       fdEdxPlane[i][j] = 0;
  88     }
  89     fTimBinPlane[i] = -1;
  90   }
  91
  92   for (k = 0; k < kMAXCLUSTERSPERTRACK; k++) {
  93     fIndex[k]       = 0;
  94     fIndexBackup[k] = 0;
  95     fdQdl[k]        = 0;
  96   }
  97
  98   for (i = 0; i < 3; i++) {
  99     fBudget[i]      = 0;
 100   }
 101
 102 }
 103
 104 //_____________________________________________________________________________
 105 AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index,
 106                          const Double_t xx[5], const Double_t cc[15],
 107                          Double_t xref, Double_t alpha)
 108   :AliKalmanTrack()
 109   ,fSeedLab(-1)
 110   ,fdEdx(0.0)
 111   ,fDE(0.0)
 112   ,fAlpha(alpha)
 113   ,fX(xref)
 114   ,fStopped(kFALSE)
 115   ,fY(xx[0])
 116   ,fZ(xx[1])
 117   ,fE(xx[2])
 118   ,fT(xx[3])
 119   ,fC(xx[4])
 120   ,fCyy(cc[0])
 121   ,fCzy(cc[1])
 122   ,fCzz(cc[2])
 123   ,fCey(cc[3])
 124   ,fCez(cc[4])
 125   ,fCee(cc[5])
 126   ,fCty(cc[6])
 127   ,fCtz(cc[7])
 128   ,fCte(cc[8])
 129   ,fCtt(cc[9])
 130   ,fCcy(cc[10])
 131   ,fCcz(cc[11])
 132   ,fCce(cc[12])
 133   ,fCct(cc[13])
 134   ,fCcc(cc[14])
 135   ,fLhElectron(0.0)
 136   ,fNWrong(0)
 137   ,fNRotate(0)
 138   ,fNCross(0)
 139   ,fNExpected(0)
 140   ,fNLast(0)
 141   ,fNExpectedLast(0)
 142   ,fNdedx(0)
 143   ,fChi2Last(1e10)
 144   ,fBackupTrack(0x0)
 145 {
 146   //
 147   // AliTRDtrack main constructor
 148   //
 149
 150   Int_t  i = 0;
 151   Int_t  j = 0;
 152   UInt_t k = 0;
 153
 154   if (fAlpha <  -TMath::Pi()) {
 155     fAlpha += 2.0 * TMath::Pi();
 156   }
 157   if (fAlpha >=  TMath::Pi()) {
 158     fAlpha -= 2.0 * TMath::Pi();
 159   }
 160
 161   SaveLocalConvConst();
 162
 163   fIndex[0] = index;
 164   SetNumberOfClusters(1);
 165
 166   for (i = 0; i < kNplane; i++) {
 167     for (j = 0; j < kNslice; j++) {
 168       fdEdxPlane[i][j] = 0;
 169     }
 170     fTimBinPlane[i] = -1;
 171   }
 172
 173   Double_t q = TMath::Abs(c->GetQ());
 174   Double_t s = fX * fC - fE;
 175   Double_t t = fT;
 176   if (s*s < 1.0) {
 177     q *= TMath::Sqrt((1-s*s)/(1+t*t));
 178   }
 179   fdQdl[0] = q;
 180
 181   for (k = 1; k < kMAXCLUSTERSPERTRACK; k++) {
 182     fdQdl[k]        = 0;
 183     fIndex[k]       = 0;
 184     fIndexBackup[k] = 0;
 185   }
 186
 187   for (i = 0; i < 3; i++) {
 188     fBudget[i]      = 0;
 189   }
 190
 191 }
 192
 193 //_____________________________________________________________________________
 194 AliTRDtrack::AliTRDtrack(const AliTRDtrack &t)
 195   :AliKalmanTrack(t)
 196   ,fSeedLab(t.fSeedLab)
 197   ,fdEdx(t.fdEdx)
 198   ,fDE(t.fDE)
 199   ,fAlpha(t.fAlpha)
 200   ,fX(t.fX)
 201   ,fStopped(t.fStopped)
 202   ,fY(t.fY)
 203   ,fZ(t.fZ)
 204   ,fE(t.fE)
 205   ,fT(t.fT)
 206   ,fC(t.fC)
 207   ,fCyy(t.fCyy)
 208   ,fCzy(t.fCzy)
 209   ,fCzz(t.fCzz)
 210   ,fCey(t.fCey)
 211   ,fCez(t.fCez)
 212   ,fCee(t.fCee)
 213   ,fCty(t.fCty)
 214   ,fCtz(t.fCtz)
 215   ,fCte(t.fCte)
 216   ,fCtt(t.fCtt)
 217   ,fCcy(t.fCcy)
 218   ,fCcz(t.fCcz)
 219   ,fCce(t.fCce)
 220   ,fCct(t.fCct)
 221   ,fCcc(t.fCcc)
 222   ,fLhElectron(0.0)
 223   ,fNWrong(t.fNWrong)
 224   ,fNRotate(t.fNRotate)
 225   ,fNCross(t.fNCross)
 226   ,fNExpected(t.fNExpected)
 227   ,fNLast(t.fNLast)
 228   ,fNExpectedLast(t.fNExpectedLast)
 229   ,fNdedx(t.fNdedx)
 230   ,fChi2Last(t.fChi2Last)
 231   ,fBackupTrack(0x0)
 232 {
 233   //
 234   // Copy constructor.
 235   //
 236
 237   Int_t  i = 0;
 238   Int_t  j = 0;
 239   UInt_t k = 0;
 240
 241   for (i = 0; i < kNplane; i++) {
 242     for (j = 0; j < kNslice; j++) {
 243       fdEdxPlane[i][j] = t.fdEdxPlane[i][j];
 244     }
 245     fTimBinPlane[i] = t.fTimBinPlane[i];
 246     fTracklets[i]   = t.fTracklets[i];
 247   }
 248
 249   Int_t n = t.GetNumberOfClusters();
 250   for (i = 0; i < n; i++) {
 251     fIndex[i]       = t.fIndex[i];
 252     fIndexBackup[i] = t.fIndex[i];
 253     fdQdl[i]        = t.fdQdl[i];
 254   }
 255   for (k = n; k < kMAXCLUSTERSPERTRACK; k++) {
 256     fdQdl[k]        = 0;
 257     fIndex[k]       = 0;
 258     fIndexBackup[k] = 0;
 259   }
 260
 261   for (i = 0; i < 6; i++) {
 262     fTracklets[i] = t.fTracklets[i];
 263   }
 264
 265   for (i = 0; i < 3; i++) {
 266     fBudget[i]    = t.fBudget[i];
 267   }
 268
 269 }
 270
 271 //_____________________________________________________________________________
 272 AliTRDtrack::AliTRDtrack(const AliKalmanTrack &t, Double_t alpha)
 273   :AliKalmanTrack(t)
 274   ,fSeedLab(-1)
 275   ,fdEdx(t.GetPIDsignal())
 276   ,fDE(0)
 277   ,fAlpha(alpha)
 278   ,fX(0)
 279   ,fStopped(kFALSE)
 280   ,fY(0)
 281   ,fZ(0)
 282   ,fE(0)
 283   ,fT(0)
 284   ,fC(0)
 285   ,fCyy(0)
 286   ,fCzy(0)
 287   ,fCzz(0)
 288   ,fCey(0)
 289   ,fCez(0)
 290   ,fCee(0)
 291   ,fCty(0)
 292   ,fCtz(0)
 293   ,fCte(0)
 294   ,fCtt(0)
 295   ,fCcy(0)
 296   ,fCcz(0)
 297   ,fCce(0)
 298   ,fCct(0)
 299   ,fCcc(0)
 300   ,fLhElectron(0.0)
 301   ,fNWrong(0)
 302   ,fNRotate(0)
 303   ,fNCross(0)
 304   ,fNExpected(0)
 305   ,fNLast(0)
 306   ,fNExpectedLast(0)
 307   ,fNdedx(0)
 308   ,fChi2Last(0.0)
 309   ,fBackupTrack(0x0)
 310 {
 311   //
 312   // Constructor from AliTPCtrack or AliITStrack .
 313   //
 314
 315   Int_t  i = 0;
 316   Int_t  j = 0;
 317   UInt_t k = 0;
 318
 319   SetChi2(0.0);
 320   SetNumberOfClusters(0);
 321
 322   for (i = 0; i < kNplane; i++) {
 323     for (j = 0; j < kNslice; j++) {
 324       fdEdxPlane[i][j] = 0.0;
 325     }
 326     fTimBinPlane[i] = -1;
 327   }
 328
 329   if      (fAlpha < -TMath::Pi()) {
 330     fAlpha += 2.0 * TMath::Pi();
 331   }
 332   else if (fAlpha >= TMath::Pi()) {
 333     fAlpha -= 2.0 * TMath::Pi();
 334   }
 335
 336   Double_t x;
 337   Double_t p[5];
 338   t.GetExternalParameters(x,p);
 339   fX = x;
 340   fY = p[0];
 341   fZ = p[1];
 342   fT = p[3];
 343   x  = GetLocalConvConst();
 344   fC = p[4] / x;
 345   fE = fC * fX - p[2];
 346
 347   // Conversion of the covariance matrix
 348   Double_t c[15];
 349   t.GetExternalCovariance(c);
 350   c[10] /= x;
 351   c[11] /= x;
 352   c[12] /= x;
 353   c[13] /= x;
 354   c[14] /= x*x;
 355
 356   Double_t c22 = fX*fX * c[14] - 2.0*fX*c[12] + c[ 5];
 357   Double_t c32 = fX    * c[13] -        c[ 8];
 358   Double_t c20 = fX    * c[10] -        c[ 3];
 359   Double_t c21 = fX    * c[11] -        c[ 4];
 360   Double_t c42 = fX    * c[14] -        c[12];
 361
 362   fCyy = c[ 0];
 363   fCzy = c[ 1];  fCzz = c[ 2];
 364   fCey = c20;    fCez = c21;    fCee = c22;
 365   fCty = c[ 6];  fCtz = c[ 7];  fCte = c32;  fCtt = c[ 9];
 366   fCcy = c[10];  fCcz = c[11];  fCce = c42;  fCct = c[13];  fCcc = c[14];
 367
 368   for (k = 0; k < kMAXCLUSTERSPERTRACK; k++) {
 369     fdQdl[k]        = 0;
 370     fIndex[k]       = 0;
 371     fIndexBackup[k] = 0;
 372   }
 373
 374   for (i = 0; i < 3; i++) {
 375     fBudget[i]      = 0;
 376   }
 377
 378 }
 379
 380 //_____________________________________________________________________________
 381 AliTRDtrack::AliTRDtrack(const AliESDtrack &t)
 382   :AliKalmanTrack()
 383   ,fSeedLab(-1)
 384   ,fdEdx(t.GetTRDsignal())
 385   ,fDE(0)
 386   ,fAlpha(t.GetAlpha())
 387   ,fX(0)
 388   ,fStopped(kFALSE)
 389   ,fY(0)
 390   ,fZ(0)
 391   ,fE(0)
 392   ,fT(0)
 393   ,fC(0)
 394   ,fCyy(1e10)
 395   ,fCzy(0)
 396   ,fCzz(1e10)
 397   ,fCey(0)
 398   ,fCez(0)
 399   ,fCee(1e10)
 400   ,fCty(0)
 401   ,fCtz(0)
 402   ,fCte(0)
 403   ,fCtt(1e10)
 404   ,fCcy(0)
 405   ,fCcz(0)
 406   ,fCce(0)
 407   ,fCct(0)
 408   ,fCcc(1e10)
 409   ,fLhElectron(0.0)
 410   ,fNWrong(0)
 411   ,fNRotate(0)
 412   ,fNCross(0)
 413   ,fNExpected(0)
 414   ,fNLast(0)
 415   ,fNExpectedLast(0)
 416   ,fNdedx(0)
 417   ,fChi2Last(0.0)
 418   ,fBackupTrack(0x0)
 419 {
 420   //
 421   // Constructor from AliESDtrack
 422   //
 423
 424   Int_t  i = 0;
 425   Int_t  j = 0;
 426   UInt_t k = 0;
 427
 428   SetLabel(t.GetLabel());
 429   SetChi2(0.);
 430   SetMass(t.GetMass());
 431   SetNumberOfClusters(t.GetTRDclusters(fIndex));
 432
 433   Int_t ncl = t.GetTRDclusters(fIndexBackup);
 434   for (k = ncl; k < kMAXCLUSTERSPERTRACK; k++) {
 435     fIndexBackup[k] = 0;
 436     fIndex[k]       = 0;
 437   }
 438
 439   for (i = 0; i < kNplane; i++) {
 440     for (j = 0; j < kNslice; j++) {
 441       fdEdxPlane[i][j] = t.GetTRDsignals(i,j);
 442     }
 443     fTimBinPlane[i] = t.GetTRDTimBin(i);
 444   }
 445
 446   if      (fAlpha <  -TMath::Pi()) {
 447     fAlpha += 2.0 * TMath::Pi();
 448   }
 449   else if (fAlpha >=  TMath::Pi()) {
 450     fAlpha -= 2.0 * TMath::Pi();
 451   }
 452
 453   // Conversion of the covariance matrix
 454   Double_t x;
 455   Double_t p[5];
 456   t.GetExternalParameters(x,p);
 457   Double_t c[15];
 458   t.GetExternalCovariance(c);
 459   if (t.GetStatus() & AliESDtrack::kTRDbackup) {
 460     t.GetOuterExternalParameters(fAlpha,x,p);
 461     t.GetOuterExternalCovariance(c);
 462     if      (fAlpha <  -TMath::Pi()) {
 463       fAlpha += 2.0 * TMath::Pi();
 464     }
 465     else if (fAlpha >=  TMath::Pi()) {
 466       fAlpha -= 2.0 * TMath::Pi();
 467     }
 468   }
 469
 470   fX = x;
 471   fY = p[0];
 472   fZ = p[1];
 473   SaveLocalConvConst();
 474   fT = p[3];
 475   x  = GetLocalConvConst();
 476   fC = p[4] / x;
 477   fE = fC*fX - p[2];
 478
 479   c[10] /= x;
 480   c[11] /= x;
 481   c[12] /= x;
 482   c[13] /= x;
 483   c[14] /= x*x;
 484
 485   Double_t c22 = fX*fX * c[14] - 2.0*fX*c[12] + c[ 5];
 486   Double_t c32 = fX    * c[13] - c[ 8];
 487   Double_t c20 = fX    * c[10] - c[ 3];
 488   Double_t c21 = fX    * c[11] - c[ 4];
 489   Double_t c42 = fX    * c[14] - c[12];
 490
 491   fCyy = c[ 0];
 492   fCzy = c[ 1];  fCzz = c[ 2];
 493   fCey = c20;    fCez = c21;    fCee = c22;
 494   fCty = c[ 6];  fCtz = c[ 7];  fCte = c32;  fCtt = c[ 9];
 495   fCcy = c[10];  fCcz = c[11];  fCce = c42;  fCct = c[13];  fCcc = c[14];
 496
 497   for (k = 0; k < kMAXCLUSTERSPERTRACK; k++) {
 498     fdQdl[k] = 0;
 499     //fIndex[k] = 0; //MI store indexes
 500   }
 501
 502   for (i = 0; i < 3; i++) {
 503     fBudget[i] = 0;
 504   }
 505   if ((t.GetStatus() & AliESDtrack::kTIME) == 0) {
 506     return;
 507   }
 508
 509   StartTimeIntegral();
 510   Double_t times[10];
 511   t.GetIntegratedTimes(times);
 512   SetIntegratedTimes(times);
 513   SetIntegratedLength(t.GetIntegratedLength());
 514
 515 }
 516
 517 //____________________________________________________________________________
 518 AliTRDtrack::~AliTRDtrack()
 519 {
 520   //
 521   // Destructor
 522   //
 523
 524   if (fBackupTrack) {
 525     delete fBackupTrack;
 526   }
 527   fBackupTrack = 0;
 528
 529 }
 530
 531 //____________________________________________________________________________
 532 AliTRDtrack &AliTRDtrack::operator=(const AliTRDtrack &t)
 533 {
 534   //
 535   // Assignment operator
 536   //
 537
 538   fLhElectron    = 0.0;
 539   fNWrong        = 0;
 540   fStopped       = 0;
 541   fNRotate       = 0;
 542   fNExpected     = 0;
 543   fNExpectedLast = 0;
 544   fNdedx         = 0;
 545   fNCross        = 0;
 546   fNLast         = 0;
 547   fChi2Last      = 0;
 548   fBackupTrack   = 0;
 549
 550   fAlpha = t.GetAlpha();
 551   if      (fAlpha <  -TMath::Pi()) {
 552     fAlpha += 2.0 * TMath::Pi();
 553   }
 554   else if (fAlpha >=  TMath::Pi()) {
 555     fAlpha -= 2.0 * TMath::Pi();
 556   }
 557
 558   return *this;
 559
 560 }
 561
 562 //____________________________________________________________________________
 563 Float_t AliTRDtrack::StatusForTOF()
 564 {
 565   //
 566   // Defines the status of the TOF extrapolation
 567   //
 568
 569   // Definition of res ????
 570   Float_t res = (0.2 + 0.8 * (fN / (fNExpected + 5.0)))
 571               * (0.4 + 0.6 * fTracklets[5].GetN() / 20.0);
 572   res *= (0.25 + 0.8 * 40.0 / (40.0 + fBudget[2]));
 573   return res;
 574
 575   // This part of the function is never reached ????
 576   // What defines these parameters ????
 577   Int_t status = 0;
 578   if (GetNumberOfClusters()                <  20)  return 0;
 579   if ((fN                                  > 110) &&
 580       (fChi2/(Float_t(fN))                 < 3.0)) return 3; // Gold
 581   if ((fNLast                              >  30) &&
 582       (fChi2Last/(Float_t(fNLast))         < 3.0)) return 3; // Gold
 583   if ((fNLast                              >  20) &&
 584       (fChi2Last/(Float_t(fNLast))         < 2.0)) return 3; // Gold
 585   if ((fNLast/(fNExpectedLast+3.0)         > 0.8) &&
 586       (fChi2Last/Float_t(fNLast)           < 5.0) &&
 587       (fNLast                              >  20)) return 2; // Silver
 588   if ((fNLast                              >   5) &&
 589       (((fNLast+1.0)/(fNExpectedLast+1.0)) > 0.8) &&
 590       (fChi2Last/(fNLast-5.0)              < 6.0)) return 1;
 591
 592   return status;
 593
 594 }
 595
 596 //_____________________________________________________________________________
 597 void AliTRDtrack::GetExternalCovariance(Double_t cc[15]) const
 598 {
 599   //
 600   // This function returns external representation of the covriance matrix.
 601   //
 602
 603   Double_t a   = GetLocalConvConst();
 604
 605   Double_t c22 = fX*fX*fCcc  - 2.0*fX*fCce+fCee;
 606   Double_t c32 = fX*fCct-fCte;
 607   Double_t c20 = fX*fCcy-fCey;
 608   Double_t c21 = fX*fCcz-fCez;
 609   Double_t c42 = fX*fCcc-fCce;
 610
 611   cc[ 0]=fCyy;
 612   cc[ 1]=fCzy;   cc[ 2]=fCzz;
 613   cc[ 3]=c20;    cc[ 4]=c21;    cc[ 5]=c22;
 614   cc[ 6]=fCty;   cc[ 7]=fCtz;   cc[ 8]=c32;   cc[ 9]=fCtt;
 615   cc[10]=fCcy*a; cc[11]=fCcz*a; cc[12]=c42*a; cc[13]=fCct*a; cc[14]=fCcc*a*a;
 616
 617 }
 618
 619 //_____________________________________________________________________________
 620 void AliTRDtrack::GetCovariance(Double_t cc[15]) const
 621 {
 622   //
 623   // Returns the track covariance matrix
 624   //
 625
 626   cc[ 0]=fCyy;
 627   cc[ 1]=fCzy; cc[ 2]=fCzz;
 628   cc[ 3]=fCey; cc[ 4]=fCez; cc[ 5]=fCee;
 629   cc[ 6]=fCcy; cc[ 7]=fCcz; cc[ 8]=fCce; cc[ 9]=fCcc;
 630   cc[10]=fCty; cc[11]=fCtz; cc[12]=fCte; cc[13]=fCct; cc[14]=fCtt;
 631
 632 }
 633
 634 //_____________________________________________________________________________
 635 Int_t AliTRDtrack::Compare(const TObject *o) const
 636 {
 637   //
 638   // Compares tracks according to their Y2 or curvature
 639   //
 640
 641   AliTRDtrack *t = (AliTRDtrack *) o;
 642   //  Double_t co=t->GetSigmaY2();
 643   //  Double_t c =GetSigmaY2();
 644
 645   Double_t co = TMath::Abs(t->GetC());
 646   Double_t c  = TMath::Abs(GetC());
 647
 648   if      (c > co) {
 649     return 1;
 650   }
 651   else if (c < co) {
 652     return -1;
 653   }
 654   return 0;
 655
 656 }
 657
 658 //_____________________________________________________________________________
 659 void AliTRDtrack::CookdEdx(Double_t low, Double_t up)
 660 {
 661   //
 662   // Calculates the truncated dE/dx within the "low" and "up" cuts.
 663   //
 664
 665   Int_t   i  = 0;
 666
 667   // Array to sort the dEdx values according to amplitude
 668   Float_t sorted[kMAXCLUSTERSPERTRACK];
 669
 670   // Number of clusters used for dedx
 671   Int_t   nc = fNdedx;
 672
 673   // Require at least 10 clusters for a dedx measurement
 674   if (nc < 10) {
 675     SetdEdx(0);
 676     return;
 677   }
 678
 679   // Lower and upper bound
 680   Int_t nl = Int_t(low * nc);
 681   Int_t nu = Int_t( up * nc);
 682
 683   // Can fdQdl be negative ????
 684   for (i = 0; i < nc; i++) {
 685     sorted[i] = TMath::Abs(fdQdl[i]);
 686   }
 687
 688   // Sort the dedx values by amplitude
 689   Int_t *index = new Int_t[nc];
 690   TMath::Sort(nc,sorted,index,kFALSE);
 691
 692   // Sum up the truncated charge between nl and nu
 693   Float_t dedx = 0.0;
 694   for (i = nl; i <= nu; i++) {
 695     dedx += sorted[index[i]];
 696   }
 697   dedx /= (nu - nl + 1.0);
 698   SetdEdx(dedx);
 699
 700   delete [] index;
 701
 702 }
 703
 704 //_____________________________________________________________________________
 705 Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho)
 706 {
 707   //
 708   // Propagates a track of particle with mass=pm to a reference plane
 709   // defined by x=xk through media of density=rho and radiationLength=x0
 710   //
 711
 712   if (xk == fX) {
 713     return 1;
 714   }
 715
 716   if (TMath::Abs(fC*xk - fE) >= 0.9) {
 717     return 0;
 718   }
 719
 720   Double_t lcc = GetLocalConvConst();
 721
 722   Double_t oldX = fX;
 723   Double_t oldY = fY;
 724   Double_t oldZ = fZ;
 725
 726   Double_t x1 = fX;
 727   Double_t x2 = x1 + (xk - x1);
 728   Double_t dx = x2 - x1;
 729   Double_t y1 = fY;
 730   Double_t z1 = fZ;
 731   Double_t c1 = fC*x1 - fE;
 732   if ((c1*c1) > 1) {
 733     return 0;
 734   }
 735
 736   Double_t r1 = TMath::Sqrt(1.0 - c1*c1);
 737   Double_t c2 = fC*x2 - fE;
 738   if ((c2*c2) > 1) {
 739     return 0;
 740   }
 741
 742   Double_t r2 = TMath::Sqrt(1.0 - c2*c2);
 743
 744   fY += dx*(c1+c2) / (r1+r2);
 745   fZ += dx*(c1+c2) / (c1*r2 + c2*r1) * fT;
 746
 747   // f = F - 1
 748   Double_t rr  =  r1+r2;
 749   Double_t cc  =  c1+c2;
 750   Double_t xx  =  x1+x2;
 751   Double_t f02 = -dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
 752   Double_t f04 =  dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
 753   Double_t cr  =  c1*r2+c2*r1;
 754   Double_t f12 = -dx*fT*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
 755   Double_t f13 =  dx*cc/cr;
 756   Double_t f14 =  dx*fT*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
 757
 758   // b = C*ft
 759   Double_t b00 = f02*fCey + f04*fCcy;
 760   Double_t b01 = f12*fCey + f14*fCcy + f13*fCty;
 761   Double_t b10 = f02*fCez + f04*fCcz;
 762   Double_t b11 = f12*fCez + f14*fCcz + f13*fCtz;
 763   Double_t b20 = f02*fCee + f04*fCce;
 764   Double_t b21 = f12*fCee + f14*fCce + f13*fCte;
 765   Double_t b30 = f02*fCte + f04*fCct;
 766   Double_t b31 = f12*fCte + f14*fCct + f13*fCtt;
 767   Double_t b40 = f02*fCce + f04*fCcc;
 768   Double_t b41 = f12*fCce + f14*fCcc + f13*fCct;
 769
 770   // a = f*b = f*C*ft
 771   Double_t a00 = f02*b20 + f04*b40;
 772   Double_t a01 = f02*b21 + f04*b41;
 773   Double_t a11 = f12*b21 + f14*b41 + f13*b31;
 774
 775   // F*C*Ft = C + (a + b + bt)
 776   fCyy += a00 + 2.0*b00;
 777   fCzy += a01 + b01 + b10;
 778   fCey += b20;
 779   fCty += b30;
 780   fCcy += b40;
 781   fCzz += a11 + 2.0*b11;
 782   fCez += b21;
 783   fCtz += b31;
 784   fCcz += b41;
 785
 786   fX = x2;
 787
 788   // Change of the magnetic field
 789   SaveLocalConvConst();
 790   cc =  fC;
 791   fC *= lcc / GetLocalConvConst();
 792   fE += fX  * (fC-cc);
 793
 794   // Multiple scattering
 795   // What is 14.1 ????
 796   Double_t d      = TMath::Sqrt((x1-fX)*(x1-fX) + (y1-fY)*(y1-fY) + (z1-fZ)*(z1-fZ));
 797   Double_t p2     = (1.0 + GetTgl()*GetTgl()) / (Get1Pt()*Get1Pt());
 798   Double_t beta2  = p2 / (p2 + GetMass()*GetMass());
 799   Double_t theta2 = 14.1*14.1 / (beta2*p2*1e6) * d / x0 * rho;
 800   Double_t ey     = fC*fX - fE;
 801   Double_t ez     = fT;
 802   Double_t xz     = fC*ez;
 803   Double_t zz1    = ez*ez + 1.0;
 804   Double_t xy     = fE + ey;
 805
 806   fCee += (2.0*ey*ez*ez*fE + 1.0 - ey*ey + ez*ez + fE*fE*ez*ez) * theta2;
 807   fCte += ez*zz1*xy*theta2;
 808   fCtt += zz1*zz1*theta2;
 809   fCce += xz*ez*xy*theta2;
 810   fCct += xz*zz1*theta2;
 811   fCcc += xz*xz*theta2;
 812
 813   // Energy losses
 814   // What is 5940.0 ???? and 0.153e-3 ????
 815   if ((5940.0*beta2 / (1.0 - beta2 + 1e-10) - beta2) < 0.0) {
 816     return 0;
 817   }
 818   Double_t dE     = 0.153e-3/beta2 * (TMath::Log(5940.0*beta2 / (1.0 - beta2 + 1e-10)) - beta2)
 819                   * d * rho;
 820   Float_t  budget = d * rho;
 821   fBudget[0] +=budget;
 822
 823   // Suspicious part - think about it ????
 824   Double_t kinE =  TMath::Sqrt(p2);
 825   if (dE > 0.8 * kinE) {
 826     dE = 0.8 * kinE;
 827   }
 828   if (dE <        0.0) {
 829     dE = 0.0;       // Not valid region for Bethe bloch
 830   }
 831   fDE += dE;
 832   if (x1 < x2) {
 833     dE = -dE;
 834   }
 835   cc  = fC;
 836   fC *= (1.0 - TMath::Sqrt(p2 + GetMass()*GetMass()) / p2 * dE);
 837   fE += fX * (fC - cc);
 838
 839   // Energy loss fluctuation
 840   // Why 0.07 ????
 841   Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE));
 842   Double_t sigmac2 = sigmade*sigmade * fC*fC * (p2 + GetMass()*GetMass()) / (p2*p2);
 843   fCcc += sigmac2;
 844   fCee += fX*fX * sigmac2;
 845
 846   // Track time measurement
 847   if (x1 < x2) {
 848     if (IsStartedTimeIntegral()) {
 849       Double_t l2 = TMath::Sqrt((fX-oldX)*(fX-oldX)
 850                               + (fY-oldY)*(fY-oldY)
 851                               + (fZ-oldZ)*(fZ-oldZ));
 852       if (TMath::Abs(l2*fC) > 0.0001){
 853         // <ake correction for curvature if neccesary
 854         l2 = 0.5 * TMath::Sqrt((fX-oldX)*(fX-oldX)
 855                              + (fY-oldY)*(fY-oldY));
 856         l2 = 2.0 * TMath::ASin(l2 * fC) / fC;
 857         l2 = TMath::Sqrt(l2*l2 + (fZ-oldZ)*(fZ-oldZ));
 858       }
 859       AddTimeStep(l2);
 860     }
 861   }
 862
 863   return 1;
 864
 865 }
 866
 867 //_____________________________________________________________________________
 868 Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq
 869                         , UInt_t index, Double_t h01)
 870 {
 871   //
 872   // Assignes a found cluster to the track and updates track information
 873   //
 874
 875   Bool_t fNoTilt = kTRUE;
 876   // What is 0.003 ????
 877   if (TMath::Abs(h01) > 0.003) {
 878     fNoTilt = kFALSE;
 879   }
 880
 881   // Add angular effect to the error contribution
 882   Float_t tangent2  = (fC*fX-fE) * (fC*fX-fE);
 883   if (tangent2 < 0.90000){
 884     tangent2 = tangent2 / (1.0 - tangent2);
 885   }
 886   // What is 0.04 ????
 887   Float_t errang    = tangent2 * 0.04;
 888   Float_t padlength = TMath::Sqrt(c->GetSigmaZ2() * 12.0);
 889
 890   Double_t r00 = c->GetSigmaY2() + errang;
 891   Double_t r01 = 0.0;
 892   Double_t r11 = c->GetSigmaZ2() * 100.0;
 893   r00 += fCyy;
 894   r01 += fCzy;
 895   r11 += fCzz;
 896   Double_t det = r00*r11 - r01*r01;
 897   Double_t tmp = r00;
 898   r00 =  r11 / det;
 899   r11 =  tmp / det;
 900   r01 = -r01 / det;
 901
 902   Double_t k00 = fCyy*r00 + fCzy*r01;
 903   Double_t k01 = fCyy*r01 + fCzy*r11;
 904   Double_t k10 = fCzy*r00 + fCzz*r01;
 905   Double_t k11 = fCzy*r01 + fCzz*r11;
 906   Double_t k20 = fCey*r00 + fCez*r01;
 907   Double_t k21 = fCey*r01 + fCez*r11;
 908   Double_t k30 = fCty*r00 + fCtz*r01;
 909   Double_t k31 = fCty*r01 + fCtz*r11;
 910   Double_t k40 = fCcy*r00 + fCcz*r01;
 911   Double_t k41 = fCcy*r01 + fCcz*r11;
 912
 913   Double_t dy  = c->GetY() - fY;
 914   Double_t dz  = c->GetZ() - fZ;
 915   Double_t cur = fC + k40*dy + k41*dz;
 916   Double_t eta = fE + k20*dy + k21*dz;
 917
 918   if (fNoTilt) {
 919
 920     if (TMath::Abs(cur*fX-eta) >= 0.9) {
 921       return 0;
 922     }
 923     fY += k00*dy + k01*dz;
 924     fZ += k10*dy + k11*dz;
 925     fE  = eta;
 926     fC  = cur;
 927
 928   }
 929   else {
 930
 931     // Empirical factor set by C.Xu in the first tilt version
 932     // Is this factor still ok ????
 933     Double_t xuFactor = 100.0;
 934     dy   = c->GetY() - fY;
 935     dz   = c->GetZ() - fZ;
 936     dy   = dy + h01*dz;
 937
 938     Float_t add = 0.0;
 939     if (TMath::Abs(dz) > padlength/2.0) {
 940       Float_t dy2  = c->GetY() - fY;
 941       Float_t sign = (dz > 0.0) ? -1.0 : 1.0;
 942       dy2 += h01 * sign * padlength/2.0;
 943       dy   = dy2;
 944       add  = 0.0;
 945     }
 946
 947     r00  = c->GetSigmaY2() + errang + add;
 948     r01  = 0.0;
 949     r11  = c->GetSigmaZ2() * xuFactor;
 950     r00 += (fCyy + 2.0*h01*fCzy + h01*h01*fCzz);
 951     r01 += (fCzy + h01*fCzz);
 952     r11 += fCzz;
 953
 954     det  =  r00*r11 - r01*r01;
 955     tmp  =  r00;
 956     r00  =  r11/det;
 957     r11  =  tmp/det;
 958     r01  = -r01/det;
 959
 960     k00  = fCyy*r00 + fCzy * (r01 + h01*r00);
 961     k01  = fCyy*r01 + fCzy * (r11 + h01*r01);
 962     k10  = fCzy*r00 + fCzz * (r01 + h01*r00);
 963     k11  = fCzy*r01 + fCzz * (r11 + h01*r01);
 964     k20  = fCey*r00 + fCez * (r01 + h01*r00);
 965     k21  = fCey*r01 + fCez * (r11 + h01*r01);
 966     k30  = fCty*r00 + fCtz * (r01 + h01*r00);
 967     k31  = fCty*r01 + fCtz * (r11 + h01*r01);
 968     k40  = fCcy*r00 + fCcz * (r01 + h01*r00);
 969     k41  = fCcy*r01 + fCcz * (r11 + h01*r01);
 970
 971     cur  = fC + k40*dy + k41*dz;
 972     eta  = fE + k20*dy + k21*dz;
 973     if (TMath::Abs(cur*fX - eta) >= 0.9) {
 974       return 0;
 975     }
 976
 977     fY  += k00*dy + k01*dz;
 978     fZ  += k10*dy + k11*dz;
 979     fE   = eta;
 980     fT  += k30*dy + k31*dz;
 981     fC   = cur;
 982
 983     k01 += h01*k00;
 984     k11 += h01*k10;
 985     k21 += h01*k20;
 986     k31 += h01*k30;
 987     k41 += h01*k40;
 988
 989   }
 990
 991   Double_t c01 = fCzy;
 992   Double_t c02 = fCey;
 993   Double_t c03 = fCty;
 994   Double_t c04 = fCcy;
 995   Double_t c12 = fCez;
 996   Double_t c13 = fCtz;
 997   Double_t c14 = fCcz;
 998
 999   fCyy -= k00*fCyy + k01*fCzy;
1000   fCzy -= k00*c01  + k01*fCzz;
1001   fCey -= k00*c02  + k01*c12;
1002   fCty -= k00*c03  + k01*c13;
1003   fCcy -= k00*c04  + k01*c14;
1004
1005   fCzz -= k10*c01  + k11*fCzz;
1006   fCez -= k10*c02  + k11*c12;
1007   fCtz -= k10*c03  + k11*c13;
1008   fCcz -= k10*c04  + k11*c14;
1009
1010   fCee -= k20*c02  + k21*c12;
1011   fCte -= k20*c03  + k21*c13;
1012   fCce -= k20*c04  + k21*c14;
1013
1014   fCtt -= k30*c03  + k31*c13;
1015   fCct -= k40*c03  + k41*c13;
1016
1017   fCcc -= k40*c04  + k41*c14;
1018
1019   Int_t n = GetNumberOfClusters();
1020   fIndex[n] = index;
1021   SetNumberOfClusters(n+1);
1022
1023   SetChi2(GetChi2()+chisq);
1024
1025   return 1;
1026
1027 }
1028
1029 //_____________________________________________________________________________
1030 Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq
1031                           , UInt_t index, Double_t h01
1032                           , Int_t /*plane*/)
1033 {
1034   //
1035   // Assignes found cluster to the track and updates track information
1036   //
1037
1038   Bool_t fNoTilt = kTRUE;
1039   if (TMath::Abs(h01) > 0.003) {
1040     fNoTilt = kFALSE;
1041   }
1042
1043   //
1044   // Add angular effect to the error contribution and make correction
1045   // Still needed ????
1046   // AliTRDclusterCorrection *corrector = AliTRDclusterCorrection::GetCorrection();
1047   //
1048
1049   Double_t tangent2 = (fC*fX-fE) * (fC*fX-fE);
1050   if (tangent2 < 0.9) {
1051     tangent2 = tangent2 / (1.0 - tangent2);
1052   }
1053   Double_t tangent  = TMath::Sqrt(tangent2);
1054   if ((fC*fX-fE) < 0.0) {
1055     tangent *= -1.0;
1056   }
1057
1058   // Where are the parameters from ????
1059   Double_t errang = tangent2 * 0.04;
1060   Double_t errsys = 0.025*0.025 * 20.0;  // Systematic error part
1061
1062   Float_t  extend = 1.0;
1063   if (c->GetNPads() == 4) extend = 2.0;
1064
1065   /////////////////////////////////////////////////////////////////////////////
1066   //
1067   // Is this still needed or will it be needed ????
1068   //
1069   //if (c->GetNPads() == 5) extend = 3.0;
1070   //if (c->GetNPads() == 6) extend = 3.0;
1071   //if (c->GetQ() < 15) {
1072   //  return 1;
1073   //}
1074   //
1075   // Will this be needed ????
1076   /*
1077   if (corrector !=0 ) {
1078   //if (0){
1079     correction = corrector->GetCorrection(plane,c->GetLocalTimeBin(),tangent);
1080     if (TMath::Abs(correction)>0){
1081       //if we have info
1082       errang     = corrector->GetSigma(plane,c->GetLocalTimeBin(),tangent);
1083       errang    *= errang;
1084       errang    += tangent2*0.04;
1085     }
1086   }
1087   */
1088   //
1089   //  Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
1090   /////////////////////////////////////////////////////////////////////////////
1091
1092   Double_t r00 = (c->GetSigmaY2() + errang + errsys) * extend;
1093   Double_t r01 = 0.0;
1094   Double_t r11 = c->GetSigmaZ2()*10000.0;
1095   r00 += fCyy;
1096   r01 += fCzy;
1097   r11 += fCzz;
1098   Double_t det =r00*r11 - r01*r01;
1099   Double_t tmp =r00;
1100   r00  =  r11 / det;
1101   r11  =  tmp / det;
1102   r01  = -r01 / det;
1103
1104   Double_t k00 = fCyy*r00 + fCzy*r01;
1105   Double_t k01 = fCyy*r01 + fCzy*r11;
1106   Double_t k10 = fCzy*r00 + fCzz*r01;
1107   Double_t k11 = fCzy*r01 + fCzz*r11;
1108   Double_t k20 = fCey*r00 + fCez*r01;
1109   Double_t k21 = fCey*r01 + fCez*r11;
1110   Double_t k30 = fCty*r00 + fCtz*r01;
1111   Double_t k31 = fCty*r01 + fCtz*r11;
1112   Double_t k40 = fCcy*r00 + fCcz*r01;
1113   Double_t k41 = fCcy*r01 + fCcz*r11;
1114
1115   Double_t dy  = c->GetY() - fY;
1116   Double_t dz  = c->GetZ() - fZ;
1117   Double_t cur = fC + k40*dy + k41*dz;
1118   Double_t eta = fE + k20*dy + k21*dz;
1119
1120   if (fNoTilt) {
1121
1122     if (TMath::Abs(cur*fX - eta) >= 0.9) {
1123       return 0;
1124     }
1125
1126     fY += k00*dy + k01*dz;
1127     fZ += k10*dy + k11*dz;
1128     fE  = eta;
1129     fC  = cur;
1130
1131   }
1132   else {
1133
1134     Double_t padlength = TMath::Sqrt(c->GetSigmaZ2() * 12.0);
1135     // Empirical factor set by C.Xu in the first tilt version
1136     Double_t xuFactor  = 1000.0;
1137
1138     dy = c->GetY() - fY;
1139     dz = c->GetZ() - fZ;
1140     //dy = dy + h01*dz + correction; // Still needed ????
1141
1142     Double_t tiltdz = dz;
1143     if (TMath::Abs(tiltdz) > padlength/2.0) {
1144       tiltdz = TMath::Sign(padlength/2.0,dz);
1145     }
1146     dy = dy + h01*tiltdz;
1147
1148     Double_t add = 0.0;
1149     if (TMath::Abs(dz) > padlength/2.0) {
1150       //Double_t dy2 = c->GetY() - fY;     // Still needed ????
1151       //Double_t sign = (dz>0) ? -1.: 1.;
1152       //dy2-=h01*sign*padlength/2.;
1153       //dy = dy2;
1154       add = 1.0;
1155     }
1156
1157     Double_t s00 = (c->GetSigmaY2() + errang) * extend + errsys + add;  // Error pad
1158     Double_t s11 = c->GetSigmaZ2() * xuFactor;                          // Error pad-row
1159     //
1160     r00  = fCyy + 2*fCzy*h01 + fCzz*h01*h01 + s00;
1161     r01  = fCzy + fCzz*h01;
1162     r11  = fCzz + s11;
1163     det  = r00*r11 - r01*r01;
1164
1165     // Inverse matrix
1166     tmp  =  r00;
1167     r00  =  r11 / det;
1168     r11  =  tmp / det;
1169     r01  = -r01 / det;
1170
1171     // K matrix
1172     k00  = fCyy*r00 + fCzy * (r01 + h01*r00);
1173     k01  = fCyy*r01 + fCzy * (r11 + h01*r01);
1174     k10  = fCzy*r00 + fCzz * (r01 + h01*r00);
1175     k11  = fCzy*r01 + fCzz * (r11 + h01*r01);
1176     k20  = fCey*r00 + fCez * (r01 + h01*r00);
1177     k21  = fCey*r01 + fCez * (r11 + h01*r01);
1178     k30  = fCty*r00 + fCtz * (r01 + h01*r00);
1179     k31  = fCty*r01 + fCtz * (r11 + h01*r01);
1180     k40  = fCcy*r00 + fCcz * (r01 + h01*r00);
1181     k41  = fCcy*r01 + fCcz * (r11 + h01*r01);
1182
1183     // Update measurement
1184     cur  = fC + k40*dy + k41*dz;
1185     eta  = fE + k20*dy + k21*dz;
1186     if (TMath::Abs(cur*fX - eta) >= 0.9) {
1187       return 0;
1188     }
1189     fY  += k00*dy + k01*dz;
1190     fZ  += k10*dy + k11*dz;
1191     fE   = eta;
1192     fT  += k30*dy + k31*dz;
1193     fC   = cur;
1194
1195     k01 += h01*k00;
1196     k11 += h01*k10;
1197     k21 += h01*k20;
1198     k31 += h01*k30;
1199     k41 += h01*k40;
1200
1201   }
1202
1203   // Update the covariance matrix
1204   Double_t oldyy = fCyy;
1205   Double_t oldzz = fCzz;
1206   //Double_t oldee = fCee;
1207   //Double_t oldcc = fCcc;
1208   Double_t oldzy = fCzy;
1209   Double_t oldey = fCey;
1210   Double_t oldty = fCty;
1211   Double_t oldcy = fCcy;
1212   Double_t oldez = fCez;
1213   Double_t oldtz = fCtz;
1214   Double_t oldcz = fCcz;
1215   //Double_t oldte = fCte;
1216   //Double_t oldce = fCce;
1217   //Double_t oldct = fCct;
1218
1219   fCyy -= k00*oldyy + k01*oldzy;
1220   fCzy -= k10*oldyy + k11*oldzy;
1221   fCey -= k20*oldyy + k21*oldzy;
1222   fCty -= k30*oldyy + k31*oldzy;
1223   fCcy -= k40*oldyy + k41*oldzy;
1224
1225   fCzz -= k10*oldzy + k11*oldzz;
1226   fCez -= k20*oldzy + k21*oldzz;
1227   fCtz -= k30*oldzy + k31*oldzz;
1228   fCcz -= k40*oldzy + k41*oldzz;
1229
1230   fCee -= k20*oldey + k21*oldez;
1231   fCte -= k30*oldey + k31*oldez;
1232   fCce -= k40*oldey + k41*oldez;
1233
1234   fCtt -= k30*oldty + k31*oldtz;
1235   fCct -= k40*oldty + k41*oldtz;
1236
1237   fCcc -= k40*oldcy + k41*oldcz;
1238
1239   Int_t n = GetNumberOfClusters();
1240   fIndex[n] = index;
1241   SetNumberOfClusters(n+1);
1242
1243   SetChi2(GetChi2() + chisq);
1244
1245   return 1;
1246
1247 }
1248
1249 //_____________________________________________________________________________
1250 Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet)
1251 {
1252   //
1253   // Assignes found tracklet to the track and updates track information
1254   //
1255
1256   Double_t r00 = (tracklet.GetTrackletSigma2());
1257   Double_t r01 = 0.0;
1258   Double_t r11 = 10000.0;
1259   r00 += fCyy;
1260   r01 += fCzy;
1261   r11 += fCzz;
1262
1263   Double_t det = r00*r11 - r01*r01;
1264   Double_t tmp = r00;
1265   r00 =  r11 / det;
1266   r11 =  tmp / det;
1267   r01 = -r01 / det;
1268
1269   Double_t dy  = tracklet.GetY() - fY;
1270   Double_t dz  = tracklet.GetZ() - fZ;
1271
1272   Double_t s00 = tracklet.GetTrackletSigma2();  // Error pad
1273   Double_t s11 = 100000.0;                      // Error pad-row
1274   Float_t  h01 = tracklet.GetTilt();
1275
1276   r00 = fCyy + fCzz*h01*h01 + s00;
1277   r01 = fCzy;
1278   r11 = fCzz + s11;
1279   det = r00*r11 - r01*r01;
1280
1281   // Inverse matrix
1282   tmp =  r00;
1283   r00 =  r11 / det;
1284   r11 =  tmp / det;
1285   r01 = -r01 / det;
1286
1287   // K matrix
1288   Double_t k00 = fCyy*r00 + fCzy*r01;
1289   Double_t k01 = fCyy*r01 + fCzy*r11;
1290   Double_t k10 = fCzy*r00 + fCzz*r01;
1291   Double_t k11 = fCzy*r01 + fCzz*r11;
1292   Double_t k20 = fCey*r00 + fCez*r01;
1293   Double_t k21 = fCey*r01 + fCez*r11;
1294   Double_t k30 = fCty*r00 + fCtz*r01;
1295   Double_t k31 = fCty*r01 + fCtz*r11;
1296   Double_t k40 = fCcy*r00 + fCcz*r01;
1297   Double_t k41 = fCcy*r01 + fCcz*r11;
1298
1299   // Update measurement
1300   Double_t cur = fC + k40*dy + k41*dz;
1301   Double_t eta = fE + k20*dy + k21*dz;
1302   if (TMath::Abs(cur*fX-eta) >= 0.9) {
1303     return 0;
1304   }
1305
1306   fY += k00*dy + k01*dz;
1307   fZ += k10*dy + k11*dz;
1308   fE  = eta;
1309   fT += k30*dy + k31*dz;
1310   fC  = cur;
1311
1312   // Update the covariance matrix
1313   Double_t oldyy = fCyy;
1314   Double_t oldzz = fCzz;
1315   //Double_t oldee = fCee;
1316   //Double_t oldcc = fCcc;
1317   Double_t oldzy = fCzy;
1318   Double_t oldey = fCey;
1319   Double_t oldty = fCty;
1320   Double_t oldcy = fCcy;
1321   Double_t oldez = fCez;
1322   Double_t oldtz = fCtz;
1323   Double_t oldcz = fCcz;
1324   //Double_t oldte = fCte;
1325   //Double_t oldce = fCce;
1326   //Double_t oldct = fCct;
1327
1328   fCyy -= k00*oldyy + k01*oldzy;
1329   fCzy -= k10*oldyy + k11*oldzy;
1330   fCey -= k20*oldyy + k21*oldzy;
1331   fCty -= k30*oldyy + k31*oldzy;
1332   fCcy -= k40*oldyy + k41*oldzy;
1333
1334   fCzz -= k10*oldzy + k11*oldzz;
1335   fCez -= k20*oldzy + k21*oldzz;
1336   fCtz -= k30*oldzy + k31*oldzz;
1337   fCcz -= k40*oldzy + k41*oldzz;
1338
1339   fCee -= k20*oldey + k21*oldez;
1340   fCte -= k30*oldey + k31*oldez;
1341   fCce -= k40*oldey + k41*oldez;
1342
1343   fCtt -= k30*oldty + k31*oldtz;
1344   fCct -= k40*oldty + k41*oldtz;
1345
1346   fCcc -= k40*oldcy + k41*oldcz;
1347
1348   return 1;
1349
1350 }
1351
1352 //_____________________________________________________________________________
1353 Int_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
1354 {
1355   //
1356   // Rotates the track parameters in the R*phi plane,
1357   // if the absolute rotation alpha is in global system.
1358   // Otherwise the rotation is relative to the current rotation angle
1359   //
1360
1361   if (absolute) {
1362     alpha -= fAlpha;
1363   }
1364   else{
1365     fNRotate++;
1366   }
1367
1368   fAlpha += alpha;
1369   if (fAlpha <  -TMath::Pi()) {
1370     fAlpha += 2.0 * TMath::Pi();
1371   }
1372   if (fAlpha >=  TMath::Pi()) {
1373     fAlpha -= 2.0 * TMath::Pi();
1374   }
1375
1376   Double_t x1 = fX;
1377   Double_t y1 = fY;
1378   Double_t ca = TMath::Cos(alpha);
1379   Double_t sa = TMath::Sin(alpha);
1380   Double_t r1 = fC*fX - fE;
1381
1382   fX =  x1*ca + y1*sa;
1383   fY = -x1*sa + y1*ca;
1384   if ((r1*r1) > 1.0) {
1385     return 0;
1386   }
1387   fE = fE*ca + (fC*y1 + TMath::Sqrt(1.0 - r1*r1)) * sa;
1388
1389   Double_t r2 = fC*fX - fE;
1390   if (TMath::Abs(r2) >= 0.9) {
1391     Int_t n = GetNumberOfClusters();
1392     if (n > 4) {
1393       AliError(Form("Rotation failed N = %d !\n",n));
1394     }
1395     return 0;
1396   }
1397
1398   if ((r2*r2) > 1.0) {
1399     return 0;
1400   }
1401   Double_t y0 = fY + TMath::Sqrt(1.0 - r2*r2) / fC;
1402   if ((fY-y0)*fC >= 0.0) {
1403     Int_t n = GetNumberOfClusters();
1404     if (n > 4) {
1405       AliError(Form("Rotation failed N = %d !\n",n));
1406     }
1407     return 0;
1408   }
1409
1410   // f = F - 1
1411   Double_t f00 = ca-1.0;
1412   Double_t f24 = (y1 - r1*x1/TMath::Sqrt(1.0 - r1*r1)) * sa;
1413   Double_t f20 = fC*sa;
1414   Double_t f22 = (ca + sa*r1/TMath::Sqrt(1.0 - r1*r1)) - 1.0;
1415
1416   // b = C*ft
1417   Double_t b00 = fCyy*f00;
1418   Double_t b02 = fCyy*f20 + fCcy*f24 + fCey*f22;
1419   Double_t b10 = fCzy*f00;
1420   Double_t b12 = fCzy*f20 + fCcz*f24 + fCez*f22;
1421   Double_t b20 = fCey*f00;
1422   Double_t b22 = fCey*f20 + fCce*f24 + fCee*f22;
1423   Double_t b30 = fCty*f00;
1424   Double_t b32 = fCty*f20 + fCct*f24 + fCte*f22;
1425   Double_t b40 = fCcy*f00;
1426   Double_t b42 = fCcy*f20 + fCcc*f24 + fCce*f22;
1427
1428   // a = f*b = f*C*ft
1429   Double_t a00 = f00*b00;
1430   Double_t a02 = f00*b02;
1431   Double_t a22 = f20*b02  + f24*b42  + f22*b22;
1432
1433   // F*C*Ft = C + (a + b + bt)
1434   fCyy += a00 + 2.0*b00;
1435   fCzy += b10;
1436   fCey += a02 +     b20 + b02;
1437   fCty += b30;
1438   fCcy += b40;
1439   fCez += b12;
1440   fCte += b32;
1441   fCee += a22 + 2.0*b22;
1442   fCce += b42;
1443
1444   return 1;
1445
1446 }
1447
1448 //_____________________________________________________________________________
1449 Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
1450 {
1451   //
1452   // Returns the track chi2
1453   //
1454
1455   Bool_t fNoTilt = kTRUE;
1456   if (TMath::Abs(h01) > 0.003) {
1457     fNoTilt = kFALSE;
1458   }
1459
1460   Double_t chi2;
1461   Double_t dy;
1462   Double_t r00;
1463   Double_t r01;
1464   Double_t r11;
1465
1466   if (fNoTilt) {
1467
1468     dy   = c->GetY() - fY;
1469     r00  = c->GetSigmaY2();
1470     chi2 = (dy*dy) / r00;
1471
1472   }
1473   else {
1474
1475     Double_t padlength = TMath::Sqrt(c->GetSigmaZ2() * 12.0);
1476
1477     r00  = c->GetSigmaY2();
1478     r01  = 0.0;
1479     r11  = c->GetSigmaZ2();
1480     r00 += fCyy;
1481     r01 += fCzy;
1482     r11 += fCzz;
1483
1484     Double_t det = r00*r11 - r01*r01;
1485     if (TMath::Abs(det) < 1.e-10) {
1486       Int_t n = GetNumberOfClusters();
1487       if (n > 4) {
1488         AliError(Form("Singular matrix N = %d!\n",n));
1489       }
1490       return 1e10;
1491     }
1492
1493     Double_t tmp = r00;
1494     r00  =  r11;
1495     r11  =  tmp;
1496     r01  = -r01;
1497     Double_t dy     = c->GetY() - fY;
1498     Double_t dz     = c->GetZ() - fZ;
1499     Double_t tiltdz = dz;
1500     if (TMath::Abs(tiltdz) > padlength/2.0) {
1501       tiltdz = TMath::Sign(padlength/2.0,dz);
1502     }
1503     dy = dy + h01*tiltdz;
1504
1505     chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz) / det;
1506
1507   }
1508
1509   return chi2;
1510
1511 }
1512
1513 //_________________________________________________________________________
1514 void AliTRDtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const
1515 {
1516   //
1517   // Returns reconstructed track momentum in the global system.
1518   //
1519
1520   Double_t pt = TMath::Abs(GetPt());
1521   Double_t r  = fC*fX - fE;
1522
1523   Double_t y0;
1524   if      (r >  1) {
1525     py =  pt;
1526     px = 0.0;
1527   }
1528   else if (r < -1) {
1529     py = -pt;
1530     px = 0.0;
1531   }
1532   else {
1533     y0 =  fY + TMath::Sqrt(1.0 - r*r) / fC;
1534     px = -pt * (fY - y0) * fC; //cos(phi);
1535     py = -pt * (fE - fX*fC);   //sin(phi);
1536   }
1537
1538   pz = pt*fT;
1539   Double_t tmp = px * TMath::Cos(fAlpha)
1540                - py * TMath::Sin(fAlpha);
1541   py = px * TMath::Sin(fAlpha)
1542      + py * TMath::Cos(fAlpha);
1543   px = tmp;
1544
1545 }
1546
1547 //_________________________________________________________________________
1548 void AliTRDtrack::GetGlobalXYZ(Double_t& x, Double_t& y, Double_t& z) const
1549 {
1550   //
1551   // Returns reconstructed track coordinates in the global system.
1552   //
1553
1554   x = fX;
1555   y = fY;
1556   z = fZ
1557 ;
1558   Double_t tmp = x * TMath::Cos(fAlpha)
1559                - y * TMath::Sin(fAlpha);
1560   y = x * TMath::Sin(fAlpha)
1561     + y * TMath::Cos(fAlpha);
1562   x = tmp;
1563
1564 }
1565
1566 //_________________________________________________________________________
1567 void AliTRDtrack::ResetCovariance()
1568 {
1569   //
1570   // Resets covariance matrix
1571   //
1572
1573   fCyy *= 10.0;
1574   fCzy  =  0.0;  fCzz *= 10.0;
1575   fCey  =  0.0;  fCez  =  0.0;  fCee *= 10.0;
1576   fCty  =  0.0;  fCtz  =  0.0;  fCte  =  0.0;  fCtt *= 10.0;
1577   fCcy  =  0.0;  fCcz  =  0.0;  fCce  =  0.0;  fCct  =  0.0;  fCcc *= 10.0;
1578
1579 }
1580
1581 //_____________________________________________________________________________
1582 void AliTRDtrack::ResetCovariance(Float_t mult)
1583 {
1584   //
1585   // Resets covariance matrix
1586   //
1587
1588   fCyy *= mult;
1589   fCzy *= 0.0;   fCzz *= 1.0;
1590   fCey *= 0.0;   fCez *= 0.0;   fCee *= mult;
1591   fCty *= 0.0;   fCtz *= 0.0;   fCte *= 0.0;   fCtt *= 1.0;
1592   fCcy *= 0.0;   fCcz *= 0.0;   fCce *= 0.0;   fCct *= 0.0;   fCcc *= mult;
1593
1594 }
1595
1596 //_____________________________________________________________________________
1597 void AliTRDtrack::MakeBackupTrack()
1598 {
1599   //
1600   // Creates a backup track
1601   //
1602
1603   if (fBackupTrack) {
1604     delete fBackupTrack;
1605   }
1606
1607   fBackupTrack = new AliTRDtrack(*this);
1608
1609 }
1610
1611 //_____________________________________________________________________________
1612 Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
1613 {
1614   //
1615   // Find prolongation at given x
1616   // Return 0 if it does not exist
1617   //
1618
1619   Double_t c1 = fC*fX - fE;
1620   if (TMath::Abs(c1) > 1.0) {
1621     return 0;
1622   }
1623
1624   Double_t r1 = TMath::Sqrt(1.0 - c1*c1);
1625   Double_t c2 = fC*xk - fE;
1626   if (TMath::Abs(c2) > 1.0) {
1627     return 0;
1628   }
1629
1630   Double_t r2 = TMath::Sqrt(1.0 - c2*c2);
1631   y = fY + (xk-fX)*(c1+c2)/(r1+r2);
1632   z = fZ + (xk-fX)*(c1+c2)/(c1*r2 + c2*r1)*fT;
1633
1634   return 1;
1635
1636 }
1637
1638 //_____________________________________________________________________________
1639 Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
1640 {
1641   //
1642   // Propagate track to a given x position
1643   // Works inside of the 20 degree segmentation
1644   // (local cooordinate frame for TRD , TPC, TOF)
1645   //
1646   // The material budget is taken from the geo manager
1647   //
1648
1649   Double_t  xyz0[3];
1650   Double_t  xyz1[3];
1651   Double_t  y;
1652   Double_t  z;
1653
1654   // Critical alpha  - cross sector indication
1655   const Double_t kAlphac  = TMath::Pi()/9.0;
1656   const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
1657
1658   // Direction +-
1659   Double_t dir = (fX > xr) ? -1.0 : 1.0;
1660
1661   for (Double_t x = fX + dir*step; dir*x < dir*xr; x += dir*step) {
1662
1663     GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);
1664     GetProlongation(x,y,z);
1665     xyz1[0] = x * TMath::Cos(fAlpha) + y * TMath::Sin(fAlpha);
1666     xyz1[1] = x * TMath::Sin(fAlpha) - y * TMath::Cos(fAlpha);
1667     xyz1[2] = z;
1668     Double_t param[7];
1669     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1670
1671     if ((param[0] > 0) &&
1672         (param[1] > 0)) {
1673       PropagateTo(x,param[1],param[0]);
1674     }
1675     if (fY >  fX*kTalphac) {
1676       Rotate(-kAlphac);
1677     }
1678     if (fY < -fX*kTalphac) {
1679       Rotate(kAlphac);
1680     }
1681
1682   }
1683
1684   PropagateTo(xr);
1685
1686   return 0;
1687
1688 }
1689
1690 //_____________________________________________________________________________
1691 Int_t AliTRDtrack::PropagateToR(Double_t r,Double_t step)
1692 {
1693   //
1694   // Propagate a track to a given radial position
1695   // The rotation is always connected to the last track position
1696   //
1697
1698   Double_t xyz0[3];
1699   Double_t xyz1[3];
1700   Double_t y;
1701   Double_t z;
1702
1703   // Direction +-
1704   Double_t radius = TMath::Sqrt(fX*fX + fY*fY);
1705   Double_t dir    = (radius > r) ? -1.0 : 1.0;
1706
1707   for (Double_t x = radius + dir*step; dir*x < dir*r; x += dir*step) {
1708     GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);
1709     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1710     Rotate(alpha,kTRUE);
1711     GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);
1712     GetProlongation(x,y,z);
1713     xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1714     xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1715     xyz1[2] = z;
1716     Double_t param[7];
1717     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1718     if (param[1] <= 0.0) {
1719       param[1] = 100000000.0;
1720     }
1721     PropagateTo(x,param[1],param[0]);
1722   }
1723
1724   GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);
1725   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
1726   Rotate(alpha,kTRUE);
1727   GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);
1728   GetProlongation(r,y,z);
1729   xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
1730   xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
1731   xyz1[2] = z;
1732   Double_t param[7];
1733   AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
1734
1735   if (param[1] <= 0.0) {
1736     param[1] = 100000000.0;
1737   }
1738   PropagateTo(r,param[1],param[0]);
1739
1740   return 0;
1741
1742 }
1743
1744 //_____________________________________________________________________________
1745 Int_t AliTRDtrack::GetSector() const
1746 {
1747   //
1748   // Return the current sector
1749   //
1750
1751   return Int_t(TVector2::Phi_0_2pi(fAlpha)
1752              / AliTRDgeometry::GetAlpha())
1753              % AliTRDgeometry::kNsect;
1754
1755 }
1756
1757 //_____________________________________________________________________________
1758 Double_t AliTRDtrack::Get1Pt() const
1759 {
1760   //
1761   // Returns the inverse Pt (1/GeV/c)
1762   // (or 1/"most probable pt", if the field is too weak)
1763   //
1764
1765   if (TMath::Abs(GetLocalConvConst()) > kVeryBigConvConst) {
1766     return 1.0 / kMostProbableMomentum
1767                / TMath::Sqrt(1.0 + GetTgl()*GetTgl());
1768   }
1769
1770   return (TMath::Sign(1.0e-9,fC) + fC) * GetLocalConvConst();
1771
1772 }
1773
1774 //_____________________________________________________________________________
1775 Double_t AliTRDtrack::GetP() const
1776 {
1777   //
1778   // Returns the total momentum
1779   //
1780
1781   return TMath::Abs(GetPt()) * TMath::Sqrt(1.0 + GetTgl()*GetTgl());
1782
1783 }
1784
1785 //_____________________________________________________________________________
1786 Double_t AliTRDtrack::GetYat(Double_t xk) const
1787 {
1788   //
1789   // This function calculates the Y-coordinate of a track at
1790   // the plane x = xk.
1791   // Needed for matching with the TOF (I.Belikov)
1792   //
1793
1794   Double_t c1 = fC*fX - fE;
1795   Double_t r1 = TMath::Sqrt(1.0 - c1*c1);
1796   Double_t c2 = fC*xk - fE;
1797   Double_t r2 = TMath::Sqrt(1.0-  c2*c2);
1798
1799   return fY + (xk-fX)*(c1+c2)/(r1+r2);
1800
1801 }
1802
1803 //_____________________________________________________________________________
1804 void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
1805 {
1806   //
1807   // The sampled energy loss
1808   //
1809
1810   Double_t s = GetSnp();
1811   Double_t t = GetTgl();
1812
1813   q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1814   fdQdl[i] = q;
1815
1816 }
1817
1818  //_____________________________________________________________________________
1819 void AliTRDtrack::SetSampledEdx(Float_t q)
1820 {
1821   //
1822   // The sampled energy loss
1823   //
1824
1825   Double_t s = GetSnp();
1826   Double_t t = GetTgl();
1827
1828   q*= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1829   fdQdl[fNdedx] = q;
1830   fNdedx++;
1831
1832 }
1833
1834 //_____________________________________________________________________________
1835 void AliTRDtrack::GetXYZ(Float_t r[3]) const
1836 {
1837   //
1838   // Returns the position of the track in the global coord. system
1839   //
1840
1841   Double_t cs = TMath::Cos(fAlpha);
1842   Double_t sn = TMath::Sin(fAlpha);
1843   r[0] = fX*cs - fY*sn;
1844   r[1] = fX*sn + fY*cs;
1845   r[2] = fZ;
1846
1847 }