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