]> git.uio.no Git - u/mrichter/AliRoot.git/blob - STEER/AliESDtrack.cxx
Next step towards supporting multiple primary vertices in ESD. Every primary vertex...
[u/mrichter/AliRoot.git] / STEER / AliESDtrack.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15 //-----------------------------------------------------------------
16 //           Implementation of the ESD track class
17 //   ESD = Event Summary Data
18 //   This is the class to deal with during the phisics analysis of data
19 //      Origin: Iouri Belikov, CERN
20 //      e-mail: Jouri.Belikov@cern.ch
21 //
22 //
23 //
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
119 #include "AliESDVertex.h"
120 #include "AliESDtrack.h"
121 #include "AliKalmanTrack.h"
122 #include "AliVTrack.h"
123 #include "AliLog.h"
124 #include "AliTrackPointArray.h"
125 #include "TPolyMarker3D.h"
126
127 ClassImp(AliESDtrack)
128
129 void SetPIDValues(Double_t * dest, const Double_t * src, Int_t n) {
130   // This function copies "n" PID weights from "scr" to "dest"
131   // and normalizes their sum to 1 thus producing conditional probabilities.
132   // The negative weights are set to 0.
133   // In case all the weights are non-positive they are replaced by
134   // uniform probabilities
135
136   if (n<=0) return;
137
138   Float_t uniform = 1./(Float_t)n;
139
140   Float_t sum = 0;
141   for (Int_t i=0; i<n; i++) 
142     if (src[i]>=0) {
143       sum+=src[i];
144       dest[i] = src[i];
145     }
146     else {
147       dest[i] = 0;
148     }
149
150   if(sum>0)
151     for (Int_t i=0; i<n; i++) dest[i] /= sum;
152   else
153     for (Int_t i=0; i<n; i++) dest[i] = uniform;
154 }
155
156 //_______________________________________________________________________
157 AliESDtrack::AliESDtrack() : 
158   AliExternalTrackParam(),
159   fCp(0),
160   fIp(0),
161   fTPCInner(0),
162   fOp(0),
163   fFriendTrack(new AliESDfriendTrack()),
164   fTPCClusterMap(159),//number of padrows
165   fTPCSharedMap(159),//number of padrows
166   fFlags(0),
167   fID(0),
168   fLabel(0),
169   fITSLabel(0),
170   fTPCLabel(0),
171   fTRDLabel(0),
172   fTOFCalChannel(0),
173   fTOFindex(-1),
174   fHMPIDqn(0),
175   fHMPIDcluIdx(-1),
176   fEMCALindex(kEMCALNoMatch),
177   fHMPIDtrkTheta(0),
178   fHMPIDtrkPhi(0),
179   fHMPIDsignal(0),
180   fTrackLength(0),
181   fdTPC(0),fzTPC(0),
182   fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
183   fCchi2TPC(0),
184   fD(0),fZ(0),
185   fCdd(0),fCdz(0),fCzz(0),
186   fCchi2(0),
187   fITSchi2(0),
188   fTPCchi2(0),
189   fTRDchi2(0),
190   fTOFchi2(0),
191   fHMPIDchi2(0),
192   fGlobalChi2(0),
193   fITSsignal(0),
194   fTPCsignal(0),
195   fTPCsignalS(0),
196   fTRDsignal(0),
197   fTRDQuality(0),
198   fTRDBudget(0),
199   fTOFsignal(0),
200   fTOFsignalToT(0),
201   fTOFsignalRaw(0),
202   fTOFsignalDz(0),
203   fHMPIDtrkX(0),
204   fHMPIDtrkY(0),
205   fHMPIDmipX(0),
206   fHMPIDmipY(0),
207   fTPCncls(0),
208   fTPCnclsF(0),
209   fTPCsignalN(0),
210   fITSncls(0),
211   fITSClusterMap(0),
212   fTRDncls(0),
213   fTRDncls0(0),
214   fTRDntracklets(0),
215   fTRDnSlices(0),
216   fTRDslices(0x0),
217   fVertexID(-2) // -2 means an orphan track 
218 {
219   //
220   // The default ESD constructor 
221   //
222   Int_t i;
223   for (i=0; i<AliPID::kSPECIES; i++) {
224     fTrackTime[i]=0.;
225     fR[i]=0.;
226     fITSr[i]=0.;
227     fTPCr[i]=0.;
228     fTRDr[i]=0.;
229     fTOFr[i]=0.;
230     fHMPIDr[i]=0.;
231   }
232   
233   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
234   for (i=0; i<3; i++)   { fV0Indexes[i]=0;}
235   for (i=0;i<kTRDnPlanes;i++) {
236     fTRDTimBin[i]=0;
237   }
238   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
239   for (i=0;i<3;i++) {fTOFLabel[i]=0;}
240   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
241   for (i=0;i<12;i++) {fITSModule[i]=-1;}
242 }
243
244 //_______________________________________________________________________
245 AliESDtrack::AliESDtrack(const AliESDtrack& track):
246   AliExternalTrackParam(track),
247   fCp(0),
248   fIp(0),
249   fTPCInner(0),
250   fOp(0),
251   fFriendTrack(0),
252   fTPCClusterMap(track.fTPCClusterMap),
253   fTPCSharedMap(track.fTPCSharedMap),
254   fFlags(track.fFlags),
255   fID(track.fID),
256   fLabel(track.fLabel),
257   fITSLabel(track.fITSLabel),
258   fTPCLabel(track.fTPCLabel),
259   fTRDLabel(track.fTRDLabel),
260   fTOFCalChannel(track.fTOFCalChannel),
261   fTOFindex(track.fTOFindex),
262   fHMPIDqn(track.fHMPIDqn),
263   fHMPIDcluIdx(track.fHMPIDcluIdx),
264   fEMCALindex(track.fEMCALindex),
265   fHMPIDtrkTheta(track.fHMPIDtrkTheta),
266   fHMPIDtrkPhi(track.fHMPIDtrkPhi),
267   fHMPIDsignal(track.fHMPIDsignal),
268   fTrackLength(track.fTrackLength),
269   fdTPC(track.fdTPC),fzTPC(track.fzTPC),
270   fCddTPC(track.fCddTPC),fCdzTPC(track.fCdzTPC),fCzzTPC(track.fCzzTPC),
271   fCchi2TPC(track.fCchi2TPC),
272   fD(track.fD),fZ(track.fZ),
273   fCdd(track.fCdd),fCdz(track.fCdz),fCzz(track.fCzz),
274   fCchi2(track.fCchi2),
275   fITSchi2(track.fITSchi2),
276   fTPCchi2(track.fTPCchi2),
277   fTRDchi2(track.fTRDchi2),
278   fTOFchi2(track.fTOFchi2),
279   fHMPIDchi2(track.fHMPIDchi2),
280   fGlobalChi2(track.fGlobalChi2),
281   fITSsignal(track.fITSsignal),
282   fTPCsignal(track.fTPCsignal),
283   fTPCsignalS(track.fTPCsignalS),
284   fTRDsignal(track.fTRDsignal),
285   fTRDQuality(track.fTRDQuality),
286   fTRDBudget(track.fTRDBudget),
287   fTOFsignal(track.fTOFsignal),
288   fTOFsignalToT(track.fTOFsignalToT),
289   fTOFsignalRaw(track.fTOFsignalRaw),
290   fTOFsignalDz(track.fTOFsignalDz),
291   fHMPIDtrkX(track.fHMPIDtrkX),
292   fHMPIDtrkY(track.fHMPIDtrkY),
293   fHMPIDmipX(track.fHMPIDmipX),
294   fHMPIDmipY(track.fHMPIDmipY),
295   fTPCncls(track.fTPCncls),
296   fTPCnclsF(track.fTPCnclsF),
297   fTPCsignalN(track.fTPCsignalN),
298   fITSncls(track.fITSncls),
299   fITSClusterMap(track.fITSClusterMap),
300   fTRDncls(track.fTRDncls),
301   fTRDncls0(track.fTRDncls0),
302   fTRDntracklets(track.fTRDntracklets),
303   fTRDnSlices(track.fTRDnSlices),
304   fTRDslices(0x0),
305   fVertexID(track.fVertexID)
306 {
307   //
308   //copy constructor
309   //
310   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i]=track.fTrackTime[i];
311   for (Int_t i=0;i<AliPID::kSPECIES;i++)  fR[i]=track.fR[i];
312   //
313   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=track.fITSr[i]; 
314   //
315   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=track.fTPCr[i]; 
316   for (Int_t i=0;i<4;i++) {fTPCPoints[i]=track.fTPCPoints[i];}
317   for (Int_t i=0; i<3;i++)   { fKinkIndexes[i]=track.fKinkIndexes[i];}
318   for (Int_t i=0; i<3;i++)   { fV0Indexes[i]=track.fV0Indexes[i];}
319   //
320   for (Int_t i=0;i<kTRDnPlanes;i++) {
321     fTRDTimBin[i]=track.fTRDTimBin[i];
322   }
323
324   if (fTRDnSlices) {
325     fTRDslices=new Double32_t[fTRDnSlices];
326     for (Int_t i=0; i<fTRDnSlices; i++) fTRDslices[i]=track.fTRDslices[i];
327   }
328
329   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i]=track.fTRDr[i]; 
330   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i]=track.fTOFr[i];
331   for (Int_t i=0;i<3;i++) fTOFLabel[i]=track.fTOFLabel[i];
332   for (Int_t i=0;i<10;i++) fTOFInfo[i]=track.fTOFInfo[i];
333   for (Int_t i=0;i<12;i++) fITSModule[i]=track.fITSModule[i];
334   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i]=track.fHMPIDr[i];
335
336   if (track.fCp) fCp=new AliExternalTrackParam(*track.fCp);
337   if (track.fIp) fIp=new AliExternalTrackParam(*track.fIp);
338   if (track.fTPCInner) fTPCInner=new AliExternalTrackParam(*track.fTPCInner);
339   if (track.fOp) fOp=new AliExternalTrackParam(*track.fOp);
340
341   if (track.fFriendTrack) fFriendTrack=new AliESDfriendTrack(*(track.fFriendTrack));
342 }
343
344 //_______________________________________________________________________
345 AliESDtrack::AliESDtrack(const AliVTrack *track) : 
346   AliExternalTrackParam(track),
347   fCp(0),
348   fIp(0),
349   fTPCInner(0),
350   fOp(0),
351   fFriendTrack(0),
352   fTPCClusterMap(159),//number of padrows
353   fTPCSharedMap(159),//number of padrows
354   fFlags(0),
355   fID(),
356   fLabel(0),
357   fITSLabel(0),
358   fTPCLabel(0),
359   fTRDLabel(0),
360   fTOFCalChannel(0),
361   fTOFindex(-1),
362   fHMPIDqn(0),
363   fHMPIDcluIdx(-1),
364   fEMCALindex(kEMCALNoMatch),
365   fHMPIDtrkTheta(0),
366   fHMPIDtrkPhi(0),
367   fHMPIDsignal(0),
368   fTrackLength(0),
369   fdTPC(0),fzTPC(0),
370   fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
371   fCchi2TPC(0),
372   fD(0),fZ(0),
373   fCdd(0),fCdz(0),fCzz(0),
374   fCchi2(0),
375   fITSchi2(0),
376   fTPCchi2(0),
377   fTRDchi2(0),
378   fTOFchi2(0),
379   fHMPIDchi2(0),
380   fGlobalChi2(0),
381   fITSsignal(0),
382   fTPCsignal(0),
383   fTPCsignalS(0),
384   fTRDsignal(0),
385   fTRDQuality(0),
386   fTRDBudget(0),
387   fTOFsignal(0),
388   fTOFsignalToT(0),
389   fTOFsignalRaw(0),
390   fTOFsignalDz(0),
391   fHMPIDtrkX(0),
392   fHMPIDtrkY(0),
393   fHMPIDmipX(0),
394   fHMPIDmipY(0),
395   fTPCncls(0),
396   fTPCnclsF(0),
397   fTPCsignalN(0),
398   fITSncls(0),
399   fITSClusterMap(0),
400   fTRDncls(0),
401   fTRDncls0(0),
402   fTRDntracklets(0),
403   fTRDnSlices(0),
404   fTRDslices(0x0),
405   fVertexID(-2)  // -2 means an orphan track
406 {
407   //
408   // ESD track from AliVTrack.
409   // This is not a copy constructor !
410   //
411
412   if (track->InheritsFrom("AliExternalTrackParam")) {
413      AliError("This is not a copy constructor. Use AliESDtrack(const AliESDtrack &) !");
414      AliWarning("Calling the default constructor...");
415      AliESDtrack();
416      return;
417   }
418
419   // Reset all the arrays
420   Int_t i;
421   for (i=0; i<AliPID::kSPECIES; i++) {
422     fTrackTime[i]=0.;
423     fR[i]=0.;
424     fITSr[i]=0.;
425     fTPCr[i]=0.;
426     fTRDr[i]=0.;
427     fTOFr[i]=0.;
428     fHMPIDr[i]=0.;
429   }
430   
431   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
432   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
433   for (i=0;i<kTRDnPlanes;i++) {
434     fTRDTimBin[i]=0;
435   }
436   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
437   for (i=0;i<3;i++) {fTOFLabel[i]=0;}
438   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
439   for (i=0;i<12;i++) {fITSModule[i]=-1;}
440
441   // Set the ID
442   SetID(track->GetID());
443
444   // Set ITS cluster map
445   fITSClusterMap=track->GetITSClusterMap();
446
447   fITSncls=0;
448   for(i=0; i<6; i++) {
449     if(HasPointOnITSLayer(i)) fITSncls++;
450   }
451
452   // Set the combined PID
453   const Double_t *pid = track->PID();
454   if(pid){
455     for (i=0; i<AliPID::kSPECIES; i++) fR[i]=pid[i];
456   }
457   // AliESD track label
458   SetLabel(track->GetLabel());
459   // Set the status
460   SetStatus(track->GetStatus());
461 }
462
463 //_______________________________________________________________________
464 AliESDtrack::AliESDtrack(TParticle * part) : 
465   AliExternalTrackParam(),
466   fCp(0),
467   fIp(0),
468   fTPCInner(0),
469   fOp(0),
470   fFriendTrack(0),
471   fTPCClusterMap(159),//number of padrows
472   fTPCSharedMap(159),//number of padrows
473   fFlags(0),
474   fID(0),
475   fLabel(0),
476   fITSLabel(0),
477   fTPCLabel(0),
478   fTRDLabel(0),
479   fTOFCalChannel(0),
480   fTOFindex(-1),
481   fHMPIDqn(0),
482   fHMPIDcluIdx(-1),
483   fEMCALindex(kEMCALNoMatch),
484   fHMPIDtrkTheta(0),
485   fHMPIDtrkPhi(0),
486   fHMPIDsignal(0),
487   fTrackLength(0),
488   fdTPC(0),fzTPC(0),
489   fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
490   fCchi2TPC(0),
491   fD(0),fZ(0),
492   fCdd(0),fCdz(0),fCzz(0),
493   fCchi2(0),
494   fITSchi2(0),
495   fTPCchi2(0),
496   fTRDchi2(0),
497   fTOFchi2(0),
498   fHMPIDchi2(0),
499   fGlobalChi2(0),
500   fITSsignal(0),
501   fTPCsignal(0),
502   fTPCsignalS(0),
503   fTRDsignal(0),
504   fTRDQuality(0),
505   fTRDBudget(0),
506   fTOFsignal(0),
507   fTOFsignalToT(0),
508   fTOFsignalRaw(0),
509   fTOFsignalDz(0),
510   fHMPIDtrkX(0),
511   fHMPIDtrkY(0),
512   fHMPIDmipX(0),
513   fHMPIDmipY(0),
514   fTPCncls(0),
515   fTPCnclsF(0),
516   fTPCsignalN(0),
517   fITSncls(0),
518   fITSClusterMap(0),
519   fTRDncls(0),
520   fTRDncls0(0),
521   fTRDntracklets(0),
522   fTRDnSlices(0),
523   fTRDslices(0x0),
524   fVertexID(-2)  // -2 means an orphan track
525 {
526   //
527   // ESD track from TParticle
528   //
529
530   // Reset all the arrays
531   Int_t i;
532   for (i=0; i<AliPID::kSPECIES; i++) {
533     fTrackTime[i]=0.;
534     fR[i]=0.;
535     fITSr[i]=0.;
536     fTPCr[i]=0.;
537     fTRDr[i]=0.;
538     fTOFr[i]=0.;
539     fHMPIDr[i]=0.;
540   }
541   
542   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
543   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
544   for (i=0;i<kTRDnPlanes;i++) {
545     fTRDTimBin[i]=0;
546   }
547   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
548   for (i=0;i<3;i++) {fTOFLabel[i]=0;}
549   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
550   for (i=0;i<12;i++) {fITSModule[i]=-1;}
551
552   // Calculate the AliExternalTrackParam content
553
554   Double_t xref;
555   Double_t alpha;
556   Double_t param[5];
557   Double_t covar[15];
558
559   // Calculate alpha: the rotation angle of the corresponding local system (TPC sector)
560   alpha = part->Phi()*180./TMath::Pi();
561   if (alpha<0) alpha+= 360.;
562   if (alpha>360) alpha -= 360.;
563
564   Int_t sector = (Int_t)(alpha/20.);
565   alpha = 10. + 20.*sector;
566   alpha /= 180;
567   alpha *= TMath::Pi();
568
569   // Covariance matrix: no errors, the parameters are exact
570   for (i=0; i<15; i++) covar[i]=0.;
571
572   // Get the vertex of origin and the momentum
573   TVector3 ver(part->Vx(),part->Vy(),part->Vz());
574   TVector3 mom(part->Px(),part->Py(),part->Pz());
575
576   // Rotate to the local coordinate system (TPC sector)
577   ver.RotateZ(-alpha);
578   mom.RotateZ(-alpha);
579
580   // X of the referense plane
581   xref = ver.X();
582
583   Int_t pdgCode = part->GetPdgCode();
584
585   Double_t charge = 
586     TDatabasePDG::Instance()->GetParticle(pdgCode)->Charge();
587
588   param[0] = ver.Y();
589   param[1] = ver.Z();
590   param[2] = TMath::Sin(mom.Phi());
591   param[3] = mom.Pz()/mom.Pt();
592   param[4] = TMath::Sign(1/mom.Pt(),charge);
593
594   // Set AliExternalTrackParam
595   Set(xref, alpha, param, covar);
596
597   // Set the PID
598   Int_t indexPID = 99;
599
600   switch (TMath::Abs(pdgCode)) {
601
602   case  11: // electron
603     indexPID = 0;
604     break;
605
606   case 13: // muon
607     indexPID = 1;
608     break;
609
610   case 211: // pion
611     indexPID = 2;
612     break;
613
614   case 321: // kaon
615     indexPID = 3;
616     break;
617
618   case 2212: // proton
619     indexPID = 4;
620     break;
621
622   default:
623     break;
624   }
625
626   // If the particle is not e,mu,pi,K or p the PID probabilities are set to 0
627   if (indexPID < AliPID::kSPECIES) {
628     fR[indexPID]=1.;
629     fITSr[indexPID]=1.;
630     fTPCr[indexPID]=1.;
631     fTRDr[indexPID]=1.;
632     fTOFr[indexPID]=1.;
633     fHMPIDr[indexPID]=1.;
634
635   }
636   // AliESD track label
637   SetLabel(part->GetUniqueID());
638
639 }
640
641 //_______________________________________________________________________
642 AliESDtrack::~AliESDtrack(){ 
643   //
644   // This is destructor according Coding Conventrions 
645   //
646   //printf("Delete track\n");
647   delete fIp; 
648   delete fTPCInner; 
649   delete fOp;
650   delete fCp; 
651   delete fFriendTrack;
652   if(fTRDnSlices)
653     delete[] fTRDslices;
654 }
655
656 AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
657   
658
659   if(&source == this) return *this;
660   AliExternalTrackParam::operator=(source);
661
662   
663   if(source.fCp){
664     // we have the trackparam: assign or copy construct
665     if(fCp)*fCp = *source.fCp;
666     else fCp = new AliExternalTrackParam(*source.fCp);
667   }
668   else{
669     // no track param delete the old one
670     if(fCp)delete fCp;
671     fCp = 0;
672   }
673
674   if(source.fIp){
675     // we have the trackparam: assign or copy construct
676     if(fIp)*fIp = *source.fIp;
677     else fIp = new AliExternalTrackParam(*source.fIp);
678   }
679   else{
680     // no track param delete the old one
681     if(fIp)delete fIp;
682     fIp = 0;
683   }
684
685
686   if(source.fTPCInner){
687     // we have the trackparam: assign or copy construct
688     if(fTPCInner) *fTPCInner = *source.fTPCInner;
689     else fTPCInner = new AliExternalTrackParam(*source.fTPCInner);
690   }
691   else{
692     // no track param delete the old one
693     if(fTPCInner)delete fTPCInner;
694     fTPCInner = 0;
695   }
696
697
698   if(source.fOp){
699     // we have the trackparam: assign or copy construct
700     if(fOp) *fOp = *source.fOp;
701     else fOp = new AliExternalTrackParam(*source.fOp);
702   }
703   else{
704     // no track param delete the old one
705     if(fOp)delete fOp;
706     fOp = 0;
707   }
708
709   // copy also the friend track 
710   // use copy constructor
711   if(source.fFriendTrack){
712     // we have the trackparam: assign or copy construct
713     delete fFriendTrack; fFriendTrack=new AliESDfriendTrack(*source.fFriendTrack);
714   }
715   else{
716     // no track param delete the old one
717     delete fFriendTrack; fFriendTrack= 0;
718   }
719
720   fTPCClusterMap = source.fTPCClusterMap; 
721   fTPCSharedMap  = source.fTPCSharedMap;  
722   // the simple stuff
723   fFlags    = source.fFlags; 
724   fID       = source.fID;             
725   fLabel    = source.fLabel;
726   fITSLabel = source.fITSLabel;
727   for(int i = 0; i< 12;++i){
728     fITSModule[i] = source.fITSModule[i];
729   }
730   fTPCLabel = source.fTPCLabel; 
731   fTRDLabel = source.fTRDLabel;
732   for(int i = 0; i< 3;++i){
733     fTOFLabel[i] = source.fTOFLabel[i];    
734   }
735   fTOFCalChannel = source.fTOFCalChannel;
736   fTOFindex      = source.fTOFindex;
737   fHMPIDqn       = source.fHMPIDqn;
738   fHMPIDcluIdx   = source.fHMPIDcluIdx; 
739   fEMCALindex    = source.fEMCALindex;
740
741   for(int i = 0; i< 3;++i){
742     fKinkIndexes[i] = source.fKinkIndexes[i]; 
743     fV0Indexes[i]   = source.fV0Indexes[i]; 
744   }
745
746   for(int i = 0; i< AliPID::kSPECIES;++i){
747     fR[i]     = source.fR[i];
748     fITSr[i]  = source.fITSr[i];
749     fTPCr[i]  = source.fTPCr[i];
750     fTRDr[i]  = source.fTRDr[i];
751     fTOFr[i]  = source.fTOFr[i];
752     fHMPIDr[i] = source.fHMPIDr[i];
753     fTrackTime[i] = source.fTrackTime[i];  
754   }
755
756   fHMPIDtrkTheta = source.fHMPIDtrkTheta;
757   fHMPIDtrkPhi   = source.fHMPIDtrkPhi;
758   fHMPIDsignal   = source.fHMPIDsignal; 
759
760   
761   fTrackLength   = source. fTrackLength;
762   fdTPC  = source.fdTPC; 
763   fzTPC  = source.fzTPC; 
764   fCddTPC = source.fCddTPC;
765   fCdzTPC = source.fCdzTPC;
766   fCzzTPC = source.fCzzTPC;
767   fCchi2TPC = source.fCchi2TPC;
768
769   fD  = source.fD; 
770   fZ  = source.fZ; 
771   fCdd = source.fCdd;
772   fCdz = source.fCdz;
773   fCzz = source.fCzz;
774   fCchi2     = source.fCchi2;
775
776   fITSchi2   = source.fITSchi2;             
777   fTPCchi2   = source.fTPCchi2;            
778   fTRDchi2   = source.fTRDchi2;      
779   fTOFchi2   = source.fTOFchi2;      
780   fHMPIDchi2 = source.fHMPIDchi2;      
781
782   fGlobalChi2 = source.fGlobalChi2;      
783
784   fITSsignal  = source.fITSsignal;     
785   fTPCsignal  = source.fTPCsignal;     
786   fTPCsignalS = source.fTPCsignalS;    
787   for(int i = 0; i< 4;++i){
788     fTPCPoints[i] = source.fTPCPoints[i];  
789   }
790   fTRDsignal = source.fTRDsignal;
791
792   for(int i = 0;i < kTRDnPlanes;++i){
793     fTRDTimBin[i] = source.fTRDTimBin[i];   
794   }
795
796   if(fTRDnSlices)
797     delete[] fTRDslices;
798   fTRDslices=0;
799   fTRDnSlices=source.fTRDnSlices;
800   if (fTRDnSlices) {
801     fTRDslices=new Double32_t[fTRDnSlices];
802     for(int j = 0;j < fTRDnSlices;++j) fTRDslices[j] = source.fTRDslices[j];
803   }
804
805   fTRDQuality =   source.fTRDQuality;     
806   fTRDBudget  =   source.fTRDBudget;      
807   fTOFsignal  =   source.fTOFsignal;     
808   fTOFsignalToT = source.fTOFsignalToT;   
809   fTOFsignalRaw = source.fTOFsignalRaw;  
810   fTOFsignalDz  = source.fTOFsignalDz;      
811   
812   for(int i = 0;i<10;++i){
813     fTOFInfo[i] = source.fTOFInfo[i];    
814   }
815
816   fHMPIDtrkX = source.fHMPIDtrkX; 
817   fHMPIDtrkY = source.fHMPIDtrkY; 
818   fHMPIDmipX = source.fHMPIDmipX;
819   fHMPIDmipY = source.fHMPIDmipY; 
820
821   fTPCncls    = source.fTPCncls;      
822   fTPCnclsF   = source.fTPCnclsF;     
823   fTPCsignalN = source.fTPCsignalN;   
824
825   fITSncls = source.fITSncls;       
826   fITSClusterMap = source.fITSClusterMap; 
827   fTRDncls   = source.fTRDncls;       
828   fTRDncls0  = source.fTRDncls0;      
829   fTRDntracklets  = source.fTRDntracklets; 
830   fVertexID = source.fVertexID;
831   return *this;
832 }
833
834
835
836 void AliESDtrack::Copy(TObject &obj) const {
837   
838   // this overwrites the virtual TOBject::Copy()
839   // to allow run time copying without casting
840   // in AliESDEvent
841
842   if(this==&obj)return;
843   AliESDtrack *robj = dynamic_cast<AliESDtrack*>(&obj);
844   if(!robj)return; // not an AliESDtrack
845   *robj = *this;
846
847 }
848
849
850
851 void AliESDtrack::AddCalibObject(TObject * object){
852   //
853   // add calib object to the list
854   //
855   if (!fFriendTrack) fFriendTrack  = new AliESDfriendTrack;
856   fFriendTrack->AddCalibObject(object);
857 }
858
859 TObject *  AliESDtrack::GetCalibObject(Int_t index){
860   //
861   // return calib objct at given position
862   //
863   if (!fFriendTrack) return 0;
864   return fFriendTrack->GetCalibObject(index);
865 }
866
867
868 Bool_t AliESDtrack::FillTPCOnlyTrack(AliESDtrack &track){
869   
870   // Fills the information of the TPC-only first reconstruction pass
871   // into the passed ESDtrack object. For consistency fTPCInner is also filled
872   // again
873
874
875
876   // For data produced before r26675
877   // RelateToVertexTPC was not properly called during reco
878   // so you'll have to call it again, before FillTPCOnlyTrack
879   //  Float_t p[2],cov[3];
880   // track->GetImpactParametersTPC(p,cov); 
881   // if(p[0]==0&&p[1]==0) // <- Default values
882   //  track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
883   
884
885   if(!fTPCInner)return kFALSE;
886
887   // fill the TPC track params to the global track parameters
888   track.Set(fTPCInner->GetX(),fTPCInner->GetAlpha(),fTPCInner->GetParameter(),fTPCInner->GetCovariance());
889   track.fD = fdTPC;
890   track.fZ = fzTPC;
891   track.fCdd = fCddTPC;
892   track.fCdz = fCdzTPC;
893   track.fCzz = fCzzTPC;
894
895   // copy the TPCinner parameters
896   if(track.fTPCInner) *track.fTPCInner = *fTPCInner;
897   else track.fTPCInner = new AliExternalTrackParam(*fTPCInner);
898   track.fdTPC   = fdTPC;
899   track.fzTPC   = fzTPC;
900   track.fCddTPC = fCddTPC;
901   track.fCdzTPC = fCdzTPC;
902   track.fCzzTPC = fCzzTPC;
903   track.fCchi2TPC = fCchi2TPC;
904
905
906   // copy all other TPC specific parameters
907
908   // replace label by TPC label
909   track.fLabel    = fTPCLabel;
910   track.fTPCLabel = fTPCLabel;
911
912   track.fTPCchi2 = fTPCchi2; 
913   track.fTPCsignal = fTPCsignal;
914   track.fTPCsignalS = fTPCsignalS;
915   for(int i = 0;i<4;++i)track.fTPCPoints[i] = fTPCPoints[i];
916
917   track.fTPCncls    = fTPCncls;     
918   track.fTPCnclsF   = fTPCnclsF;     
919   track.fTPCsignalN =  fTPCsignalN;
920
921   // PID 
922   for(int i=0;i<AliPID::kSPECIES;++i){
923     track.fTPCr[i] = fTPCr[i];
924     // combined PID is TPC only!
925     track.fR[i] = fTPCr[i];
926   }
927   track.fTPCClusterMap = fTPCClusterMap;
928   track.fTPCSharedMap = fTPCSharedMap;
929
930
931   // reset the flags
932   track.fFlags = kTPCin;
933   track.fID    = fID;
934
935   track.fFlags |= fFlags & kTPCpid; //copy the TPCpid status flag
936  
937   for (Int_t i=0;i<3;i++) track.fKinkIndexes[i] = fKinkIndexes[i];
938   
939   return kTRUE;
940     
941 }
942
943 //_______________________________________________________________________
944 void AliESDtrack::MakeMiniESDtrack(){
945   // Resets everything except
946   // fFlags: Reconstruction status flags 
947   // fLabel: Track label
948   // fID:  Unique ID of the track
949   // Impact parameter information
950   // fR[AliPID::kSPECIES]: combined "detector response probability"
951   // Running track parameters in the base class (AliExternalTrackParam)
952   
953   fTrackLength = 0;
954
955   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i] = 0;
956
957   // Reset track parameters constrained to the primary vertex
958   delete fCp;fCp = 0;
959
960   // Reset track parameters at the inner wall of TPC
961   delete fIp;fIp = 0;
962   delete fTPCInner;fTPCInner=0;
963   // Reset track parameters at the inner wall of the TRD
964   delete fOp;fOp = 0;
965
966
967   // Reset ITS track related information
968   fITSchi2 = 0;
969   fITSncls = 0;       
970   fITSClusterMap=0;
971   fITSsignal = 0;     
972   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=0; 
973   fITSLabel = 0;       
974
975   // Reset TPC related track information
976   fTPCchi2 = 0;       
977   fTPCncls = 0;       
978   fTPCnclsF = 0;       
979   fTPCClusterMap = 0;  
980   fTPCSharedMap = 0;  
981   fTPCsignal= 0;      
982   fTPCsignalS= 0;      
983   fTPCsignalN= 0;      
984   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0; 
985   fTPCLabel=0;       
986   for (Int_t i=0;i<4;i++) fTPCPoints[i] = 0;
987   for (Int_t i=0; i<3;i++)   fKinkIndexes[i] = 0;
988   for (Int_t i=0; i<3;i++)   fV0Indexes[i] = 0;
989
990   // Reset TRD related track information
991   fTRDchi2 = 0;        
992   fTRDncls = 0;       
993   fTRDncls0 = 0;       
994   fTRDsignal = 0;      
995   for (Int_t i=0;i<kTRDnPlanes;i++) {
996     fTRDTimBin[i]  = 0;
997   }
998   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0; 
999   fTRDLabel = 0;       
1000   fTRDQuality  = 0;
1001   fTRDntracklets = 0;
1002   if(fTRDnSlices)
1003     delete[] fTRDslices;
1004   fTRDslices=0x0;
1005   fTRDnSlices=0;
1006   fTRDBudget  = 0;
1007
1008   // Reset TOF related track information
1009   fTOFchi2 = 0;        
1010   fTOFindex = -1;       
1011   fTOFsignal = 0;      
1012   fTOFCalChannel = 0;
1013   fTOFsignalToT = 0;
1014   fTOFsignalRaw = 0;
1015   fTOFsignalDz = 0;
1016   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i] = 0;
1017   for (Int_t i=0;i<3;i++) fTOFLabel[i] = 0;
1018   for (Int_t i=0;i<10;i++) fTOFInfo[i] = 0;
1019
1020   // Reset HMPID related track information
1021   fHMPIDchi2 = 0;     
1022   fHMPIDqn = 0;     
1023   fHMPIDcluIdx = -1;     
1024   fHMPIDsignal = 0;     
1025   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i] = 0;
1026   fHMPIDtrkTheta = 0;     
1027   fHMPIDtrkPhi = 0;      
1028   fHMPIDtrkX = 0;     
1029   fHMPIDtrkY = 0;      
1030   fHMPIDmipX = 0;
1031   fHMPIDmipY = 0;
1032   fEMCALindex = kEMCALNoMatch;
1033
1034   // reset global track chi2
1035   fGlobalChi2 = 0;
1036
1037   fVertexID = -2; // an orphan track
1038
1039   delete fFriendTrack; fFriendTrack = 0;
1040
1041 //_______________________________________________________________________
1042 Double_t AliESDtrack::GetMass() const {
1043   // Returns the mass of the most probable particle type
1044   Float_t max=0.;
1045   Int_t k=-1;
1046   for (Int_t i=0; i<AliPID::kSPECIES; i++) {
1047     if (fR[i]>max) {k=i; max=fR[i];}
1048   }
1049   if (k==0) { // dE/dx "crossing points" in the TPC
1050      Double_t p=GetP();
1051      if ((p>0.38)&&(p<0.48))
1052         if (fR[0]<fR[3]*10.) return AliPID::ParticleMass(AliPID::kKaon);
1053      if ((p>0.75)&&(p<0.85))
1054         if (fR[0]<fR[4]*10.) return AliPID::ParticleMass(AliPID::kProton);
1055      return 0.00051;
1056   }
1057   if (k==1) return AliPID::ParticleMass(AliPID::kMuon); 
1058   if (k==2||k==-1) return AliPID::ParticleMass(AliPID::kPion);
1059   if (k==3) return AliPID::ParticleMass(AliPID::kKaon);
1060   if (k==4) return AliPID::ParticleMass(AliPID::kProton);
1061   AliWarning("Undefined mass !");
1062   return AliPID::ParticleMass(AliPID::kPion);
1063 }
1064
1065 //______________________________________________________________________________
1066 Double_t AliESDtrack::E() const
1067 {
1068   // Returns the energy of the particle given its assumed mass.
1069   // Assumes the pion mass if the particle can't be identified properly.
1070   
1071   Double_t m = M();
1072   Double_t p = P();
1073   return TMath::Sqrt(p*p + m*m);
1074 }
1075
1076 //______________________________________________________________________________
1077 Double_t AliESDtrack::Y() const
1078 {
1079   // Returns the rapidity of a particle given its assumed mass.
1080   // Assumes the pion mass if the particle can't be identified properly.
1081   
1082   Double_t e = E();
1083   Double_t pz = Pz();
1084   if (e != TMath::Abs(pz)) { // energy was not equal to pz
1085     return 0.5*TMath::Log((e+pz)/(e-pz));
1086   } else { // energy was equal to pz
1087     return -999.;
1088   }
1089 }
1090
1091 //_______________________________________________________________________
1092 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
1093   //
1094   // This function updates track's running parameters 
1095   //
1096   Int_t *index=0;
1097   Bool_t rc=kTRUE;
1098
1099   SetStatus(flags);
1100   fLabel=t->GetLabel();
1101
1102   if (t->IsStartedTimeIntegral()) {
1103     SetStatus(kTIME);
1104     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
1105     SetIntegratedLength(t->GetIntegratedLength());
1106   }
1107
1108   Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1109   
1110   switch (flags) {
1111     
1112   case kITSin: case kITSout: case kITSrefit:
1113     fITSClusterMap=0;
1114     fITSncls=t->GetNumberOfClusters();
1115     index=fFriendTrack->GetITSindices(); 
1116     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
1117         index[i]=t->GetClusterIndex(i);
1118         if (i<fITSncls) {
1119            Int_t l=(index[i] & 0xf0000000) >> 28;
1120            SETBIT(fITSClusterMap,l);                 
1121         }
1122     }
1123     fITSchi2=t->GetChi2();
1124     fITSsignal=t->GetPIDsignal();
1125     fITSLabel = t->GetLabel();
1126     // keep in fOp the parameters outside ITS for ITS stand-alone tracks 
1127     if (flags==kITSout) { 
1128       if (!fOp) fOp=new AliExternalTrackParam(*t);
1129       else 
1130         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1131     }      
1132     break;
1133     
1134   case kTPCin: case kTPCrefit:
1135     fTPCLabel = t->GetLabel();
1136     if (flags==kTPCin)  fTPCInner=new AliExternalTrackParam(*t);
1137     if (!fIp) fIp=new AliExternalTrackParam(*t);
1138     else 
1139       fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1140   case kTPCout:
1141     index=fFriendTrack->GetTPCindices(); 
1142     if (flags & kTPCout){
1143       if (!fOp) fOp=new AliExternalTrackParam(*t);
1144       else 
1145         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1146     }
1147     fTPCncls=t->GetNumberOfClusters();    
1148     fTPCchi2=t->GetChi2();
1149     
1150      {//prevrow must be declared in separate namespace, otherwise compiler cries:
1151       //"jump to case label crosses initialization of `Int_t prevrow'"
1152        Int_t prevrow = -1;
1153        //       for (Int_t i=0;i<fTPCncls;i++) 
1154        for (Int_t i=0;i<AliESDfriendTrack::kMaxTPCcluster;i++) 
1155         {
1156           index[i]=t->GetClusterIndex(i);
1157           Int_t idx = index[i];
1158
1159           if (idx<0) continue; 
1160
1161           // Piotr's Cluster Map for HBT  
1162           // ### please change accordingly if cluster array is changing 
1163           // to "New TPC Tracking" style (with gaps in array) 
1164           Int_t sect = (idx&0xff000000)>>24;
1165           Int_t row = (idx&0x00ff0000)>>16;
1166           if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
1167
1168           fTPCClusterMap.SetBitNumber(row,kTRUE);
1169
1170           //Fill the gap between previous row and this row with 0 bits
1171           //In case  ###  pleas change it as well - just set bit 0 in case there 
1172           //is no associated clusters for current "i"
1173           if (prevrow < 0) 
1174            {
1175              prevrow = row;//if previous bit was not assigned yet == this is the first one
1176            }
1177           else
1178            { //we don't know the order (inner to outer or reverse)
1179              //just to be save in case it is going to change
1180              Int_t n = 0, m = 0;
1181              if (prevrow < row)
1182               {
1183                 n = prevrow;
1184                 m = row;
1185               }
1186              else
1187               {
1188                 n = row;
1189                 m = prevrow;
1190               }
1191
1192              for (Int_t j = n+1; j < m; j++)
1193               {
1194                 fTPCClusterMap.SetBitNumber(j,kFALSE);
1195               }
1196              prevrow = row; 
1197            }
1198           // End Of Piotr's Cluster Map for HBT
1199         }
1200      }
1201     fTPCsignal=t->GetPIDsignal();
1202     break;
1203
1204   case kTRDin: case kTRDrefit:
1205     break;
1206   case kTRDout:
1207     index     = fFriendTrack->GetTRDindices();
1208     fTRDLabel = t->GetLabel(); 
1209     fTRDchi2  = t->GetChi2();
1210     fTRDncls  = t->GetNumberOfClusters();
1211     for (Int_t i=0;i<6;i++) index[i]=t->GetTrackletIndex(i);
1212     
1213     fTRDsignal=t->GetPIDsignal();
1214     break;
1215   case kTRDbackup:
1216     if (!fOp) fOp=new AliExternalTrackParam(*t);
1217     else 
1218       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1219     fTRDncls0 = t->GetNumberOfClusters(); 
1220     break;
1221   case kTOFin: 
1222     break;
1223   case kTOFout: 
1224     break;
1225   case kTRDStop:
1226     break;
1227   default: 
1228     AliError("Wrong flag !");
1229     return kFALSE;
1230   }
1231
1232   return rc;
1233 }
1234
1235 //_______________________________________________________________________
1236 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
1237   //---------------------------------------------------------------------
1238   // This function returns external representation of the track parameters
1239   //---------------------------------------------------------------------
1240   x=GetX();
1241   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
1242 }
1243
1244 //_______________________________________________________________________
1245 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
1246   //---------------------------------------------------------------------
1247   // This function returns external representation of the cov. matrix
1248   //---------------------------------------------------------------------
1249   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
1250 }
1251
1252 //_______________________________________________________________________
1253 Bool_t AliESDtrack::GetConstrainedExternalParameters
1254                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1255   //---------------------------------------------------------------------
1256   // This function returns the constrained external track parameters
1257   //---------------------------------------------------------------------
1258   if (!fCp) return kFALSE;
1259   alpha=fCp->GetAlpha();
1260   x=fCp->GetX();
1261   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
1262   return kTRUE;
1263 }
1264
1265 //_______________________________________________________________________
1266 Bool_t 
1267 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
1268   //---------------------------------------------------------------------
1269   // This function returns the constrained external cov. matrix
1270   //---------------------------------------------------------------------
1271   if (!fCp) return kFALSE;
1272   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
1273   return kTRUE;
1274 }
1275
1276 Bool_t
1277 AliESDtrack::GetInnerExternalParameters
1278                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1279   //---------------------------------------------------------------------
1280   // This function returns external representation of the track parameters 
1281   // at the inner layer of TPC
1282   //---------------------------------------------------------------------
1283   if (!fIp) return kFALSE;
1284   alpha=fIp->GetAlpha();
1285   x=fIp->GetX();
1286   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
1287   return kTRUE;
1288 }
1289
1290 Bool_t 
1291 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
1292  //---------------------------------------------------------------------
1293  // This function returns external representation of the cov. matrix 
1294  // at the inner layer of TPC
1295  //---------------------------------------------------------------------
1296   if (!fIp) return kFALSE;
1297   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
1298   return kTRUE;
1299 }
1300
1301 void 
1302 AliESDtrack::SetOuterParam(const AliExternalTrackParam *p, ULong_t flags) {
1303   //
1304   // This is a direct setter for the outer track parameters
1305   //
1306   SetStatus(flags);
1307   if (fOp) delete fOp;
1308   fOp=new AliExternalTrackParam(*p);
1309 }
1310
1311 Bool_t 
1312 AliESDtrack::GetOuterExternalParameters
1313                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1314   //---------------------------------------------------------------------
1315   // This function returns external representation of the track parameters 
1316   // at the inner layer of TRD
1317   //---------------------------------------------------------------------
1318   if (!fOp) return kFALSE;
1319   alpha=fOp->GetAlpha();
1320   x=fOp->GetX();
1321   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
1322   return kTRUE;
1323 }
1324
1325 Bool_t 
1326 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
1327  //---------------------------------------------------------------------
1328  // This function returns external representation of the cov. matrix 
1329  // at the inner layer of TRD
1330  //---------------------------------------------------------------------
1331   if (!fOp) return kFALSE;
1332   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
1333   return kTRUE;
1334 }
1335
1336 Int_t AliESDtrack::GetNcls(Int_t idet) const
1337 {
1338   // Get number of clusters by subdetector index
1339   //
1340   Int_t ncls = 0;
1341   switch(idet){
1342   case 0:
1343     ncls = fITSncls;
1344     break;
1345   case 1:
1346     ncls = fTPCncls;
1347     break;
1348   case 2:
1349     ncls = fTRDncls;
1350     break;
1351   case 3:
1352     if (fTOFindex != -1)
1353       ncls = 1;
1354     break;
1355   case 4: //PHOS
1356     break;
1357   case 5: //HMPID
1358     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1359       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1360         ncls = 1;
1361       }
1362     }    
1363     break;
1364   default:
1365     break;
1366   }
1367   return ncls;
1368 }
1369
1370 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
1371 {
1372   // Get cluster index array by subdetector index
1373   //
1374   Int_t ncls = 0;
1375   switch(idet){
1376   case 0:
1377     ncls = GetITSclusters(idx);
1378     break;
1379   case 1:
1380     ncls = GetTPCclusters(idx);
1381     break;
1382   case 2:
1383     ncls = GetTRDclusters(idx);
1384     break;
1385   case 3:
1386     if (fTOFindex != -1) {
1387       idx[0] = fTOFindex;
1388       ncls = 1;
1389     }
1390     break;
1391   case 4: //PHOS
1392     break;
1393   case 5:
1394     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1395       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1396         idx[0] = GetHMPIDcluIdx();
1397         ncls = 1;
1398       }
1399     }    
1400     break;
1401   case 6: //EMCAL
1402     break;
1403   default:
1404     break;
1405   }
1406   return ncls;
1407 }
1408
1409 //_______________________________________________________________________
1410 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
1411   // Returns the array with integrated times for each particle hypothesis
1412   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
1413 }
1414
1415 //_______________________________________________________________________
1416 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
1417   // Sets the array with integrated times for each particle hypotesis
1418   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
1419 }
1420
1421 //_______________________________________________________________________
1422 void AliESDtrack::SetITSpid(const Double_t *p) {
1423   // Sets values for the probability of each particle type (in ITS)
1424   SetPIDValues(fITSr,p,AliPID::kSPECIES);
1425   SetStatus(AliESDtrack::kITSpid);
1426 }
1427
1428 //_______________________________________________________________________
1429 void AliESDtrack::GetITSpid(Double_t *p) const {
1430   // Gets the probability of each particle type (in ITS)
1431   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
1432 }
1433
1434 //_______________________________________________________________________
1435 Char_t AliESDtrack::GetITSclusters(Int_t *idx) const {
1436   //---------------------------------------------------------------------
1437   // This function returns indices of the assgined ITS clusters 
1438   //---------------------------------------------------------------------
1439   if (idx!=0) {
1440      Int_t *index=fFriendTrack->GetITSindices();
1441      for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) {
1442          if ( (i>=fITSncls) && (i<6) ) idx[i]=-1;
1443          else idx[i]=index[i];
1444      }
1445   }
1446   return fITSncls;
1447 }
1448
1449 //_______________________________________________________________________
1450 Bool_t AliESDtrack::GetITSModuleIndexInfo(Int_t ilayer,Int_t &idet,Int_t &status,
1451                                          Float_t &xloc,Float_t &zloc) const {
1452   //----------------------------------------------------------------------
1453   // This function encodes in the module number also the status of cluster association
1454   // "status" can have the following values: 
1455   // 1 "found" (cluster is associated), 
1456   // 2 "dead" (module is dead from OCDB), 
1457   // 3 "skipped" (module or layer forced to be skipped),
1458   // 4 "outinz" (track out of z acceptance), 
1459   // 5 "nocls" (no clusters in the road), 
1460   // 6 "norefit" (cluster rejected during refit), 
1461   // 7 "deadzspd" (holes in z in SPD)
1462   // Also given are the coordinates of the crossing point of track and module
1463   // (in the local module ref. system)
1464   // WARNING: THIS METHOD HAS TO BE SYNCHRONIZED WITH AliITStrackV2::GetModuleIndexInfo()!
1465   //----------------------------------------------------------------------
1466
1467   if(fITSModule[ilayer]==-1) {
1468     idet = -1;
1469     status=0;
1470     xloc=-99.; zloc=-99.;
1471     return kFALSE;
1472   }
1473
1474   Int_t module = fITSModule[ilayer];
1475
1476   idet = Int_t(module/1000000);
1477
1478   module -= idet*1000000;
1479
1480   status = Int_t(module/100000);
1481
1482   module -= status*100000;
1483
1484   Int_t signs = Int_t(module/10000);
1485
1486   module-=signs*10000;
1487
1488   Int_t xInt = Int_t(module/100);
1489   module -= xInt*100;
1490
1491   Int_t zInt = module;
1492
1493   if(signs==1) { xInt*=1; zInt*=1; }
1494   if(signs==2) { xInt*=1; zInt*=-1; }
1495   if(signs==3) { xInt*=-1; zInt*=1; }
1496   if(signs==4) { xInt*=-1; zInt*=-1; }
1497
1498   xloc = 0.1*(Float_t)xInt;
1499   zloc = 0.1*(Float_t)zInt;
1500
1501   if(status==4) idet = -1;
1502
1503   return kTRUE;
1504 }
1505
1506 //_______________________________________________________________________
1507 UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
1508   //---------------------------------------------------------------------
1509   // This function returns indices of the assgined ITS clusters 
1510   //---------------------------------------------------------------------
1511   if (idx!=0) {
1512     Int_t *index=fFriendTrack->GetTPCindices();
1513     for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
1514   }
1515   return fTPCncls;
1516 }
1517
1518 Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
1519   //
1520   // GetDensity of the clusters on given region between row0 and row1
1521   // Dead zone effect takin into acoount
1522   //
1523   Int_t good  = 0;
1524   Int_t found = 0;
1525   //  
1526   Int_t *index=fFriendTrack->GetTPCindices();
1527   for (Int_t i=row0;i<=row1;i++){     
1528     Int_t idx = index[i];
1529     if (idx!=-1)  good++;             // track outside of dead zone
1530     if (idx>0)    found++;
1531   }
1532   Float_t density=0.5;
1533   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
1534   return density;
1535 }
1536
1537 //_______________________________________________________________________
1538 void AliESDtrack::SetTPCpid(const Double_t *p) {  
1539   // Sets values for the probability of each particle type (in TPC)
1540   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
1541   SetStatus(AliESDtrack::kTPCpid);
1542 }
1543
1544 //_______________________________________________________________________
1545 void AliESDtrack::GetTPCpid(Double_t *p) const {
1546   // Gets the probability of each particle type (in TPC)
1547   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
1548 }
1549
1550 //_______________________________________________________________________
1551 UChar_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
1552   //---------------------------------------------------------------------
1553   // This function returns indices of the assgined TRD clusters 
1554   //---------------------------------------------------------------------
1555   if (idx!=0) {
1556      Int_t *index=fFriendTrack->GetTRDindices();
1557      for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
1558   }
1559   return fTRDncls;
1560 }
1561
1562 //_______________________________________________________________________
1563 UChar_t AliESDtrack::GetTRDtracklets(Int_t *idx) const {
1564   //---------------------------------------------------------------------
1565   // This function returns indices of the assigned TRD tracklets 
1566   //---------------------------------------------------------------------
1567   if (idx!=0) {
1568      Int_t *index=fFriendTrack->GetTRDindices();
1569      for (Int_t i=0; i<6/*AliESDfriendTrack::kMaxTRDcluster*/; i++) idx[i]=index[i];
1570   }
1571   return fTRDncls;
1572 }
1573
1574 //_______________________________________________________________________
1575 void AliESDtrack::SetTRDpid(const Double_t *p) {  
1576   // Sets values for the probability of each particle type (in TRD)
1577   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
1578   SetStatus(AliESDtrack::kTRDpid);
1579 }
1580
1581 //_______________________________________________________________________
1582 void AliESDtrack::GetTRDpid(Double_t *p) const {
1583   // Gets the probability of each particle type (in TRD)
1584   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
1585 }
1586
1587 //_______________________________________________________________________
1588 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
1589 {
1590   // Sets the probability of particle type iSpecies to p (in TRD)
1591   fTRDr[iSpecies] = p;
1592 }
1593
1594 Double_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
1595 {
1596   // Returns the probability of particle type iSpecies (in TRD)
1597   return fTRDr[iSpecies];
1598 }
1599
1600 void  AliESDtrack::SetNumberOfTRDslices(Int_t n) {
1601   //Sets the number of slices used for PID 
1602   if (fTRDnSlices != 0) return;
1603   fTRDnSlices=kTRDnPlanes*n;
1604   fTRDslices=new Double32_t[fTRDnSlices];
1605   for (Int_t i=0; i<fTRDnSlices; i++) fTRDslices[i]=-1.;
1606 }
1607
1608 void  AliESDtrack::SetTRDslice(Double_t q, Int_t plane, Int_t slice) {
1609   //Sets the charge q in the slice of the plane
1610   Int_t ns=GetNumberOfTRDslices();
1611   if (ns==0) {
1612     AliError("No TRD slices allocated for this track !");
1613     return;
1614   }
1615
1616   if ((plane<0) || (plane>=kTRDnPlanes)) {
1617     AliError("Wrong TRD plane !");
1618     return;
1619   }
1620   if ((slice<0) || (slice>=ns)) {
1621     AliError("Wrong TRD slice !");
1622     return;
1623   }
1624   Int_t n=plane*ns + slice;
1625   fTRDslices[n]=q;
1626 }
1627
1628 Double_t  AliESDtrack::GetTRDslice(Int_t plane, Int_t slice) const {
1629   //Gets the charge from the slice of the plane
1630   Int_t ns=GetNumberOfTRDslices();
1631   if (ns==0) {
1632     //AliError("No TRD slices allocated for this track !");
1633     return -1.;
1634   }
1635
1636   if ((plane<0) || (plane>=kTRDnPlanes)) {
1637     AliError("Wrong TRD plane !");
1638     return -1.;
1639   }
1640   if ((slice<-1) || (slice>=ns)) {
1641     //AliError("Wrong TRD slice !");  
1642     return -1.;
1643   }
1644
1645   if (slice==-1) {
1646     Double_t q=0.;
1647     for (Int_t i=0; i<ns; i++) q+=fTRDslices[plane*ns + i];
1648     return q/ns;
1649   }
1650
1651   return fTRDslices[plane*ns + slice];
1652 }
1653
1654
1655 //_______________________________________________________________________
1656 void AliESDtrack::SetTOFpid(const Double_t *p) {  
1657   // Sets the probability of each particle type (in TOF)
1658   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
1659   SetStatus(AliESDtrack::kTOFpid);
1660 }
1661
1662 //_______________________________________________________________________
1663 void AliESDtrack::SetTOFLabel(const Int_t *p) {  
1664   // Sets  (in TOF)
1665   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
1666 }
1667
1668 //_______________________________________________________________________
1669 void AliESDtrack::GetTOFpid(Double_t *p) const {
1670   // Gets probabilities of each particle type (in TOF)
1671   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
1672 }
1673
1674 //_______________________________________________________________________
1675 void AliESDtrack::GetTOFLabel(Int_t *p) const {
1676   // Gets (in TOF)
1677   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
1678 }
1679
1680 //_______________________________________________________________________
1681 void AliESDtrack::GetTOFInfo(Float_t *info) const {
1682   // Gets (in TOF)
1683   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
1684 }
1685
1686 //_______________________________________________________________________
1687 void AliESDtrack::SetTOFInfo(Float_t*info) {
1688   // Gets (in TOF)
1689   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
1690 }
1691
1692
1693
1694 //_______________________________________________________________________
1695 void AliESDtrack::SetHMPIDpid(const Double_t *p) {  
1696   // Sets the probability of each particle type (in HMPID)
1697   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
1698   SetStatus(AliESDtrack::kHMPIDpid);
1699 }
1700
1701 //_______________________________________________________________________
1702 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
1703   // Gets probabilities of each particle type (in HMPID)
1704   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
1705 }
1706
1707
1708
1709 //_______________________________________________________________________
1710 void AliESDtrack::SetESDpid(const Double_t *p) {  
1711   // Sets the probability of each particle type for the ESD track
1712   SetPIDValues(fR,p,AliPID::kSPECIES);
1713   SetStatus(AliESDtrack::kESDpid);
1714 }
1715
1716 //_______________________________________________________________________
1717 void AliESDtrack::GetESDpid(Double_t *p) const {
1718   // Gets probability of each particle type for the ESD track
1719   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
1720 }
1721
1722 //_______________________________________________________________________
1723 Bool_t AliESDtrack::RelateToVertexTPC(const AliESDVertex *vtx, 
1724 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
1725   //
1726   // Try to relate the TPC-only track parameters to the vertex "vtx", 
1727   // if the (rough) transverse impact parameter is not bigger then "maxd". 
1728   //            Magnetic field is "b" (kG).
1729   //
1730   // a) The TPC-only paramters are extapolated to the DCA to the vertex.
1731   // b) The impact parameters and their covariance matrix are calculated.
1732   // c) An attempt to constrain the TPC-only params to the vertex is done.
1733   //    The constrained params are returned via "cParam".
1734   //
1735   // In the case of success, the returned value is kTRUE
1736   // otherwise, it's kFALSE)
1737   // 
1738
1739   if (!fTPCInner) return kFALSE;
1740   if (!vtx) return kFALSE;
1741
1742   Double_t dz[2],cov[3];
1743   if (!fTPCInner->PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
1744
1745   fdTPC = dz[0];
1746   fzTPC = dz[1];  
1747   fCddTPC = cov[0];
1748   fCdzTPC = cov[1];
1749   fCzzTPC = cov[2];
1750   
1751   Double_t covar[6]; vtx->GetCovMatrix(covar);
1752   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
1753   Double_t c[3]={covar[2],0.,covar[5]};
1754
1755   Double_t chi2=GetPredictedChi2(p,c);
1756   if (chi2>kVeryBig) return kFALSE;
1757
1758   fCchi2TPC=chi2;
1759
1760   if (!cParam) return kTRUE;
1761
1762   *cParam = *fTPCInner;
1763   if (!cParam->Update(p,c)) return kFALSE;
1764
1765   return kTRUE;
1766 }
1767
1768 //_______________________________________________________________________
1769 Bool_t AliESDtrack::RelateToVertex(const AliESDVertex *vtx, 
1770 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
1771   //
1772   // Try to relate this track to the vertex "vtx", 
1773   // if the (rough) transverse impact parameter is not bigger then "maxd". 
1774   //            Magnetic field is "b" (kG).
1775   //
1776   // a) The track gets extapolated to the DCA to the vertex.
1777   // b) The impact parameters and their covariance matrix are calculated.
1778   // c) An attempt to constrain this track to the vertex is done.
1779   //    The constrained params are returned via "cParam".
1780   //
1781   // In the case of success, the returned value is kTRUE
1782   // (otherwise, it's kFALSE)
1783   //  
1784
1785   if (!vtx) return kFALSE;
1786
1787   Double_t dz[2],cov[3];
1788   if (!PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
1789
1790   fD = dz[0];
1791   fZ = dz[1];  
1792   fCdd = cov[0];
1793   fCdz = cov[1];
1794   fCzz = cov[2];
1795   
1796   Double_t covar[6]; vtx->GetCovMatrix(covar);
1797   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
1798   Double_t c[3]={covar[2],0.,covar[5]};
1799
1800   Double_t chi2=GetPredictedChi2(p,c);
1801   if (chi2>kVeryBig) return kFALSE;
1802
1803   fCchi2=chi2;
1804
1805
1806   //--- Could now these lines be removed ? ---
1807   delete fCp;
1808   fCp=new AliExternalTrackParam(*this);  
1809
1810   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
1811   //----------------------------------------
1812
1813   fVertexID = vtx->GetID();
1814
1815   if (!cParam) return kTRUE;
1816
1817   *cParam = *this;
1818   if (!cParam->Update(p,c)) return kFALSE; 
1819
1820   return kTRUE;
1821 }
1822
1823 //_______________________________________________________________________
1824 void AliESDtrack::Print(Option_t *) const {
1825   // Prints info on the track
1826   AliExternalTrackParam::Print();
1827   printf("ESD track info\n") ; 
1828   Double_t p[AliPID::kSPECIESN] ; 
1829   Int_t index = 0 ; 
1830   if( IsOn(kITSpid) ){
1831     printf("From ITS: ") ; 
1832     GetITSpid(p) ; 
1833     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1834       printf("%f, ", p[index]) ;
1835     printf("\n           signal = %f\n", GetITSsignal()) ;
1836   } 
1837   if( IsOn(kTPCpid) ){
1838     printf("From TPC: ") ; 
1839     GetTPCpid(p) ; 
1840     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1841       printf("%f, ", p[index]) ;
1842     printf("\n           signal = %f\n", GetTPCsignal()) ;
1843   }
1844   if( IsOn(kTRDpid) ){
1845     printf("From TRD: ") ; 
1846     GetTRDpid(p) ; 
1847     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1848       printf("%f, ", p[index]) ;
1849       printf("\n           signal = %f\n", GetTRDsignal()) ;
1850   }
1851   if( IsOn(kTOFpid) ){
1852     printf("From TOF: ") ; 
1853     GetTOFpid(p) ; 
1854     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1855       printf("%f, ", p[index]) ;
1856     printf("\n           signal = %f\n", GetTOFsignal()) ;
1857   }
1858   if( IsOn(kHMPIDpid) ){
1859     printf("From HMPID: ") ; 
1860     GetHMPIDpid(p) ; 
1861     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1862       printf("%f, ", p[index]) ;
1863     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
1864   }
1865
1866
1867
1868 //
1869 // Draw functionality
1870 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
1871 //
1872 void AliESDtrack::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
1873   //
1874   // Fill points in the polymarker
1875   //
1876   TObjArray arrayRef;
1877   arrayRef.AddLast(new AliExternalTrackParam(*this));
1878   if (fIp) arrayRef.AddLast(new AliExternalTrackParam(*fIp));
1879   if (fOp) arrayRef.AddLast(new AliExternalTrackParam(*fOp));
1880   //
1881   Double_t mpos[3]={0,0,0};
1882   Int_t entries=arrayRef.GetEntries();
1883   for (Int_t i=0;i<entries;i++){
1884     Double_t pos[3];
1885     ((AliExternalTrackParam*)arrayRef.At(i))->GetXYZ(pos);
1886     mpos[0]+=pos[0]/entries;
1887     mpos[1]+=pos[1]/entries;
1888     mpos[2]+=pos[2]/entries;    
1889   }
1890   // Rotate to the mean position
1891   //
1892   Float_t fi= TMath::ATan2(mpos[1],mpos[0]);
1893   for (Int_t i=0;i<entries;i++){
1894     Bool_t res = ((AliExternalTrackParam*)arrayRef.At(i))->Rotate(fi);
1895     if (!res) delete arrayRef.RemoveAt(i);
1896   }
1897   Int_t counter=0;
1898   for (Double_t r=minR; r<maxR; r+=stepR){
1899     Double_t sweight=0;
1900     Double_t mlpos[3]={0,0,0};
1901     for (Int_t i=0;i<entries;i++){
1902       Double_t point[3]={0,0,0};
1903       AliExternalTrackParam *param = ((AliExternalTrackParam*)arrayRef.At(i));
1904       if (!param) continue;
1905       if (param->GetXYZAt(r,magF,point)){
1906         Double_t weight = 1./(10.+(r-param->GetX())*(r-param->GetX()));
1907         sweight+=weight;
1908         mlpos[0]+=point[0]*weight;
1909         mlpos[1]+=point[1]*weight;
1910         mlpos[2]+=point[2]*weight;
1911       }
1912     }
1913     if (sweight>0){
1914       mlpos[0]/=sweight;
1915       mlpos[1]/=sweight;
1916       mlpos[2]/=sweight;      
1917       pol->SetPoint(counter,mlpos[0],mlpos[1], mlpos[2]);
1918       printf("xyz\t%f\t%f\t%f\n",mlpos[0], mlpos[1],mlpos[2]);
1919       counter++;
1920     }
1921   }
1922 }