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