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