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 TPCinner parameters
 974   if(track.fTPCInner) *track.fTPCInner = *fTPCInner;
 975   else track.fTPCInner = new AliExternalTrackParam(*fTPCInner);
 976   track.fdTPC   = fdTPC;
 977   track.fzTPC   = fzTPC;
 978   track.fCddTPC = fCddTPC;
 979   track.fCdzTPC = fCdzTPC;
 980   track.fCzzTPC = fCzzTPC;
 981   track.fCchi2TPC = fCchi2TPC;
 982
 983
 984   // copy all other TPC specific parameters
 985
 986   // replace label by TPC label
 987   track.fLabel    = fTPCLabel;
 988   track.fTPCLabel = fTPCLabel;
 989
 990   track.fTPCchi2 = fTPCchi2;
 991   track.fTPCchi2Iter1 = fTPCchi2Iter1;
 992   track.fTPCsignal = fTPCsignal;
 993   track.fTPCsignalS = fTPCsignalS;
 994   for(int i = 0;i<4;++i)track.fTPCPoints[i] = fTPCPoints[i];
 995
 996   track.fTPCncls    = fTPCncls;
 997   track.fTPCnclsF   = fTPCnclsF;
 998   track.fTPCsignalN =  fTPCsignalN;
 999   track.fTPCnclsIter1    = fTPCnclsIter1;
1000   track.fTPCnclsFIter1   = fTPCnclsFIter1;
1001
1002   // PID
1003   for(int i=0;i<AliPID::kSPECIES;++i){
1004     track.fTPCr[i] = fTPCr[i];
1005     // combined PID is TPC only!
1006     track.fR[i] = fTPCr[i];
1007   }
1008   track.fTPCClusterMap = fTPCClusterMap;
1009   track.fTPCSharedMap = fTPCSharedMap;
1010
1011
1012   // reset the flags
1013   track.fFlags = kTPCin;
1014   track.fID    = fID;
1015
1016   track.fFlags |= fFlags & kTPCpid; //copy the TPCpid status flag
1017
1018   for (Int_t i=0;i<3;i++) track.fKinkIndexes[i] = fKinkIndexes[i];
1019
1020   return kTRUE;
1021
1022 }
1023
1024 //_______________________________________________________________________
1025 void AliESDtrack::MakeMiniESDtrack(){
1026   // Resets everything except
1027   // fFlags: Reconstruction status flags
1028   // fLabel: Track label
1029   // fID:  Unique ID of the track
1030   // Impact parameter information
1031   // fR[AliPID::kSPECIES]: combined "detector response probability"
1032   // Running track parameters in the base class (AliExternalTrackParam)
1033
1034   fTrackLength = 0;
1035
1036   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i] = 0;
1037
1038   // Reset track parameters constrained to the primary vertex
1039   delete fCp;fCp = 0;
1040
1041   // Reset track parameters at the inner wall of TPC
1042   delete fIp;fIp = 0;
1043   delete fTPCInner;fTPCInner=0;
1044   // Reset track parameters at the inner wall of the TRD
1045   delete fOp;fOp = 0;
1046   // Reset track parameters at the HMPID
1047   delete fHMPIDp;fHMPIDp = 0;
1048
1049
1050   // Reset ITS track related information
1051   fITSchi2 = 0;
1052   fITSncls = 0;
1053   fITSClusterMap=0;
1054   fITSSharedMap=0;
1055   fITSsignal = 0;
1056   for (Int_t i=0;i<4;i++) fITSdEdxSamples[i] = 0.;
1057   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=0;
1058   fITSLabel = 0;
1059
1060   // Reset TPC related track information
1061   fTPCchi2 = 0;
1062   fTPCchi2Iter1 = 0;
1063   fTPCncls = 0;
1064   fTPCnclsF = 0;
1065   fTPCnclsIter1 = 0;
1066   fTPCnclsFIter1 = 0;
1067   fTPCClusterMap = 0;
1068   fTPCSharedMap = 0;
1069   fTPCsignal= 0;
1070   fTPCsignalS= 0;
1071   fTPCsignalN= 0;
1072   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0;
1073   fTPCLabel=0;
1074   for (Int_t i=0;i<4;i++) fTPCPoints[i] = 0;
1075   for (Int_t i=0; i<3;i++)   fKinkIndexes[i] = 0;
1076   for (Int_t i=0; i<3;i++)   fV0Indexes[i] = 0;
1077
1078   // Reset TRD related track information
1079   fTRDchi2 = 0;
1080   fTRDncls = 0;
1081   fTRDncls0 = 0;
1082   fTRDsignal = 0;
1083   for (Int_t i=0;i<kTRDnPlanes;i++) {
1084     fTRDTimBin[i]  = 0;
1085   }
1086   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0;
1087   fTRDLabel = 0;
1088   fTRDQuality  = 0;
1089   fTRDntracklets = 0;
1090   if(fTRDnSlices)
1091     delete[] fTRDslices;
1092   fTRDslices=0x0;
1093   fTRDnSlices=0;
1094   fTRDBudget  = 0;
1095
1096   // Reset TOF related track information
1097   fTOFchi2 = 0;
1098   fTOFindex = -1;
1099   fTOFsignal = 99999;
1100   fTOFCalChannel = -1;
1101   fTOFsignalToT = 99999;
1102   fTOFsignalRaw = 99999;
1103   fTOFsignalDz = 999;
1104   fTOFsignalDx = 999;
1105   fTOFdeltaBC = 999;
1106   fTOFl0l1 = 999;
1107   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i] = 0;
1108   for (Int_t i=0;i<3;i++) fTOFLabel[i] = -1;
1109   for (Int_t i=0;i<10;i++) fTOFInfo[i] = 0;
1110
1111   // Reset HMPID related track information
1112   fHMPIDchi2 = 0;
1113   fHMPIDqn = 0;
1114   fHMPIDcluIdx = -1;
1115   fHMPIDsignal = 0;
1116   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i] = 0;
1117   fHMPIDtrkTheta = 0;
1118   fHMPIDtrkPhi = 0;
1119   fHMPIDtrkX = 0;
1120   fHMPIDtrkY = 0;
1121   fHMPIDmipX = 0;
1122   fHMPIDmipY = 0;
1123   fCaloIndex = kEMCALNoMatch;
1124
1125   // reset global track chi2
1126   fGlobalChi2 = 0;
1127
1128   fVertexID = -2; // an orphan track
1129
1130   delete fFriendTrack; fFriendTrack = 0;
1131 }
1132 //_______________________________________________________________________
1133 Double_t AliESDtrack::GetMass() const {
1134   // Returns the mass of the most probable particle type
1135
1136   Int_t i;
1137   for (i=0; i<AliPID::kSPECIES-1; i++) {
1138       if (fR[i] != fR[i+1]) break;
1139   }
1140   // If all the probabilities are equal, return the pion mass
1141   if (i == AliPID::kSPECIES-1) return AliPID::ParticleMass(AliPID::kPion);
1142
1143   Float_t max=0.;
1144   Int_t k=-1;
1145   for (i=0; i<AliPID::kSPECIES; i++) {
1146     if (fR[i]>max) {k=i; max=fR[i];}
1147   }
1148   if (k==0) { // dE/dx "crossing points" in the TPC
1149      Double_t p=GetP();
1150      if ((p>0.38)&&(p<0.48))
1151         if (fR[0]<fR[3]*10.) return AliPID::ParticleMass(AliPID::kKaon);
1152      if ((p>0.75)&&(p<0.85))
1153         if (fR[0]<fR[4]*10.) return AliPID::ParticleMass(AliPID::kProton);
1154      return 0.00051;
1155   }
1156   if (k==1) return AliPID::ParticleMass(AliPID::kMuon);
1157   if (k==2||k==-1) return AliPID::ParticleMass(AliPID::kPion);
1158   if (k==3) return AliPID::ParticleMass(AliPID::kKaon);
1159   if (k==4) return AliPID::ParticleMass(AliPID::kProton);
1160   AliWarning("Undefined mass !");
1161   return AliPID::ParticleMass(AliPID::kPion);
1162 }
1163
1164 //______________________________________________________________________________
1165 Double_t AliESDtrack::M() const
1166 {
1167   // Returns the assumed mass
1168   // (the pion mass, if the particle can't be identified properly).
1169
1170   AliWarning("This is the ESD mass. Use it with care !");
1171   return GetMass();
1172 }
1173
1174 //______________________________________________________________________________
1175 Double_t AliESDtrack::E() const
1176 {
1177   // Returns the energy of the particle given its assumed mass.
1178   // Assumes the pion mass if the particle can't be identified properly.
1179
1180   Double_t m = M();
1181   Double_t p = P();
1182   return TMath::Sqrt(p*p + m*m);
1183 }
1184
1185 //______________________________________________________________________________
1186 Double_t AliESDtrack::Y() const
1187 {
1188   // Returns the rapidity of a particle given its assumed mass.
1189   // Assumes the pion mass if the particle can't be identified properly.
1190
1191   Double_t e = E();
1192   Double_t pz = Pz();
1193   if (e != TMath::Abs(pz)) { // energy was not equal to pz
1194     return 0.5*TMath::Log((e+pz)/(e-pz));
1195   } else { // energy was equal to pz
1196     return -999.;
1197   }
1198 }
1199
1200 //_______________________________________________________________________
1201 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
1202   //
1203   // This function updates track's running parameters
1204   //
1205   Bool_t rc=kTRUE;
1206
1207   SetStatus(flags);
1208   fLabel=t->GetLabel();
1209
1210   if (t->IsStartedTimeIntegral()) {
1211     SetStatus(kTIME);
1212     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
1213     SetIntegratedLength(t->GetIntegratedLength());
1214   }
1215
1216   Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1217   if (flags==kITSout) fFriendTrack->SetITSOut(*t);
1218   if (flags==kTPCout) fFriendTrack->SetTPCOut(*t);
1219   if (flags==kTRDrefit) fFriendTrack->SetTRDIn(*t);
1220
1221   switch (flags) {
1222
1223   case kITSin: case kITSout: case kITSrefit:
1224     {
1225     fITSClusterMap=0;
1226     fITSncls=t->GetNumberOfClusters();
1227     Int_t* indexITS = new Int_t[AliESDfriendTrack::kMaxITScluster];
1228     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
1229         indexITS[i]=t->GetClusterIndex(i);
1230
1231         if (i<fITSncls) {
1232           Int_t l=(indexITS[i] & 0xf0000000) >> 28;
1233            SETBIT(fITSClusterMap,l);
1234         }
1235     }
1236     fFriendTrack->SetITSIndices(indexITS,AliESDfriendTrack::kMaxITScluster);
1237     delete [] indexITS;
1238
1239     fITSchi2=t->GetChi2();
1240     fITSsignal=t->GetPIDsignal();
1241     fITSLabel = t->GetLabel();
1242     // keep in fOp the parameters outside ITS for ITS stand-alone tracks
1243     if (flags==kITSout) {
1244       if (!fOp) fOp=new AliExternalTrackParam(*t);
1245       else
1246         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1247     }
1248     }
1249     break;
1250
1251   case kTPCin: case kTPCrefit:
1252     {
1253     fTPCLabel = t->GetLabel();
1254     if (flags==kTPCin)  {
1255       fTPCInner=new AliExternalTrackParam(*t);
1256       fTPCnclsIter1=t->GetNumberOfClusters();
1257       fTPCchi2Iter1=t->GetChi2();
1258     }
1259     if (!fIp) fIp=new AliExternalTrackParam(*t);
1260     else
1261       fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1262     }
1263   case kTPCout:
1264     {
1265     Int_t* indexTPC = new Int_t[AliESDfriendTrack::kMaxTPCcluster];
1266     if (flags & kTPCout){
1267       if (!fOp) fOp=new AliExternalTrackParam(*t);
1268       else
1269         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1270     }
1271     fTPCncls=t->GetNumberOfClusters();
1272     fTPCchi2=t->GetChi2();
1273
1274      {//prevrow must be declared in separate namespace, otherwise compiler cries:
1275       //"jump to case label crosses initialization of `Int_t prevrow'"
1276        Int_t prevrow = -1;
1277        //       for (Int_t i=0;i<fTPCncls;i++)
1278        for (Int_t i=0;i<AliESDfriendTrack::kMaxTPCcluster;i++)
1279         {
1280           indexTPC[i]=t->GetClusterIndex(i);
1281           Int_t idx = indexTPC[i];
1282
1283           if (idx<0) continue;
1284
1285           // Piotr's Cluster Map for HBT
1286           // ### please change accordingly if cluster array is changing
1287           // to "New TPC Tracking" style (with gaps in array)
1288           Int_t sect = (idx&0xff000000)>>24;
1289           Int_t row = (idx&0x00ff0000)>>16;
1290           if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
1291
1292           fTPCClusterMap.SetBitNumber(row,kTRUE);
1293
1294           //Fill the gap between previous row and this row with 0 bits
1295           //In case  ###  pleas change it as well - just set bit 0 in case there
1296           //is no associated clusters for current "i"
1297           if (prevrow < 0)
1298            {
1299              prevrow = row;//if previous bit was not assigned yet == this is the first one
1300            }
1301           else
1302            { //we don't know the order (inner to outer or reverse)
1303              //just to be save in case it is going to change
1304              Int_t n = 0, m = 0;
1305              if (prevrow < row)
1306               {
1307                 n = prevrow;
1308                 m = row;
1309               }
1310              else
1311               {
1312                 n = row;
1313                 m = prevrow;
1314               }
1315
1316              for (Int_t j = n+1; j < m; j++)
1317               {
1318                 fTPCClusterMap.SetBitNumber(j,kFALSE);
1319               }
1320              prevrow = row;
1321            }
1322           // End Of Piotr's Cluster Map for HBT
1323         }
1324         fFriendTrack->SetTPCIndices(indexTPC,AliESDfriendTrack::kMaxTPCcluster);
1325         delete [] indexTPC;
1326
1327      }
1328     fTPCsignal=t->GetPIDsignal();
1329     }
1330     break;
1331
1332   case kTRDin: case kTRDrefit:
1333     break;
1334   case kTRDout:
1335     {
1336     fTRDLabel = t->GetLabel();
1337     fTRDchi2  = t->GetChi2();
1338     fTRDncls  = t->GetNumberOfClusters();
1339       Int_t* indexTRD = new Int_t[AliESDfriendTrack::kMaxTRDcluster];
1340       for (Int_t i=0;i<AliESDfriendTrack::kMaxTRDcluster;i++) indexTRD[i]=-2;
1341       for (Int_t i=0;i<6;i++) indexTRD[i]=t->GetTrackletIndex(i);
1342       fFriendTrack->SetTRDIndices(indexTRD,AliESDfriendTrack::kMaxTRDcluster);
1343       delete [] indexTRD;
1344
1345
1346     fTRDsignal=t->GetPIDsignal();
1347     }
1348     break;
1349   case kTRDbackup:
1350     if (!fOp) fOp=new AliExternalTrackParam(*t);
1351     else
1352       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1353     fTRDncls0 = t->GetNumberOfClusters();
1354     break;
1355   case kTOFin:
1356     break;
1357   case kTOFout:
1358     break;
1359   case kTRDStop:
1360     break;
1361   case kHMPIDout:
1362   if (!fHMPIDp) fHMPIDp=new AliExternalTrackParam(*t);
1363     else
1364       fHMPIDp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1365     break;
1366   default:
1367     AliError("Wrong flag !");
1368     return kFALSE;
1369   }
1370
1371   return rc;
1372 }
1373
1374 //_______________________________________________________________________
1375 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
1376   //---------------------------------------------------------------------
1377   // This function returns external representation of the track parameters
1378   //---------------------------------------------------------------------
1379   x=GetX();
1380   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
1381 }
1382
1383 //_______________________________________________________________________
1384 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
1385   //---------------------------------------------------------------------
1386   // This function returns external representation of the cov. matrix
1387   //---------------------------------------------------------------------
1388   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
1389 }
1390
1391 //_______________________________________________________________________
1392 Bool_t AliESDtrack::GetConstrainedExternalParameters
1393                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1394   //---------------------------------------------------------------------
1395   // This function returns the constrained external track parameters
1396   //---------------------------------------------------------------------
1397   if (!fCp) return kFALSE;
1398   alpha=fCp->GetAlpha();
1399   x=fCp->GetX();
1400   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
1401   return kTRUE;
1402 }
1403
1404 //_______________________________________________________________________
1405 Bool_t
1406 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
1407   //---------------------------------------------------------------------
1408   // This function returns the constrained external cov. matrix
1409   //---------------------------------------------------------------------
1410   if (!fCp) return kFALSE;
1411   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
1412   return kTRUE;
1413 }
1414
1415 Bool_t
1416 AliESDtrack::GetInnerExternalParameters
1417                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1418   //---------------------------------------------------------------------
1419   // This function returns external representation of the track parameters
1420   // at the inner layer of TPC
1421   //---------------------------------------------------------------------
1422   if (!fIp) return kFALSE;
1423   alpha=fIp->GetAlpha();
1424   x=fIp->GetX();
1425   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
1426   return kTRUE;
1427 }
1428
1429 Bool_t
1430 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
1431  //---------------------------------------------------------------------
1432  // This function returns external representation of the cov. matrix
1433  // at the inner layer of TPC
1434  //---------------------------------------------------------------------
1435   if (!fIp) return kFALSE;
1436   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
1437   return kTRUE;
1438 }
1439
1440 void
1441 AliESDtrack::SetOuterParam(const AliExternalTrackParam *p, ULong_t flags) {
1442   //
1443   // This is a direct setter for the outer track parameters
1444   //
1445   SetStatus(flags);
1446   if (fOp) delete fOp;
1447   fOp=new AliExternalTrackParam(*p);
1448 }
1449
1450 void
1451 AliESDtrack::SetOuterHmpParam(const AliExternalTrackParam *p, ULong_t flags) {
1452   //
1453   // This is a direct setter for the outer track parameters
1454   //
1455   SetStatus(flags);
1456   if (fHMPIDp) delete fHMPIDp;
1457   fHMPIDp=new AliExternalTrackParam(*p);
1458 }
1459
1460 Bool_t
1461 AliESDtrack::GetOuterExternalParameters
1462                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1463   //---------------------------------------------------------------------
1464   // This function returns external representation of the track parameters
1465   // at the inner layer of TRD
1466   //---------------------------------------------------------------------
1467   if (!fOp) return kFALSE;
1468   alpha=fOp->GetAlpha();
1469   x=fOp->GetX();
1470   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
1471   return kTRUE;
1472 }
1473
1474 Bool_t
1475 AliESDtrack::GetOuterHmpExternalParameters
1476                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1477   //---------------------------------------------------------------------
1478   // This function returns external representation of the track parameters
1479   // at the inner layer of TRD
1480   //---------------------------------------------------------------------
1481   if (!fHMPIDp) return kFALSE;
1482   alpha=fHMPIDp->GetAlpha();
1483   x=fHMPIDp->GetX();
1484   for (Int_t i=0; i<5; i++) p[i]=fHMPIDp->GetParameter()[i];
1485   return kTRUE;
1486 }
1487
1488 Bool_t
1489 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
1490  //---------------------------------------------------------------------
1491  // This function returns external representation of the cov. matrix
1492  // at the inner layer of TRD
1493  //---------------------------------------------------------------------
1494   if (!fOp) return kFALSE;
1495   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
1496   return kTRUE;
1497 }
1498
1499 Bool_t
1500 AliESDtrack::GetOuterHmpExternalCovariance(Double_t cov[15]) const {
1501  //---------------------------------------------------------------------
1502  // This function returns external representation of the cov. matrix
1503  // at the inner layer of TRD
1504  //---------------------------------------------------------------------
1505   if (!fHMPIDp) return kFALSE;
1506   for (Int_t i=0; i<15; i++) cov[i]=fHMPIDp->GetCovariance()[i];
1507   return kTRUE;
1508 }
1509
1510 Int_t AliESDtrack::GetNcls(Int_t idet) const
1511 {
1512   // Get number of clusters by subdetector index
1513   //
1514   Int_t ncls = 0;
1515   switch(idet){
1516   case 0:
1517     ncls = fITSncls;
1518     break;
1519   case 1:
1520     ncls = fTPCncls;
1521     break;
1522   case 2:
1523     ncls = fTRDncls;
1524     break;
1525   case 3:
1526     if (fTOFindex != -1)
1527       ncls = 1;
1528     break;
1529   case 4: //PHOS
1530     break;
1531   case 5: //HMPID
1532     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1533       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1534         ncls = 1;
1535       }
1536     }
1537     break;
1538   default:
1539     break;
1540   }
1541   return ncls;
1542 }
1543
1544 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
1545 {
1546   // Get cluster index array by subdetector index
1547   //
1548   Int_t ncls = 0;
1549   switch(idet){
1550   case 0:
1551     ncls = GetITSclusters(idx);
1552     break;
1553   case 1:
1554     ncls = GetTPCclusters(idx);
1555     break;
1556   case 2:
1557     ncls = GetTRDclusters(idx);
1558     break;
1559   case 3:
1560     if (fTOFindex != -1) {
1561       idx[0] = fTOFindex;
1562       ncls = 1;
1563     }
1564     break;
1565   case 4: //PHOS
1566     break;
1567   case 5:
1568     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1569       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1570         idx[0] = GetHMPIDcluIdx();
1571         ncls = 1;
1572       }
1573     }
1574     break;
1575   case 6: //EMCAL
1576     break;
1577   default:
1578     break;
1579   }
1580   return ncls;
1581 }
1582
1583 //_______________________________________________________________________
1584 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
1585   // Returns the array with integrated times for each particle hypothesis
1586   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
1587 }
1588
1589 //_______________________________________________________________________
1590 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
1591   // Sets the array with integrated times for each particle hypotesis
1592   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
1593 }
1594
1595 //_______________________________________________________________________
1596 void AliESDtrack::SetITSpid(const Double_t *p) {
1597   // Sets values for the probability of each particle type (in ITS)
1598   SetPIDValues(fITSr,p,AliPID::kSPECIES);
1599   SetStatus(AliESDtrack::kITSpid);
1600 }
1601
1602 //_______________________________________________________________________
1603 void AliESDtrack::GetITSpid(Double_t *p) const {
1604   // Gets the probability of each particle type (in ITS)
1605   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
1606 }
1607
1608 //_______________________________________________________________________
1609 Char_t AliESDtrack::GetITSclusters(Int_t *idx) const {
1610   //---------------------------------------------------------------------
1611   // This function returns indices of the assgined ITS clusters
1612   //---------------------------------------------------------------------
1613   if (idx) {
1614     Int_t *index=fFriendTrack->GetITSindices();
1615     for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) {
1616       if ( (i>=fITSncls) && (i<6) ) idx[i]=-1;
1617       else {
1618         if (index) {
1619           idx[i]=index[i];
1620         }
1621         else idx[i]= -2;
1622       }
1623     }
1624   }
1625   return fITSncls;
1626 }
1627
1628 //_______________________________________________________________________
1629 Bool_t AliESDtrack::GetITSModuleIndexInfo(Int_t ilayer,Int_t &idet,Int_t &status,
1630                                          Float_t &xloc,Float_t &zloc) const {
1631   //----------------------------------------------------------------------
1632   // This function encodes in the module number also the status of cluster association
1633   // "status" can have the following values:
1634   // 1 "found" (cluster is associated),
1635   // 2 "dead" (module is dead from OCDB),
1636   // 3 "skipped" (module or layer forced to be skipped),
1637   // 4 "outinz" (track out of z acceptance),
1638   // 5 "nocls" (no clusters in the road),
1639   // 6 "norefit" (cluster rejected during refit),
1640   // 7 "deadzspd" (holes in z in SPD)
1641   // Also given are the coordinates of the crossing point of track and module
1642   // (in the local module ref. system)
1643   // WARNING: THIS METHOD HAS TO BE SYNCHRONIZED WITH AliITStrackV2::GetModuleIndexInfo()!
1644   //----------------------------------------------------------------------
1645
1646   if(fITSModule[ilayer]==-1) {
1647     idet = -1;
1648     status=0;
1649     xloc=-99.; zloc=-99.;
1650     return kFALSE;
1651   }
1652
1653   Int_t module = fITSModule[ilayer];
1654
1655   idet = Int_t(module/1000000);
1656
1657   module -= idet*1000000;
1658
1659   status = Int_t(module/100000);
1660
1661   module -= status*100000;
1662
1663   Int_t signs = Int_t(module/10000);
1664
1665   module-=signs*10000;
1666
1667   Int_t xInt = Int_t(module/100);
1668   module -= xInt*100;
1669
1670   Int_t zInt = module;
1671
1672   if(signs==1) { xInt*=1; zInt*=1; }
1673   if(signs==2) { xInt*=1; zInt*=-1; }
1674   if(signs==3) { xInt*=-1; zInt*=1; }
1675   if(signs==4) { xInt*=-1; zInt*=-1; }
1676
1677   xloc = 0.1*(Float_t)xInt;
1678   zloc = 0.1*(Float_t)zInt;
1679
1680   if(status==4) idet = -1;
1681
1682   return kTRUE;
1683 }
1684
1685 //_______________________________________________________________________
1686 UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
1687   //---------------------------------------------------------------------
1688   // This function returns indices of the assgined ITS clusters
1689   //---------------------------------------------------------------------
1690   if (idx) {
1691     Int_t *index=fFriendTrack->GetTPCindices();
1692
1693     if (index){
1694       for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
1695     }
1696     else {
1697       for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=-2;
1698     }
1699   }
1700   return fTPCncls;
1701 }
1702
1703 //_______________________________________________________________________
1704 Float_t AliESDtrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
1705 {
1706   //
1707   // TPC cluster information
1708   // type 0: get fraction of found/findable clusters with neighbourhood definition
1709   //      1: findable clusters with neighbourhood definition
1710   //      2: found clusters
1711   //
1712   // definition of findable clusters:
1713   //            a cluster is defined as findable if there is another cluster
1714   //           within +- nNeighbours pad rows. The idea is to overcome threshold
1715   //           effects with a very simple algorithm.
1716   //
1717
1718   if (type==2) return fTPCClusterMap.CountBits();
1719
1720   Int_t found=0;
1721   Int_t findable=0;
1722   Int_t last=-nNeighbours;
1723
1724   for (Int_t i=row0; i<row1; ++i){
1725     //look to current row
1726     if (fTPCClusterMap[i]) {
1727       last=i;
1728       ++found;
1729       ++findable;
1730       continue;
1731     }
1732     //look to nNeighbours before
1733     if ((i-last)<=nNeighbours) {
1734       ++findable;
1735       continue;
1736     }
1737     //look to nNeighbours after
1738     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
1739       if (fTPCClusterMap[j]){
1740         ++findable;
1741         break;
1742       }
1743     }
1744   }
1745   if (type==1) return findable;
1746
1747   if (type==0){
1748     Float_t fraction=0;
1749     if (findable>0)
1750       fraction=(Float_t)found/(Float_t)findable;
1751     else
1752       fraction=0;
1753     return fraction;
1754   }
1755   return 0;  // undefined type - default value
1756 }
1757
1758 //_______________________________________________________________________
1759 Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
1760   //
1761   // GetDensity of the clusters on given region between row0 and row1
1762   // Dead zone effect takin into acoount
1763   //
1764   Int_t good  = 0;
1765   Int_t found = 0;
1766   //
1767   Int_t *index=fFriendTrack->GetTPCindices();
1768   for (Int_t i=row0;i<=row1;i++){
1769     Int_t idx = index[i];
1770     if (idx!=-1)  good++;             // track outside of dead zone
1771     if (idx>0)    found++;
1772   }
1773   Float_t density=0.5;
1774   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
1775   return density;
1776 }
1777
1778 //_______________________________________________________________________
1779 void AliESDtrack::SetTPCpid(const Double_t *p) {
1780   // Sets values for the probability of each particle type (in TPC)
1781   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
1782   SetStatus(AliESDtrack::kTPCpid);
1783 }
1784
1785 //_______________________________________________________________________
1786 void AliESDtrack::GetTPCpid(Double_t *p) const {
1787   // Gets the probability of each particle type (in TPC)
1788   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
1789 }
1790
1791 //_______________________________________________________________________
1792 UChar_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
1793   //---------------------------------------------------------------------
1794   // This function returns indices of the assgined TRD clusters
1795   //---------------------------------------------------------------------
1796   if (idx) {
1797     Int_t *index=fFriendTrack->GetTRDindices();
1798
1799     if (index) {
1800       for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
1801     }
1802     else {
1803       for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=-2;
1804     }
1805   }
1806   return fTRDncls;
1807 }
1808
1809 //_______________________________________________________________________
1810 UChar_t AliESDtrack::GetTRDtracklets(Int_t *idx) const {
1811 //
1812 // This function returns the number of TRD tracklets used in tracking
1813 // and it fills the indices of these tracklets in the array "idx" as they
1814 // are registered in the TRD track list.
1815 //
1816 // Caution :
1817 //   1. The idx array has to be allocated with a size >= AliESDtrack::kTRDnPlanes
1818 //   2. The idx array store not only the index but also the layer of the tracklet.
1819 //      Therefore tracks with TRD gaps contain default values for indices [-1]
1820
1821   if (!idx) return GetTRDntracklets();
1822   Int_t *index=fFriendTrack->GetTRDindices();
1823   Int_t n = 0;
1824   for (Int_t i=0; i<kTRDnPlanes; i++){
1825     if (index){
1826       if(index[i]>=0) n++;
1827       idx[i]=index[i];
1828     }
1829     else idx[i] = -2;
1830   }
1831   return n;
1832 }
1833
1834 //_______________________________________________________________________
1835 void AliESDtrack::SetTRDpid(const Double_t *p) {
1836   // Sets values for the probability of each particle type (in TRD)
1837   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
1838   SetStatus(AliESDtrack::kTRDpid);
1839 }
1840
1841 //_______________________________________________________________________
1842 void AliESDtrack::GetTRDpid(Double_t *p) const {
1843   // Gets the probability of each particle type (in TRD)
1844   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
1845 }
1846
1847 //_______________________________________________________________________
1848 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
1849 {
1850   // Sets the probability of particle type iSpecies to p (in TRD)
1851   fTRDr[iSpecies] = p;
1852 }
1853
1854 Double_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
1855 {
1856   // Returns the probability of particle type iSpecies (in TRD)
1857   return fTRDr[iSpecies];
1858 }
1859
1860 //____________________________________________________
1861 Int_t AliESDtrack::GetNumberOfTRDslices() const
1862 {
1863   // built in backward compatibility
1864   Int_t idx = fTRDnSlices - (kTRDnPlanes<<1);
1865   return idx<18 ? fTRDnSlices/kTRDnPlanes : idx/kTRDnPlanes;
1866 }
1867
1868 //____________________________________________________
1869 Double_t AliESDtrack::GetTRDmomentum(Int_t plane, Double_t *sp) const
1870 {
1871 //Returns momentum estimation and optional its error (sp)
1872 // in TRD layer "plane".
1873
1874   if (!fTRDnSlices) {
1875     AliDebug(2, "No TRD info allocated for this track.");
1876     return -1.;
1877   }
1878   if ((plane<0) || (plane>=kTRDnPlanes)) {
1879     AliWarning(Form("Request for TRD plane[%d] outside range.", plane));
1880     return -1.;
1881   }
1882
1883   Int_t idx = fTRDnSlices-(kTRDnPlanes<<1)+plane;
1884   // Protection for backward compatibility
1885   if(idx<(GetNumberOfTRDslices()*kTRDnPlanes)) return -1.;
1886
1887   if(sp) (*sp) = fTRDslices[idx+kTRDnPlanes];
1888   return fTRDslices[idx];
1889 }
1890
1891 //____________________________________________________
1892 Double_t  AliESDtrack::GetTRDslice(Int_t plane, Int_t slice) const {
1893   //Gets the charge from the slice of the plane
1894
1895   if(!fTRDslices) {
1896     //AliError("No TRD slices allocated for this track !");
1897     return -1.;
1898   }
1899   if ((plane<0) || (plane>=kTRDnPlanes)) {
1900     AliError("Info for TRD plane not available !");
1901     return -1.;
1902   }
1903   Int_t ns=GetNumberOfTRDslices();
1904   if ((slice<-1) || (slice>=ns)) {
1905     //AliError("Wrong TRD slice !");
1906     return -1.;
1907   }
1908
1909   if(slice>=0) return fTRDslices[plane*ns + slice];
1910
1911   // return average of the dEdx measurements
1912   Double_t q=0.; Double32_t *s = &fTRDslices[plane*ns];
1913   for (Int_t i=0; i<ns; i++, s++) if((*s)>0.) q+=(*s);
1914   return q/ns;
1915 }
1916
1917 //____________________________________________________
1918 void  AliESDtrack::SetNumberOfTRDslices(Int_t n) {
1919   //Sets the number of slices used for PID
1920   if (fTRDnSlices) return;
1921
1922   fTRDnSlices=n;
1923   fTRDslices=new Double32_t[fTRDnSlices];
1924
1925   // set-up correctly the allocated memory
1926   memset(fTRDslices, 0, n*sizeof(Double32_t));
1927   for (Int_t i=GetNumberOfTRDslices(); i--;) fTRDslices[i]=-1.;
1928 }
1929
1930 //____________________________________________________
1931 void  AliESDtrack::SetTRDslice(Double_t q, Int_t plane, Int_t slice) {
1932   //Sets the charge q in the slice of the plane
1933   if(!fTRDslices) {
1934     AliError("No TRD slices allocated for this track !");
1935     return;
1936   }
1937   if ((plane<0) || (plane>=kTRDnPlanes)) {
1938     AliError("Info for TRD plane not allocated !");
1939     return;
1940   }
1941   Int_t ns=GetNumberOfTRDslices();
1942   if ((slice<0) || (slice>=ns)) {
1943     AliError("Wrong TRD slice !");
1944     return;
1945   }
1946   Int_t n=plane*ns + slice;
1947   fTRDslices[n]=q;
1948 }
1949
1950
1951 //____________________________________________________
1952 void AliESDtrack::SetTRDmomentum(Double_t p, Int_t plane, Double_t *sp)
1953 {
1954   if(!fTRDslices) {
1955     AliError("No TRD slices allocated for this track !");
1956     return;
1957   }
1958   if ((plane<0) || (plane>=kTRDnPlanes)) {
1959     AliError("Info for TRD plane not allocated !");
1960     return;
1961   }
1962
1963   Int_t idx = fTRDnSlices-(kTRDnPlanes<<1)+plane;
1964   // Protection for backward compatibility
1965   if(idx<GetNumberOfTRDslices()*kTRDnPlanes) return;
1966
1967   if(sp) fTRDslices[idx+kTRDnPlanes] = (*sp);
1968   fTRDslices[idx] = p;
1969 }
1970
1971
1972 //_______________________________________________________________________
1973 void AliESDtrack::SetTOFpid(const Double_t *p) {
1974   // Sets the probability of each particle type (in TOF)
1975   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
1976   SetStatus(AliESDtrack::kTOFpid);
1977 }
1978
1979 //_______________________________________________________________________
1980 void AliESDtrack::SetTOFLabel(const Int_t *p) {
1981   // Sets  (in TOF)
1982   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
1983 }
1984
1985 //_______________________________________________________________________
1986 void AliESDtrack::GetTOFpid(Double_t *p) const {
1987   // Gets probabilities of each particle type (in TOF)
1988   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
1989 }
1990
1991 //_______________________________________________________________________
1992 void AliESDtrack::GetTOFLabel(Int_t *p) const {
1993   // Gets (in TOF)
1994   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
1995 }
1996
1997 //_______________________________________________________________________
1998 void AliESDtrack::GetTOFInfo(Float_t *info) const {
1999   // Gets (in TOF)
2000   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
2001 }
2002
2003 //_______________________________________________________________________
2004 void AliESDtrack::SetTOFInfo(Float_t*info) {
2005   // Gets (in TOF)
2006   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
2007 }
2008
2009
2010
2011 //_______________________________________________________________________
2012 void AliESDtrack::SetHMPIDpid(const Double_t *p) {
2013   // Sets the probability of each particle type (in HMPID)
2014   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
2015   SetStatus(AliESDtrack::kHMPIDpid);
2016 }
2017
2018 //_______________________________________________________________________
2019 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
2020   // Gets probabilities of each particle type (in HMPID)
2021   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
2022 }
2023
2024
2025
2026 //_______________________________________________________________________
2027 void AliESDtrack::SetESDpid(const Double_t *p) {
2028   // Sets the probability of each particle type for the ESD track
2029   SetPIDValues(fR,p,AliPID::kSPECIES);
2030   SetStatus(AliESDtrack::kESDpid);
2031 }
2032
2033 //_______________________________________________________________________
2034 void AliESDtrack::GetESDpid(Double_t *p) const {
2035   // Gets probability of each particle type for the ESD track
2036   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
2037 }
2038
2039 //_______________________________________________________________________
2040 Bool_t AliESDtrack::RelateToVertexTPC(const AliESDVertex *vtx,
2041 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
2042   //
2043   // Try to relate the TPC-only track parameters to the vertex "vtx",
2044   // if the (rough) transverse impact parameter is not bigger then "maxd".
2045   //            Magnetic field is "b" (kG).
2046   //
2047   // a) The TPC-only paramters are extapolated to the DCA to the vertex.
2048   // b) The impact parameters and their covariance matrix are calculated.
2049   // c) An attempt to constrain the TPC-only params to the vertex is done.
2050   //    The constrained params are returned via "cParam".
2051   //
2052   // In the case of success, the returned value is kTRUE
2053   // otherwise, it's kFALSE)
2054   //
2055
2056   if (!fTPCInner) return kFALSE;
2057   if (!vtx) return kFALSE;
2058
2059   Double_t dz[2],cov[3];
2060   if (!fTPCInner->PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
2061
2062   fdTPC = dz[0];
2063   fzTPC = dz[1];
2064   fCddTPC = cov[0];
2065   fCdzTPC = cov[1];
2066   fCzzTPC = cov[2];
2067
2068   Double_t covar[6]; vtx->GetCovMatrix(covar);
2069   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2070   Double_t c[3]={covar[2],0.,covar[5]};
2071
2072   Double_t chi2=GetPredictedChi2(p,c);
2073   if (chi2>kVeryBig) return kFALSE;
2074
2075   fCchi2TPC=chi2;
2076
2077   if (!cParam) return kTRUE;
2078
2079   *cParam = *fTPCInner;
2080   if (!cParam->Update(p,c)) return kFALSE;
2081
2082   return kTRUE;
2083 }
2084
2085 //_______________________________________________________________________
2086 Bool_t AliESDtrack::RelateToVertexTPCBxByBz(const AliESDVertex *vtx,
2087 Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam) {
2088   //
2089   // Try to relate the TPC-only track parameters to the vertex "vtx",
2090   // if the (rough) transverse impact parameter is not bigger then "maxd".
2091   //
2092   // All three components of the magnetic field ,"b[3]" (kG),
2093   // are taken into account.
2094   //
2095   // a) The TPC-only paramters are extapolated to the DCA to the vertex.
2096   // b) The impact parameters and their covariance matrix are calculated.
2097   // c) An attempt to constrain the TPC-only params to the vertex is done.
2098   //    The constrained params are returned via "cParam".
2099   //
2100   // In the case of success, the returned value is kTRUE
2101   // otherwise, it's kFALSE)
2102   //
2103
2104   if (!fTPCInner) return kFALSE;
2105   if (!vtx) return kFALSE;
2106
2107   Double_t dz[2],cov[3];
2108   if (!fTPCInner->PropagateToDCABxByBz(vtx, b, maxd, dz, cov)) return kFALSE;
2109
2110   fdTPC = dz[0];
2111   fzTPC = dz[1];
2112   fCddTPC = cov[0];
2113   fCdzTPC = cov[1];
2114   fCzzTPC = cov[2];
2115
2116   Double_t covar[6]; vtx->GetCovMatrix(covar);
2117   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2118   Double_t c[3]={covar[2],0.,covar[5]};
2119
2120   Double_t chi2=GetPredictedChi2(p,c);
2121   if (chi2>kVeryBig) return kFALSE;
2122
2123   fCchi2TPC=chi2;
2124
2125   if (!cParam) return kTRUE;
2126
2127   *cParam = *fTPCInner;
2128   if (!cParam->Update(p,c)) return kFALSE;
2129
2130   return kTRUE;
2131 }
2132
2133 //_______________________________________________________________________
2134 Bool_t AliESDtrack::RelateToVertex(const AliESDVertex *vtx,
2135 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
2136   //
2137   // Try to relate this track to the vertex "vtx",
2138   // if the (rough) transverse impact parameter is not bigger then "maxd".
2139   //            Magnetic field is "b" (kG).
2140   //
2141   // a) The track gets extapolated to the DCA to the vertex.
2142   // b) The impact parameters and their covariance matrix are calculated.
2143   // c) An attempt to constrain this track to the vertex is done.
2144   //    The constrained params are returned via "cParam".
2145   //
2146   // In the case of success, the returned value is kTRUE
2147   // (otherwise, it's kFALSE)
2148   //
2149
2150   if (!vtx) return kFALSE;
2151
2152   Double_t dz[2],cov[3];
2153   if (!PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
2154
2155   fD = dz[0];
2156   fZ = dz[1];
2157   fCdd = cov[0];
2158   fCdz = cov[1];
2159   fCzz = cov[2];
2160
2161   Double_t covar[6]; vtx->GetCovMatrix(covar);
2162   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2163   Double_t c[3]={covar[2],0.,covar[5]};
2164
2165   Double_t chi2=GetPredictedChi2(p,c);
2166   if (chi2>kVeryBig) return kFALSE;
2167
2168   fCchi2=chi2;
2169
2170
2171   //--- Could now these lines be removed ? ---
2172   delete fCp;
2173   fCp=new AliExternalTrackParam(*this);
2174
2175   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
2176   //----------------------------------------
2177
2178   fVertexID = vtx->GetID();
2179
2180   if (!cParam) return kTRUE;
2181
2182   *cParam = *this;
2183   if (!cParam->Update(p,c)) return kFALSE;
2184
2185   return kTRUE;
2186 }
2187
2188 //_______________________________________________________________________
2189 Bool_t AliESDtrack::RelateToVertexBxByBz(const AliESDVertex *vtx,
2190 Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam) {
2191   //
2192   // Try to relate this track to the vertex "vtx",
2193   // if the (rough) transverse impact parameter is not bigger then "maxd".
2194   //            Magnetic field is "b" (kG).
2195   //
2196   // a) The track gets extapolated to the DCA to the vertex.
2197   // b) The impact parameters and their covariance matrix are calculated.
2198   // c) An attempt to constrain this track to the vertex is done.
2199   //    The constrained params are returned via "cParam".
2200   //
2201   // In the case of success, the returned value is kTRUE
2202   // (otherwise, it's kFALSE)
2203   //
2204
2205   if (!vtx) return kFALSE;
2206
2207   Double_t dz[2],cov[3];
2208   if (!PropagateToDCABxByBz(vtx, b, maxd, dz, cov)) return kFALSE;
2209
2210   fD = dz[0];
2211   fZ = dz[1];
2212   fCdd = cov[0];
2213   fCdz = cov[1];
2214   fCzz = cov[2];
2215
2216   Double_t covar[6]; vtx->GetCovMatrix(covar);
2217   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2218   Double_t c[3]={covar[2],0.,covar[5]};
2219
2220   Double_t chi2=GetPredictedChi2(p,c);
2221   if (chi2>kVeryBig) return kFALSE;
2222
2223   fCchi2=chi2;
2224
2225
2226   //--- Could now these lines be removed ? ---
2227   delete fCp;
2228   fCp=new AliExternalTrackParam(*this);
2229
2230   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
2231   //----------------------------------------
2232
2233   fVertexID = vtx->GetID();
2234
2235   if (!cParam) return kTRUE;
2236
2237   *cParam = *this;
2238   if (!cParam->Update(p,c)) return kFALSE;
2239
2240   return kTRUE;
2241 }
2242
2243 //_______________________________________________________________________
2244 void AliESDtrack::Print(Option_t *) const {
2245   // Prints info on the track
2246   AliExternalTrackParam::Print();
2247   printf("ESD track info\n") ;
2248   Double_t p[AliPID::kSPECIESN] ;
2249   Int_t index = 0 ;
2250   if( IsOn(kITSpid) ){
2251     printf("From ITS: ") ;
2252     GetITSpid(p) ;
2253     for(index = 0 ; index < AliPID::kSPECIES; index++)
2254       printf("%f, ", p[index]) ;
2255     printf("\n           signal = %f\n", GetITSsignal()) ;
2256   }
2257   if( IsOn(kTPCpid) ){
2258     printf("From TPC: ") ;
2259     GetTPCpid(p) ;
2260     for(index = 0 ; index < AliPID::kSPECIES; index++)
2261       printf("%f, ", p[index]) ;
2262     printf("\n           signal = %f\n", GetTPCsignal()) ;
2263   }
2264   if( IsOn(kTRDpid) ){
2265     printf("From TRD: ") ;
2266     GetTRDpid(p) ;
2267     for(index = 0 ; index < AliPID::kSPECIES; index++)
2268       printf("%f, ", p[index]) ;
2269       printf("\n           signal = %f\n", GetTRDsignal()) ;
2270   }
2271   if( IsOn(kTOFpid) ){
2272     printf("From TOF: ") ;
2273     GetTOFpid(p) ;
2274     for(index = 0 ; index < AliPID::kSPECIES; index++)
2275       printf("%f, ", p[index]) ;
2276     printf("\n           signal = %f\n", GetTOFsignal()) ;
2277   }
2278   if( IsOn(kHMPIDpid) ){
2279     printf("From HMPID: ") ;
2280     GetHMPIDpid(p) ;
2281     for(index = 0 ; index < AliPID::kSPECIES; index++)
2282       printf("%f, ", p[index]) ;
2283     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
2284   }
2285 }
2286
2287
2288 //
2289 // Draw functionality
2290 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
2291 //
2292 void AliESDtrack::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
2293   //
2294   // Fill points in the polymarker
2295   //
2296   TObjArray arrayRef;
2297   arrayRef.AddLast(new AliExternalTrackParam(*this));
2298   if (fIp) arrayRef.AddLast(new AliExternalTrackParam(*fIp));
2299   if (fOp) arrayRef.AddLast(new AliExternalTrackParam(*fOp));
2300   if (fHMPIDp) arrayRef.AddLast(new AliExternalTrackParam(*fHMPIDp));
2301   //
2302   Double_t mpos[3]={0,0,0};
2303   Int_t entries=arrayRef.GetEntries();
2304   for (Int_t i=0;i<entries;i++){
2305     Double_t pos[3];
2306     ((AliExternalTrackParam*)arrayRef.At(i))->GetXYZ(pos);
2307     mpos[0]+=pos[0]/entries;
2308     mpos[1]+=pos[1]/entries;
2309     mpos[2]+=pos[2]/entries;
2310   }
2311   // Rotate to the mean position
2312   //
2313   Float_t fi= TMath::ATan2(mpos[1],mpos[0]);
2314   for (Int_t i=0;i<entries;i++){
2315     Bool_t res = ((AliExternalTrackParam*)arrayRef.At(i))->Rotate(fi);
2316     if (!res) delete arrayRef.RemoveAt(i);
2317   }
2318   Int_t counter=0;
2319   for (Double_t r=minR; r<maxR; r+=stepR){
2320     Double_t sweight=0;
2321     Double_t mlpos[3]={0,0,0};
2322     for (Int_t i=0;i<entries;i++){
2323       Double_t point[3]={0,0,0};
2324       AliExternalTrackParam *param = ((AliExternalTrackParam*)arrayRef.At(i));
2325       if (!param) continue;
2326       if (param->GetXYZAt(r,magF,point)){
2327         Double_t weight = 1./(10.+(r-param->GetX())*(r-param->GetX()));
2328         sweight+=weight;
2329         mlpos[0]+=point[0]*weight;
2330         mlpos[1]+=point[1]*weight;
2331         mlpos[2]+=point[2]*weight;
2332       }
2333     }
2334     if (sweight>0){
2335       mlpos[0]/=sweight;
2336       mlpos[1]/=sweight;
2337       mlpos[2]/=sweight;
2338       pol->SetPoint(counter,mlpos[0],mlpos[1], mlpos[2]);
2339       printf("xyz\t%f\t%f\t%f\n",mlpos[0], mlpos[1],mlpos[2]);
2340       counter++;
2341     }
2342   }
2343 }
2344
2345 //_______________________________________________________________________
2346 void AliESDtrack::SetITSdEdxSamples(const Double_t s[4]) {
2347   //
2348   // Store the dE/dx samples measured by the two SSD and two SDD layers.
2349   // These samples are corrected for the track segment length.
2350   //
2351   for (Int_t i=0; i<4; i++) fITSdEdxSamples[i]=s[i];
2352 }
2353
2354 //_______________________________________________________________________
2355 void AliESDtrack::GetITSdEdxSamples(Double_t *s) const {
2356   //
2357   // Get the dE/dx samples measured by the two SSD and two SDD layers.
2358   // These samples are corrected for the track segment length.
2359   //
2360   for (Int_t i=0; i<4; i++) s[i]=fITSdEdxSamples[i];
2361 }
2362
2363
2364 UShort_t   AliESDtrack::GetTPCnclsS(Int_t i0,Int_t i1) const{
2365   //
2366   // get number of shared TPC clusters
2367   //
2368   return  fTPCSharedMap.CountBits(i0)-fTPCSharedMap.CountBits(i1);
2369 }
2370
2371 UShort_t   AliESDtrack::GetTPCncls(Int_t i0,Int_t i1) const{
2372   //
2373   // get number of TPC clusters
2374   //
2375   return  fTPCClusterMap.CountBits(i0)-fTPCClusterMap.CountBits(i1);
2376 }