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