Update in the interface for HLT and for visualization
[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 "AliESDtrack.h"
19 #include "AliTracker.h"
20
21 #include "AliTRDtrackV1.h"
22 #include "AliTRDcluster.h"
23 #include "AliTRDcalibDB.h"
24 #include "AliTRDrecoParam.h"
25
26 ClassImp(AliTRDtrackV1)
27
28 ///////////////////////////////////////////////////////////////////////////////
29 //                                                                           //
30 //  Represents a reconstructed TRD track                                     //
31 //  Local TRD Kalman track                                                   //
32 //                                                                           //
33 //  Authors:                                                                 //
34 //    Alex Bercuci <A.Bercuci@gsi.de>                                        //
35 //    Markus Fasel <M.Fasel@gsi.de>                                          //
36 //                                                                           //
37 ///////////////////////////////////////////////////////////////////////////////
38
39 //_______________________________________________________________
40 AliTRDtrackV1::AliTRDtrackV1() : AliKalmanTrack()
41   ,fPIDquality(0)
42   ,fDE(0.)
43   ,fBackupTrack(0x0)
44 {
45   //
46   // Default constructor
47   //
48   for(int i =0; i<3; i++) fBudget[i] = 0.;
49   
50   Float_t pid = 1./AliPID::kSPECIES;
51   for(int is =0; is<AliPID::kSPECIES; is++) fPID[is] = pid;
52
53   for(int ip=0; ip<kNplane; ip++){
54     fTrackletIndex[ip] = 0xffff;
55     fTracklet[ip]      = 0x0;
56   }
57 }
58
59 //_______________________________________________________________
60 AliTRDtrackV1::AliTRDtrackV1(const AliTRDtrackV1 &ref) : AliKalmanTrack(ref)
61   ,fPIDquality(ref.fPIDquality)
62   ,fDE(ref.fDE)
63   ,fBackupTrack(0x0)
64 {
65   //
66   // Copy constructor
67   //
68
69   for(int ip=0; ip<kNplane; ip++){ 
70     fTrackletIndex[ip] = ref.fTrackletIndex[ip];
71     fTracklet[ip]      = ref.fTracklet[ip];
72   }
73
74   for (Int_t i = 0; i < 3;i++) fBudget[i]      = ref.fBudget[i];
75
76         for(Int_t is = 0; is<AliPID::kSPECIES; is++) fPID[is] = ref.fPID[is];
77   
78   AliKalmanTrack::SetNumberOfClusters(ref.GetNumberOfClusters());
79 }
80
81 //_______________________________________________________________
82 AliTRDtrackV1::AliTRDtrackV1(const AliESDtrack &t) : AliKalmanTrack()
83   ,fPIDquality(0)
84   ,fDE(0.)
85   ,fBackupTrack(0x0)
86 {
87   //
88   // Constructor from AliESDtrack
89   //
90
91   SetLabel(t.GetLabel());
92   SetChi2(0.0);
93   SetMass(t.GetMass());
94   AliKalmanTrack::SetNumberOfClusters(t.GetTRDncls()); 
95   Int_t ti[kNplane]; t.GetTRDtracklets(&ti[0]);
96   for(int ip=0; ip<kNplane; ip++){ 
97     fTrackletIndex[ip] = ti[ip] < 0 ? 0xffff : ti[ip];
98     fTracklet[ip]      = 0x0;
99   }
100   for(int i =0; i<3; i++) fBudget[i] = 0.;
101   
102   Float_t pid = 1./AliPID::kSPECIES;
103   for(int is =0; is<AliPID::kSPECIES; is++) fPID[is] = pid;
104
105   const AliExternalTrackParam *par = &t;
106   if (t.GetStatus() & AliESDtrack::kTRDbackup) { 
107     par = t.GetOuterParam();
108     if (!par) {
109       AliError("No backup info!"); 
110       par = &t;
111     }
112   }
113   Set(par->GetX() 
114      ,par->GetAlpha()
115      ,par->GetParameter()
116      ,par->GetCovariance());
117
118   if(t.GetStatus() & AliESDtrack::kTIME) {
119     StartTimeIntegral();
120     Double_t times[10]; 
121     t.GetIntegratedTimes(times); 
122     SetIntegratedTimes(times);
123     SetIntegratedLength(t.GetIntegratedLength());
124   }
125 }
126
127 //_______________________________________________________________
128 AliTRDtrackV1::AliTRDtrackV1(AliTRDseedV1 *trklts, const Double_t p[5], const Double_t cov[15]
129              , Double_t x, Double_t alpha) : AliKalmanTrack()
130   ,fPIDquality(0)
131   ,fDE(0.)
132   ,fBackupTrack(0x0)
133 {
134   //
135   // The stand alone tracking constructor
136   // TEMPORARY !!!!!!!!!!!
137   // to check :
138   // 1. covariance matrix
139   // 2. dQdl calculation
140   //
141
142   Double_t cnv   = 1.0 / (GetBz() * kB2C);
143
144   Double_t pp[5] = { p[0]    
145                    , p[1]
146                    , x*p[4] - p[2]
147                    , p[3]
148                    , p[4]*cnv      };
149
150   Double_t c22 = x*x*cov[14] - 2*x*cov[12] + cov[ 5];
151   Double_t c32 =   x*cov[13] -     cov[ 8];
152   Double_t c20 =   x*cov[10] -     cov[ 3];
153   Double_t c21 =   x*cov[11] -     cov[ 4];
154   Double_t c42 =   x*cov[14] -     cov[12];
155
156   Double_t cc[15] = { cov[ 0]
157                     , cov[ 1],     cov[ 2]
158                     , c20,         c21,         c22
159                     , cov[ 6],     cov[ 7],     c32,     cov[ 9]
160                     , cov[10]*cnv, cov[11]*cnv, c42*cnv, cov[13]*cnv, cov[14]*cnv*cnv };
161   
162         Set(x,alpha,pp,cc);
163         for(int iplane=0; iplane<kNplane; iplane++){
164                 fTrackletIndex[iplane] = 0xffff;
165                 fTracklet[iplane] = &trklts[iplane];
166         }
167         SetNumberOfClusters();
168 }
169
170 //_______________________________________________________________
171 AliTRDtrackV1::~AliTRDtrackV1()
172 {
173   if(fBackupTrack) {
174     delete fBackupTrack;
175   }
176   fBackupTrack = 0x0;
177
178   if(TestBit(1)){
179     for(Int_t ip=0; ip<kNplane; ip++){
180       if(fTracklet[ip]) delete fTracklet[ip];
181       fTracklet[ip] = 0x0;
182     }
183   }
184 }
185         
186 //_______________________________________________________________
187 Bool_t AliTRDtrackV1::CookLabel(Float_t wrong)
188 {
189   // set MC label for this track
190   
191   Int_t s[kMAXCLUSTERSPERTRACK][2];
192   for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
193     s[i][0] = -1;
194     s[i][1] =  0;
195   }
196   
197   Bool_t labelAdded;
198   Int_t label;
199   AliTRDcluster *c    = 0x0;
200   for (Int_t ip = 0; ip < kNplane; ip++) {
201     if(fTrackletIndex[ip] == 0xffff) continue;
202     for (Int_t ic = 0; ic < AliTRDseed::knTimebins; ic++) {
203       if(!(c = fTracklet[ip]->GetClusters(ic))) continue;
204       for (Int_t k = 0; k < 3; k++) { 
205         label      = c->GetLabel(k);
206         labelAdded = kFALSE; 
207         Int_t j = 0;
208         if (label >= 0) {
209           while ((!labelAdded) && (j < kMAXCLUSTERSPERTRACK)) {
210             if ((s[j][0] == label) || 
211                 (s[j][1] ==     0)) {
212               s[j][0] = label; 
213               s[j][1]++; 
214               labelAdded = kTRUE;
215             }
216             j++;
217           }
218         }
219       }
220     }
221   }
222   
223   Int_t max = 0;
224   label = -123456789;
225   for (Int_t i = 0; i < kMAXCLUSTERSPERTRACK; i++) {
226     if (s[i][1] <= max) continue;
227     max   = s[i][1]; 
228     label = s[i][0];
229   }
230   
231   if ((1. - Float_t(max)/GetNumberOfClusters()) > wrong) label = -label;
232   
233   SetLabel(label); 
234   
235   return kTRUE;
236 }
237
238 //_______________________________________________________________
239 Bool_t AliTRDtrackV1::CookPID()
240 {
241   //
242   // Cook the PID information
243   //
244   
245   // Reset the a priori probabilities
246   Double_t pid = 1. / AliPID::kSPECIES;
247   for(int ispec=0; ispec<AliPID::kSPECIES; ispec++) {
248     fPID[ispec] = pid;  
249   }
250   fPIDquality = 0;
251   
252   // steer PID calculation @ tracklet level
253   Double_t *prob = 0x0;
254   for(int ip=0; ip<kNplane; ip++){
255     if(fTrackletIndex[ip] == 0xffff) continue;
256     if(!fTracklet[ip]->IsOK()) continue;
257     if(!(prob = fTracklet[ip]->GetProbability())) return kFALSE;
258     
259     Int_t nspec = 0; // quality check of tracklet dEdx
260     for(int ispec=0; ispec<AliPID::kSPECIES; ispec++){
261       if(prob[ispec] < 0.) continue;
262       fPID[ispec] *= prob[ispec];
263       nspec++;
264     }
265     if(!nspec) continue;
266     
267     fPIDquality++;
268   }
269   
270   // no tracklet found for PID calculations
271   if(!fPIDquality) return kTRUE;
272   
273   // slot for PID calculation @ track level
274   
275   
276   // normalize probabilities
277   Double_t probTotal = 0.0;
278   for (Int_t is = 0; is < AliPID::kSPECIES; is++) probTotal += fPID[is];
279   
280   
281   if (probTotal <= 0.0) {
282     AliWarning("The total probability over all species <= 0. This may be caused by some error in the reference data.");
283     return kFALSE;
284   }
285   
286   for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++) fPID[iSpecies] /= probTotal;
287   
288   return kTRUE;
289 }
290
291 //_____________________________________________________________________________
292 Double_t AliTRDtrackV1::GetBz() const 
293 {
294   //
295   // Returns Bz component of the magnetic field (kG)
296   //
297
298   if (AliTracker::UniformField()) return AliTracker::GetBz();
299
300   Double_t r[3]; 
301   GetXYZ(r);
302   return AliTracker::GetBz(r);
303 }
304
305 //_______________________________________________________________
306 Int_t  AliTRDtrackV1::GetClusterIndex(Int_t id) const
307 {
308   Int_t n = 0;
309   for(Int_t ip=0; ip<kNplane; ip++){
310     if(!fTracklet[ip]) continue;
311     if(n+fTracklet[ip]->GetN() < id){ 
312       n+=fTracklet[ip]->GetN();
313       continue;
314     }
315     for(Int_t ic=34; ic>=0; ic--){
316       if(!fTracklet[ip]->GetClusters(ic)) continue;
317       n++;
318       if(n<id) continue;
319       return fTracklet[ip]->GetIndexes(ic);
320     }
321   }
322   return -1;
323 }
324
325 //_______________________________________________________________
326 Double_t AliTRDtrackV1::GetPredictedChi2(const AliTRDseedV1 *trklt) const
327 {
328   //
329   // Get the predicted chi2
330   //
331   
332   Double_t x      = trklt->GetX0();
333   Double_t p[2]   = { trklt->GetYat(x)
334                     , trklt->GetZat(x) };
335   Double_t cov[3];
336   trklt->GetCovAt(x, cov);
337   
338   return AliExternalTrackParam::GetPredictedChi2(p, cov);
339 }
340
341         
342 //_____________________________________________________________________________
343 void AliTRDtrackV1::MakeBackupTrack()
344 {
345   //
346   // Creates a backup track
347   //
348
349   if(fBackupTrack) {
350     fBackupTrack->~AliTRDtrackV1();
351     new(fBackupTrack) AliTRDtrackV1((AliTRDtrackV1&)(*this));
352   }
353   fBackupTrack = new AliTRDtrackV1((AliTRDtrackV1&)(*this));
354 }
355
356 //_____________________________________________________________________________
357 Int_t AliTRDtrackV1::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
358 {
359   //
360   // Find a prolongation at given x
361   // Return 0 if it does not exist
362   //  
363
364   Double_t bz = GetBz();
365   if (!AliExternalTrackParam::GetYAt(xk,bz,y)) return 0;
366   if (!AliExternalTrackParam::GetZAt(xk,bz,z)) return 0;
367
368   return 1;  
369
370 }
371
372 //_____________________________________________________________________________
373 Bool_t AliTRDtrackV1::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
374 {
375   //
376   // Propagates this track to a reference plane defined by "xk" [cm] 
377   // correcting for the mean crossed material.
378   //
379   // "xx0"  - thickness/rad.length [units of the radiation length] 
380   // "xrho" - thickness*density    [g/cm^2] 
381   // 
382
383   if (xk == GetX()) {
384     return kTRUE;
385   }
386
387   Double_t oldX = GetX();
388   Double_t oldY = GetY();
389   Double_t oldZ = GetZ();
390
391   Double_t bz   = GetBz();
392
393   if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
394     return kFALSE;
395   }
396
397   Double_t x = GetX();
398   Double_t y = GetY();
399   Double_t z = GetZ();
400
401   if (oldX < xk) {
402     xrho = -xrho;
403     if (IsStartedTimeIntegral()) {
404       Double_t l2  = TMath::Sqrt((x-oldX)*(x-oldX) 
405                                + (y-oldY)*(y-oldY) 
406                                + (z-oldZ)*(z-oldZ));
407       Double_t crv = AliExternalTrackParam::GetC(bz);
408       if (TMath::Abs(l2*crv) > 0.0001) {
409         // Make correction for curvature if neccesary
410         l2 = 0.5 * TMath::Sqrt((x-oldX)*(x-oldX) 
411                              + (y-oldY)*(y-oldY));
412         l2 = 2.0 * TMath::ASin(l2 * crv) / crv;
413         l2 = TMath::Sqrt(l2*l2 + (z-oldZ)*(z-oldZ));
414       }
415       AddTimeStep(l2);
416     }
417   }
418
419   if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0,xrho,GetMass())) { 
420     return kFALSE;
421   }
422
423   {
424
425     // Energy losses
426     Double_t p2    = (1.0 + GetTgl()*GetTgl()) / (GetSigned1Pt()*GetSigned1Pt());
427     Double_t beta2 = p2 / (p2 + GetMass()*GetMass());
428     if ((beta2 < 1.0e-10) || 
429         ((5940.0 * beta2/(1.0 - beta2 + 1.0e-10) - beta2) < 0.0)) {
430       return kFALSE;
431     }
432
433     Double_t dE    = 0.153e-3 / beta2 
434                    * (TMath::Log(5940.0 * beta2/(1.0 - beta2 + 1.0e-10)) - beta2)
435                    * xrho;
436     fBudget[0] += xrho;
437
438     /*
439     // Suspicious part - think about it ?
440     Double_t kinE =  TMath::Sqrt(p2);
441     if (dE > 0.8*kinE) dE = 0.8 * kinE;  //      
442     if (dE < 0)        dE = 0.0;         // Not valid region for Bethe bloch 
443     */
444  
445     fDE += dE;
446
447     /*
448     // Suspicious ! I.B.
449     Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE));   // Energy loss fluctuation 
450     Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
451     fCcc += sigmac2;
452     fCee += fX*fX * sigmac2;  
453     */
454
455   }
456
457   return kTRUE;
458 }
459
460 //_____________________________________________________________________________
461 Int_t   AliTRDtrackV1::PropagateToR(Double_t r,Double_t step)
462 {
463   //
464   // Propagate track to the radial position
465   // Rotation always connected to the last track position
466   //
467
468   Double_t xyz0[3];
469   Double_t xyz1[3];
470   Double_t y;
471   Double_t z; 
472
473   Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
474   // Direction +-
475   Double_t dir    = (radius > r) ? -1.0 : 1.0;   
476
477   for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
478
479     GetXYZ(xyz0);       
480     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
481     Rotate(alpha,kTRUE);
482     GetXYZ(xyz0);       
483     GetProlongation(x,y,z);
484     xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha); 
485     xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
486     xyz1[2] = z;
487     Double_t param[7];
488     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
489     if (param[1] <= 0) {
490       param[1] = 100000000;
491     }
492     PropagateTo(x,param[1],param[0]*param[4]);
493
494   } 
495
496   GetXYZ(xyz0); 
497   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
498   Rotate(alpha,kTRUE);
499   GetXYZ(xyz0); 
500   GetProlongation(r,y,z);
501   xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha); 
502   xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
503   xyz1[2] = z;
504   Double_t param[7];
505   AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
506
507   if (param[1] <= 0) {
508     param[1] = 100000000;
509   }
510   PropagateTo(r,param[1],param[0]*param[4]);
511
512   return 0;
513
514 }
515
516 //_____________________________________________________________________________
517 Bool_t AliTRDtrackV1::Rotate(Double_t alpha, Bool_t absolute)
518 {
519   //
520   // Rotates track parameters in R*phi plane
521   // if absolute rotation alpha is in global system
522   // otherwise alpha rotation is relative to the current rotation angle
523   //  
524
525   if (absolute) alpha -= GetAlpha();
526   //else fNRotate++;
527
528   return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
529 }
530
531 //___________________________________________________________
532 void AliTRDtrackV1::SetNumberOfClusters() 
533 {
534 // Calculate the number of clusters attached to this track
535         
536   Int_t ncls = 0;
537   for(int ip=0; ip<kNplane; ip++){
538     if(fTracklet[ip] && fTrackletIndex[ip] != 0xffff) ncls += fTracklet[ip]->GetN();
539   }
540   AliKalmanTrack::SetNumberOfClusters(ncls);    
541 }
542
543         
544 //_______________________________________________________________
545 void AliTRDtrackV1::SetOwner()
546 {
547   //
548   // Toggle ownership of tracklets
549   //
550
551   if(TestBit(1)) return;
552   for (Int_t ip = 0; ip < kNplane; ip++) {
553     if(fTrackletIndex[ip] == 0xffff) continue;
554     fTracklet[ip] = new AliTRDseedV1(*fTracklet[ip]);
555     fTracklet[ip]->SetOwner();
556   }
557   SetBit(1);
558 }
559
560 //_______________________________________________________________
561 void AliTRDtrackV1::SetTracklet(AliTRDseedV1 *trklt, Int_t index)
562 {
563   //
564   // Set the tracklets
565   //
566   Int_t plane = trklt->GetPlane();
567   fTracklet[plane]      = trklt;
568   fTrackletIndex[plane] = index;
569 }
570
571 //_______________________________________________________________
572 Bool_t  AliTRDtrackV1::Update(AliTRDseedV1 *trklt, Double_t chisq)
573 {
574   //
575   // Update track and tracklet parameters 
576   //
577   
578   Double_t x      = trklt->GetX0();
579   Double_t p[2]   = { trklt->GetYat(x)
580                     , trklt->GetZat(x) };
581   Double_t cov[3];
582   trklt->GetCovAt(x, cov);
583   
584   
585   if(!AliExternalTrackParam::Update(p, cov)) return kFALSE;
586   
587   AliTRDcluster *c = 0x0;
588   Int_t ic = 0; while(!(c = trklt->GetClusters(ic))) ic++;
589   AliTracker::FillResiduals(this, p, cov, c->GetVolumeId());
590   
591   
592   // Register info to track
593   SetNumberOfClusters();
594   SetChi2(GetChi2() + chisq);
595   
596   // update tracklet
597   trklt->SetMomentum(GetP());
598   trklt->SetSnp(GetSnp());
599   trklt->SetTgl(GetTgl());
600   return kTRUE;
601 }
602
603 //_______________________________________________________________
604 void AliTRDtrackV1::UpdateESDtrack(AliESDtrack *track)
605 {
606   //
607   // Update the TRD PID information in the ESD track
608   //
609
610   track->SetNumberOfTRDslices(kNslice);
611
612   for (Int_t ip = 0; ip < kNplane; ip++) {
613     if(fTrackletIndex[ip] == 0xffff) continue;
614     fTracklet[ip]->CookdEdx(kNslice);
615     Float_t *dedx = fTracklet[ip]->GetdEdx();
616     for (Int_t js = 0; js < kNslice; js++) track->SetTRDslice(dedx[js], ip, js);
617   }
618
619   // copy PID to ESD
620   track->SetTRDpid(fPID);
621   track->SetTRDpidQuality(fPIDquality);
622 }