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