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
 914   for (Int_t i=0;i<3;i++) track.fKinkIndexes[i] = fKinkIndexes[i];
 915
 916   return kTRUE;
 917
 918 }
 919
 920 //_______________________________________________________________________
 921 void AliESDtrack::MakeMiniESDtrack(){
 922   // Resets everything except
 923   // fFlags: Reconstruction status flags
 924   // fLabel: Track label
 925   // fID:  Unique ID of the track
 926   // Impact parameter information
 927   // fR[AliPID::kSPECIES]: combined "detector response probability"
 928   // Running track parameters in the base class (AliExternalTrackParam)
 929
 930   fTrackLength = 0;
 931
 932   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i] = 0;
 933
 934   // Reset track parameters constrained to the primary vertex
 935   delete fCp;fCp = 0;
 936
 937   // Reset track parameters at the inner wall of TPC
 938   delete fIp;fIp = 0;
 939   delete fTPCInner;fTPCInner=0;
 940   // Reset track parameters at the inner wall of the TRD
 941   delete fOp;fOp = 0;
 942
 943
 944   // Reset ITS track related information
 945   fITSchi2 = 0;
 946   fITSncls = 0;
 947   fITSClusterMap=0;
 948   fITSsignal = 0;
 949   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=0;
 950   fITSLabel = 0;
 951
 952   // Reset TPC related track information
 953   fTPCchi2 = 0;
 954   fTPCncls = 0;
 955   fTPCnclsF = 0;
 956   fTPCClusterMap = 0;
 957   fTPCSharedMap = 0;
 958   fTPCsignal= 0;
 959   fTPCsignalS= 0;
 960   fTPCsignalN= 0;
 961   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0;
 962   fTPCLabel=0;
 963   for (Int_t i=0;i<4;i++) fTPCPoints[i] = 0;
 964   for (Int_t i=0; i<3;i++)   fKinkIndexes[i] = 0;
 965   for (Int_t i=0; i<3;i++)   fV0Indexes[i] = 0;
 966
 967   // Reset TRD related track information
 968   fTRDchi2 = 0;
 969   fTRDncls = 0;
 970   fTRDncls0 = 0;
 971   fTRDsignal = 0;
 972   for (Int_t i=0;i<kTRDnPlanes;i++) {
 973     fTRDTimBin[i]  = 0;
 974   }
 975   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0;
 976   fTRDLabel = 0;
 977   fTRDQuality  = 0;
 978   fTRDpidQuality = 0;
 979   if(fTRDnSlices)
 980     delete[] fTRDslices;
 981   fTRDslices=0x0;
 982   fTRDnSlices=0;
 983   fTRDBudget  = 0;
 984
 985   // Reset TOF related track information
 986   fTOFchi2 = 0;
 987   fTOFindex = -1;
 988   fTOFsignal = 0;
 989   fTOFCalChannel = 0;
 990   fTOFsignalToT = 0;
 991   fTOFsignalRaw = 0;
 992   fTOFsignalDz = 0;
 993   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i] = 0;
 994   for (Int_t i=0;i<3;i++) fTOFLabel[i] = 0;
 995   for (Int_t i=0;i<10;i++) fTOFInfo[i] = 0;
 996
 997   // Reset HMPID related track information
 998   fHMPIDchi2 = 0;
 999   fHMPIDqn = 0;
1000   fHMPIDcluIdx = -1;
1001   fHMPIDsignal = 0;
1002   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i] = 0;
1003   fHMPIDtrkTheta = 0;
1004   fHMPIDtrkPhi = 0;
1005   fHMPIDtrkX = 0;
1006   fHMPIDtrkY = 0;
1007   fHMPIDmipX = 0;
1008   fHMPIDmipY = 0;
1009   fEMCALindex = kEMCALNoMatch;
1010
1011   delete fFriendTrack; fFriendTrack = 0;
1012 }
1013 //_______________________________________________________________________
1014 Double_t AliESDtrack::GetMass() const {
1015   // Returns the mass of the most probable particle type
1016   Float_t max=0.;
1017   Int_t k=-1;
1018   for (Int_t i=0; i<AliPID::kSPECIES; i++) {
1019     if (fR[i]>max) {k=i; max=fR[i];}
1020   }
1021   if (k==0) { // dE/dx "crossing points" in the TPC
1022      Double_t p=GetP();
1023      if ((p>0.38)&&(p<0.48))
1024         if (fR[0]<fR[3]*10.) return AliPID::ParticleMass(AliPID::kKaon);
1025      if ((p>0.75)&&(p<0.85))
1026         if (fR[0]<fR[4]*10.) return AliPID::ParticleMass(AliPID::kProton);
1027      return 0.00051;
1028   }
1029   if (k==1) return AliPID::ParticleMass(AliPID::kMuon);
1030   if (k==2||k==-1) return AliPID::ParticleMass(AliPID::kPion);
1031   if (k==3) return AliPID::ParticleMass(AliPID::kKaon);
1032   if (k==4) return AliPID::ParticleMass(AliPID::kProton);
1033   AliWarning("Undefined mass !");
1034   return AliPID::ParticleMass(AliPID::kPion);
1035 }
1036
1037 //______________________________________________________________________________
1038 Double_t AliESDtrack::E() const
1039 {
1040   // Returns the energy of the particle given its assumed mass.
1041   // Assumes the pion mass if the particle can't be identified properly.
1042
1043   Double_t m = M();
1044   Double_t p = P();
1045   return TMath::Sqrt(p*p + m*m);
1046 }
1047
1048 //______________________________________________________________________________
1049 Double_t AliESDtrack::Y() const
1050 {
1051   // Returns the rapidity of a particle given its assumed mass.
1052   // Assumes the pion mass if the particle can't be identified properly.
1053
1054   Double_t e = E();
1055   Double_t pz = Pz();
1056   if (e != TMath::Abs(pz)) { // energy was not equal to pz
1057     return 0.5*TMath::Log((e+pz)/(e-pz));
1058   } else { // energy was equal to pz
1059     return -999.;
1060   }
1061 }
1062
1063 //_______________________________________________________________________
1064 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
1065   //
1066   // This function updates track's running parameters
1067   //
1068   Int_t *index=0;
1069   Bool_t rc=kTRUE;
1070
1071   SetStatus(flags);
1072   fLabel=t->GetLabel();
1073
1074   if (t->IsStartedTimeIntegral()) {
1075     SetStatus(kTIME);
1076     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
1077     SetIntegratedLength(t->GetIntegratedLength());
1078   }
1079
1080   Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1081
1082   switch (flags) {
1083
1084   case kITSin: case kITSout: case kITSrefit:
1085     fITSClusterMap=0;
1086     fITSncls=t->GetNumberOfClusters();
1087     index=fFriendTrack->GetITSindices();
1088     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
1089         index[i]=t->GetClusterIndex(i);
1090         if (i<fITSncls) {
1091            Int_t l=(index[i] & 0xf0000000) >> 28;
1092            SETBIT(fITSClusterMap,l);
1093         }
1094     }
1095     fITSchi2=t->GetChi2();
1096     fITSsignal=t->GetPIDsignal();
1097     fITSLabel = t->GetLabel();
1098     // keep in fOp the parameters outside ITS for ITS stand-alone tracks
1099     if (flags==kITSout) {
1100       if (!fOp) fOp=new AliExternalTrackParam(*t);
1101       else
1102         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1103     }
1104     break;
1105
1106   case kTPCin: case kTPCrefit:
1107     fTPCLabel = t->GetLabel();
1108     if (flags==kTPCin)  fTPCInner=new AliExternalTrackParam(*t);
1109     if (!fIp) fIp=new AliExternalTrackParam(*t);
1110     else
1111       fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1112   case kTPCout:
1113     index=fFriendTrack->GetTPCindices();
1114     if (flags & kTPCout){
1115       if (!fOp) fOp=new AliExternalTrackParam(*t);
1116       else
1117         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1118     }
1119     fTPCncls=t->GetNumberOfClusters();
1120     fTPCchi2=t->GetChi2();
1121
1122      {//prevrow must be declared in separate namespace, otherwise compiler cries:
1123       //"jump to case label crosses initialization of `Int_t prevrow'"
1124        Int_t prevrow = -1;
1125        //       for (Int_t i=0;i<fTPCncls;i++)
1126        for (Int_t i=0;i<AliESDfriendTrack::kMaxTPCcluster;i++)
1127         {
1128           index[i]=t->GetClusterIndex(i);
1129           Int_t idx = index[i];
1130
1131           if (idx<0) continue;
1132
1133           // Piotr's Cluster Map for HBT
1134           // ### please change accordingly if cluster array is changing
1135           // to "New TPC Tracking" style (with gaps in array)
1136           Int_t sect = (idx&0xff000000)>>24;
1137           Int_t row = (idx&0x00ff0000)>>16;
1138           if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
1139
1140           fTPCClusterMap.SetBitNumber(row,kTRUE);
1141
1142           //Fill the gap between previous row and this row with 0 bits
1143           //In case  ###  pleas change it as well - just set bit 0 in case there
1144           //is no associated clusters for current "i"
1145           if (prevrow < 0)
1146            {
1147              prevrow = row;//if previous bit was not assigned yet == this is the first one
1148            }
1149           else
1150            { //we don't know the order (inner to outer or reverse)
1151              //just to be save in case it is going to change
1152              Int_t n = 0, m = 0;
1153              if (prevrow < row)
1154               {
1155                 n = prevrow;
1156                 m = row;
1157               }
1158              else
1159               {
1160                 n = row;
1161                 m = prevrow;
1162               }
1163
1164              for (Int_t j = n+1; j < m; j++)
1165               {
1166                 fTPCClusterMap.SetBitNumber(j,kFALSE);
1167               }
1168              prevrow = row;
1169            }
1170           // End Of Piotr's Cluster Map for HBT
1171         }
1172      }
1173     fTPCsignal=t->GetPIDsignal();
1174     break;
1175
1176   case kTRDout: case kTRDin: case kTRDrefit:
1177     index     = fFriendTrack->GetTRDindices();
1178     fTRDLabel = t->GetLabel();
1179     fTRDchi2  = t->GetChi2();
1180     fTRDncls  = t->GetNumberOfClusters();
1181     for (Int_t i=0;i<6;i++) index[i]=t->GetTrackletIndex(i);
1182
1183     fTRDsignal=t->GetPIDsignal();
1184     break;
1185   case kTRDbackup:
1186     if (!fOp) fOp=new AliExternalTrackParam(*t);
1187     else
1188       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1189     fTRDncls0 = t->GetNumberOfClusters();
1190     break;
1191   case kTOFin:
1192     break;
1193   case kTOFout:
1194     break;
1195   case kTRDStop:
1196     break;
1197   default:
1198     AliError("Wrong flag !");
1199     return kFALSE;
1200   }
1201
1202   return rc;
1203 }
1204
1205 //_______________________________________________________________________
1206 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
1207   //---------------------------------------------------------------------
1208   // This function returns external representation of the track parameters
1209   //---------------------------------------------------------------------
1210   x=GetX();
1211   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
1212 }
1213
1214 //_______________________________________________________________________
1215 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
1216   //---------------------------------------------------------------------
1217   // This function returns external representation of the cov. matrix
1218   //---------------------------------------------------------------------
1219   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
1220 }
1221
1222 //_______________________________________________________________________
1223 Bool_t AliESDtrack::GetConstrainedExternalParameters
1224                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1225   //---------------------------------------------------------------------
1226   // This function returns the constrained external track parameters
1227   //---------------------------------------------------------------------
1228   if (!fCp) return kFALSE;
1229   alpha=fCp->GetAlpha();
1230   x=fCp->GetX();
1231   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
1232   return kTRUE;
1233 }
1234
1235 //_______________________________________________________________________
1236 Bool_t
1237 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
1238   //---------------------------------------------------------------------
1239   // This function returns the constrained external cov. matrix
1240   //---------------------------------------------------------------------
1241   if (!fCp) return kFALSE;
1242   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
1243   return kTRUE;
1244 }
1245
1246 Bool_t
1247 AliESDtrack::GetInnerExternalParameters
1248                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1249   //---------------------------------------------------------------------
1250   // This function returns external representation of the track parameters
1251   // at the inner layer of TPC
1252   //---------------------------------------------------------------------
1253   if (!fIp) return kFALSE;
1254   alpha=fIp->GetAlpha();
1255   x=fIp->GetX();
1256   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
1257   return kTRUE;
1258 }
1259
1260 Bool_t
1261 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
1262  //---------------------------------------------------------------------
1263  // This function returns external representation of the cov. matrix
1264  // at the inner layer of TPC
1265  //---------------------------------------------------------------------
1266   if (!fIp) return kFALSE;
1267   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
1268   return kTRUE;
1269 }
1270
1271 void
1272 AliESDtrack::SetOuterParam(const AliExternalTrackParam *p, ULong_t flags) {
1273   //
1274   // This is a direct setter for the outer track parameters
1275   //
1276   SetStatus(flags);
1277   if (fOp) delete fOp;
1278   fOp=new AliExternalTrackParam(*p);
1279 }
1280
1281 Bool_t
1282 AliESDtrack::GetOuterExternalParameters
1283                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1284   //---------------------------------------------------------------------
1285   // This function returns external representation of the track parameters
1286   // at the inner layer of TRD
1287   //---------------------------------------------------------------------
1288   if (!fOp) return kFALSE;
1289   alpha=fOp->GetAlpha();
1290   x=fOp->GetX();
1291   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
1292   return kTRUE;
1293 }
1294
1295 Bool_t
1296 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
1297  //---------------------------------------------------------------------
1298  // This function returns external representation of the cov. matrix
1299  // at the inner layer of TRD
1300  //---------------------------------------------------------------------
1301   if (!fOp) return kFALSE;
1302   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
1303   return kTRUE;
1304 }
1305
1306 Int_t AliESDtrack::GetNcls(Int_t idet) const
1307 {
1308   // Get number of clusters by subdetector index
1309   //
1310   Int_t ncls = 0;
1311   switch(idet){
1312   case 0:
1313     ncls = fITSncls;
1314     break;
1315   case 1:
1316     ncls = fTPCncls;
1317     break;
1318   case 2:
1319     ncls = fTRDncls;
1320     break;
1321   case 3:
1322     if (fTOFindex != -1)
1323       ncls = 1;
1324     break;
1325   case 4: //PHOS
1326     break;
1327   case 5: //HMPID
1328     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1329       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1330         ncls = 1;
1331       }
1332     }
1333     break;
1334   default:
1335     break;
1336   }
1337   return ncls;
1338 }
1339
1340 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
1341 {
1342   // Get cluster index array by subdetector index
1343   //
1344   Int_t ncls = 0;
1345   switch(idet){
1346   case 0:
1347     ncls = GetITSclusters(idx);
1348     break;
1349   case 1:
1350     ncls = GetTPCclusters(idx);
1351     break;
1352   case 2:
1353     ncls = GetTRDclusters(idx);
1354     break;
1355   case 3:
1356     if (fTOFindex != -1) {
1357       idx[0] = fTOFindex;
1358       ncls = 1;
1359     }
1360     break;
1361   case 4: //PHOS
1362     break;
1363   case 5:
1364     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1365       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1366         idx[0] = GetHMPIDcluIdx();
1367         ncls = 1;
1368       }
1369     }
1370     break;
1371   case 6: //EMCAL
1372     break;
1373   default:
1374     break;
1375   }
1376   return ncls;
1377 }
1378
1379 //_______________________________________________________________________
1380 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
1381   // Returns the array with integrated times for each particle hypothesis
1382   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
1383 }
1384
1385 //_______________________________________________________________________
1386 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
1387   // Sets the array with integrated times for each particle hypotesis
1388   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
1389 }
1390
1391 //_______________________________________________________________________
1392 void AliESDtrack::SetITSpid(const Double_t *p) {
1393   // Sets values for the probability of each particle type (in ITS)
1394   SetPIDValues(fITSr,p,AliPID::kSPECIES);
1395   SetStatus(AliESDtrack::kITSpid);
1396 }
1397
1398 //_______________________________________________________________________
1399 void AliESDtrack::GetITSpid(Double_t *p) const {
1400   // Gets the probability of each particle type (in ITS)
1401   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
1402 }
1403
1404 //_______________________________________________________________________
1405 Char_t AliESDtrack::GetITSclusters(Int_t *idx) const {
1406   //---------------------------------------------------------------------
1407   // This function returns indices of the assgined ITS clusters
1408   //---------------------------------------------------------------------
1409   if (idx!=0) {
1410      Int_t *index=fFriendTrack->GetITSindices();
1411      for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) {
1412          if ( (i>=fITSncls) && (i<6) ) idx[i]=-1;
1413          else idx[i]=index[i];
1414      }
1415   }
1416   return fITSncls;
1417 }
1418
1419 //_______________________________________________________________________
1420 Bool_t AliESDtrack::GetITSModuleIndexInfo(Int_t ilayer,Int_t &idet,Int_t &status,
1421                                          Float_t &xloc,Float_t &zloc) const {
1422   //----------------------------------------------------------------------
1423   // This function encodes in the module number also the status of cluster association
1424   // "status" can have the following values:
1425   // 1 "found" (cluster is associated),
1426   // 2 "dead" (module is dead from OCDB),
1427   // 3 "skipped" (module or layer forced to be skipped),
1428   // 4 "outinz" (track out of z acceptance),
1429   // 5 "nocls" (no clusters in the road),
1430   // 6 "norefit" (cluster rejected during refit),
1431   // 7 "deadzspd" (holes in z in SPD)
1432   // Also given are the coordinates of the crossing point of track and module
1433   // (in the local module ref. system)
1434   // WARNING: THIS METHOD HAS TO BE SYNCHRONIZED WITH AliITStrackV2::GetModuleIndexInfo()!
1435   //----------------------------------------------------------------------
1436
1437   if(fITSModule[ilayer]==-1) {
1438     idet = -1;
1439     status=0;
1440     xloc=-99.; zloc=-99.;
1441     return kFALSE;
1442   }
1443
1444   Int_t module = fITSModule[ilayer];
1445
1446   idet = Int_t(module/1000000);
1447
1448   module -= idet*1000000;
1449
1450   status = Int_t(module/100000);
1451
1452   module -= status*100000;
1453
1454   Int_t signs = Int_t(module/10000);
1455
1456   module-=signs*10000;
1457
1458   Int_t xInt = Int_t(module/100);
1459   module -= xInt*100;
1460
1461   Int_t zInt = module;
1462
1463   if(signs==1) { xInt*=1; zInt*=1; }
1464   if(signs==2) { xInt*=1; zInt*=-1; }
1465   if(signs==3) { xInt*=-1; zInt*=1; }
1466   if(signs==4) { xInt*=-1; zInt*=-1; }
1467
1468   xloc = 0.1*(Float_t)xInt;
1469   zloc = 0.1*(Float_t)zInt;
1470
1471   if(status==4) idet = -1;
1472
1473   return kTRUE;
1474 }
1475
1476 //_______________________________________________________________________
1477 UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
1478   //---------------------------------------------------------------------
1479   // This function returns indices of the assgined ITS clusters
1480   //---------------------------------------------------------------------
1481   if (idx!=0) {
1482     Int_t *index=fFriendTrack->GetTPCindices();
1483     for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
1484   }
1485   return fTPCncls;
1486 }
1487
1488 Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
1489   //
1490   // GetDensity of the clusters on given region between row0 and row1
1491   // Dead zone effect takin into acoount
1492   //
1493   Int_t good  = 0;
1494   Int_t found = 0;
1495   //
1496   Int_t *index=fFriendTrack->GetTPCindices();
1497   for (Int_t i=row0;i<=row1;i++){
1498     Int_t idx = index[i];
1499     if (idx!=-1)  good++;             // track outside of dead zone
1500     if (idx>0)    found++;
1501   }
1502   Float_t density=0.5;
1503   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
1504   return density;
1505 }
1506
1507 //_______________________________________________________________________
1508 void AliESDtrack::SetTPCpid(const Double_t *p) {
1509   // Sets values for the probability of each particle type (in TPC)
1510   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
1511   SetStatus(AliESDtrack::kTPCpid);
1512 }
1513
1514 //_______________________________________________________________________
1515 void AliESDtrack::GetTPCpid(Double_t *p) const {
1516   // Gets the probability of each particle type (in TPC)
1517   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
1518 }
1519
1520 //_______________________________________________________________________
1521 UChar_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
1522   //---------------------------------------------------------------------
1523   // This function returns indices of the assgined TRD clusters
1524   //---------------------------------------------------------------------
1525   if (idx!=0) {
1526      Int_t *index=fFriendTrack->GetTRDindices();
1527      for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
1528   }
1529   return fTRDncls;
1530 }
1531
1532 //_______________________________________________________________________
1533 UChar_t AliESDtrack::GetTRDtracklets(Int_t *idx) const {
1534   //---------------------------------------------------------------------
1535   // This function returns indices of the assigned TRD tracklets
1536   //---------------------------------------------------------------------
1537   if (idx!=0) {
1538      Int_t *index=fFriendTrack->GetTRDindices();
1539      for (Int_t i=0; i<6/*AliESDfriendTrack::kMaxTRDcluster*/; i++) idx[i]=index[i];
1540   }
1541   return fTRDncls;
1542 }
1543
1544 //_______________________________________________________________________
1545 void AliESDtrack::SetTRDpid(const Double_t *p) {
1546   // Sets values for the probability of each particle type (in TRD)
1547   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
1548   SetStatus(AliESDtrack::kTRDpid);
1549 }
1550
1551 //_______________________________________________________________________
1552 void AliESDtrack::GetTRDpid(Double_t *p) const {
1553   // Gets the probability of each particle type (in TRD)
1554   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
1555 }
1556
1557 //_______________________________________________________________________
1558 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
1559 {
1560   // Sets the probability of particle type iSpecies to p (in TRD)
1561   fTRDr[iSpecies] = p;
1562 }
1563
1564 Double_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
1565 {
1566   // Returns the probability of particle type iSpecies (in TRD)
1567   return fTRDr[iSpecies];
1568 }
1569
1570 void  AliESDtrack::SetNumberOfTRDslices(Int_t n) {
1571   //Sets the number of slices used for PID
1572   if (fTRDnSlices != 0) return;
1573   fTRDnSlices=kTRDnPlanes*n;
1574   fTRDslices=new Double32_t[fTRDnSlices];
1575   for (Int_t i=0; i<fTRDnSlices; i++) fTRDslices[i]=-1.;
1576 }
1577
1578 void  AliESDtrack::SetTRDslice(Double_t q, Int_t plane, Int_t slice) {
1579   //Sets the charge q in the slice of the plane
1580   Int_t ns=GetNumberOfTRDslices();
1581   if (ns==0) {
1582     AliError("No TRD slices allocated for this track !");
1583     return;
1584   }
1585
1586   if ((plane<0) || (plane>=kTRDnPlanes)) {
1587     AliError("Wrong TRD plane !");
1588     return;
1589   }
1590   if ((slice<0) || (slice>=ns)) {
1591     AliError("Wrong TRD slice !");
1592     return;
1593   }
1594   Int_t n=plane*ns + slice;
1595   fTRDslices[n]=q;
1596 }
1597
1598 Double_t  AliESDtrack::GetTRDslice(Int_t plane, Int_t slice) const {
1599   //Gets the charge from the slice of the plane
1600   Int_t ns=GetNumberOfTRDslices();
1601   if (ns==0) {
1602     //AliError("No TRD slices allocated for this track !");
1603     return -1.;
1604   }
1605
1606   if ((plane<0) || (plane>=kTRDnPlanes)) {
1607     AliError("Wrong TRD plane !");
1608     return -1.;
1609   }
1610   if ((slice<-1) || (slice>=ns)) {
1611     //AliError("Wrong TRD slice !");
1612     return -1.;
1613   }
1614
1615   if (slice==-1) {
1616     Double_t q=0.;
1617     for (Int_t i=0; i<ns; i++) q+=fTRDslices[plane*ns + i];
1618     return q/ns;
1619   }
1620
1621   return fTRDslices[plane*ns + slice];
1622 }
1623
1624
1625 //_______________________________________________________________________
1626 void AliESDtrack::SetTOFpid(const Double_t *p) {
1627   // Sets the probability of each particle type (in TOF)
1628   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
1629   SetStatus(AliESDtrack::kTOFpid);
1630 }
1631
1632 //_______________________________________________________________________
1633 void AliESDtrack::SetTOFLabel(const Int_t *p) {
1634   // Sets  (in TOF)
1635   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
1636 }
1637
1638 //_______________________________________________________________________
1639 void AliESDtrack::GetTOFpid(Double_t *p) const {
1640   // Gets probabilities of each particle type (in TOF)
1641   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
1642 }
1643
1644 //_______________________________________________________________________
1645 void AliESDtrack::GetTOFLabel(Int_t *p) const {
1646   // Gets (in TOF)
1647   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
1648 }
1649
1650 //_______________________________________________________________________
1651 void AliESDtrack::GetTOFInfo(Float_t *info) const {
1652   // Gets (in TOF)
1653   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
1654 }
1655
1656 //_______________________________________________________________________
1657 void AliESDtrack::SetTOFInfo(Float_t*info) {
1658   // Gets (in TOF)
1659   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
1660 }
1661
1662
1663
1664 //_______________________________________________________________________
1665 void AliESDtrack::SetHMPIDpid(const Double_t *p) {
1666   // Sets the probability of each particle type (in HMPID)
1667   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
1668   SetStatus(AliESDtrack::kHMPIDpid);
1669 }
1670
1671 //_______________________________________________________________________
1672 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
1673   // Gets probabilities of each particle type (in HMPID)
1674   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
1675 }
1676
1677
1678
1679 //_______________________________________________________________________
1680 void AliESDtrack::SetESDpid(const Double_t *p) {
1681   // Sets the probability of each particle type for the ESD track
1682   SetPIDValues(fR,p,AliPID::kSPECIES);
1683   SetStatus(AliESDtrack::kESDpid);
1684 }
1685
1686 //_______________________________________________________________________
1687 void AliESDtrack::GetESDpid(Double_t *p) const {
1688   // Gets probability of each particle type for the ESD track
1689   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
1690 }
1691
1692 //_______________________________________________________________________
1693 Bool_t AliESDtrack::RelateToVertexTPC(const AliESDVertex *vtx,
1694 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
1695   //
1696   // Try to relate the TPC-only track parameters to the vertex "vtx",
1697   // if the (rough) transverse impact parameter is not bigger then "maxd".
1698   //            Magnetic field is "b" (kG).
1699   //
1700   // a) The TPC-only paramters are extapolated to the DCA to the vertex.
1701   // b) The impact parameters and their covariance matrix are calculated.
1702   // c) An attempt to constrain the TPC-only params to the vertex is done.
1703   //    The constrained params are returned via "cParam".
1704   //
1705   // In the case of success, the returned value is kTRUE
1706   // otherwise, it's kFALSE)
1707   //
1708
1709   if (!fTPCInner) return kFALSE;
1710   if (!vtx) return kFALSE;
1711
1712   Double_t dz[2],cov[3];
1713   if (!fTPCInner->PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
1714
1715   fdTPC = dz[0];
1716   fzTPC = dz[1];
1717   fCddTPC = cov[0];
1718   fCdzTPC = cov[1];
1719   fCzzTPC = cov[2];
1720
1721   Double_t covar[6]; vtx->GetCovMatrix(covar);
1722   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
1723   Double_t c[3]={covar[2],0.,covar[5]};
1724
1725   Double_t chi2=GetPredictedChi2(p,c);
1726   if (chi2>kVeryBig) return kFALSE;
1727
1728   fCchi2TPC=chi2;
1729
1730   if (!cParam) return kTRUE;
1731
1732   *cParam = *fTPCInner;
1733   if (!cParam->Update(p,c)) return kFALSE;
1734
1735   return kTRUE;
1736 }
1737
1738 //_______________________________________________________________________
1739 Bool_t AliESDtrack::RelateToVertex(const AliESDVertex *vtx,
1740 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
1741   //
1742   // Try to relate this track to the vertex "vtx",
1743   // if the (rough) transverse impact parameter is not bigger then "maxd".
1744   //            Magnetic field is "b" (kG).
1745   //
1746   // a) The track gets extapolated to the DCA to the vertex.
1747   // b) The impact parameters and their covariance matrix are calculated.
1748   // c) An attempt to constrain this track to the vertex is done.
1749   //    The constrained params are returned via "cParam".
1750   //
1751   // In the case of success, the returned value is kTRUE
1752   // (otherwise, it's kFALSE)
1753   //
1754
1755   if (!vtx) return kFALSE;
1756
1757   Double_t dz[2],cov[3];
1758   if (!PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
1759
1760   fD = dz[0];
1761   fZ = dz[1];
1762   fCdd = cov[0];
1763   fCdz = cov[1];
1764   fCzz = cov[2];
1765
1766   Double_t covar[6]; vtx->GetCovMatrix(covar);
1767   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
1768   Double_t c[3]={covar[2],0.,covar[5]};
1769
1770   Double_t chi2=GetPredictedChi2(p,c);
1771   if (chi2>kVeryBig) return kFALSE;
1772
1773   fCchi2=chi2;
1774
1775
1776   //--- Could now these lines be removed ? ---
1777   delete fCp;
1778   fCp=new AliExternalTrackParam(*this);
1779
1780   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
1781   //----------------------------------------
1782
1783
1784   if (!cParam) return kTRUE;
1785
1786   *cParam = *this;
1787   if (!cParam->Update(p,c)) return kFALSE;
1788
1789   return kTRUE;
1790 }
1791
1792 //_______________________________________________________________________
1793 void AliESDtrack::Print(Option_t *) const {
1794   // Prints info on the track
1795   AliExternalTrackParam::Print();
1796   printf("ESD track info\n") ;
1797   Double_t p[AliPID::kSPECIESN] ;
1798   Int_t index = 0 ;
1799   if( IsOn(kITSpid) ){
1800     printf("From ITS: ") ;
1801     GetITSpid(p) ;
1802     for(index = 0 ; index < AliPID::kSPECIES; index++)
1803       printf("%f, ", p[index]) ;
1804     printf("\n           signal = %f\n", GetITSsignal()) ;
1805   }
1806   if( IsOn(kTPCpid) ){
1807     printf("From TPC: ") ;
1808     GetTPCpid(p) ;
1809     for(index = 0 ; index < AliPID::kSPECIES; index++)
1810       printf("%f, ", p[index]) ;
1811     printf("\n           signal = %f\n", GetTPCsignal()) ;
1812   }
1813   if( IsOn(kTRDpid) ){
1814     printf("From TRD: ") ;
1815     GetTRDpid(p) ;
1816     for(index = 0 ; index < AliPID::kSPECIES; index++)
1817       printf("%f, ", p[index]) ;
1818       printf("\n           signal = %f\n", GetTRDsignal()) ;
1819   }
1820   if( IsOn(kTOFpid) ){
1821     printf("From TOF: ") ;
1822     GetTOFpid(p) ;
1823     for(index = 0 ; index < AliPID::kSPECIES; index++)
1824       printf("%f, ", p[index]) ;
1825     printf("\n           signal = %f\n", GetTOFsignal()) ;
1826   }
1827   if( IsOn(kHMPIDpid) ){
1828     printf("From HMPID: ") ;
1829     GetHMPIDpid(p) ;
1830     for(index = 0 ; index < AliPID::kSPECIES; index++)
1831       printf("%f, ", p[index]) ;
1832     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
1833   }
1834 }
1835
1836
1837 //
1838 // Draw functionality
1839 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
1840 //
1841 void AliESDtrack::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
1842   //
1843   // Fill points in the polymarker
1844   //
1845   TObjArray arrayRef;
1846   arrayRef.AddLast(new AliExternalTrackParam(*this));
1847   if (fIp) arrayRef.AddLast(new AliExternalTrackParam(*fIp));
1848   if (fOp) arrayRef.AddLast(new AliExternalTrackParam(*fOp));
1849   //
1850   Double_t mpos[3]={0,0,0};
1851   Int_t entries=arrayRef.GetEntries();
1852   for (Int_t i=0;i<entries;i++){
1853     Double_t pos[3];
1854     ((AliExternalTrackParam*)arrayRef.At(i))->GetXYZ(pos);
1855     mpos[0]+=pos[0]/entries;
1856     mpos[1]+=pos[1]/entries;
1857     mpos[2]+=pos[2]/entries;
1858   }
1859   // Rotate to the mean position
1860   //
1861   Float_t fi= TMath::ATan2(mpos[1],mpos[0]);
1862   for (Int_t i=0;i<entries;i++){
1863     Bool_t res = ((AliExternalTrackParam*)arrayRef.At(i))->Rotate(fi);
1864     if (!res) delete arrayRef.RemoveAt(i);
1865   }
1866   Int_t counter=0;
1867   for (Double_t r=minR; r<maxR; r+=stepR){
1868     Double_t sweight=0;
1869     Double_t mlpos[3]={0,0,0};
1870     for (Int_t i=0;i<entries;i++){
1871       Double_t point[3]={0,0,0};
1872       AliExternalTrackParam *param = ((AliExternalTrackParam*)arrayRef.At(i));
1873       if (!param) continue;
1874       if (param->GetXYZAt(r,magF,point)){
1875         Double_t weight = 1./(10.+(r-param->GetX())*(r-param->GetX()));
1876         sweight+=weight;
1877         mlpos[0]+=point[0]*weight;
1878         mlpos[1]+=point[1]*weight;
1879         mlpos[2]+=point[2]*weight;
1880       }
1881     }
1882     if (sweight>0){
1883       mlpos[0]/=sweight;
1884       mlpos[1]/=sweight;
1885       mlpos[2]/=sweight;
1886       pol->SetPoint(counter,mlpos[0],mlpos[1], mlpos[2]);
1887       printf("xyz\t%f\t%f\t%f\n",mlpos[0], mlpos[1],mlpos[2]);
1888       counter++;
1889     }
1890   }
1891 }