Fix for #88344
[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   if (fFriendTrack) delete fFriendTrack;
743   if (fTPCdEdxInfo) delete fTPCdEdxInfo;
744   fFriendTrack=NULL;
745   if(fTRDnSlices)
746     delete[] fTRDslices;
747 }
748
749 AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
750   
751
752   if(&source == this) return *this;
753   AliExternalTrackParam::operator=(source);
754
755   
756   if(source.fCp){
757     // we have the trackparam: assign or copy construct
758     if(fCp)*fCp = *source.fCp;
759     else fCp = new AliExternalTrackParam(*source.fCp);
760   }
761   else{
762     // no track param delete the old one
763     if(fCp)delete fCp;
764     fCp = 0;
765   }
766
767   if(source.fIp){
768     // we have the trackparam: assign or copy construct
769     if(fIp)*fIp = *source.fIp;
770     else fIp = new AliExternalTrackParam(*source.fIp);
771   }
772   else{
773     // no track param delete the old one
774     if(fIp)delete fIp;
775     fIp = 0;
776   }
777
778
779   if(source.fTPCInner){
780     // we have the trackparam: assign or copy construct
781     if(fTPCInner) *fTPCInner = *source.fTPCInner;
782     else fTPCInner = new AliExternalTrackParam(*source.fTPCInner);
783   }
784   else{
785     // no track param delete the old one
786     if(fTPCInner)delete fTPCInner;
787     fTPCInner = 0;
788   }
789
790   if(source.fTPCdEdxInfo) {
791     if(fTPCdEdxInfo) *fTPCdEdxInfo = *source.fTPCdEdxInfo;
792     fTPCdEdxInfo = new AliTPCdEdxInfo(*source.fTPCdEdxInfo);
793   }
794
795   if(source.fOp){
796     // we have the trackparam: assign or copy construct
797     if(fOp) *fOp = *source.fOp;
798     else fOp = new AliExternalTrackParam(*source.fOp);
799   }
800   else{
801     // no track param delete the old one
802     if(fOp)delete fOp;
803     fOp = 0;
804   }
805
806   
807   if(source.fHMPIDp){
808     // we have the trackparam: assign or copy construct
809     if(fHMPIDp) *fHMPIDp = *source.fHMPIDp;
810     else fHMPIDp = new AliExternalTrackParam(*source.fHMPIDp);
811   }
812   else{
813     // no track param delete the old one
814     if(fHMPIDp)delete fHMPIDp;
815     fHMPIDp = 0;
816   }
817
818   // copy also the friend track 
819   // use copy constructor
820   if(source.fFriendTrack){
821     // we have the trackparam: assign or copy construct
822     delete fFriendTrack; fFriendTrack=new AliESDfriendTrack(*source.fFriendTrack);
823   }
824   else{
825     // no track param delete the old one
826     delete fFriendTrack; fFriendTrack= 0;
827   }
828
829   fTPCFitMap = source.fTPCFitMap; 
830   fTPCClusterMap = source.fTPCClusterMap; 
831   fTPCSharedMap  = source.fTPCSharedMap;  
832   // the simple stuff
833   fFlags    = source.fFlags; 
834   fID       = source.fID;             
835   fLabel    = source.fLabel;
836   fITSLabel = source.fITSLabel;
837   for(int i = 0; i< 12;++i){
838     fITSModule[i] = source.fITSModule[i];
839   }
840   fTPCLabel = source.fTPCLabel; 
841   fTRDLabel = source.fTRDLabel;
842   for(int i = 0; i< 3;++i){
843     fTOFLabel[i] = source.fTOFLabel[i];    
844   }
845   fTOFCalChannel = source.fTOFCalChannel;
846   fTOFindex      = source.fTOFindex;
847   fHMPIDqn       = source.fHMPIDqn;
848   fHMPIDcluIdx   = source.fHMPIDcluIdx; 
849   fCaloIndex    = source.fCaloIndex;
850   for (int i=kNITSchi2Std;i--;) fITSchi2Std[i] = source.fITSchi2Std[i];
851   for(int i = 0; i< 3;++i){
852     fKinkIndexes[i] = source.fKinkIndexes[i]; 
853     fV0Indexes[i]   = source.fV0Indexes[i]; 
854   }
855
856   for(int i = 0; i< AliPID::kSPECIES;++i){
857     fR[i]     = source.fR[i];
858     fITSr[i]  = source.fITSr[i];
859     fTPCr[i]  = source.fTPCr[i];
860     fTRDr[i]  = source.fTRDr[i];
861     fTOFr[i]  = source.fTOFr[i];
862     fHMPIDr[i] = source.fHMPIDr[i];
863     fTrackTime[i] = source.fTrackTime[i];  
864   }
865
866   fHMPIDtrkTheta = source.fHMPIDtrkTheta;
867   fHMPIDtrkPhi   = source.fHMPIDtrkPhi;
868   fHMPIDsignal   = source.fHMPIDsignal; 
869
870   
871   fTrackLength   = source. fTrackLength;
872   fdTPC  = source.fdTPC; 
873   fzTPC  = source.fzTPC; 
874   fCddTPC = source.fCddTPC;
875   fCdzTPC = source.fCdzTPC;
876   fCzzTPC = source.fCzzTPC;
877   fCchi2TPC = source.fCchi2TPC;
878
879   fD  = source.fD; 
880   fZ  = source.fZ; 
881   fCdd = source.fCdd;
882   fCdz = source.fCdz;
883   fCzz = source.fCzz;
884   fCchi2     = source.fCchi2;
885
886   fITSchi2   = source.fITSchi2;             
887   fTPCchi2   = source.fTPCchi2;            
888   fTPCchi2Iter1   = source.fTPCchi2Iter1;            
889   fTRDchi2   = source.fTRDchi2;      
890   fTOFchi2   = source.fTOFchi2;      
891   fHMPIDchi2 = source.fHMPIDchi2;      
892
893   fGlobalChi2 = source.fGlobalChi2;      
894
895   fITSsignal  = source.fITSsignal;     
896   for (Int_t i=0;i<4;i++) {fITSdEdxSamples[i]=source.fITSdEdxSamples[i];}
897   fTPCsignal  = source.fTPCsignal;     
898   fTPCsignalS = source.fTPCsignalS;    
899   for(int i = 0; i< 4;++i){
900     fTPCPoints[i] = source.fTPCPoints[i];  
901   }
902   fTRDsignal = source.fTRDsignal;
903
904   for(int i = 0;i < kTRDnPlanes;++i){
905     fTRDTimBin[i] = source.fTRDTimBin[i];   
906   }
907
908   if(fTRDnSlices)
909     delete[] fTRDslices;
910   fTRDslices=0;
911   fTRDnSlices=source.fTRDnSlices;
912   if (fTRDnSlices) {
913     fTRDslices=new Double32_t[fTRDnSlices];
914     for(int j = 0;j < fTRDnSlices;++j) fTRDslices[j] = source.fTRDslices[j];
915   }
916
917   fTRDQuality =   source.fTRDQuality;     
918   fTRDBudget  =   source.fTRDBudget;      
919   fTOFsignal  =   source.fTOFsignal;     
920   fTOFsignalToT = source.fTOFsignalToT;   
921   fTOFsignalRaw = source.fTOFsignalRaw;  
922   fTOFsignalDz  = source.fTOFsignalDz;      
923   fTOFsignalDx  = source.fTOFsignalDx;      
924   fTOFdeltaBC   = source.fTOFdeltaBC;
925   fTOFl0l1      = source.fTOFl0l1;
926  
927   for(int i = 0;i<10;++i){
928     fTOFInfo[i] = source.fTOFInfo[i];    
929   }
930
931   fHMPIDtrkX = source.fHMPIDtrkX; 
932   fHMPIDtrkY = source.fHMPIDtrkY; 
933   fHMPIDmipX = source.fHMPIDmipX;
934   fHMPIDmipY = source.fHMPIDmipY; 
935
936   fTPCncls    = source.fTPCncls;      
937   fTPCnclsF   = source.fTPCnclsF;     
938   fTPCsignalN = source.fTPCsignalN;   
939   fTPCnclsIter1    = source.fTPCnclsIter1;      
940   fTPCnclsFIter1   = source.fTPCnclsFIter1;     
941
942   fITSncls = source.fITSncls;       
943   fITSClusterMap = source.fITSClusterMap; 
944   fITSSharedMap = source.fITSSharedMap; 
945   fTRDncls   = source.fTRDncls;       
946   fTRDncls0  = source.fTRDncls0;      
947   fTRDntracklets  = source.fTRDntracklets; 
948   fVertexID = source.fVertexID;
949   return *this;
950 }
951
952
953
954 void AliESDtrack::Copy(TObject &obj) const {
955   
956   // this overwrites the virtual TOBject::Copy()
957   // to allow run time copying without casting
958   // in AliESDEvent
959
960   if(this==&obj)return;
961   AliESDtrack *robj = dynamic_cast<AliESDtrack*>(&obj);
962   if(!robj)return; // not an AliESDtrack
963   *robj = *this;
964
965 }
966
967
968
969 void AliESDtrack::AddCalibObject(TObject * object){
970   //
971   // add calib object to the list
972   //
973   if (!fFriendTrack) fFriendTrack  = new AliESDfriendTrack;
974   if (!fFriendTrack) return;
975   fFriendTrack->AddCalibObject(object);
976 }
977
978 TObject *  AliESDtrack::GetCalibObject(Int_t index){
979   //
980   // return calib objct at given position
981   //
982   if (!fFriendTrack) return 0;
983   return fFriendTrack->GetCalibObject(index);
984 }
985
986
987 Bool_t AliESDtrack::FillTPCOnlyTrack(AliESDtrack &track){
988   
989   // Fills the information of the TPC-only first reconstruction pass
990   // into the passed ESDtrack object. For consistency fTPCInner is also filled
991   // again
992
993
994
995   // For data produced before r26675
996   // RelateToVertexTPC was not properly called during reco
997   // so you'll have to call it again, before FillTPCOnlyTrack
998   //  Float_t p[2],cov[3];
999   // track->GetImpactParametersTPC(p,cov); 
1000   // if(p[0]==0&&p[1]==0) // <- Default values
1001   //  track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
1002   
1003
1004   if(!fTPCInner)return kFALSE;
1005
1006   // fill the TPC track params to the global track parameters
1007   track.Set(fTPCInner->GetX(),fTPCInner->GetAlpha(),fTPCInner->GetParameter(),fTPCInner->GetCovariance());
1008   track.fD = fdTPC;
1009   track.fZ = fzTPC;
1010   track.fCdd = fCddTPC;
1011   track.fCdz = fCdzTPC;
1012   track.fCzz = fCzzTPC;
1013
1014   // copy the inner params
1015   if(track.fIp) *track.fIp = *fIp;
1016   else track.fIp = new AliExternalTrackParam(*fIp);
1017
1018   // copy the TPCinner parameters
1019   if(track.fTPCInner) *track.fTPCInner = *fTPCInner;
1020   else track.fTPCInner = new AliExternalTrackParam(*fTPCInner);
1021   track.fdTPC   = fdTPC;
1022   track.fzTPC   = fzTPC;
1023   track.fCddTPC = fCddTPC;
1024   track.fCdzTPC = fCdzTPC;
1025   track.fCzzTPC = fCzzTPC;
1026   track.fCchi2TPC = fCchi2TPC;
1027
1028   // copy all other TPC specific parameters
1029
1030   // replace label by TPC label
1031   track.fLabel    = fTPCLabel;
1032   track.fTPCLabel = fTPCLabel;
1033
1034   track.fTPCchi2 = fTPCchi2; 
1035   track.fTPCchi2Iter1 = fTPCchi2Iter1; 
1036   track.fTPCsignal = fTPCsignal;
1037   track.fTPCsignalS = fTPCsignalS;
1038   for(int i = 0;i<4;++i)track.fTPCPoints[i] = fTPCPoints[i];
1039
1040   track.fTPCncls    = fTPCncls;     
1041   track.fTPCnclsF   = fTPCnclsF;     
1042   track.fTPCsignalN =  fTPCsignalN;
1043   track.fTPCnclsIter1    = fTPCnclsIter1;     
1044   track.fTPCnclsFIter1   = fTPCnclsFIter1;     
1045
1046   // PID 
1047   for(int i=0;i<AliPID::kSPECIES;++i){
1048     track.fTPCr[i] = fTPCr[i];
1049     // combined PID is TPC only!
1050     track.fR[i] = fTPCr[i];
1051   }
1052   track.fTPCFitMap = fTPCFitMap;
1053   track.fTPCClusterMap = fTPCClusterMap;
1054   track.fTPCSharedMap = fTPCSharedMap;
1055
1056
1057   // reset the flags
1058   track.fFlags = kTPCin;
1059   track.fID    = fID;
1060
1061   track.fFlags |= fFlags & kTPCpid; //copy the TPCpid status flag
1062  
1063   for (Int_t i=0;i<3;i++) track.fKinkIndexes[i] = fKinkIndexes[i];
1064   
1065   return kTRUE;
1066     
1067 }
1068
1069 //_______________________________________________________________________
1070 void AliESDtrack::MakeMiniESDtrack(){
1071   // Resets everything except
1072   // fFlags: Reconstruction status flags 
1073   // fLabel: Track label
1074   // fID:  Unique ID of the track
1075   // Impact parameter information
1076   // fR[AliPID::kSPECIES]: combined "detector response probability"
1077   // Running track parameters in the base class (AliExternalTrackParam)
1078   
1079   fTrackLength = 0;
1080
1081   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i] = 0;
1082
1083   // Reset track parameters constrained to the primary vertex
1084   delete fCp;fCp = 0;
1085
1086   // Reset track parameters at the inner wall of TPC
1087   delete fIp;fIp = 0;
1088   delete fTPCInner;fTPCInner=0;
1089   // Reset track parameters at the inner wall of the TRD
1090   delete fOp;fOp = 0;
1091   // Reset track parameters at the HMPID
1092   delete fHMPIDp;fHMPIDp = 0;
1093
1094
1095   // Reset ITS track related information
1096   fITSchi2 = 0;
1097   fITSncls = 0;       
1098   fITSClusterMap=0;
1099   fITSSharedMap=0;
1100   fITSsignal = 0;     
1101   for (Int_t i=0;i<4;i++) fITSdEdxSamples[i] = 0.;
1102   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=0; 
1103   fITSLabel = 0;       
1104
1105   // Reset TPC related track information
1106   fTPCchi2 = 0;       
1107   fTPCchi2Iter1 = 0;       
1108   fTPCncls = 0;       
1109   fTPCnclsF = 0;       
1110   fTPCnclsIter1 = 0;       
1111   fTPCnclsFIter1 = 0;  
1112   fTPCFitMap = 0;       
1113   fTPCClusterMap = 0;  
1114   fTPCSharedMap = 0;  
1115   fTPCsignal= 0;      
1116   fTPCsignalS= 0;      
1117   fTPCsignalN= 0;      
1118   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0; 
1119   fTPCLabel=0;       
1120   for (Int_t i=0;i<4;i++) fTPCPoints[i] = 0;
1121   for (Int_t i=0; i<3;i++)   fKinkIndexes[i] = 0;
1122   for (Int_t i=0; i<3;i++)   fV0Indexes[i] = 0;
1123
1124   // Reset TRD related track information
1125   fTRDchi2 = 0;        
1126   fTRDncls = 0;       
1127   fTRDncls0 = 0;       
1128   fTRDsignal = 0;      
1129   for (Int_t i=0;i<kTRDnPlanes;i++) {
1130     fTRDTimBin[i]  = 0;
1131   }
1132   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0; 
1133   fTRDLabel = 0;       
1134   fTRDQuality  = 0;
1135   fTRDntracklets = 0;
1136   if(fTRDnSlices)
1137     delete[] fTRDslices;
1138   fTRDslices=0x0;
1139   fTRDnSlices=0;
1140   fTRDBudget  = 0;
1141
1142   // Reset TOF related track information
1143   fTOFchi2 = 0;        
1144   fTOFindex = -1;       
1145   fTOFsignal = 99999;      
1146   fTOFCalChannel = -1;
1147   fTOFsignalToT = 99999;
1148   fTOFsignalRaw = 99999;
1149   fTOFsignalDz = 999;
1150   fTOFsignalDx = 999;
1151   fTOFdeltaBC = 999;
1152   fTOFl0l1 = 999;
1153   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i] = 0;
1154   for (Int_t i=0;i<3;i++) fTOFLabel[i] = -1;
1155   for (Int_t i=0;i<10;i++) fTOFInfo[i] = 0;
1156
1157   // Reset HMPID related track information
1158   fHMPIDchi2 = 0;     
1159   fHMPIDqn = 0;     
1160   fHMPIDcluIdx = -1;     
1161   fHMPIDsignal = 0;     
1162   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i] = 0;
1163   fHMPIDtrkTheta = 0;     
1164   fHMPIDtrkPhi = 0;      
1165   fHMPIDtrkX = 0;     
1166   fHMPIDtrkY = 0;      
1167   fHMPIDmipX = 0;
1168   fHMPIDmipY = 0;
1169   fCaloIndex = kEMCALNoMatch;
1170
1171   // reset global track chi2
1172   fGlobalChi2 = 0;
1173
1174   fVertexID = -2; // an orphan track
1175
1176   delete fFriendTrack; fFriendTrack = 0;
1177
1178
1179 //_______________________________________________________________________
1180 Int_t AliESDtrack::GetPID() const 
1181 {
1182   // Returns the particle most probable id
1183   Int_t i;
1184   for (i=0; i<AliPID::kSPECIES-1; i++) if (fR[i] != fR[i+1]) break;
1185   //
1186   if (i == AliPID::kSPECIES-1) return AliPID::kPion;  // If all the probabilities are equal, return the pion mass
1187   //
1188   Float_t max=0.;
1189   Int_t k=-1;
1190   for (i=0; i<AliPID::kSPECIES; i++) if (fR[i]>max) {k=i; max=fR[i];}
1191   //
1192   if (k==0) { // dE/dx "crossing points" in the TPC
1193     Double_t p=GetP();
1194     if ((p>0.38)&&(p<0.48))
1195       if (fR[0]<fR[3]*10.) return AliPID::kKaon;
1196     if ((p>0.75)&&(p<0.85))
1197       if (fR[0]<fR[4]*10.) return AliPID::kProton;
1198     return AliPID::kElectron;
1199   }
1200   if (k==1) return AliPID::kMuon; 
1201   if (k==2||k==-1) return AliPID::kPion;
1202   if (k==3) return AliPID::kKaon;
1203   if (k==4) return AliPID::kProton;
1204   AliWarning("Undefined PID !");
1205   return AliPID::kPion;
1206 }
1207
1208 //_______________________________________________________________________
1209 Int_t AliESDtrack::GetTOFBunchCrossing(Double_t b) const 
1210 {
1211   // Returns the number of bunch crossings after trigger (assuming 25ns spacing)
1212   const double kSpacing = 25e3; // min interbanch spacing
1213   const double kShift = 0;
1214   Int_t bcid = kTOFBCNA; // defualt one
1215   if (!IsOn(kTOFout) || !IsOn(kESDpid)) return bcid; // no info
1216   //
1217   double tdif = fTOFsignal;
1218   if (IsOn(kTIME)) { // integrated time info is there
1219     int pid = GetPID();
1220     tdif -= fTrackTime[pid];
1221   }
1222   else { // assume integrated time info from TOF radius and momentum
1223     const double kRTOF = 385.;
1224     const double kCSpeed = 3.e-2; // cm/ps
1225     double p = GetP();
1226     if (p<0.01) return bcid;
1227     double m = GetMass();
1228     double curv = GetC(b);
1229     double path = TMath::Abs(curv)>kAlmost0 ? // account for curvature
1230       2./curv*TMath::ASin(kRTOF*curv/2.)*TMath::Sqrt(1.+GetTgl()*GetTgl()) : kRTOF;
1231     tdif -= path/kCSpeed*TMath::Sqrt(1.+m*m/(p*p));
1232   }
1233   bcid = TMath::Nint((tdif - kShift)/kSpacing);
1234   return bcid;
1235 }
1236
1237 //______________________________________________________________________________
1238 Double_t AliESDtrack::M() const
1239 {
1240   // Returns the assumed mass
1241   // (the pion mass, if the particle can't be identified properly).
1242   static Bool_t printerr=kTRUE;
1243   if (printerr) {
1244      AliWarning("WARNING !!! ... THIS WILL BE PRINTED JUST ONCE !!!");
1245      printerr = kFALSE;
1246      AliWarning("This is the ESD mass. Use it with care !"); 
1247   }
1248   return GetMass(); 
1249 }
1250   
1251 //______________________________________________________________________________
1252 Double_t AliESDtrack::E() const
1253 {
1254   // Returns the energy of the particle given its assumed mass.
1255   // Assumes the pion mass if the particle can't be identified properly.
1256   
1257   Double_t m = M();
1258   Double_t p = P();
1259   return TMath::Sqrt(p*p + m*m);
1260 }
1261
1262 //______________________________________________________________________________
1263 Double_t AliESDtrack::Y() const
1264 {
1265   // Returns the rapidity of a particle given its assumed mass.
1266   // Assumes the pion mass if the particle can't be identified properly.
1267   
1268   Double_t e = E();
1269   Double_t pz = Pz();
1270   if (e != TMath::Abs(pz)) { // energy was not equal to pz
1271     return 0.5*TMath::Log((e+pz)/(e-pz));
1272   } else { // energy was equal to pz
1273     return -999.;
1274   }
1275 }
1276
1277 //_______________________________________________________________________
1278 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
1279   //
1280   // This function updates track's running parameters 
1281   //
1282   Bool_t rc=kTRUE;
1283
1284   SetStatus(flags);
1285   fLabel=t->GetLabel();
1286
1287   if (t->IsStartedTimeIntegral()) {
1288     SetStatus(kTIME);
1289     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
1290     SetIntegratedLength(t->GetIntegratedLength());
1291   }
1292
1293   Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1294   if (fFriendTrack) {
1295   if (flags==kITSout) fFriendTrack->SetITSOut(*t);
1296   if (flags==kTPCout) fFriendTrack->SetTPCOut(*t);
1297   if (flags==kTRDrefit) fFriendTrack->SetTRDIn(*t);
1298   }
1299   
1300   switch (flags) {
1301     
1302   case kITSin: 
1303     fITSchi2Std[0] = t->GetChi2();
1304     //
1305   case kITSout: 
1306     fITSchi2Std[1] = t->GetChi2();
1307   case kITSrefit:
1308     {
1309     fITSchi2Std[2] = t->GetChi2();
1310     fITSClusterMap=0;
1311     fITSncls=t->GetNumberOfClusters();
1312     if (fFriendTrack) {
1313     Int_t* indexITS = new Int_t[AliESDfriendTrack::kMaxITScluster];
1314     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
1315         indexITS[i]=t->GetClusterIndex(i);
1316
1317         if (i<fITSncls) {
1318           Int_t l=(indexITS[i] & 0xf0000000) >> 28;
1319            SETBIT(fITSClusterMap,l);                 
1320         }
1321     }
1322     fFriendTrack->SetITSIndices(indexITS,AliESDfriendTrack::kMaxITScluster);
1323     delete [] indexITS;
1324     }
1325
1326     fITSchi2=t->GetChi2();
1327     fITSsignal=t->GetPIDsignal();
1328     fITSLabel = t->GetLabel();
1329     // keep in fOp the parameters outside ITS for ITS stand-alone tracks 
1330     if (flags==kITSout) { 
1331       if (!fOp) fOp=new AliExternalTrackParam(*t);
1332       else 
1333         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1334     }   
1335     }
1336     break;
1337     
1338   case kTPCin: case kTPCrefit:
1339     {
1340     fTPCLabel = t->GetLabel();
1341     if (flags==kTPCin)  {
1342       fTPCInner=new AliExternalTrackParam(*t); 
1343       fTPCnclsIter1=t->GetNumberOfClusters();    
1344       fTPCchi2Iter1=t->GetChi2();
1345     }
1346     if (!fIp) fIp=new AliExternalTrackParam(*t);
1347     else 
1348       fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1349     }
1350     // Intentionally no break statement; need to set general TPC variables as well
1351   case kTPCout:
1352     {
1353     if (flags & kTPCout){
1354       if (!fOp) fOp=new AliExternalTrackParam(*t);
1355       else 
1356         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1357     }
1358     fTPCncls=t->GetNumberOfClusters();    
1359     fTPCchi2=t->GetChi2();
1360     
1361     if (fFriendTrack) {  // Copy cluster indices
1362       Int_t* indexTPC = new Int_t[AliESDfriendTrack::kMaxTPCcluster];
1363       for (Int_t i=0;i<AliESDfriendTrack::kMaxTPCcluster;i++)         
1364         indexTPC[i]=t->GetClusterIndex(i);
1365       fFriendTrack->SetTPCIndices(indexTPC,AliESDfriendTrack::kMaxTPCcluster);
1366       delete [] indexTPC;
1367     }
1368     fTPCsignal=t->GetPIDsignal();
1369     }
1370     break;
1371
1372   case kTRDin: case kTRDrefit:
1373     break;
1374   case kTRDout:
1375     {
1376     fTRDLabel = t->GetLabel(); 
1377     fTRDchi2  = t->GetChi2();
1378     fTRDncls  = t->GetNumberOfClusters();
1379     if (fFriendTrack) {
1380       Int_t* indexTRD = new Int_t[AliESDfriendTrack::kMaxTRDcluster];
1381       for (Int_t i=0;i<AliESDfriendTrack::kMaxTRDcluster;i++) indexTRD[i]=-2;
1382       for (Int_t i=0;i<6;i++) indexTRD[i]=t->GetTrackletIndex(i);
1383       fFriendTrack->SetTRDIndices(indexTRD,AliESDfriendTrack::kMaxTRDcluster);
1384       delete [] indexTRD;
1385     }    
1386     
1387     fTRDsignal=t->GetPIDsignal();
1388     }
1389     break;
1390   case kTRDbackup:
1391     if (!fOp) fOp=new AliExternalTrackParam(*t);
1392     else 
1393       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1394     fTRDncls0 = t->GetNumberOfClusters(); 
1395     break;
1396   case kTOFin: 
1397     break;
1398   case kTOFout: 
1399     break;
1400   case kTRDStop:
1401     break;
1402   case kHMPIDout:
1403   if (!fHMPIDp) fHMPIDp=new AliExternalTrackParam(*t);
1404     else 
1405       fHMPIDp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1406     break;
1407   default: 
1408     AliError("Wrong flag !");
1409     return kFALSE;
1410   }
1411
1412   return rc;
1413 }
1414
1415 //_______________________________________________________________________
1416 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
1417   //---------------------------------------------------------------------
1418   // This function returns external representation of the track parameters
1419   //---------------------------------------------------------------------
1420   x=GetX();
1421   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
1422 }
1423
1424 //_______________________________________________________________________
1425 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
1426   //---------------------------------------------------------------------
1427   // This function returns external representation of the cov. matrix
1428   //---------------------------------------------------------------------
1429   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
1430 }
1431
1432 //_______________________________________________________________________
1433 Bool_t AliESDtrack::GetConstrainedExternalParameters
1434                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1435   //---------------------------------------------------------------------
1436   // This function returns the constrained external track parameters
1437   //---------------------------------------------------------------------
1438   if (!fCp) return kFALSE;
1439   alpha=fCp->GetAlpha();
1440   x=fCp->GetX();
1441   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
1442   return kTRUE;
1443 }
1444
1445 //_______________________________________________________________________
1446 Bool_t 
1447 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
1448   //---------------------------------------------------------------------
1449   // This function returns the constrained external cov. matrix
1450   //---------------------------------------------------------------------
1451   if (!fCp) return kFALSE;
1452   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
1453   return kTRUE;
1454 }
1455
1456 Bool_t
1457 AliESDtrack::GetInnerExternalParameters
1458                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1459   //---------------------------------------------------------------------
1460   // This function returns external representation of the track parameters 
1461   // at the inner layer of TPC
1462   //---------------------------------------------------------------------
1463   if (!fIp) return kFALSE;
1464   alpha=fIp->GetAlpha();
1465   x=fIp->GetX();
1466   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
1467   return kTRUE;
1468 }
1469
1470 Bool_t 
1471 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
1472  //---------------------------------------------------------------------
1473  // This function returns external representation of the cov. matrix 
1474  // at the inner layer of TPC
1475  //---------------------------------------------------------------------
1476   if (!fIp) return kFALSE;
1477   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
1478   return kTRUE;
1479 }
1480
1481 void 
1482 AliESDtrack::SetOuterParam(const AliExternalTrackParam *p, ULong_t flags) {
1483   //
1484   // This is a direct setter for the outer track parameters
1485   //
1486   SetStatus(flags);
1487   if (fOp) delete fOp;
1488   fOp=new AliExternalTrackParam(*p);
1489 }
1490
1491 void 
1492 AliESDtrack::SetOuterHmpParam(const AliExternalTrackParam *p, ULong_t flags) {
1493   //
1494   // This is a direct setter for the outer track parameters
1495   //
1496   SetStatus(flags);
1497   if (fHMPIDp) delete fHMPIDp;
1498   fHMPIDp=new AliExternalTrackParam(*p);
1499 }
1500
1501 Bool_t 
1502 AliESDtrack::GetOuterExternalParameters
1503                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1504   //---------------------------------------------------------------------
1505   // This function returns external representation of the track parameters 
1506   // at the inner layer of TRD
1507   //---------------------------------------------------------------------
1508   if (!fOp) return kFALSE;
1509   alpha=fOp->GetAlpha();
1510   x=fOp->GetX();
1511   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
1512   return kTRUE;
1513 }
1514
1515 Bool_t 
1516 AliESDtrack::GetOuterHmpExternalParameters
1517                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1518   //---------------------------------------------------------------------
1519   // This function returns external representation of the track parameters 
1520   // at the inner layer of TRD
1521   //---------------------------------------------------------------------
1522   if (!fHMPIDp) return kFALSE;
1523   alpha=fHMPIDp->GetAlpha();
1524   x=fHMPIDp->GetX();
1525   for (Int_t i=0; i<5; i++) p[i]=fHMPIDp->GetParameter()[i];
1526   return kTRUE;
1527 }
1528
1529 Bool_t 
1530 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
1531  //---------------------------------------------------------------------
1532  // This function returns external representation of the cov. matrix 
1533  // at the inner layer of TRD
1534  //---------------------------------------------------------------------
1535   if (!fOp) return kFALSE;
1536   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
1537   return kTRUE;
1538 }
1539
1540 Bool_t 
1541 AliESDtrack::GetOuterHmpExternalCovariance(Double_t cov[15]) const {
1542  //---------------------------------------------------------------------
1543  // This function returns external representation of the cov. matrix 
1544  // at the inner layer of TRD
1545  //---------------------------------------------------------------------
1546   if (!fHMPIDp) return kFALSE;
1547   for (Int_t i=0; i<15; i++) cov[i]=fHMPIDp->GetCovariance()[i];
1548   return kTRUE;
1549 }
1550
1551 Int_t AliESDtrack::GetNcls(Int_t idet) const
1552 {
1553   // Get number of clusters by subdetector index
1554   //
1555   Int_t ncls = 0;
1556   switch(idet){
1557   case 0:
1558     ncls = fITSncls;
1559     break;
1560   case 1:
1561     ncls = fTPCncls;
1562     break;
1563   case 2:
1564     ncls = fTRDncls;
1565     break;
1566   case 3:
1567     if (fTOFindex != -1)
1568       ncls = 1;
1569     break;
1570   case 4: //PHOS
1571     break;
1572   case 5: //HMPID
1573     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1574       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1575         ncls = 1;
1576       }
1577     }    
1578     break;
1579   default:
1580     break;
1581   }
1582   return ncls;
1583 }
1584
1585 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
1586 {
1587   // Get cluster index array by subdetector index
1588   //
1589   Int_t ncls = 0;
1590   switch(idet){
1591   case 0:
1592     ncls = GetITSclusters(idx);
1593     break;
1594   case 1:
1595     ncls = GetTPCclusters(idx);
1596     break;
1597   case 2:
1598     ncls = GetTRDclusters(idx);
1599     break;
1600   case 3:
1601     if (fTOFindex != -1) {
1602       idx[0] = fTOFindex;
1603       ncls = 1;
1604     }
1605     break;
1606   case 4: //PHOS
1607     break;
1608   case 5:
1609     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1610       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1611         idx[0] = GetHMPIDcluIdx();
1612         ncls = 1;
1613       }
1614     }    
1615     break;
1616   case 6: //EMCAL
1617     break;
1618   default:
1619     break;
1620   }
1621   return ncls;
1622 }
1623
1624 //_______________________________________________________________________
1625 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
1626   // Returns the array with integrated times for each particle hypothesis
1627   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
1628 }
1629
1630 //_______________________________________________________________________
1631 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
1632   // Sets the array with integrated times for each particle hypotesis
1633   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
1634 }
1635
1636 //_______________________________________________________________________
1637 void AliESDtrack::SetITSpid(const Double_t *p) {
1638   // Sets values for the probability of each particle type (in ITS)
1639   SetPIDValues(fITSr,p,AliPID::kSPECIES);
1640   SetStatus(AliESDtrack::kITSpid);
1641 }
1642
1643 //_______________________________________________________________________
1644 void AliESDtrack::GetITSpid(Double_t *p) const {
1645   // Gets the probability of each particle type (in ITS)
1646   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
1647 }
1648
1649 //_______________________________________________________________________
1650 Char_t AliESDtrack::GetITSclusters(Int_t *idx) const {
1651   //---------------------------------------------------------------------
1652   // This function returns indices of the assgined ITS clusters 
1653   //---------------------------------------------------------------------
1654   if (idx && fFriendTrack) {
1655     Int_t *index=fFriendTrack->GetITSindices();
1656     for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) {
1657       if ( (i>=fITSncls) && (i<6) ) idx[i]=-1;
1658       else {
1659         if (index) {
1660           idx[i]=index[i];
1661         }
1662         else idx[i]= -2;
1663       }
1664     }
1665   }
1666   return fITSncls;
1667 }
1668
1669 //_______________________________________________________________________
1670 Bool_t AliESDtrack::GetITSModuleIndexInfo(Int_t ilayer,Int_t &idet,Int_t &status,
1671                                          Float_t &xloc,Float_t &zloc) const {
1672   //----------------------------------------------------------------------
1673   // This function encodes in the module number also the status of cluster association
1674   // "status" can have the following values: 
1675   // 1 "found" (cluster is associated), 
1676   // 2 "dead" (module is dead from OCDB), 
1677   // 3 "skipped" (module or layer forced to be skipped),
1678   // 4 "outinz" (track out of z acceptance), 
1679   // 5 "nocls" (no clusters in the road), 
1680   // 6 "norefit" (cluster rejected during refit), 
1681   // 7 "deadzspd" (holes in z in SPD)
1682   // Also given are the coordinates of the crossing point of track and module
1683   // (in the local module ref. system)
1684   // WARNING: THIS METHOD HAS TO BE SYNCHRONIZED WITH AliITStrackV2::GetModuleIndexInfo()!
1685   //----------------------------------------------------------------------
1686
1687   if(fITSModule[ilayer]==-1) {
1688     idet = -1;
1689     status=0;
1690     xloc=-99.; zloc=-99.;
1691     return kFALSE;
1692   }
1693
1694   Int_t module = fITSModule[ilayer];
1695
1696   idet = Int_t(module/1000000);
1697
1698   module -= idet*1000000;
1699
1700   status = Int_t(module/100000);
1701
1702   module -= status*100000;
1703
1704   Int_t signs = Int_t(module/10000);
1705
1706   module-=signs*10000;
1707
1708   Int_t xInt = Int_t(module/100);
1709   module -= xInt*100;
1710
1711   Int_t zInt = module;
1712
1713   if(signs==1) { xInt*=1; zInt*=1; }
1714   if(signs==2) { xInt*=1; zInt*=-1; }
1715   if(signs==3) { xInt*=-1; zInt*=1; }
1716   if(signs==4) { xInt*=-1; zInt*=-1; }
1717
1718   xloc = 0.1*(Float_t)xInt;
1719   zloc = 0.1*(Float_t)zInt;
1720
1721   if(status==4) idet = -1;
1722
1723   return kTRUE;
1724 }
1725
1726 //_______________________________________________________________________
1727 UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
1728   //---------------------------------------------------------------------
1729   // This function returns indices of the assgined ITS clusters 
1730   //---------------------------------------------------------------------
1731   if (idx && fFriendTrack) {
1732     Int_t *index=fFriendTrack->GetTPCindices();
1733
1734     if (index){
1735       for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
1736     }
1737     else {
1738       for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=-2;
1739     }
1740   }
1741   return fTPCncls;
1742 }
1743
1744 //_______________________________________________________________________
1745 Float_t AliESDtrack::GetTPCCrossedRows() const
1746 {
1747   // This function calls GetTPCClusterInfo with some default parameters which are used in the track selection and caches the outcome
1748   // because GetTPCClusterInfo is quite time-consuming
1749   
1750   if (fCacheNCrossedRows > -1)
1751     return fCacheNCrossedRows;
1752   
1753   fCacheNCrossedRows = GetTPCClusterInfo(2, 1);
1754   return fCacheNCrossedRows;
1755 }
1756
1757 //_______________________________________________________________________
1758 Float_t AliESDtrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1759 {
1760   //
1761   // TPC cluster information
1762   // type 0: get fraction of found/findable clusters with neighbourhood definition
1763   //      1: findable clusters with neighbourhood definition
1764   //      2: found clusters
1765   //
1766   // definition of findable clusters:
1767   //            a cluster is defined as findable if there is another cluster
1768   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1769   //           effects with a very simple algorithm.
1770   //
1771
1772   if (type==2) return fTPCClusterMap.CountBits();
1773   
1774   Int_t found=0;
1775   Int_t findable=0;
1776   Int_t last=-nNeighbours;
1777   
1778   Int_t upperBound=fTPCClusterMap.GetNbits();
1779   if (upperBound>row1) upperBound=row1;
1780   for (Int_t i=row0; i<upperBound; ++i){
1781     //look to current row
1782     if (fTPCClusterMap[i]) {
1783       last=i;
1784       ++found;
1785       ++findable;
1786       continue;
1787     }
1788     //look to nNeighbours before
1789     if ((i-last)<=nNeighbours) {
1790       ++findable;
1791       continue;
1792     }
1793     //look to nNeighbours after
1794     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1795       if (fTPCClusterMap[j]){
1796         ++findable;
1797         break;
1798       }
1799     }
1800   }
1801   if (type==1) return findable;
1802   
1803   if (type==0){
1804     Float_t fraction=0;
1805     if (findable>0) 
1806       fraction=(Float_t)found/(Float_t)findable;
1807     else 
1808       fraction=0;
1809     return fraction;
1810   }  
1811   return 0;  // undefined type - default value
1812 }
1813
1814 //_______________________________________________________________________
1815 Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
1816   //
1817   // GetDensity of the clusters on given region between row0 and row1
1818   // Dead zone effect takin into acoount
1819   //
1820   if (!fFriendTrack) return 0.0;
1821   Int_t good  = 0;
1822   Int_t found = 0;
1823   //  
1824   Int_t *index=fFriendTrack->GetTPCindices();
1825   for (Int_t i=row0;i<=row1;i++){     
1826     Int_t idx = index[i];
1827     if (idx!=-1)  good++;             // track outside of dead zone
1828     if (idx>0)    found++;
1829   }
1830   Float_t density=0.5;
1831   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
1832   return density;
1833 }
1834
1835 //_______________________________________________________________________
1836 void AliESDtrack::SetTPCpid(const Double_t *p) {  
1837   // Sets values for the probability of each particle type (in TPC)
1838   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
1839   SetStatus(AliESDtrack::kTPCpid);
1840 }
1841
1842 //_______________________________________________________________________
1843 void AliESDtrack::GetTPCpid(Double_t *p) const {
1844   // Gets the probability of each particle type (in TPC)
1845   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
1846 }
1847
1848 //_______________________________________________________________________
1849 UChar_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
1850   //---------------------------------------------------------------------
1851   // This function returns indices of the assgined TRD clusters 
1852   //---------------------------------------------------------------------
1853   if (idx && fFriendTrack) {
1854     Int_t *index=fFriendTrack->GetTRDindices();
1855
1856     if (index) {
1857       for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
1858     }
1859     else {
1860       for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=-2;
1861     }
1862   }
1863   return fTRDncls;
1864 }
1865
1866 //_______________________________________________________________________
1867 UChar_t AliESDtrack::GetTRDtracklets(Int_t *idx) const {
1868 //
1869 // This function returns the number of TRD tracklets used in tracking
1870 // and it fills the indices of these tracklets in the array "idx" as they 
1871 // are registered in the TRD track list. 
1872 // 
1873 // Caution :
1874 //   1. The idx array has to be allocated with a size >= AliESDtrack::kTRDnPlanes
1875 //   2. The idx array store not only the index but also the layer of the tracklet. 
1876 //      Therefore tracks with TRD gaps contain default values for indices [-1] 
1877
1878   if (!fFriendTrack) return 0;
1879   if (!idx) return GetTRDntracklets();
1880   Int_t *index=fFriendTrack->GetTRDindices();
1881   Int_t n = 0;
1882   for (Int_t i=0; i<kTRDnPlanes; i++){ 
1883     if (index){
1884       if(index[i]>=0) n++;
1885       idx[i]=index[i];
1886     }
1887     else idx[i] = -2;
1888   }
1889   return n;
1890 }
1891
1892 //_______________________________________________________________________
1893 void AliESDtrack::SetTRDpid(const Double_t *p) {  
1894   // Sets values for the probability of each particle type (in TRD)
1895   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
1896   SetStatus(AliESDtrack::kTRDpid);
1897 }
1898
1899 //_______________________________________________________________________
1900 void AliESDtrack::GetTRDpid(Double_t *p) const {
1901   // Gets the probability of each particle type (in TRD)
1902   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
1903 }
1904
1905 //_______________________________________________________________________
1906 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
1907 {
1908   // Sets the probability of particle type iSpecies to p (in TRD)
1909   fTRDr[iSpecies] = p;
1910 }
1911
1912 Double_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
1913 {
1914   // Returns the probability of particle type iSpecies (in TRD)
1915   return fTRDr[iSpecies];
1916 }
1917
1918 //____________________________________________________
1919 Int_t AliESDtrack::GetNumberOfTRDslices() const 
1920 {
1921   // built in backward compatibility
1922   Int_t idx = fTRDnSlices - (kTRDnPlanes<<1);
1923   return idx<18 ? fTRDnSlices/kTRDnPlanes : idx/kTRDnPlanes;
1924 }
1925
1926 //____________________________________________________
1927 Double_t AliESDtrack::GetTRDmomentum(Int_t plane, Double_t *sp) const
1928 {
1929 //Returns momentum estimation and optional its error (sp)
1930 // in TRD layer "plane".
1931
1932   if (!fTRDnSlices) {
1933     AliDebug(2, "No TRD info allocated for this track.");
1934     return -1.;
1935   }
1936   if ((plane<0) || (plane>=kTRDnPlanes)) {
1937     AliWarning(Form("Request for TRD plane[%d] outside range.", plane)); 
1938     return -1.;
1939   }
1940
1941   Int_t idx = fTRDnSlices-(kTRDnPlanes<<1)+plane;
1942   // Protection for backward compatibility
1943   if(idx<(GetNumberOfTRDslices()*kTRDnPlanes)) return -1.;
1944
1945   if(sp) (*sp) = fTRDslices[idx+kTRDnPlanes];
1946   return fTRDslices[idx];
1947 }
1948
1949 //____________________________________________________
1950 Double_t  AliESDtrack::GetTRDslice(Int_t plane, Int_t slice) const {
1951   //Gets the charge from the slice of the plane
1952
1953   if(!fTRDslices) {
1954     //AliError("No TRD slices allocated for this track !");
1955     return -1.;
1956   }
1957   if ((plane<0) || (plane>=kTRDnPlanes)) {
1958     AliError("Info for TRD plane not available !");
1959     return -1.;
1960   }
1961   Int_t ns=GetNumberOfTRDslices();
1962   if ((slice<-1) || (slice>=ns)) {
1963     //AliError("Wrong TRD slice !");  
1964     return -1.;
1965   }
1966
1967   if(slice>=0) return fTRDslices[plane*ns + slice];
1968
1969   // return average of the dEdx measurements
1970   Double_t q=0.; Double32_t *s = &fTRDslices[plane*ns];
1971   for (Int_t i=0; i<ns; i++, s++) if((*s)>0.) q+=(*s);
1972   return q/ns;
1973 }
1974
1975 //____________________________________________________
1976 void  AliESDtrack::SetNumberOfTRDslices(Int_t n) {
1977   //Sets the number of slices used for PID 
1978   if (fTRDnSlices) return;
1979
1980   fTRDnSlices=n;
1981   fTRDslices=new Double32_t[fTRDnSlices];
1982   
1983   // set-up correctly the allocated memory
1984   memset(fTRDslices, 0, n*sizeof(Double32_t));
1985   for (Int_t i=GetNumberOfTRDslices(); i--;) fTRDslices[i]=-1.;
1986 }
1987
1988 //____________________________________________________
1989 void  AliESDtrack::SetTRDslice(Double_t q, Int_t plane, Int_t slice) {
1990   //Sets the charge q in the slice of the plane
1991   if(!fTRDslices) {
1992     AliError("No TRD slices allocated for this track !");
1993     return;
1994   }
1995   if ((plane<0) || (plane>=kTRDnPlanes)) {
1996     AliError("Info for TRD plane not allocated !");
1997     return;
1998   }
1999   Int_t ns=GetNumberOfTRDslices();
2000   if ((slice<0) || (slice>=ns)) {
2001     AliError("Wrong TRD slice !");
2002     return;
2003   }
2004   Int_t n=plane*ns + slice;
2005   fTRDslices[n]=q;
2006 }
2007
2008
2009 //____________________________________________________
2010 void AliESDtrack::SetTRDmomentum(Double_t p, Int_t plane, Double_t *sp)
2011 {
2012   if(!fTRDslices) {
2013     AliError("No TRD slices allocated for this track !");
2014     return;
2015   }
2016   if ((plane<0) || (plane>=kTRDnPlanes)) {
2017     AliError("Info for TRD plane not allocated !");
2018     return;
2019   }
2020
2021   Int_t idx = fTRDnSlices-(kTRDnPlanes<<1)+plane;
2022   // Protection for backward compatibility
2023   if(idx<GetNumberOfTRDslices()*kTRDnPlanes) return;
2024
2025   if(sp) fTRDslices[idx+kTRDnPlanes] = (*sp);
2026   fTRDslices[idx] = p;
2027 }
2028
2029
2030 //_______________________________________________________________________
2031 void AliESDtrack::SetTOFpid(const Double_t *p) {  
2032   // Sets the probability of each particle type (in TOF)
2033   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
2034   SetStatus(AliESDtrack::kTOFpid);
2035 }
2036
2037 //_______________________________________________________________________
2038 void AliESDtrack::SetTOFLabel(const Int_t *p) {  
2039   // Sets  (in TOF)
2040   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
2041 }
2042
2043 //_______________________________________________________________________
2044 void AliESDtrack::GetTOFpid(Double_t *p) const {
2045   // Gets probabilities of each particle type (in TOF)
2046   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
2047 }
2048
2049 //_______________________________________________________________________
2050 void AliESDtrack::GetTOFLabel(Int_t *p) const {
2051   // Gets (in TOF)
2052   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
2053 }
2054
2055 //_______________________________________________________________________
2056 void AliESDtrack::GetTOFInfo(Float_t *info) const {
2057   // Gets (in TOF)
2058   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
2059 }
2060
2061 //_______________________________________________________________________
2062 void AliESDtrack::SetTOFInfo(Float_t*info) {
2063   // Gets (in TOF)
2064   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
2065 }
2066
2067
2068
2069 //_______________________________________________________________________
2070 void AliESDtrack::SetHMPIDpid(const Double_t *p) {  
2071   // Sets the probability of each particle type (in HMPID)
2072   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
2073   SetStatus(AliESDtrack::kHMPIDpid);
2074 }
2075
2076 //_______________________________________________________________________
2077 void  AliESDtrack::SetTPCdEdxInfo(AliTPCdEdxInfo * dEdxInfo){ 
2078   if(fTPCdEdxInfo) delete fTPCdEdxInfo;
2079   fTPCdEdxInfo = dEdxInfo; 
2080 }
2081
2082 //_______________________________________________________________________
2083 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
2084   // Gets probabilities of each particle type (in HMPID)
2085   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
2086 }
2087
2088
2089
2090 //_______________________________________________________________________
2091 void AliESDtrack::SetESDpid(const Double_t *p) {  
2092   // Sets the probability of each particle type for the ESD track
2093   SetPIDValues(fR,p,AliPID::kSPECIES);
2094   SetStatus(AliESDtrack::kESDpid);
2095 }
2096
2097 //_______________________________________________________________________
2098 void AliESDtrack::GetESDpid(Double_t *p) const {
2099   // Gets probability of each particle type for the ESD track
2100   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
2101 }
2102
2103 //_______________________________________________________________________
2104 Bool_t AliESDtrack::RelateToVertexTPC(const AliESDVertex *vtx, 
2105 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
2106   //
2107   // Try to relate the TPC-only track parameters to the vertex "vtx", 
2108   // if the (rough) transverse impact parameter is not bigger then "maxd". 
2109   //            Magnetic field is "b" (kG).
2110   //
2111   // a) The TPC-only paramters are extapolated to the DCA to the vertex.
2112   // b) The impact parameters and their covariance matrix are calculated.
2113   // c) An attempt to constrain the TPC-only params to the vertex is done.
2114   //    The constrained params are returned via "cParam".
2115   //
2116   // In the case of success, the returned value is kTRUE
2117   // otherwise, it's kFALSE)
2118   // 
2119
2120   if (!fTPCInner) return kFALSE;
2121   if (!vtx) return kFALSE;
2122
2123   Double_t dz[2],cov[3];
2124   if (!fTPCInner->PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
2125
2126   fdTPC = dz[0];
2127   fzTPC = dz[1];  
2128   fCddTPC = cov[0];
2129   fCdzTPC = cov[1];
2130   fCzzTPC = cov[2];
2131   
2132   Double_t covar[6]; vtx->GetCovMatrix(covar);
2133   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2134   Double_t c[3]={covar[2],0.,covar[5]};
2135
2136   Double_t chi2=GetPredictedChi2(p,c);
2137   if (chi2>kVeryBig) return kFALSE;
2138
2139   fCchi2TPC=chi2;
2140
2141   if (!cParam) return kTRUE;
2142
2143   *cParam = *fTPCInner;
2144   if (!cParam->Update(p,c)) return kFALSE;
2145
2146   return kTRUE;
2147 }
2148
2149 //_______________________________________________________________________
2150 Bool_t AliESDtrack::RelateToVertexTPCBxByBz(const AliESDVertex *vtx, 
2151 Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam) {
2152   //
2153   // Try to relate the TPC-only track parameters to the vertex "vtx", 
2154   // if the (rough) transverse impact parameter is not bigger then "maxd". 
2155   //
2156   // All three components of the magnetic field ,"b[3]" (kG), 
2157   // are taken into account.
2158   //
2159   // a) The TPC-only paramters are extapolated to the DCA to the vertex.
2160   // b) The impact parameters and their covariance matrix are calculated.
2161   // c) An attempt to constrain the TPC-only params to the vertex is done.
2162   //    The constrained params are returned via "cParam".
2163   //
2164   // In the case of success, the returned value is kTRUE
2165   // otherwise, it's kFALSE)
2166   // 
2167
2168   if (!fTPCInner) return kFALSE;
2169   if (!vtx) return kFALSE;
2170
2171   Double_t dz[2],cov[3];
2172   if (!fTPCInner->PropagateToDCABxByBz(vtx, b, maxd, dz, cov)) return kFALSE;
2173
2174   fdTPC = dz[0];
2175   fzTPC = dz[1];  
2176   fCddTPC = cov[0];
2177   fCdzTPC = cov[1];
2178   fCzzTPC = cov[2];
2179   
2180   Double_t covar[6]; vtx->GetCovMatrix(covar);
2181   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2182   Double_t c[3]={covar[2],0.,covar[5]};
2183
2184   Double_t chi2=GetPredictedChi2(p,c);
2185   if (chi2>kVeryBig) return kFALSE;
2186
2187   fCchi2TPC=chi2;
2188
2189   if (!cParam) return kTRUE;
2190
2191   *cParam = *fTPCInner;
2192   if (!cParam->Update(p,c)) return kFALSE;
2193
2194   return kTRUE;
2195 }
2196
2197 //_______________________________________________________________________
2198 Bool_t AliESDtrack::RelateToVertex(const AliESDVertex *vtx, 
2199 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
2200   //
2201   // Try to relate this track to the vertex "vtx", 
2202   // if the (rough) transverse impact parameter is not bigger then "maxd". 
2203   //            Magnetic field is "b" (kG).
2204   //
2205   // a) The track gets extapolated to the DCA to the vertex.
2206   // b) The impact parameters and their covariance matrix are calculated.
2207   // c) An attempt to constrain this track to the vertex is done.
2208   //    The constrained params are returned via "cParam".
2209   //
2210   // In the case of success, the returned value is kTRUE
2211   // (otherwise, it's kFALSE)
2212   //  
2213
2214   if (!vtx) return kFALSE;
2215
2216   Double_t dz[2],cov[3];
2217   if (!PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
2218
2219   fD = dz[0];
2220   fZ = dz[1];  
2221   fCdd = cov[0];
2222   fCdz = cov[1];
2223   fCzz = cov[2];
2224   
2225   Double_t covar[6]; vtx->GetCovMatrix(covar);
2226   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2227   Double_t c[3]={covar[2],0.,covar[5]};
2228
2229   Double_t chi2=GetPredictedChi2(p,c);
2230   if (chi2>kVeryBig) return kFALSE;
2231
2232   fCchi2=chi2;
2233
2234
2235   //--- Could now these lines be removed ? ---
2236   delete fCp;
2237   fCp=new AliExternalTrackParam(*this);  
2238
2239   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
2240   //----------------------------------------
2241
2242   fVertexID = vtx->GetID();
2243
2244   if (!cParam) return kTRUE;
2245
2246   *cParam = *this;
2247   if (!cParam->Update(p,c)) return kFALSE; 
2248
2249   return kTRUE;
2250 }
2251
2252 //_______________________________________________________________________
2253 Bool_t AliESDtrack::RelateToVertexBxByBz(const AliESDVertex *vtx, 
2254 Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam) {
2255   //
2256   // Try to relate this track to the vertex "vtx", 
2257   // if the (rough) transverse impact parameter is not bigger then "maxd". 
2258   //            Magnetic field is "b" (kG).
2259   //
2260   // a) The track gets extapolated to the DCA to the vertex.
2261   // b) The impact parameters and their covariance matrix are calculated.
2262   // c) An attempt to constrain this track to the vertex is done.
2263   //    The constrained params are returned via "cParam".
2264   //
2265   // In the case of success, the returned value is kTRUE
2266   // (otherwise, it's kFALSE)
2267   //  
2268
2269   if (!vtx) return kFALSE;
2270
2271   Double_t dz[2],cov[3];
2272   if (!PropagateToDCABxByBz(vtx, b, maxd, dz, cov)) return kFALSE;
2273
2274   fD = dz[0];
2275   fZ = dz[1];  
2276   fCdd = cov[0];
2277   fCdz = cov[1];
2278   fCzz = cov[2];
2279   
2280   Double_t covar[6]; vtx->GetCovMatrix(covar);
2281   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2282   Double_t c[3]={covar[2],0.,covar[5]};
2283
2284   Double_t chi2=GetPredictedChi2(p,c);
2285   if (chi2>kVeryBig) return kFALSE;
2286
2287   fCchi2=chi2;
2288
2289
2290   //--- Could now these lines be removed ? ---
2291   delete fCp;
2292   fCp=new AliExternalTrackParam(*this);  
2293
2294   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
2295   //----------------------------------------
2296
2297   fVertexID = vtx->GetID();
2298
2299   if (!cParam) return kTRUE;
2300
2301   *cParam = *this;
2302   if (!cParam->Update(p,c)) return kFALSE; 
2303
2304   return kTRUE;
2305 }
2306
2307 //_______________________________________________________________________
2308 void AliESDtrack::Print(Option_t *) const {
2309   // Prints info on the track
2310   AliExternalTrackParam::Print();
2311   printf("ESD track info\n") ; 
2312   Double_t p[AliPID::kSPECIESN] ; 
2313   Int_t index = 0 ; 
2314   if( IsOn(kITSpid) ){
2315     printf("From ITS: ") ; 
2316     GetITSpid(p) ; 
2317     for(index = 0 ; index < AliPID::kSPECIES; index++) 
2318       printf("%f, ", p[index]) ;
2319     printf("\n           signal = %f\n", GetITSsignal()) ;
2320   } 
2321   if( IsOn(kTPCpid) ){
2322     printf("From TPC: ") ; 
2323     GetTPCpid(p) ; 
2324     for(index = 0 ; index < AliPID::kSPECIES; index++) 
2325       printf("%f, ", p[index]) ;
2326     printf("\n           signal = %f\n", GetTPCsignal()) ;
2327   }
2328   if( IsOn(kTRDpid) ){
2329     printf("From TRD: ") ; 
2330     GetTRDpid(p) ; 
2331     for(index = 0 ; index < AliPID::kSPECIES; index++) 
2332       printf("%f, ", p[index]) ;
2333       printf("\n           signal = %f\n", GetTRDsignal()) ;
2334   }
2335   if( IsOn(kTOFpid) ){
2336     printf("From TOF: ") ; 
2337     GetTOFpid(p) ; 
2338     for(index = 0 ; index < AliPID::kSPECIES; index++) 
2339       printf("%f, ", p[index]) ;
2340     printf("\n           signal = %f\n", GetTOFsignal()) ;
2341   }
2342   if( IsOn(kHMPIDpid) ){
2343     printf("From HMPID: ") ; 
2344     GetHMPIDpid(p) ; 
2345     for(index = 0 ; index < AliPID::kSPECIES; index++) 
2346       printf("%f, ", p[index]) ;
2347     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
2348   }
2349
2350
2351
2352 //
2353 // Draw functionality
2354 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
2355 //
2356 void AliESDtrack::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
2357   //
2358   // Fill points in the polymarker
2359   //
2360   TObjArray arrayRef;
2361   arrayRef.AddLast(new AliExternalTrackParam(*this));
2362   if (fIp) arrayRef.AddLast(new AliExternalTrackParam(*fIp));
2363   if (fOp) arrayRef.AddLast(new AliExternalTrackParam(*fOp));
2364   if (fHMPIDp) arrayRef.AddLast(new AliExternalTrackParam(*fHMPIDp));
2365   //
2366   Double_t mpos[3]={0,0,0};
2367   Int_t entries=arrayRef.GetEntries();
2368   for (Int_t i=0;i<entries;i++){
2369     Double_t pos[3];
2370     ((AliExternalTrackParam*)arrayRef.At(i))->GetXYZ(pos);
2371     mpos[0]+=pos[0]/entries;
2372     mpos[1]+=pos[1]/entries;
2373     mpos[2]+=pos[2]/entries;    
2374   }
2375   // Rotate to the mean position
2376   //
2377   Float_t fi= TMath::ATan2(mpos[1],mpos[0]);
2378   for (Int_t i=0;i<entries;i++){
2379     Bool_t res = ((AliExternalTrackParam*)arrayRef.At(i))->Rotate(fi);
2380     if (!res) delete arrayRef.RemoveAt(i);
2381   }
2382   Int_t counter=0;
2383   for (Double_t r=minR; r<maxR; r+=stepR){
2384     Double_t sweight=0;
2385     Double_t mlpos[3]={0,0,0};
2386     for (Int_t i=0;i<entries;i++){
2387       Double_t point[3]={0,0,0};
2388       AliExternalTrackParam *param = ((AliExternalTrackParam*)arrayRef.At(i));
2389       if (!param) continue;
2390       if (param->GetXYZAt(r,magF,point)){
2391         Double_t weight = 1./(10.+(r-param->GetX())*(r-param->GetX()));
2392         sweight+=weight;
2393         mlpos[0]+=point[0]*weight;
2394         mlpos[1]+=point[1]*weight;
2395         mlpos[2]+=point[2]*weight;
2396       }
2397     }
2398     if (sweight>0){
2399       mlpos[0]/=sweight;
2400       mlpos[1]/=sweight;
2401       mlpos[2]/=sweight;      
2402       pol->SetPoint(counter,mlpos[0],mlpos[1], mlpos[2]);
2403       //      printf("xyz\t%f\t%f\t%f\n",mlpos[0], mlpos[1],mlpos[2]);
2404       counter++;
2405     }
2406   }
2407 }
2408
2409 //_______________________________________________________________________
2410 void AliESDtrack::SetITSdEdxSamples(const Double_t s[4]) {
2411   //
2412   // Store the dE/dx samples measured by the two SSD and two SDD layers.
2413   // These samples are corrected for the track segment length. 
2414   //
2415   for (Int_t i=0; i<4; i++) fITSdEdxSamples[i]=s[i];
2416 }
2417
2418 //_______________________________________________________________________
2419 void AliESDtrack::GetITSdEdxSamples(Double_t *s) const {
2420   //
2421   // Get the dE/dx samples measured by the two SSD and two SDD layers.  
2422   // These samples are corrected for the track segment length.
2423   //
2424   for (Int_t i=0; i<4; i++) s[i]=fITSdEdxSamples[i];
2425 }
2426
2427
2428 UShort_t   AliESDtrack::GetTPCnclsS(Int_t i0,Int_t i1) const{
2429   //
2430   // get number of shared TPC clusters
2431   //
2432   return  fTPCSharedMap.CountBits(i0)-fTPCSharedMap.CountBits(i1);
2433 }
2434
2435 UShort_t   AliESDtrack::GetTPCncls(Int_t i0,Int_t i1) const{
2436   //
2437   // get number of TPC clusters
2438   //
2439   return  fTPCClusterMap.CountBits(i0)-fTPCClusterMap.CountBits(i1);
2440 }
2441
2442 //____________________________________________________________________
2443 Double_t AliESDtrack::GetChi2TPCConstrainedVsGlobal(const AliESDVertex* vtx) const
2444 {
2445   // Calculates the chi2 between the TPC track (TPCinner) constrained to the primary vertex and the global track
2446   //
2447   // Returns -1 in case the calculation failed
2448   //
2449   // Value is cached as a non-persistent member.
2450   //
2451   // Code adapted from original code by GSI group (Jacek, Marian, Michael)
2452   
2453   // cache, ignoring that a different vertex might be passed
2454   if (fCacheChi2TPCConstrainedVsGlobalVertex == vtx)
2455     return fCacheChi2TPCConstrainedVsGlobal;
2456   
2457   fCacheChi2TPCConstrainedVsGlobal = -1;
2458   fCacheChi2TPCConstrainedVsGlobalVertex = vtx;
2459   
2460   Double_t x[3];
2461   GetXYZ(x);
2462   Double_t b[3];
2463   AliTrackerBase::GetBxByBz(x,b);
2464
2465   if (!fTPCInner)  { 
2466     AliWarning("Could not get TPC Inner Param.");
2467     return fCacheChi2TPCConstrainedVsGlobal;
2468   }
2469   
2470   // clone for constraining
2471   AliExternalTrackParam* tpcInnerC = new AliExternalTrackParam(*fTPCInner);
2472   if (!tpcInnerC) { 
2473     AliWarning("Clone of TPCInnerParam failed.");
2474     return fCacheChi2TPCConstrainedVsGlobal;  
2475   }
2476   
2477   // transform to the track reference frame 
2478   Bool_t isOK = tpcInnerC->Rotate(GetAlpha());
2479   isOK &= tpcInnerC->PropagateTo(GetX(), b[2]);
2480   if (!isOK) { 
2481     delete tpcInnerC;
2482     tpcInnerC = 0; 
2483     AliWarning("Rotation/Propagation of track failed.") ; 
2484     return fCacheChi2TPCConstrainedVsGlobal;    
2485   }  
2486
2487   // constrain TPCinner 
2488   isOK = tpcInnerC->ConstrainToVertex(vtx, b);
2489   
2490   // transform to the track reference frame 
2491   isOK &= tpcInnerC->Rotate(GetAlpha());
2492   isOK &= tpcInnerC->PropagateTo(GetX(), b[2]);
2493
2494   if (!isOK) {
2495     AliWarning("ConstrainTPCInner failed.") ;
2496     delete tpcInnerC;
2497     tpcInnerC = 0; 
2498     return fCacheChi2TPCConstrainedVsGlobal;  
2499   }
2500   
2501   // calculate chi2 between vi and vj vectors
2502   // with covi and covj covariance matrices
2503   // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)
2504   TMatrixD deltaT(5,1);
2505   TMatrixD delta(1,5);
2506   TMatrixD covarM(5,5);
2507
2508   for (Int_t ipar=0; ipar<5; ipar++) {
2509     deltaT(ipar,0) = tpcInnerC->GetParameter()[ipar] - GetParameter()[ipar];
2510     delta(0,ipar) = tpcInnerC->GetParameter()[ipar] - GetParameter()[ipar];
2511
2512     for (Int_t jpar=0; jpar<5; jpar++) {
2513       Int_t index = GetIndex(ipar,jpar);
2514       covarM(ipar,jpar) = GetCovariance()[index]+tpcInnerC->GetCovariance()[index];
2515     }
2516   }
2517   // chi2 distance TPC constrained and TPC+ITS
2518   TMatrixD covarMInv = covarM.Invert();
2519   TMatrixD mat2 = covarMInv*deltaT;
2520   TMatrixD chi2 = delta*mat2; 
2521   
2522   delete tpcInnerC; 
2523   tpcInnerC = 0;
2524   
2525   fCacheChi2TPCConstrainedVsGlobal = chi2(0,0);
2526   return fCacheChi2TPCConstrainedVsGlobal;
2527 }