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