TRD/AliTRDtrackV1.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 "AliLog.h"
  19 #include "AliESDtrack.h"
  20 #include "AliTracker.h"
  21
  22 #include "AliTRDtrackV1.h"
  23 #include "AliTRDcluster.h"
  24 #include "AliTRDcalibDB.h"
  25 #include "AliTRDReconstructor.h"
  26 #include "AliTRDrecoParam.h"
  27
  28 ClassImp(AliTRDtrackV1)
  29
  30 ///////////////////////////////////////////////////////////////////////////////
  31 //                                                                           //
  32 //  Represents a reconstructed TRD track                                     //
  33 //  Local TRD Kalman track                                                   //
  34 //                                                                           //
  35 //  Authors:                                                                 //
  36 //    Alex Bercuci <A.Bercuci@gsi.de>                                        //
  37 //    Markus Fasel <M.Fasel@gsi.de>                                          //
  38 //                                                                           //
  39 ///////////////////////////////////////////////////////////////////////////////
  40
  41 //_______________________________________________________________
  42 AliTRDtrackV1::AliTRDtrackV1() : AliKalmanTrack()
  43   ,fPIDquality(0)
  44   ,fDE(0.)
  45   ,fReconstructor(0x0)
  46   ,fBackupTrack(0x0)
  47 {
  48   //
  49   // Default constructor
  50   //
  51   //printf("AliTRDtrackV1::AliTRDtrackV1()\n");
  52
  53   for(int i =0; i<3; i++) fBudget[i] = 0.;
  54
  55   Float_t pid = 1./AliPID::kSPECIES;
  56   for(int is =0; is<AliPID::kSPECIES; is++) fPID[is] = pid;
  57
  58   for(int ip=0; ip<kNplane; ip++){
  59     fTrackletIndex[ip] = 0xffff;
  60     fTracklet[ip]      = 0x0;
  61   }
  62 }
  63
  64 //_______________________________________________________________
  65 AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref) : AliKalmanTrack(ref)
  66   ,fPIDquality(ref.fPIDquality)
  67   ,fDE(ref.fDE)
  68   ,fReconstructor(ref.fReconstructor)
  69   ,fBackupTrack(0x0)
  70 {
  71   //
  72   // Copy constructor
  73   //
  74
  75   //printf("AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &)\n");
  76   SetBit(kOwner, kFALSE);
  77   for(int ip=0; ip<kNplane; ip++){
  78     fTrackletIndex[ip] = ref.fTrackletIndex[ip];
  79     fTracklet[ip]      = ref.fTracklet[ip];
  80   }
  81
  82   for (Int_t i = 0; i < 3;i++) fBudget[i]      = ref.fBudget[i];
  83
  84         for(Int_t is = 0; is<AliPID::kSPECIES; is++) fPID[is] = ref.fPID[is];
  85
  86   AliKalmanTrack::SetNumberOfClusters(ref.GetNumberOfClusters());
  87 }
  88
  89 //_______________________________________________________________
  90 AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) : AliKalmanTrack()
  91   ,fPIDquality(0)
  92   ,fDE(0.)
  93   ,fReconstructor(0x0)
  94   ,fBackupTrack(0x0)
  95 {
  96   //
  97   // Constructor from AliESDtrack
  98   //
  99
 100   SetLabel(t.GetLabel());
 101   SetChi2(0.0);
 102   SetMass(t.GetMass());
 103   AliKalmanTrack::SetNumberOfClusters(t.GetTRDncls());
 104   Int_t ti[kNplane]; t.GetTRDtracklets(&ti[0]);
 105   for(int ip=0; ip<kNplane; ip++){
 106     fTrackletIndex[ip] = ti[ip] < 0 ? 0xffff : ti[ip];
 107     fTracklet[ip]      = 0x0;
 108   }
 109   for(int i =0; i<3; i++) fBudget[i] = 0.;
 110
 111   Float_t pid = 1./AliPID::kSPECIES;
 112   for(int is =0; is<AliPID::kSPECIES; is++) fPID[is] = pid;
 113
 114   const AliExternalTrackParam *par = &t;
 115   if (t.GetStatus() & AliESDtrack::kTRDbackup) {
 116     par = t.GetOuterParam();
 117     if (!par) {
 118       AliError("No backup info!");
 119       par = &t;
 120     }
 121   }
 122   Set(par->GetX()
 123      ,par->GetAlpha()
 124      ,par->GetParameter()
 125      ,par->GetCovariance());
 126
 127   if(t.GetStatus() & AliESDtrack::kTIME) {
 128     StartTimeIntegral();
 129     Double_t times[10];
 130     t.GetIntegratedTimes(times);
 131     SetIntegratedTimes(times);
 132     SetIntegratedLength(t.GetIntegratedLength());
 133   }
 134 }
 135
 136 //_______________________________________________________________
 137 AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Double_t cov[15]
 138              , Double_t x, Double_t alpha) : AliKalmanTrack()
 139   ,fPIDquality(0)
 140   ,fDE(0.)
 141   ,fReconstructor(0x0)
 142   ,fBackupTrack(0x0)
 143 {
 144   //
 145   // The stand alone tracking constructor
 146   // TEMPORARY !!!!!!!!!!!
 147   // to check :
 148   // 1. covariance matrix
 149   // 2. dQdl calculation
 150   //
 151
 152   Double_t cnv = GetBz() < 1.e-5 ? 1.e5 : 1.0 / (GetBz() * kB2C);
 153   //  Double_t cnv   = 1.0 / (GetBz() * kB2C);
 154
 155   Double_t pp[5] = { p[0]
 156                    , p[1]
 157                    , p[2]
 158                    , p[3]
 159                    , p[4]*cnv      };
 160
 161   Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[ 5];
 162   Double_t c32 =   x*cov[13] -     cov[ 8];
 163   Double_t c20 =   x*cov[10] -     cov[ 3];
 164   Double_t c21 =   x*cov[11] -     cov[ 4];
 165   Double_t c42 =   x*cov[14] -     cov[12];
 166
 167   Double_t cc[15] = { cov[ 0]
 168                     , cov[ 1],     cov[ 2]
 169                     , c20,         c21,         c22
 170                     , cov[ 6],     cov[ 7],     c32,     cov[ 9]
 171                     , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
 172
 173   Double_t mostProbablePt=AliExternalTrackParam::GetMostProbablePt();
 174   Double_t p0=TMath::Sign(1/mostProbablePt,pp[4]);
 175   Double_t w0=cc[14]/(cc[14] + p0*p0), w1=p0*p0/(cc[14] + p0*p0);
 176   pp[4] = w0*p0 + w1*pp[4];
 177   cc[10]*=w1; cc[11]*=w1; cc[12]*=w1; cc[13]*=w1; cc[14]*=w1;
 178
 179         Set(x,alpha,pp,cc);
 180   Int_t ncls = 0;
 181         for(int iplane=0; iplane<kNplane; iplane++){
 182     fTrackletIndex[iplane] = 0xffff;
 183                 if(!trklts[iplane].IsOK()) fTracklet[iplane] = 0x0;
 184     else{
 185       fTracklet[iplane] = &trklts[iplane];
 186       ncls += fTracklet[iplane]->GetN();
 187     }
 188         }
 189   AliKalmanTrack::SetNumberOfClusters(ncls);
 190   for(int i =0; i<3; i++) fBudget[i] = 0.;
 191
 192   Float_t pid = 1./AliPID::kSPECIES;
 193   for(int is =0; is<AliPID::kSPECIES; is++) fPID[is] = pid;
 194
 195 }
 196
 197 //_______________________________________________________________
 198 AliTRDtrackV1::~AliTRDtrackV1()
 199 {
 200   //AliInfo("");
 201   //printf("I-AliTRDtrackV1::~AliTRDtrackV1() : Owner[%s]\n", TestBit(kOwner)?"YES":"NO");
 202
 203   if(fBackupTrack) delete fBackupTrack;
 204   fBackupTrack = 0x0;
 205
 206   for(Int_t ip=0; ip<kNplane; ip++){
 207     if(TestBit(kOwner) && fTracklet[ip]) delete fTracklet[ip];
 208     fTracklet[ip] = 0x0;
 209     fTrackletIndex[ip] = 0xffff;
 210   }
 211 }
 212
 213 //_______________________________________________________________
 214 Bool_t AliTRDtrackV1::CookLabel(Float_t wrong)
 215 {
 216   // set MC label for this track
 217
 218   Int_t s[kMAXCLUSTERSPERTRACK][2];
 219   for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
 220     s[i][0] = -1;
 221     s[i][1] =  0;
 222   }
 223
 224   Bool_t labelAdded;
 225   Int_t label;
 226   AliTRDcluster *c    = 0x0;
 227   for (Int_t ip = 0; ip < kNplane; ip++) {
 228     if(fTrackletIndex[ip] == 0xffff) continue;
 229     for (Int_t ic = 0; ic < AliTRDseed::knTimebins; ic++) {
 230       if(!(c = fTracklet[ip]->GetClusters(ic))) continue;
 231       for (Int_t k = 0; k < 3; k++) {
 232         label      = c->GetLabel(k);
 233         labelAdded = kFALSE;
 234         Int_t j = 0;
 235         if (label >= 0) {
 236           while ((!labelAdded) && (j < kMAXCLUSTERSPERTRACK)) {
 237             if ((s[j][0] == label) ||
 238                 (s[j][1] ==     0)) {
 239               s[j][0] = label;
 240               s[j][1]++;
 241               labelAdded = kTRUE;
 242             }
 243             j++;
 244           }
 245         }
 246       }
 247     }
 248   }
 249
 250   Int_t max = 0;
 251   label = -123456789;
 252   for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
 253     if (s[i][1] <= max) continue;
 254     max   = s[i][1];
 255     label = s[i][0];
 256   }
 257
 258   if ((1. - Float_t(max)/GetNumberOfClusters()) > wrong) label = -label;
 259
 260   SetLabel(label);
 261
 262   return kTRUE;
 263 }
 264
 265 //_______________________________________________________________
 266 Bool_t AliTRDtrackV1::CookPID()
 267 {
 268   //
 269   // Cook the PID information
 270   //
 271
 272   // Reset the a priori probabilities
 273   Double_t pid = 1. / AliPID::kSPECIES;
 274   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
 275     fPID[ispec] = pid;
 276   }
 277   fPIDquality = 0;
 278
 279   // steer PID calculation @ tracklet level
 280   Double_t *prob = 0x0;
 281   for(int ip=0; ip<kNplane; ip++){
 282     if(fTrackletIndex[ip] == 0xffff) continue;
 283     if(!fTracklet[ip]->IsOK()) continue;
 284     if(!(prob = fTracklet[ip]->GetProbability())) return kFALSE;
 285
 286     Int_t nspec = 0; // quality check of tracklet dEdx
 287     for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
 288       if(prob[ispec] < 0.) continue;
 289       fPID[ispec] *= prob[ispec];
 290       nspec++;
 291     }
 292     if(!nspec) continue;
 293
 294     fPIDquality++;
 295   }
 296
 297   // no tracklet found for PID calculations
 298   if(!fPIDquality) return kTRUE;
 299
 300   // slot for PID calculation @ track level
 301
 302
 303   // normalize probabilities
 304   Double_t probTotal = 0.0;
 305   for (Int_t is = 0; is < AliPID::kSPECIES; is++) probTotal += fPID[is];
 306
 307
 308   if (probTotal <= 0.0) {
 309     AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference data.");
 310     return kFALSE;
 311   }
 312
 313   for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) fPID[iSpecies] /= probTotal;
 314
 315   return kTRUE;
 316 }
 317
 318 //_____________________________________________________________________________
 319 Double_t AliTRDtrackV1::GetBz() const
 320 {
 321   //
 322   // Returns Bz component of the magnetic field (kG)
 323   //
 324
 325   if (AliTracker::UniformField()) return AliTracker::GetBz();
 326
 327   Double_t r[3];
 328   GetXYZ(r);
 329   return AliTracker::GetBz(r);
 330 }
 331
 332 //_______________________________________________________________
 333 AliTRDcluster* AliTRDtrackV1::GetCluster(Int_t id)
 334 {
 335   Int_t n = 0;
 336   for(Int_t ip=0; ip<kNplane; ip++){
 337     if(!fTracklet[ip]) continue;
 338     if(n+fTracklet[ip]->GetN() <= id){
 339       n+=fTracklet[ip]->GetN();
 340       continue;
 341     }
 342     AliTRDcluster *c = 0x0;
 343     for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
 344       if(!(c = fTracklet[ip]->GetClusters(ic))) continue;
 345
 346       if(n<id){n++; continue;}
 347       return c;
 348     }
 349   }
 350   return 0x0;
 351 }
 352
 353 //_______________________________________________________________
 354 Int_t  AliTRDtrackV1::GetClusterIndex(Int_t id) const
 355 {
 356   Int_t n = 0;
 357   for(Int_t ip=0; ip<kNplane; ip++){
 358     if(!fTracklet[ip]) continue;
 359     if(n+fTracklet[ip]->GetN() <= id){
 360       n+=fTracklet[ip]->GetN();
 361       continue;
 362     }
 363     AliTRDcluster *c = 0x0;
 364     for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
 365       if(!(c = fTracklet[ip]->GetClusters(ic))) continue;
 366
 367       if(n<id){n++; continue;}
 368       return fTracklet[ip]->GetIndexes(ic);
 369     }
 370   }
 371   return -1;
 372 }
 373
 374 //_______________________________________________________________
 375 Double_t AliTRDtrackV1::GetPredictedChi2(const AliTRDseedV1 *trklt) const
 376 {
 377   //
 378   // Get the predicted chi2
 379   //
 380
 381   Double_t x      = trklt->GetX0();
 382   Double_t p[2]   = { trklt->GetYat(x)
 383                     , trklt->GetZat(x) };
 384   Double_t cov[3];
 385   trklt->GetCovAt(x, cov);
 386
 387   return AliExternalTrackParam::GetPredictedChi2(p, cov);
 388 }
 389
 390 //_______________________________________________________________
 391 Bool_t AliTRDtrackV1::IsEqual(const TObject *o) const
 392 {
 393   if (!o) return kFALSE;
 394   const AliTRDtrackV1 *inTrack = dynamic_cast<const AliTRDtrackV1*>(o);
 395   if (!inTrack) return kFALSE;
 396
 397   if ( fPIDquality != inTrack->GetPIDquality() ) return kFALSE;
 398
 399   for(Int_t i = 0; i < AliPID::kSPECIES; i++){
 400     if ( fPID[i] != inTrack->GetPID(i) ) return kFALSE;
 401   }
 402
 403   for (Int_t i = 0; i < 3; i++){
 404     if ( fBudget[i] != inTrack->GetBudget(i) ) return kFALSE;
 405   }
 406   if ( fDE != inTrack->GetEdep() ) return kFALSE;
 407   if ( fFakeRatio != inTrack->GetFakeRatio() ) return kFALSE;
 408   if ( fChi2 != inTrack->GetChi2() ) return kFALSE;
 409   if ( fMass != inTrack->GetMass() ) return kFALSE;
 410   if ( fLab != inTrack->GetLabel() ) return kFALSE;
 411   if ( fN != inTrack->GetNumberOfClusters() ) return kFALSE;
 412   if ( AliKalmanTrack::GetIntegratedLength() != inTrack->GetIntegratedLength() ) return kFALSE;
 413
 414   if ( GetX() != inTrack->GetX() ) return kFALSE;
 415   if ( GetAlpha() != inTrack->GetAlpha() ) return kFALSE;
 416   const Double_t *inP = inTrack->GetParameter();
 417   const Double_t *curP = GetParameter();
 418   for (Int_t i = 0; i < 5; i++){
 419     if ( curP[i] != inP[i]) return kFALSE;
 420   }
 421   const Double_t *inC = inTrack->GetCovariance();
 422   const Double_t *curC = GetCovariance();
 423   for (Int_t i = 0; i < 15; i++){
 424     if ( curC[i] != inC[i]) return kFALSE;
 425   }
 426
 427   for (Int_t iTracklet = 0; iTracklet < kNplane; iTracklet++){
 428     AliTRDseedV1 *curTracklet = fTracklet[iTracklet];
 429     AliTRDseedV1 *inTracklet = inTrack->GetTracklet(iTracklet);
 430     if (curTracklet && inTracklet){
 431       if (! curTracklet->IsEqual(inTracklet) ) {
 432         curTracklet->Print();
 433         inTracklet->Print();
 434         return kFALSE;
 435       }
 436     } else {
 437       // if one tracklet exists, and corresponding
 438       // in other track doesn't - return kFALSE
 439       if(inTracklet || curTracklet) return kFALSE;
 440     }
 441   }
 442
 443   return kTRUE;
 444 }
 445
 446 //_______________________________________________________________
 447 Bool_t AliTRDtrackV1::IsElectron() const
 448 {
 449   if(GetPID(0) > fReconstructor->GetRecoParam()->GetPIDThreshold(GetP())) return kTRUE;
 450   return kFALSE;
 451 }
 452
 453
 454 //_____________________________________________________________________________
 455 void AliTRDtrackV1::MakeBackupTrack()
 456 {
 457   //
 458   // Creates a backup track
 459   //
 460
 461   if(fBackupTrack) {
 462     fBackupTrack->~AliTRDtrackV1();
 463     new(fBackupTrack) AliTRDtrackV1((AliTRDtrackV1&)(*this));
 464     return;
 465   }
 466   fBackupTrack = new AliTRDtrackV1((AliTRDtrackV1&)(*this));
 467 }
 468
 469 //_____________________________________________________________________________
 470 Int_t AliTRDtrackV1::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
 471 {
 472   //
 473   // Find a prolongation at given x
 474   // Return 0 if it does not exist
 475   //
 476
 477   Double_t bz = GetBz();
 478   if (!AliExternalTrackParam::GetYAt(xk,bz,y)) return 0;
 479   if (!AliExternalTrackParam::GetZAt(xk,bz,z)) return 0;
 480
 481   return 1;
 482
 483 }
 484
 485 //_____________________________________________________________________________
 486 Bool_t AliTRDtrackV1::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
 487 {
 488   //
 489   // Propagates this track to a reference plane defined by "xk" [cm]
 490   // correcting for the mean crossed material.
 491   //
 492   // "xx0"  - thickness/rad.length [units of the radiation length]
 493   // "xrho" - thickness*density    [g/cm^2]
 494   //
 495
 496   if (xk == GetX()) {
 497     return kTRUE;
 498   }
 499
 500   Double_t oldX = GetX();
 501   Double_t oldY = GetY();
 502   Double_t oldZ = GetZ();
 503
 504   Double_t bz   = GetBz();
 505
 506   if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
 507     return kFALSE;
 508   }
 509
 510   Double_t x = GetX();
 511   Double_t y = GetY();
 512   Double_t z = GetZ();
 513
 514   if (oldX < xk) {
 515     xrho = -xrho;
 516     if (IsStartedTimeIntegral()) {
 517       Double_t l2  = TMath::Sqrt((x-oldX)*(x-oldX)
 518                                + (y-oldY)*(y-oldY)
 519                                + (z-oldZ)*(z-oldZ));
 520       Double_t crv = AliExternalTrackParam::GetC(bz);
 521       if (TMath::Abs(l2*crv) > 0.0001) {
 522         // Make correction for curvature if neccesary
 523         l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX)
 524                              + (y-oldY)*(y-oldY));
 525         l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
 526         l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ));
 527       }
 528       AddTimeStep(l2);
 529     }
 530   }
 531
 532   if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0, xrho, GetMass())) return kFALSE;
 533
 534
 535   {
 536
 537     // Energy losses
 538     Double_t p2    = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt());
 539     Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
 540     if ((beta2 < 1.0e-10) ||
 541         ((5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0)) {
 542       return kFALSE;
 543     }
 544
 545     Double_t dE    = 0.153e-3 / beta2
 546                    * (TMath::Log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
 547                    * xrho;
 548     fBudget[0] += xrho;
 549
 550     /*
 551     // Suspicious part - think about it ?
 552     Double_t kinE =  TMath::Sqrt(p2);
 553     if (dE > 0.8*kinE) dE = 0.8 * kinE;  //
 554     if (dE < 0)        dE = 0.0;         // Not valid region for Bethe bloch
 555     */
 556
 557     fDE += dE;
 558
 559     /*
 560     // Suspicious ! I.B.
 561     Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE));   // Energy loss fluctuation
 562     Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
 563     fCcc += sigmac2;
 564     fCee += fX*fX * sigmac2;
 565     */
 566
 567   }
 568
 569   return kTRUE;
 570 }
 571
 572 //_____________________________________________________________________________
 573 Int_t   AliTRDtrackV1::PropagateToR(Double_t r,Double_t step)
 574 {
 575   //
 576   // Propagate track to the radial position
 577   // Rotation always connected to the last track position
 578   //
 579
 580   Double_t xyz0[3];
 581   Double_t xyz1[3];
 582   Double_t y;
 583   Double_t z;
 584
 585   Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
 586   // Direction +-
 587   Double_t dir    = (radius > r) ? -1.0 : 1.0;
 588
 589   for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
 590
 591     GetXYZ(xyz0);
 592     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 593     Rotate(alpha,kTRUE);
 594     GetXYZ(xyz0);
 595     GetProlongation(x,y,z);
 596     xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
 597     xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
 598     xyz1[2] = z;
 599     Double_t param[7];
 600     if(AliTracker::MeanMaterialBudget(xyz0,xyz1,param)<=0.) return -1;
 601     if (param[1] <= 0) {
 602       param[1] = 100000000;
 603     }
 604     PropagateTo(x,param[1],param[0]*param[4]);
 605
 606   }
 607
 608   GetXYZ(xyz0);
 609   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 610   Rotate(alpha,kTRUE);
 611   GetXYZ(xyz0);
 612   GetProlongation(r,y,z);
 613   xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
 614   xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
 615   xyz1[2] = z;
 616   Double_t param[7];
 617   if(AliTracker::MeanMaterialBudget(xyz0,xyz1,param) <= 0.) return -1;
 618
 619   if (param[1] <= 0) {
 620     param[1] = 100000000;
 621   }
 622   PropagateTo(r,param[1],param[0]*param[4]);
 623
 624   return 0;
 625
 626 }
 627
 628 //_____________________________________________________________________________
 629 void AliTRDtrackV1::Print(Option_t *o) const
 630 {
 631   AliInfo(Form("PID q[%d] [%4.1f %4.1f %4.1f %4.1f %4.1f]", fPIDquality, 1.E2*fPID[0], 1.E2*fPID[1], 1.E2*fPID[2], 1.E2*fPID[3], 1.E2*fPID[4]));
 632   AliInfo(Form("Material[%5.2f %5.2f %5.2f]", fBudget[0], fBudget[1], fBudget[2]));
 633
 634   AliInfo(Form("x[%7.2f] t[%7.4f] alpha[%f] mass[%f]", GetX(), GetIntegratedLength(), GetAlpha(), fMass));
 635   AliInfo(Form("Ntr[%1d] Ncl[%3d] lab[%3d]", GetNumberOfTracklets(), fN, fLab));
 636
 637   if(strcmp(o, "a")!=0) return;
 638   printf("|X| = (");
 639   const Double_t *curP = GetParameter();
 640   for (Int_t i = 0; i < 5; i++) printf("%7.2f ", curP[i]);
 641   printf(")\n");
 642
 643   printf("|V| = \n");
 644   const Double_t *curC = GetCovariance();
 645   for (Int_t i = 0, j=4, k=0; i<15; i++, k++){
 646     printf("%7.2f ", curC[i]);
 647     if(k==j){
 648       printf("\n");
 649       k=-1; j--;
 650     }
 651   }
 652
 653   for(Int_t ip=0; ip<kNplane; ip++){
 654     if(!fTracklet[ip]) continue;
 655     fTracklet[ip]->Print(o);
 656   }
 657 }
 658
 659
 660 //_____________________________________________________________________________
 661 Bool_t AliTRDtrackV1::Rotate(Double_t alpha, Bool_t absolute)
 662 {
 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
 669   if (absolute) alpha -= GetAlpha();
 670   //else fNRotate++;
 671
 672   return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
 673 }
 674
 675 //___________________________________________________________
 676 void AliTRDtrackV1::SetNumberOfClusters()
 677 {
 678 // Calculate the number of clusters attached to this track
 679
 680   Int_t ncls = 0;
 681   for(int ip=0; ip<kNplane; ip++){
 682     if(fTracklet[ip] && fTrackletIndex[ip] != 0xffff) ncls += fTracklet[ip]->GetN();
 683   }
 684   AliKalmanTrack::SetNumberOfClusters(ncls);
 685 }
 686
 687
 688 //_______________________________________________________________
 689 void AliTRDtrackV1::SetOwner()
 690 {
 691   //
 692   // Toggle ownership of tracklets
 693   //
 694
 695   if(TestBit(kOwner)) return;
 696   for (Int_t ip = 0; ip < kNplane; ip++) {
 697     if(fTrackletIndex[ip] == 0xffff) continue;
 698     fTracklet[ip] = new AliTRDseedV1(*fTracklet[ip]);
 699     fTracklet[ip]->SetOwner();
 700   }
 701   SetBit(kOwner);
 702 }
 703
 704 //_______________________________________________________________
 705 void AliTRDtrackV1::SetTracklet(AliTRDseedV1 *trklt, Int_t index)
 706 {
 707   //
 708   // Set the tracklets
 709   //
 710   Int_t plane = trklt->GetPlane();
 711
 712   fTracklet[plane]      = trklt;
 713   fTrackletIndex[plane] = index;
 714 }
 715
 716 //_______________________________________________________________
 717 void AliTRDtrackV1::UnsetTracklet(Int_t plane)
 718 {
 719   if(plane<0 && plane >= kNplane) return;
 720   fTrackletIndex[plane] = 0xffff;
 721   fTracklet[plane] = 0x0;
 722 }
 723
 724
 725 //_______________________________________________________________
 726 Bool_t  AliTRDtrackV1::Update(AliTRDseedV1 *trklt, Double_t chisq)
 727 {
 728   //
 729   // Update track and tracklet parameters
 730   //
 731
 732   Double_t x      = GetX();
 733   Double_t p[2]   = { trklt->GetYat(x)
 734                     , trklt->GetZat(x) };
 735   Double_t cov[3]/*, covR[3], cov0[3]*/;
 736
 737 //   printf("\tD[%3d] Ly[%d] Trk: x[%f] y[%f] z[%f]\n", trklt->GetDetector(), trklt->GetPlane(), GetX(), GetY(), GetZ());
 738 // //
 739 //   Double_t xref = trklt->GetXref();
 740 //   trklt->GetCovAt(xref, covR);
 741 //   printf("xr=%5.3f y=%f+-%f z=%f+-%f (covYZ=%f)\n", xref, trklt->GetYat(xref), TMath::Sqrt(covR[0]), trklt->GetZat(xref), TMath::Sqrt(covR[2]), covR[1]);
 742 //
 743 //   Double_t x0 = trklt->GetX0();
 744 //   trklt->GetCovAt(x0, cov0);
 745 //   printf("x0=%5.3f y=%f+-%f z=%f+-%f (covYZ=%f)\n", x0, trklt->GetYat(x0), TMath::Sqrt(cov0[0]), trklt->GetZat(x0), TMath::Sqrt(cov0[2]), cov0[1]);
 746 //
 747 //   trklt->GetCovAt(x, cov);
 748 //   printf("x =%5.3f y=%f+-%f z=%f+-%f (covYZ=%f)\n", x, p[0], TMath::Sqrt(cov[0]), p[1], TMath::Sqrt(cov[2]), cov[1]);
 749 //
 750 //   const Double_t *cc = GetCovariance();
 751 //   printf("yklm[0] = %f +- %f\n", GetY(), TMath::Sqrt(cc[0]));
 752
 753   trklt->GetCovAt(x, cov);
 754   if(!AliExternalTrackParam::Update(p, cov)) return kFALSE;
 755 //   cc = GetCovariance();
 756 //   printf("yklm[1] = %f +- %f\n", GetY(), TMath::Sqrt(cc[0]));
 757
 758   AliTRDcluster *c = 0x0;
 759   Int_t ic = 0; while(!(c = trklt->GetClusters(ic))) ic++;
 760   AliTracker::FillResiduals(this, p, cov, c->GetVolumeId());
 761
 762   // Register info to track
 763   SetNumberOfClusters();
 764   SetChi2(GetChi2() + chisq);
 765   return kTRUE;
 766 }
 767
 768 //_______________________________________________________________
 769 void AliTRDtrackV1::UpdateESDtrack(AliESDtrack *track)
 770 {
 771   //
 772   // Update the TRD PID information in the ESD track
 773   //
 774
 775   Int_t nslices = fReconstructor->IsEightSlices() ? (Int_t)AliTRDReconstructor::kNNslices : (Int_t)AliTRDReconstructor::kLQslices;
 776   track->SetNumberOfTRDslices(nslices);
 777
 778   for (Int_t ip = 0; ip < kNplane; ip++) {
 779     if(fTrackletIndex[ip] == 0xffff) continue;
 780     if(!fTracklet[ip]->IsOK()) continue;
 781     fTracklet[ip]->CookdEdx(nslices);
 782     Float_t *dedx = fTracklet[ip]->GetdEdx();
 783     for (Int_t js = 0; js < nslices; js++, dedx++) track->SetTRDslice(*dedx, ip, js);
 784   }
 785
 786   // copy PID to ESD
 787   if(!fPIDquality) return;
 788   track->SetTRDpid(fPID);
 789   track->SetTRDpidQuality(fPIDquality);
 790 }