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 #include <TDatabasePDG.h>
 119
 120 #include "AliESDVertex.h"
 121 #include "AliESDtrack.h"
 122 #include "AliESDEvent.h"
 123 #include "AliKalmanTrack.h"
 124 #include "AliVTrack.h"
 125 #include "AliLog.h"
 126 #include "AliTrackPointArray.h"
 127 #include "TPolyMarker3D.h"
 128
 129 ClassImp(AliESDtrack)
 130
 131 void SetPIDValues(Double_t * dest, const Double_t * src, Int_t n) {
 132   // This function copies "n" PID weights from "scr" to "dest"
 133   // and normalizes their sum to 1 thus producing conditional probabilities.
 134   // The negative weights are set to 0.
 135   // In case all the weights are non-positive they are replaced by
 136   // uniform probabilities
 137
 138   if (n<=0) return;
 139
 140   Float_t uniform = 1./(Float_t)n;
 141
 142   Float_t sum = 0;
 143   for (Int_t i=0; i<n; i++)
 144     if (src[i]>=0) {
 145       sum+=src[i];
 146       dest[i] = src[i];
 147     }
 148     else {
 149       dest[i] = 0;
 150     }
 151
 152   if(sum>0)
 153     for (Int_t i=0; i<n; i++) dest[i] /= sum;
 154   else
 155     for (Int_t i=0; i<n; i++) dest[i] = uniform;
 156 }
 157
 158 //_______________________________________________________________________
 159 AliESDtrack::AliESDtrack() :
 160   AliExternalTrackParam(),
 161   fCp(0),
 162   fIp(0),
 163   fTPCInner(0),
 164   fOp(0),
 165   fHMPIDp(0),
 166   fFriendTrack(NULL),
 167   fTPCClusterMap(159),//number of padrows
 168   fTPCSharedMap(159),//number of padrows
 169   fFlags(0),
 170   fID(0),
 171   fLabel(0),
 172   fITSLabel(0),
 173   fTPCLabel(0),
 174   fTRDLabel(0),
 175   fTOFCalChannel(-1),
 176   fTOFindex(-1),
 177   fHMPIDqn(0),
 178   fHMPIDcluIdx(-1),
 179   fCaloIndex(kEMCALNoMatch),
 180   fHMPIDtrkTheta(0),
 181   fHMPIDtrkPhi(0),
 182   fHMPIDsignal(0),
 183   fTrackLength(0),
 184   fdTPC(0),fzTPC(0),
 185   fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
 186   fCchi2TPC(0),
 187   fD(0),fZ(0),
 188   fCdd(0),fCdz(0),fCzz(0),
 189   fCchi2(0),
 190   fITSchi2(0),
 191   fTPCchi2(0),
 192   fTPCchi2Iter1(0),
 193   fTRDchi2(0),
 194   fTOFchi2(0),
 195   fHMPIDchi2(0),
 196   fGlobalChi2(0),
 197   fITSsignal(0),
 198   fTPCsignal(0),
 199   fTPCsignalS(0),
 200   fTRDsignal(0),
 201   fTRDQuality(0),
 202   fTRDBudget(0),
 203   fTOFsignal(99999),
 204   fTOFsignalToT(99999),
 205   fTOFsignalRaw(99999),
 206   fTOFsignalDz(999),
 207   fTOFsignalDx(999),
 208   fTOFdeltaBC(999),
 209   fTOFl0l1(999),
 210   fCaloDx(0),
 211   fCaloDz(0),
 212   fHMPIDtrkX(0),
 213   fHMPIDtrkY(0),
 214   fHMPIDmipX(0),
 215   fHMPIDmipY(0),
 216   fTPCncls(0),
 217   fTPCnclsF(0),
 218   fTPCsignalN(0),
 219   fTPCnclsIter1(0),
 220   fTPCnclsFIter1(0),
 221   fITSncls(0),
 222   fITSClusterMap(0),
 223   fITSSharedMap(0),
 224   fTRDncls(0),
 225   fTRDncls0(0),
 226   fTRDntracklets(0),
 227   fTRDnSlices(0),
 228   fTRDslices(0x0),
 229   fVertexID(-2),// -2 means an orphan track
 230   fESDEvent(0)
 231 {
 232   //
 233   // The default ESD constructor
 234   //
 235   if (!OnlineMode()) fFriendTrack=new AliESDfriendTrack();
 236
 237   Int_t i;
 238   for (i=0; i<AliPID::kSPECIES; i++) {
 239     fTrackTime[i]=0.;
 240     fR[i]=0.;
 241     fITSr[i]=0.;
 242     fTPCr[i]=0.;
 243     fTRDr[i]=0.;
 244     fTOFr[i]=0.;
 245     fHMPIDr[i]=0.;
 246   }
 247
 248   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
 249   for (i=0; i<3; i++)   { fV0Indexes[i]=0;}
 250   for (i=0;i<kTRDnPlanes;i++) {
 251     fTRDTimBin[i]=0;
 252   }
 253   for (i=0;i<4;i++) {fITSdEdxSamples[i]=0.;}
 254   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
 255   for (i=0;i<3;i++) {fTOFLabel[i]=-1;}
 256   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
 257   for (i=0;i<12;i++) {fITSModule[i]=-1;}
 258 }
 259
 260 bool AliESDtrack::fgkOnlineMode=false;
 261
 262 //_______________________________________________________________________
 263 AliESDtrack::AliESDtrack(const AliESDtrack& track):
 264   AliExternalTrackParam(track),
 265   fCp(0),
 266   fIp(0),
 267   fTPCInner(0),
 268   fOp(0),
 269   fHMPIDp(0),
 270   fFriendTrack(0),
 271   fTPCClusterMap(track.fTPCClusterMap),
 272   fTPCSharedMap(track.fTPCSharedMap),
 273   fFlags(track.fFlags),
 274   fID(track.fID),
 275   fLabel(track.fLabel),
 276   fITSLabel(track.fITSLabel),
 277   fTPCLabel(track.fTPCLabel),
 278   fTRDLabel(track.fTRDLabel),
 279   fTOFCalChannel(track.fTOFCalChannel),
 280   fTOFindex(track.fTOFindex),
 281   fHMPIDqn(track.fHMPIDqn),
 282   fHMPIDcluIdx(track.fHMPIDcluIdx),
 283   fCaloIndex(track.fCaloIndex),
 284   fHMPIDtrkTheta(track.fHMPIDtrkTheta),
 285   fHMPIDtrkPhi(track.fHMPIDtrkPhi),
 286   fHMPIDsignal(track.fHMPIDsignal),
 287   fTrackLength(track.fTrackLength),
 288   fdTPC(track.fdTPC),fzTPC(track.fzTPC),
 289   fCddTPC(track.fCddTPC),fCdzTPC(track.fCdzTPC),fCzzTPC(track.fCzzTPC),
 290   fCchi2TPC(track.fCchi2TPC),
 291   fD(track.fD),fZ(track.fZ),
 292   fCdd(track.fCdd),fCdz(track.fCdz),fCzz(track.fCzz),
 293   fCchi2(track.fCchi2),
 294   fITSchi2(track.fITSchi2),
 295   fTPCchi2(track.fTPCchi2),
 296   fTPCchi2Iter1(track.fTPCchi2Iter1),
 297   fTRDchi2(track.fTRDchi2),
 298   fTOFchi2(track.fTOFchi2),
 299   fHMPIDchi2(track.fHMPIDchi2),
 300   fGlobalChi2(track.fGlobalChi2),
 301   fITSsignal(track.fITSsignal),
 302   fTPCsignal(track.fTPCsignal),
 303   fTPCsignalS(track.fTPCsignalS),
 304   fTRDsignal(track.fTRDsignal),
 305   fTRDQuality(track.fTRDQuality),
 306   fTRDBudget(track.fTRDBudget),
 307   fTOFsignal(track.fTOFsignal),
 308   fTOFsignalToT(track.fTOFsignalToT),
 309   fTOFsignalRaw(track.fTOFsignalRaw),
 310   fTOFsignalDz(track.fTOFsignalDz),
 311   fTOFsignalDx(track.fTOFsignalDx),
 312   fTOFdeltaBC(track.fTOFdeltaBC),
 313   fTOFl0l1(track.fTOFl0l1),
 314   fCaloDx(track.fCaloDx),
 315   fCaloDz(track.fCaloDz),
 316   fHMPIDtrkX(track.fHMPIDtrkX),
 317   fHMPIDtrkY(track.fHMPIDtrkY),
 318   fHMPIDmipX(track.fHMPIDmipX),
 319   fHMPIDmipY(track.fHMPIDmipY),
 320   fTPCncls(track.fTPCncls),
 321   fTPCnclsF(track.fTPCnclsF),
 322   fTPCsignalN(track.fTPCsignalN),
 323   fTPCnclsIter1(track.fTPCnclsIter1),
 324   fTPCnclsFIter1(track.fTPCnclsIter1),
 325   fITSncls(track.fITSncls),
 326   fITSClusterMap(track.fITSClusterMap),
 327   fITSSharedMap(track.fITSSharedMap),
 328   fTRDncls(track.fTRDncls),
 329   fTRDncls0(track.fTRDncls0),
 330   fTRDntracklets(track.fTRDntracklets),
 331   fTRDnSlices(track.fTRDnSlices),
 332   fTRDslices(0x0),
 333   fVertexID(track.fVertexID),
 334   fESDEvent(track.fESDEvent)
 335 {
 336   //
 337   //copy constructor
 338   //
 339   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i]=track.fTrackTime[i];
 340   for (Int_t i=0;i<AliPID::kSPECIES;i++)  fR[i]=track.fR[i];
 341   //
 342   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=track.fITSr[i];
 343   //
 344   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=track.fTPCr[i];
 345   for (Int_t i=0;i<4;i++) {fITSdEdxSamples[i]=track.fITSdEdxSamples[i];}
 346   for (Int_t i=0;i<4;i++) {fTPCPoints[i]=track.fTPCPoints[i];}
 347   for (Int_t i=0; i<3;i++)   { fKinkIndexes[i]=track.fKinkIndexes[i];}
 348   for (Int_t i=0; i<3;i++)   { fV0Indexes[i]=track.fV0Indexes[i];}
 349   //
 350   for (Int_t i=0;i<kTRDnPlanes;i++) {
 351     fTRDTimBin[i]=track.fTRDTimBin[i];
 352   }
 353
 354   if (fTRDnSlices) {
 355     fTRDslices=new Double32_t[fTRDnSlices];
 356     for (Int_t i=0; i<fTRDnSlices; i++) fTRDslices[i]=track.fTRDslices[i];
 357   }
 358
 359   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i]=track.fTRDr[i];
 360   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i]=track.fTOFr[i];
 361   for (Int_t i=0;i<3;i++) fTOFLabel[i]=track.fTOFLabel[i];
 362   for (Int_t i=0;i<10;i++) fTOFInfo[i]=track.fTOFInfo[i];
 363   for (Int_t i=0;i<12;i++) fITSModule[i]=track.fITSModule[i];
 364   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i]=track.fHMPIDr[i];
 365
 366   if (track.fCp) fCp=new AliExternalTrackParam(*track.fCp);
 367   if (track.fIp) fIp=new AliExternalTrackParam(*track.fIp);
 368   if (track.fTPCInner) fTPCInner=new AliExternalTrackParam(*track.fTPCInner);
 369   if (track.fOp) fOp=new AliExternalTrackParam(*track.fOp);
 370   if (track.fHMPIDp) fHMPIDp=new AliExternalTrackParam(*track.fHMPIDp);
 371
 372   if (track.fFriendTrack) fFriendTrack=new AliESDfriendTrack(*(track.fFriendTrack));
 373 }
 374
 375 //_______________________________________________________________________
 376 AliESDtrack::AliESDtrack(const AliVTrack *track) :
 377   AliExternalTrackParam(track),
 378   fCp(0),
 379   fIp(0),
 380   fTPCInner(0),
 381   fOp(0),
 382   fHMPIDp(0),
 383   fFriendTrack(0),
 384   fTPCClusterMap(159),//number of padrows
 385   fTPCSharedMap(159),//number of padrows
 386   fFlags(0),
 387   fID(),
 388   fLabel(0),
 389   fITSLabel(0),
 390   fTPCLabel(0),
 391   fTRDLabel(0),
 392   fTOFCalChannel(-1),
 393   fTOFindex(-1),
 394   fHMPIDqn(0),
 395   fHMPIDcluIdx(-1),
 396   fCaloIndex(kEMCALNoMatch),
 397   fHMPIDtrkTheta(0),
 398   fHMPIDtrkPhi(0),
 399   fHMPIDsignal(0),
 400   fTrackLength(0),
 401   fdTPC(0),fzTPC(0),
 402   fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
 403   fCchi2TPC(0),
 404   fD(0),fZ(0),
 405   fCdd(0),fCdz(0),fCzz(0),
 406   fCchi2(0),
 407   fITSchi2(0),
 408   fTPCchi2(0),
 409   fTPCchi2Iter1(0),
 410   fTRDchi2(0),
 411   fTOFchi2(0),
 412   fHMPIDchi2(0),
 413   fGlobalChi2(0),
 414   fITSsignal(0),
 415   fTPCsignal(0),
 416   fTPCsignalS(0),
 417   fTRDsignal(0),
 418   fTRDQuality(0),
 419   fTRDBudget(0),
 420   fTOFsignal(99999),
 421   fTOFsignalToT(99999),
 422   fTOFsignalRaw(99999),
 423   fTOFsignalDz(999),
 424   fTOFsignalDx(999),
 425   fTOFdeltaBC(999),
 426   fTOFl0l1(999),
 427   fCaloDx(0),
 428   fCaloDz(0),
 429   fHMPIDtrkX(0),
 430   fHMPIDtrkY(0),
 431   fHMPIDmipX(0),
 432   fHMPIDmipY(0),
 433   fTPCncls(0),
 434   fTPCnclsF(0),
 435   fTPCsignalN(0),
 436   fTPCnclsIter1(0),
 437   fTPCnclsFIter1(0),
 438   fITSncls(0),
 439   fITSClusterMap(0),
 440   fITSSharedMap(0),
 441   fTRDncls(0),
 442   fTRDncls0(0),
 443   fTRDntracklets(0),
 444   fTRDnSlices(0),
 445   fTRDslices(0x0),
 446   fVertexID(-2),  // -2 means an orphan track
 447   fESDEvent(0)
 448 {
 449   //
 450   // ESD track from AliVTrack.
 451   // This is not a copy constructor !
 452   //
 453
 454   if (track->InheritsFrom("AliExternalTrackParam")) {
 455      AliError("This is not a copy constructor. Use AliESDtrack(const AliESDtrack &) !");
 456      AliWarning("Calling the default constructor...");
 457      AliESDtrack();
 458      return;
 459   }
 460
 461   // Reset all the arrays
 462   Int_t i;
 463   for (i=0; i<AliPID::kSPECIES; i++) {
 464     fTrackTime[i]=0.;
 465     fR[i]=0.;
 466     fITSr[i]=0.;
 467     fTPCr[i]=0.;
 468     fTRDr[i]=0.;
 469     fTOFr[i]=0.;
 470     fHMPIDr[i]=0.;
 471   }
 472
 473   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
 474   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
 475   for (i=0;i<kTRDnPlanes;i++) {
 476     fTRDTimBin[i]=0;
 477   }
 478   for (i=0;i<4;i++) {fITSdEdxSamples[i]=0.;}
 479   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
 480   for (i=0;i<3;i++) {fTOFLabel[i]=-1;}
 481   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
 482   for (i=0;i<12;i++) {fITSModule[i]=-1;}
 483
 484   // Set the ID
 485   SetID(track->GetID());
 486
 487   // Set ITS cluster map
 488   fITSClusterMap=track->GetITSClusterMap();
 489   fITSSharedMap=0;
 490
 491   fITSncls=0;
 492   for(i=0; i<6; i++) {
 493     if(HasPointOnITSLayer(i)) fITSncls++;
 494   }
 495
 496   // Set TPC ncls
 497   fTPCncls=track->GetTPCNcls();
 498
 499
 500   // Set the combined PID
 501   const Double_t *pid = track->PID();
 502   if(pid){
 503     for (i=0; i<AliPID::kSPECIES; i++) fR[i]=pid[i];
 504   }
 505   // AliESD track label
 506   SetLabel(track->GetLabel());
 507   // Set the status
 508   SetStatus(track->GetStatus());
 509 }
 510
 511 //_______________________________________________________________________
 512 AliESDtrack::AliESDtrack(TParticle * part) :
 513   AliExternalTrackParam(),
 514   fCp(0),
 515   fIp(0),
 516   fTPCInner(0),
 517   fOp(0),
 518   fHMPIDp(0),
 519   fFriendTrack(0),
 520   fTPCClusterMap(159),//number of padrows
 521   fTPCSharedMap(159),//number of padrows
 522   fFlags(0),
 523   fID(0),
 524   fLabel(0),
 525   fITSLabel(0),
 526   fTPCLabel(0),
 527   fTRDLabel(0),
 528   fTOFCalChannel(-1),
 529   fTOFindex(-1),
 530   fHMPIDqn(0),
 531   fHMPIDcluIdx(-1),
 532   fCaloIndex(kEMCALNoMatch),
 533   fHMPIDtrkTheta(0),
 534   fHMPIDtrkPhi(0),
 535   fHMPIDsignal(0),
 536   fTrackLength(0),
 537   fdTPC(0),fzTPC(0),
 538   fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
 539   fCchi2TPC(0),
 540   fD(0),fZ(0),
 541   fCdd(0),fCdz(0),fCzz(0),
 542   fCchi2(0),
 543   fITSchi2(0),
 544   fTPCchi2(0),
 545   fTPCchi2Iter1(0),
 546   fTRDchi2(0),
 547   fTOFchi2(0),
 548   fHMPIDchi2(0),
 549   fGlobalChi2(0),
 550   fITSsignal(0),
 551   fTPCsignal(0),
 552   fTPCsignalS(0),
 553   fTRDsignal(0),
 554   fTRDQuality(0),
 555   fTRDBudget(0),
 556   fTOFsignal(99999),
 557   fTOFsignalToT(99999),
 558   fTOFsignalRaw(99999),
 559   fTOFsignalDz(999),
 560   fTOFsignalDx(999),
 561   fTOFdeltaBC(999),
 562   fTOFl0l1(999),
 563   fCaloDx(0),
 564   fCaloDz(0),
 565   fHMPIDtrkX(0),
 566   fHMPIDtrkY(0),
 567   fHMPIDmipX(0),
 568   fHMPIDmipY(0),
 569   fTPCncls(0),
 570   fTPCnclsF(0),
 571   fTPCsignalN(0),
 572   fTPCnclsIter1(0),
 573   fTPCnclsFIter1(0),
 574   fITSncls(0),
 575   fITSClusterMap(0),
 576   fITSSharedMap(0),
 577   fTRDncls(0),
 578   fTRDncls0(0),
 579   fTRDntracklets(0),
 580   fTRDnSlices(0),
 581   fTRDslices(0x0),
 582   fVertexID(-2),  // -2 means an orphan track
 583   fESDEvent(0)
 584 {
 585   //
 586   // ESD track from TParticle
 587   //
 588
 589   // Reset all the arrays
 590   Int_t i;
 591   for (i=0; i<AliPID::kSPECIES; i++) {
 592     fTrackTime[i]=0.;
 593     fR[i]=0.;
 594     fITSr[i]=0.;
 595     fTPCr[i]=0.;
 596     fTRDr[i]=0.;
 597     fTOFr[i]=0.;
 598     fHMPIDr[i]=0.;
 599   }
 600
 601   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
 602   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
 603   for (i=0;i<kTRDnPlanes;i++) {
 604     fTRDTimBin[i]=0;
 605   }
 606   for (i=0;i<4;i++) {fITSdEdxSamples[i]=0.;}
 607   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
 608   for (i=0;i<3;i++) {fTOFLabel[i]=-1;}
 609   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
 610   for (i=0;i<12;i++) {fITSModule[i]=-1;}
 611
 612   // Calculate the AliExternalTrackParam content
 613
 614   Double_t xref;
 615   Double_t alpha;
 616   Double_t param[5];
 617   Double_t covar[15];
 618
 619   // Calculate alpha: the rotation angle of the corresponding local system (TPC sector)
 620   alpha = part->Phi()*180./TMath::Pi();
 621   if (alpha<0) alpha+= 360.;
 622   if (alpha>360) alpha -= 360.;
 623
 624   Int_t sector = (Int_t)(alpha/20.);
 625   alpha = 10. + 20.*sector;
 626   alpha /= 180;
 627   alpha *= TMath::Pi();
 628
 629   // Covariance matrix: no errors, the parameters are exact
 630   for (i=0; i<15; i++) covar[i]=0.;
 631
 632   // Get the vertex of origin and the momentum
 633   TVector3 ver(part->Vx(),part->Vy(),part->Vz());
 634   TVector3 mom(part->Px(),part->Py(),part->Pz());
 635
 636   // Rotate to the local coordinate system (TPC sector)
 637   ver.RotateZ(-alpha);
 638   mom.RotateZ(-alpha);
 639
 640   // X of the referense plane
 641   xref = ver.X();
 642
 643   Int_t pdgCode = part->GetPdgCode();
 644
 645   Double_t charge =
 646     TDatabasePDG::Instance()->GetParticle(pdgCode)->Charge();
 647
 648   param[0] = ver.Y();
 649   param[1] = ver.Z();
 650   param[2] = TMath::Sin(mom.Phi());
 651   param[3] = mom.Pz()/mom.Pt();
 652   param[4] = TMath::Sign(1/mom.Pt(),charge);
 653
 654   // Set AliExternalTrackParam
 655   Set(xref, alpha, param, covar);
 656
 657   // Set the PID
 658   Int_t indexPID = 99;
 659
 660   switch (TMath::Abs(pdgCode)) {
 661
 662   case  11: // electron
 663     indexPID = 0;
 664     break;
 665
 666   case 13: // muon
 667     indexPID = 1;
 668     break;
 669
 670   case 211: // pion
 671     indexPID = 2;
 672     break;
 673
 674   case 321: // kaon
 675     indexPID = 3;
 676     break;
 677
 678   case 2212: // proton
 679     indexPID = 4;
 680     break;
 681
 682   default:
 683     break;
 684   }
 685
 686   // If the particle is not e,mu,pi,K or p the PID probabilities are set to 0
 687   if (indexPID < AliPID::kSPECIES) {
 688     fR[indexPID]=1.;
 689     fITSr[indexPID]=1.;
 690     fTPCr[indexPID]=1.;
 691     fTRDr[indexPID]=1.;
 692     fTOFr[indexPID]=1.;
 693     fHMPIDr[indexPID]=1.;
 694
 695   }
 696   // AliESD track label
 697   SetLabel(part->GetUniqueID());
 698
 699 }
 700
 701 //_______________________________________________________________________
 702 AliESDtrack::~AliESDtrack(){
 703   //
 704   // This is destructor according Coding Conventrions
 705   //
 706   //printf("Delete track\n");
 707   delete fIp;
 708   delete fTPCInner;
 709   delete fOp;
 710   delete fHMPIDp;
 711   delete fCp;
 712   if (fFriendTrack) delete fFriendTrack;
 713   fFriendTrack=NULL;
 714   if(fTRDnSlices)
 715     delete[] fTRDslices;
 716 }
 717
 718 AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
 719
 720
 721   if(&source == this) return *this;
 722   AliExternalTrackParam::operator=(source);
 723
 724
 725   if(source.fCp){
 726     // we have the trackparam: assign or copy construct
 727     if(fCp)*fCp = *source.fCp;
 728     else fCp = new AliExternalTrackParam(*source.fCp);
 729   }
 730   else{
 731     // no track param delete the old one
 732     if(fCp)delete fCp;
 733     fCp = 0;
 734   }
 735
 736   if(source.fIp){
 737     // we have the trackparam: assign or copy construct
 738     if(fIp)*fIp = *source.fIp;
 739     else fIp = new AliExternalTrackParam(*source.fIp);
 740   }
 741   else{
 742     // no track param delete the old one
 743     if(fIp)delete fIp;
 744     fIp = 0;
 745   }
 746
 747
 748   if(source.fTPCInner){
 749     // we have the trackparam: assign or copy construct
 750     if(fTPCInner) *fTPCInner = *source.fTPCInner;
 751     else fTPCInner = new AliExternalTrackParam(*source.fTPCInner);
 752   }
 753   else{
 754     // no track param delete the old one
 755     if(fTPCInner)delete fTPCInner;
 756     fTPCInner = 0;
 757   }
 758
 759
 760   if(source.fOp){
 761     // we have the trackparam: assign or copy construct
 762     if(fOp) *fOp = *source.fOp;
 763     else fOp = new AliExternalTrackParam(*source.fOp);
 764   }
 765   else{
 766     // no track param delete the old one
 767     if(fOp)delete fOp;
 768     fOp = 0;
 769   }
 770
 771
 772   if(source.fHMPIDp){
 773     // we have the trackparam: assign or copy construct
 774     if(fHMPIDp) *fHMPIDp = *source.fHMPIDp;
 775     else fHMPIDp = new AliExternalTrackParam(*source.fHMPIDp);
 776   }
 777   else{
 778     // no track param delete the old one
 779     if(fHMPIDp)delete fHMPIDp;
 780     fHMPIDp = 0;
 781   }
 782
 783
 784   // copy also the friend track
 785   // use copy constructor
 786   if(source.fFriendTrack){
 787     // we have the trackparam: assign or copy construct
 788     delete fFriendTrack; fFriendTrack=new AliESDfriendTrack(*source.fFriendTrack);
 789   }
 790   else{
 791     // no track param delete the old one
 792     delete fFriendTrack; fFriendTrack= 0;
 793   }
 794
 795   fTPCClusterMap = source.fTPCClusterMap;
 796   fTPCSharedMap  = source.fTPCSharedMap;
 797   // the simple stuff
 798   fFlags    = source.fFlags;
 799   fID       = source.fID;
 800   fLabel    = source.fLabel;
 801   fITSLabel = source.fITSLabel;
 802   for(int i = 0; i< 12;++i){
 803     fITSModule[i] = source.fITSModule[i];
 804   }
 805   fTPCLabel = source.fTPCLabel;
 806   fTRDLabel = source.fTRDLabel;
 807   for(int i = 0; i< 3;++i){
 808     fTOFLabel[i] = source.fTOFLabel[i];
 809   }
 810   fTOFCalChannel = source.fTOFCalChannel;
 811   fTOFindex      = source.fTOFindex;
 812   fHMPIDqn       = source.fHMPIDqn;
 813   fHMPIDcluIdx   = source.fHMPIDcluIdx;
 814   fCaloIndex    = source.fCaloIndex;
 815
 816   for(int i = 0; i< 3;++i){
 817     fKinkIndexes[i] = source.fKinkIndexes[i];
 818     fV0Indexes[i]   = source.fV0Indexes[i];
 819   }
 820
 821   for(int i = 0; i< AliPID::kSPECIES;++i){
 822     fR[i]     = source.fR[i];
 823     fITSr[i]  = source.fITSr[i];
 824     fTPCr[i]  = source.fTPCr[i];
 825     fTRDr[i]  = source.fTRDr[i];
 826     fTOFr[i]  = source.fTOFr[i];
 827     fHMPIDr[i] = source.fHMPIDr[i];
 828     fTrackTime[i] = source.fTrackTime[i];
 829   }
 830
 831   fHMPIDtrkTheta = source.fHMPIDtrkTheta;
 832   fHMPIDtrkPhi   = source.fHMPIDtrkPhi;
 833   fHMPIDsignal   = source.fHMPIDsignal;
 834
 835
 836   fTrackLength   = source. fTrackLength;
 837   fdTPC  = source.fdTPC;
 838   fzTPC  = source.fzTPC;
 839   fCddTPC = source.fCddTPC;
 840   fCdzTPC = source.fCdzTPC;
 841   fCzzTPC = source.fCzzTPC;
 842   fCchi2TPC = source.fCchi2TPC;
 843
 844   fD  = source.fD;
 845   fZ  = source.fZ;
 846   fCdd = source.fCdd;
 847   fCdz = source.fCdz;
 848   fCzz = source.fCzz;
 849   fCchi2     = source.fCchi2;
 850
 851   fITSchi2   = source.fITSchi2;
 852   fTPCchi2   = source.fTPCchi2;
 853   fTPCchi2Iter1   = source.fTPCchi2Iter1;
 854   fTRDchi2   = source.fTRDchi2;
 855   fTOFchi2   = source.fTOFchi2;
 856   fHMPIDchi2 = source.fHMPIDchi2;
 857
 858   fGlobalChi2 = source.fGlobalChi2;
 859
 860   fITSsignal  = source.fITSsignal;
 861   for (Int_t i=0;i<4;i++) {fITSdEdxSamples[i]=source.fITSdEdxSamples[i];}
 862   fTPCsignal  = source.fTPCsignal;
 863   fTPCsignalS = source.fTPCsignalS;
 864   for(int i = 0; i< 4;++i){
 865     fTPCPoints[i] = source.fTPCPoints[i];
 866   }
 867   fTRDsignal = source.fTRDsignal;
 868
 869   for(int i = 0;i < kTRDnPlanes;++i){
 870     fTRDTimBin[i] = source.fTRDTimBin[i];
 871   }
 872
 873   if(fTRDnSlices)
 874     delete[] fTRDslices;
 875   fTRDslices=0;
 876   fTRDnSlices=source.fTRDnSlices;
 877   if (fTRDnSlices) {
 878     fTRDslices=new Double32_t[fTRDnSlices];
 879     for(int j = 0;j < fTRDnSlices;++j) fTRDslices[j] = source.fTRDslices[j];
 880   }
 881
 882   fTRDQuality =   source.fTRDQuality;
 883   fTRDBudget  =   source.fTRDBudget;
 884   fTOFsignal  =   source.fTOFsignal;
 885   fTOFsignalToT = source.fTOFsignalToT;
 886   fTOFsignalRaw = source.fTOFsignalRaw;
 887   fTOFsignalDz  = source.fTOFsignalDz;
 888   fTOFsignalDx  = source.fTOFsignalDx;
 889   fTOFdeltaBC   = source.fTOFdeltaBC;
 890   fTOFl0l1      = source.fTOFl0l1;
 891
 892   for(int i = 0;i<10;++i){
 893     fTOFInfo[i] = source.fTOFInfo[i];
 894   }
 895
 896   fHMPIDtrkX = source.fHMPIDtrkX;
 897   fHMPIDtrkY = source.fHMPIDtrkY;
 898   fHMPIDmipX = source.fHMPIDmipX;
 899   fHMPIDmipY = source.fHMPIDmipY;
 900
 901   fTPCncls    = source.fTPCncls;
 902   fTPCnclsF   = source.fTPCnclsF;
 903   fTPCsignalN = source.fTPCsignalN;
 904   fTPCnclsIter1    = source.fTPCnclsIter1;
 905   fTPCnclsFIter1   = source.fTPCnclsFIter1;
 906
 907   fITSncls = source.fITSncls;
 908   fITSClusterMap = source.fITSClusterMap;
 909   fITSSharedMap = source.fITSSharedMap;
 910   fTRDncls   = source.fTRDncls;
 911   fTRDncls0  = source.fTRDncls0;
 912   fTRDntracklets  = source.fTRDntracklets;
 913   fVertexID = source.fVertexID;
 914   return *this;
 915 }
 916
 917
 918
 919 void AliESDtrack::Copy(TObject &obj) const {
 920
 921   // this overwrites the virtual TOBject::Copy()
 922   // to allow run time copying without casting
 923   // in AliESDEvent
 924
 925   if(this==&obj)return;
 926   AliESDtrack *robj = dynamic_cast<AliESDtrack*>(&obj);
 927   if(!robj)return; // not an AliESDtrack
 928   *robj = *this;
 929
 930 }
 931
 932
 933
 934 void AliESDtrack::AddCalibObject(TObject * object){
 935   //
 936   // add calib object to the list
 937   //
 938   if (!fFriendTrack) fFriendTrack  = new AliESDfriendTrack;
 939   if (!fFriendTrack) return;
 940   fFriendTrack->AddCalibObject(object);
 941 }
 942
 943 TObject *  AliESDtrack::GetCalibObject(Int_t index){
 944   //
 945   // return calib objct at given position
 946   //
 947   if (!fFriendTrack) return 0;
 948   return fFriendTrack->GetCalibObject(index);
 949 }
 950
 951
 952 Bool_t AliESDtrack::FillTPCOnlyTrack(AliESDtrack &track){
 953
 954   // Fills the information of the TPC-only first reconstruction pass
 955   // into the passed ESDtrack object. For consistency fTPCInner is also filled
 956   // again
 957
 958
 959
 960   // For data produced before r26675
 961   // RelateToVertexTPC was not properly called during reco
 962   // so you'll have to call it again, before FillTPCOnlyTrack
 963   //  Float_t p[2],cov[3];
 964   // track->GetImpactParametersTPC(p,cov);
 965   // if(p[0]==0&&p[1]==0) // <- Default values
 966   //  track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
 967
 968
 969   if(!fTPCInner)return kFALSE;
 970
 971   // fill the TPC track params to the global track parameters
 972   track.Set(fTPCInner->GetX(),fTPCInner->GetAlpha(),fTPCInner->GetParameter(),fTPCInner->GetCovariance());
 973   track.fD = fdTPC;
 974   track.fZ = fzTPC;
 975   track.fCdd = fCddTPC;
 976   track.fCdz = fCdzTPC;
 977   track.fCzz = fCzzTPC;
 978
 979   // copy the inner params
 980   if(track.fIp) *track.fIp = *fIp;
 981   else track.fIp = new AliExternalTrackParam(*fIp);
 982
 983   // copy the TPCinner parameters
 984   if(track.fTPCInner) *track.fTPCInner = *fTPCInner;
 985   else track.fTPCInner = new AliExternalTrackParam(*fTPCInner);
 986   track.fdTPC   = fdTPC;
 987   track.fzTPC   = fzTPC;
 988   track.fCddTPC = fCddTPC;
 989   track.fCdzTPC = fCdzTPC;
 990   track.fCzzTPC = fCzzTPC;
 991   track.fCchi2TPC = fCchi2TPC;
 992
 993   // copy all other TPC specific parameters
 994
 995   // replace label by TPC label
 996   track.fLabel    = fTPCLabel;
 997   track.fTPCLabel = fTPCLabel;
 998
 999   track.fTPCchi2 = fTPCchi2;
1000   track.fTPCchi2Iter1 = fTPCchi2Iter1;
1001   track.fTPCsignal = fTPCsignal;
1002   track.fTPCsignalS = fTPCsignalS;
1003   for(int i = 0;i<4;++i)track.fTPCPoints[i] = fTPCPoints[i];
1004
1005   track.fTPCncls    = fTPCncls;
1006   track.fTPCnclsF   = fTPCnclsF;
1007   track.fTPCsignalN =  fTPCsignalN;
1008   track.fTPCnclsIter1    = fTPCnclsIter1;
1009   track.fTPCnclsFIter1   = fTPCnclsFIter1;
1010
1011   // PID
1012   for(int i=0;i<AliPID::kSPECIES;++i){
1013     track.fTPCr[i] = fTPCr[i];
1014     // combined PID is TPC only!
1015     track.fR[i] = fTPCr[i];
1016   }
1017   track.fTPCClusterMap = fTPCClusterMap;
1018   track.fTPCSharedMap = fTPCSharedMap;
1019
1020
1021   // reset the flags
1022   track.fFlags = kTPCin;
1023   track.fID    = fID;
1024
1025   track.fFlags |= fFlags & kTPCpid; //copy the TPCpid status flag
1026
1027   for (Int_t i=0;i<3;i++) track.fKinkIndexes[i] = fKinkIndexes[i];
1028
1029   return kTRUE;
1030
1031 }
1032
1033 //_______________________________________________________________________
1034 void AliESDtrack::MakeMiniESDtrack(){
1035   // Resets everything except
1036   // fFlags: Reconstruction status flags
1037   // fLabel: Track label
1038   // fID:  Unique ID of the track
1039   // Impact parameter information
1040   // fR[AliPID::kSPECIES]: combined "detector response probability"
1041   // Running track parameters in the base class (AliExternalTrackParam)
1042
1043   fTrackLength = 0;
1044
1045   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i] = 0;
1046
1047   // Reset track parameters constrained to the primary vertex
1048   delete fCp;fCp = 0;
1049
1050   // Reset track parameters at the inner wall of TPC
1051   delete fIp;fIp = 0;
1052   delete fTPCInner;fTPCInner=0;
1053   // Reset track parameters at the inner wall of the TRD
1054   delete fOp;fOp = 0;
1055   // Reset track parameters at the HMPID
1056   delete fHMPIDp;fHMPIDp = 0;
1057
1058
1059   // Reset ITS track related information
1060   fITSchi2 = 0;
1061   fITSncls = 0;
1062   fITSClusterMap=0;
1063   fITSSharedMap=0;
1064   fITSsignal = 0;
1065   for (Int_t i=0;i<4;i++) fITSdEdxSamples[i] = 0.;
1066   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=0;
1067   fITSLabel = 0;
1068
1069   // Reset TPC related track information
1070   fTPCchi2 = 0;
1071   fTPCchi2Iter1 = 0;
1072   fTPCncls = 0;
1073   fTPCnclsF = 0;
1074   fTPCnclsIter1 = 0;
1075   fTPCnclsFIter1 = 0;
1076   fTPCClusterMap = 0;
1077   fTPCSharedMap = 0;
1078   fTPCsignal= 0;
1079   fTPCsignalS= 0;
1080   fTPCsignalN= 0;
1081   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0;
1082   fTPCLabel=0;
1083   for (Int_t i=0;i<4;i++) fTPCPoints[i] = 0;
1084   for (Int_t i=0; i<3;i++)   fKinkIndexes[i] = 0;
1085   for (Int_t i=0; i<3;i++)   fV0Indexes[i] = 0;
1086
1087   // Reset TRD related track information
1088   fTRDchi2 = 0;
1089   fTRDncls = 0;
1090   fTRDncls0 = 0;
1091   fTRDsignal = 0;
1092   for (Int_t i=0;i<kTRDnPlanes;i++) {
1093     fTRDTimBin[i]  = 0;
1094   }
1095   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0;
1096   fTRDLabel = 0;
1097   fTRDQuality  = 0;
1098   fTRDntracklets = 0;
1099   if(fTRDnSlices)
1100     delete[] fTRDslices;
1101   fTRDslices=0x0;
1102   fTRDnSlices=0;
1103   fTRDBudget  = 0;
1104
1105   // Reset TOF related track information
1106   fTOFchi2 = 0;
1107   fTOFindex = -1;
1108   fTOFsignal = 99999;
1109   fTOFCalChannel = -1;
1110   fTOFsignalToT = 99999;
1111   fTOFsignalRaw = 99999;
1112   fTOFsignalDz = 999;
1113   fTOFsignalDx = 999;
1114   fTOFdeltaBC = 999;
1115   fTOFl0l1 = 999;
1116   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i] = 0;
1117   for (Int_t i=0;i<3;i++) fTOFLabel[i] = -1;
1118   for (Int_t i=0;i<10;i++) fTOFInfo[i] = 0;
1119
1120   // Reset HMPID related track information
1121   fHMPIDchi2 = 0;
1122   fHMPIDqn = 0;
1123   fHMPIDcluIdx = -1;
1124   fHMPIDsignal = 0;
1125   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i] = 0;
1126   fHMPIDtrkTheta = 0;
1127   fHMPIDtrkPhi = 0;
1128   fHMPIDtrkX = 0;
1129   fHMPIDtrkY = 0;
1130   fHMPIDmipX = 0;
1131   fHMPIDmipY = 0;
1132   fCaloIndex = kEMCALNoMatch;
1133
1134   // reset global track chi2
1135   fGlobalChi2 = 0;
1136
1137   fVertexID = -2; // an orphan track
1138
1139   delete fFriendTrack; fFriendTrack = 0;
1140 }
1141
1142 //_______________________________________________________________________
1143 Int_t AliESDtrack::GetPID() const
1144 {
1145   // Returns the particle most probable id
1146   Int_t i;
1147   for (i=0; i<AliPID::kSPECIES-1; i++) if (fR[i] != fR[i+1]) break;
1148   //
1149   if (i == AliPID::kSPECIES-1) return AliPID::kPion;  // If all the probabilities are equal, return the pion mass
1150   //
1151   Float_t max=0.;
1152   Int_t k=-1;
1153   for (i=0; i<AliPID::kSPECIES; i++) if (fR[i]>max) {k=i; max=fR[i];}
1154   //
1155   if (k==0) { // dE/dx "crossing points" in the TPC
1156     Double_t p=GetP();
1157     if ((p>0.38)&&(p<0.48))
1158       if (fR[0]<fR[3]*10.) return AliPID::kKaon;
1159     if ((p>0.75)&&(p<0.85))
1160       if (fR[0]<fR[4]*10.) return AliPID::kProton;
1161     return AliPID::kElectron;
1162   }
1163   if (k==1) return AliPID::kMuon;
1164   if (k==2||k==-1) return AliPID::kPion;
1165   if (k==3) return AliPID::kKaon;
1166   if (k==4) return AliPID::kProton;
1167   AliWarning("Undefined PID !");
1168   return AliPID::kPion;
1169 }
1170
1171 //_______________________________________________________________________
1172 Int_t AliESDtrack::GetTOFBunchCrossing() const
1173 {
1174   // Returns the number of bunch crossings after trigger (assuming 25ns spacing)
1175   const UInt_t kAskBits = kESDpid | kTOFout | kTIME;
1176   const double kSpacing = 25e3; // min interbanch spacing
1177   const double kShift = 0;
1178   Int_t bcid = -1; // defualt one
1179   if ( (GetStatus()&kAskBits) != kAskBits) return bcid;
1180   int pid = GetPID();
1181   double tdif = fTOFsignal - fTrackTime[pid];
1182   bcid = TMath::Nint((tdif - kShift)/kSpacing);
1183   return bcid;
1184 }
1185
1186 //______________________________________________________________________________
1187 Double_t AliESDtrack::M() const
1188 {
1189   // Returns the assumed mass
1190   // (the pion mass, if the particle can't be identified properly).
1191   static Bool_t printerr=kTRUE;
1192   if (printerr) {
1193      AliWarning("WARNING !!! ... THIS WILL BE PRINTED JUST ONCE !!!");
1194      printerr = kFALSE;
1195      AliWarning("This is the ESD mass. Use it with care !");
1196   }
1197   return GetMass();
1198 }
1199
1200 //______________________________________________________________________________
1201 Double_t AliESDtrack::E() const
1202 {
1203   // Returns the energy of the particle given its assumed mass.
1204   // Assumes the pion mass if the particle can't be identified properly.
1205
1206   Double_t m = M();
1207   Double_t p = P();
1208   return TMath::Sqrt(p*p + m*m);
1209 }
1210
1211 //______________________________________________________________________________
1212 Double_t AliESDtrack::Y() const
1213 {
1214   // Returns the rapidity of a particle given its assumed mass.
1215   // Assumes the pion mass if the particle can't be identified properly.
1216
1217   Double_t e = E();
1218   Double_t pz = Pz();
1219   if (e != TMath::Abs(pz)) { // energy was not equal to pz
1220     return 0.5*TMath::Log((e+pz)/(e-pz));
1221   } else { // energy was equal to pz
1222     return -999.;
1223   }
1224 }
1225
1226 //_______________________________________________________________________
1227 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
1228   //
1229   // This function updates track's running parameters
1230   //
1231   Bool_t rc=kTRUE;
1232
1233   SetStatus(flags);
1234   fLabel=t->GetLabel();
1235
1236   if (t->IsStartedTimeIntegral()) {
1237     SetStatus(kTIME);
1238     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
1239     SetIntegratedLength(t->GetIntegratedLength());
1240   }
1241
1242   Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1243   if (fFriendTrack) {
1244   if (flags==kITSout) fFriendTrack->SetITSOut(*t);
1245   if (flags==kTPCout) fFriendTrack->SetTPCOut(*t);
1246   if (flags==kTRDrefit) fFriendTrack->SetTRDIn(*t);
1247   }
1248
1249   switch (flags) {
1250
1251   case kITSin: case kITSout: case kITSrefit:
1252     {
1253     fITSClusterMap=0;
1254     fITSncls=t->GetNumberOfClusters();
1255     if (fFriendTrack) {
1256     Int_t* indexITS = new Int_t[AliESDfriendTrack::kMaxITScluster];
1257     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
1258         indexITS[i]=t->GetClusterIndex(i);
1259
1260         if (i<fITSncls) {
1261           Int_t l=(indexITS[i] & 0xf0000000) >> 28;
1262            SETBIT(fITSClusterMap,l);
1263         }
1264     }
1265     fFriendTrack->SetITSIndices(indexITS,AliESDfriendTrack::kMaxITScluster);
1266     delete [] indexITS;
1267     }
1268
1269     fITSchi2=t->GetChi2();
1270     fITSsignal=t->GetPIDsignal();
1271     fITSLabel = t->GetLabel();
1272     // keep in fOp the parameters outside ITS for ITS stand-alone tracks
1273     if (flags==kITSout) {
1274       if (!fOp) fOp=new AliExternalTrackParam(*t);
1275       else
1276         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1277     }
1278     }
1279     break;
1280
1281   case kTPCin: case kTPCrefit:
1282     {
1283     fTPCLabel = t->GetLabel();
1284     if (flags==kTPCin)  {
1285       fTPCInner=new AliExternalTrackParam(*t);
1286       fTPCnclsIter1=t->GetNumberOfClusters();
1287       fTPCchi2Iter1=t->GetChi2();
1288     }
1289     if (!fIp) fIp=new AliExternalTrackParam(*t);
1290     else
1291       fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1292     }
1293   case kTPCout:
1294     {
1295     if (flags & kTPCout){
1296       if (!fOp) fOp=new AliExternalTrackParam(*t);
1297       else
1298         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1299     }
1300     fTPCncls=t->GetNumberOfClusters();
1301     fTPCchi2=t->GetChi2();
1302
1303     if (fFriendTrack)
1304      {//prevrow must be declared in separate namespace, otherwise compiler cries:
1305       //"jump to case label crosses initialization of `Int_t prevrow'"
1306        Int_t* indexTPC = new Int_t[AliESDfriendTrack::kMaxTPCcluster];
1307        Int_t prevrow = -1;
1308        //       for (Int_t i=0;i<fTPCncls;i++)
1309        for (Int_t i=0;i<AliESDfriendTrack::kMaxTPCcluster;i++)
1310         {
1311           indexTPC[i]=t->GetClusterIndex(i);
1312           Int_t idx = indexTPC[i];
1313
1314           if (idx<0) continue;
1315
1316           // Piotr's Cluster Map for HBT
1317           // ### please change accordingly if cluster array is changing
1318           // to "New TPC Tracking" style (with gaps in array)
1319           Int_t sect = (idx&0xff000000)>>24;
1320           Int_t row = (idx&0x00ff0000)>>16;
1321           if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
1322
1323           fTPCClusterMap.SetBitNumber(row,kTRUE);
1324
1325           //Fill the gap between previous row and this row with 0 bits
1326           //In case  ###  pleas change it as well - just set bit 0 in case there
1327           //is no associated clusters for current "i"
1328           if (prevrow < 0)
1329            {
1330              prevrow = row;//if previous bit was not assigned yet == this is the first one
1331            }
1332           else
1333            { //we don't know the order (inner to outer or reverse)
1334              //just to be save in case it is going to change
1335              Int_t n = 0, m = 0;
1336              if (prevrow < row)
1337               {
1338                 n = prevrow;
1339                 m = row;
1340               }
1341              else
1342               {
1343                 n = row;
1344                 m = prevrow;
1345               }
1346
1347              for (Int_t j = n+1; j < m; j++)
1348               {
1349                 fTPCClusterMap.SetBitNumber(j,kFALSE);
1350               }
1351              prevrow = row;
1352            }
1353           // End Of Piotr's Cluster Map for HBT
1354         }
1355         fFriendTrack->SetTPCIndices(indexTPC,AliESDfriendTrack::kMaxTPCcluster);
1356         delete [] indexTPC;
1357
1358      }
1359     fTPCsignal=t->GetPIDsignal();
1360     }
1361     break;
1362
1363   case kTRDin: case kTRDrefit:
1364     break;
1365   case kTRDout:
1366     {
1367     fTRDLabel = t->GetLabel();
1368     fTRDchi2  = t->GetChi2();
1369     fTRDncls  = t->GetNumberOfClusters();
1370     if (fFriendTrack) {
1371       Int_t* indexTRD = new Int_t[AliESDfriendTrack::kMaxTRDcluster];
1372       for (Int_t i=0;i<AliESDfriendTrack::kMaxTRDcluster;i++) indexTRD[i]=-2;
1373       for (Int_t i=0;i<6;i++) indexTRD[i]=t->GetTrackletIndex(i);
1374       fFriendTrack->SetTRDIndices(indexTRD,AliESDfriendTrack::kMaxTRDcluster);
1375       delete [] indexTRD;
1376     }
1377
1378     fTRDsignal=t->GetPIDsignal();
1379     }
1380     break;
1381   case kTRDbackup:
1382     if (!fOp) fOp=new AliExternalTrackParam(*t);
1383     else
1384       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1385     fTRDncls0 = t->GetNumberOfClusters();
1386     break;
1387   case kTOFin:
1388     break;
1389   case kTOFout:
1390     break;
1391   case kTRDStop:
1392     break;
1393   case kHMPIDout:
1394   if (!fHMPIDp) fHMPIDp=new AliExternalTrackParam(*t);
1395     else
1396       fHMPIDp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1397     break;
1398   default:
1399     AliError("Wrong flag !");
1400     return kFALSE;
1401   }
1402
1403   return rc;
1404 }
1405
1406 //_______________________________________________________________________
1407 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
1408   //---------------------------------------------------------------------
1409   // This function returns external representation of the track parameters
1410   //---------------------------------------------------------------------
1411   x=GetX();
1412   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
1413 }
1414
1415 //_______________________________________________________________________
1416 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
1417   //---------------------------------------------------------------------
1418   // This function returns external representation of the cov. matrix
1419   //---------------------------------------------------------------------
1420   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
1421 }
1422
1423 //_______________________________________________________________________
1424 Bool_t AliESDtrack::GetConstrainedExternalParameters
1425                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1426   //---------------------------------------------------------------------
1427   // This function returns the constrained external track parameters
1428   //---------------------------------------------------------------------
1429   if (!fCp) return kFALSE;
1430   alpha=fCp->GetAlpha();
1431   x=fCp->GetX();
1432   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
1433   return kTRUE;
1434 }
1435
1436 //_______________________________________________________________________
1437 Bool_t
1438 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
1439   //---------------------------------------------------------------------
1440   // This function returns the constrained external cov. matrix
1441   //---------------------------------------------------------------------
1442   if (!fCp) return kFALSE;
1443   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
1444   return kTRUE;
1445 }
1446
1447 Bool_t
1448 AliESDtrack::GetInnerExternalParameters
1449                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1450   //---------------------------------------------------------------------
1451   // This function returns external representation of the track parameters
1452   // at the inner layer of TPC
1453   //---------------------------------------------------------------------
1454   if (!fIp) return kFALSE;
1455   alpha=fIp->GetAlpha();
1456   x=fIp->GetX();
1457   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
1458   return kTRUE;
1459 }
1460
1461 Bool_t
1462 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
1463  //---------------------------------------------------------------------
1464  // This function returns external representation of the cov. matrix
1465  // at the inner layer of TPC
1466  //---------------------------------------------------------------------
1467   if (!fIp) return kFALSE;
1468   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
1469   return kTRUE;
1470 }
1471
1472 void
1473 AliESDtrack::SetOuterParam(const AliExternalTrackParam *p, ULong_t flags) {
1474   //
1475   // This is a direct setter for the outer track parameters
1476   //
1477   SetStatus(flags);
1478   if (fOp) delete fOp;
1479   fOp=new AliExternalTrackParam(*p);
1480 }
1481
1482 void
1483 AliESDtrack::SetOuterHmpParam(const AliExternalTrackParam *p, ULong_t flags) {
1484   //
1485   // This is a direct setter for the outer track parameters
1486   //
1487   SetStatus(flags);
1488   if (fHMPIDp) delete fHMPIDp;
1489   fHMPIDp=new AliExternalTrackParam(*p);
1490 }
1491
1492 Bool_t
1493 AliESDtrack::GetOuterExternalParameters
1494                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1495   //---------------------------------------------------------------------
1496   // This function returns external representation of the track parameters
1497   // at the inner layer of TRD
1498   //---------------------------------------------------------------------
1499   if (!fOp) return kFALSE;
1500   alpha=fOp->GetAlpha();
1501   x=fOp->GetX();
1502   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
1503   return kTRUE;
1504 }
1505
1506 Bool_t
1507 AliESDtrack::GetOuterHmpExternalParameters
1508                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1509   //---------------------------------------------------------------------
1510   // This function returns external representation of the track parameters
1511   // at the inner layer of TRD
1512   //---------------------------------------------------------------------
1513   if (!fHMPIDp) return kFALSE;
1514   alpha=fHMPIDp->GetAlpha();
1515   x=fHMPIDp->GetX();
1516   for (Int_t i=0; i<5; i++) p[i]=fHMPIDp->GetParameter()[i];
1517   return kTRUE;
1518 }
1519
1520 Bool_t
1521 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
1522  //---------------------------------------------------------------------
1523  // This function returns external representation of the cov. matrix
1524  // at the inner layer of TRD
1525  //---------------------------------------------------------------------
1526   if (!fOp) return kFALSE;
1527   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
1528   return kTRUE;
1529 }
1530
1531 Bool_t
1532 AliESDtrack::GetOuterHmpExternalCovariance(Double_t cov[15]) const {
1533  //---------------------------------------------------------------------
1534  // This function returns external representation of the cov. matrix
1535  // at the inner layer of TRD
1536  //---------------------------------------------------------------------
1537   if (!fHMPIDp) return kFALSE;
1538   for (Int_t i=0; i<15; i++) cov[i]=fHMPIDp->GetCovariance()[i];
1539   return kTRUE;
1540 }
1541
1542 Int_t AliESDtrack::GetNcls(Int_t idet) const
1543 {
1544   // Get number of clusters by subdetector index
1545   //
1546   Int_t ncls = 0;
1547   switch(idet){
1548   case 0:
1549     ncls = fITSncls;
1550     break;
1551   case 1:
1552     ncls = fTPCncls;
1553     break;
1554   case 2:
1555     ncls = fTRDncls;
1556     break;
1557   case 3:
1558     if (fTOFindex != -1)
1559       ncls = 1;
1560     break;
1561   case 4: //PHOS
1562     break;
1563   case 5: //HMPID
1564     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1565       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1566         ncls = 1;
1567       }
1568     }
1569     break;
1570   default:
1571     break;
1572   }
1573   return ncls;
1574 }
1575
1576 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
1577 {
1578   // Get cluster index array by subdetector index
1579   //
1580   Int_t ncls = 0;
1581   switch(idet){
1582   case 0:
1583     ncls = GetITSclusters(idx);
1584     break;
1585   case 1:
1586     ncls = GetTPCclusters(idx);
1587     break;
1588   case 2:
1589     ncls = GetTRDclusters(idx);
1590     break;
1591   case 3:
1592     if (fTOFindex != -1) {
1593       idx[0] = fTOFindex;
1594       ncls = 1;
1595     }
1596     break;
1597   case 4: //PHOS
1598     break;
1599   case 5:
1600     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1601       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1602         idx[0] = GetHMPIDcluIdx();
1603         ncls = 1;
1604       }
1605     }
1606     break;
1607   case 6: //EMCAL
1608     break;
1609   default:
1610     break;
1611   }
1612   return ncls;
1613 }
1614
1615 //_______________________________________________________________________
1616 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
1617   // Returns the array with integrated times for each particle hypothesis
1618   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
1619 }
1620
1621 //_______________________________________________________________________
1622 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
1623   // Sets the array with integrated times for each particle hypotesis
1624   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
1625 }
1626
1627 //_______________________________________________________________________
1628 void AliESDtrack::SetITSpid(const Double_t *p) {
1629   // Sets values for the probability of each particle type (in ITS)
1630   SetPIDValues(fITSr,p,AliPID::kSPECIES);
1631   SetStatus(AliESDtrack::kITSpid);
1632 }
1633
1634 //_______________________________________________________________________
1635 void AliESDtrack::GetITSpid(Double_t *p) const {
1636   // Gets the probability of each particle type (in ITS)
1637   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
1638 }
1639
1640 //_______________________________________________________________________
1641 Char_t AliESDtrack::GetITSclusters(Int_t *idx) const {
1642   //---------------------------------------------------------------------
1643   // This function returns indices of the assgined ITS clusters
1644   //---------------------------------------------------------------------
1645   if (idx && fFriendTrack) {
1646     Int_t *index=fFriendTrack->GetITSindices();
1647     for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) {
1648       if ( (i>=fITSncls) && (i<6) ) idx[i]=-1;
1649       else {
1650         if (index) {
1651           idx[i]=index[i];
1652         }
1653         else idx[i]= -2;
1654       }
1655     }
1656   }
1657   return fITSncls;
1658 }
1659
1660 //_______________________________________________________________________
1661 Bool_t AliESDtrack::GetITSModuleIndexInfo(Int_t ilayer,Int_t &idet,Int_t &status,
1662                                          Float_t &xloc,Float_t &zloc) const {
1663   //----------------------------------------------------------------------
1664   // This function encodes in the module number also the status of cluster association
1665   // "status" can have the following values:
1666   // 1 "found" (cluster is associated),
1667   // 2 "dead" (module is dead from OCDB),
1668   // 3 "skipped" (module or layer forced to be skipped),
1669   // 4 "outinz" (track out of z acceptance),
1670   // 5 "nocls" (no clusters in the road),
1671   // 6 "norefit" (cluster rejected during refit),
1672   // 7 "deadzspd" (holes in z in SPD)
1673   // Also given are the coordinates of the crossing point of track and module
1674   // (in the local module ref. system)
1675   // WARNING: THIS METHOD HAS TO BE SYNCHRONIZED WITH AliITStrackV2::GetModuleIndexInfo()!
1676   //----------------------------------------------------------------------
1677
1678   if(fITSModule[ilayer]==-1) {
1679     idet = -1;
1680     status=0;
1681     xloc=-99.; zloc=-99.;
1682     return kFALSE;
1683   }
1684
1685   Int_t module = fITSModule[ilayer];
1686
1687   idet = Int_t(module/1000000);
1688
1689   module -= idet*1000000;
1690
1691   status = Int_t(module/100000);
1692
1693   module -= status*100000;
1694
1695   Int_t signs = Int_t(module/10000);
1696
1697   module-=signs*10000;
1698
1699   Int_t xInt = Int_t(module/100);
1700   module -= xInt*100;
1701
1702   Int_t zInt = module;
1703
1704   if(signs==1) { xInt*=1; zInt*=1; }
1705   if(signs==2) { xInt*=1; zInt*=-1; }
1706   if(signs==3) { xInt*=-1; zInt*=1; }
1707   if(signs==4) { xInt*=-1; zInt*=-1; }
1708
1709   xloc = 0.1*(Float_t)xInt;
1710   zloc = 0.1*(Float_t)zInt;
1711
1712   if(status==4) idet = -1;
1713
1714   return kTRUE;
1715 }
1716
1717 //_______________________________________________________________________
1718 UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
1719   //---------------------------------------------------------------------
1720   // This function returns indices of the assgined ITS clusters
1721   //---------------------------------------------------------------------
1722   if (idx && fFriendTrack) {
1723     Int_t *index=fFriendTrack->GetTPCindices();
1724
1725     if (index){
1726       for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
1727     }
1728     else {
1729       for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=-2;
1730     }
1731   }
1732   return fTPCncls;
1733 }
1734
1735 //_______________________________________________________________________
1736 Float_t AliESDtrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1737 {
1738   //
1739   // TPC cluster information
1740   // type 0: get fraction of found/findable clusters with neighbourhood definition
1741   //      1: findable clusters with neighbourhood definition
1742   //      2: found clusters
1743   //
1744   // definition of findable clusters:
1745   //            a cluster is defined as findable if there is another cluster
1746   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1747   //           effects with a very simple algorithm.
1748   //
1749
1750   if (type==2) return fTPCClusterMap.CountBits();
1751
1752   Int_t found=0;
1753   Int_t findable=0;
1754   Int_t last=-nNeighbours;
1755
1756   Int_t upperBound=fTPCClusterMap.GetNbits();
1757   if (upperBound>row1) upperBound=row1;
1758   for (Int_t i=row0; i<upperBound; ++i){
1759     //look to current row
1760     if (fTPCClusterMap[i]) {
1761       last=i;
1762       ++found;
1763       ++findable;
1764       continue;
1765     }
1766     //look to nNeighbours before
1767     if ((i-last)<=nNeighbours) {
1768       ++findable;
1769       continue;
1770     }
1771     //look to nNeighbours after
1772     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1773       if (fTPCClusterMap[j]){
1774         ++findable;
1775         break;
1776       }
1777     }
1778   }
1779   if (type==1) return findable;
1780
1781   if (type==0){
1782     Float_t fraction=0;
1783     if (findable>0)
1784       fraction=(Float_t)found/(Float_t)findable;
1785     else
1786       fraction=0;
1787     return fraction;
1788   }
1789   return 0;  // undefined type - default value
1790 }
1791
1792 //_______________________________________________________________________
1793 Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
1794   //
1795   // GetDensity of the clusters on given region between row0 and row1
1796   // Dead zone effect takin into acoount
1797   //
1798   if (!fFriendTrack) return 0.0;
1799   Int_t good  = 0;
1800   Int_t found = 0;
1801   //
1802   Int_t *index=fFriendTrack->GetTPCindices();
1803   for (Int_t i=row0;i<=row1;i++){
1804     Int_t idx = index[i];
1805     if (idx!=-1)  good++;             // track outside of dead zone
1806     if (idx>0)    found++;
1807   }
1808   Float_t density=0.5;
1809   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
1810   return density;
1811 }
1812
1813 //_______________________________________________________________________
1814 void AliESDtrack::SetTPCpid(const Double_t *p) {
1815   // Sets values for the probability of each particle type (in TPC)
1816   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
1817   SetStatus(AliESDtrack::kTPCpid);
1818 }
1819
1820 //_______________________________________________________________________
1821 void AliESDtrack::GetTPCpid(Double_t *p) const {
1822   // Gets the probability of each particle type (in TPC)
1823   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
1824 }
1825
1826 //_______________________________________________________________________
1827 UChar_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
1828   //---------------------------------------------------------------------
1829   // This function returns indices of the assgined TRD clusters
1830   //---------------------------------------------------------------------
1831   if (idx && fFriendTrack) {
1832     Int_t *index=fFriendTrack->GetTRDindices();
1833
1834     if (index) {
1835       for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
1836     }
1837     else {
1838       for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=-2;
1839     }
1840   }
1841   return fTRDncls;
1842 }
1843
1844 //_______________________________________________________________________
1845 UChar_t AliESDtrack::GetTRDtracklets(Int_t *idx) const {
1846 //
1847 // This function returns the number of TRD tracklets used in tracking
1848 // and it fills the indices of these tracklets in the array "idx" as they
1849 // are registered in the TRD track list.
1850 //
1851 // Caution :
1852 //   1. The idx array has to be allocated with a size >= AliESDtrack::kTRDnPlanes
1853 //   2. The idx array store not only the index but also the layer of the tracklet.
1854 //      Therefore tracks with TRD gaps contain default values for indices [-1]
1855
1856   if (!fFriendTrack) return 0;
1857   if (!idx) return GetTRDntracklets();
1858   Int_t *index=fFriendTrack->GetTRDindices();
1859   Int_t n = 0;
1860   for (Int_t i=0; i<kTRDnPlanes; i++){
1861     if (index){
1862       if(index[i]>=0) n++;
1863       idx[i]=index[i];
1864     }
1865     else idx[i] = -2;
1866   }
1867   return n;
1868 }
1869
1870 //_______________________________________________________________________
1871 void AliESDtrack::SetTRDpid(const Double_t *p) {
1872   // Sets values for the probability of each particle type (in TRD)
1873   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
1874   SetStatus(AliESDtrack::kTRDpid);
1875 }
1876
1877 //_______________________________________________________________________
1878 void AliESDtrack::GetTRDpid(Double_t *p) const {
1879   // Gets the probability of each particle type (in TRD)
1880   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
1881 }
1882
1883 //_______________________________________________________________________
1884 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
1885 {
1886   // Sets the probability of particle type iSpecies to p (in TRD)
1887   fTRDr[iSpecies] = p;
1888 }
1889
1890 Double_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
1891 {
1892   // Returns the probability of particle type iSpecies (in TRD)
1893   return fTRDr[iSpecies];
1894 }
1895
1896 //____________________________________________________
1897 Int_t AliESDtrack::GetNumberOfTRDslices() const
1898 {
1899   // built in backward compatibility
1900   Int_t idx = fTRDnSlices - (kTRDnPlanes<<1);
1901   return idx<18 ? fTRDnSlices/kTRDnPlanes : idx/kTRDnPlanes;
1902 }
1903
1904 //____________________________________________________
1905 Double_t AliESDtrack::GetTRDmomentum(Int_t plane, Double_t *sp) const
1906 {
1907 //Returns momentum estimation and optional its error (sp)
1908 // in TRD layer "plane".
1909
1910   if (!fTRDnSlices) {
1911     AliDebug(2, "No TRD info allocated for this track.");
1912     return -1.;
1913   }
1914   if ((plane<0) || (plane>=kTRDnPlanes)) {
1915     AliWarning(Form("Request for TRD plane[%d] outside range.", plane));
1916     return -1.;
1917   }
1918
1919   Int_t idx = fTRDnSlices-(kTRDnPlanes<<1)+plane;
1920   // Protection for backward compatibility
1921   if(idx<(GetNumberOfTRDslices()*kTRDnPlanes)) return -1.;
1922
1923   if(sp) (*sp) = fTRDslices[idx+kTRDnPlanes];
1924   return fTRDslices[idx];
1925 }
1926
1927 //____________________________________________________
1928 Double_t  AliESDtrack::GetTRDslice(Int_t plane, Int_t slice) const {
1929   //Gets the charge from the slice of the plane
1930
1931   if(!fTRDslices) {
1932     //AliError("No TRD slices allocated for this track !");
1933     return -1.;
1934   }
1935   if ((plane<0) || (plane>=kTRDnPlanes)) {
1936     AliError("Info for TRD plane not available !");
1937     return -1.;
1938   }
1939   Int_t ns=GetNumberOfTRDslices();
1940   if ((slice<-1) || (slice>=ns)) {
1941     //AliError("Wrong TRD slice !");
1942     return -1.;
1943   }
1944
1945   if(slice>=0) return fTRDslices[plane*ns + slice];
1946
1947   // return average of the dEdx measurements
1948   Double_t q=0.; Double32_t *s = &fTRDslices[plane*ns];
1949   for (Int_t i=0; i<ns; i++, s++) if((*s)>0.) q+=(*s);
1950   return q/ns;
1951 }
1952
1953 //____________________________________________________
1954 void  AliESDtrack::SetNumberOfTRDslices(Int_t n) {
1955   //Sets the number of slices used for PID
1956   if (fTRDnSlices) return;
1957
1958   fTRDnSlices=n;
1959   fTRDslices=new Double32_t[fTRDnSlices];
1960
1961   // set-up correctly the allocated memory
1962   memset(fTRDslices, 0, n*sizeof(Double32_t));
1963   for (Int_t i=GetNumberOfTRDslices(); i--;) fTRDslices[i]=-1.;
1964 }
1965
1966 //____________________________________________________
1967 void  AliESDtrack::SetTRDslice(Double_t q, Int_t plane, Int_t slice) {
1968   //Sets the charge q in the slice of the plane
1969   if(!fTRDslices) {
1970     AliError("No TRD slices allocated for this track !");
1971     return;
1972   }
1973   if ((plane<0) || (plane>=kTRDnPlanes)) {
1974     AliError("Info for TRD plane not allocated !");
1975     return;
1976   }
1977   Int_t ns=GetNumberOfTRDslices();
1978   if ((slice<0) || (slice>=ns)) {
1979     AliError("Wrong TRD slice !");
1980     return;
1981   }
1982   Int_t n=plane*ns + slice;
1983   fTRDslices[n]=q;
1984 }
1985
1986
1987 //____________________________________________________
1988 void AliESDtrack::SetTRDmomentum(Double_t p, Int_t plane, Double_t *sp)
1989 {
1990   if(!fTRDslices) {
1991     AliError("No TRD slices allocated for this track !");
1992     return;
1993   }
1994   if ((plane<0) || (plane>=kTRDnPlanes)) {
1995     AliError("Info for TRD plane not allocated !");
1996     return;
1997   }
1998
1999   Int_t idx = fTRDnSlices-(kTRDnPlanes<<1)+plane;
2000   // Protection for backward compatibility
2001   if(idx<GetNumberOfTRDslices()*kTRDnPlanes) return;
2002
2003   if(sp) fTRDslices[idx+kTRDnPlanes] = (*sp);
2004   fTRDslices[idx] = p;
2005 }
2006
2007
2008 //_______________________________________________________________________
2009 void AliESDtrack::SetTOFpid(const Double_t *p) {
2010   // Sets the probability of each particle type (in TOF)
2011   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
2012   SetStatus(AliESDtrack::kTOFpid);
2013 }
2014
2015 //_______________________________________________________________________
2016 void AliESDtrack::SetTOFLabel(const Int_t *p) {
2017   // Sets  (in TOF)
2018   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
2019 }
2020
2021 //_______________________________________________________________________
2022 void AliESDtrack::GetTOFpid(Double_t *p) const {
2023   // Gets probabilities of each particle type (in TOF)
2024   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
2025 }
2026
2027 //_______________________________________________________________________
2028 void AliESDtrack::GetTOFLabel(Int_t *p) const {
2029   // Gets (in TOF)
2030   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
2031 }
2032
2033 //_______________________________________________________________________
2034 void AliESDtrack::GetTOFInfo(Float_t *info) const {
2035   // Gets (in TOF)
2036   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
2037 }
2038
2039 //_______________________________________________________________________
2040 void AliESDtrack::SetTOFInfo(Float_t*info) {
2041   // Gets (in TOF)
2042   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
2043 }
2044
2045
2046
2047 //_______________________________________________________________________
2048 void AliESDtrack::SetHMPIDpid(const Double_t *p) {
2049   // Sets the probability of each particle type (in HMPID)
2050   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
2051   SetStatus(AliESDtrack::kHMPIDpid);
2052 }
2053
2054 //_______________________________________________________________________
2055 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
2056   // Gets probabilities of each particle type (in HMPID)
2057   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
2058 }
2059
2060
2061
2062 //_______________________________________________________________________
2063 void AliESDtrack::SetESDpid(const Double_t *p) {
2064   // Sets the probability of each particle type for the ESD track
2065   SetPIDValues(fR,p,AliPID::kSPECIES);
2066   SetStatus(AliESDtrack::kESDpid);
2067 }
2068
2069 //_______________________________________________________________________
2070 void AliESDtrack::GetESDpid(Double_t *p) const {
2071   // Gets probability of each particle type for the ESD track
2072   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
2073 }
2074
2075 //_______________________________________________________________________
2076 Bool_t AliESDtrack::RelateToVertexTPC(const AliESDVertex *vtx,
2077 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
2078   //
2079   // Try to relate the TPC-only track parameters to the vertex "vtx",
2080   // if the (rough) transverse impact parameter is not bigger then "maxd".
2081   //            Magnetic field is "b" (kG).
2082   //
2083   // a) The TPC-only paramters are extapolated to the DCA to the vertex.
2084   // b) The impact parameters and their covariance matrix are calculated.
2085   // c) An attempt to constrain the TPC-only params to the vertex is done.
2086   //    The constrained params are returned via "cParam".
2087   //
2088   // In the case of success, the returned value is kTRUE
2089   // otherwise, it's kFALSE)
2090   //
2091
2092   if (!fTPCInner) return kFALSE;
2093   if (!vtx) return kFALSE;
2094
2095   Double_t dz[2],cov[3];
2096   if (!fTPCInner->PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
2097
2098   fdTPC = dz[0];
2099   fzTPC = dz[1];
2100   fCddTPC = cov[0];
2101   fCdzTPC = cov[1];
2102   fCzzTPC = cov[2];
2103
2104   Double_t covar[6]; vtx->GetCovMatrix(covar);
2105   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2106   Double_t c[3]={covar[2],0.,covar[5]};
2107
2108   Double_t chi2=GetPredictedChi2(p,c);
2109   if (chi2>kVeryBig) return kFALSE;
2110
2111   fCchi2TPC=chi2;
2112
2113   if (!cParam) return kTRUE;
2114
2115   *cParam = *fTPCInner;
2116   if (!cParam->Update(p,c)) return kFALSE;
2117
2118   return kTRUE;
2119 }
2120
2121 //_______________________________________________________________________
2122 Bool_t AliESDtrack::RelateToVertexTPCBxByBz(const AliESDVertex *vtx,
2123 Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam) {
2124   //
2125   // Try to relate the TPC-only track parameters to the vertex "vtx",
2126   // if the (rough) transverse impact parameter is not bigger then "maxd".
2127   //
2128   // All three components of the magnetic field ,"b[3]" (kG),
2129   // are taken into account.
2130   //
2131   // a) The TPC-only paramters are extapolated to the DCA to the vertex.
2132   // b) The impact parameters and their covariance matrix are calculated.
2133   // c) An attempt to constrain the TPC-only params to the vertex is done.
2134   //    The constrained params are returned via "cParam".
2135   //
2136   // In the case of success, the returned value is kTRUE
2137   // otherwise, it's kFALSE)
2138   //
2139
2140   if (!fTPCInner) return kFALSE;
2141   if (!vtx) return kFALSE;
2142
2143   Double_t dz[2],cov[3];
2144   if (!fTPCInner->PropagateToDCABxByBz(vtx, b, maxd, dz, cov)) return kFALSE;
2145
2146   fdTPC = dz[0];
2147   fzTPC = dz[1];
2148   fCddTPC = cov[0];
2149   fCdzTPC = cov[1];
2150   fCzzTPC = cov[2];
2151
2152   Double_t covar[6]; vtx->GetCovMatrix(covar);
2153   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2154   Double_t c[3]={covar[2],0.,covar[5]};
2155
2156   Double_t chi2=GetPredictedChi2(p,c);
2157   if (chi2>kVeryBig) return kFALSE;
2158
2159   fCchi2TPC=chi2;
2160
2161   if (!cParam) return kTRUE;
2162
2163   *cParam = *fTPCInner;
2164   if (!cParam->Update(p,c)) return kFALSE;
2165
2166   return kTRUE;
2167 }
2168
2169 //_______________________________________________________________________
2170 Bool_t AliESDtrack::RelateToVertex(const AliESDVertex *vtx,
2171 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
2172   //
2173   // Try to relate this track to the vertex "vtx",
2174   // if the (rough) transverse impact parameter is not bigger then "maxd".
2175   //            Magnetic field is "b" (kG).
2176   //
2177   // a) The track gets extapolated to the DCA to the vertex.
2178   // b) The impact parameters and their covariance matrix are calculated.
2179   // c) An attempt to constrain this track to the vertex is done.
2180   //    The constrained params are returned via "cParam".
2181   //
2182   // In the case of success, the returned value is kTRUE
2183   // (otherwise, it's kFALSE)
2184   //
2185
2186   if (!vtx) return kFALSE;
2187
2188   Double_t dz[2],cov[3];
2189   if (!PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
2190
2191   fD = dz[0];
2192   fZ = dz[1];
2193   fCdd = cov[0];
2194   fCdz = cov[1];
2195   fCzz = cov[2];
2196
2197   Double_t covar[6]; vtx->GetCovMatrix(covar);
2198   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2199   Double_t c[3]={covar[2],0.,covar[5]};
2200
2201   Double_t chi2=GetPredictedChi2(p,c);
2202   if (chi2>kVeryBig) return kFALSE;
2203
2204   fCchi2=chi2;
2205
2206
2207   //--- Could now these lines be removed ? ---
2208   delete fCp;
2209   fCp=new AliExternalTrackParam(*this);
2210
2211   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
2212   //----------------------------------------
2213
2214   fVertexID = vtx->GetID();
2215
2216   if (!cParam) return kTRUE;
2217
2218   *cParam = *this;
2219   if (!cParam->Update(p,c)) return kFALSE;
2220
2221   return kTRUE;
2222 }
2223
2224 //_______________________________________________________________________
2225 Bool_t AliESDtrack::RelateToVertexBxByBz(const AliESDVertex *vtx,
2226 Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam) {
2227   //
2228   // Try to relate this track to the vertex "vtx",
2229   // if the (rough) transverse impact parameter is not bigger then "maxd".
2230   //            Magnetic field is "b" (kG).
2231   //
2232   // a) The track gets extapolated to the DCA to the vertex.
2233   // b) The impact parameters and their covariance matrix are calculated.
2234   // c) An attempt to constrain this track to the vertex is done.
2235   //    The constrained params are returned via "cParam".
2236   //
2237   // In the case of success, the returned value is kTRUE
2238   // (otherwise, it's kFALSE)
2239   //
2240
2241   if (!vtx) return kFALSE;
2242
2243   Double_t dz[2],cov[3];
2244   if (!PropagateToDCABxByBz(vtx, b, maxd, dz, cov)) return kFALSE;
2245
2246   fD = dz[0];
2247   fZ = dz[1];
2248   fCdd = cov[0];
2249   fCdz = cov[1];
2250   fCzz = cov[2];
2251
2252   Double_t covar[6]; vtx->GetCovMatrix(covar);
2253   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2254   Double_t c[3]={covar[2],0.,covar[5]};
2255
2256   Double_t chi2=GetPredictedChi2(p,c);
2257   if (chi2>kVeryBig) return kFALSE;
2258
2259   fCchi2=chi2;
2260
2261
2262   //--- Could now these lines be removed ? ---
2263   delete fCp;
2264   fCp=new AliExternalTrackParam(*this);
2265
2266   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
2267   //----------------------------------------
2268
2269   fVertexID = vtx->GetID();
2270
2271   if (!cParam) return kTRUE;
2272
2273   *cParam = *this;
2274   if (!cParam->Update(p,c)) return kFALSE;
2275
2276   return kTRUE;
2277 }
2278
2279 //_______________________________________________________________________
2280 void AliESDtrack::Print(Option_t *) const {
2281   // Prints info on the track
2282   AliExternalTrackParam::Print();
2283   printf("ESD track info\n") ;
2284   Double_t p[AliPID::kSPECIESN] ;
2285   Int_t index = 0 ;
2286   if( IsOn(kITSpid) ){
2287     printf("From ITS: ") ;
2288     GetITSpid(p) ;
2289     for(index = 0 ; index < AliPID::kSPECIES; index++)
2290       printf("%f, ", p[index]) ;
2291     printf("\n           signal = %f\n", GetITSsignal()) ;
2292   }
2293   if( IsOn(kTPCpid) ){
2294     printf("From TPC: ") ;
2295     GetTPCpid(p) ;
2296     for(index = 0 ; index < AliPID::kSPECIES; index++)
2297       printf("%f, ", p[index]) ;
2298     printf("\n           signal = %f\n", GetTPCsignal()) ;
2299   }
2300   if( IsOn(kTRDpid) ){
2301     printf("From TRD: ") ;
2302     GetTRDpid(p) ;
2303     for(index = 0 ; index < AliPID::kSPECIES; index++)
2304       printf("%f, ", p[index]) ;
2305       printf("\n           signal = %f\n", GetTRDsignal()) ;
2306   }
2307   if( IsOn(kTOFpid) ){
2308     printf("From TOF: ") ;
2309     GetTOFpid(p) ;
2310     for(index = 0 ; index < AliPID::kSPECIES; index++)
2311       printf("%f, ", p[index]) ;
2312     printf("\n           signal = %f\n", GetTOFsignal()) ;
2313   }
2314   if( IsOn(kHMPIDpid) ){
2315     printf("From HMPID: ") ;
2316     GetHMPIDpid(p) ;
2317     for(index = 0 ; index < AliPID::kSPECIES; index++)
2318       printf("%f, ", p[index]) ;
2319     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
2320   }
2321 }
2322
2323
2324 //
2325 // Draw functionality
2326 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
2327 //
2328 void AliESDtrack::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
2329   //
2330   // Fill points in the polymarker
2331   //
2332   TObjArray arrayRef;
2333   arrayRef.AddLast(new AliExternalTrackParam(*this));
2334   if (fIp) arrayRef.AddLast(new AliExternalTrackParam(*fIp));
2335   if (fOp) arrayRef.AddLast(new AliExternalTrackParam(*fOp));
2336   if (fHMPIDp) arrayRef.AddLast(new AliExternalTrackParam(*fHMPIDp));
2337   //
2338   Double_t mpos[3]={0,0,0};
2339   Int_t entries=arrayRef.GetEntries();
2340   for (Int_t i=0;i<entries;i++){
2341     Double_t pos[3];
2342     ((AliExternalTrackParam*)arrayRef.At(i))->GetXYZ(pos);
2343     mpos[0]+=pos[0]/entries;
2344     mpos[1]+=pos[1]/entries;
2345     mpos[2]+=pos[2]/entries;
2346   }
2347   // Rotate to the mean position
2348   //
2349   Float_t fi= TMath::ATan2(mpos[1],mpos[0]);
2350   for (Int_t i=0;i<entries;i++){
2351     Bool_t res = ((AliExternalTrackParam*)arrayRef.At(i))->Rotate(fi);
2352     if (!res) delete arrayRef.RemoveAt(i);
2353   }
2354   Int_t counter=0;
2355   for (Double_t r=minR; r<maxR; r+=stepR){
2356     Double_t sweight=0;
2357     Double_t mlpos[3]={0,0,0};
2358     for (Int_t i=0;i<entries;i++){
2359       Double_t point[3]={0,0,0};
2360       AliExternalTrackParam *param = ((AliExternalTrackParam*)arrayRef.At(i));
2361       if (!param) continue;
2362       if (param->GetXYZAt(r,magF,point)){
2363         Double_t weight = 1./(10.+(r-param->GetX())*(r-param->GetX()));
2364         sweight+=weight;
2365         mlpos[0]+=point[0]*weight;
2366         mlpos[1]+=point[1]*weight;
2367         mlpos[2]+=point[2]*weight;
2368       }
2369     }
2370     if (sweight>0){
2371       mlpos[0]/=sweight;
2372       mlpos[1]/=sweight;
2373       mlpos[2]/=sweight;
2374       pol->SetPoint(counter,mlpos[0],mlpos[1], mlpos[2]);
2375       printf("xyz\t%f\t%f\t%f\n",mlpos[0], mlpos[1],mlpos[2]);
2376       counter++;
2377     }
2378   }
2379 }
2380
2381 //_______________________________________________________________________
2382 void AliESDtrack::SetITSdEdxSamples(const Double_t s[4]) {
2383   //
2384   // Store the dE/dx samples measured by the two SSD and two SDD layers.
2385   // These samples are corrected for the track segment length.
2386   //
2387   for (Int_t i=0; i<4; i++) fITSdEdxSamples[i]=s[i];
2388 }
2389
2390 //_______________________________________________________________________
2391 void AliESDtrack::GetITSdEdxSamples(Double_t *s) const {
2392   //
2393   // Get the dE/dx samples measured by the two SSD and two SDD layers.
2394   // These samples are corrected for the track segment length.
2395   //
2396   for (Int_t i=0; i<4; i++) s[i]=fITSdEdxSamples[i];
2397 }
2398
2399
2400 UShort_t   AliESDtrack::GetTPCnclsS(Int_t i0,Int_t i1) const{
2401   //
2402   // get number of shared TPC clusters
2403   //
2404   return  fTPCSharedMap.CountBits(i0)-fTPCSharedMap.CountBits(i1);
2405 }
2406
2407 UShort_t   AliESDtrack::GetTPCncls(Int_t i0,Int_t i1) const{
2408   //
2409   // get number of TPC clusters
2410   //
2411   return  fTPCClusterMap.CountBits(i0)-fTPCClusterMap.CountBits(i1);
2412 }