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