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 AliTRDcluster* AliTRDtrackV1::GetCluster(Int_t id)
 320 {
 321   Int_t n = 0;
 322   for(Int_t ip=0; ip<kNplane; ip++){
 323     if(!fTracklet[ip]) continue;
 324     if(n+fTracklet[ip]->GetN() <= id){
 325       n+=fTracklet[ip]->GetN();
 326       continue;
 327     }
 328     AliTRDcluster *c = 0x0;
 329     for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
 330       if(!(c = fTracklet[ip]->GetClusters(ic))) continue;
 331
 332       if(n<id){n++; continue;}
 333       return c;
 334     }
 335   }
 336   return 0x0;
 337 }
 338
 339 //_______________________________________________________________
 340 Int_t  AliTRDtrackV1::GetClusterIndex(Int_t id) const
 341 {
 342   Int_t n = 0;
 343   for(Int_t ip=0; ip<kNplane; ip++){
 344     if(!fTracklet[ip]) continue;
 345     if(n+fTracklet[ip]->GetN() <= id){
 346       n+=fTracklet[ip]->GetN();
 347       continue;
 348     }
 349     AliTRDcluster *c = 0x0;
 350     for(Int_t ic=AliTRDseed::knTimebins-1; ic>=0; ic--){
 351       if(!(c = fTracklet[ip]->GetClusters(ic))) continue;
 352
 353       if(n<id){n++; continue;}
 354       return fTracklet[ip]->GetIndexes(ic);
 355     }
 356   }
 357   return -1;
 358 }
 359
 360 //_______________________________________________________________
 361 Double_t AliTRDtrackV1::GetPredictedChi2(const AliTRDseedV1 *trklt) const
 362 {
 363   //
 364   // Get the predicted chi2
 365   //
 366
 367   Double_t x      = trklt->GetX0();
 368   Double_t p[2]   = { trklt->GetYat(x)
 369                     , trklt->GetZat(x) };
 370   Double_t cov[3];
 371   trklt->GetCovAt(x, cov);
 372
 373   return AliExternalTrackParam::GetPredictedChi2(p, cov);
 374 }
 375
 376 //_______________________________________________________________
 377 Bool_t AliTRDtrackV1::IsEqual(const TObject *o) const
 378 {
 379   if (!o) return kFALSE;
 380   const AliTRDtrackV1 *inTrack = dynamic_cast<const AliTRDtrackV1*>(o);
 381   if (!inTrack) return kFALSE;
 382
 383   if ( fPIDquality != inTrack->GetPIDquality() ) return kFALSE;
 384
 385   for(Int_t i = 0; i < AliPID::kSPECIES; i++){
 386     if ( fPID[i] != inTrack->GetPID(i) ) return kFALSE;
 387   }
 388
 389   for (Int_t i = 0; i < 3; i++){
 390     if ( fBudget[i] != inTrack->GetBudget(i) ) return kFALSE;
 391   }
 392   if ( fDE != inTrack->GetEdep() ) return kFALSE;
 393   if ( fFakeRatio != inTrack->GetFakeRatio() ) return kFALSE;
 394   if ( fChi2 != inTrack->GetChi2() ) return kFALSE;
 395   if ( fMass != inTrack->GetMass() ) return kFALSE;
 396   if ( fLab != inTrack->GetLabel() ) return kFALSE;
 397   if ( fN != inTrack->GetNumberOfClusters() ) return kFALSE;
 398   if ( AliKalmanTrack::GetIntegratedLength() != inTrack->GetIntegratedLength() ) return kFALSE;
 399
 400   if ( GetX() != inTrack->GetX() ) return kFALSE;
 401   if ( GetAlpha() != inTrack->GetAlpha() ) return kFALSE;
 402   const Double_t *inP = inTrack->GetParameter();
 403   const Double_t *curP = GetParameter();
 404   for (Int_t i = 0; i < 5; i++){
 405     if ( curP[i] != inP[i]) return kFALSE;
 406   }
 407   const Double_t *inC = inTrack->GetCovariance();
 408   const Double_t *curC = GetCovariance();
 409   for (Int_t i = 0; i < 15; i++){
 410     if ( curC[i] != inC[i]) return kFALSE;
 411   }
 412
 413   for (Int_t iTracklet = 0; iTracklet < kNplane; iTracklet++){
 414     AliTRDseedV1 *curTracklet = fTracklet[iTracklet];
 415     AliTRDseedV1 *inTracklet = inTrack->GetTracklet(iTracklet);
 416     if (curTracklet && inTracklet){
 417       if (! curTracklet->IsEqual(inTracklet) ) {
 418         curTracklet->Print();
 419         inTracklet->Print();
 420         return kFALSE;
 421       }
 422     } else {
 423       // if one tracklet exists, and corresponding
 424       // in other track doesn't - return kFALSE
 425       if(inTracklet || curTracklet) return kFALSE;
 426     }
 427   }
 428
 429   return kTRUE;
 430 }
 431
 432 //_______________________________________________________________
 433 Bool_t AliTRDtrackV1::IsElectron() const
 434 {
 435   if(GetPID(0) > fReconstructor->GetRecoParam()->GetPIDThreshold(GetP())) return kTRUE;
 436   return kFALSE;
 437 }
 438
 439
 440 //_____________________________________________________________________________
 441 void AliTRDtrackV1::MakeBackupTrack()
 442 {
 443   //
 444   // Creates a backup track
 445   //
 446
 447   if(fBackupTrack) {
 448     fBackupTrack->~AliTRDtrackV1();
 449     new(fBackupTrack) AliTRDtrackV1((AliTRDtrackV1&)(*this));
 450     return;
 451   }
 452   fBackupTrack = new AliTRDtrackV1((AliTRDtrackV1&)(*this));
 453 }
 454
 455 //_____________________________________________________________________________
 456 Int_t AliTRDtrackV1::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
 457 {
 458   //
 459   // Find a prolongation at given x
 460   // Return 0 if it does not exist
 461   //
 462
 463   Double_t bz = GetBz();
 464   if (!AliExternalTrackParam::GetYAt(xk,bz,y)) return 0;
 465   if (!AliExternalTrackParam::GetZAt(xk,bz,z)) return 0;
 466
 467   return 1;
 468
 469 }
 470
 471 //_____________________________________________________________________________
 472 Bool_t AliTRDtrackV1::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
 473 {
 474   //
 475   // Propagates this track to a reference plane defined by "xk" [cm]
 476   // correcting for the mean crossed material.
 477   //
 478   // "xx0"  - thickness/rad.length [units of the radiation length]
 479   // "xrho" - thickness*density    [g/cm^2]
 480   //
 481
 482   if (xk == GetX()) {
 483     return kTRUE;
 484   }
 485
 486   Double_t oldX = GetX();
 487   Double_t oldY = GetY();
 488   Double_t oldZ = GetZ();
 489
 490   Double_t bz   = GetBz();
 491
 492   if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
 493     return kFALSE;
 494   }
 495
 496   Double_t x = GetX();
 497   Double_t y = GetY();
 498   Double_t z = GetZ();
 499
 500   if (oldX < xk) {
 501     xrho = -xrho;
 502     if (IsStartedTimeIntegral()) {
 503       Double_t l2  = TMath::Sqrt((x-oldX)*(x-oldX)
 504                                + (y-oldY)*(y-oldY)
 505                                + (z-oldZ)*(z-oldZ));
 506       Double_t crv = AliExternalTrackParam::GetC(bz);
 507       if (TMath::Abs(l2*crv) > 0.0001) {
 508         // Make correction for curvature if neccesary
 509         l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX)
 510                              + (y-oldY)*(y-oldY));
 511         l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
 512         l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ));
 513       }
 514       AddTimeStep(l2);
 515     }
 516   }
 517
 518   if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0, xrho, GetMass())) return kFALSE;
 519
 520
 521   {
 522
 523     // Energy losses
 524     Double_t p2    = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt());
 525     Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
 526     if ((beta2 < 1.0e-10) ||
 527         ((5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0)) {
 528       return kFALSE;
 529     }
 530
 531     Double_t dE    = 0.153e-3 / beta2
 532                    * (TMath::Log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
 533                    * xrho;
 534     fBudget[0] += xrho;
 535
 536     /*
 537     // Suspicious part - think about it ?
 538     Double_t kinE =  TMath::Sqrt(p2);
 539     if (dE > 0.8*kinE) dE = 0.8 * kinE;  //
 540     if (dE < 0)        dE = 0.0;         // Not valid region for Bethe bloch
 541     */
 542
 543     fDE += dE;
 544
 545     /*
 546     // Suspicious ! I.B.
 547     Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE));   // Energy loss fluctuation
 548     Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
 549     fCcc += sigmac2;
 550     fCee += fX*fX * sigmac2;
 551     */
 552
 553   }
 554
 555   return kTRUE;
 556 }
 557
 558 //_____________________________________________________________________________
 559 Int_t   AliTRDtrackV1::PropagateToR(Double_t r,Double_t step)
 560 {
 561   //
 562   // Propagate track to the radial position
 563   // Rotation always connected to the last track position
 564   //
 565
 566   Double_t xyz0[3];
 567   Double_t xyz1[3];
 568   Double_t y;
 569   Double_t z;
 570
 571   Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
 572   // Direction +-
 573   Double_t dir    = (radius > r) ? -1.0 : 1.0;
 574
 575   for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
 576
 577     GetXYZ(xyz0);
 578     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 579     Rotate(alpha,kTRUE);
 580     GetXYZ(xyz0);
 581     GetProlongation(x,y,z);
 582     xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha);
 583     xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
 584     xyz1[2] = z;
 585     Double_t param[7];
 586     if(AliTracker::MeanMaterialBudget(xyz0,xyz1,param)<=0.) return -1;
 587     if (param[1] <= 0) {
 588       param[1] = 100000000;
 589     }
 590     PropagateTo(x,param[1],param[0]*param[4]);
 591
 592   }
 593
 594   GetXYZ(xyz0);
 595   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 596   Rotate(alpha,kTRUE);
 597   GetXYZ(xyz0);
 598   GetProlongation(r,y,z);
 599   xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha);
 600   xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
 601   xyz1[2] = z;
 602   Double_t param[7];
 603   if(AliTracker::MeanMaterialBudget(xyz0,xyz1,param) <= 0.) return -1;
 604
 605   if (param[1] <= 0) {
 606     param[1] = 100000000;
 607   }
 608   PropagateTo(r,param[1],param[0]*param[4]);
 609
 610   return 0;
 611
 612 }
 613
 614 //_____________________________________________________________________________
 615 void AliTRDtrackV1::Print(Option_t *o) const
 616 {
 617   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]));
 618   AliInfo(Form("Material[%5.2f %5.2f %5.2f]", fBudget[0], fBudget[1], fBudget[2]));
 619
 620   AliInfo(Form("x[%7.2f] t[%7.4f] alpha[%f] mass[%f]", GetX(), GetIntegratedLength(), GetAlpha(), fMass));
 621   AliInfo(Form("Ntr[%1d] Ncl[%3d] lab[%3d]", GetNumberOfTracklets(), fN, fLab));
 622
 623   if(strcmp(o, "a")!=0) return;
 624   printf("|X| = (");
 625   const Double_t *curP = GetParameter();
 626   for (Int_t i = 0; i < 5; i++) printf("%7.2f ", curP[i]);
 627   printf(")\n");
 628
 629   printf("|V| = \n");
 630   const Double_t *curC = GetCovariance();
 631   for (Int_t i = 0, j=4, k=0; i<15; i++, k++){
 632     printf("%7.2f ", curC[i]);
 633     if(k==j){
 634       printf("\n");
 635       k=-1; j--;
 636     }
 637   }
 638
 639   for(Int_t ip=0; ip<kNplane; ip++){
 640     if(!fTracklet[ip]) continue;
 641     fTracklet[ip]->Print(o);
 642   }
 643 }
 644
 645
 646 //_____________________________________________________________________________
 647 Bool_t AliTRDtrackV1::Rotate(Double_t alpha, Bool_t absolute)
 648 {
 649   //
 650   // Rotates track parameters in R*phi plane
 651   // if absolute rotation alpha is in global system
 652   // otherwise alpha rotation is relative to the current rotation angle
 653   //
 654
 655   if (absolute) alpha -= GetAlpha();
 656   //else fNRotate++;
 657
 658   return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
 659 }
 660
 661 //___________________________________________________________
 662 void AliTRDtrackV1::SetNumberOfClusters()
 663 {
 664 // Calculate the number of clusters attached to this track
 665
 666   Int_t ncls = 0;
 667   for(int ip=0; ip<kNplane; ip++){
 668     if(fTracklet[ip] && fTrackletIndex[ip] != 0xffff) ncls += fTracklet[ip]->GetN();
 669   }
 670   AliKalmanTrack::SetNumberOfClusters(ncls);
 671 }
 672
 673
 674 //_______________________________________________________________
 675 void AliTRDtrackV1::SetOwner()
 676 {
 677   //
 678   // Toggle ownership of tracklets
 679   //
 680
 681   if(TestBit(kOwner)) return;
 682   for (Int_t ip = 0; ip < kNplane; ip++) {
 683     if(fTrackletIndex[ip] == 0xffff) continue;
 684     fTracklet[ip] = new AliTRDseedV1(*fTracklet[ip]);
 685     fTracklet[ip]->SetOwner();
 686   }
 687   SetBit(kOwner);
 688 }
 689
 690 //_______________________________________________________________
 691 void AliTRDtrackV1::SetTracklet(AliTRDseedV1 *trklt, Int_t index)
 692 {
 693   //
 694   // Set the tracklets
 695   //
 696   Int_t plane = trklt->GetPlane();
 697
 698   fTracklet[plane]      = trklt;
 699   fTrackletIndex[plane] = index;
 700 }
 701
 702 //_______________________________________________________________
 703 void AliTRDtrackV1::UnsetTracklet(Int_t plane)
 704 {
 705   if(plane<0 && plane >= kNplane) return;
 706   fTrackletIndex[plane] = 0xffff;
 707   fTracklet[plane] = 0x0;
 708 }
 709
 710
 711 //_______________________________________________________________
 712 Bool_t  AliTRDtrackV1::Update(AliTRDseedV1 *trklt, Double_t chisq)
 713 {
 714   //
 715   // Update track and tracklet parameters
 716   //
 717
 718   Double_t x      = GetX();
 719   Double_t p[2]   = { trklt->GetYat(x)
 720                     , trklt->GetZat(x) };
 721   Double_t cov[3]/*, covR[3], cov0[3]*/;
 722
 723 //   printf("\tD[%3d] Ly[%d] Trk: x[%f] y[%f] z[%f]\n", trklt->GetDetector(), trklt->GetPlane(), GetX(), GetY(), GetZ());
 724 // //
 725 //   Double_t xref = trklt->GetXref();
 726 //   trklt->GetCovAt(xref, covR);
 727 //   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]);
 728 //
 729 //   Double_t x0 = trklt->GetX0();
 730 //   trklt->GetCovAt(x0, cov0);
 731 //   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]);
 732 //
 733 //   trklt->GetCovAt(x, cov);
 734 //   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]);
 735 //
 736 //   const Double_t *cc = GetCovariance();
 737 //   printf("yklm[0] = %f +- %f\n", GetY(), TMath::Sqrt(cc[0]));
 738
 739   trklt->GetCovAt(x, cov);
 740   if(!AliExternalTrackParam::Update(p, cov)) return kFALSE;
 741 //   cc = GetCovariance();
 742 //   printf("yklm[1] = %f +- %f\n", GetY(), TMath::Sqrt(cc[0]));
 743
 744   AliTRDcluster *c = 0x0;
 745   Int_t ic = 0; while(!(c = trklt->GetClusters(ic))) ic++;
 746   AliTracker::FillResiduals(this, p, cov, c->GetVolumeId());
 747
 748   // Register info to track
 749   SetNumberOfClusters();
 750   SetChi2(GetChi2() + chisq);
 751   return kTRUE;
 752 }
 753
 754 //_______________________________________________________________
 755 void AliTRDtrackV1::UpdateESDtrack(AliESDtrack *track)
 756 {
 757   //
 758   // Update the TRD PID information in the ESD track
 759   //
 760
 761   Int_t nslices = fReconstructor->IsEightSlices() ? (Int_t)AliTRDpidUtil::kNNslices : (Int_t)AliTRDpidUtil::kLQslices;
 762   track->SetNumberOfTRDslices(nslices);
 763
 764   for (Int_t ip = 0; ip < kNplane; ip++) {
 765     if(fTrackletIndex[ip] == 0xffff) continue;
 766     if(!fTracklet[ip]->IsOK()) continue;
 767     fTracklet[ip]->CookdEdx(nslices);
 768     Float_t *dedx = fTracklet[ip]->GetdEdx();
 769     for (Int_t js = 0; js < nslices; js++, dedx++) track->SetTRDslice(*dedx, ip, js);
 770   }
 771
 772   // copy PID to ESD
 773   if(!fPIDquality) return;
 774   track->SetTRDpid(fPID);
 775   track->SetTRDpidQuality(fPIDquality);
 776 }