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