]> git.uio.no Git - u/mrichter/AliRoot.git/blob - STEER/AliESDtrack.cxx
upgraded the version number
[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 kTRDout: case kTRDin: case kTRDrefit:
1199     index     = fFriendTrack->GetTRDindices();
1200     fTRDLabel = t->GetLabel(); 
1201     fTRDchi2  = t->GetChi2();
1202     fTRDncls  = t->GetNumberOfClusters();
1203     for (Int_t i=0;i<6;i++) index[i]=t->GetTrackletIndex(i);
1204     
1205     fTRDsignal=t->GetPIDsignal();
1206     break;
1207   case kTRDbackup:
1208     if (!fOp) fOp=new AliExternalTrackParam(*t);
1209     else 
1210       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1211     fTRDncls0 = t->GetNumberOfClusters(); 
1212     break;
1213   case kTOFin: 
1214     break;
1215   case kTOFout: 
1216     break;
1217   case kTRDStop:
1218     break;
1219   default: 
1220     AliError("Wrong flag !");
1221     return kFALSE;
1222   }
1223
1224   return rc;
1225 }
1226
1227 //_______________________________________________________________________
1228 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
1229   //---------------------------------------------------------------------
1230   // This function returns external representation of the track parameters
1231   //---------------------------------------------------------------------
1232   x=GetX();
1233   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
1234 }
1235
1236 //_______________________________________________________________________
1237 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
1238   //---------------------------------------------------------------------
1239   // This function returns external representation of the cov. matrix
1240   //---------------------------------------------------------------------
1241   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
1242 }
1243
1244 //_______________________________________________________________________
1245 Bool_t AliESDtrack::GetConstrainedExternalParameters
1246                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1247   //---------------------------------------------------------------------
1248   // This function returns the constrained external track parameters
1249   //---------------------------------------------------------------------
1250   if (!fCp) return kFALSE;
1251   alpha=fCp->GetAlpha();
1252   x=fCp->GetX();
1253   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
1254   return kTRUE;
1255 }
1256
1257 //_______________________________________________________________________
1258 Bool_t 
1259 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
1260   //---------------------------------------------------------------------
1261   // This function returns the constrained external cov. matrix
1262   //---------------------------------------------------------------------
1263   if (!fCp) return kFALSE;
1264   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
1265   return kTRUE;
1266 }
1267
1268 Bool_t
1269 AliESDtrack::GetInnerExternalParameters
1270                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1271   //---------------------------------------------------------------------
1272   // This function returns external representation of the track parameters 
1273   // at the inner layer of TPC
1274   //---------------------------------------------------------------------
1275   if (!fIp) return kFALSE;
1276   alpha=fIp->GetAlpha();
1277   x=fIp->GetX();
1278   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
1279   return kTRUE;
1280 }
1281
1282 Bool_t 
1283 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
1284  //---------------------------------------------------------------------
1285  // This function returns external representation of the cov. matrix 
1286  // at the inner layer of TPC
1287  //---------------------------------------------------------------------
1288   if (!fIp) return kFALSE;
1289   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
1290   return kTRUE;
1291 }
1292
1293 void 
1294 AliESDtrack::SetOuterParam(const AliExternalTrackParam *p, ULong_t flags) {
1295   //
1296   // This is a direct setter for the outer track parameters
1297   //
1298   SetStatus(flags);
1299   if (fOp) delete fOp;
1300   fOp=new AliExternalTrackParam(*p);
1301 }
1302
1303 Bool_t 
1304 AliESDtrack::GetOuterExternalParameters
1305                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1306   //---------------------------------------------------------------------
1307   // This function returns external representation of the track parameters 
1308   // at the inner layer of TRD
1309   //---------------------------------------------------------------------
1310   if (!fOp) return kFALSE;
1311   alpha=fOp->GetAlpha();
1312   x=fOp->GetX();
1313   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
1314   return kTRUE;
1315 }
1316
1317 Bool_t 
1318 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
1319  //---------------------------------------------------------------------
1320  // This function returns external representation of the cov. matrix 
1321  // at the inner layer of TRD
1322  //---------------------------------------------------------------------
1323   if (!fOp) return kFALSE;
1324   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
1325   return kTRUE;
1326 }
1327
1328 Int_t AliESDtrack::GetNcls(Int_t idet) const
1329 {
1330   // Get number of clusters by subdetector index
1331   //
1332   Int_t ncls = 0;
1333   switch(idet){
1334   case 0:
1335     ncls = fITSncls;
1336     break;
1337   case 1:
1338     ncls = fTPCncls;
1339     break;
1340   case 2:
1341     ncls = fTRDncls;
1342     break;
1343   case 3:
1344     if (fTOFindex != -1)
1345       ncls = 1;
1346     break;
1347   case 4: //PHOS
1348     break;
1349   case 5: //HMPID
1350     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1351       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1352         ncls = 1;
1353       }
1354     }    
1355     break;
1356   default:
1357     break;
1358   }
1359   return ncls;
1360 }
1361
1362 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
1363 {
1364   // Get cluster index array by subdetector index
1365   //
1366   Int_t ncls = 0;
1367   switch(idet){
1368   case 0:
1369     ncls = GetITSclusters(idx);
1370     break;
1371   case 1:
1372     ncls = GetTPCclusters(idx);
1373     break;
1374   case 2:
1375     ncls = GetTRDclusters(idx);
1376     break;
1377   case 3:
1378     if (fTOFindex != -1) {
1379       idx[0] = fTOFindex;
1380       ncls = 1;
1381     }
1382     break;
1383   case 4: //PHOS
1384     break;
1385   case 5:
1386     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1387       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1388         idx[0] = GetHMPIDcluIdx();
1389         ncls = 1;
1390       }
1391     }    
1392     break;
1393   case 6: //EMCAL
1394     break;
1395   default:
1396     break;
1397   }
1398   return ncls;
1399 }
1400
1401 //_______________________________________________________________________
1402 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
1403   // Returns the array with integrated times for each particle hypothesis
1404   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
1405 }
1406
1407 //_______________________________________________________________________
1408 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
1409   // Sets the array with integrated times for each particle hypotesis
1410   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
1411 }
1412
1413 //_______________________________________________________________________
1414 void AliESDtrack::SetITSpid(const Double_t *p) {
1415   // Sets values for the probability of each particle type (in ITS)
1416   SetPIDValues(fITSr,p,AliPID::kSPECIES);
1417   SetStatus(AliESDtrack::kITSpid);
1418 }
1419
1420 //_______________________________________________________________________
1421 void AliESDtrack::GetITSpid(Double_t *p) const {
1422   // Gets the probability of each particle type (in ITS)
1423   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
1424 }
1425
1426 //_______________________________________________________________________
1427 Char_t AliESDtrack::GetITSclusters(Int_t *idx) const {
1428   //---------------------------------------------------------------------
1429   // This function returns indices of the assgined ITS clusters 
1430   //---------------------------------------------------------------------
1431   if (idx!=0) {
1432      Int_t *index=fFriendTrack->GetITSindices();
1433      for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) {
1434          if ( (i>=fITSncls) && (i<6) ) idx[i]=-1;
1435          else idx[i]=index[i];
1436      }
1437   }
1438   return fITSncls;
1439 }
1440
1441 //_______________________________________________________________________
1442 Bool_t AliESDtrack::GetITSModuleIndexInfo(Int_t ilayer,Int_t &idet,Int_t &status,
1443                                          Float_t &xloc,Float_t &zloc) const {
1444   //----------------------------------------------------------------------
1445   // This function encodes in the module number also the status of cluster association
1446   // "status" can have the following values: 
1447   // 1 "found" (cluster is associated), 
1448   // 2 "dead" (module is dead from OCDB), 
1449   // 3 "skipped" (module or layer forced to be skipped),
1450   // 4 "outinz" (track out of z acceptance), 
1451   // 5 "nocls" (no clusters in the road), 
1452   // 6 "norefit" (cluster rejected during refit), 
1453   // 7 "deadzspd" (holes in z in SPD)
1454   // Also given are the coordinates of the crossing point of track and module
1455   // (in the local module ref. system)
1456   // WARNING: THIS METHOD HAS TO BE SYNCHRONIZED WITH AliITStrackV2::GetModuleIndexInfo()!
1457   //----------------------------------------------------------------------
1458
1459   if(fITSModule[ilayer]==-1) {
1460     idet = -1;
1461     status=0;
1462     xloc=-99.; zloc=-99.;
1463     return kFALSE;
1464   }
1465
1466   Int_t module = fITSModule[ilayer];
1467
1468   idet = Int_t(module/1000000);
1469
1470   module -= idet*1000000;
1471
1472   status = Int_t(module/100000);
1473
1474   module -= status*100000;
1475
1476   Int_t signs = Int_t(module/10000);
1477
1478   module-=signs*10000;
1479
1480   Int_t xInt = Int_t(module/100);
1481   module -= xInt*100;
1482
1483   Int_t zInt = module;
1484
1485   if(signs==1) { xInt*=1; zInt*=1; }
1486   if(signs==2) { xInt*=1; zInt*=-1; }
1487   if(signs==3) { xInt*=-1; zInt*=1; }
1488   if(signs==4) { xInt*=-1; zInt*=-1; }
1489
1490   xloc = 0.1*(Float_t)xInt;
1491   zloc = 0.1*(Float_t)zInt;
1492
1493   if(status==4) idet = -1;
1494
1495   return kTRUE;
1496 }
1497
1498 //_______________________________________________________________________
1499 UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
1500   //---------------------------------------------------------------------
1501   // This function returns indices of the assgined ITS clusters 
1502   //---------------------------------------------------------------------
1503   if (idx!=0) {
1504     Int_t *index=fFriendTrack->GetTPCindices();
1505     for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
1506   }
1507   return fTPCncls;
1508 }
1509
1510 Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
1511   //
1512   // GetDensity of the clusters on given region between row0 and row1
1513   // Dead zone effect takin into acoount
1514   //
1515   Int_t good  = 0;
1516   Int_t found = 0;
1517   //  
1518   Int_t *index=fFriendTrack->GetTPCindices();
1519   for (Int_t i=row0;i<=row1;i++){     
1520     Int_t idx = index[i];
1521     if (idx!=-1)  good++;             // track outside of dead zone
1522     if (idx>0)    found++;
1523   }
1524   Float_t density=0.5;
1525   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
1526   return density;
1527 }
1528
1529 //_______________________________________________________________________
1530 void AliESDtrack::SetTPCpid(const Double_t *p) {  
1531   // Sets values for the probability of each particle type (in TPC)
1532   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
1533   SetStatus(AliESDtrack::kTPCpid);
1534 }
1535
1536 //_______________________________________________________________________
1537 void AliESDtrack::GetTPCpid(Double_t *p) const {
1538   // Gets the probability of each particle type (in TPC)
1539   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
1540 }
1541
1542 //_______________________________________________________________________
1543 UChar_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
1544   //---------------------------------------------------------------------
1545   // This function returns indices of the assgined TRD clusters 
1546   //---------------------------------------------------------------------
1547   if (idx!=0) {
1548      Int_t *index=fFriendTrack->GetTRDindices();
1549      for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
1550   }
1551   return fTRDncls;
1552 }
1553
1554 //_______________________________________________________________________
1555 UChar_t AliESDtrack::GetTRDtracklets(Int_t *idx) const {
1556   //---------------------------------------------------------------------
1557   // This function returns indices of the assigned TRD tracklets 
1558   //---------------------------------------------------------------------
1559   if (idx!=0) {
1560      Int_t *index=fFriendTrack->GetTRDindices();
1561      for (Int_t i=0; i<6/*AliESDfriendTrack::kMaxTRDcluster*/; i++) idx[i]=index[i];
1562   }
1563   return fTRDncls;
1564 }
1565
1566 //_______________________________________________________________________
1567 void AliESDtrack::SetTRDpid(const Double_t *p) {  
1568   // Sets values for the probability of each particle type (in TRD)
1569   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
1570   SetStatus(AliESDtrack::kTRDpid);
1571 }
1572
1573 //_______________________________________________________________________
1574 void AliESDtrack::GetTRDpid(Double_t *p) const {
1575   // Gets the probability of each particle type (in TRD)
1576   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
1577 }
1578
1579 //_______________________________________________________________________
1580 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
1581 {
1582   // Sets the probability of particle type iSpecies to p (in TRD)
1583   fTRDr[iSpecies] = p;
1584 }
1585
1586 Double_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
1587 {
1588   // Returns the probability of particle type iSpecies (in TRD)
1589   return fTRDr[iSpecies];
1590 }
1591
1592 void  AliESDtrack::SetNumberOfTRDslices(Int_t n) {
1593   //Sets the number of slices used for PID 
1594   if (fTRDnSlices != 0) return;
1595   fTRDnSlices=kTRDnPlanes*n;
1596   fTRDslices=new Double32_t[fTRDnSlices];
1597   for (Int_t i=0; i<fTRDnSlices; i++) fTRDslices[i]=-1.;
1598 }
1599
1600 void  AliESDtrack::SetTRDslice(Double_t q, Int_t plane, Int_t slice) {
1601   //Sets the charge q in the slice of the plane
1602   Int_t ns=GetNumberOfTRDslices();
1603   if (ns==0) {
1604     AliError("No TRD slices allocated for this track !");
1605     return;
1606   }
1607
1608   if ((plane<0) || (plane>=kTRDnPlanes)) {
1609     AliError("Wrong TRD plane !");
1610     return;
1611   }
1612   if ((slice<0) || (slice>=ns)) {
1613     AliError("Wrong TRD slice !");
1614     return;
1615   }
1616   Int_t n=plane*ns + slice;
1617   fTRDslices[n]=q;
1618 }
1619
1620 Double_t  AliESDtrack::GetTRDslice(Int_t plane, Int_t slice) const {
1621   //Gets the charge from the slice of the plane
1622   Int_t ns=GetNumberOfTRDslices();
1623   if (ns==0) {
1624     //AliError("No TRD slices allocated for this track !");
1625     return -1.;
1626   }
1627
1628   if ((plane<0) || (plane>=kTRDnPlanes)) {
1629     AliError("Wrong TRD plane !");
1630     return -1.;
1631   }
1632   if ((slice<-1) || (slice>=ns)) {
1633     //AliError("Wrong TRD slice !");  
1634     return -1.;
1635   }
1636
1637   if (slice==-1) {
1638     Double_t q=0.;
1639     for (Int_t i=0; i<ns; i++) q+=fTRDslices[plane*ns + i];
1640     return q/ns;
1641   }
1642
1643   return fTRDslices[plane*ns + slice];
1644 }
1645
1646
1647 //_______________________________________________________________________
1648 void AliESDtrack::SetTOFpid(const Double_t *p) {  
1649   // Sets the probability of each particle type (in TOF)
1650   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
1651   SetStatus(AliESDtrack::kTOFpid);
1652 }
1653
1654 //_______________________________________________________________________
1655 void AliESDtrack::SetTOFLabel(const Int_t *p) {  
1656   // Sets  (in TOF)
1657   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
1658 }
1659
1660 //_______________________________________________________________________
1661 void AliESDtrack::GetTOFpid(Double_t *p) const {
1662   // Gets probabilities of each particle type (in TOF)
1663   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
1664 }
1665
1666 //_______________________________________________________________________
1667 void AliESDtrack::GetTOFLabel(Int_t *p) const {
1668   // Gets (in TOF)
1669   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
1670 }
1671
1672 //_______________________________________________________________________
1673 void AliESDtrack::GetTOFInfo(Float_t *info) const {
1674   // Gets (in TOF)
1675   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
1676 }
1677
1678 //_______________________________________________________________________
1679 void AliESDtrack::SetTOFInfo(Float_t*info) {
1680   // Gets (in TOF)
1681   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
1682 }
1683
1684
1685
1686 //_______________________________________________________________________
1687 void AliESDtrack::SetHMPIDpid(const Double_t *p) {  
1688   // Sets the probability of each particle type (in HMPID)
1689   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
1690   SetStatus(AliESDtrack::kHMPIDpid);
1691 }
1692
1693 //_______________________________________________________________________
1694 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
1695   // Gets probabilities of each particle type (in HMPID)
1696   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
1697 }
1698
1699
1700
1701 //_______________________________________________________________________
1702 void AliESDtrack::SetESDpid(const Double_t *p) {  
1703   // Sets the probability of each particle type for the ESD track
1704   SetPIDValues(fR,p,AliPID::kSPECIES);
1705   SetStatus(AliESDtrack::kESDpid);
1706 }
1707
1708 //_______________________________________________________________________
1709 void AliESDtrack::GetESDpid(Double_t *p) const {
1710   // Gets probability of each particle type for the ESD track
1711   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
1712 }
1713
1714 //_______________________________________________________________________
1715 Bool_t AliESDtrack::RelateToVertexTPC(const AliESDVertex *vtx, 
1716 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
1717   //
1718   // Try to relate the TPC-only track parameters to the vertex "vtx", 
1719   // if the (rough) transverse impact parameter is not bigger then "maxd". 
1720   //            Magnetic field is "b" (kG).
1721   //
1722   // a) The TPC-only paramters are extapolated to the DCA to the vertex.
1723   // b) The impact parameters and their covariance matrix are calculated.
1724   // c) An attempt to constrain the TPC-only params to the vertex is done.
1725   //    The constrained params are returned via "cParam".
1726   //
1727   // In the case of success, the returned value is kTRUE
1728   // otherwise, it's kFALSE)
1729   // 
1730
1731   if (!fTPCInner) return kFALSE;
1732   if (!vtx) return kFALSE;
1733
1734   Double_t dz[2],cov[3];
1735   if (!fTPCInner->PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
1736
1737   fdTPC = dz[0];
1738   fzTPC = dz[1];  
1739   fCddTPC = cov[0];
1740   fCdzTPC = cov[1];
1741   fCzzTPC = cov[2];
1742   
1743   Double_t covar[6]; vtx->GetCovMatrix(covar);
1744   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
1745   Double_t c[3]={covar[2],0.,covar[5]};
1746
1747   Double_t chi2=GetPredictedChi2(p,c);
1748   if (chi2>kVeryBig) return kFALSE;
1749
1750   fCchi2TPC=chi2;
1751
1752   if (!cParam) return kTRUE;
1753
1754   *cParam = *fTPCInner;
1755   if (!cParam->Update(p,c)) return kFALSE;
1756
1757   return kTRUE;
1758 }
1759
1760 //_______________________________________________________________________
1761 Bool_t AliESDtrack::RelateToVertex(const AliESDVertex *vtx, 
1762 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
1763   //
1764   // Try to relate this track to the vertex "vtx", 
1765   // if the (rough) transverse impact parameter is not bigger then "maxd". 
1766   //            Magnetic field is "b" (kG).
1767   //
1768   // a) The track gets extapolated to the DCA to the vertex.
1769   // b) The impact parameters and their covariance matrix are calculated.
1770   // c) An attempt to constrain this track to the vertex is done.
1771   //    The constrained params are returned via "cParam".
1772   //
1773   // In the case of success, the returned value is kTRUE
1774   // (otherwise, it's kFALSE)
1775   //  
1776
1777   if (!vtx) return kFALSE;
1778
1779   Double_t dz[2],cov[3];
1780   if (!PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
1781
1782   fD = dz[0];
1783   fZ = dz[1];  
1784   fCdd = cov[0];
1785   fCdz = cov[1];
1786   fCzz = cov[2];
1787   
1788   Double_t covar[6]; vtx->GetCovMatrix(covar);
1789   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
1790   Double_t c[3]={covar[2],0.,covar[5]};
1791
1792   Double_t chi2=GetPredictedChi2(p,c);
1793   if (chi2>kVeryBig) return kFALSE;
1794
1795   fCchi2=chi2;
1796
1797
1798   //--- Could now these lines be removed ? ---
1799   delete fCp;
1800   fCp=new AliExternalTrackParam(*this);  
1801
1802   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
1803   //----------------------------------------
1804
1805
1806   if (!cParam) return kTRUE;
1807
1808   *cParam = *this;
1809   if (!cParam->Update(p,c)) return kFALSE; 
1810
1811   return kTRUE;
1812 }
1813
1814 //_______________________________________________________________________
1815 void AliESDtrack::Print(Option_t *) const {
1816   // Prints info on the track
1817   AliExternalTrackParam::Print();
1818   printf("ESD track info\n") ; 
1819   Double_t p[AliPID::kSPECIESN] ; 
1820   Int_t index = 0 ; 
1821   if( IsOn(kITSpid) ){
1822     printf("From ITS: ") ; 
1823     GetITSpid(p) ; 
1824     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1825       printf("%f, ", p[index]) ;
1826     printf("\n           signal = %f\n", GetITSsignal()) ;
1827   } 
1828   if( IsOn(kTPCpid) ){
1829     printf("From TPC: ") ; 
1830     GetTPCpid(p) ; 
1831     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1832       printf("%f, ", p[index]) ;
1833     printf("\n           signal = %f\n", GetTPCsignal()) ;
1834   }
1835   if( IsOn(kTRDpid) ){
1836     printf("From TRD: ") ; 
1837     GetTRDpid(p) ; 
1838     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1839       printf("%f, ", p[index]) ;
1840       printf("\n           signal = %f\n", GetTRDsignal()) ;
1841   }
1842   if( IsOn(kTOFpid) ){
1843     printf("From TOF: ") ; 
1844     GetTOFpid(p) ; 
1845     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1846       printf("%f, ", p[index]) ;
1847     printf("\n           signal = %f\n", GetTOFsignal()) ;
1848   }
1849   if( IsOn(kHMPIDpid) ){
1850     printf("From HMPID: ") ; 
1851     GetHMPIDpid(p) ; 
1852     for(index = 0 ; index < AliPID::kSPECIES; index++) 
1853       printf("%f, ", p[index]) ;
1854     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
1855   }
1856
1857
1858
1859 //
1860 // Draw functionality
1861 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
1862 //
1863 void AliESDtrack::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
1864   //
1865   // Fill points in the polymarker
1866   //
1867   TObjArray arrayRef;
1868   arrayRef.AddLast(new AliExternalTrackParam(*this));
1869   if (fIp) arrayRef.AddLast(new AliExternalTrackParam(*fIp));
1870   if (fOp) arrayRef.AddLast(new AliExternalTrackParam(*fOp));
1871   //
1872   Double_t mpos[3]={0,0,0};
1873   Int_t entries=arrayRef.GetEntries();
1874   for (Int_t i=0;i<entries;i++){
1875     Double_t pos[3];
1876     ((AliExternalTrackParam*)arrayRef.At(i))->GetXYZ(pos);
1877     mpos[0]+=pos[0]/entries;
1878     mpos[1]+=pos[1]/entries;
1879     mpos[2]+=pos[2]/entries;    
1880   }
1881   // Rotate to the mean position
1882   //
1883   Float_t fi= TMath::ATan2(mpos[1],mpos[0]);
1884   for (Int_t i=0;i<entries;i++){
1885     Bool_t res = ((AliExternalTrackParam*)arrayRef.At(i))->Rotate(fi);
1886     if (!res) delete arrayRef.RemoveAt(i);
1887   }
1888   Int_t counter=0;
1889   for (Double_t r=minR; r<maxR; r+=stepR){
1890     Double_t sweight=0;
1891     Double_t mlpos[3]={0,0,0};
1892     for (Int_t i=0;i<entries;i++){
1893       Double_t point[3]={0,0,0};
1894       AliExternalTrackParam *param = ((AliExternalTrackParam*)arrayRef.At(i));
1895       if (!param) continue;
1896       if (param->GetXYZAt(r,magF,point)){
1897         Double_t weight = 1./(10.+(r-param->GetX())*(r-param->GetX()));
1898         sweight+=weight;
1899         mlpos[0]+=point[0]*weight;
1900         mlpos[1]+=point[1]*weight;
1901         mlpos[2]+=point[2]*weight;
1902       }
1903     }
1904     if (sweight>0){
1905       mlpos[0]/=sweight;
1906       mlpos[1]/=sweight;
1907       mlpos[2]/=sweight;      
1908       pol->SetPoint(counter,mlpos[0],mlpos[1], mlpos[2]);
1909       printf("xyz\t%f\t%f\t%f\n",mlpos[0], mlpos[1],mlpos[2]);
1910       counter++;
1911     }
1912   }
1913 }