]> git.uio.no Git - u/mrichter/AliRoot.git/blob - STEER/ESD/AliESDtrack.cxx
7b80878247d859d8ffea6cddb7043be73a05f476
[u/mrichter/AliRoot.git] / STEER / ESD / 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 //
24 //  What do you need to know before starting analysis
25 //  (by Marian Ivanov: marian.ivanov@cern.ch)
26 //
27 //
28 //   AliESDtrack:
29 //   1.  What is the AliESDtrack
30 //   2.  What informations do we store
31 //   3.  How to use the information for analysis
32 //   
33 //
34 //   1.AliESDtrack is the container of the information about the track/particle
35 //     reconstructed during Barrel Tracking.
36 //     The track information is propagated from one tracking detector to 
37 //     other using the functionality of AliESDtrack - Current parameters.  
38 //
39 //     No global fit model is used.
40 //     Barrel tracking use Kalman filtering technique, it gives optimal local 
41 //     track parameters at given point under certian assumptions.
42 //             
43 //     Kalman filter take into account additional effect which are 
44 //     difficult to handle using global fit.
45 //     Effects:
46 //        a.) Multiple scattering
47 //        b.) Energy loss
48 //        c.) Non homogenous magnetic field
49 //
50 //     In general case, following barrel detectors are contributing to 
51 //     the Kalman track information:
52 //         a. TPC
53 //         b. ITS
54 //         c. TRD
55 //
56 //      In general 3 reconstruction itteration are performed:
57 //         1. Find tracks   - sequence TPC->ITS
58 //         2. PropagateBack - sequence ITS->TPC->TRD -> Outer PID detectors
59 //         3. Refit invward - sequence TRD->TPC->ITS
60 //      The current tracks are updated after each detector (see bellow).
61 //      In specical cases a track  sanpshots are stored.   
62 // 
63 //
64 //      For some type of analysis (+visualization) track local parameters at 
65 //      different position are neccesary. A snapshots during the track 
66 //      propagation are created.
67 //      (See AliExternalTrackParam class for desctiption of variables and 
68 //      functionality)
69 //      Snapshots:
70 //      a. Current parameters - class itself (AliExternalTrackParam)
71 //         Contributors: general case TRD->TPC->ITS
72 //         Preferable usage:  Decission  - primary or secondary track
73 //         NOTICE - By default the track parameters are stored at the DCA point
74 //                  to the primary vertex. optimal for primary tracks, 
75 //                  far from optimal for secondary tracks.
76 //      b. Constrained parameters - Kalman information updated with 
77 //         the Primary vertex information 
78 //         Contributors: general case TRD->TPC->ITS
79 //         Preferable usage: Use only for tracks selected as primary
80 //         NOTICE - not real constrain - taken as additional measurement 
81 //         with corresponding error
82 //         Function:  
83 //       const AliExternalTrackParam *GetConstrainedParam() const {return fCp;}
84 //      c. Inner parameters -  Track parameters at inner wall of the TPC 
85 //         Contributors: general case TRD->TPC
86 //         function:
87 //           const AliExternalTrackParam *GetInnerParam() const { return fIp;}
88 //
89 //      d. TPCinnerparam  - contributors - TPC only
90 //         Contributors:  TPC
91 //         Preferable usage: Requested for HBT study 
92 //                         (smaller correlations as using also ITS information)
93 //         NOTICE - the track parameters are propagated to the DCA to  
94 //         to primary vertex
95 //         Optimal for primary, far from optimal for secondary tracks
96 //         Function:
97 //    const AliExternalTrackParam *GetTPCInnerParam() const {return fTPCInner;}
98 //     
99 //      e. Outer parameters - 
100 //           Contributors-  general case - ITS-> TPC -> TRD
101 //           The last point - Outer parameters radius is determined
102 //           e.a) Local inclination angle bigger than threshold - 
103 //                Low momenta tracks 
104 //           e.a) Catastrofic energy losss in material
105 //           e.b) Not further improvement (no space points)
106 //           Usage:             
107 //              a.) Tracking: Starting parameter for Refit inward 
108 //              b.) Visualization
109 //              c.) QA
110 //         NOTICE: Should be not used for the physic analysis
111 //         Function:
112 //            const AliExternalTrackParam *GetOuterParam() const { return fOp;}
113 //
114 //-----------------------------------------------------------------
115
116 #include <TMath.h>
117 #include <TParticle.h>
118 #include <TDatabasePDG.h>
119 #include <TMatrixD.h>
120
121 #include "AliESDVertex.h"
122 #include "AliESDtrack.h"
123 #include "AliESDEvent.h"
124 #include "AliKalmanTrack.h"
125 #include "AliVTrack.h"
126 #include "AliLog.h"
127 #include "AliTrackPointArray.h"
128 #include "TPolyMarker3D.h"
129 #include "AliTrackerBase.h"
130 #include "AliTPCdEdxInfo.h"
131
132 ClassImp(AliESDtrack)
133
134 void SetPIDValues(Double_t * dest, const Double_t * src, Int_t n) {
135   // This function copies "n" PID weights from "scr" to "dest"
136   // and normalizes their sum to 1 thus producing conditional probabilities.
137   // The negative weights are set to 0.
138   // In case all the weights are non-positive they are replaced by
139   // uniform probabilities
140
141   if (n<=0) return;
142
143   Float_t uniform = 1./(Float_t)n;
144
145   Float_t sum = 0;
146   for (Int_t i=0; i<n; i++) 
147     if (src[i]>=0) {
148       sum+=src[i];
149       dest[i] = src[i];
150     }
151     else {
152       dest[i] = 0;
153     }
154
155   if(sum>0)
156     for (Int_t i=0; i<n; i++) dest[i] /= sum;
157   else
158     for (Int_t i=0; i<n; i++) dest[i] = uniform;
159 }
160
161 //_______________________________________________________________________
162 AliESDtrack::AliESDtrack() : 
163   AliExternalTrackParam(),
164   fCp(0),
165   fIp(0),
166   fTPCInner(0),
167   fOp(0),
168   fHMPIDp(0),  
169   fFriendTrack(NULL),
170   fTPCFitMap(159),//number of padrows
171   fTPCClusterMap(159),//number of padrows
172   fTPCSharedMap(159),//number of padrows
173   fFlags(0),
174   fID(0),
175   fLabel(0),
176   fITSLabel(0),
177   fTPCLabel(0),
178   fTRDLabel(0),
179   fTOFCalChannel(-1),
180   fTOFindex(-1),
181   fHMPIDqn(0),
182   fHMPIDcluIdx(-1),
183   fCaloIndex(kEMCALNoMatch),
184   fHMPIDtrkTheta(0),
185   fHMPIDtrkPhi(0),
186   fHMPIDsignal(0),
187   fTrackLength(0),
188   fdTPC(0),fzTPC(0),
189   fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
190   fCchi2TPC(0),
191   fD(0),fZ(0),
192   fCdd(0),fCdz(0),fCzz(0),
193   fCchi2(0),
194   fITSchi2(0),
195   fTPCchi2(0),
196   fTPCchi2Iter1(0),
197   fTRDchi2(0),
198   fTOFchi2(0),
199   fHMPIDchi2(0),
200   fGlobalChi2(0),
201   fITSsignal(0),
202   fTPCsignal(0),
203   fTPCsignalS(0),
204   fTPCdEdxInfo(0),
205   fTRDsignal(0),
206   fTRDQuality(0),
207   fTRDBudget(0),
208   fTOFsignal(99999),
209   fTOFsignalToT(99999),
210   fTOFsignalRaw(99999),
211   fTOFsignalDz(999),
212   fTOFsignalDx(999),
213   fTOFdeltaBC(999),
214   fTOFl0l1(999),
215   fCaloDx(0),
216   fCaloDz(0),
217   fHMPIDtrkX(0),
218   fHMPIDtrkY(0),
219   fHMPIDmipX(0),
220   fHMPIDmipY(0),
221   fTPCncls(0),
222   fTPCnclsF(0),
223   fTPCsignalN(0),
224   fTPCnclsIter1(0),
225   fTPCnclsFIter1(0),
226   fITSncls(0),
227   fITSClusterMap(0),
228   fITSSharedMap(0),
229   fTRDncls(0),
230   fTRDncls0(0),
231   fTRDntracklets(0),
232   fTRDnSlices(0),
233   fTRDslices(0x0),
234   fVertexID(-2),// -2 means an orphan track 
235   fESDEvent(0),
236   fCacheNCrossedRows(-10),
237   fCacheChi2TPCConstrainedVsGlobal(-10),
238   fCacheChi2TPCConstrainedVsGlobalVertex(0)
239 {
240   //
241   // The default ESD constructor 
242   //
243   if (!OnlineMode()) fFriendTrack=new AliESDfriendTrack();
244
245   Int_t i;
246   for (i=kNITSchi2Std;i--;) fITSchi2Std[i] = 0;
247   for (i=0; i<AliPID::kSPECIES; i++) {
248     fTrackTime[i]=0.;
249     fR[i]=0.;
250     fITSr[i]=0.;
251     fTPCr[i]=0.;
252     fTRDr[i]=0.;
253     fTOFr[i]=0.;
254     fHMPIDr[i]=0.;
255   }
256   
257   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
258   for (i=0; i<3; i++)   { fV0Indexes[i]=0;}
259   for (i=0;i<kTRDnPlanes;i++) {
260     fTRDTimBin[i]=0;
261   }
262   for (i=0;i<4;i++) {fITSdEdxSamples[i]=0.;}
263   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
264   for (i=0;i<3;i++) {fTOFLabel[i]=-1;}
265   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
266   for (i=0;i<12;i++) {fITSModule[i]=-1;}
267 }
268
269 bool AliESDtrack::fgkOnlineMode=false;
270
271 //_______________________________________________________________________
272 AliESDtrack::AliESDtrack(const AliESDtrack& track):
273   AliExternalTrackParam(track),
274   fCp(0),
275   fIp(0),
276   fTPCInner(0),
277   fOp(0),
278   fHMPIDp(0),  
279   fFriendTrack(0),
280   fTPCFitMap(track.fTPCFitMap),
281   fTPCClusterMap(track.fTPCClusterMap),
282   fTPCSharedMap(track.fTPCSharedMap),
283   fFlags(track.fFlags),
284   fID(track.fID),
285   fLabel(track.fLabel),
286   fITSLabel(track.fITSLabel),
287   fTPCLabel(track.fTPCLabel),
288   fTRDLabel(track.fTRDLabel),
289   fTOFCalChannel(track.fTOFCalChannel),
290   fTOFindex(track.fTOFindex),
291   fHMPIDqn(track.fHMPIDqn),
292   fHMPIDcluIdx(track.fHMPIDcluIdx),
293   fCaloIndex(track.fCaloIndex),
294   fHMPIDtrkTheta(track.fHMPIDtrkTheta),
295   fHMPIDtrkPhi(track.fHMPIDtrkPhi),
296   fHMPIDsignal(track.fHMPIDsignal),
297   fTrackLength(track.fTrackLength),
298   fdTPC(track.fdTPC),fzTPC(track.fzTPC),
299   fCddTPC(track.fCddTPC),fCdzTPC(track.fCdzTPC),fCzzTPC(track.fCzzTPC),
300   fCchi2TPC(track.fCchi2TPC),
301   fD(track.fD),fZ(track.fZ),
302   fCdd(track.fCdd),fCdz(track.fCdz),fCzz(track.fCzz),
303   fCchi2(track.fCchi2),
304   fITSchi2(track.fITSchi2),
305   fTPCchi2(track.fTPCchi2),
306   fTPCchi2Iter1(track.fTPCchi2Iter1),
307   fTRDchi2(track.fTRDchi2),
308   fTOFchi2(track.fTOFchi2),
309   fHMPIDchi2(track.fHMPIDchi2),
310   fGlobalChi2(track.fGlobalChi2),
311   fITSsignal(track.fITSsignal),
312   fTPCsignal(track.fTPCsignal),
313   fTPCsignalS(track.fTPCsignalS),
314   fTPCdEdxInfo(0),
315   fTRDsignal(track.fTRDsignal),
316   fTRDQuality(track.fTRDQuality),
317   fTRDBudget(track.fTRDBudget),
318   fTOFsignal(track.fTOFsignal),
319   fTOFsignalToT(track.fTOFsignalToT),
320   fTOFsignalRaw(track.fTOFsignalRaw),
321   fTOFsignalDz(track.fTOFsignalDz),
322   fTOFsignalDx(track.fTOFsignalDx),
323   fTOFdeltaBC(track.fTOFdeltaBC),
324   fTOFl0l1(track.fTOFl0l1),
325   fCaloDx(track.fCaloDx),
326   fCaloDz(track.fCaloDz),
327   fHMPIDtrkX(track.fHMPIDtrkX),
328   fHMPIDtrkY(track.fHMPIDtrkY),
329   fHMPIDmipX(track.fHMPIDmipX),
330   fHMPIDmipY(track.fHMPIDmipY),
331   fTPCncls(track.fTPCncls),
332   fTPCnclsF(track.fTPCnclsF),
333   fTPCsignalN(track.fTPCsignalN),
334   fTPCnclsIter1(track.fTPCnclsIter1),
335   fTPCnclsFIter1(track.fTPCnclsIter1),
336   fITSncls(track.fITSncls),
337   fITSClusterMap(track.fITSClusterMap),
338   fITSSharedMap(track.fITSSharedMap),
339   fTRDncls(track.fTRDncls),
340   fTRDncls0(track.fTRDncls0),
341   fTRDntracklets(track.fTRDntracklets),
342   fTRDnSlices(track.fTRDnSlices),
343   fTRDslices(0x0),
344   fVertexID(track.fVertexID),
345   fESDEvent(track.fESDEvent),
346   fCacheNCrossedRows(track.fCacheNCrossedRows),
347   fCacheChi2TPCConstrainedVsGlobal(track.fCacheChi2TPCConstrainedVsGlobal),
348   fCacheChi2TPCConstrainedVsGlobalVertex(track.fCacheChi2TPCConstrainedVsGlobalVertex)
349 {
350   //
351   //copy constructor
352   //
353   for (Int_t i=kNITSchi2Std;i--;) fITSchi2Std[i] = track.fTrackTime[i];
354   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i]=track.fTrackTime[i];
355   for (Int_t i=0;i<AliPID::kSPECIES;i++)  fR[i]=track.fR[i];
356   //
357   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=track.fITSr[i]; 
358   //
359   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=track.fTPCr[i]; 
360   for (Int_t i=0;i<4;i++) {fITSdEdxSamples[i]=track.fITSdEdxSamples[i];}
361   for (Int_t i=0;i<4;i++) {fTPCPoints[i]=track.fTPCPoints[i];}
362   for (Int_t i=0; i<3;i++)   { fKinkIndexes[i]=track.fKinkIndexes[i];}
363   for (Int_t i=0; i<3;i++)   { fV0Indexes[i]=track.fV0Indexes[i];}
364   //
365   for (Int_t i=0;i<kTRDnPlanes;i++) {
366     fTRDTimBin[i]=track.fTRDTimBin[i];
367   }
368
369   if (fTRDnSlices) {
370     fTRDslices=new Double32_t[fTRDnSlices];
371     for (Int_t i=0; i<fTRDnSlices; i++) fTRDslices[i]=track.fTRDslices[i];
372   }
373
374
375   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i]=track.fTRDr[i]; 
376   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i]=track.fTOFr[i];
377   for (Int_t i=0;i<3;i++) fTOFLabel[i]=track.fTOFLabel[i];
378   for (Int_t i=0;i<10;i++) fTOFInfo[i]=track.fTOFInfo[i];
379   for (Int_t i=0;i<12;i++) fITSModule[i]=track.fITSModule[i];
380   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i]=track.fHMPIDr[i];
381
382   if (track.fCp) fCp=new AliExternalTrackParam(*track.fCp);
383   if (track.fIp) fIp=new AliExternalTrackParam(*track.fIp);
384   if (track.fTPCInner) fTPCInner=new AliExternalTrackParam(*track.fTPCInner);
385   if (track.fOp) fOp=new AliExternalTrackParam(*track.fOp);
386   if (track.fHMPIDp) fHMPIDp=new AliExternalTrackParam(*track.fHMPIDp);
387   if (track.fTPCdEdxInfo) fTPCdEdxInfo = new AliTPCdEdxInfo(*track.fTPCdEdxInfo);
388
389   
390   if (track.fFriendTrack) fFriendTrack=new AliESDfriendTrack(*(track.fFriendTrack));
391 }
392
393 //_______________________________________________________________________
394 AliESDtrack::AliESDtrack(const AliVTrack *track) : 
395   AliExternalTrackParam(track),
396   fCp(0),
397   fIp(0),
398   fTPCInner(0),
399   fOp(0),
400   fHMPIDp(0),  
401   fFriendTrack(0),
402   fTPCFitMap(159),//number of padrows
403   fTPCClusterMap(159),//number of padrows
404   fTPCSharedMap(159),//number of padrows
405   fFlags(0),
406   fID(),
407   fLabel(0),
408   fITSLabel(0),
409   fTPCLabel(0),
410   fTRDLabel(0),
411   fTOFCalChannel(-1),
412   fTOFindex(-1),
413   fHMPIDqn(0),
414   fHMPIDcluIdx(-1),
415   fCaloIndex(kEMCALNoMatch),
416   fHMPIDtrkTheta(0),
417   fHMPIDtrkPhi(0),
418   fHMPIDsignal(0),
419   fTrackLength(0),
420   fdTPC(0),fzTPC(0),
421   fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
422   fCchi2TPC(0),
423   fD(0),fZ(0),
424   fCdd(0),fCdz(0),fCzz(0),
425   fCchi2(0),
426   fITSchi2(0),
427   fTPCchi2(0),
428   fTPCchi2Iter1(0),
429   fTRDchi2(0),
430   fTOFchi2(0),
431   fHMPIDchi2(0),
432   fGlobalChi2(0),
433   fITSsignal(0),
434   fTPCsignal(0),
435   fTPCsignalS(0),
436   fTPCdEdxInfo(0),
437   fTRDsignal(0),
438   fTRDQuality(0),
439   fTRDBudget(0),
440   fTOFsignal(99999),
441   fTOFsignalToT(99999),
442   fTOFsignalRaw(99999),
443   fTOFsignalDz(999),
444   fTOFsignalDx(999),
445   fTOFdeltaBC(999),
446   fTOFl0l1(999),
447   fCaloDx(0),
448   fCaloDz(0),
449   fHMPIDtrkX(0),
450   fHMPIDtrkY(0),
451   fHMPIDmipX(0),
452   fHMPIDmipY(0),
453   fTPCncls(0),
454   fTPCnclsF(0),
455   fTPCsignalN(0),
456   fTPCnclsIter1(0),
457   fTPCnclsFIter1(0),
458   fITSncls(0),
459   fITSClusterMap(0),
460   fITSSharedMap(0),
461   fTRDncls(0),
462   fTRDncls0(0),
463   fTRDntracklets(0),
464   fTRDnSlices(0),
465   fTRDslices(0x0),
466   fVertexID(-2),  // -2 means an orphan track
467   fESDEvent(0),
468   fCacheNCrossedRows(-10),
469   fCacheChi2TPCConstrainedVsGlobal(-10),
470   fCacheChi2TPCConstrainedVsGlobalVertex(0)
471 {
472   //
473   // ESD track from AliVTrack.
474   // This is not a copy constructor !
475   //
476
477   if (track->InheritsFrom("AliExternalTrackParam")) {
478      AliError("This is not a copy constructor. Use AliESDtrack(const AliESDtrack &) !");
479      AliWarning("Calling the default constructor...");
480      AliESDtrack();
481      return;
482   }
483
484   // Reset all the arrays
485   Int_t i;
486   for (i=kNITSchi2Std;i--;) fITSchi2Std[i] = 0;
487   for (i=0; i<AliPID::kSPECIES; i++) {
488     fTrackTime[i]=0.;
489     fR[i]=0.;
490     fITSr[i]=0.;
491     fTPCr[i]=0.;
492     fTRDr[i]=0.;
493     fTOFr[i]=0.;
494     fHMPIDr[i]=0.;
495   }
496   
497   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
498   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
499   for (i=0;i<kTRDnPlanes;i++) {
500     fTRDTimBin[i]=0;
501   }
502   for (i=0;i<4;i++) {fITSdEdxSamples[i]=0.;}
503   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
504   for (i=0;i<3;i++) {fTOFLabel[i]=-1;}
505   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
506   for (i=0;i<12;i++) {fITSModule[i]=-1;}
507
508   // Set the ID
509   SetID(track->GetID());
510
511   // Set ITS cluster map
512   fITSClusterMap=track->GetITSClusterMap();
513   fITSSharedMap=0;
514
515   fITSncls=0;
516   for(i=0; i<6; i++) {
517     if(HasPointOnITSLayer(i)) fITSncls++;
518   }
519
520   // Set TPC ncls 
521   fTPCncls=track->GetTPCNcls();
522
523
524   // Set the combined PID
525   const Double_t *pid = track->PID();
526   if(pid){
527     for (i=0; i<AliPID::kSPECIES; i++) fR[i]=pid[i];
528   }
529   // AliESD track label
530   SetLabel(track->GetLabel());
531   // Set the status
532   SetStatus(track->GetStatus());
533 }
534
535 //_______________________________________________________________________
536 AliESDtrack::AliESDtrack(TParticle * part) : 
537   AliExternalTrackParam(),
538   fCp(0),
539   fIp(0),
540   fTPCInner(0),
541   fOp(0),
542   fHMPIDp(0),  
543   fFriendTrack(0),
544   fTPCFitMap(159),//number of padrows
545   fTPCClusterMap(159),//number of padrows
546   fTPCSharedMap(159),//number of padrows
547   fFlags(0),
548   fID(0),
549   fLabel(0),
550   fITSLabel(0),
551   fTPCLabel(0),
552   fTRDLabel(0),
553   fTOFCalChannel(-1),
554   fTOFindex(-1),
555   fHMPIDqn(0),
556   fHMPIDcluIdx(-1),
557   fCaloIndex(kEMCALNoMatch),
558   fHMPIDtrkTheta(0),
559   fHMPIDtrkPhi(0),
560   fHMPIDsignal(0),
561   fTrackLength(0),
562   fdTPC(0),fzTPC(0),
563   fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
564   fCchi2TPC(0),
565   fD(0),fZ(0),
566   fCdd(0),fCdz(0),fCzz(0),
567   fCchi2(0),
568   fITSchi2(0),
569   fTPCchi2(0),
570   fTPCchi2Iter1(0),  
571   fTRDchi2(0),
572   fTOFchi2(0),
573   fHMPIDchi2(0),
574   fGlobalChi2(0),
575   fITSsignal(0),
576   fTPCsignal(0),
577   fTPCsignalS(0),
578   fTPCdEdxInfo(0),
579   fTRDsignal(0),
580   fTRDQuality(0),
581   fTRDBudget(0),
582   fTOFsignal(99999),
583   fTOFsignalToT(99999),
584   fTOFsignalRaw(99999),
585   fTOFsignalDz(999),
586   fTOFsignalDx(999),
587   fTOFdeltaBC(999),
588   fTOFl0l1(999),
589   fCaloDx(0),
590   fCaloDz(0),
591   fHMPIDtrkX(0),
592   fHMPIDtrkY(0),
593   fHMPIDmipX(0),
594   fHMPIDmipY(0),
595   fTPCncls(0),
596   fTPCnclsF(0),
597   fTPCsignalN(0),
598   fTPCnclsIter1(0),
599   fTPCnclsFIter1(0),
600   fITSncls(0),
601   fITSClusterMap(0),
602   fITSSharedMap(0),
603   fTRDncls(0),
604   fTRDncls0(0),
605   fTRDntracklets(0),
606   fTRDnSlices(0),
607   fTRDslices(0x0),
608   fVertexID(-2),  // -2 means an orphan track
609   fESDEvent(0),
610   fCacheNCrossedRows(-10),
611   fCacheChi2TPCConstrainedVsGlobal(-10),
612   fCacheChi2TPCConstrainedVsGlobalVertex(0)  
613 {
614   //
615   // ESD track from TParticle
616   //
617
618   // Reset all the arrays
619   Int_t i;
620   for (i=kNITSchi2Std;i--;) fITSchi2Std[i] = 0;
621   for (i=0; i<AliPID::kSPECIES; i++) {
622     fTrackTime[i]=0.;
623     fR[i]=0.;
624     fITSr[i]=0.;
625     fTPCr[i]=0.;
626     fTRDr[i]=0.;
627     fTOFr[i]=0.;
628     fHMPIDr[i]=0.;
629   }
630   
631   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
632   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
633   for (i=0;i<kTRDnPlanes;i++) {
634     fTRDTimBin[i]=0;
635   }
636   for (i=0;i<4;i++) {fITSdEdxSamples[i]=0.;}
637   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
638   for (i=0;i<3;i++) {fTOFLabel[i]=-1;}
639   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
640   for (i=0;i<12;i++) {fITSModule[i]=-1;}
641
642   // Calculate the AliExternalTrackParam content
643
644   Double_t xref;
645   Double_t alpha;
646   Double_t param[5];
647   Double_t covar[15];
648
649   // Calculate alpha: the rotation angle of the corresponding local system (TPC sector)
650   alpha = part->Phi()*180./TMath::Pi();
651   if (alpha<0) alpha+= 360.;
652   if (alpha>360) alpha -= 360.;
653
654   Int_t sector = (Int_t)(alpha/20.);
655   alpha = 10. + 20.*sector;
656   alpha /= 180;
657   alpha *= TMath::Pi();
658
659   // Covariance matrix: no errors, the parameters are exact
660   for (i=0; i<15; i++) covar[i]=0.;
661
662   // Get the vertex of origin and the momentum
663   TVector3 ver(part->Vx(),part->Vy(),part->Vz());
664   TVector3 mom(part->Px(),part->Py(),part->Pz());
665
666   // Rotate to the local coordinate system (TPC sector)
667   ver.RotateZ(-alpha);
668   mom.RotateZ(-alpha);
669
670   // X of the referense plane
671   xref = ver.X();
672
673   Int_t pdgCode = part->GetPdgCode();
674
675   Double_t charge = 
676     TDatabasePDG::Instance()->GetParticle(pdgCode)->Charge();
677
678   param[0] = ver.Y();
679   param[1] = ver.Z();
680   param[2] = TMath::Sin(mom.Phi());
681   param[3] = mom.Pz()/mom.Pt();
682   param[4] = TMath::Sign(1/mom.Pt(),charge);
683
684   // Set AliExternalTrackParam
685   Set(xref, alpha, param, covar);
686
687   // Set the PID
688   Int_t indexPID = 99;
689
690   switch (TMath::Abs(pdgCode)) {
691
692   case  11: // electron
693     indexPID = 0;
694     break;
695
696   case 13: // muon
697     indexPID = 1;
698     break;
699
700   case 211: // pion
701     indexPID = 2;
702     break;
703
704   case 321: // kaon
705     indexPID = 3;
706     break;
707
708   case 2212: // proton
709     indexPID = 4;
710     break;
711
712   default:
713     break;
714   }
715
716   // If the particle is not e,mu,pi,K or p the PID probabilities are set to 0
717   if (indexPID < AliPID::kSPECIES) {
718     fR[indexPID]=1.;
719     fITSr[indexPID]=1.;
720     fTPCr[indexPID]=1.;
721     fTRDr[indexPID]=1.;
722     fTOFr[indexPID]=1.;
723     fHMPIDr[indexPID]=1.;
724
725   }
726   // AliESD track label
727   SetLabel(part->GetUniqueID());
728
729 }
730
731 //_______________________________________________________________________
732 AliESDtrack::~AliESDtrack(){ 
733   //
734   // This is destructor according Coding Conventrions 
735   //
736   //printf("Delete track\n");
737   delete fIp; 
738   delete fTPCInner; 
739   delete fOp;
740   delete fHMPIDp;
741   delete fCp; 
742   delete fFriendTrack;
743   delete fTPCdEdxInfo;
744   if(fTRDnSlices)
745     delete[] fTRDslices;
746
747   //Reset cached values - needed for TClonesArray in AliESDInputHandler
748   fCacheNCrossedRows = -10.;
749   fCacheChi2TPCConstrainedVsGlobal = -10.;
750   if(fCacheChi2TPCConstrainedVsGlobalVertex) fCacheChi2TPCConstrainedVsGlobalVertex = 0;
751
752 }
753
754 AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
755   
756
757   if(&source == this) return *this;
758   AliExternalTrackParam::operator=(source);
759
760   
761   if(source.fCp){
762     // we have the trackparam: assign or copy construct
763     if(fCp)*fCp = *source.fCp;
764     else fCp = new AliExternalTrackParam(*source.fCp);
765   }
766   else{
767     // no track param delete the old one
768     delete fCp;
769     fCp = 0;
770   }
771
772   if(source.fIp){
773     // we have the trackparam: assign or copy construct
774     if(fIp)*fIp = *source.fIp;
775     else fIp = new AliExternalTrackParam(*source.fIp);
776   }
777   else{
778     // no track param delete the old one
779     delete fIp;
780     fIp = 0;
781   }
782
783
784   if(source.fTPCInner){
785     // we have the trackparam: assign or copy construct
786     if(fTPCInner) *fTPCInner = *source.fTPCInner;
787     else fTPCInner = new AliExternalTrackParam(*source.fTPCInner);
788   }
789   else{
790     // no track param delete the old one
791     delete fTPCInner;
792     fTPCInner = 0;
793   }
794
795   if(source.fTPCdEdxInfo) {
796     if(fTPCdEdxInfo) *fTPCdEdxInfo = *source.fTPCdEdxInfo;
797     fTPCdEdxInfo = new AliTPCdEdxInfo(*source.fTPCdEdxInfo);
798   }
799
800   if(source.fOp){
801     // we have the trackparam: assign or copy construct
802     if(fOp) *fOp = *source.fOp;
803     else fOp = new AliExternalTrackParam(*source.fOp);
804   }
805   else{
806     // no track param delete the old one
807     delete fOp;
808     fOp = 0;
809   }
810
811   
812   if(source.fHMPIDp){
813     // we have the trackparam: assign or copy construct
814     if(fHMPIDp) *fHMPIDp = *source.fHMPIDp;
815     else fHMPIDp = new AliExternalTrackParam(*source.fHMPIDp);
816   }
817   else{
818     // no track param delete the old one
819     delete fHMPIDp;
820     fHMPIDp = 0;
821   }
822
823   // copy also the friend track 
824   // use copy constructor
825   if(source.fFriendTrack){
826     // we have the trackparam: assign or copy construct
827     delete fFriendTrack; fFriendTrack=new AliESDfriendTrack(*source.fFriendTrack);
828   }
829   else{
830     // no track param delete the old one
831     delete fFriendTrack; fFriendTrack= 0;
832   }
833
834   fTPCFitMap = source.fTPCFitMap; 
835   fTPCClusterMap = source.fTPCClusterMap; 
836   fTPCSharedMap  = source.fTPCSharedMap;  
837   // the simple stuff
838   fFlags    = source.fFlags; 
839   fID       = source.fID;             
840   fLabel    = source.fLabel;
841   fITSLabel = source.fITSLabel;
842   for(int i = 0; i< 12;++i){
843     fITSModule[i] = source.fITSModule[i];
844   }
845   fTPCLabel = source.fTPCLabel; 
846   fTRDLabel = source.fTRDLabel;
847   for(int i = 0; i< 3;++i){
848     fTOFLabel[i] = source.fTOFLabel[i];    
849   }
850   fTOFCalChannel = source.fTOFCalChannel;
851   fTOFindex      = source.fTOFindex;
852   fHMPIDqn       = source.fHMPIDqn;
853   fHMPIDcluIdx   = source.fHMPIDcluIdx; 
854   fCaloIndex    = source.fCaloIndex;
855   for (int i=kNITSchi2Std;i--;) fITSchi2Std[i] = source.fITSchi2Std[i];
856   for(int i = 0; i< 3;++i){
857     fKinkIndexes[i] = source.fKinkIndexes[i]; 
858     fV0Indexes[i]   = source.fV0Indexes[i]; 
859   }
860
861   for(int i = 0; i< AliPID::kSPECIES;++i){
862     fR[i]     = source.fR[i];
863     fITSr[i]  = source.fITSr[i];
864     fTPCr[i]  = source.fTPCr[i];
865     fTRDr[i]  = source.fTRDr[i];
866     fTOFr[i]  = source.fTOFr[i];
867     fHMPIDr[i] = source.fHMPIDr[i];
868     fTrackTime[i] = source.fTrackTime[i];  
869   }
870
871   fHMPIDtrkTheta = source.fHMPIDtrkTheta;
872   fHMPIDtrkPhi   = source.fHMPIDtrkPhi;
873   fHMPIDsignal   = source.fHMPIDsignal; 
874
875   
876   fTrackLength   = source. fTrackLength;
877   fdTPC  = source.fdTPC; 
878   fzTPC  = source.fzTPC; 
879   fCddTPC = source.fCddTPC;
880   fCdzTPC = source.fCdzTPC;
881   fCzzTPC = source.fCzzTPC;
882   fCchi2TPC = source.fCchi2TPC;
883
884   fD  = source.fD; 
885   fZ  = source.fZ; 
886   fCdd = source.fCdd;
887   fCdz = source.fCdz;
888   fCzz = source.fCzz;
889   fCchi2     = source.fCchi2;
890
891   fITSchi2   = source.fITSchi2;             
892   fTPCchi2   = source.fTPCchi2;            
893   fTPCchi2Iter1   = source.fTPCchi2Iter1;            
894   fTRDchi2   = source.fTRDchi2;      
895   fTOFchi2   = source.fTOFchi2;      
896   fHMPIDchi2 = source.fHMPIDchi2;      
897
898   fGlobalChi2 = source.fGlobalChi2;      
899
900   fITSsignal  = source.fITSsignal;     
901   for (Int_t i=0;i<4;i++) {fITSdEdxSamples[i]=source.fITSdEdxSamples[i];}
902   fTPCsignal  = source.fTPCsignal;     
903   fTPCsignalS = source.fTPCsignalS;    
904   for(int i = 0; i< 4;++i){
905     fTPCPoints[i] = source.fTPCPoints[i];  
906   }
907   fTRDsignal = source.fTRDsignal;
908
909   for(int i = 0;i < kTRDnPlanes;++i){
910     fTRDTimBin[i] = source.fTRDTimBin[i];   
911   }
912
913   if(fTRDnSlices)
914     delete[] fTRDslices;
915   fTRDslices=0;
916   fTRDnSlices=source.fTRDnSlices;
917   if (fTRDnSlices) {
918     fTRDslices=new Double32_t[fTRDnSlices];
919     for(int j = 0;j < fTRDnSlices;++j) fTRDslices[j] = source.fTRDslices[j];
920   }
921
922   fTRDQuality =   source.fTRDQuality;     
923   fTRDBudget  =   source.fTRDBudget;      
924   fTOFsignal  =   source.fTOFsignal;     
925   fTOFsignalToT = source.fTOFsignalToT;   
926   fTOFsignalRaw = source.fTOFsignalRaw;  
927   fTOFsignalDz  = source.fTOFsignalDz;      
928   fTOFsignalDx  = source.fTOFsignalDx;      
929   fTOFdeltaBC   = source.fTOFdeltaBC;
930   fTOFl0l1      = source.fTOFl0l1;
931  
932   for(int i = 0;i<10;++i){
933     fTOFInfo[i] = source.fTOFInfo[i];    
934   }
935
936   fHMPIDtrkX = source.fHMPIDtrkX; 
937   fHMPIDtrkY = source.fHMPIDtrkY; 
938   fHMPIDmipX = source.fHMPIDmipX;
939   fHMPIDmipY = source.fHMPIDmipY; 
940
941   fTPCncls    = source.fTPCncls;      
942   fTPCnclsF   = source.fTPCnclsF;     
943   fTPCsignalN = source.fTPCsignalN;   
944   fTPCnclsIter1    = source.fTPCnclsIter1;      
945   fTPCnclsFIter1   = source.fTPCnclsFIter1;     
946
947   fITSncls = source.fITSncls;       
948   fITSClusterMap = source.fITSClusterMap; 
949   fITSSharedMap = source.fITSSharedMap; 
950   fTRDncls   = source.fTRDncls;       
951   fTRDncls0  = source.fTRDncls0;      
952   fTRDntracklets  = source.fTRDntracklets; 
953   fVertexID = source.fVertexID;
954
955   fCacheNCrossedRows = source.fCacheNCrossedRows;
956   fCacheChi2TPCConstrainedVsGlobal = source.fCacheChi2TPCConstrainedVsGlobal;
957   fCacheChi2TPCConstrainedVsGlobalVertex = source.fCacheChi2TPCConstrainedVsGlobalVertex;
958
959   return *this;
960 }
961
962
963
964 void AliESDtrack::Copy(TObject &obj) const {
965   
966   // this overwrites the virtual TOBject::Copy()
967   // to allow run time copying without casting
968   // in AliESDEvent
969
970   if(this==&obj)return;
971   AliESDtrack *robj = dynamic_cast<AliESDtrack*>(&obj);
972   if(!robj)return; // not an AliESDtrack
973   *robj = *this;
974
975 }
976
977
978
979 void AliESDtrack::AddCalibObject(TObject * object){
980   //
981   // add calib object to the list
982   //
983   if (!fFriendTrack) fFriendTrack  = new AliESDfriendTrack;
984   if (!fFriendTrack) return;
985   fFriendTrack->AddCalibObject(object);
986 }
987
988 TObject *  AliESDtrack::GetCalibObject(Int_t index){
989   //
990   // return calib objct at given position
991   //
992   if (!fFriendTrack) return 0;
993   return fFriendTrack->GetCalibObject(index);
994 }
995
996
997 Bool_t AliESDtrack::FillTPCOnlyTrack(AliESDtrack &track){
998   
999   // Fills the information of the TPC-only first reconstruction pass
1000   // into the passed ESDtrack object. For consistency fTPCInner is also filled
1001   // again
1002
1003
1004
1005   // For data produced before r26675
1006   // RelateToVertexTPC was not properly called during reco
1007   // so you'll have to call it again, before FillTPCOnlyTrack
1008   //  Float_t p[2],cov[3];
1009   // track->GetImpactParametersTPC(p,cov); 
1010   // if(p[0]==0&&p[1]==0) // <- Default values
1011   //  track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
1012   
1013
1014   if(!fTPCInner)return kFALSE;
1015
1016   // fill the TPC track params to the global track parameters
1017   track.Set(fTPCInner->GetX(),fTPCInner->GetAlpha(),fTPCInner->GetParameter(),fTPCInner->GetCovariance());
1018   track.fD = fdTPC;
1019   track.fZ = fzTPC;
1020   track.fCdd = fCddTPC;
1021   track.fCdz = fCdzTPC;
1022   track.fCzz = fCzzTPC;
1023
1024   // copy the inner params
1025   if(track.fIp) *track.fIp = *fIp;
1026   else track.fIp = new AliExternalTrackParam(*fIp);
1027
1028   // copy the TPCinner parameters
1029   if(track.fTPCInner) *track.fTPCInner = *fTPCInner;
1030   else track.fTPCInner = new AliExternalTrackParam(*fTPCInner);
1031   track.fdTPC   = fdTPC;
1032   track.fzTPC   = fzTPC;
1033   track.fCddTPC = fCddTPC;
1034   track.fCdzTPC = fCdzTPC;
1035   track.fCzzTPC = fCzzTPC;
1036   track.fCchi2TPC = fCchi2TPC;
1037
1038   // copy all other TPC specific parameters
1039
1040   // replace label by TPC label
1041   track.fLabel    = fTPCLabel;
1042   track.fTPCLabel = fTPCLabel;
1043
1044   track.fTPCchi2 = fTPCchi2; 
1045   track.fTPCchi2Iter1 = fTPCchi2Iter1; 
1046   track.fTPCsignal = fTPCsignal;
1047   track.fTPCsignalS = fTPCsignalS;
1048   for(int i = 0;i<4;++i)track.fTPCPoints[i] = fTPCPoints[i];
1049
1050   track.fTPCncls    = fTPCncls;     
1051   track.fTPCnclsF   = fTPCnclsF;     
1052   track.fTPCsignalN =  fTPCsignalN;
1053   track.fTPCnclsIter1    = fTPCnclsIter1;     
1054   track.fTPCnclsFIter1   = fTPCnclsFIter1;     
1055
1056   // PID 
1057   for(int i=0;i<AliPID::kSPECIES;++i){
1058     track.fTPCr[i] = fTPCr[i];
1059     // combined PID is TPC only!
1060     track.fR[i] = fTPCr[i];
1061   }
1062   track.fTPCFitMap = fTPCFitMap;
1063   track.fTPCClusterMap = fTPCClusterMap;
1064   track.fTPCSharedMap = fTPCSharedMap;
1065
1066
1067   // reset the flags
1068   track.fFlags = kTPCin;
1069   track.fID    = fID;
1070
1071   track.fFlags |= fFlags & kTPCpid; //copy the TPCpid status flag
1072  
1073   for (Int_t i=0;i<3;i++) track.fKinkIndexes[i] = fKinkIndexes[i];
1074   
1075   return kTRUE;
1076     
1077 }
1078
1079 //_______________________________________________________________________
1080 void AliESDtrack::MakeMiniESDtrack(){
1081   // Resets everything except
1082   // fFlags: Reconstruction status flags 
1083   // fLabel: Track label
1084   // fID:  Unique ID of the track
1085   // Impact parameter information
1086   // fR[AliPID::kSPECIES]: combined "detector response probability"
1087   // Running track parameters in the base class (AliExternalTrackParam)
1088   
1089   fTrackLength = 0;
1090
1091   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i] = 0;
1092
1093   // Reset track parameters constrained to the primary vertex
1094   delete fCp;fCp = 0;
1095
1096   // Reset track parameters at the inner wall of TPC
1097   delete fIp;fIp = 0;
1098   delete fTPCInner;fTPCInner=0;
1099   // Reset track parameters at the inner wall of the TRD
1100   delete fOp;fOp = 0;
1101   // Reset track parameters at the HMPID
1102   delete fHMPIDp;fHMPIDp = 0;
1103
1104
1105   // Reset ITS track related information
1106   fITSchi2 = 0;
1107   fITSncls = 0;       
1108   fITSClusterMap=0;
1109   fITSSharedMap=0;
1110   fITSsignal = 0;     
1111   for (Int_t i=0;i<4;i++) fITSdEdxSamples[i] = 0.;
1112   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=0; 
1113   fITSLabel = 0;       
1114
1115   // Reset TPC related track information
1116   fTPCchi2 = 0;       
1117   fTPCchi2Iter1 = 0;       
1118   fTPCncls = 0;       
1119   fTPCnclsF = 0;       
1120   fTPCnclsIter1 = 0;       
1121   fTPCnclsFIter1 = 0;  
1122   fTPCFitMap = 0;       
1123   fTPCClusterMap = 0;  
1124   fTPCSharedMap = 0;  
1125   fTPCsignal= 0;      
1126   fTPCsignalS= 0;      
1127   fTPCsignalN= 0;      
1128   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0; 
1129   fTPCLabel=0;       
1130   for (Int_t i=0;i<4;i++) fTPCPoints[i] = 0;
1131   for (Int_t i=0; i<3;i++)   fKinkIndexes[i] = 0;
1132   for (Int_t i=0; i<3;i++)   fV0Indexes[i] = 0;
1133
1134   // Reset TRD related track information
1135   fTRDchi2 = 0;        
1136   fTRDncls = 0;       
1137   fTRDncls0 = 0;       
1138   fTRDsignal = 0;      
1139   for (Int_t i=0;i<kTRDnPlanes;i++) {
1140     fTRDTimBin[i]  = 0;
1141   }
1142   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0; 
1143   fTRDLabel = 0;       
1144   fTRDQuality  = 0;
1145   fTRDntracklets = 0;
1146   if(fTRDnSlices)
1147     delete[] fTRDslices;
1148   fTRDslices=0x0;
1149   fTRDnSlices=0;
1150   fTRDBudget  = 0;
1151
1152   // Reset TOF related track information
1153   fTOFchi2 = 0;        
1154   fTOFindex = -1;       
1155   fTOFsignal = 99999;      
1156   fTOFCalChannel = -1;
1157   fTOFsignalToT = 99999;
1158   fTOFsignalRaw = 99999;
1159   fTOFsignalDz = 999;
1160   fTOFsignalDx = 999;
1161   fTOFdeltaBC = 999;
1162   fTOFl0l1 = 999;
1163   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i] = 0;
1164   for (Int_t i=0;i<3;i++) fTOFLabel[i] = -1;
1165   for (Int_t i=0;i<10;i++) fTOFInfo[i] = 0;
1166
1167   // Reset HMPID related track information
1168   fHMPIDchi2 = 0;     
1169   fHMPIDqn = 0;     
1170   fHMPIDcluIdx = -1;     
1171   fHMPIDsignal = 0;     
1172   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i] = 0;
1173   fHMPIDtrkTheta = 0;     
1174   fHMPIDtrkPhi = 0;      
1175   fHMPIDtrkX = 0;     
1176   fHMPIDtrkY = 0;      
1177   fHMPIDmipX = 0;
1178   fHMPIDmipY = 0;
1179   fCaloIndex = kEMCALNoMatch;
1180
1181   // reset global track chi2
1182   fGlobalChi2 = 0;
1183
1184   fVertexID = -2; // an orphan track
1185
1186   delete fFriendTrack; fFriendTrack = 0;
1187
1188
1189 //_______________________________________________________________________
1190 Int_t AliESDtrack::GetPID() const 
1191 {
1192   // Returns the particle most probable id
1193   Int_t i;
1194   for (i=0; i<AliPID::kSPECIES-1; i++) if (fR[i] != fR[i+1]) break;
1195   //
1196   if (i == AliPID::kSPECIES-1) return AliPID::kPion;  // If all the probabilities are equal, return the pion mass
1197   //
1198   Float_t max=0.;
1199   Int_t k=-1;
1200   for (i=0; i<AliPID::kSPECIES; i++) if (fR[i]>max) {k=i; max=fR[i];}
1201   //
1202   if (k==0) { // dE/dx "crossing points" in the TPC
1203     Double_t p=GetP();
1204     if ((p>0.38)&&(p<0.48))
1205       if (fR[0]<fR[3]*10.) return AliPID::kKaon;
1206     if ((p>0.75)&&(p<0.85))
1207       if (fR[0]<fR[4]*10.) return AliPID::kProton;
1208     return AliPID::kElectron;
1209   }
1210   if (k==1) return AliPID::kMuon; 
1211   if (k==2||k==-1) return AliPID::kPion;
1212   if (k==3) return AliPID::kKaon;
1213   if (k==4) return AliPID::kProton;
1214   AliWarning("Undefined PID !");
1215   return AliPID::kPion;
1216 }
1217
1218 //_______________________________________________________________________
1219 Int_t AliESDtrack::GetTOFBunchCrossing(Double_t b) const 
1220 {
1221   // Returns the number of bunch crossings after trigger (assuming 25ns spacing)
1222   const double kSpacing = 25e3; // min interbanch spacing
1223   const double kShift = 0;
1224   Int_t bcid = kTOFBCNA; // defualt one
1225   if (!IsOn(kTOFout) || !IsOn(kESDpid)) return bcid; // no info
1226   //
1227   double tdif = fTOFsignal;
1228   if (IsOn(kTIME)) { // integrated time info is there
1229     int pid = GetPID();
1230     tdif -= fTrackTime[pid];
1231   }
1232   else { // assume integrated time info from TOF radius and momentum
1233     const double kRTOF = 385.;
1234     const double kCSpeed = 3.e-2; // cm/ps
1235     double p = GetP();
1236     if (p<0.01) return bcid;
1237     double m = GetMass();
1238     double curv = GetC(b);
1239     double path = TMath::Abs(curv)>kAlmost0 ? // account for curvature
1240       2./curv*TMath::ASin(kRTOF*curv/2.)*TMath::Sqrt(1.+GetTgl()*GetTgl()) : kRTOF;
1241     tdif -= path/kCSpeed*TMath::Sqrt(1.+m*m/(p*p));
1242   }
1243   bcid = TMath::Nint((tdif - kShift)/kSpacing);
1244   return bcid;
1245 }
1246
1247 //______________________________________________________________________________
1248 Double_t AliESDtrack::M() const
1249 {
1250   // Returns the assumed mass
1251   // (the pion mass, if the particle can't be identified properly).
1252   static Bool_t printerr=kTRUE;
1253   if (printerr) {
1254      AliWarning("WARNING !!! ... THIS WILL BE PRINTED JUST ONCE !!!");
1255      printerr = kFALSE;
1256      AliWarning("This is the ESD mass. Use it with care !"); 
1257   }
1258   return GetMass(); 
1259 }
1260   
1261 //______________________________________________________________________________
1262 Double_t AliESDtrack::E() const
1263 {
1264   // Returns the energy of the particle given its assumed mass.
1265   // Assumes the pion mass if the particle can't be identified properly.
1266   
1267   Double_t m = M();
1268   Double_t p = P();
1269   return TMath::Sqrt(p*p + m*m);
1270 }
1271
1272 //______________________________________________________________________________
1273 Double_t AliESDtrack::Y() const
1274 {
1275   // Returns the rapidity of a particle given its assumed mass.
1276   // Assumes the pion mass if the particle can't be identified properly.
1277   
1278   Double_t e = E();
1279   Double_t pz = Pz();
1280   if (e != TMath::Abs(pz)) { // energy was not equal to pz
1281     return 0.5*TMath::Log((e+pz)/(e-pz));
1282   } else { // energy was equal to pz
1283     return -999.;
1284   }
1285 }
1286
1287 //_______________________________________________________________________
1288 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
1289   //
1290   // This function updates track's running parameters 
1291   //
1292   Bool_t rc=kTRUE;
1293
1294   SetStatus(flags);
1295   fLabel=t->GetLabel();
1296
1297   if (t->IsStartedTimeIntegral()) {
1298     SetStatus(kTIME);
1299     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
1300     SetIntegratedLength(t->GetIntegratedLength());
1301   }
1302
1303   Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1304   if (fFriendTrack) {
1305   if (flags==kITSout) fFriendTrack->SetITSOut(*t);
1306   if (flags==kTPCout) fFriendTrack->SetTPCOut(*t);
1307   if (flags==kTRDrefit) fFriendTrack->SetTRDIn(*t);
1308   }
1309   
1310   switch (flags) {
1311     
1312   case kITSin: 
1313     fITSchi2Std[0] = t->GetChi2();
1314     //
1315   case kITSout: 
1316     fITSchi2Std[1] = t->GetChi2();
1317   case kITSrefit:
1318     {
1319     fITSchi2Std[2] = t->GetChi2();
1320     fITSClusterMap=0;
1321     fITSncls=t->GetNumberOfClusters();
1322     if (fFriendTrack) {
1323     Int_t* indexITS = new Int_t[AliESDfriendTrack::kMaxITScluster];
1324     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
1325         indexITS[i]=t->GetClusterIndex(i);
1326
1327         if (i<fITSncls) {
1328           Int_t l=(indexITS[i] & 0xf0000000) >> 28;
1329            SETBIT(fITSClusterMap,l);                 
1330         }
1331     }
1332     fFriendTrack->SetITSIndices(indexITS,AliESDfriendTrack::kMaxITScluster);
1333     delete [] indexITS;
1334     }
1335
1336     fITSchi2=t->GetChi2();
1337     fITSsignal=t->GetPIDsignal();
1338     fITSLabel = t->GetLabel();
1339     // keep in fOp the parameters outside ITS for ITS stand-alone tracks 
1340     if (flags==kITSout) { 
1341       if (!fOp) fOp=new AliExternalTrackParam(*t);
1342       else 
1343         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1344     }   
1345     }
1346     break;
1347     
1348   case kTPCin: case kTPCrefit:
1349     {
1350     fTPCLabel = t->GetLabel();
1351     if (flags==kTPCin)  {
1352       fTPCInner=new AliExternalTrackParam(*t); 
1353       fTPCnclsIter1=t->GetNumberOfClusters();    
1354       fTPCchi2Iter1=t->GetChi2();
1355     }
1356     if (!fIp) fIp=new AliExternalTrackParam(*t);
1357     else 
1358       fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1359     }
1360     // Intentionally no break statement; need to set general TPC variables as well
1361   case kTPCout:
1362     {
1363     if (flags & kTPCout){
1364       if (!fOp) fOp=new AliExternalTrackParam(*t);
1365       else 
1366         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1367     }
1368     fTPCncls=t->GetNumberOfClusters();    
1369     fTPCchi2=t->GetChi2();
1370     
1371     if (fFriendTrack) {  // Copy cluster indices
1372       Int_t* indexTPC = new Int_t[AliESDfriendTrack::kMaxTPCcluster];
1373       for (Int_t i=0;i<AliESDfriendTrack::kMaxTPCcluster;i++)         
1374         indexTPC[i]=t->GetClusterIndex(i);
1375       fFriendTrack->SetTPCIndices(indexTPC,AliESDfriendTrack::kMaxTPCcluster);
1376       delete [] indexTPC;
1377     }
1378     fTPCsignal=t->GetPIDsignal();
1379     }
1380     break;
1381
1382   case kTRDin: case kTRDrefit:
1383     break;
1384   case kTRDout:
1385     {
1386     fTRDLabel = t->GetLabel(); 
1387     fTRDchi2  = t->GetChi2();
1388     fTRDncls  = t->GetNumberOfClusters();
1389     if (fFriendTrack) {
1390       Int_t* indexTRD = new Int_t[AliESDfriendTrack::kMaxTRDcluster];
1391       for (Int_t i=0;i<AliESDfriendTrack::kMaxTRDcluster;i++) indexTRD[i]=-2;
1392       for (Int_t i=0;i<6;i++) indexTRD[i]=t->GetTrackletIndex(i);
1393       fFriendTrack->SetTRDIndices(indexTRD,AliESDfriendTrack::kMaxTRDcluster);
1394       delete [] indexTRD;
1395     }    
1396     
1397     fTRDsignal=t->GetPIDsignal();
1398     }
1399     break;
1400   case kTRDbackup:
1401     if (!fOp) fOp=new AliExternalTrackParam(*t);
1402     else 
1403       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1404     fTRDncls0 = t->GetNumberOfClusters(); 
1405     break;
1406   case kTOFin: 
1407     break;
1408   case kTOFout: 
1409     break;
1410   case kTRDStop:
1411     break;
1412   case kHMPIDout:
1413   if (!fHMPIDp) fHMPIDp=new AliExternalTrackParam(*t);
1414     else 
1415       fHMPIDp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1416     break;
1417   default: 
1418     AliError("Wrong flag !");
1419     return kFALSE;
1420   }
1421
1422   return rc;
1423 }
1424
1425 //_______________________________________________________________________
1426 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
1427   //---------------------------------------------------------------------
1428   // This function returns external representation of the track parameters
1429   //---------------------------------------------------------------------
1430   x=GetX();
1431   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
1432 }
1433
1434 //_______________________________________________________________________
1435 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
1436   //---------------------------------------------------------------------
1437   // This function returns external representation of the cov. matrix
1438   //---------------------------------------------------------------------
1439   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
1440 }
1441
1442 //_______________________________________________________________________
1443 Bool_t AliESDtrack::GetConstrainedExternalParameters
1444                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1445   //---------------------------------------------------------------------
1446   // This function returns the constrained external track parameters
1447   //---------------------------------------------------------------------
1448   if (!fCp) return kFALSE;
1449   alpha=fCp->GetAlpha();
1450   x=fCp->GetX();
1451   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
1452   return kTRUE;
1453 }
1454
1455 //_______________________________________________________________________
1456 Bool_t 
1457 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
1458   //---------------------------------------------------------------------
1459   // This function returns the constrained external cov. matrix
1460   //---------------------------------------------------------------------
1461   if (!fCp) return kFALSE;
1462   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
1463   return kTRUE;
1464 }
1465
1466 Bool_t
1467 AliESDtrack::GetInnerExternalParameters
1468                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1469   //---------------------------------------------------------------------
1470   // This function returns external representation of the track parameters 
1471   // at the inner layer of TPC
1472   //---------------------------------------------------------------------
1473   if (!fIp) return kFALSE;
1474   alpha=fIp->GetAlpha();
1475   x=fIp->GetX();
1476   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
1477   return kTRUE;
1478 }
1479
1480 Bool_t 
1481 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
1482  //---------------------------------------------------------------------
1483  // This function returns external representation of the cov. matrix 
1484  // at the inner layer of TPC
1485  //---------------------------------------------------------------------
1486   if (!fIp) return kFALSE;
1487   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
1488   return kTRUE;
1489 }
1490
1491 void 
1492 AliESDtrack::SetOuterParam(const AliExternalTrackParam *p, ULong_t flags) {
1493   //
1494   // This is a direct setter for the outer track parameters
1495   //
1496   SetStatus(flags);
1497   if (fOp) delete fOp;
1498   fOp=new AliExternalTrackParam(*p);
1499 }
1500
1501 void 
1502 AliESDtrack::SetOuterHmpParam(const AliExternalTrackParam *p, ULong_t flags) {
1503   //
1504   // This is a direct setter for the outer track parameters
1505   //
1506   SetStatus(flags);
1507   if (fHMPIDp) delete fHMPIDp;
1508   fHMPIDp=new AliExternalTrackParam(*p);
1509 }
1510
1511 Bool_t 
1512 AliESDtrack::GetOuterExternalParameters
1513                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1514   //---------------------------------------------------------------------
1515   // This function returns external representation of the track parameters 
1516   // at the inner layer of TRD
1517   //---------------------------------------------------------------------
1518   if (!fOp) return kFALSE;
1519   alpha=fOp->GetAlpha();
1520   x=fOp->GetX();
1521   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
1522   return kTRUE;
1523 }
1524
1525 Bool_t 
1526 AliESDtrack::GetOuterHmpExternalParameters
1527                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1528   //---------------------------------------------------------------------
1529   // This function returns external representation of the track parameters 
1530   // at the inner layer of TRD
1531   //---------------------------------------------------------------------
1532   if (!fHMPIDp) return kFALSE;
1533   alpha=fHMPIDp->GetAlpha();
1534   x=fHMPIDp->GetX();
1535   for (Int_t i=0; i<5; i++) p[i]=fHMPIDp->GetParameter()[i];
1536   return kTRUE;
1537 }
1538
1539 Bool_t 
1540 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
1541  //---------------------------------------------------------------------
1542  // This function returns external representation of the cov. matrix 
1543  // at the inner layer of TRD
1544  //---------------------------------------------------------------------
1545   if (!fOp) return kFALSE;
1546   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
1547   return kTRUE;
1548 }
1549
1550 Bool_t 
1551 AliESDtrack::GetOuterHmpExternalCovariance(Double_t cov[15]) const {
1552  //---------------------------------------------------------------------
1553  // This function returns external representation of the cov. matrix 
1554  // at the inner layer of TRD
1555  //---------------------------------------------------------------------
1556   if (!fHMPIDp) return kFALSE;
1557   for (Int_t i=0; i<15; i++) cov[i]=fHMPIDp->GetCovariance()[i];
1558   return kTRUE;
1559 }
1560
1561 Int_t AliESDtrack::GetNcls(Int_t idet) const
1562 {
1563   // Get number of clusters by subdetector index
1564   //
1565   Int_t ncls = 0;
1566   switch(idet){
1567   case 0:
1568     ncls = fITSncls;
1569     break;
1570   case 1:
1571     ncls = fTPCncls;
1572     break;
1573   case 2:
1574     ncls = fTRDncls;
1575     break;
1576   case 3:
1577     if (fTOFindex != -1)
1578       ncls = 1;
1579     break;
1580   case 4: //PHOS
1581     break;
1582   case 5: //HMPID
1583     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1584       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1585         ncls = 1;
1586       }
1587     }    
1588     break;
1589   default:
1590     break;
1591   }
1592   return ncls;
1593 }
1594
1595 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
1596 {
1597   // Get cluster index array by subdetector index
1598   //
1599   Int_t ncls = 0;
1600   switch(idet){
1601   case 0:
1602     ncls = GetITSclusters(idx);
1603     break;
1604   case 1:
1605     ncls = GetTPCclusters(idx);
1606     break;
1607   case 2:
1608     ncls = GetTRDclusters(idx);
1609     break;
1610   case 3:
1611     if (fTOFindex != -1) {
1612       idx[0] = fTOFindex;
1613       ncls = 1;
1614     }
1615     break;
1616   case 4: //PHOS
1617     break;
1618   case 5:
1619     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1620       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1621         idx[0] = GetHMPIDcluIdx();
1622         ncls = 1;
1623       }
1624     }    
1625     break;
1626   case 6: //EMCAL
1627     break;
1628   default:
1629     break;
1630   }
1631   return ncls;
1632 }
1633
1634 //_______________________________________________________________________
1635 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
1636   // Returns the array with integrated times for each particle hypothesis
1637   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
1638 }
1639
1640 //_______________________________________________________________________
1641 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
1642   // Sets the array with integrated times for each particle hypotesis
1643   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
1644 }
1645
1646 //_______________________________________________________________________
1647 void AliESDtrack::SetITSpid(const Double_t *p) {
1648   // Sets values for the probability of each particle type (in ITS)
1649   SetPIDValues(fITSr,p,AliPID::kSPECIES);
1650   SetStatus(AliESDtrack::kITSpid);
1651 }
1652
1653 //_______________________________________________________________________
1654 void AliESDtrack::GetITSpid(Double_t *p) const {
1655   // Gets the probability of each particle type (in ITS)
1656   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
1657 }
1658
1659 //_______________________________________________________________________
1660 Char_t AliESDtrack::GetITSclusters(Int_t *idx) const {
1661   //---------------------------------------------------------------------
1662   // This function returns indices of the assgined ITS clusters 
1663   //---------------------------------------------------------------------
1664   if (idx && fFriendTrack) {
1665     Int_t *index=fFriendTrack->GetITSindices();
1666     for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) {
1667       if ( (i>=fITSncls) && (i<6) ) idx[i]=-1;
1668       else {
1669         if (index) {
1670           idx[i]=index[i];
1671         }
1672         else idx[i]= -2;
1673       }
1674     }
1675   }
1676   return fITSncls;
1677 }
1678
1679 //_______________________________________________________________________
1680 Bool_t AliESDtrack::GetITSModuleIndexInfo(Int_t ilayer,Int_t &idet,Int_t &status,
1681                                          Float_t &xloc,Float_t &zloc) const {
1682   //----------------------------------------------------------------------
1683   // This function encodes in the module number also the status of cluster association
1684   // "status" can have the following values: 
1685   // 1 "found" (cluster is associated), 
1686   // 2 "dead" (module is dead from OCDB), 
1687   // 3 "skipped" (module or layer forced to be skipped),
1688   // 4 "outinz" (track out of z acceptance), 
1689   // 5 "nocls" (no clusters in the road), 
1690   // 6 "norefit" (cluster rejected during refit), 
1691   // 7 "deadzspd" (holes in z in SPD)
1692   // Also given are the coordinates of the crossing point of track and module
1693   // (in the local module ref. system)
1694   // WARNING: THIS METHOD HAS TO BE SYNCHRONIZED WITH AliITStrackV2::GetModuleIndexInfo()!
1695   //----------------------------------------------------------------------
1696
1697   if(fITSModule[ilayer]==-1) {
1698     idet = -1;
1699     status=0;
1700     xloc=-99.; zloc=-99.;
1701     return kFALSE;
1702   }
1703
1704   Int_t module = fITSModule[ilayer];
1705
1706   idet = Int_t(module/1000000);
1707
1708   module -= idet*1000000;
1709
1710   status = Int_t(module/100000);
1711
1712   module -= status*100000;
1713
1714   Int_t signs = Int_t(module/10000);
1715
1716   module-=signs*10000;
1717
1718   Int_t xInt = Int_t(module/100);
1719   module -= xInt*100;
1720
1721   Int_t zInt = module;
1722
1723   if(signs==1) { xInt*=1; zInt*=1; }
1724   if(signs==2) { xInt*=1; zInt*=-1; }
1725   if(signs==3) { xInt*=-1; zInt*=1; }
1726   if(signs==4) { xInt*=-1; zInt*=-1; }
1727
1728   xloc = 0.1*(Float_t)xInt;
1729   zloc = 0.1*(Float_t)zInt;
1730
1731   if(status==4) idet = -1;
1732
1733   return kTRUE;
1734 }
1735
1736 //_______________________________________________________________________
1737 UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
1738   //---------------------------------------------------------------------
1739   // This function returns indices of the assgined ITS clusters 
1740   //---------------------------------------------------------------------
1741   if (idx && fFriendTrack) {
1742     Int_t *index=fFriendTrack->GetTPCindices();
1743
1744     if (index){
1745       for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
1746     }
1747     else {
1748       for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=-2;
1749     }
1750   }
1751   return fTPCncls;
1752 }
1753
1754 //_______________________________________________________________________
1755 Float_t AliESDtrack::GetTPCCrossedRows() const
1756 {
1757   // This function calls GetTPCClusterInfo with some default parameters which are used in the track selection and caches the outcome
1758   // because GetTPCClusterInfo is quite time-consuming
1759   
1760   if (fCacheNCrossedRows > -1)
1761     return fCacheNCrossedRows;
1762   
1763   fCacheNCrossedRows = GetTPCClusterInfo(2, 1);
1764   return fCacheNCrossedRows;
1765 }
1766
1767 //_______________________________________________________________________
1768 Float_t AliESDtrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1, Int_t bitType ) const
1769 {
1770   //
1771   // TPC cluster information
1772   // type 0: get fraction of found/findable clusters with neighbourhood definition
1773   //      1: findable clusters with neighbourhood definition
1774   //      2: found clusters
1775   // bitType:
1776   //      0 - all cluster used
1777   //      1 - clusters  used for the kalman update
1778   // definition of findable clusters:
1779   //            a cluster is defined as findable if there is another cluster
1780   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1781   //           effects with a very simple algorithm.
1782   //
1783
1784   
1785   Int_t found=0;
1786   Int_t findable=0;
1787   Int_t last=-nNeighbours;
1788   const TBits & clusterMap = (bitType%2==0) ? fTPCClusterMap : fTPCFitMap;
1789   
1790   Int_t upperBound=clusterMap.GetNbits();
1791   if (upperBound>row1) upperBound=row1;
1792   for (Int_t i=row0; i<upperBound; ++i){
1793     //look to current row
1794     if (clusterMap[i]) {
1795       last=i;
1796       ++found;
1797       ++findable;
1798       continue;
1799     }
1800     //look to nNeighbours before
1801     if ((i-last)<=nNeighbours) {
1802       ++findable;
1803       continue;
1804     }
1805     //look to nNeighbours after
1806     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1807       if (clusterMap[j]){
1808         ++findable;
1809         break;
1810       }
1811     }
1812   }
1813   if (type==2) return found;
1814   if (type==1) return findable;
1815   
1816   if (type==0){
1817     Float_t fraction=0;
1818     if (findable>0) 
1819       fraction=(Float_t)found/(Float_t)findable;
1820     else 
1821       fraction=0;
1822     return fraction;
1823   }  
1824   return 0;  // undefined type - default value
1825 }
1826
1827 //_______________________________________________________________________
1828 Float_t AliESDtrack::GetTPCClusterDensity(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1, Int_t bitType ) const
1829 {
1830   //
1831   // TPC cluster density -  only rows where signal before and after given row are used
1832   //                     -  slower function
1833   // type 0: get fraction of found/findable clusters with neighbourhood definition
1834   //      1: findable clusters with neighbourhood definition
1835   //      2: found clusters
1836   // bitType:
1837   //      0 - all cluster used
1838   //      1 - clusters  used for the kalman update
1839   // definition of findable clusters:
1840   //            a cluster is defined as findable if there is another cluster
1841   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1842   //           effects with a very simple algorithm.
1843   //  
1844   Int_t found=0;
1845   Int_t findable=0;
1846   //  Int_t last=-nNeighbours;
1847   const TBits & clusterMap = (bitType%2==0) ? fTPCClusterMap : fTPCFitMap;
1848   Int_t upperBound=clusterMap.GetNbits();
1849   if (upperBound>row1) upperBound=row1;
1850   for (Int_t i=row0; i<upperBound; ++i){
1851     Bool_t isUp=kFALSE;
1852     Bool_t isDown=kFALSE;
1853     for (Int_t idelta=1; idelta<=nNeighbours; idelta++){
1854       if (i-idelta>=0 && clusterMap[i-idelta]) isDown=kTRUE;
1855       if (i+idelta<upperBound && clusterMap[i+idelta]) isUp=kTRUE;
1856     }
1857     if (isUp&&isDown){
1858       ++findable;
1859       if (clusterMap[i]) ++found;
1860     }
1861   }
1862   if (type==2) return found;
1863   if (type==1) return findable;
1864   
1865   if (type==0){
1866     Float_t fraction=0;
1867     if (findable>0) 
1868       fraction=(Float_t)found/(Float_t)findable;
1869     else 
1870       fraction=0;
1871     return fraction;
1872   }  
1873   return 0;  // undefined type - default value
1874 }
1875
1876
1877
1878
1879 //_______________________________________________________________________
1880 Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
1881   //
1882   // GetDensity of the clusters on given region between row0 and row1
1883   // Dead zone effect takin into acoount
1884   //
1885   if (!fFriendTrack) return 0.0;
1886   Int_t good  = 0;
1887   Int_t found = 0;
1888   //  
1889   Int_t *index=fFriendTrack->GetTPCindices();
1890   for (Int_t i=row0;i<=row1;i++){     
1891     Int_t idx = index[i];
1892     if (idx!=-1)  good++;             // track outside of dead zone
1893     if (idx>0)    found++;
1894   }
1895   Float_t density=0.5;
1896   if (good>TMath::Max((row1-row0)*0.5,0.0)) density = Float_t(found)/Float_t(good);
1897   return density;
1898 }
1899
1900 //_______________________________________________________________________
1901 void AliESDtrack::SetTPCpid(const Double_t *p) {  
1902   // Sets values for the probability of each particle type (in TPC)
1903   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
1904   SetStatus(AliESDtrack::kTPCpid);
1905 }
1906
1907 //_______________________________________________________________________
1908 void AliESDtrack::GetTPCpid(Double_t *p) const {
1909   // Gets the probability of each particle type (in TPC)
1910   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
1911 }
1912
1913 //_______________________________________________________________________
1914 UChar_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
1915   //---------------------------------------------------------------------
1916   // This function returns indices of the assgined TRD clusters 
1917   //---------------------------------------------------------------------
1918   if (idx && fFriendTrack) {
1919     Int_t *index=fFriendTrack->GetTRDindices();
1920
1921     if (index) {
1922       for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
1923     }
1924     else {
1925       for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=-2;
1926     }
1927   }
1928   return fTRDncls;
1929 }
1930
1931 //_______________________________________________________________________
1932 UChar_t AliESDtrack::GetTRDtracklets(Int_t *idx) const {
1933 //
1934 // This function returns the number of TRD tracklets used in tracking
1935 // and it fills the indices of these tracklets in the array "idx" as they 
1936 // are registered in the TRD track list. 
1937 // 
1938 // Caution :
1939 //   1. The idx array has to be allocated with a size >= AliESDtrack::kTRDnPlanes
1940 //   2. The idx array store not only the index but also the layer of the tracklet. 
1941 //      Therefore tracks with TRD gaps contain default values for indices [-1] 
1942
1943   if (!fFriendTrack) return 0;
1944   if (!idx) return GetTRDntracklets();
1945   Int_t *index=fFriendTrack->GetTRDindices();
1946   Int_t n = 0;
1947   for (Int_t i=0; i<kTRDnPlanes; i++){ 
1948     if (index){
1949       if(index[i]>=0) n++;
1950       idx[i]=index[i];
1951     }
1952     else idx[i] = -2;
1953   }
1954   return n;
1955 }
1956
1957 //_______________________________________________________________________
1958 void AliESDtrack::SetTRDpid(const Double_t *p) {  
1959   // Sets values for the probability of each particle type (in TRD)
1960   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
1961   SetStatus(AliESDtrack::kTRDpid);
1962 }
1963
1964 //_______________________________________________________________________
1965 void AliESDtrack::GetTRDpid(Double_t *p) const {
1966   // Gets the probability of each particle type (in TRD)
1967   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
1968 }
1969
1970 //_______________________________________________________________________
1971 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
1972 {
1973   // Sets the probability of particle type iSpecies to p (in TRD)
1974   fTRDr[iSpecies] = p;
1975 }
1976
1977 Double_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
1978 {
1979   // Returns the probability of particle type iSpecies (in TRD)
1980   return fTRDr[iSpecies];
1981 }
1982
1983 //____________________________________________________
1984 Int_t AliESDtrack::GetNumberOfTRDslices() const 
1985 {
1986   // built in backward compatibility
1987   Int_t idx = fTRDnSlices - (kTRDnPlanes<<1);
1988   return idx<18 ? fTRDnSlices/kTRDnPlanes : idx/kTRDnPlanes;
1989 }
1990
1991 //____________________________________________________
1992 Double_t AliESDtrack::GetTRDmomentum(Int_t plane, Double_t *sp) const
1993 {
1994 //Returns momentum estimation and optional its error (sp)
1995 // in TRD layer "plane".
1996
1997   if (!fTRDnSlices) {
1998     AliDebug(2, "No TRD info allocated for this track.");
1999     return -1.;
2000   }
2001   if ((plane<0) || (plane>=kTRDnPlanes)) {
2002     AliWarning(Form("Request for TRD plane[%d] outside range.", plane)); 
2003     return -1.;
2004   }
2005
2006   Int_t idx = fTRDnSlices-(kTRDnPlanes<<1)+plane;
2007   // Protection for backward compatibility
2008   if(idx<(GetNumberOfTRDslices()*kTRDnPlanes)) return -1.;
2009
2010   if(sp) (*sp) = fTRDslices[idx+kTRDnPlanes];
2011   return fTRDslices[idx];
2012 }
2013
2014 //____________________________________________________
2015 Double_t  AliESDtrack::GetTRDslice(Int_t plane, Int_t slice) const {
2016   //Gets the charge from the slice of the plane
2017
2018   if(!fTRDslices) {
2019     //AliError("No TRD slices allocated for this track !");
2020     return -1.;
2021   }
2022   if ((plane<0) || (plane>=kTRDnPlanes)) {
2023     AliError("Info for TRD plane not available !");
2024     return -1.;
2025   }
2026   Int_t ns=GetNumberOfTRDslices();
2027   if ((slice<-1) || (slice>=ns)) {
2028     //AliError("Wrong TRD slice !");  
2029     return -1.;
2030   }
2031
2032   if(slice>=0) return fTRDslices[plane*ns + slice];
2033
2034   // return average of the dEdx measurements
2035   Double_t q=0.; Double32_t *s = &fTRDslices[plane*ns];
2036   for (Int_t i=0; i<ns; i++, s++) if((*s)>0.) q+=(*s);
2037   return q/ns;
2038 }
2039
2040 //____________________________________________________
2041 void  AliESDtrack::SetNumberOfTRDslices(Int_t n) {
2042   //Sets the number of slices used for PID 
2043   if (fTRDnSlices) return;
2044
2045   fTRDnSlices=n;
2046   fTRDslices=new Double32_t[fTRDnSlices];
2047   
2048   // set-up correctly the allocated memory
2049   memset(fTRDslices, 0, n*sizeof(Double32_t));
2050   for (Int_t i=GetNumberOfTRDslices(); i--;) fTRDslices[i]=-1.;
2051 }
2052
2053 //____________________________________________________
2054 void  AliESDtrack::SetTRDslice(Double_t q, Int_t plane, Int_t slice) {
2055   //Sets the charge q in the slice of the plane
2056   if(!fTRDslices) {
2057     AliError("No TRD slices allocated for this track !");
2058     return;
2059   }
2060   if ((plane<0) || (plane>=kTRDnPlanes)) {
2061     AliError("Info for TRD plane not allocated !");
2062     return;
2063   }
2064   Int_t ns=GetNumberOfTRDslices();
2065   if ((slice<0) || (slice>=ns)) {
2066     AliError("Wrong TRD slice !");
2067     return;
2068   }
2069   Int_t n=plane*ns + slice;
2070   fTRDslices[n]=q;
2071 }
2072
2073
2074 //____________________________________________________
2075 void AliESDtrack::SetTRDmomentum(Double_t p, Int_t plane, Double_t *sp)
2076 {
2077   if(!fTRDslices) {
2078     AliError("No TRD slices allocated for this track !");
2079     return;
2080   }
2081   if ((plane<0) || (plane>=kTRDnPlanes)) {
2082     AliError("Info for TRD plane not allocated !");
2083     return;
2084   }
2085
2086   Int_t idx = fTRDnSlices-(kTRDnPlanes<<1)+plane;
2087   // Protection for backward compatibility
2088   if(idx<GetNumberOfTRDslices()*kTRDnPlanes) return;
2089
2090   if(sp) fTRDslices[idx+kTRDnPlanes] = (*sp);
2091   fTRDslices[idx] = p;
2092 }
2093
2094
2095 //_______________________________________________________________________
2096 void AliESDtrack::SetTOFpid(const Double_t *p) {  
2097   // Sets the probability of each particle type (in TOF)
2098   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
2099   SetStatus(AliESDtrack::kTOFpid);
2100 }
2101
2102 //_______________________________________________________________________
2103 void AliESDtrack::SetTOFLabel(const Int_t *p) {  
2104   // Sets  (in TOF)
2105   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
2106 }
2107
2108 //_______________________________________________________________________
2109 void AliESDtrack::GetTOFpid(Double_t *p) const {
2110   // Gets probabilities of each particle type (in TOF)
2111   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
2112 }
2113
2114 //_______________________________________________________________________
2115 void AliESDtrack::GetTOFLabel(Int_t *p) const {
2116   // Gets (in TOF)
2117   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
2118 }
2119
2120 //_______________________________________________________________________
2121 void AliESDtrack::GetTOFInfo(Float_t *info) const {
2122   // Gets (in TOF)
2123   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
2124 }
2125
2126 //_______________________________________________________________________
2127 void AliESDtrack::SetTOFInfo(Float_t*info) {
2128   // Gets (in TOF)
2129   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
2130 }
2131
2132
2133
2134 //_______________________________________________________________________
2135 void AliESDtrack::SetHMPIDpid(const Double_t *p) {  
2136   // Sets the probability of each particle type (in HMPID)
2137   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
2138   SetStatus(AliESDtrack::kHMPIDpid);
2139 }
2140
2141 //_______________________________________________________________________
2142 void  AliESDtrack::SetTPCdEdxInfo(AliTPCdEdxInfo * dEdxInfo){ 
2143   if(fTPCdEdxInfo) delete fTPCdEdxInfo;
2144   fTPCdEdxInfo = dEdxInfo; 
2145 }
2146
2147 //_______________________________________________________________________
2148 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
2149   // Gets probabilities of each particle type (in HMPID)
2150   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
2151 }
2152
2153
2154
2155 //_______________________________________________________________________
2156 void AliESDtrack::SetESDpid(const Double_t *p) {  
2157   // Sets the probability of each particle type for the ESD track
2158   SetPIDValues(fR,p,AliPID::kSPECIES);
2159   SetStatus(AliESDtrack::kESDpid);
2160 }
2161
2162 //_______________________________________________________________________
2163 void AliESDtrack::GetESDpid(Double_t *p) const {
2164   // Gets probability of each particle type for the ESD track
2165   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
2166 }
2167
2168 //_______________________________________________________________________
2169 Bool_t AliESDtrack::RelateToVertexTPC(const AliESDVertex *vtx, 
2170 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
2171   //
2172   // Try to relate the TPC-only track parameters to the vertex "vtx", 
2173   // if the (rough) transverse impact parameter is not bigger then "maxd". 
2174   //            Magnetic field is "b" (kG).
2175   //
2176   // a) The TPC-only paramters are extapolated to the DCA to the vertex.
2177   // b) The impact parameters and their covariance matrix are calculated.
2178   // c) An attempt to constrain the TPC-only params to the vertex is done.
2179   //    The constrained params are returned via "cParam".
2180   //
2181   // In the case of success, the returned value is kTRUE
2182   // otherwise, it's kFALSE)
2183   // 
2184
2185   if (!fTPCInner) return kFALSE;
2186   if (!vtx) return kFALSE;
2187
2188   Double_t dz[2],cov[3];
2189   if (!fTPCInner->PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
2190
2191   fdTPC = dz[0];
2192   fzTPC = dz[1];  
2193   fCddTPC = cov[0];
2194   fCdzTPC = cov[1];
2195   fCzzTPC = cov[2];
2196   
2197   Double_t covar[6]; vtx->GetCovMatrix(covar);
2198   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2199   Double_t c[3]={covar[2],0.,covar[5]};
2200
2201   Double_t chi2=GetPredictedChi2(p,c);
2202   if (chi2>kVeryBig) return kFALSE;
2203
2204   fCchi2TPC=chi2;
2205
2206   if (!cParam) return kTRUE;
2207
2208   *cParam = *fTPCInner;
2209   if (!cParam->Update(p,c)) return kFALSE;
2210
2211   return kTRUE;
2212 }
2213
2214 //_______________________________________________________________________
2215 Bool_t AliESDtrack::RelateToVertexTPCBxByBz(const AliESDVertex *vtx, 
2216 Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam) {
2217   //
2218   // Try to relate the TPC-only track parameters to the vertex "vtx", 
2219   // if the (rough) transverse impact parameter is not bigger then "maxd". 
2220   //
2221   // All three components of the magnetic field ,"b[3]" (kG), 
2222   // are taken into account.
2223   //
2224   // a) The TPC-only paramters are extapolated to the DCA to the vertex.
2225   // b) The impact parameters and their covariance matrix are calculated.
2226   // c) An attempt to constrain the TPC-only params to the vertex is done.
2227   //    The constrained params are returned via "cParam".
2228   //
2229   // In the case of success, the returned value is kTRUE
2230   // otherwise, it's kFALSE)
2231   // 
2232
2233   if (!fTPCInner) return kFALSE;
2234   if (!vtx) return kFALSE;
2235
2236   Double_t dz[2],cov[3];
2237   if (!fTPCInner->PropagateToDCABxByBz(vtx, b, maxd, dz, cov)) return kFALSE;
2238
2239   fdTPC = dz[0];
2240   fzTPC = dz[1];  
2241   fCddTPC = cov[0];
2242   fCdzTPC = cov[1];
2243   fCzzTPC = cov[2];
2244   
2245   Double_t covar[6]; vtx->GetCovMatrix(covar);
2246   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2247   Double_t c[3]={covar[2],0.,covar[5]};
2248
2249   Double_t chi2=GetPredictedChi2(p,c);
2250   if (chi2>kVeryBig) return kFALSE;
2251
2252   fCchi2TPC=chi2;
2253
2254   if (!cParam) return kTRUE;
2255
2256   *cParam = *fTPCInner;
2257   if (!cParam->Update(p,c)) return kFALSE;
2258
2259   return kTRUE;
2260 }
2261
2262 //_______________________________________________________________________
2263 Bool_t AliESDtrack::RelateToVertex(const AliESDVertex *vtx, 
2264 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
2265   //
2266   // Try to relate this track to the vertex "vtx", 
2267   // if the (rough) transverse impact parameter is not bigger then "maxd". 
2268   //            Magnetic field is "b" (kG).
2269   //
2270   // a) The track gets extapolated to the DCA to the vertex.
2271   // b) The impact parameters and their covariance matrix are calculated.
2272   // c) An attempt to constrain this track to the vertex is done.
2273   //    The constrained params are returned via "cParam".
2274   //
2275   // In the case of success, the returned value is kTRUE
2276   // (otherwise, it's kFALSE)
2277   //  
2278
2279   if (!vtx) return kFALSE;
2280
2281   Double_t dz[2],cov[3];
2282   if (!PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
2283
2284   fD = dz[0];
2285   fZ = dz[1];  
2286   fCdd = cov[0];
2287   fCdz = cov[1];
2288   fCzz = cov[2];
2289   
2290   Double_t covar[6]; vtx->GetCovMatrix(covar);
2291   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2292   Double_t c[3]={covar[2],0.,covar[5]};
2293
2294   Double_t chi2=GetPredictedChi2(p,c);
2295   if (chi2>kVeryBig) return kFALSE;
2296
2297   fCchi2=chi2;
2298
2299
2300   //--- Could now these lines be removed ? ---
2301   delete fCp;
2302   fCp=new AliExternalTrackParam(*this);  
2303
2304   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
2305   //----------------------------------------
2306
2307   fVertexID = vtx->GetID();
2308
2309   if (!cParam) return kTRUE;
2310
2311   *cParam = *this;
2312   if (!cParam->Update(p,c)) return kFALSE; 
2313
2314   return kTRUE;
2315 }
2316
2317 //_______________________________________________________________________
2318 Bool_t AliESDtrack::RelateToVertexBxByBz(const AliESDVertex *vtx, 
2319 Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam) {
2320   //
2321   // Try to relate this track to the vertex "vtx", 
2322   // if the (rough) transverse impact parameter is not bigger then "maxd". 
2323   //            Magnetic field is "b" (kG).
2324   //
2325   // a) The track gets extapolated to the DCA to the vertex.
2326   // b) The impact parameters and their covariance matrix are calculated.
2327   // c) An attempt to constrain this track to the vertex is done.
2328   //    The constrained params are returned via "cParam".
2329   //
2330   // In the case of success, the returned value is kTRUE
2331   // (otherwise, it's kFALSE)
2332   //  
2333
2334   if (!vtx) return kFALSE;
2335
2336   Double_t dz[2],cov[3];
2337   if (!PropagateToDCABxByBz(vtx, b, maxd, dz, cov)) return kFALSE;
2338
2339   fD = dz[0];
2340   fZ = dz[1];  
2341   fCdd = cov[0];
2342   fCdz = cov[1];
2343   fCzz = cov[2];
2344   
2345   Double_t covar[6]; vtx->GetCovMatrix(covar);
2346   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2347   Double_t c[3]={covar[2],0.,covar[5]};
2348
2349   Double_t chi2=GetPredictedChi2(p,c);
2350   if (chi2>kVeryBig) return kFALSE;
2351
2352   fCchi2=chi2;
2353
2354
2355   //--- Could now these lines be removed ? ---
2356   delete fCp;
2357   fCp=new AliExternalTrackParam(*this);  
2358
2359   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
2360   //----------------------------------------
2361
2362   fVertexID = vtx->GetID();
2363
2364   if (!cParam) return kTRUE;
2365
2366   *cParam = *this;
2367   if (!cParam->Update(p,c)) return kFALSE; 
2368
2369   return kTRUE;
2370 }
2371
2372 //_______________________________________________________________________
2373 void AliESDtrack::Print(Option_t *) const {
2374   // Prints info on the track
2375   AliExternalTrackParam::Print();
2376   printf("ESD track info\n") ; 
2377   Double_t p[AliPID::kSPECIESN] ; 
2378   Int_t index = 0 ; 
2379   if( IsOn(kITSpid) ){
2380     printf("From ITS: ") ; 
2381     GetITSpid(p) ; 
2382     for(index = 0 ; index < AliPID::kSPECIES; index++) 
2383       printf("%f, ", p[index]) ;
2384     printf("\n           signal = %f\n", GetITSsignal()) ;
2385   } 
2386   if( IsOn(kTPCpid) ){
2387     printf("From TPC: ") ; 
2388     GetTPCpid(p) ; 
2389     for(index = 0 ; index < AliPID::kSPECIES; index++) 
2390       printf("%f, ", p[index]) ;
2391     printf("\n           signal = %f\n", GetTPCsignal()) ;
2392   }
2393   if( IsOn(kTRDpid) ){
2394     printf("From TRD: ") ; 
2395     GetTRDpid(p) ; 
2396     for(index = 0 ; index < AliPID::kSPECIES; index++) 
2397       printf("%f, ", p[index]) ;
2398       printf("\n           signal = %f\n", GetTRDsignal()) ;
2399   }
2400   if( IsOn(kTOFpid) ){
2401     printf("From TOF: ") ; 
2402     GetTOFpid(p) ; 
2403     for(index = 0 ; index < AliPID::kSPECIES; index++) 
2404       printf("%f, ", p[index]) ;
2405     printf("\n           signal = %f\n", GetTOFsignal()) ;
2406   }
2407   if( IsOn(kHMPIDpid) ){
2408     printf("From HMPID: ") ; 
2409     GetHMPIDpid(p) ; 
2410     for(index = 0 ; index < AliPID::kSPECIES; index++) 
2411       printf("%f, ", p[index]) ;
2412     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
2413   }
2414
2415
2416
2417 //
2418 // Draw functionality
2419 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
2420 //
2421 void AliESDtrack::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
2422   //
2423   // Fill points in the polymarker
2424   //
2425   TObjArray arrayRef;
2426   arrayRef.AddLast(new AliExternalTrackParam(*this));
2427   if (fIp) arrayRef.AddLast(new AliExternalTrackParam(*fIp));
2428   if (fOp) arrayRef.AddLast(new AliExternalTrackParam(*fOp));
2429   if (fHMPIDp) arrayRef.AddLast(new AliExternalTrackParam(*fHMPIDp));
2430   //
2431   Double_t mpos[3]={0,0,0};
2432   Int_t entries=arrayRef.GetEntries();
2433   for (Int_t i=0;i<entries;i++){
2434     Double_t pos[3];
2435     ((AliExternalTrackParam*)arrayRef.At(i))->GetXYZ(pos);
2436     mpos[0]+=pos[0]/entries;
2437     mpos[1]+=pos[1]/entries;
2438     mpos[2]+=pos[2]/entries;    
2439   }
2440   // Rotate to the mean position
2441   //
2442   Float_t fi= TMath::ATan2(mpos[1],mpos[0]);
2443   for (Int_t i=0;i<entries;i++){
2444     Bool_t res = ((AliExternalTrackParam*)arrayRef.At(i))->Rotate(fi);
2445     if (!res) delete arrayRef.RemoveAt(i);
2446   }
2447   Int_t counter=0;
2448   for (Double_t r=minR; r<maxR; r+=stepR){
2449     Double_t sweight=0;
2450     Double_t mlpos[3]={0,0,0};
2451     for (Int_t i=0;i<entries;i++){
2452       Double_t point[3]={0,0,0};
2453       AliExternalTrackParam *param = ((AliExternalTrackParam*)arrayRef.At(i));
2454       if (!param) continue;
2455       if (param->GetXYZAt(r,magF,point)){
2456         Double_t weight = 1./(10.+(r-param->GetX())*(r-param->GetX()));
2457         sweight+=weight;
2458         mlpos[0]+=point[0]*weight;
2459         mlpos[1]+=point[1]*weight;
2460         mlpos[2]+=point[2]*weight;
2461       }
2462     }
2463     if (sweight>0){
2464       mlpos[0]/=sweight;
2465       mlpos[1]/=sweight;
2466       mlpos[2]/=sweight;      
2467       pol->SetPoint(counter,mlpos[0],mlpos[1], mlpos[2]);
2468       //      printf("xyz\t%f\t%f\t%f\n",mlpos[0], mlpos[1],mlpos[2]);
2469       counter++;
2470     }
2471   }
2472 }
2473
2474 //_______________________________________________________________________
2475 void AliESDtrack::SetITSdEdxSamples(const Double_t s[4]) {
2476   //
2477   // Store the dE/dx samples measured by the two SSD and two SDD layers.
2478   // These samples are corrected for the track segment length. 
2479   //
2480   for (Int_t i=0; i<4; i++) fITSdEdxSamples[i]=s[i];
2481 }
2482
2483 //_______________________________________________________________________
2484 void AliESDtrack::GetITSdEdxSamples(Double_t *s) const {
2485   //
2486   // Get the dE/dx samples measured by the two SSD and two SDD layers.  
2487   // These samples are corrected for the track segment length.
2488   //
2489   for (Int_t i=0; i<4; i++) s[i]=fITSdEdxSamples[i];
2490 }
2491
2492
2493 UShort_t   AliESDtrack::GetTPCnclsS(Int_t i0,Int_t i1) const{
2494   //
2495   // get number of shared TPC clusters
2496   //
2497   return  fTPCSharedMap.CountBits(i0)-fTPCSharedMap.CountBits(i1);
2498 }
2499
2500 UShort_t   AliESDtrack::GetTPCncls(Int_t i0,Int_t i1) const{
2501   //
2502   // get number of TPC clusters
2503   //
2504   return  fTPCClusterMap.CountBits(i0)-fTPCClusterMap.CountBits(i1);
2505 }
2506
2507 //____________________________________________________________________
2508 Double_t AliESDtrack::GetChi2TPCConstrainedVsGlobal(const AliESDVertex* vtx) const
2509 {
2510   // Calculates the chi2 between the TPC track (TPCinner) constrained to the primary vertex and the global track
2511   //
2512   // Returns -1 in case the calculation failed
2513   //
2514   // Value is cached as a non-persistent member.
2515   //
2516   // Code adapted from original code by GSI group (Jacek, Marian, Michael)
2517   
2518   // cache, ignoring that a different vertex might be passed
2519   if (fCacheChi2TPCConstrainedVsGlobalVertex == vtx)
2520     return fCacheChi2TPCConstrainedVsGlobal;
2521   
2522   fCacheChi2TPCConstrainedVsGlobal = -1;
2523   fCacheChi2TPCConstrainedVsGlobalVertex = vtx;
2524   
2525   Double_t x[3];
2526   GetXYZ(x);
2527   Double_t b[3];
2528   AliTrackerBase::GetBxByBz(x,b);
2529
2530   if (!fTPCInner)  { 
2531     AliWarning("Could not get TPC Inner Param.");
2532     return fCacheChi2TPCConstrainedVsGlobal;
2533   }
2534   
2535   // clone for constraining
2536   AliExternalTrackParam* tpcInnerC = new AliExternalTrackParam(*fTPCInner);
2537   if (!tpcInnerC) { 
2538     AliWarning("Clone of TPCInnerParam failed.");
2539     return fCacheChi2TPCConstrainedVsGlobal;  
2540   }
2541   
2542   // transform to the track reference frame 
2543   Bool_t isOK = tpcInnerC->Rotate(GetAlpha());
2544   isOK &= tpcInnerC->PropagateTo(GetX(), b[2]);
2545   if (!isOK) { 
2546     delete tpcInnerC;
2547     tpcInnerC = 0; 
2548     AliWarning("Rotation/Propagation of track failed.") ; 
2549     return fCacheChi2TPCConstrainedVsGlobal;    
2550   }  
2551
2552   // constrain TPCinner 
2553   isOK = tpcInnerC->ConstrainToVertex(vtx, b);
2554   
2555   // transform to the track reference frame 
2556   isOK &= tpcInnerC->Rotate(GetAlpha());
2557   isOK &= tpcInnerC->PropagateTo(GetX(), b[2]);
2558
2559   if (!isOK) {
2560     AliWarning("ConstrainTPCInner failed.") ;
2561     delete tpcInnerC;
2562     tpcInnerC = 0; 
2563     return fCacheChi2TPCConstrainedVsGlobal;  
2564   }
2565   
2566   // calculate chi2 between vi and vj vectors
2567   // with covi and covj covariance matrices
2568   // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)
2569   TMatrixD deltaT(5,1);
2570   TMatrixD delta(1,5);
2571   TMatrixD covarM(5,5);
2572
2573   for (Int_t ipar=0; ipar<5; ipar++) {
2574     deltaT(ipar,0) = tpcInnerC->GetParameter()[ipar] - GetParameter()[ipar];
2575     delta(0,ipar) = tpcInnerC->GetParameter()[ipar] - GetParameter()[ipar];
2576
2577     for (Int_t jpar=0; jpar<5; jpar++) {
2578       Int_t index = GetIndex(ipar,jpar);
2579       covarM(ipar,jpar) = GetCovariance()[index]+tpcInnerC->GetCovariance()[index];
2580     }
2581   }
2582   // chi2 distance TPC constrained and TPC+ITS
2583   TMatrixD covarMInv = covarM.Invert();
2584   TMatrixD mat2 = covarMInv*deltaT;
2585   TMatrixD chi2 = delta*mat2; 
2586   
2587   delete tpcInnerC; 
2588   tpcInnerC = 0;
2589   
2590   fCacheChi2TPCConstrainedVsGlobal = chi2(0,0);
2591   return fCacheChi2TPCConstrainedVsGlobal;
2592 }