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   // Array to sort the dEdx values according to amplitude
 403   Float_t sorted[kMAXCLUSTERSPERTRACK];
 404   fdEdx = 0.0;
 405
 406   // Require at least 10 clusters for a dedx measurement
 407   if (fNdedx < 10) return;
 408
 409   // Can fdQdl be negative ????
 410   for (Int_t i = 0; i < fNdedx; i++) {
 411     sorted[i] = TMath::Abs(fdQdl[i]);
 412   }
 413   // Sort the dedx values by amplitude
 414   Int_t *index = new Int_t[fNdedx];
 415   TMath::Sort(fNdedx, sorted, index, kFALSE);
 416
 417   // Sum up the truncated charge between lower and upper bounds
 418   Int_t nl = Int_t(low * fNdedx);
 419   Int_t nu = Int_t( up * fNdedx);
 420   for (Int_t i = nl; i <= nu; i++) {
 421     fdEdx += sorted[index[i]];
 422   }
 423   fdEdx /= (nu - nl + 1.0);
 424
 425   delete[] index;
 426
 427 }
 428
 429 //_____________________________________________________________________________
 430 Bool_t AliTRDtrack::PropagateTo(Double_t xk, Double_t x0, Double_t rho)
 431 {
 432   //
 433   // Propagates a track of particle with mass=pm to a reference plane
 434   // defined by x=xk through media of density=rho and radiationLength=x0
 435   //
 436
 437   if (xk == GetX()) {
 438     return kTRUE;
 439   }
 440
 441   Double_t oldX = GetX();
 442   Double_t oldY = GetY();
 443   Double_t oldZ = GetZ();
 444
 445   Double_t bz   = GetBz();
 446
 447   if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
 448     return kFALSE;
 449   }
 450
 451   Double_t x = GetX();
 452   Double_t y = GetY();
 453   Double_t z = GetZ();
 454   Double_t d = TMath::Sqrt((x-oldX)*(x-oldX)
 455                          + (y-oldY)*(y-oldY)
 456                          + (z-oldZ)*(z-oldZ));
 457   if (oldX < xk) {
 458     if (IsStartedTimeIntegral()) {
 459       Double_t l2  = d;
 460       Double_t crv = GetC();
 461       if (TMath::Abs(l2*crv) > 0.0001) {
 462         // Make correction for curvature if neccesary
 463         l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX) + (y-oldY)*(y-oldY));
 464         l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
 465         l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ));
 466       }
 467       AddTimeStep(l2);
 468     }
 469   }
 470
 471   Double_t ll = (oldX < xk) ? -d : d;
 472   if (!AliExternalTrackParam::CorrectForMaterial(ll*rho/x0,x0,GetMass())) {
 473     return kFALSE;
 474   }
 475
 476   {
 477
 478     // Energy losses************************
 479     Double_t p2    = (1.0 + GetTgl()*GetTgl()) / (Get1Pt()*Get1Pt());
 480     Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
 481     if (beta2<1.0e-10 || (5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0) {
 482       return kFALSE;
 483     }
 484
 485     Double_t dE    = 0.153e-3 / beta2
 486                    * (log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
 487                    * d * rho;
 488     Float_t budget = d * rho;
 489     fBudget[0] += budget;
 490
 491     /*
 492     // Suspicious part - think about it ?
 493     Double_t kinE =  TMath::Sqrt(p2);
 494     if (dE > 0.8*kinE) dE = 0.8 * kinE;  //
 495     if (dE < 0)        dE = 0.0;         // Not valid region for Bethe bloch
 496     */
 497
 498     fDE += dE;
 499
 500     /*
 501     // Suspicious ! I.B.
 502     Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE));   // Energy loss fluctuation
 503     Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
 504     fCcc += sigmac2;
 505     fCee += fX*fX * sigmac2;
 506     */
 507
 508   }
 509
 510   return kTRUE;
 511
 512 }
 513
 514 //_____________________________________________________________________________
 515 Bool_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, Int_t index
 516                          , Double_t h01)
 517 {
 518   //
 519   // Assignes found cluster to the track and updates track information
 520   //
 521
 522   Bool_t fNoTilt = kTRUE;
 523   if (TMath::Abs(h01) > 0.003) {
 524     fNoTilt = kFALSE;
 525   }
 526
 527   // Add angular effect to the error contribution -  MI
 528   Float_t tangent2 = GetSnp()*GetSnp();
 529   if (tangent2 < 0.90000) {
 530     tangent2 = tangent2 / (1.0 - tangent2);
 531   }
 532   //Float_t errang = tangent2 * 0.04;
 533
 534   Double_t p[2]   = {c->GetY(), c->GetZ() };
 535   //Double_t cov[3] = {c->GetSigmaY2()+errang, 0.0, c->GetSigmaZ2()*100.0 };
 536   Double_t sy2    = c->GetSigmaY2() * 4.0;
 537   Double_t sz2    = c->GetSigmaZ2() * 4.0;
 538   Double_t cov[3] = {sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
 539
 540   if (!AliExternalTrackParam::Update(p,cov)) {
 541     return kFALSE;
 542   }
 543
 544   Int_t n = GetNumberOfClusters();
 545   fIndex[n] = index;
 546   SetNumberOfClusters(n+1);
 547
 548   SetChi2(GetChi2()+chisq);
 549
 550   return kTRUE;
 551
 552 }
 553
 554 //_____________________________________________________________________________
 555 Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, Int_t index
 556                           , Double_t h01, Int_t /*plane*/)
 557 {
 558   //
 559   // Assignes found cluster to the track and updates track information
 560   //
 561
 562   Bool_t fNoTilt = kTRUE;
 563   if (TMath::Abs(h01) > 0.003) {
 564     fNoTilt = kFALSE;
 565   }
 566
 567   // Add angular effect to the error contribution and make correction  -  MI
 568   Double_t tangent2 = GetSnp()*GetSnp();
 569   if (tangent2 < 0.90000){
 570     tangent2 = tangent2 / (1.0-tangent2);
 571   }
 572   Double_t tangent = TMath::Sqrt(tangent2);
 573   if (GetSnp() < 0) {
 574     tangent *= -1;
 575   }
 576
 577   //
 578   // Is the following still needed ????
 579   //
 580
 581   //  Double_t correction = 0*plane;
 582   /*
 583   Double_t errang = tangent2*0.04;  //
 584   Double_t errsys =0.025*0.025*20;  //systematic error part
 585
 586   Float_t extend =1;
 587   if (c->GetNPads()==4) extend=2;
 588   */
 589   //if (c->GetNPads()==5)  extend=3;
 590   //if (c->GetNPads()==6)  extend=3;
 591   //if (c->GetQ()<15) return 1;
 592
 593   /*
 594   if (corrector!=0){
 595   //if (0){
 596     correction = corrector->GetCorrection(plane,c->GetLocalTimeBin(),tangent);
 597     if (TMath::Abs(correction)>0){
 598       //if we have info
 599       errang     = corrector->GetSigma(plane,c->GetLocalTimeBin(),tangent);
 600       errang    *= errang;
 601       errang    += tangent2*0.04;
 602     }
 603   }
 604   */
 605   //
 606   //Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
 607   /*
 608     {
 609       Double_t dy=c->GetY() - GetY(), dz=c->GetZ() - GetZ();
 610       printf("%e %e %e %e\n",dy,dz,padlength/2,h01);
 611     }
 612   */
 613
 614   Double_t p[2]   = { c->GetY(), c->GetZ() };
 615   //Double_t cov[3]={ (c->GetSigmaY2()+errang+errsys)*extend, 0.0, c->GetSigmaZ2()*10000.0 };
 616   Double_t sy2    = c->GetSigmaY2() * 4.0;
 617   Double_t sz2    = c->GetSigmaZ2() * 4.0;
 618   Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
 619
 620   if (!AliExternalTrackParam::Update(p,cov)) {
 621     return kFALSE;
 622   }
 623
 624   Int_t n = GetNumberOfClusters();
 625   fIndex[n] = index;
 626   SetNumberOfClusters(n+1);
 627   SetChi2(GetChi2() + chisq);
 628
 629   return kTRUE;
 630
 631 }
 632
 633 // //_____________________________________________________________________________
 634 // Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet)
 635 // {
 636 //   //
 637 //   // Assignes found tracklet to the track and updates track information
 638 //   //
 639 //   // Can this be removed ????
 640 //   //
 641 //
 642 //   Double_t r00=(tracklet.GetTrackletSigma2()), r01=0., r11= 10000.;
 643 //   r00+=fCyy; r01+=fCzy; r11+=fCzz;
 644 //   //
 645 //   Double_t det=r00*r11 - r01*r01;
 646 //   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 647 //   //
 648
 649 //   Double_t dy=tracklet.GetY() - fY, dz=tracklet.GetZ() - fZ;
 650
 651
 652 //   Double_t s00 = tracklet.GetTrackletSigma2();  // error pad
 653 //   Double_t s11 = 100000;   // error pad-row
 654 //   Float_t  h01 = tracklet.GetTilt();
 655 //   //
 656 //   //  r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00;
 657 //   r00 = fCyy + fCzz*h01*h01+s00;
 658 //   //  r01 = fCzy + fCzz*h01;
 659 //   r01 = fCzy ;
 660 //   r11 = fCzz + s11;
 661 //   det = r00*r11 - r01*r01;
 662 //   // inverse matrix
 663 //   tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 664
 665 //   Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11;
 666 //   Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11;
 667 //   Double_t k20=fCey*r00+fCez*r01, k21=fCey*r01+fCez*r11;
 668 //   Double_t k30=fCty*r00+fCtz*r01, k31=fCty*r01+fCtz*r11;
 669 //   Double_t k40=fCcy*r00+fCcz*r01, k41=fCcy*r01+fCcz*r11;
 670
 671 //   // K matrix
 672 // //   k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01);
 673 // //   k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01);
 674 // //   k20=fCey*r00+fCez*(r01+h01*r00),k21=fCey*r01+fCez*(r11+h01*r01);
 675 // //   k30=fCty*r00+fCtz*(r01+h01*r00),k31=fCty*r01+fCtz*(r11+h01*r01);
 676 // //   k40=fCcy*r00+fCcz*(r01+h01*r00),k41=fCcy*r01+fCcz*(r11+h01*r01);
 677 //   //
 678 //   //Update measurement
 679 //   Double_t cur=fC + k40*dy + k41*dz, eta=fE + k20*dy + k21*dz;
 680 //   //  cur=fC + k40*dy + k41*dz; eta=fE + k20*dy + k21*dz;
 681 //   if (TMath::Abs(cur*fX-eta) >= 0.90000) {
 682 //     //Int_t n=GetNumberOfClusters();
 683 //     //      if (n>4) cerr<<n<<" AliTRDtrack warning: Filtering failed !\n";
 684 //     return 0;
 685 //   }
 686 // //   k01+=h01*k00;
 687 // //   k11+=h01*k10;
 688 // //   k21+=h01*k20;
 689 // //   k31+=h01*k30;
 690 // //   k41+=h01*k40;
 691
 692
 693 //   fY += k00*dy + k01*dz;
 694 //   fZ += k10*dy + k11*dz;
 695 //   fE  = eta;
 696 //   fT += k30*dy + k31*dz;
 697 //   fC  = cur;
 698
 699
 700 //   //Update covariance
 701 //   //
 702 //   //
 703 //   Double_t oldyy = fCyy, oldzz = fCzz; //, oldee=fCee, oldcc =fCcc;
 704 //   Double_t oldzy = fCzy, oldey = fCey, oldty=fCty, oldcy =fCcy;
 705 //   Double_t oldez = fCez, oldtz = fCtz, oldcz=fCcz;
 706 //   //Double_t oldte = fCte, oldce = fCce;
 707 //   //Double_t oldct = fCct;
 708
 709 //   fCyy-=k00*oldyy+k01*oldzy;
 710 //   fCzy-=k10*oldyy+k11*oldzy;
 711 //   fCey-=k20*oldyy+k21*oldzy;
 712 //   fCty-=k30*oldyy+k31*oldzy;
 713 //   fCcy-=k40*oldyy+k41*oldzy;
 714 //   //
 715 //   fCzz-=k10*oldzy+k11*oldzz;
 716 //   fCez-=k20*oldzy+k21*oldzz;
 717 //   fCtz-=k30*oldzy+k31*oldzz;
 718 //   fCcz-=k40*oldzy+k41*oldzz;
 719 //   //
 720 //   fCee-=k20*oldey+k21*oldez;
 721 //   fCte-=k30*oldey+k31*oldez;
 722 //   fCce-=k40*oldey+k41*oldez;
 723 //   //
 724 //   fCtt-=k30*oldty+k31*oldtz;
 725 //   fCct-=k40*oldty+k41*oldtz;
 726 //   //
 727 //   fCcc-=k40*oldcy+k41*oldcz;
 728 //   //
 729
 730 //   //Int_t n=GetNumberOfClusters();
 731 //   //fIndex[n]=index;
 732 //   //SetNumberOfClusters(n+1);
 733
 734 //   //SetChi2(GetChi2()+chisq);
 735 //   //  cerr<<"in update: fIndex["<<fN<<"] = "<<index<<endl;
 736
 737 //   return 1;
 738
 739 // }
 740
 741 //_____________________________________________________________________________
 742 Bool_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
 743 {
 744   //
 745   // Rotates track parameters in R*phi plane
 746   // if absolute rotation alpha is in global system
 747   // otherwise alpha rotation is relative to the current rotation angle
 748   //
 749
 750   if (absolute) {
 751     alpha -= GetAlpha();
 752   }
 753   else{
 754     fNRotate++;
 755   }
 756
 757   return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
 758
 759 }
 760
 761 //_____________________________________________________________________________
 762 Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
 763 {
 764   //
 765   // Returns the track chi2
 766   //
 767
 768   Double_t p[2]   = { c->GetY(), c->GetZ() };
 769   Double_t sy2    = c->GetSigmaY2() * 4.0;
 770   Double_t sz2    = c->GetSigmaZ2() * 4.0;
 771   Double_t cov[3] = { sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2 };
 772
 773   return AliExternalTrackParam::GetPredictedChi2(p,cov);
 774
 775   //
 776   // Can the following be removed ????
 777   //
 778   /*
 779   Bool_t fNoTilt = kTRUE;
 780   if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE;
 781
 782   return (c->GetY() - GetY())*(c->GetY() - GetY())/c->GetSigmaY2();
 783   */
 784
 785   /*
 786   Double_t chi2, dy, r00, r01, r11;
 787
 788   if(fNoTilt) {
 789     dy=c->GetY() - fY;
 790     r00=c->GetSigmaY2();
 791     chi2 = (dy*dy)/r00;
 792   }
 793   else {
 794     Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12);
 795     //
 796     r00=c->GetSigmaY2(); r01=0.; r11=c->GetSigmaZ2();
 797     r00+=fCyy; r01+=fCzy; r11+=fCzz;
 798
 799     Double_t det=r00*r11 - r01*r01;
 800     if (TMath::Abs(det) < 1.e-10) {
 801       Int_t n=GetNumberOfClusters();
 802       if (n>4) cerr<<n<<" AliTRDtrack warning: Singular matrix !\n";
 803       return 1e10;
 804     }
 805     Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
 806     Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
 807     Double_t tiltdz = dz;
 808     if (TMath::Abs(tiltdz)>padlength/2.) {
 809       tiltdz = TMath::Sign(padlength/2,dz);
 810     }
 811     //    dy=dy+h01*dz;
 812     dy=dy+h01*tiltdz;
 813
 814     chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
 815   }
 816
 817   return chi2;
 818   */
 819
 820 }
 821
 822 //_____________________________________________________________________________
 823 void AliTRDtrack::MakeBackupTrack()
 824 {
 825   //
 826   // Creates a backup track
 827   //
 828
 829   if (fBackupTrack) {
 830     delete fBackupTrack;
 831   }
 832   fBackupTrack = new AliTRDtrack(*this);
 833
 834 }
 835
 836 //_____________________________________________________________________________
 837 Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
 838 {
 839   //
 840   // Find a prolongation at given x
 841   // Return 0 if it does not exist
 842   //
 843
 844   Double_t bz = GetBz();
 845
 846   if (!AliExternalTrackParam::GetYAt(xk,bz,y)) {
 847     return 0;
 848   }
 849   if (!AliExternalTrackParam::GetZAt(xk,bz,z)) {
 850     return 0;
 851   }
 852
 853   return 1;
 854
 855 }
 856
 857 //_____________________________________________________________________________
 858 Int_t AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
 859 {
 860   //
 861   // Propagate track to given x  position
 862   // Works inside of the 20 degree segmentation (local cooordinate frame for TRD , TPC, TOF)
 863   //
 864   // Material budget from geo manager
 865   //
 866
 867   Double_t xyz0[3];
 868   Double_t xyz1[3];
 869   Double_t y;
 870   Double_t z;
 871
 872   const Double_t kAlphac  = TMath::Pi()/9.0;
 873   const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
 874
 875   // Critical alpha  - cross sector indication
 876   Double_t dir = (GetX()>xr) ? -1.0 : 1.0;
 877
 878   // Direction +-
 879   for (Double_t x = GetX()+dir*step; dir*x < dir*xr; x += dir*step) {
 880
 881     GetXYZ(xyz0);
 882     GetProlongation(x,y,z);
 883     xyz1[0] = x * TMath::Cos(GetAlpha()) + y * TMath::Sin(GetAlpha());
 884     xyz1[1] = x * TMath::Sin(GetAlpha()) - y * TMath::Cos(GetAlpha());
 885     xyz1[2] = z;
 886     Double_t param[7];
 887     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
 888
 889     if ((param[0] > 0) &&
 890         (param[1] > 0)) {
 891       PropagateTo(x,param[1],param[0]);
 892     }
 893
 894     if (GetY() >  GetX()*kTalphac) {
 895       Rotate(-kAlphac);
 896     }
 897     if (GetY() < -GetX()*kTalphac) {
 898       Rotate( kAlphac);
 899     }
 900
 901   }
 902
 903   PropagateTo(xr);
 904
 905   return 0;
 906
 907 }
 908
 909 //_____________________________________________________________________________
 910 Int_t   AliTRDtrack::PropagateToR(Double_t r,Double_t step)
 911 {
 912   //
 913   // Propagate track to the radial position
 914   // Rotation always connected to the last track position
 915   //
 916
 917   Double_t xyz0[3];
 918   Double_t xyz1[3];
 919   Double_t y;
 920   Double_t z;
 921
 922   Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
 923   // Direction +-
 924   Double_t dir    = (radius>r) ? -1.0 : 1.0;
 925
 926   for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
 927
 928     GetXYZ(xyz0);
 929     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 930     Rotate(alpha,kTRUE);
 931     GetXYZ(xyz0);
 932     GetProlongation(x,y,z);
 933     xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
 934     xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
 935     xyz1[2] = z;
 936     Double_t param[7];
 937     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
 938     if (param[1] <= 0) {
 939       param[1] = 100000000;
 940     }
 941     PropagateTo(x,param[1],param[0]);
 942
 943   }
 944
 945   GetXYZ(xyz0);
 946   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 947   Rotate(alpha,kTRUE);
 948   GetXYZ(xyz0);
 949   GetProlongation(r,y,z);
 950   xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
 951   xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
 952   xyz1[2] = z;
 953   Double_t param[7];
 954   AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
 955
 956   if (param[1] <= 0) {
 957     param[1] = 100000000;
 958   }
 959   PropagateTo(r,param[1],param[0]);
 960
 961   return 0;
 962
 963 }
 964
 965 //_____________________________________________________________________________
 966 Int_t AliTRDtrack::GetSector() const
 967 {
 968   //
 969   // Return the current sector
 970   //
 971
 972   return Int_t(TVector2::Phi_0_2pi(GetAlpha()) / AliTRDgeometry::GetAlpha())
 973          % AliTRDgeometry::kNsect;
 974
 975 }
 976
 977 //_____________________________________________________________________________
 978 void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
 979 {
 980   //
 981   // The sampled energy loss
 982   //
 983
 984   Double_t s = GetSnp();
 985   Double_t t = GetTgl();
 986   q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
 987   fdQdl[i] = q;
 988
 989 }
 990
 991 //_____________________________________________________________________________
 992 void AliTRDtrack::SetSampledEdx(Float_t q)
 993 {
 994   //
 995   // The sampled energy loss
 996   //
 997
 998   Double_t s = GetSnp();
 999   Double_t t = GetTgl();
1000   q *= TMath::Sqrt((1.0 - s*s) / (1.0 + t*t));
1001   fdQdl[fNdedx] = q;
1002   fNdedx++;
1003
1004 }
1005
1006 //_____________________________________________________________________________
1007 Double_t AliTRDtrack::GetBz() const
1008 {
1009   //
1010   // Returns Bz component of the magnetic field (kG)
1011   //
1012
1013   if (AliTracker::UniformField()) {
1014     return AliTracker::GetBz();
1015   }
1016   Double_t r[3];
1017   GetXYZ(r);
1018
1019   return AliTracker::GetBz(r);
1020
1021 }