226afd1f76e674a6435960ce8163207f23e5fd26
[u/mrichter/AliRoot.git] / STEER / AliESDtrack.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 //           Implementation of the ESD track class
17 //   ESD = Event Summary Data
18 //   This is the class to deal with during the phisics analysis of data
19 //      Origin: Iouri Belikov, CERN
20 //      e-mail: Jouri.Belikov@cern.ch
21 //-----------------------------------------------------------------
22
23 #include <TMath.h>
24 #include <TParticle.h>
25
26 #include "AliESDVertex.h"
27 #include "AliESDtrack.h"
28 #include "AliKalmanTrack.h"
29 #include "AliLog.h"
30 #include "AliTrackPointArray.h"
31
32 ClassImp(AliESDtrack)
33
34 void SetPIDValues(Float_t * dest, const Double_t * src, Int_t n) {
35   // This function copies "n" PID weights from "scr" to "dest"
36   // and normalizes their sum to 1 thus producing conditional probabilities.
37   // The negative weights are set to 0.
38   // In case all the weights are non-positive they are replaced by
39   // uniform probabilities
40
41   if (n<=0) return;
42
43   Float_t uniform = 1./(Float_t)n;
44
45   Float_t sum = 0;
46   for (Int_t i=0; i<n; i++) 
47     if (src[i]>=0) {
48       sum+=src[i];
49       dest[i] = src[i];
50     }
51     else {
52       dest[i] = 0;
53     }
54
55   if(sum>0)
56     for (Int_t i=0; i<n; i++) dest[i] /= sum;
57   else
58     for (Int_t i=0; i<n; i++) dest[i] = uniform;
59 }
60
61 //_______________________________________________________________________
62 AliESDtrack::AliESDtrack() : 
63   AliExternalTrackParam(),
64   fFlags(0),
65   fLabel(0),
66   fID(0),
67   fTrackLength(0),
68   fD(0),fZ(0),
69   fCdd(0),fCdz(0),fCzz(0),
70   fStopVertex(0),
71   fCp(0),
72   fCchi2(1e10),
73   fIp(0),
74   fTPCInner(0),
75   fOp(0),
76   fITSchi2(0),
77   fITSncls(0),
78   fITSClusterMap(0),
79   fITSsignal(0),
80   fITSLabel(0),
81   fTPCchi2(0),
82   fTPCncls(0),
83   fTPCnclsF(0),
84   fTPCClusterMap(159),//number of padrows
85   fTPCSharedMap(159),//number of padrows
86   fTPCsignal(0),
87   fTPCsignalN(0),
88   fTPCsignalS(0),
89   fTPCLabel(0),
90   fTRDchi2(0),
91   fTRDncls(0),
92   fTRDncls0(0),
93   fTRDsignal(0),
94   fTRDLabel(0),
95   fTRDQuality(0),
96   fTRDBudget(0),
97   fTOFchi2(0),
98   fTOFindex(0),
99   fTOFCalChannel(-1),
100   fTOFsignal(-1),
101   fTOFsignalToT(0),
102   fTOFsignalRaw(0),
103   fTOFsignalDz(0),
104   fHMPIDchi2(1e10),
105   fHMPIDqn(-1),
106   fHMPIDcluIdx(-1),
107   fHMPIDsignal(-1),
108   fHMPIDtrkTheta(-1),
109   fHMPIDtrkPhi(-1),
110   fHMPIDtrkX(-1),
111   fHMPIDtrkY(-1),
112   fHMPIDmipX(-1),
113   fHMPIDmipY(-1),
114   fEMCALindex(kEMCALNoMatch),
115   fFriendTrack(new AliESDfriendTrack())
116 {
117   //
118   // The default ESD constructor 
119   //
120   Int_t i, j;
121   for (i=0; i<AliPID::kSPECIES; i++) {
122     fTrackTime[i]=0.;
123     fR[i]=1.;
124     fITSr[i]=1.;
125     fTPCr[i]=1.;
126     fTRDr[i]=1.;
127     fTOFr[i]=1.;
128     fHMPIDr[i]=1.;
129   }
130   
131   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
132   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
133   for (i=0;i<kNPlane;i++) {
134     for (j=0;j<kNSlice;j++) {
135       fTRDsignals[i][j]=0.; 
136     }
137     fTRDTimBin[i]=-1;
138   }
139   for (i=0;i<4;i++) {fTPCPoints[i]=-1;}
140   for (i=0;i<3;i++) {fTOFLabel[i]=-1;}
141   for (i=0;i<10;i++) {fTOFInfo[i]=-1;}
142 }
143
144 //_______________________________________________________________________
145 AliESDtrack::AliESDtrack(const AliESDtrack& track):
146   AliExternalTrackParam(track),
147   fFlags(track.fFlags),
148   fLabel(track.fLabel),
149   fID(track.fID),
150   fTrackLength(track.fTrackLength),
151   fD(track.fD),fZ(track.fZ),
152   fCdd(track.fCdd),fCdz(track.fCdz),fCzz(track.fCzz),
153   fStopVertex(track.fStopVertex),
154   fCp(0),
155   fCchi2(track.fCchi2),
156   fIp(0),
157   fTPCInner(0),
158   fOp(0),
159   fITSchi2(track.fITSchi2),
160   fITSncls(track.fITSncls),
161   fITSClusterMap(track.fITSClusterMap),
162   fITSsignal(track.fITSsignal),
163   fITSLabel(track.fITSLabel),
164   fTPCchi2(track.fTPCchi2),
165   fTPCncls(track.fTPCncls),
166   fTPCnclsF(track.fTPCnclsF),
167   fTPCClusterMap(track.fTPCClusterMap),
168   fTPCSharedMap(track.fTPCSharedMap),
169   fTPCsignal(track.fTPCsignal),
170   fTPCsignalN(track.fTPCsignalN),
171   fTPCsignalS(track.fTPCsignalS),
172   fTPCLabel(track.fTPCLabel),
173   fTRDchi2(track.fTRDchi2),
174   fTRDncls(track.fTRDncls),
175   fTRDncls0(track.fTRDncls0),
176   fTRDsignal(track.fTRDsignal),
177   fTRDLabel(track.fTRDLabel),
178   fTRDQuality(track.fTRDQuality),
179   fTRDBudget(track.fTRDBudget),
180   fTOFchi2(track.fTOFchi2),
181   fTOFindex(track.fTOFindex),
182   fTOFCalChannel(track.fTOFCalChannel),
183   fTOFsignal(track.fTOFsignal),
184   fTOFsignalToT(track.fTOFsignalToT),
185   fTOFsignalRaw(track.fTOFsignalRaw),
186   fTOFsignalDz(track.fTOFsignalDz),
187   fHMPIDchi2(track.fHMPIDchi2),
188   fHMPIDqn(track.fHMPIDqn),
189   fHMPIDcluIdx(track.fHMPIDcluIdx),
190   fHMPIDsignal(track.fHMPIDsignal),
191   fHMPIDtrkTheta(track.fHMPIDtrkTheta),
192   fHMPIDtrkPhi(track.fHMPIDtrkPhi),
193   fHMPIDtrkX(track.fHMPIDtrkX),
194   fHMPIDtrkY(track.fHMPIDtrkY),
195   fHMPIDmipX(track.fHMPIDmipX),
196   fHMPIDmipY(track.fHMPIDmipY),
197   fEMCALindex(track.fEMCALindex),
198   fFriendTrack(0)
199 {
200   //
201   //copy constructor
202   //
203   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i]=track.fTrackTime[i];
204   for (Int_t i=0;i<AliPID::kSPECIES;i++)  fR[i]=track.fR[i];
205   //
206   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=track.fITSr[i]; 
207   //
208   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=track.fTPCr[i]; 
209   for (Int_t i=0;i<4;i++) {fTPCPoints[i]=track.fTPCPoints[i];}
210   for (Int_t i=0; i<3;i++)   { fKinkIndexes[i]=track.fKinkIndexes[i];}
211   for (Int_t i=0; i<3;i++)   { fV0Indexes[i]=track.fV0Indexes[i];}
212   //
213   for (Int_t i=0;i<kNPlane;i++) {
214     for (Int_t j=0;j<kNSlice;j++) {
215       fTRDsignals[i][j]=track.fTRDsignals[i][j]; 
216     }
217     fTRDTimBin[i]=track.fTRDTimBin[i];
218   }
219   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i]=track.fTRDr[i]; 
220   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i]=track.fTOFr[i];
221   for (Int_t i=0;i<3;i++) fTOFLabel[i]=track.fTOFLabel[i];
222   for (Int_t i=0;i<10;i++) fTOFInfo[i]=track.fTOFInfo[i];
223   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i]=track.fHMPIDr[i];
224
225   if (track.fCp) fCp=new AliExternalTrackParam(*track.fCp);
226   if (track.fIp) fIp=new AliExternalTrackParam(*track.fIp);
227   if (track.fTPCInner) fTPCInner=new AliExternalTrackParam(*track.fTPCInner);
228   if (track.fOp) fOp=new AliExternalTrackParam(*track.fOp);
229
230   if (track.fFriendTrack) fFriendTrack=new AliESDfriendTrack(*(track.fFriendTrack));
231 }
232
233 //_______________________________________________________________________
234 AliESDtrack::AliESDtrack(TParticle * part) : 
235   AliExternalTrackParam(),
236   fFlags(0),
237   fLabel(0),
238   fID(0),
239   fTrackLength(0),
240   fD(0),fZ(0),
241   fCdd(0),fCdz(0),fCzz(0),
242   fStopVertex(0),
243   fCp(0),
244   fCchi2(1e10),
245   fIp(0),
246   fTPCInner(0),
247   fOp(0),
248   fITSchi2(0),
249   fITSncls(0),
250   fITSClusterMap(0),
251   fITSsignal(0),
252   fITSLabel(0),
253   fTPCchi2(0),
254   fTPCncls(0),
255   fTPCnclsF(0),
256   fTPCClusterMap(159),//number of padrows
257   fTPCSharedMap(159),//number of padrows
258   fTPCsignal(0),
259   fTPCsignalN(0),
260   fTPCsignalS(0),
261   fTPCLabel(0),
262   fTRDchi2(0),
263   fTRDncls(0),
264   fTRDncls0(0),
265   fTRDsignal(0),
266   fTRDLabel(0),
267   fTRDQuality(0),
268   fTRDBudget(0),
269   fTOFchi2(0),
270   fTOFindex(0),
271   fTOFCalChannel(-1),
272   fTOFsignal(-1),
273   fTOFsignalToT(0),
274   fTOFsignalRaw(0),
275   fTOFsignalDz(0),
276   fHMPIDchi2(1e10),
277   fHMPIDqn(-1),
278   fHMPIDcluIdx(-1),
279   fHMPIDsignal(-1),
280   fHMPIDtrkTheta(-1),
281   fHMPIDtrkPhi(-1),
282   fHMPIDtrkX(-1),
283   fHMPIDtrkY(-1),
284   fHMPIDmipX(-1),
285   fHMPIDmipY(-1),
286   fEMCALindex(kEMCALNoMatch),
287   fFriendTrack(0)
288 {
289   //
290   // ESD track from TParticle
291   //
292
293   // Reset all the arrays
294   Int_t i, j;
295   for (i=0; i<AliPID::kSPECIES; i++) {
296     fTrackTime[i]=0.;
297     fR[i]=0.;
298     fITSr[i]=0.;
299     fTPCr[i]=0.;
300     fTRDr[i]=0.;
301     fTOFr[i]=0.;
302     fHMPIDr[i]=0.;
303   }
304   
305   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
306   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
307   for (i=0;i<kNPlane;i++) {
308     for (j=0;j<kNSlice;j++) {
309       fTRDsignals[i][j]=0.; 
310     }
311     fTRDTimBin[i]=-1;
312   }
313   for (i=0;i<4;i++) {fTPCPoints[i]=-1;}
314   for (i=0;i<3;i++) {fTOFLabel[i]=-1;}
315   for (i=0;i<10;i++) {fTOFInfo[i]=-1;}
316
317   // Calculate the AliExternalTrackParam content
318
319   Double_t xref;
320   Double_t alpha;
321   Double_t param[5];
322   Double_t covar[15];
323
324   // Calculate alpha: the rotation angle of the corresponding local system (TPC sector)
325   alpha = part->Phi()*180./TMath::Pi();
326   if (alpha<0) alpha+= 360.;
327   if (alpha>360) alpha -= 360.;
328
329   Int_t sector = (Int_t)(alpha/20.);
330   alpha = 10. + 20.*sector;
331   alpha /= 180;
332   alpha *= TMath::Pi();
333
334   // Covariance matrix: no errors, the parameters are exact
335   for (i=0; i<15; i++) covar[i]=0.;
336
337   // Get the vertex of origin and the momentum
338   TVector3 ver(part->Vx(),part->Vy(),part->Vz());
339   TVector3 mom(part->Px(),part->Py(),part->Pz());
340
341   // Rotate to the local coordinate system (TPC sector)
342   ver.RotateZ(-alpha);
343   mom.RotateZ(-alpha);
344
345   // X of the referense plane
346   xref = ver.X();
347
348   Int_t pdgCode = part->GetPdgCode();
349
350   Double_t charge = 
351     TDatabasePDG::Instance()->GetParticle(pdgCode)->Charge();
352
353   param[0] = ver.Y();
354   param[1] = ver.Z();
355   param[2] = TMath::Sin(mom.Phi());
356   param[3] = mom.Pz()/mom.Pt();
357   param[4] = TMath::Sign(1/mom.Pt(),charge);
358
359   // Set AliExternalTrackParam
360   Set(xref, alpha, param, covar);
361
362   // Set the PID
363   Int_t indexPID = 99;
364
365   switch (TMath::Abs(pdgCode)) {
366
367   case  11: // electron
368     indexPID = 0;
369     break;
370
371   case 13: // muon
372     indexPID = 1;
373     break;
374
375   case 211: // pion
376     indexPID = 2;
377     break;
378
379   case 321: // kaon
380     indexPID = 3;
381     break;
382
383   case 2212: // proton
384     indexPID = 4;
385     break;
386
387   default:
388     break;
389   }
390
391   // If the particle is not e,mu,pi,K or p the PID probabilities are set to 0
392   if (indexPID < AliPID::kSPECIES) {
393     fR[indexPID]=1.;
394     fITSr[indexPID]=1.;
395     fTPCr[indexPID]=1.;
396     fTRDr[indexPID]=1.;
397     fTOFr[indexPID]=1.;
398     fHMPIDr[indexPID]=1.;
399
400   }
401   // AliESD track label
402   SetLabel(part->GetUniqueID());
403
404 }
405
406 //_______________________________________________________________________
407 AliESDtrack::~AliESDtrack(){ 
408   //
409   // This is destructor according Coding Conventrions 
410   //
411   //printf("Delete track\n");
412   delete fIp; 
413   delete fTPCInner; 
414   delete fOp;
415   delete fCp; 
416   delete fFriendTrack;
417 }
418
419 void AliESDtrack::AddCalibObject(TObject * object){
420   //
421   // add calib object to the list
422   //
423   if (!fFriendTrack) fFriendTrack  = new AliESDfriendTrack;
424   fFriendTrack->AddCalibObject(object);
425 }
426
427 TObject *  AliESDtrack::GetCalibObject(Int_t index){
428   //
429   // return calib objct at given position
430   //
431   if (!fFriendTrack) return 0;
432   return fFriendTrack->GetCalibObject(index);
433 }
434
435
436 //_______________________________________________________________________
437 void AliESDtrack::MakeMiniESDtrack(){
438   // Resets everything except
439   // fFlags: Reconstruction status flags 
440   // fLabel: Track label
441   // fID:  Unique ID of the track
442   // fD: Impact parameter in XY-plane
443   // fZ: Impact parameter in Z 
444   // fR[AliPID::kSPECIES]: combined "detector response probability"
445   // Running track parameters in the base class (AliExternalTrackParam)
446   
447   fTrackLength = 0;
448   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i] = 0;
449   fStopVertex = 0;
450
451   // Reset track parameters constrained to the primary vertex
452   fCp = 0;
453   fCchi2 = 0;
454
455   // Reset track parameters at the inner wall of TPC
456   fIp = 0;
457   fTPCInner=0;
458   // Reset track parameters at the inner wall of the TRD
459   fOp = 0;
460
461   // Reset ITS track related information
462   fITSchi2 = 0;
463   fITSncls = 0;       
464   fITSClusterMap=0;
465   fITSsignal = 0;     
466   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=0; 
467   fITSLabel = 0;       
468
469   // Reset TPC related track information
470   fTPCchi2 = 0;       
471   fTPCncls = 0;       
472   fTPCnclsF = 0;       
473   fTPCClusterMap = 0;  
474   fTPCSharedMap = 0;  
475   fTPCsignal= 0;      
476   fTPCsignalS= 0;      
477   fTPCsignalN= 0;      
478   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0; 
479   fTPCLabel=0;       
480   for (Int_t i=0;i<4;i++) fTPCPoints[i] = 0;
481   for (Int_t i=0; i<3;i++)   fKinkIndexes[i] = 0;
482   for (Int_t i=0; i<3;i++)   fV0Indexes[i] = 0;
483
484   // Reset TRD related track information
485   fTRDchi2 = 0;        
486   fTRDncls = 0;       
487   fTRDncls0 = 0;       
488   fTRDsignal = 0;      
489   for (Int_t i=0;i<kNPlane;i++) {
490     for (Int_t j=0;j<kNSlice;j++) {
491       fTRDsignals[i][j] = 0; 
492     }
493     fTRDTimBin[i]  = 0;
494   }
495   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0; 
496   fTRDLabel = 0;       
497   fTRDQuality  = 0;
498   fTRDBudget  = 0;
499
500   // Reset TOF related track information
501   fTOFchi2 = 0;        
502   fTOFindex = 0;       
503   fTOFsignal = 0;      
504   fTOFCalChannel = -1;
505   fTOFsignalToT = 0;
506   fTOFsignalRaw = 0;
507   fTOFsignalDz = 0;
508   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i] = 0;
509   for (Int_t i=0;i<3;i++) fTOFLabel[i] = 0;
510   for (Int_t i=0;i<10;i++) fTOFInfo[i] = 0;
511
512   // Reset HMPID related track information
513   fHMPIDchi2 = 0;     
514   fHMPIDqn = -1;     
515   fHMPIDcluIdx = -1;     
516   fHMPIDsignal = 0;     
517   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i] = 0;
518   fHMPIDtrkTheta = -1;     
519   fHMPIDtrkPhi = -1;      
520   fHMPIDtrkX = -1;     
521   fHMPIDtrkY = -1;      
522   fHMPIDmipX = -1;
523   fHMPIDmipY = -1;
524   fEMCALindex = kEMCALNoMatch;
525
526   delete fFriendTrack; fFriendTrack = 0;
527
528 //_______________________________________________________________________
529 Double_t AliESDtrack::GetMass() const {
530   // Returns the mass of the most probable particle type
531   Float_t max=0.;
532   Int_t k=-1;
533   for (Int_t i=0; i<AliPID::kSPECIES; i++) {
534     if (fR[i]>max) {k=i; max=fR[i];}
535   }
536   if (k==0) { // dE/dx "crossing points" in the TPC
537      Double_t p=GetP();
538      if ((p>0.38)&&(p<0.48))
539         if (fR[0]<fR[3]*10.) return AliPID::ParticleMass(AliPID::kKaon);
540      if ((p>0.75)&&(p<0.85))
541         if (fR[0]<fR[4]*10.) return AliPID::ParticleMass(AliPID::kProton);
542      return 0.00051;
543   }
544   if (k==1) return AliPID::ParticleMass(AliPID::kMuon); 
545   if (k==2||k==-1) return AliPID::ParticleMass(AliPID::kPion);
546   if (k==3) return AliPID::ParticleMass(AliPID::kKaon);
547   if (k==4) return AliPID::ParticleMass(AliPID::kProton);
548   AliWarning("Undefined mass !");
549   return AliPID::ParticleMass(AliPID::kPion);
550 }
551
552 //______________________________________________________________________________
553 Double_t AliESDtrack::E() const
554 {
555   // Returns the energy of the particle given its assumed mass.
556   // Assumes the pion mass if the particle can't be identified properly.
557   
558   Double_t m = M();
559   Double_t p = P();
560   return TMath::Sqrt(p*p + m*m);
561 }
562
563 //______________________________________________________________________________
564 Double_t AliESDtrack::Y() const
565 {
566   // Returns the rapidity of a particle given its assumed mass.
567   // Assumes the pion mass if the particle can't be identified properly.
568   
569   Double_t e = E();
570   Double_t pz = Pz();
571   if (e != TMath::Abs(pz)) { // energy was not equal to pz
572     return 0.5*TMath::Log((e+pz)/(e-pz));
573   } else { // energy was equal to pz
574     return -999.;
575   }
576 }
577
578 //_______________________________________________________________________
579 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
580   //
581   // This function updates track's running parameters 
582   //
583   Int_t *index=0;
584   Bool_t rc=kTRUE;
585
586   SetStatus(flags);
587   fLabel=t->GetLabel();
588
589   if (t->IsStartedTimeIntegral()) {
590     SetStatus(kTIME);
591     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
592     SetIntegratedLength(t->GetIntegratedLength());
593   }
594
595   Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
596   
597   switch (flags) {
598     
599   case kITSin: case kITSout: case kITSrefit:
600     fITSncls=t->GetNumberOfClusters();
601     index=fFriendTrack->GetITSindices(); 
602     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
603         index[i]=t->GetClusterIndex(i);
604         if (i<fITSncls) {
605            Int_t l=(index[i] & 0xf0000000) >> 28;
606            SETBIT(fITSClusterMap,l);                 
607         }
608     }
609     fITSchi2=t->GetChi2();
610     fITSsignal=t->GetPIDsignal();
611     fITSLabel = t->GetLabel();
612     break;
613     
614   case kTPCin: case kTPCrefit:
615     fTPCLabel = t->GetLabel();
616     if (flags==kTPCin)  fTPCInner=new AliExternalTrackParam(*t);
617     if (!fIp) fIp=new AliExternalTrackParam(*t);
618     else 
619       fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
620   case kTPCout:
621     index=fFriendTrack->GetTPCindices(); 
622     if (flags & kTPCout){
623       if (!fOp) fOp=new AliExternalTrackParam(*t);
624       else 
625         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
626     }
627     fTPCncls=t->GetNumberOfClusters();    
628     fTPCchi2=t->GetChi2();
629     
630      {//prevrow must be declared in separate namespace, otherwise compiler cries:
631       //"jump to case label crosses initialization of `Int_t prevrow'"
632        Int_t prevrow = -1;
633        //       for (Int_t i=0;i<fTPCncls;i++) 
634        for (Int_t i=0;i<AliESDfriendTrack::kMaxTPCcluster;i++) 
635         {
636           index[i]=t->GetClusterIndex(i);
637           Int_t idx = index[i];
638
639           if (idx<0) continue; 
640
641           // Piotr's Cluster Map for HBT  
642           // ### please change accordingly if cluster array is changing 
643           // to "New TPC Tracking" style (with gaps in array) 
644           Int_t sect = (idx&0xff000000)>>24;
645           Int_t row = (idx&0x00ff0000)>>16;
646           if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
647
648           fTPCClusterMap.SetBitNumber(row,kTRUE);
649
650           //Fill the gap between previous row and this row with 0 bits
651           //In case  ###  pleas change it as well - just set bit 0 in case there 
652           //is no associated clusters for current "i"
653           if (prevrow < 0) 
654            {
655              prevrow = row;//if previous bit was not assigned yet == this is the first one
656            }
657           else
658            { //we don't know the order (inner to outer or reverse)
659              //just to be save in case it is going to change
660              Int_t n = 0, m = 0;
661              if (prevrow < row)
662               {
663                 n = prevrow;
664                 m = row;
665               }
666              else
667               {
668                 n = row;
669                 m = prevrow;
670               }
671
672              for (Int_t j = n+1; j < m; j++)
673               {
674                 fTPCClusterMap.SetBitNumber(j,kFALSE);
675               }
676              prevrow = row; 
677            }
678           // End Of Piotr's Cluster Map for HBT
679         }
680      }
681     fTPCsignal=t->GetPIDsignal();
682     break;
683
684   case kTRDout: case kTRDin: case kTRDrefit:
685     index=fFriendTrack->GetTRDindices();
686     fTRDLabel = t->GetLabel(); 
687     fTRDncls=t->GetNumberOfClusters();
688     fTRDchi2=t->GetChi2();
689     for (Int_t i=0;i<fTRDncls;i++) index[i]=t->GetClusterIndex(i);
690     fTRDsignal=t->GetPIDsignal();
691     break;
692   case kTRDbackup:
693     if (!fOp) fOp=new AliExternalTrackParam(*t);
694     else 
695       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
696     fTRDncls0 = t->GetNumberOfClusters(); 
697     break;
698   case kTOFin: 
699     break;
700   case kTOFout: 
701     break;
702   case kTRDStop:
703     break;
704   default: 
705     AliError("Wrong flag !");
706     return kFALSE;
707   }
708
709   return rc;
710 }
711
712 //_______________________________________________________________________
713 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
714   //---------------------------------------------------------------------
715   // This function returns external representation of the track parameters
716   //---------------------------------------------------------------------
717   x=GetX();
718   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
719 }
720
721 //_______________________________________________________________________
722 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
723   //---------------------------------------------------------------------
724   // This function returns external representation of the cov. matrix
725   //---------------------------------------------------------------------
726   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
727 }
728
729 //_______________________________________________________________________
730 Bool_t AliESDtrack::GetConstrainedExternalParameters
731                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
732   //---------------------------------------------------------------------
733   // This function returns the constrained external track parameters
734   //---------------------------------------------------------------------
735   if (!fCp) return kFALSE;
736   alpha=fCp->GetAlpha();
737   x=fCp->GetX();
738   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
739   return kTRUE;
740 }
741
742 //_______________________________________________________________________
743 Bool_t 
744 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
745   //---------------------------------------------------------------------
746   // This function returns the constrained external cov. matrix
747   //---------------------------------------------------------------------
748   if (!fCp) return kFALSE;
749   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
750   return kTRUE;
751 }
752
753 Bool_t
754 AliESDtrack::GetInnerExternalParameters
755                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
756   //---------------------------------------------------------------------
757   // This function returns external representation of the track parameters 
758   // at the inner layer of TPC
759   //---------------------------------------------------------------------
760   if (!fIp) return kFALSE;
761   alpha=fIp->GetAlpha();
762   x=fIp->GetX();
763   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
764   return kTRUE;
765 }
766
767 Bool_t 
768 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
769  //---------------------------------------------------------------------
770  // This function returns external representation of the cov. matrix 
771  // at the inner layer of TPC
772  //---------------------------------------------------------------------
773   if (!fIp) return kFALSE;
774   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
775   return kTRUE;
776 }
777
778 Bool_t 
779 AliESDtrack::GetOuterExternalParameters
780                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
781   //---------------------------------------------------------------------
782   // This function returns external representation of the track parameters 
783   // at the inner layer of TRD
784   //---------------------------------------------------------------------
785   if (!fOp) return kFALSE;
786   alpha=fOp->GetAlpha();
787   x=fOp->GetX();
788   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
789   return kTRUE;
790 }
791
792 Bool_t 
793 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
794  //---------------------------------------------------------------------
795  // This function returns external representation of the cov. matrix 
796  // at the inner layer of TRD
797  //---------------------------------------------------------------------
798   if (!fOp) return kFALSE;
799   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
800   return kTRUE;
801 }
802
803 Int_t AliESDtrack::GetNcls(Int_t idet) const
804 {
805   // Get number of clusters by subdetector index
806   //
807   Int_t ncls = 0;
808   switch(idet){
809   case 0:
810     ncls = fITSncls;
811     break;
812   case 1:
813     ncls = fTPCncls;
814     break;
815   case 2:
816     ncls = fTRDncls;
817     break;
818   case 3:
819     if (fTOFindex != 0)
820       ncls = 1;
821     break;
822   default:
823     break;
824   }
825   return ncls;
826 }
827
828 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
829 {
830   // Get cluster index array by subdetector index
831   //
832   Int_t ncls = 0;
833   switch(idet){
834   case 0:
835     ncls = GetITSclusters(idx);
836     break;
837   case 1:
838     ncls = GetTPCclusters(idx);
839     break;
840   case 2:
841     ncls = GetTRDclusters(idx);
842     break;
843   case 3:
844     if (fTOFindex != 0) {
845       idx[0] = GetTOFcluster();
846       ncls = 1;
847     }
848     break;
849   default:
850     break;
851   }
852   return ncls;
853 }
854
855 //_______________________________________________________________________
856 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
857   // Returns the array with integrated times for each particle hypothesis
858   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
859 }
860
861 //_______________________________________________________________________
862 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
863   // Sets the array with integrated times for each particle hypotesis
864   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
865 }
866
867 //_______________________________________________________________________
868 void AliESDtrack::SetITSpid(const Double_t *p) {
869   // Sets values for the probability of each particle type (in ITS)
870   SetPIDValues(fITSr,p,AliPID::kSPECIES);
871   SetStatus(AliESDtrack::kITSpid);
872 }
873
874 //_______________________________________________________________________
875 void AliESDtrack::GetITSpid(Double_t *p) const {
876   // Gets the probability of each particle type (in ITS)
877   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
878 }
879
880 //_______________________________________________________________________
881 Int_t AliESDtrack::GetITSclusters(Int_t *idx) const {
882   //---------------------------------------------------------------------
883   // This function returns indices of the assgined ITS clusters 
884   //---------------------------------------------------------------------
885   if (idx!=0) {
886      Int_t *index=fFriendTrack->GetITSindices();
887      for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) idx[i]=index[i];
888   }
889   return fITSncls;
890 }
891
892 //_______________________________________________________________________
893 Int_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
894   //---------------------------------------------------------------------
895   // This function returns indices of the assgined ITS clusters 
896   //---------------------------------------------------------------------
897   if (idx!=0) {
898     Int_t *index=fFriendTrack->GetTPCindices();
899     for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
900   }
901   return fTPCncls;
902 }
903
904 Float_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
905   //
906   // GetDensity of the clusters on given region between row0 and row1
907   // Dead zone effect takin into acoount
908   //
909   Int_t good  = 0;
910   Int_t found = 0;
911   //  
912   Int_t *index=fFriendTrack->GetTPCindices();
913   for (Int_t i=row0;i<=row1;i++){     
914     Int_t idx = index[i];
915     if (idx!=-1)  good++;             // track outside of dead zone
916     if (idx>0)    found++;
917   }
918   Float_t density=0.5;
919   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
920   return density;
921 }
922
923 //_______________________________________________________________________
924 void AliESDtrack::SetTPCpid(const Double_t *p) {  
925   // Sets values for the probability of each particle type (in TPC)
926   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
927   SetStatus(AliESDtrack::kTPCpid);
928 }
929
930 //_______________________________________________________________________
931 void AliESDtrack::GetTPCpid(Double_t *p) const {
932   // Gets the probability of each particle type (in TPC)
933   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
934 }
935
936 //_______________________________________________________________________
937 Int_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
938   //---------------------------------------------------------------------
939   // This function returns indices of the assgined TRD clusters 
940   //---------------------------------------------------------------------
941   if (idx!=0) {
942      Int_t *index=fFriendTrack->GetTRDindices();
943      for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
944   }
945   return fTRDncls;
946 }
947
948 //_______________________________________________________________________
949 void AliESDtrack::SetTRDpid(const Double_t *p) {  
950   // Sets values for the probability of each particle type (in TRD)
951   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
952   SetStatus(AliESDtrack::kTRDpid);
953 }
954
955 //_______________________________________________________________________
956 void AliESDtrack::GetTRDpid(Double_t *p) const {
957   // Gets the probability of each particle type (in TRD)
958   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
959 }
960
961 //_______________________________________________________________________
962 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
963 {
964   // Sets the probability of particle type iSpecies to p (in TRD)
965   fTRDr[iSpecies] = p;
966 }
967
968 Float_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
969 {
970   // Returns the probability of particle type iSpecies (in TRD)
971   return fTRDr[iSpecies];
972 }
973
974 //_______________________________________________________________________
975 void AliESDtrack::SetTOFpid(const Double_t *p) {  
976   // Sets the probability of each particle type (in TOF)
977   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
978   SetStatus(AliESDtrack::kTOFpid);
979 }
980
981 //_______________________________________________________________________
982 void AliESDtrack::SetTOFLabel(const Int_t *p) {  
983   // Sets  (in TOF)
984   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
985 }
986
987 //_______________________________________________________________________
988 void AliESDtrack::GetTOFpid(Double_t *p) const {
989   // Gets probabilities of each particle type (in TOF)
990   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
991 }
992
993 //_______________________________________________________________________
994 void AliESDtrack::GetTOFLabel(Int_t *p) const {
995   // Gets (in TOF)
996   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
997 }
998
999 //_______________________________________________________________________
1000 void AliESDtrack::GetTOFInfo(Float_t *info) const {
1001   // Gets (in TOF)
1002   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
1003 }
1004
1005 //_______________________________________________________________________
1006 void AliESDtrack::SetTOFInfo(Float_t*info) {
1007   // Gets (in TOF)
1008   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
1009 }
1010
1011
1012
1013 //_______________________________________________________________________
1014 void AliESDtrack::SetHMPIDpid(const Double_t *p) {  
1015   // Sets the probability of each particle type (in HMPID)
1016   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
1017   SetStatus(AliESDtrack::kHMPIDpid);
1018 }
1019
1020 //_______________________________________________________________________
1021 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
1022   // Gets probabilities of each particle type (in HMPID)
1023   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
1024 }
1025
1026
1027
1028 //_______________________________________________________________________
1029 void AliESDtrack::SetESDpid(const Double_t *p) {  
1030   // Sets the probability of each particle type for the ESD track
1031   SetPIDValues(fR,p,AliPID::kSPECIES);
1032   SetStatus(AliESDtrack::kESDpid);
1033 }
1034
1035 //_______________________________________________________________________
1036 void AliESDtrack::GetESDpid(Double_t *p) const {
1037   // Gets probability of each particle type for the ESD track
1038   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
1039 }
1040
1041 //_______________________________________________________________________
1042 Bool_t AliESDtrack::RelateToVertex
1043 (const AliESDVertex *vtx, Double_t b, Double_t maxd) {
1044   //
1045   // Try to relate this track to the vertex "vtx", 
1046   // if the (rough) transverse impact parameter is not bigger then "maxd". 
1047   //            Magnetic field is "b" (kG).
1048   //
1049   // a) The track gets extapolated to the DCA to the vertex.
1050   // b) The impact parameters and their covariance matrix are calculated.
1051   // c) An attempt to constrain this track to the vertex is done.
1052   //
1053   //    In the case of success, the returned value is kTRUE
1054   //    (otherwise, it's kFALSE)
1055   //  
1056
1057   if (!vtx) return kFALSE;
1058
1059   Double_t alpha=GetAlpha();
1060   Double_t sn=TMath::Sin(alpha), cs=TMath::Cos(alpha);
1061   Double_t x=GetX(), y=GetParameter()[0], snp=GetParameter()[2];
1062   Double_t xv= vtx->GetXv()*cs + vtx->GetYv()*sn;
1063   Double_t yv=-vtx->GetXv()*sn + vtx->GetYv()*cs, zv=vtx->GetZv();
1064   x-=xv; y-=yv;
1065
1066   //Estimate the impact parameter neglecting the track curvature
1067   Double_t d=TMath::Abs(x*snp - y*TMath::Sqrt(1.- snp*snp));
1068   if (d > maxd) return kFALSE; 
1069
1070   //Propagate to the DCA
1071   Double_t crv=kB2C*b*GetParameter()[4];
1072   if (TMath::Abs(b) < kAlmost0Field) crv=0.;
1073
1074   Double_t tgfv=-(crv*x - snp)/(crv*y + TMath::Sqrt(1.-snp*snp));
1075   sn=tgfv/TMath::Sqrt(1.+ tgfv*tgfv);
1076   if (TMath::Abs(tgfv)>0.) cs = sn/tgfv;
1077   else cs=1.;
1078
1079   x = xv*cs + yv*sn;
1080   yv=-xv*sn + yv*cs; xv=x;
1081
1082   if (!Propagate(alpha+TMath::ASin(sn),xv,b)) return kFALSE;
1083
1084   fD = GetParameter()[0] - yv;
1085   fZ = GetParameter()[1] - zv;
1086   
1087   Double_t cov[6]; vtx->GetCovMatrix(cov);
1088
1089   //***** Improvements by A.Dainese
1090   alpha=GetAlpha(); sn=TMath::Sin(alpha); cs=TMath::Cos(alpha);
1091   Double_t s2ylocvtx = cov[0]*sn*sn + cov[2]*cs*cs - 2.*cov[1]*cs*sn;
1092   fCdd = GetCovariance()[0] + s2ylocvtx;   // neglecting correlations
1093   fCdz = GetCovariance()[1];               // between (x,y) and z
1094   fCzz = GetCovariance()[2] + cov[5];      // in vertex's covariance matrix
1095   //*****
1096
1097   {//Try to constrain 
1098     Double_t p[2]={yv,zv}, c[3]={cov[2],0.,cov[5]};
1099     Double_t chi2=GetPredictedChi2(p,c);
1100
1101     if (chi2>77.) return kFALSE;
1102
1103     AliExternalTrackParam tmp(*this);
1104     if (!tmp.Update(p,c)) return kFALSE;
1105
1106     fCchi2=chi2;
1107     if (!fCp) fCp=new AliExternalTrackParam();
1108     new (fCp) AliExternalTrackParam(tmp);
1109   }
1110
1111   return kTRUE;
1112 }
1113
1114 //_______________________________________________________________________
1115 void AliESDtrack::Print(Option_t *) const {
1116   // Prints info on the track
1117   
1118   printf("ESD track info\n") ; 
1119   Double_t p[AliPID::kSPECIESN] ; 
1120   Int_t index = 0 ; 
1121   if( IsOn(kITSpid) ){
1122     printf("From ITS: ") ; 
1123     GetITSpid(p) ; 
1124     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1125       printf("%f, ", p[index]) ;
1126     printf("\n           signal = %f\n", GetITSsignal()) ;
1127   } 
1128   if( IsOn(kTPCpid) ){
1129     printf("From TPC: ") ; 
1130     GetTPCpid(p) ; 
1131     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1132       printf("%f, ", p[index]) ;
1133     printf("\n           signal = %f\n", GetTPCsignal()) ;
1134   }
1135   if( IsOn(kTRDpid) ){
1136     printf("From TRD: ") ; 
1137     GetTRDpid(p) ; 
1138     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1139       printf("%f, ", p[index]) ;
1140     printf("\n           signal = %f\n", GetTRDsignal()) ;
1141   }
1142   if( IsOn(kTOFpid) ){
1143     printf("From TOF: ") ; 
1144     GetTOFpid(p) ; 
1145     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1146       printf("%f, ", p[index]) ;
1147     printf("\n           signal = %f\n", GetTOFsignal()) ;
1148   }
1149   if( IsOn(kHMPIDpid) ){
1150     printf("From HMPID: ") ; 
1151     GetHMPIDpid(p) ; 
1152     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1153       printf("%f, ", p[index]) ;
1154     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
1155   }
1156
1157