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