Change from Int_t/SHort_t to smaller data types, the change from Float_t/Double_t...
[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(Double_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   fCp(0),
65   fIp(0),
66   fTPCInner(0),
67   fOp(0),
68   fFriendTrack(new AliESDfriendTrack()),
69   fTPCClusterMap(159),//number of padrows
70   fTPCSharedMap(159),//number of padrows
71   fFlags(0),
72   fID(0),
73   fLabel(0),
74   fITSLabel(0),
75   fTPCLabel(0),
76   fTRDLabel(0),
77   fTOFCalChannel(0),
78   fTOFindex(0),
79   fHMPIDqn(0),
80   fHMPIDcluIdx(0),
81   fEMCALindex(kEMCALNoMatch),
82   fHMPIDtrkTheta(0),
83   fHMPIDtrkPhi(0),
84   fHMPIDsignal(0),
85   fTrackLength(0),
86   fD(0),fZ(0),
87   fCdd(0),fCdz(0),fCzz(0),
88   fCchi2(0),
89   fITSchi2(0),
90   fTPCchi2(0),
91   fTRDchi2(0),
92   fTOFchi2(0),
93   fHMPIDchi2(0),
94   fITSsignal(0),
95   fTPCsignal(0),
96   fTPCsignalS(0),
97   fTRDsignal(0),
98   fTRDQuality(0),
99   fTRDBudget(0),
100   fTOFsignal(0),
101   fTOFsignalToT(0),
102   fTOFsignalRaw(0),
103   fTOFsignalDz(0),
104   fHMPIDtrkX(0),
105   fHMPIDtrkY(0),
106   fHMPIDmipX(0),
107   fHMPIDmipY(0),
108   fTPCncls(0),
109   fTPCnclsF(0),
110   fTPCsignalN(0),
111   fITSncls(0),
112   fITSClusterMap(0),
113   fTRDncls(0),
114   fTRDncls0(0),
115   fTRDpidQuality(0)
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]=0.;
124     fITSr[i]=0.;
125     fTPCr[i]=0.;
126     fTRDr[i]=0.;
127     fTOFr[i]=0.;
128     fHMPIDr[i]=0.;
129   }
130   
131   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
132   for (i=0; i<3; i++)   { fV0Indexes[i]=0;}
133   for (i=0;i<kNPlane;i++) {
134     for (j=0;j<kNSlice;j++) {
135       fTRDsignals[i][j]=0.; 
136     }
137     fTRDTimBin[i]=0;
138   }
139   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
140   for (i=0;i<3;i++) {fTOFLabel[i]=0;}
141   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
142 }
143
144 //_______________________________________________________________________
145 AliESDtrack::AliESDtrack(const AliESDtrack& track):
146   AliExternalTrackParam(track),
147   fCp(0),
148   fIp(0),
149   fTPCInner(0),
150   fOp(0),
151   fFriendTrack(0),
152   fTPCClusterMap(track.fTPCClusterMap),
153   fTPCSharedMap(track.fTPCSharedMap),
154   fFlags(track.fFlags),
155   fID(track.fID),
156   fLabel(track.fLabel),
157   fITSLabel(track.fITSLabel),
158   fTPCLabel(track.fTPCLabel),
159   fTRDLabel(track.fTRDLabel),
160   fTOFCalChannel(track.fTOFCalChannel),
161   fTOFindex(track.fTOFindex),
162   fHMPIDqn(track.fHMPIDqn),
163   fHMPIDcluIdx(track.fHMPIDcluIdx),
164   fEMCALindex(track.fEMCALindex),
165   fHMPIDtrkTheta(track.fHMPIDtrkTheta),
166   fHMPIDtrkPhi(track.fHMPIDtrkPhi),
167   fHMPIDsignal(track.fHMPIDsignal),
168   fTrackLength(track.fTrackLength),
169   fD(track.fD),fZ(track.fZ),
170   fCdd(track.fCdd),fCdz(track.fCdz),fCzz(track.fCzz),
171   fCchi2(track.fCchi2),
172   fITSchi2(track.fITSchi2),
173   fTPCchi2(track.fTPCchi2),
174   fTRDchi2(track.fTRDchi2),
175   fTOFchi2(track.fTOFchi2),
176   fHMPIDchi2(track.fHMPIDchi2),
177   fITSsignal(track.fITSsignal),
178   fTPCsignal(track.fTPCsignal),
179   fTPCsignalS(track.fTPCsignalS),
180   fTRDsignal(track.fTRDsignal),
181   fTRDQuality(track.fTRDQuality),
182   fTRDBudget(track.fTRDBudget),
183   fTOFsignal(track.fTOFsignal),
184   fTOFsignalToT(track.fTOFsignalToT),
185   fTOFsignalRaw(track.fTOFsignalRaw),
186   fTOFsignalDz(track.fTOFsignalDz),
187   fHMPIDtrkX(track.fHMPIDtrkX),
188   fHMPIDtrkY(track.fHMPIDtrkY),
189   fHMPIDmipX(track.fHMPIDmipX),
190   fHMPIDmipY(track.fHMPIDmipY),
191   fTPCncls(track.fTPCncls),
192   fTPCnclsF(track.fTPCnclsF),
193   fTPCsignalN(track.fTPCsignalN),
194   fITSncls(track.fITSncls),
195   fITSClusterMap(track.fITSClusterMap),
196   fTRDncls(track.fTRDncls),
197   fTRDncls0(track.fTRDncls0),
198   fTRDpidQuality(track.fTRDpidQuality)
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   fCp(0),
237   fIp(0),
238   fTPCInner(0),
239   fOp(0),
240   fFriendTrack(0),
241   fTPCClusterMap(159),//number of padrows
242   fTPCSharedMap(159),//number of padrows
243   fFlags(0),
244   fID(0),
245   fLabel(0),
246   fITSLabel(0),
247   fTPCLabel(0),
248   fTRDLabel(0),
249   fTOFCalChannel(0),
250   fTOFindex(0),
251   fHMPIDqn(0),
252   fHMPIDcluIdx(0),
253   fEMCALindex(kEMCALNoMatch),
254   fHMPIDtrkTheta(0),
255   fHMPIDtrkPhi(0),
256   fHMPIDsignal(0),
257   fTrackLength(0),
258   fD(0),fZ(0),
259   fCdd(0),fCdz(0),fCzz(0),
260   fCchi2(0),
261   fITSchi2(0),
262   fTPCchi2(0),
263   fTRDchi2(0),
264   fTOFchi2(0),
265   fHMPIDchi2(0),
266   fITSsignal(0),
267   fTPCsignal(0),
268   fTPCsignalS(0),
269   fTRDsignal(0),
270   fTRDQuality(0),
271   fTRDBudget(0),
272   fTOFsignal(0),
273   fTOFsignalToT(0),
274   fTOFsignalRaw(0),
275   fTOFsignalDz(0),
276   fHMPIDtrkX(0),
277   fHMPIDtrkY(0),
278   fHMPIDmipX(0),
279   fHMPIDmipY(0),
280   fTPCncls(0),
281   fTPCnclsF(0),
282   fTPCsignalN(0),
283   fITSncls(0),
284   fITSClusterMap(0),
285   fTRDncls(0),
286   fTRDncls0(0),
287   fTRDpidQuality(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]=0;
312   }
313   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
314   for (i=0;i<3;i++) {fTOFLabel[i]=0;}
315   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
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
449   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i] = 0;
450
451   // Reset track parameters constrained to the primary vertex
452   delete fCp;fCp = 0;
453   fCchi2 = 0;
454
455   // Reset track parameters at the inner wall of TPC
456   delete fIp;fIp = 0;
457   delete fTPCInner;fTPCInner=0;
458   // Reset track parameters at the inner wall of the TRD
459   delete fOp;fOp = 0;
460
461
462   // Reset ITS track related information
463   fITSchi2 = 0;
464   fITSncls = 0;       
465   fITSClusterMap=0;
466   fITSsignal = 0;     
467   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=0; 
468   fITSLabel = 0;       
469
470   // Reset TPC related track information
471   fTPCchi2 = 0;       
472   fTPCncls = 0;       
473   fTPCnclsF = 0;       
474   fTPCClusterMap = 0;  
475   fTPCSharedMap = 0;  
476   fTPCsignal= 0;      
477   fTPCsignalS= 0;      
478   fTPCsignalN= 0;      
479   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0; 
480   fTPCLabel=0;       
481   for (Int_t i=0;i<4;i++) fTPCPoints[i] = 0;
482   for (Int_t i=0; i<3;i++)   fKinkIndexes[i] = 0;
483   for (Int_t i=0; i<3;i++)   fV0Indexes[i] = 0;
484
485   // Reset TRD related track information
486   fTRDchi2 = 0;        
487   fTRDncls = 0;       
488   fTRDncls0 = 0;       
489   fTRDsignal = 0;      
490   for (Int_t i=0;i<kNPlane;i++) {
491     for (Int_t j=0;j<kNSlice;j++) {
492       fTRDsignals[i][j] = 0; 
493     }
494     fTRDTimBin[i]  = 0;
495   }
496   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0; 
497   fTRDLabel = 0;       
498   fTRDQuality  = 0;
499   fTRDpidQuality = 0;
500   fTRDBudget  = 0;
501
502   // Reset TOF related track information
503   fTOFchi2 = 0;        
504   fTOFindex = 0;       
505   fTOFsignal = 0;      
506   fTOFCalChannel = 0;
507   fTOFsignalToT = 0;
508   fTOFsignalRaw = 0;
509   fTOFsignalDz = 0;
510   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i] = 0;
511   for (Int_t i=0;i<3;i++) fTOFLabel[i] = 0;
512   for (Int_t i=0;i<10;i++) fTOFInfo[i] = 0;
513
514   // Reset HMPID related track information
515   fHMPIDchi2 = 0;     
516   fHMPIDqn = 0;     
517   fHMPIDcluIdx = 0;     
518   fHMPIDsignal = 0;     
519   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i] = 0;
520   fHMPIDtrkTheta = 0;     
521   fHMPIDtrkPhi = 0;      
522   fHMPIDtrkX = 0;     
523   fHMPIDtrkY = 0;      
524   fHMPIDmipX = 0;
525   fHMPIDmipY = 0;
526   fEMCALindex = kEMCALNoMatch;
527
528   delete fFriendTrack; fFriendTrack = 0;
529
530 //_______________________________________________________________________
531 Double_t AliESDtrack::GetMass() const {
532   // Returns the mass of the most probable particle type
533   Float_t max=0.;
534   Int_t k=-1;
535   for (Int_t i=0; i<AliPID::kSPECIES; i++) {
536     if (fR[i]>max) {k=i; max=fR[i];}
537   }
538   if (k==0) { // dE/dx "crossing points" in the TPC
539      Double_t p=GetP();
540      if ((p>0.38)&&(p<0.48))
541         if (fR[0]<fR[3]*10.) return AliPID::ParticleMass(AliPID::kKaon);
542      if ((p>0.75)&&(p<0.85))
543         if (fR[0]<fR[4]*10.) return AliPID::ParticleMass(AliPID::kProton);
544      return 0.00051;
545   }
546   if (k==1) return AliPID::ParticleMass(AliPID::kMuon); 
547   if (k==2||k==-1) return AliPID::ParticleMass(AliPID::kPion);
548   if (k==3) return AliPID::ParticleMass(AliPID::kKaon);
549   if (k==4) return AliPID::ParticleMass(AliPID::kProton);
550   AliWarning("Undefined mass !");
551   return AliPID::ParticleMass(AliPID::kPion);
552 }
553
554 //______________________________________________________________________________
555 Double_t AliESDtrack::E() const
556 {
557   // Returns the energy of the particle given its assumed mass.
558   // Assumes the pion mass if the particle can't be identified properly.
559   
560   Double_t m = M();
561   Double_t p = P();
562   return TMath::Sqrt(p*p + m*m);
563 }
564
565 //______________________________________________________________________________
566 Double_t AliESDtrack::Y() const
567 {
568   // Returns the rapidity of a particle given its assumed mass.
569   // Assumes the pion mass if the particle can't be identified properly.
570   
571   Double_t e = E();
572   Double_t pz = Pz();
573   if (e != TMath::Abs(pz)) { // energy was not equal to pz
574     return 0.5*TMath::Log((e+pz)/(e-pz));
575   } else { // energy was equal to pz
576     return -999.;
577   }
578 }
579
580 //_______________________________________________________________________
581 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
582   //
583   // This function updates track's running parameters 
584   //
585   Int_t *index=0;
586   Bool_t rc=kTRUE;
587
588   SetStatus(flags);
589   fLabel=t->GetLabel();
590
591   if (t->IsStartedTimeIntegral()) {
592     SetStatus(kTIME);
593     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
594     SetIntegratedLength(t->GetIntegratedLength());
595   }
596
597   Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
598   
599   switch (flags) {
600     
601   case kITSin: case kITSout: case kITSrefit:
602     fITSncls=t->GetNumberOfClusters();
603     index=fFriendTrack->GetITSindices(); 
604     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
605         index[i]=t->GetClusterIndex(i);
606         if (i<fITSncls) {
607            Int_t l=(index[i] & 0xf0000000) >> 28;
608            SETBIT(fITSClusterMap,l);                 
609         }
610     }
611     fITSchi2=t->GetChi2();
612     fITSsignal=t->GetPIDsignal();
613     fITSLabel = t->GetLabel();
614     break;
615     
616   case kTPCin: case kTPCrefit:
617     fTPCLabel = t->GetLabel();
618     if (flags==kTPCin)  fTPCInner=new AliExternalTrackParam(*t);
619     if (!fIp) fIp=new AliExternalTrackParam(*t);
620     else 
621       fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
622   case kTPCout:
623     index=fFriendTrack->GetTPCindices(); 
624     if (flags & kTPCout){
625       if (!fOp) fOp=new AliExternalTrackParam(*t);
626       else 
627         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
628     }
629     fTPCncls=t->GetNumberOfClusters();    
630     fTPCchi2=t->GetChi2();
631     
632      {//prevrow must be declared in separate namespace, otherwise compiler cries:
633       //"jump to case label crosses initialization of `Int_t prevrow'"
634        Int_t prevrow = -1;
635        //       for (Int_t i=0;i<fTPCncls;i++) 
636        for (Int_t i=0;i<AliESDfriendTrack::kMaxTPCcluster;i++) 
637         {
638           index[i]=t->GetClusterIndex(i);
639           Int_t idx = index[i];
640
641           if (idx<0) continue; 
642
643           // Piotr's Cluster Map for HBT  
644           // ### please change accordingly if cluster array is changing 
645           // to "New TPC Tracking" style (with gaps in array) 
646           Int_t sect = (idx&0xff000000)>>24;
647           Int_t row = (idx&0x00ff0000)>>16;
648           if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
649
650           fTPCClusterMap.SetBitNumber(row,kTRUE);
651
652           //Fill the gap between previous row and this row with 0 bits
653           //In case  ###  pleas change it as well - just set bit 0 in case there 
654           //is no associated clusters for current "i"
655           if (prevrow < 0) 
656            {
657              prevrow = row;//if previous bit was not assigned yet == this is the first one
658            }
659           else
660            { //we don't know the order (inner to outer or reverse)
661              //just to be save in case it is going to change
662              Int_t n = 0, m = 0;
663              if (prevrow < row)
664               {
665                 n = prevrow;
666                 m = row;
667               }
668              else
669               {
670                 n = row;
671                 m = prevrow;
672               }
673
674              for (Int_t j = n+1; j < m; j++)
675               {
676                 fTPCClusterMap.SetBitNumber(j,kFALSE);
677               }
678              prevrow = row; 
679            }
680           // End Of Piotr's Cluster Map for HBT
681         }
682      }
683     fTPCsignal=t->GetPIDsignal();
684     break;
685
686   case kTRDout: case kTRDin: case kTRDrefit:
687     index=fFriendTrack->GetTRDindices();
688     fTRDLabel = t->GetLabel(); 
689     fTRDncls=t->GetNumberOfClusters();
690     fTRDchi2=t->GetChi2();
691     for (Int_t i=0;i<fTRDncls;i++) index[i]=t->GetClusterIndex(i);
692     fTRDsignal=t->GetPIDsignal();
693     break;
694   case kTRDbackup:
695     if (!fOp) fOp=new AliExternalTrackParam(*t);
696     else 
697       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
698     fTRDncls0 = t->GetNumberOfClusters(); 
699     break;
700   case kTOFin: 
701     break;
702   case kTOFout: 
703     break;
704   case kTRDStop:
705     break;
706   default: 
707     AliError("Wrong flag !");
708     return kFALSE;
709   }
710
711   return rc;
712 }
713
714 //_______________________________________________________________________
715 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
716   //---------------------------------------------------------------------
717   // This function returns external representation of the track parameters
718   //---------------------------------------------------------------------
719   x=GetX();
720   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
721 }
722
723 //_______________________________________________________________________
724 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
725   //---------------------------------------------------------------------
726   // This function returns external representation of the cov. matrix
727   //---------------------------------------------------------------------
728   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
729 }
730
731 //_______________________________________________________________________
732 Bool_t AliESDtrack::GetConstrainedExternalParameters
733                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
734   //---------------------------------------------------------------------
735   // This function returns the constrained external track parameters
736   //---------------------------------------------------------------------
737   if (!fCp) return kFALSE;
738   alpha=fCp->GetAlpha();
739   x=fCp->GetX();
740   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
741   return kTRUE;
742 }
743
744 //_______________________________________________________________________
745 Bool_t 
746 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
747   //---------------------------------------------------------------------
748   // This function returns the constrained external cov. matrix
749   //---------------------------------------------------------------------
750   if (!fCp) return kFALSE;
751   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
752   return kTRUE;
753 }
754
755 Bool_t
756 AliESDtrack::GetInnerExternalParameters
757                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
758   //---------------------------------------------------------------------
759   // This function returns external representation of the track parameters 
760   // at the inner layer of TPC
761   //---------------------------------------------------------------------
762   if (!fIp) return kFALSE;
763   alpha=fIp->GetAlpha();
764   x=fIp->GetX();
765   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
766   return kTRUE;
767 }
768
769 Bool_t 
770 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
771  //---------------------------------------------------------------------
772  // This function returns external representation of the cov. matrix 
773  // at the inner layer of TPC
774  //---------------------------------------------------------------------
775   if (!fIp) return kFALSE;
776   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
777   return kTRUE;
778 }
779
780 Bool_t 
781 AliESDtrack::GetOuterExternalParameters
782                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
783   //---------------------------------------------------------------------
784   // This function returns external representation of the track parameters 
785   // at the inner layer of TRD
786   //---------------------------------------------------------------------
787   if (!fOp) return kFALSE;
788   alpha=fOp->GetAlpha();
789   x=fOp->GetX();
790   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
791   return kTRUE;
792 }
793
794 Bool_t 
795 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
796  //---------------------------------------------------------------------
797  // This function returns external representation of the cov. matrix 
798  // at the inner layer of TRD
799  //---------------------------------------------------------------------
800   if (!fOp) return kFALSE;
801   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
802   return kTRUE;
803 }
804
805 Int_t AliESDtrack::GetNcls(Int_t idet) const
806 {
807   // Get number of clusters by subdetector index
808   //
809   Int_t ncls = 0;
810   switch(idet){
811   case 0:
812     ncls = fITSncls;
813     break;
814   case 1:
815     ncls = fTPCncls;
816     break;
817   case 2:
818     ncls = fTRDncls;
819     break;
820   case 3:
821     if (fTOFindex != 0)
822       ncls = 1;
823     break;
824   default:
825     break;
826   }
827   return ncls;
828 }
829
830 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
831 {
832   // Get cluster index array by subdetector index
833   //
834   Int_t ncls = 0;
835   switch(idet){
836   case 0:
837     ncls = GetITSclusters(idx);
838     break;
839   case 1:
840     ncls = GetTPCclusters(idx);
841     break;
842   case 2:
843     ncls = GetTRDclusters(idx);
844     break;
845   case 3:
846     if (fTOFindex != 0) {
847       idx[0] = GetTOFcluster();
848       ncls = 1;
849     }
850     break;
851   default:
852     break;
853   }
854   return ncls;
855 }
856
857 //_______________________________________________________________________
858 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
859   // Returns the array with integrated times for each particle hypothesis
860   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
861 }
862
863 //_______________________________________________________________________
864 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
865   // Sets the array with integrated times for each particle hypotesis
866   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
867 }
868
869 //_______________________________________________________________________
870 void AliESDtrack::SetITSpid(const Double_t *p) {
871   // Sets values for the probability of each particle type (in ITS)
872   SetPIDValues(fITSr,p,AliPID::kSPECIES);
873   SetStatus(AliESDtrack::kITSpid);
874 }
875
876 //_______________________________________________________________________
877 void AliESDtrack::GetITSpid(Double_t *p) const {
878   // Gets the probability of each particle type (in ITS)
879   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
880 }
881
882 //_______________________________________________________________________
883 Char_t AliESDtrack::GetITSclusters(Int_t *idx) const {
884   //---------------------------------------------------------------------
885   // This function returns indices of the assgined ITS clusters 
886   //---------------------------------------------------------------------
887   if (idx!=0) {
888      Int_t *index=fFriendTrack->GetITSindices();
889      for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) idx[i]=index[i];
890   }
891   return fITSncls;
892 }
893
894 //_______________________________________________________________________
895 UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
896   //---------------------------------------------------------------------
897   // This function returns indices of the assgined ITS clusters 
898   //---------------------------------------------------------------------
899   if (idx!=0) {
900     Int_t *index=fFriendTrack->GetTPCindices();
901     for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
902   }
903   return fTPCncls;
904 }
905
906 Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
907   //
908   // GetDensity of the clusters on given region between row0 and row1
909   // Dead zone effect takin into acoount
910   //
911   Int_t good  = 0;
912   Int_t found = 0;
913   //  
914   Int_t *index=fFriendTrack->GetTPCindices();
915   for (Int_t i=row0;i<=row1;i++){     
916     Int_t idx = index[i];
917     if (idx!=-1)  good++;             // track outside of dead zone
918     if (idx>0)    found++;
919   }
920   Float_t density=0.5;
921   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
922   return density;
923 }
924
925 //_______________________________________________________________________
926 void AliESDtrack::SetTPCpid(const Double_t *p) {  
927   // Sets values for the probability of each particle type (in TPC)
928   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
929   SetStatus(AliESDtrack::kTPCpid);
930 }
931
932 //_______________________________________________________________________
933 void AliESDtrack::GetTPCpid(Double_t *p) const {
934   // Gets the probability of each particle type (in TPC)
935   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
936 }
937
938 //_______________________________________________________________________
939 UChar_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
940   //---------------------------------------------------------------------
941   // This function returns indices of the assgined TRD clusters 
942   //---------------------------------------------------------------------
943   if (idx!=0) {
944      Int_t *index=fFriendTrack->GetTRDindices();
945      for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
946   }
947   return fTRDncls;
948 }
949
950 //_______________________________________________________________________
951 void AliESDtrack::SetTRDpid(const Double_t *p) {  
952   // Sets values for the probability of each particle type (in TRD)
953   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
954   SetStatus(AliESDtrack::kTRDpid);
955 }
956
957 //_______________________________________________________________________
958 void AliESDtrack::GetTRDpid(Double_t *p) const {
959   // Gets the probability of each particle type (in TRD)
960   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
961 }
962
963 //_______________________________________________________________________
964 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
965 {
966   // Sets the probability of particle type iSpecies to p (in TRD)
967   fTRDr[iSpecies] = p;
968 }
969
970 Double_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
971 {
972   // Returns the probability of particle type iSpecies (in TRD)
973   return fTRDr[iSpecies];
974 }
975
976 //_______________________________________________________________________
977 void AliESDtrack::SetTOFpid(const Double_t *p) {  
978   // Sets the probability of each particle type (in TOF)
979   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
980   SetStatus(AliESDtrack::kTOFpid);
981 }
982
983 //_______________________________________________________________________
984 void AliESDtrack::SetTOFLabel(const Int_t *p) {  
985   // Sets  (in TOF)
986   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
987 }
988
989 //_______________________________________________________________________
990 void AliESDtrack::GetTOFpid(Double_t *p) const {
991   // Gets probabilities of each particle type (in TOF)
992   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
993 }
994
995 //_______________________________________________________________________
996 void AliESDtrack::GetTOFLabel(Int_t *p) const {
997   // Gets (in TOF)
998   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
999 }
1000
1001 //_______________________________________________________________________
1002 void AliESDtrack::GetTOFInfo(Float_t *info) const {
1003   // Gets (in TOF)
1004   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
1005 }
1006
1007 //_______________________________________________________________________
1008 void AliESDtrack::SetTOFInfo(Float_t*info) {
1009   // Gets (in TOF)
1010   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
1011 }
1012
1013
1014
1015 //_______________________________________________________________________
1016 void AliESDtrack::SetHMPIDpid(const Double_t *p) {  
1017   // Sets the probability of each particle type (in HMPID)
1018   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
1019   SetStatus(AliESDtrack::kHMPIDpid);
1020 }
1021
1022 //_______________________________________________________________________
1023 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
1024   // Gets probabilities of each particle type (in HMPID)
1025   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
1026 }
1027
1028
1029
1030 //_______________________________________________________________________
1031 void AliESDtrack::SetESDpid(const Double_t *p) {  
1032   // Sets the probability of each particle type for the ESD track
1033   SetPIDValues(fR,p,AliPID::kSPECIES);
1034   SetStatus(AliESDtrack::kESDpid);
1035 }
1036
1037 //_______________________________________________________________________
1038 void AliESDtrack::GetESDpid(Double_t *p) const {
1039   // Gets probability of each particle type for the ESD track
1040   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
1041 }
1042
1043 //_______________________________________________________________________
1044 Bool_t AliESDtrack::RelateToVertex
1045 (const AliESDVertex *vtx, Double_t b, Double_t maxd) {
1046   //
1047   // Try to relate this track to the vertex "vtx", 
1048   // if the (rough) transverse impact parameter is not bigger then "maxd". 
1049   //            Magnetic field is "b" (kG).
1050   //
1051   // a) The track gets extapolated to the DCA to the vertex.
1052   // b) The impact parameters and their covariance matrix are calculated.
1053   // c) An attempt to constrain this track to the vertex is done.
1054   //
1055   //    In the case of success, the returned value is kTRUE
1056   //    (otherwise, it's kFALSE)
1057   //  
1058
1059   if (!vtx) return kFALSE;
1060
1061   Double_t alpha=GetAlpha();
1062   Double_t sn=TMath::Sin(alpha), cs=TMath::Cos(alpha);
1063   Double_t x=GetX(), y=GetParameter()[0], snp=GetParameter()[2];
1064   Double_t xv= vtx->GetXv()*cs + vtx->GetYv()*sn;
1065   Double_t yv=-vtx->GetXv()*sn + vtx->GetYv()*cs, zv=vtx->GetZv();
1066   x-=xv; y-=yv;
1067
1068   //Estimate the impact parameter neglecting the track curvature
1069   Double_t d=TMath::Abs(x*snp - y*TMath::Sqrt(1.- snp*snp));
1070   if (d > maxd) return kFALSE; 
1071
1072   //Propagate to the DCA
1073   Double_t crv=kB2C*b*GetParameter()[4];
1074   if (TMath::Abs(b) < kAlmost0Field) crv=0.;
1075
1076   Double_t tgfv=-(crv*x - snp)/(crv*y + TMath::Sqrt(1.-snp*snp));
1077   sn=tgfv/TMath::Sqrt(1.+ tgfv*tgfv);
1078   if (TMath::Abs(tgfv)>0.) cs = sn/tgfv;
1079   else cs=1.;
1080
1081   x = xv*cs + yv*sn;
1082   yv=-xv*sn + yv*cs; xv=x;
1083
1084   if (!Propagate(alpha+TMath::ASin(sn),xv,b)) return kFALSE;
1085
1086   fD = GetParameter()[0] - yv;
1087   fZ = GetParameter()[1] - zv;
1088   
1089   Double_t cov[6]; vtx->GetCovMatrix(cov);
1090
1091   //***** Improvements by A.Dainese
1092   alpha=GetAlpha(); sn=TMath::Sin(alpha); cs=TMath::Cos(alpha);
1093   Double_t s2ylocvtx = cov[0]*sn*sn + cov[2]*cs*cs - 2.*cov[1]*cs*sn;
1094   fCdd = GetCovariance()[0] + s2ylocvtx;   // neglecting correlations
1095   fCdz = GetCovariance()[1];               // between (x,y) and z
1096   fCzz = GetCovariance()[2] + cov[5];      // in vertex's covariance matrix
1097   //*****
1098
1099   {//Try to constrain 
1100     Double_t p[2]={yv,zv}, c[3]={cov[2],0.,cov[5]};
1101     Double_t chi2=GetPredictedChi2(p,c);
1102
1103     if (chi2>77.) return kFALSE;
1104
1105     AliExternalTrackParam tmp(*this);
1106     if (!tmp.Update(p,c)) return kFALSE;
1107
1108     fCchi2=chi2;
1109     if (!fCp) fCp=new AliExternalTrackParam();
1110     new (fCp) AliExternalTrackParam(tmp);
1111   }
1112
1113   return kTRUE;
1114 }
1115
1116 //_______________________________________________________________________
1117 void AliESDtrack::Print(Option_t *) const {
1118   // Prints info on the track
1119   
1120   printf("ESD track info\n") ; 
1121   Double_t p[AliPID::kSPECIESN] ; 
1122   Int_t index = 0 ; 
1123   if( IsOn(kITSpid) ){
1124     printf("From ITS: ") ; 
1125     GetITSpid(p) ; 
1126     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1127       printf("%f, ", p[index]) ;
1128     printf("\n           signal = %f\n", GetITSsignal()) ;
1129   } 
1130   if( IsOn(kTPCpid) ){
1131     printf("From TPC: ") ; 
1132     GetTPCpid(p) ; 
1133     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1134       printf("%f, ", p[index]) ;
1135     printf("\n           signal = %f\n", GetTPCsignal()) ;
1136   }
1137   if( IsOn(kTRDpid) ){
1138     printf("From TRD: ") ; 
1139     GetTRDpid(p) ; 
1140     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1141       printf("%f, ", p[index]) ;
1142     printf("\n           signal = %f\n", GetTRDsignal()) ;
1143   }
1144   if( IsOn(kTOFpid) ){
1145     printf("From TOF: ") ; 
1146     GetTOFpid(p) ; 
1147     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1148       printf("%f, ", p[index]) ;
1149     printf("\n           signal = %f\n", GetTOFsignal()) ;
1150   }
1151   if( IsOn(kHMPIDpid) ){
1152     printf("From HMPID: ") ; 
1153     GetHMPIDpid(p) ; 
1154     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1155       printf("%f, ", p[index]) ;
1156     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
1157   }
1158
1159