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   if (GetLabel()==277)
 676     printf("Rotate %e %e %e %e\n",GetAlpha(),GetX(),GetY(),GetZ());
 677
 678   return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
 679 }
 680
 681 //_____________________________________________________________________________
 682 Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
 683 {
 684   //
 685   // Returns the track chi2
 686   //
 687
 688   Double_t p[2]={c->GetY(), c->GetZ()};
 689   Double_t sy2=c->GetSigmaY2()*4;
 690   Double_t sz2=c->GetSigmaZ2()*4;
 691   Double_t cov[3]={sy2 + h01*h01*sz2, h01*(sy2-sz2), sz2 + h01*h01*sy2};
 692
 693   return AliExternalTrackParam::GetPredictedChi2(p,cov);
 694
 695   /*
 696   Bool_t fNoTilt = kTRUE;
 697   if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE;
 698
 699   return (c->GetY() - GetY())*(c->GetY() - GetY())/c->GetSigmaY2();
 700   */
 701
 702   /*
 703   Double_t chi2, dy, r00, r01, r11;
 704
 705   if(fNoTilt) {
 706     dy=c->GetY() - fY;
 707     r00=c->GetSigmaY2();
 708     chi2 = (dy*dy)/r00;
 709   }
 710   else {
 711     Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12);
 712     //
 713     r00=c->GetSigmaY2(); r01=0.; r11=c->GetSigmaZ2();
 714     r00+=fCyy; r01+=fCzy; r11+=fCzz;
 715
 716     Double_t det=r00*r11 - r01*r01;
 717     if (TMath::Abs(det) < 1.e-10) {
 718       Int_t n=GetNumberOfClusters();
 719       if (n>4) cerr<<n<<" AliTRDtrack warning: Singular matrix !\n";
 720       return 1e10;
 721     }
 722     Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
 723     Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
 724     Double_t tiltdz = dz;
 725     if (TMath::Abs(tiltdz)>padlength/2.) {
 726       tiltdz = TMath::Sign(padlength/2,dz);
 727     }
 728     //    dy=dy+h01*dz;
 729     dy=dy+h01*tiltdz;
 730
 731     chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
 732   }
 733
 734   return chi2;
 735   */
 736 }
 737
 738 //_____________________________________________________________________________
 739 void AliTRDtrack::MakeBackupTrack()
 740 {
 741   //
 742   // Creates a backup track
 743   //
 744
 745   if (fBackupTrack) delete fBackupTrack;
 746   fBackupTrack = new AliTRDtrack(*this);
 747
 748 }
 749
 750 //_____________________________________________________________________________
 751 Int_t AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
 752 {
 753   //
 754   // Find prolongation at given x
 755   // return 0 if not exist
 756
 757   Double_t bz=GetBz();
 758
 759   if (!AliExternalTrackParam::GetYAt(xk,bz,y)) return 0;
 760   if (!AliExternalTrackParam::GetZAt(xk,bz,z)) return 0;
 761
 762   return 1;
 763 }
 764
 765 //_____________________________________________________________________________
 766 Int_t   AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
 767 {
 768   //
 769   // Propagate track to given x  position
 770   // works inside of the 20 degree segmentation (local cooordinate frame for TRD , TPC, TOF)
 771   //
 772   // material budget from geo manager
 773   //
 774   Double_t  xyz0[3], xyz1[3],y,z;
 775   const Double_t kAlphac  = TMath::Pi()/9.;
 776   const Double_t kTalphac = TMath::Tan(kAlphac*0.5);
 777   // critical alpha  - cross sector indication
 778   //
 779   Double_t dir = (GetX()>xr) ? -1.:1.;
 780   // direction +-
 781   for (Double_t x=GetX()+dir*step;dir*x<dir*xr;x+=dir*step){
 782     //
 783     GetXYZ(xyz0);
 784     GetProlongation(x,y,z);
 785     xyz1[0] = x*TMath::Cos(GetAlpha())+y*TMath::Sin(GetAlpha());
 786     xyz1[1] = x*TMath::Sin(GetAlpha())-y*TMath::Cos(GetAlpha());
 787     xyz1[2] = z;
 788     Double_t param[7];
 789     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
 790     //
 791     if (param[0]>0&&param[1]>0) PropagateTo(x,param[1],param[0]);
 792     if (GetY()>GetX()*kTalphac){
 793       Rotate(-kAlphac);
 794     }
 795     if (GetY()<-GetX()*kTalphac){
 796       Rotate(kAlphac);
 797     }
 798   }
 799   //
 800   PropagateTo(xr);
 801
 802   return 0;
 803
 804 }
 805
 806 //_____________________________________________________________________________
 807 Int_t   AliTRDtrack::PropagateToR(Double_t r,Double_t step)
 808 {
 809   //
 810   // propagate track to the radial position
 811   // rotation always connected to the last track position
 812   //
 813   Double_t  xyz0[3], xyz1[3],y,z;
 814   Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
 815   Double_t dir = (radius>r) ? -1.:1.;   // direction +-
 816   //
 817   for (Double_t x=radius+dir*step;dir*x<dir*r;x+=dir*step){
 818     GetXYZ(xyz0);
 819     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 820     Rotate(alpha,kTRUE);
 821     GetXYZ(xyz0);
 822     GetProlongation(x,y,z);
 823     xyz1[0] = x*TMath::Cos(alpha)+y*TMath::Sin(alpha);
 824     xyz1[1] = x*TMath::Sin(alpha)-y*TMath::Cos(alpha);
 825     xyz1[2] = z;
 826     Double_t param[7];
 827     AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
 828     if (param[1]<=0) param[1] =100000000;
 829     PropagateTo(x,param[1],param[0]);
 830   }
 831   GetXYZ(xyz0);
 832   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 833   Rotate(alpha,kTRUE);
 834   GetXYZ(xyz0);
 835   GetProlongation(r,y,z);
 836   xyz1[0] = r*TMath::Cos(alpha)+y*TMath::Sin(alpha);
 837   xyz1[1] = r*TMath::Sin(alpha)-y*TMath::Cos(alpha);
 838   xyz1[2] = z;
 839   Double_t param[7];
 840   AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
 841   //
 842   if (param[1]<=0) param[1] =100000000;
 843   PropagateTo(r,param[1],param[0]);
 844
 845   return 0;
 846
 847 }
 848
 849 //_____________________________________________________________________________
 850 Int_t AliTRDtrack::GetSector() const
 851 {
 852   //
 853   // Return the current sector
 854   //
 855
 856   return Int_t(TVector2::Phi_0_2pi(GetAlpha())
 857              / AliTRDgeometry::GetAlpha())
 858              % AliTRDgeometry::kNsect;
 859
 860 }
 861
 862 //_____________________________________________________________________________
 863 void AliTRDtrack::SetSampledEdx(Float_t q, Int_t i)
 864 {
 865   //
 866   // The sampled energy loss
 867   //
 868
 869   Double_t s = GetSnp();
 870   Double_t t = GetTgl();
 871   q *= TMath::Sqrt((1-s*s)/(1+t*t));
 872   fdQdl[i] = q;
 873
 874 }
 875
 876  //_____________________________________________________________________________
 877 void AliTRDtrack::SetSampledEdx(Float_t q)
 878 {
 879   //
 880   // The sampled energy loss
 881   //
 882
 883   Double_t s = GetSnp();
 884   Double_t t = GetTgl();
 885   q*= TMath::Sqrt((1-s*s)/(1+t*t));
 886   fdQdl[fNdedx] = q;
 887   fNdedx++;
 888
 889 }
 890
 891 Double_t AliTRDtrack::GetBz() const {
 892   //
 893   // returns Bz component of the magnetic field (kG)
 894   //
 895   if (AliTracker::UniformField()) return AliTracker::GetBz();
 896   Double_t r[3]; GetXYZ(r);
 897   return AliTracker::GetBz(r);
 898 }
 899
 900