STEER/ESD/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 //     Content:
  37 //        a.) Track parameters and covariance - AliExternalTrackParam  - snapshots along trajectory at differnt tracking steps
  38 //            current, fIp, fTPCinner, fCp, fOp, fHMPIDp, + friendTrack (fITSout, fTPCout, fTRDin)
  39 //        b.) Flags - per detector status bits
  40 //        c.) Track fit quality -chi2, number of clusters
  41 //        d.) Detector PID information
  42 //        d.) Different detector specific information (e.g number of tracklets in TRD, TOF cluster descriptor ...)
  43 //
  44 //
  45 //     The track information is propagated from one tracking detector to
  46 //     other using the functionality of AliESDtrack (AliExternalTrackParam - current parameters)
  47 //
  48 //     Barrel tracking uses Kalman filtering technique (no global fit model is used). Kalman provides optimal local
  49 //     track parameters estimate at given position under certian assumptions.
  50 //     Following approximations were used:
  51 //           a.) gaussian Multiple scattering
  52 //           b.) gaussian energy loss
  53 //           c.) gaussian error of the space point residuals
  54 //
  55 //     Kalman filter take into account following effects which are
  56 //     difficult to handle using global fit:
  57 //        a.) Multiple scattering
  58 //        b.) Energy loss
  59 //        c.) Non homogenous magnetic field
  60 //
  61 //     In general case, following barrel detectors are contributing to
  62 //     the Kalman track information:
  63 //         a. TPC
  64 //         b. ITS
  65 //         c. TRD
  66 //
  67 //      Track findind/fitting procedure is done in 3 steps:
  68 //         1. Cluster2Track(in)   - inward sequence TPC->ITS
  69 //         2. PropagateBack(out)   - outward sequence ITS->TPC->TRD -> Outer PID detectors
  70 //         3. RefitInward(refit)     - inward sequence TRD->TPC->ITS
  71 //      After each recosntruction step detector status is updated in the data member fFlags
  72 //      fFlags|=k<DetectorName><step> where step={1:in,2:out,3:refit,}
  73 //      For some of detectors a special flags were implemented. Declaration of list of all flags can be find in $ALICE_ROOT/STEER/STEERBase/AliVTrack.h
  74 //
  75 //
  76 //      The current track parameter is updated after each detector (see bellow).
  77 //      In specical cases a track  snapshots (AliExternalTrackParam) are stored
  78 //
  79 //
  80 //
  81 //      For some type of analysis (+visualization) track local parameters at
  82 //      different position are neccesary. A snapshots during the track
  83 //      propagation are created and stored either in track itself (for analysis purposes) or assiciated friend track (for calibration and debugging purposes)
  84 //      (See AliExternalTrackParam class for desctiption of variables and functionality)
  85 //      Snapshots:
  86 //      a.)  Current parameters (AliESDtrack itself)
  87 //         Description:  as obtained in the last succesfull tracking step
  88 //         Contributors: best case TRD->TPC->ITS after RefitInward
  89 //         Recomended way to get track parameters. It includes all of the information.
  90 //         NOTICE - By default the track parameters are stored at the DCA point to the primary vertex.
  91 //                  Optimal for primary tracks, far from optimal for deeply secondary tracks.
  92 //                  To get optimal track parameters at the secondary vertex, OnTheFly V0s with associated track parameters should be used
  93 //
  94 //      b.) Constrained parameters (fCp)
  95 //         Description:
  96 //            Kalman track  updated with the Primary vertex information with corresponding error (soft  constraint - see http://en.wikipedia.org/wiki/Constraint_(mathematics)#Hard_and_soft_constraints)
  97 //         Function:
  98 //            const AliExternalTrackParam *GetConstrainedParam() const {return fCp;}
  99 //         Contributors: best case TRD->TPC->ITS after RefitInward
 100 //         Recommended usage: Use only for tracks selected as primary (check GetConstrainedChi2())
 101 //         NOTICE - not real constraint only soft constraint
 102 //
 103 //      c.) Inner parameters (fIp) -
 104 //         Description: Track parameters at inner wall of the TPC
 105 //         Function:
 106 //           const AliExternalTrackParam *GetInnerParam() const { return fIp;}
 107 //         Contributors: general case TRD->TPC (during RefitInward)
 108 //         Recomended usage: To provide momenta for the TPC PID and to estimate quality of the track determination for further
 109 //
 110 //      d.)  TPCinnerParam (fTPCinner):
 111 //         Description: TPC only parameters at DCA to the primary vertex (corrected for the material between TPC and vertex)
 112 //         Function:
 113 //                const AliExternalTrackParam *GetTPCInnerParam() const {return fTPCInner;}
 114 //         Contributors:  TPC only from in step  1 (Cluster2Tracks)
 115 //         Recomended usage: Requested for HBT study
 116 //                           (smaller correlations as using also ITS information)
 117 //         NOTICE: Optimal for primary, far from optimal for secondary tracks (similar to track parameters a.)
 118 //                 ! We should always use the c.) fIp in case of the TPC PID analysis, or undo material budget correction!
 119 //
 120 //      e.) Outer parameters - (fOp)
 121 //         Description: track parameters during PropagateBack in the last sucessfull propagation
 122 //                  Reason to generate backup  OuterParameters
 123 //                        a.) Local inclination angle bigger than threshold -
 124 //                            Low momenta tracks
 125 //                        b.) Catastrofic (large relative>~20 %)energy loss in material outside of the TPC
 126 //                        c.) No additional  space points contributing to track
 127 //         Function:
 128 //            const AliExternalTrackParam *GetOuterParam() const { return fOp;}
 129 //         Contributors:  general case - ITS-> TPC -> TRD ( during PropagateBack )
 130 //         Recomended usage:
 131 //              a.) Tracking: Starting parameter for Refit inward
 132 //              b.) Visualization
 133 //              c.) QA
 134 //         NOTICE: Should be not used for the physic analysis
 135 //         Function:
 136 //            const AliExternalTrackParam *GetOuterParam() const { return fOp;}
 137 //
 138 //-----------------------------------------------------------------
 139
 140 #include <TMath.h>
 141 #include <TParticle.h>
 142 #include <TDatabasePDG.h>
 143 #include <TMatrixD.h>
 144
 145 #include "AliESDVertex.h"
 146 #include "AliESDtrack.h"
 147 #include "AliESDEvent.h"
 148 #include "AliKalmanTrack.h"
 149 #include "AliVTrack.h"
 150 #include "AliLog.h"
 151 #include "AliTrackPointArray.h"
 152 #include "TPolyMarker3D.h"
 153 #include "AliTrackerBase.h"
 154 #include "AliTPCdEdxInfo.h"
 155 #include "AliDetectorPID.h"
 156 #include "TTreeStream.h"
 157 #include "TObjArray.h"
 158
 159 ClassImp(AliESDtrack)
 160
 161 void SetPIDValues(Double_t * dest, const Double_t * src, Int_t n) {
 162   // This function copies "n" PID weights from "scr" to "dest"
 163   // and normalizes their sum to 1 thus producing conditional probabilities.
 164   // The negative weights are set to 0.
 165   // In case all the weights are non-positive they are replaced by
 166   // uniform probabilities
 167
 168   if (n<=0) return;
 169
 170   Float_t uniform = 1./(Float_t)n;
 171
 172   Float_t sum = 0;
 173   for (Int_t i=0; i<n; i++)
 174     if (src[i]>=0) {
 175       sum+=src[i];
 176       dest[i] = src[i];
 177     }
 178     else {
 179       dest[i] = 0;
 180     }
 181
 182   if(sum>0)
 183     for (Int_t i=0; i<n; i++) dest[i] /= sum;
 184   else
 185     for (Int_t i=0; i<n; i++) dest[i] = uniform;
 186 }
 187
 188 //_______________________________________________________________________
 189 AliESDtrack::AliESDtrack() :
 190   AliExternalTrackParam(),
 191   fCp(0),
 192   fIp(0),
 193   fTPCInner(0),
 194   fOp(0),
 195   fHMPIDp(0),
 196   fFriendTrack(NULL),
 197   fTPCFitMap(159),//number of padrows
 198   fTPCClusterMap(159),//number of padrows
 199   fTPCSharedMap(159),//number of padrows
 200   fFlags(0),
 201   fID(0),
 202   fLabel(0),
 203   fITSLabel(0),
 204   fTPCLabel(0),
 205   fTRDLabel(0),
 206   fTOFLabel(NULL),
 207   fTOFCalChannel(-1),
 208   fTOFindex(-1),
 209   fHMPIDqn(0),
 210   fHMPIDcluIdx(-1),
 211   fCaloIndex(kEMCALNoMatch),
 212   fR(0),
 213   fITSr(0),
 214   fTPCr(0),
 215   fTRDr(0),
 216   fTOFr(0),
 217   fHMPIDr(0),
 218   fHMPIDtrkTheta(0),
 219   fHMPIDtrkPhi(0),
 220   fHMPIDsignal(0),
 221   fTrackTime(0),
 222   fTrackLength(0),
 223   fdTPC(0),fzTPC(0),
 224   fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
 225   fCchi2TPC(0),
 226   fD(0),fZ(0),
 227   fCdd(0),fCdz(0),fCzz(0),
 228   fCchi2(0),
 229   fITSchi2(0),
 230   fTPCchi2(0),
 231   fTPCchi2Iter1(0),
 232   fTRDchi2(0),
 233   fTOFchi2(0),
 234   fHMPIDchi2(0),
 235   fGlobalChi2(0),
 236   fITSsignal(0),
 237   fTPCsignal(0),
 238   fTPCsignalTuned(0),
 239   fTPCsignalS(0),
 240   fTPCdEdxInfo(0),
 241   fTRDsignal(0),
 242   fTRDQuality(0),
 243   fTRDBudget(0),
 244   fTOFsignal(99999),
 245   fTOFsignalTuned(99999),
 246   fTOFsignalToT(99999),
 247   fTOFsignalRaw(99999),
 248   fTOFsignalDz(999),
 249   fTOFsignalDx(999),
 250   fTOFdeltaBC(999),
 251   fTOFl0l1(999),
 252   fCaloDx(0),
 253   fCaloDz(0),
 254   fHMPIDtrkX(0),
 255   fHMPIDtrkY(0),
 256   fHMPIDmipX(0),
 257   fHMPIDmipY(0),
 258   fTPCncls(0),
 259   fTPCnclsF(0),
 260   fTPCsignalN(0),
 261   fTPCnclsIter1(0),
 262   fTPCnclsFIter1(0),
 263   fITSncls(0),
 264   fITSClusterMap(0),
 265   fITSSharedMap(0),
 266   fTRDncls(0),
 267   fTRDncls0(0),
 268   fTRDntracklets(0),
 269   fTRDNchamberdEdx(0),
 270   fTRDNclusterdEdx(0),
 271   fTRDnSlices(0),
 272   fTRDslices(0x0),
 273   fVertexID(-2),// -2 means an orphan track
 274   fPIDForTracking(AliPID::kPion),
 275   fESDEvent(0),
 276   fCacheNCrossedRows(-10),
 277   fCacheChi2TPCConstrainedVsGlobal(-10),
 278   fCacheChi2TPCConstrainedVsGlobalVertex(0),
 279   fDetectorPID(0x0),
 280   fTrackPhiOnEMCal(-999),
 281   fTrackEtaOnEMCal(-999),
 282   fTrackPtOnEMCal(-999),
 283   fNtofClusters(0),
 284   fTOFcluster(NULL)
 285 {
 286   //
 287   // The default ESD constructor
 288   //
 289   if (!OnlineMode()) fFriendTrack=new AliESDfriendTrack();
 290
 291   Int_t i;
 292   for (i=kNITSchi2Std;i--;) fITSchi2Std[i] = 0;
 293
 294   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
 295   for (i=0; i<3; i++)   { fV0Indexes[i]=0;}
 296   for (i=0;i<kTRDnPlanes;i++) {
 297     fTRDTimBin[i]=0;
 298   }
 299   for (i=0;i<4;i++) {fITSdEdxSamples[i]=0.;}
 300   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
 301   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
 302   for (i=0;i<12;i++) {fITSModule[i]=-1;}
 303 }
 304
 305 bool AliESDtrack::fgkOnlineMode=false;
 306
 307 //_______________________________________________________________________
 308 AliESDtrack::AliESDtrack(const AliESDtrack& track):
 309   AliExternalTrackParam(track),
 310   fCp(0),
 311   fIp(0),
 312   fTPCInner(0),
 313   fOp(0),
 314   fHMPIDp(0),
 315   fFriendTrack(0),
 316   fTPCFitMap(track.fTPCFitMap),
 317   fTPCClusterMap(track.fTPCClusterMap),
 318   fTPCSharedMap(track.fTPCSharedMap),
 319   fFlags(track.fFlags),
 320   fID(track.fID),
 321   fLabel(track.fLabel),
 322   fITSLabel(track.fITSLabel),
 323   fTPCLabel(track.fTPCLabel),
 324   fTRDLabel(track.fTRDLabel),
 325   fTOFLabel(NULL),
 326   fTOFCalChannel(track.fTOFCalChannel),
 327   fTOFindex(track.fTOFindex),
 328   fHMPIDqn(track.fHMPIDqn),
 329   fHMPIDcluIdx(track.fHMPIDcluIdx),
 330   fCaloIndex(track.fCaloIndex),
 331   fR(0),
 332   fITSr(0),
 333   fTPCr(0),
 334   fTRDr(0),
 335   fTOFr(0),
 336   fHMPIDr(0),
 337   fHMPIDtrkTheta(track.fHMPIDtrkTheta),
 338   fHMPIDtrkPhi(track.fHMPIDtrkPhi),
 339   fHMPIDsignal(track.fHMPIDsignal),
 340   fTrackTime(NULL),
 341   fTrackLength(track.fTrackLength),
 342   fdTPC(track.fdTPC),fzTPC(track.fzTPC),
 343   fCddTPC(track.fCddTPC),fCdzTPC(track.fCdzTPC),fCzzTPC(track.fCzzTPC),
 344   fCchi2TPC(track.fCchi2TPC),
 345   fD(track.fD),fZ(track.fZ),
 346   fCdd(track.fCdd),fCdz(track.fCdz),fCzz(track.fCzz),
 347   fCchi2(track.fCchi2),
 348   fITSchi2(track.fITSchi2),
 349   fTPCchi2(track.fTPCchi2),
 350   fTPCchi2Iter1(track.fTPCchi2Iter1),
 351   fTRDchi2(track.fTRDchi2),
 352   fTOFchi2(track.fTOFchi2),
 353   fHMPIDchi2(track.fHMPIDchi2),
 354   fGlobalChi2(track.fGlobalChi2),
 355   fITSsignal(track.fITSsignal),
 356   fTPCsignal(track.fTPCsignal),
 357   fTPCsignalTuned(track.fTPCsignalTuned),
 358   fTPCsignalS(track.fTPCsignalS),
 359   fTPCdEdxInfo(0),
 360   fTRDsignal(track.fTRDsignal),
 361   fTRDQuality(track.fTRDQuality),
 362   fTRDBudget(track.fTRDBudget),
 363   fTOFsignal(track.fTOFsignal),
 364   fTOFsignalTuned(track.fTOFsignalTuned),
 365   fTOFsignalToT(track.fTOFsignalToT),
 366   fTOFsignalRaw(track.fTOFsignalRaw),
 367   fTOFsignalDz(track.fTOFsignalDz),
 368   fTOFsignalDx(track.fTOFsignalDx),
 369   fTOFdeltaBC(track.fTOFdeltaBC),
 370   fTOFl0l1(track.fTOFl0l1),
 371   fCaloDx(track.fCaloDx),
 372   fCaloDz(track.fCaloDz),
 373   fHMPIDtrkX(track.fHMPIDtrkX),
 374   fHMPIDtrkY(track.fHMPIDtrkY),
 375   fHMPIDmipX(track.fHMPIDmipX),
 376   fHMPIDmipY(track.fHMPIDmipY),
 377   fTPCncls(track.fTPCncls),
 378   fTPCnclsF(track.fTPCnclsF),
 379   fTPCsignalN(track.fTPCsignalN),
 380   fTPCnclsIter1(track.fTPCnclsIter1),
 381   fTPCnclsFIter1(track.fTPCnclsIter1),
 382   fITSncls(track.fITSncls),
 383   fITSClusterMap(track.fITSClusterMap),
 384   fITSSharedMap(track.fITSSharedMap),
 385   fTRDncls(track.fTRDncls),
 386   fTRDncls0(track.fTRDncls0),
 387   fTRDntracklets(track.fTRDntracklets),
 388   fTRDNchamberdEdx(track.fTRDNchamberdEdx),
 389   fTRDNclusterdEdx(track.fTRDNclusterdEdx),
 390   fTRDnSlices(track.fTRDnSlices),
 391   fTRDslices(0x0),
 392   fVertexID(track.fVertexID),
 393   fPIDForTracking(AliPID::kPion),
 394   fESDEvent(track.fESDEvent),
 395   fCacheNCrossedRows(track.fCacheNCrossedRows),
 396   fCacheChi2TPCConstrainedVsGlobal(track.fCacheChi2TPCConstrainedVsGlobal),
 397   fCacheChi2TPCConstrainedVsGlobalVertex(track.fCacheChi2TPCConstrainedVsGlobalVertex),
 398   fDetectorPID(0x0),
 399   fTrackPhiOnEMCal(track.fTrackPhiOnEMCal),
 400   fTrackEtaOnEMCal(track.fTrackEtaOnEMCal),
 401   fTrackPtOnEMCal(track.fTrackPtOnEMCal),
 402   fNtofClusters(track.fNtofClusters),
 403   fTOFcluster(NULL)
 404 {
 405   //
 406   //copy constructor
 407   //
 408   for (Int_t i=kNITSchi2Std;i--;) fITSchi2Std[i] = track.fITSchi2Std[i];
 409
 410   if(track.fTrackTime){
 411     fTrackTime = new Double32_t[AliPID::kSPECIESC];
 412     for (Int_t i=0;i<AliPID::kSPECIESC;i++) fTrackTime[i]=track.fTrackTime[i];
 413   }
 414
 415   if (track.fR) {
 416     fR = new Double32_t[AliPID::kSPECIES];
 417     for (Int_t i=AliPID::kSPECIES;i--;)  fR[i]=track.fR[i];
 418   }
 419   if (track.fITSr) {
 420     fITSr = new Double32_t[AliPID::kSPECIES];
 421     for (Int_t i=AliPID::kSPECIES;i--;)  fITSr[i]=track.fITSr[i];
 422   }
 423   //
 424   if (track.fTPCr) {
 425     fTPCr = new Double32_t[AliPID::kSPECIES];
 426     for (Int_t i=AliPID::kSPECIES;i--;) fTPCr[i]=track.fTPCr[i];
 427   }
 428
 429   for (Int_t i=0;i<4;i++) {fITSdEdxSamples[i]=track.fITSdEdxSamples[i];}
 430   for (Int_t i=0;i<4;i++) {fTPCPoints[i]=track.fTPCPoints[i];}
 431   for (Int_t i=0; i<3;i++)   { fKinkIndexes[i]=track.fKinkIndexes[i];}
 432   for (Int_t i=0; i<3;i++)   { fV0Indexes[i]=track.fV0Indexes[i];}
 433   //
 434   for (Int_t i=0;i<kTRDnPlanes;i++) {
 435     fTRDTimBin[i]=track.fTRDTimBin[i];
 436   }
 437
 438   if (fTRDnSlices) {
 439     fTRDslices=new Double32_t[fTRDnSlices];
 440     for (Int_t i=0; i<fTRDnSlices; i++) fTRDslices[i]=track.fTRDslices[i];
 441   }
 442
 443   if (track.fDetectorPID) fDetectorPID = new AliDetectorPID(*track.fDetectorPID);
 444
 445   if (track.fTRDr) {
 446     fTRDr  = new Double32_t[AliPID::kSPECIES];
 447     for (Int_t i=AliPID::kSPECIES;i--;) fTRDr[i]=track.fTRDr[i];
 448   }
 449
 450   if (track.fTOFr) {
 451     fTOFr = new Double32_t[AliPID::kSPECIES];
 452     for (Int_t i=AliPID::kSPECIES;i--;) fTOFr[i]=track.fTOFr[i];
 453   }
 454
 455   if(track.fTOFLabel){
 456     if(!fTOFLabel) fTOFLabel = new Int_t[3];
 457     for (Int_t i=0;i<3;i++) fTOFLabel[i]=track.fTOFLabel[i];
 458   }
 459
 460   for (Int_t i=0;i<10;i++) fTOFInfo[i]=track.fTOFInfo[i];
 461   for (Int_t i=0;i<12;i++) fITSModule[i]=track.fITSModule[i];
 462
 463   if (track.fHMPIDr) {
 464     fHMPIDr = new Double32_t[AliPID::kSPECIES];
 465     for (Int_t i=AliPID::kSPECIES;i--;) fHMPIDr[i]=track.fHMPIDr[i];
 466   }
 467
 468   if (track.fCp) fCp=new AliExternalTrackParam(*track.fCp);
 469   if (track.fIp) fIp=new AliExternalTrackParam(*track.fIp);
 470   if (track.fTPCInner) fTPCInner=new AliExternalTrackParam(*track.fTPCInner);
 471   if (track.fOp) fOp=new AliExternalTrackParam(*track.fOp);
 472   if (track.fHMPIDp) fHMPIDp=new AliExternalTrackParam(*track.fHMPIDp);
 473   if (track.fTPCdEdxInfo) fTPCdEdxInfo = new AliTPCdEdxInfo(*track.fTPCdEdxInfo);
 474
 475
 476   if (track.fFriendTrack) fFriendTrack=new AliESDfriendTrack(*(track.fFriendTrack));
 477
 478   if(fNtofClusters > 0){
 479     fTOFcluster = new Int_t[fNtofClusters];
 480         for(Int_t i=0;i < fNtofClusters;i++) fTOFcluster[i] = track.fTOFcluster[i];
 481   }
 482 }
 483
 484 //_______________________________________________________________________
 485 AliESDtrack::AliESDtrack(const AliVTrack *track) :
 486   AliExternalTrackParam(track),
 487   fCp(0),
 488   fIp(0),
 489   fTPCInner(0),
 490   fOp(0),
 491   fHMPIDp(0),
 492   fFriendTrack(0),
 493   fTPCFitMap(159),//number of padrows
 494   fTPCClusterMap(159),//number of padrows
 495   fTPCSharedMap(159),//number of padrows
 496   fFlags(0),
 497   fID(),
 498   fLabel(0),
 499   fITSLabel(0),
 500   fTPCLabel(0),
 501   fTRDLabel(0),
 502   fTOFLabel(NULL),
 503   fTOFCalChannel(-1),
 504   fTOFindex(-1),
 505   fHMPIDqn(0),
 506   fHMPIDcluIdx(-1),
 507   fCaloIndex(kEMCALNoMatch),
 508   fR(0),
 509   fITSr(0),
 510   fTPCr(0),
 511   fTRDr(0),
 512   fTOFr(0),
 513   fHMPIDr(0),
 514   fHMPIDtrkTheta(0),
 515   fHMPIDtrkPhi(0),
 516   fHMPIDsignal(0),
 517   fTrackTime(NULL),
 518   fTrackLength(0),
 519   fdTPC(0),fzTPC(0),
 520   fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
 521   fCchi2TPC(0),
 522   fD(0),fZ(0),
 523   fCdd(0),fCdz(0),fCzz(0),
 524   fCchi2(0),
 525   fITSchi2(0),
 526   fTPCchi2(0),
 527   fTPCchi2Iter1(0),
 528   fTRDchi2(0),
 529   fTOFchi2(0),
 530   fHMPIDchi2(0),
 531   fGlobalChi2(0),
 532   fITSsignal(0),
 533   fTPCsignal(0),
 534   fTPCsignalTuned(0),
 535   fTPCsignalS(0),
 536   fTPCdEdxInfo(0),
 537   fTRDsignal(0),
 538   fTRDQuality(0),
 539   fTRDBudget(0),
 540   fTOFsignal(99999),
 541   fTOFsignalTuned(99999),
 542   fTOFsignalToT(99999),
 543   fTOFsignalRaw(99999),
 544   fTOFsignalDz(999),
 545   fTOFsignalDx(999),
 546   fTOFdeltaBC(999),
 547   fTOFl0l1(999),
 548   fCaloDx(0),
 549   fCaloDz(0),
 550   fHMPIDtrkX(0),
 551   fHMPIDtrkY(0),
 552   fHMPIDmipX(0),
 553   fHMPIDmipY(0),
 554   fTPCncls(0),
 555   fTPCnclsF(0),
 556   fTPCsignalN(0),
 557   fTPCnclsIter1(0),
 558   fTPCnclsFIter1(0),
 559   fITSncls(0),
 560   fITSClusterMap(0),
 561   fITSSharedMap(0),
 562   fTRDncls(0),
 563   fTRDncls0(0),
 564   fTRDntracklets(0),
 565   fTRDNchamberdEdx(0),
 566   fTRDNclusterdEdx(0),
 567   fTRDnSlices(0),
 568   fTRDslices(0x0),
 569   fVertexID(-2),  // -2 means an orphan track
 570   fPIDForTracking(track->GetPIDForTracking()),
 571   fESDEvent(0),
 572   fCacheNCrossedRows(-10),
 573   fCacheChi2TPCConstrainedVsGlobal(-10),
 574   fCacheChi2TPCConstrainedVsGlobalVertex(0),
 575   fDetectorPID(0x0),
 576   fTrackPhiOnEMCal(-999),
 577   fTrackEtaOnEMCal(-999),
 578   fTrackPtOnEMCal(-999),
 579   fNtofClusters(0),
 580   fTOFcluster(NULL)
 581 {
 582   //
 583   // ESD track from AliVTrack.
 584   // This is not a copy constructor !
 585   //
 586
 587   if (track->InheritsFrom("AliExternalTrackParam")) {
 588      AliError("This is not a copy constructor. Use AliESDtrack(const AliESDtrack &) !");
 589      AliWarning("Calling the default constructor...");
 590      AliESDtrack();
 591      return;
 592   }
 593
 594   // Reset all the arrays
 595   Int_t i;
 596   for (i=kNITSchi2Std;i--;) fITSchi2Std[i] = 0;
 597
 598   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
 599   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
 600   for (i=0;i<kTRDnPlanes;i++) {
 601     fTRDTimBin[i]=0;
 602   }
 603   for (i=0;i<4;i++) {fITSdEdxSamples[i]=0.;}
 604   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
 605   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
 606   for (i=0;i<12;i++) {fITSModule[i]=-1;}
 607
 608
 609   // Set ITS cluster map
 610   fITSClusterMap=track->GetITSClusterMap();
 611   fITSSharedMap=0;
 612
 613   fITSncls=0;
 614   for(i=0; i<6; i++) {
 615     if(HasPointOnITSLayer(i)) fITSncls++;
 616   }
 617
 618   // Set TPC ncls
 619   fTPCncls=track->GetTPCNcls();
 620   fTPCnclsF=track->GetTPCNclsF();
 621   // TPC cluster maps
 622   const TBits* bmap = track->GetTPCClusterMapPtr();
 623   if (bmap) SetTPCClusterMap(*bmap);
 624   bmap = GetTPCFitMapPtr();
 625   if (bmap) SetTPCFitMap(*bmap);
 626   bmap = GetTPCSharedMapPtr();
 627   if (bmap) SetTPCSharedMap(*bmap);
 628   //
 629   // Set the combined PID
 630   const Double_t *pid = track->PID();
 631   if(pid) {
 632     fR = new Double32_t[AliPID::kSPECIES];
 633     for (i=AliPID::kSPECIES; i--;) fR[i]=pid[i];
 634   }
 635   //
 636   // calo matched cluster id
 637   SetEMCALcluster(track->GetEMCALcluster());
 638   // AliESD track label
 639   //
 640   // PID info
 641   fITSsignal = track->GetITSsignal();
 642   double itsdEdx[4];
 643   track->GetITSdEdxSamples(itsdEdx);
 644   SetITSdEdxSamples(itsdEdx);
 645   //
 646   SetTPCsignal(track->GetTPCsignal(),fTPCsignalS,track->GetTPCsignalN()); // No signalS in AODPi
 647   AliTPCdEdxInfo * dEdxInfo = track->GetTPCdEdxInfo();
 648   if (dEdxInfo) SetTPCdEdxInfo(new AliTPCdEdxInfo(*dEdxInfo));
 649   //
 650   SetTRDsignal(track->GetTRDsignal());
 651   int ntrdsl = track->GetNumberOfTRDslices();
 652   if (ntrdsl>0) {
 653     SetNumberOfTRDslices((ntrdsl+2)*kTRDnPlanes);
 654     for (int ipl=kTRDnPlanes;ipl--;){
 655       for (int isl=ntrdsl;isl--;) SetTRDslice(track->GetTRDslice(ipl,isl),ipl,isl);
 656       Double_t sp, p = track->GetTRDmomentum(ipl, &sp);
 657       SetTRDmomentum(p, ipl, &sp);
 658     }
 659   }
 660   //
 661   fTRDncls = track->GetTRDncls();
 662   fTRDntracklets &= 0xff & track->GetTRDntrackletsPID();
 663   fTRDchi2 = track->GetTRDchi2();
 664   //
 665   SetTOFsignal(track->GetTOFsignal());
 666   Double_t expt[AliPID::kSPECIESC];
 667   track->GetIntegratedTimes(expt,AliPID::kSPECIESC);
 668   SetIntegratedTimes(expt);
 669   //
 670   SetTrackPhiEtaPtOnEMCal(track->GetTrackPhiOnEMCal(),track->GetTrackEtaOnEMCal(),track->GetTrackPtOnEMCal());
 671   //
 672   SetLabel(track->GetLabel());
 673   // Set the status
 674   SetStatus(track->GetStatus());
 675   //
 676   // Set the ID
 677   SetID(track->GetID());
 678   //
 679 }
 680
 681 //_______________________________________________________________________
 682 AliESDtrack::AliESDtrack(TParticle * part) :
 683   AliExternalTrackParam(),
 684   fCp(0),
 685   fIp(0),
 686   fTPCInner(0),
 687   fOp(0),
 688   fHMPIDp(0),
 689   fFriendTrack(0),
 690   fTPCFitMap(159),//number of padrows
 691   fTPCClusterMap(159),//number of padrows
 692   fTPCSharedMap(159),//number of padrows
 693   fFlags(0),
 694   fID(0),
 695   fLabel(0),
 696   fITSLabel(0),
 697   fTPCLabel(0),
 698   fTRDLabel(0),
 699   fTOFLabel(NULL),
 700   fTOFCalChannel(-1),
 701   fTOFindex(-1),
 702   fHMPIDqn(0),
 703   fHMPIDcluIdx(-1),
 704   fCaloIndex(kEMCALNoMatch),
 705   fR(0),
 706   fITSr(0),
 707   fTPCr(0),
 708   fTRDr(0),
 709   fTOFr(0),
 710   fHMPIDr(0),
 711   fHMPIDtrkTheta(0),
 712   fHMPIDtrkPhi(0),
 713   fHMPIDsignal(0),
 714   fTrackTime(NULL),
 715   fTrackLength(0),
 716   fdTPC(0),fzTPC(0),
 717   fCddTPC(0),fCdzTPC(0),fCzzTPC(0),
 718   fCchi2TPC(0),
 719   fD(0),fZ(0),
 720   fCdd(0),fCdz(0),fCzz(0),
 721   fCchi2(0),
 722   fITSchi2(0),
 723   fTPCchi2(0),
 724   fTPCchi2Iter1(0),
 725   fTRDchi2(0),
 726   fTOFchi2(0),
 727   fHMPIDchi2(0),
 728   fGlobalChi2(0),
 729   fITSsignal(0),
 730   fTPCsignal(0),
 731   fTPCsignalTuned(0),
 732   fTPCsignalS(0),
 733   fTPCdEdxInfo(0),
 734   fTRDsignal(0),
 735   fTRDQuality(0),
 736   fTRDBudget(0),
 737   fTOFsignal(99999),
 738   fTOFsignalTuned(99999),
 739   fTOFsignalToT(99999),
 740   fTOFsignalRaw(99999),
 741   fTOFsignalDz(999),
 742   fTOFsignalDx(999),
 743   fTOFdeltaBC(999),
 744   fTOFl0l1(999),
 745   fCaloDx(0),
 746   fCaloDz(0),
 747   fHMPIDtrkX(0),
 748   fHMPIDtrkY(0),
 749   fHMPIDmipX(0),
 750   fHMPIDmipY(0),
 751   fTPCncls(0),
 752   fTPCnclsF(0),
 753   fTPCsignalN(0),
 754   fTPCnclsIter1(0),
 755   fTPCnclsFIter1(0),
 756   fITSncls(0),
 757   fITSClusterMap(0),
 758   fITSSharedMap(0),
 759   fTRDncls(0),
 760   fTRDncls0(0),
 761   fTRDntracklets(0),
 762   fTRDNchamberdEdx(0),
 763   fTRDNclusterdEdx(0),
 764   fTRDnSlices(0),
 765   fTRDslices(0x0),
 766   fVertexID(-2),  // -2 means an orphan track
 767   fPIDForTracking(AliPID::kPion),
 768   fESDEvent(0),
 769   fCacheNCrossedRows(-10),
 770   fCacheChi2TPCConstrainedVsGlobal(-10),
 771   fCacheChi2TPCConstrainedVsGlobalVertex(0),
 772   fDetectorPID(0x0),
 773   fTrackPhiOnEMCal(-999),
 774   fTrackEtaOnEMCal(-999),
 775   fTrackPtOnEMCal(-999),
 776   fNtofClusters(0),
 777   fTOFcluster(NULL)
 778 {
 779   //
 780   // ESD track from TParticle
 781   //
 782
 783   // Reset all the arrays
 784   Int_t i;
 785   for (i=kNITSchi2Std;i--;) fITSchi2Std[i] = 0;
 786
 787   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
 788   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
 789   for (i=0;i<kTRDnPlanes;i++) {
 790     fTRDTimBin[i]=0;
 791   }
 792   for (i=0;i<4;i++) {fITSdEdxSamples[i]=0.;}
 793   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
 794   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
 795   for (i=0;i<12;i++) {fITSModule[i]=-1;}
 796
 797   // Calculate the AliExternalTrackParam content
 798
 799   Double_t xref;
 800   Double_t alpha;
 801   Double_t param[5];
 802   Double_t covar[15];
 803
 804   // Calculate alpha: the rotation angle of the corresponding local system (TPC sector)
 805   alpha = part->Phi()*180./TMath::Pi();
 806   if (alpha<0) alpha+= 360.;
 807   if (alpha>360) alpha -= 360.;
 808
 809   Int_t sector = (Int_t)(alpha/20.);
 810   alpha = 10. + 20.*sector;
 811   alpha /= 180;
 812   alpha *= TMath::Pi();
 813
 814   // Covariance matrix: no errors, the parameters are exact
 815   for (i=0; i<15; i++) covar[i]=0.;
 816
 817   // Get the vertex of origin and the momentum
 818   TVector3 ver(part->Vx(),part->Vy(),part->Vz());
 819   TVector3 mom(part->Px(),part->Py(),part->Pz());
 820
 821   // Rotate to the local coordinate system (TPC sector)
 822   ver.RotateZ(-alpha);
 823   mom.RotateZ(-alpha);
 824
 825   // X of the referense plane
 826   xref = ver.X();
 827
 828   Int_t pdgCode = part->GetPdgCode();
 829
 830   Double_t charge =
 831     TDatabasePDG::Instance()->GetParticle(pdgCode)->Charge();
 832
 833   param[0] = ver.Y();
 834   param[1] = ver.Z();
 835   param[2] = TMath::Sin(mom.Phi());
 836   param[3] = mom.Pz()/mom.Pt();
 837   param[4] = TMath::Sign(1/mom.Pt(),charge);
 838
 839   // Set AliExternalTrackParam
 840   Set(xref, alpha, param, covar);
 841
 842   // Set the PID
 843   Int_t indexPID = 99;
 844   if (pdgCode<0) pdgCode = -pdgCode;
 845   for (i=0;i<AliPID::kSPECIESC;i++) if (pdgCode==AliPID::ParticleCode(i)) {indexPID = i; break;}
 846
 847   if (indexPID < AliPID::kSPECIESC) fPIDForTracking = indexPID;
 848
 849   // AliESD track label
 850   SetLabel(part->GetUniqueID());
 851
 852 }
 853
 854 //_______________________________________________________________________
 855 AliESDtrack::~AliESDtrack(){
 856   //
 857   // This is destructor according Coding Conventrions
 858   //
 859   //printf("Delete track\n");
 860   delete fIp;
 861   delete fTPCInner;
 862   delete fOp;
 863   delete fHMPIDp;
 864   delete fCp;
 865   delete fFriendTrack;
 866   delete fTPCdEdxInfo;
 867   if(fTRDnSlices)
 868     delete[] fTRDslices;
 869
 870   //Reset cached values - needed for TClonesArray in AliESDInputHandler
 871   fCacheNCrossedRows = -10.;
 872   fCacheChi2TPCConstrainedVsGlobal = -10.;
 873   if(fCacheChi2TPCConstrainedVsGlobalVertex) fCacheChi2TPCConstrainedVsGlobalVertex = 0;
 874
 875   if(fTOFcluster)
 876     delete[] fTOFcluster;
 877   fTOFcluster = NULL;
 878   fNtofClusters=0;
 879
 880   delete fDetectorPID;
 881
 882   delete[] fR;
 883   delete[] fITSr;
 884   delete[] fTPCr;
 885   delete[] fTRDr;
 886   delete[] fTOFr;
 887   delete[] fHMPIDr;
 888   //
 889   if(fTrackTime) delete[] fTrackTime;
 890   if(fTOFLabel) delete[] fTOFLabel;
 891 }
 892
 893 //_______________________________________________________________________
 894 AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source)
 895 {
 896   // == operator
 897
 898   if(&source == this) return *this;
 899   AliExternalTrackParam::operator=(source);
 900
 901
 902   if(source.fCp){
 903     // we have the trackparam: assign or copy construct
 904     if(fCp)*fCp = *source.fCp;
 905     else fCp = new AliExternalTrackParam(*source.fCp);
 906   }
 907   else{
 908     // no track param delete the old one
 909     delete fCp;
 910     fCp = 0;
 911   }
 912
 913   if(source.fIp){
 914     // we have the trackparam: assign or copy construct
 915     if(fIp)*fIp = *source.fIp;
 916     else fIp = new AliExternalTrackParam(*source.fIp);
 917   }
 918   else{
 919     // no track param delete the old one
 920     delete fIp;
 921     fIp = 0;
 922   }
 923
 924
 925   if(source.fTPCInner){
 926     // we have the trackparam: assign or copy construct
 927     if(fTPCInner) *fTPCInner = *source.fTPCInner;
 928     else fTPCInner = new AliExternalTrackParam(*source.fTPCInner);
 929   }
 930   else{
 931     // no track param delete the old one
 932     delete fTPCInner;
 933     fTPCInner = 0;
 934   }
 935
 936   if(source.fTPCdEdxInfo) {
 937     if(fTPCdEdxInfo) *fTPCdEdxInfo = *source.fTPCdEdxInfo;
 938     fTPCdEdxInfo = new AliTPCdEdxInfo(*source.fTPCdEdxInfo);
 939   }
 940
 941   if(source.fOp){
 942     // we have the trackparam: assign or copy construct
 943     if(fOp) *fOp = *source.fOp;
 944     else fOp = new AliExternalTrackParam(*source.fOp);
 945   }
 946   else{
 947     // no track param delete the old one
 948     delete fOp;
 949     fOp = 0;
 950   }
 951
 952
 953   if(source.fHMPIDp){
 954     // we have the trackparam: assign or copy construct
 955     if(fHMPIDp) *fHMPIDp = *source.fHMPIDp;
 956     else fHMPIDp = new AliExternalTrackParam(*source.fHMPIDp);
 957   }
 958   else{
 959     // no track param delete the old one
 960     delete fHMPIDp;
 961     fHMPIDp = 0;
 962   }
 963
 964   // copy also the friend track
 965   // use copy constructor
 966   if(source.fFriendTrack){
 967     // we have the trackparam: assign or copy construct
 968     delete fFriendTrack; fFriendTrack=new AliESDfriendTrack(*source.fFriendTrack);
 969   }
 970   else{
 971     // no track param delete the old one
 972     delete fFriendTrack; fFriendTrack= 0;
 973   }
 974
 975   fTPCFitMap = source.fTPCFitMap;
 976   fTPCClusterMap = source.fTPCClusterMap;
 977   fTPCSharedMap  = source.fTPCSharedMap;
 978   // the simple stuff
 979   fFlags    = source.fFlags;
 980   fID       = source.fID;
 981   fLabel    = source.fLabel;
 982   fITSLabel = source.fITSLabel;
 983   for(int i = 0; i< 12;++i){
 984     fITSModule[i] = source.fITSModule[i];
 985   }
 986   fTPCLabel = source.fTPCLabel;
 987   fTRDLabel = source.fTRDLabel;
 988   if(source.fTOFLabel){
 989     if(!fTOFLabel) fTOFLabel = new Int_t[3];
 990     for(int i = 0; i< 3;++i){
 991       fTOFLabel[i] = source.fTOFLabel[i];
 992     }
 993   }
 994   fTOFCalChannel = source.fTOFCalChannel;
 995   fTOFindex      = source.fTOFindex;
 996   fHMPIDqn       = source.fHMPIDqn;
 997   fHMPIDcluIdx   = source.fHMPIDcluIdx;
 998   fCaloIndex    = source.fCaloIndex;
 999   for (int i=kNITSchi2Std;i--;) fITSchi2Std[i] = source.fITSchi2Std[i];
1000   for(int i = 0; i< 3;++i){
1001     fKinkIndexes[i] = source.fKinkIndexes[i];
1002     fV0Indexes[i]   = source.fV0Indexes[i];
1003   }
1004
1005   if (source.fR) {
1006     if (!fR) fR = new Double32_t[AliPID::kSPECIES];
1007     for (Int_t i=AliPID::kSPECIES;i--;)  fR[i]=source.fR[i];
1008   }
1009   else {delete[] fR; fR = 0;}
1010
1011   if (source.fITSr) {
1012     if (!fITSr) fITSr = new Double32_t[AliPID::kSPECIES];
1013     for (Int_t i=AliPID::kSPECIES;i--;)  fITSr[i]=source.fITSr[i];
1014   }
1015   else {delete[] fITSr; fITSr = 0;}
1016   //
1017   if (source.fTPCr) {
1018     if (!fTPCr) fTPCr = new Double32_t[AliPID::kSPECIES];
1019     for (Int_t i=AliPID::kSPECIES;i--;) fTPCr[i]=source.fTPCr[i];
1020   }
1021   else {delete[] fTPCr; fTPCr = 0;}
1022
1023   if (source.fTRDr) {
1024     if (!fTRDr) fTRDr  = new Double32_t[AliPID::kSPECIES];
1025     for (Int_t i=AliPID::kSPECIES;i--;) fTRDr[i]=source.fTRDr[i];
1026   }
1027   else {delete[] fTRDr; fTRDr = 0;}
1028
1029   if (source.fTOFr) {
1030     if (!fTOFr) fTOFr = new Double32_t[AliPID::kSPECIES];
1031     for (Int_t i=AliPID::kSPECIES;i--;) fTOFr[i]=source.fTOFr[i];
1032   }
1033   else {delete[] fTOFr; fTOFr = 0;}
1034
1035   if (source.fHMPIDr) {
1036     if (!fHMPIDr) fHMPIDr = new Double32_t[AliPID::kSPECIES];
1037     for (Int_t i=AliPID::kSPECIES;i--;) fHMPIDr[i]=source.fHMPIDr[i];
1038   }
1039   else {delete[] fHMPIDr; fHMPIDr = 0;}
1040
1041   fPIDForTracking = source.fPIDForTracking;
1042
1043   fHMPIDtrkTheta = source.fHMPIDtrkTheta;
1044   fHMPIDtrkPhi   = source.fHMPIDtrkPhi;
1045   fHMPIDsignal   = source.fHMPIDsignal;
1046
1047
1048   if(fTrackTime){
1049     delete[] fTrackTime;
1050   }
1051   if(source.fTrackTime){
1052     fTrackTime = new Double32_t[AliPID::kSPECIESC];
1053     for(Int_t i=0;i < AliPID::kSPECIESC;i++)
1054       fTrackTime[i] = source.fTrackTime[i];
1055   }
1056
1057   fTrackLength   = source. fTrackLength;
1058   fdTPC  = source.fdTPC;
1059   fzTPC  = source.fzTPC;
1060   fCddTPC = source.fCddTPC;
1061   fCdzTPC = source.fCdzTPC;
1062   fCzzTPC = source.fCzzTPC;
1063   fCchi2TPC = source.fCchi2TPC;
1064
1065   fD  = source.fD;
1066   fZ  = source.fZ;
1067   fCdd = source.fCdd;
1068   fCdz = source.fCdz;
1069   fCzz = source.fCzz;
1070   fCchi2     = source.fCchi2;
1071
1072   fITSchi2   = source.fITSchi2;
1073   fTPCchi2   = source.fTPCchi2;
1074   fTPCchi2Iter1   = source.fTPCchi2Iter1;
1075   fTRDchi2   = source.fTRDchi2;
1076   fTOFchi2   = source.fTOFchi2;
1077   fHMPIDchi2 = source.fHMPIDchi2;
1078
1079   fGlobalChi2 = source.fGlobalChi2;
1080
1081   fITSsignal  = source.fITSsignal;
1082   for (Int_t i=0;i<4;i++) {fITSdEdxSamples[i]=source.fITSdEdxSamples[i];}
1083   fTPCsignal  = source.fTPCsignal;
1084   fTPCsignalTuned  = source.fTPCsignalTuned;
1085   fTPCsignalS = source.fTPCsignalS;
1086   for(int i = 0; i< 4;++i){
1087     fTPCPoints[i] = source.fTPCPoints[i];
1088   }
1089   fTRDsignal = source.fTRDsignal;
1090   fTRDNchamberdEdx = source.fTRDNchamberdEdx;
1091   fTRDNclusterdEdx = source.fTRDNclusterdEdx;
1092
1093   for(int i = 0;i < kTRDnPlanes;++i){
1094     fTRDTimBin[i] = source.fTRDTimBin[i];
1095   }
1096
1097   if(fTRDnSlices)
1098     delete[] fTRDslices;
1099   fTRDslices=0;
1100   fTRDnSlices=source.fTRDnSlices;
1101   if (fTRDnSlices) {
1102     fTRDslices=new Double32_t[fTRDnSlices];
1103     for(int j = 0;j < fTRDnSlices;++j) fTRDslices[j] = source.fTRDslices[j];
1104   }
1105
1106   fTRDQuality =   source.fTRDQuality;
1107   fTRDBudget  =   source.fTRDBudget;
1108   fTOFsignal  =   source.fTOFsignal;
1109   fTOFsignalTuned  = source.fTOFsignalTuned;
1110   fTOFsignalToT = source.fTOFsignalToT;
1111   fTOFsignalRaw = source.fTOFsignalRaw;
1112   fTOFsignalDz  = source.fTOFsignalDz;
1113   fTOFsignalDx  = source.fTOFsignalDx;
1114   fTOFdeltaBC   = source.fTOFdeltaBC;
1115   fTOFl0l1      = source.fTOFl0l1;
1116
1117   for(int i = 0;i<10;++i){
1118     fTOFInfo[i] = source.fTOFInfo[i];
1119   }
1120
1121   fHMPIDtrkX = source.fHMPIDtrkX;
1122   fHMPIDtrkY = source.fHMPIDtrkY;
1123   fHMPIDmipX = source.fHMPIDmipX;
1124   fHMPIDmipY = source.fHMPIDmipY;
1125
1126   fTPCncls    = source.fTPCncls;
1127   fTPCnclsF   = source.fTPCnclsF;
1128   fTPCsignalN = source.fTPCsignalN;
1129   fTPCnclsIter1    = source.fTPCnclsIter1;
1130   fTPCnclsFIter1   = source.fTPCnclsFIter1;
1131
1132   fITSncls = source.fITSncls;
1133   fITSClusterMap = source.fITSClusterMap;
1134   fITSSharedMap = source.fITSSharedMap;
1135   fTRDncls   = source.fTRDncls;
1136   fTRDncls0  = source.fTRDncls0;
1137   fTRDntracklets  = source.fTRDntracklets;
1138   fVertexID = source.fVertexID;
1139   fPIDForTracking = source.fPIDForTracking;
1140
1141   fCacheNCrossedRows = source.fCacheNCrossedRows;
1142   fCacheChi2TPCConstrainedVsGlobal = source.fCacheChi2TPCConstrainedVsGlobal;
1143   fCacheChi2TPCConstrainedVsGlobalVertex = source.fCacheChi2TPCConstrainedVsGlobalVertex;
1144
1145   delete fDetectorPID;
1146   fDetectorPID=0x0;
1147   if (source.fDetectorPID) fDetectorPID = new AliDetectorPID(*source.fDetectorPID);
1148
1149   fTrackPhiOnEMCal= source.fTrackPhiOnEMCal;
1150   fTrackEtaOnEMCal= source.fTrackEtaOnEMCal;
1151   fTrackPtOnEMCal= source.fTrackPtOnEMCal;
1152
1153   if(fTOFcluster){
1154     delete[] fTOFcluster;
1155   }
1156   fNtofClusters = source.fNtofClusters;
1157   if(fNtofClusters > 0){
1158     fTOFcluster = new Int_t[fNtofClusters];
1159         for(Int_t i=0;i < fNtofClusters;i++) fTOFcluster[i] = source.fTOFcluster[i];
1160   }
1161
1162   return *this;
1163 }
1164
1165
1166
1167 void AliESDtrack::Copy(TObject &obj) const {
1168
1169   // this overwrites the virtual TOBject::Copy()
1170   // to allow run time copying without casting
1171   // in AliESDEvent
1172
1173   if(this==&obj)return;
1174   AliESDtrack *robj = dynamic_cast<AliESDtrack*>(&obj);
1175   if(!robj)return; // not an AliESDtrack
1176   *robj = *this;
1177
1178 }
1179
1180
1181
1182 void AliESDtrack::AddCalibObject(TObject * object){
1183   //
1184   // add calib object to the list
1185   //
1186   if (!fFriendTrack) fFriendTrack  = new AliESDfriendTrack;
1187   if (!fFriendTrack) return;
1188   fFriendTrack->AddCalibObject(object);
1189 }
1190
1191 TObject *  AliESDtrack::GetCalibObject(Int_t index){
1192   //
1193   // return calib objct at given position
1194   //
1195   if (!fFriendTrack) return 0;
1196   return fFriendTrack->GetCalibObject(index);
1197 }
1198
1199
1200 Bool_t AliESDtrack::FillTPCOnlyTrack(AliESDtrack &track){
1201
1202   // Fills the information of the TPC-only first reconstruction pass
1203   // into the passed ESDtrack object. For consistency fTPCInner is also filled
1204   // again
1205
1206
1207
1208   // For data produced before r26675
1209   // RelateToVertexTPC was not properly called during reco
1210   // so you'll have to call it again, before FillTPCOnlyTrack
1211   //  Float_t p[2],cov[3];
1212   // track->GetImpactParametersTPC(p,cov);
1213   // if(p[0]==0&&p[1]==0) // <- Default values
1214   //  track->RelateToVertexTPC(esd->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
1215
1216
1217   if(!fTPCInner)return kFALSE;
1218
1219   // fill the TPC track params to the global track parameters
1220   track.Set(fTPCInner->GetX(),fTPCInner->GetAlpha(),fTPCInner->GetParameter(),fTPCInner->GetCovariance());
1221   track.fD = fdTPC;
1222   track.fZ = fzTPC;
1223   track.fCdd = fCddTPC;
1224   track.fCdz = fCdzTPC;
1225   track.fCzz = fCzzTPC;
1226
1227   // copy the inner params
1228   if(track.fIp) *track.fIp = *fIp;
1229   else track.fIp = new AliExternalTrackParam(*fIp);
1230
1231   // copy the TPCinner parameters
1232   if(track.fTPCInner) *track.fTPCInner = *fTPCInner;
1233   else track.fTPCInner = new AliExternalTrackParam(*fTPCInner);
1234   track.fdTPC   = fdTPC;
1235   track.fzTPC   = fzTPC;
1236   track.fCddTPC = fCddTPC;
1237   track.fCdzTPC = fCdzTPC;
1238   track.fCzzTPC = fCzzTPC;
1239   track.fCchi2TPC = fCchi2TPC;
1240
1241   // copy all other TPC specific parameters
1242
1243   // replace label by TPC label
1244   track.fLabel    = fTPCLabel;
1245   track.fTPCLabel = fTPCLabel;
1246
1247   track.fTPCchi2 = fTPCchi2;
1248   track.fTPCchi2Iter1 = fTPCchi2Iter1;
1249   track.fTPCsignal = fTPCsignal;
1250   track.fTPCsignalTuned = fTPCsignalTuned;
1251   track.fTPCsignalS = fTPCsignalS;
1252   for(int i = 0;i<4;++i)track.fTPCPoints[i] = fTPCPoints[i];
1253
1254   track.fTPCncls    = fTPCncls;
1255   track.fTPCnclsF   = fTPCnclsF;
1256   track.fTPCsignalN =  fTPCsignalN;
1257   track.fTPCnclsIter1    = fTPCnclsIter1;
1258   track.fTPCnclsFIter1   = fTPCnclsFIter1;
1259
1260   // PID
1261   if (fTPCr) {
1262     if (!track.fTPCr) track.fTPCr = new Double32_t[AliPID::kSPECIES];
1263     for(int i=AliPID::kSPECIES;i--;) track.fTPCr[i] = fTPCr[i];
1264   }
1265   //
1266   if (fR) {
1267     if (!track.fR) track.fR = new Double32_t[AliPID::kSPECIES];
1268     for(int i=AliPID::kSPECIES;i--;) track.fR[i] = fR[i];
1269   }
1270
1271   track.fTPCFitMap = fTPCFitMap;
1272   track.fTPCClusterMap = fTPCClusterMap;
1273   track.fTPCSharedMap = fTPCSharedMap;
1274
1275
1276   // reset the flags
1277   track.fFlags = kTPCin;
1278   track.fID    = fID;
1279
1280   track.fFlags |= fFlags & kTPCpid; //copy the TPCpid status flag
1281
1282   for (Int_t i=0;i<3;i++) track.fKinkIndexes[i] = fKinkIndexes[i];
1283
1284   return kTRUE;
1285
1286 }
1287
1288 //_______________________________________________________________________
1289 void AliESDtrack::MakeMiniESDtrack(){
1290   // Resets everything except
1291   // fFlags: Reconstruction status flags
1292   // fLabel: Track label
1293   // fID:  Unique ID of the track
1294   // Impact parameter information
1295   // fR[AliPID::kSPECIES]: combined "detector response probability"
1296   // Running track parameters in the base class (AliExternalTrackParam)
1297
1298   fTrackLength = 0;
1299
1300   if(fTrackTime)
1301     for (Int_t i=0;i<AliPID::kSPECIESC;i++) fTrackTime[i] = 0;
1302
1303   // Reset track parameters constrained to the primary vertex
1304   delete fCp;fCp = 0;
1305
1306   // Reset track parameters at the inner wall of TPC
1307   delete fIp;fIp = 0;
1308   delete fTPCInner;fTPCInner=0;
1309   // Reset track parameters at the inner wall of the TRD
1310   delete fOp;fOp = 0;
1311   // Reset track parameters at the HMPID
1312   delete fHMPIDp;fHMPIDp = 0;
1313
1314
1315   // Reset ITS track related information
1316   fITSchi2 = 0;
1317   fITSncls = 0;
1318   fITSClusterMap=0;
1319   fITSSharedMap=0;
1320   fITSsignal = 0;
1321   for (Int_t i=0;i<4;i++) fITSdEdxSamples[i] = 0.;
1322   if (fITSr) for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=0;
1323   fITSLabel = 0;
1324
1325   // Reset TPC related track information
1326   fTPCchi2 = 0;
1327   fTPCchi2Iter1 = 0;
1328   fTPCncls = 0;
1329   fTPCnclsF = 0;
1330   fTPCnclsIter1 = 0;
1331   fTPCnclsFIter1 = 0;
1332   fTPCFitMap = 0;
1333   fTPCClusterMap = 0;
1334   fTPCSharedMap = 0;
1335   fTPCsignal= 0;
1336   fTPCsignalTuned= 0;
1337   fTPCsignalS= 0;
1338   fTPCsignalN= 0;
1339   if (fTPCr) for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i] = 0;
1340   fTPCLabel=0;
1341   for (Int_t i=0;i<4;i++) fTPCPoints[i] = 0;
1342   for (Int_t i=0; i<3;i++)   fKinkIndexes[i] = 0;
1343   for (Int_t i=0; i<3;i++)   fV0Indexes[i] = 0;
1344
1345   // Reset TRD related track information
1346   fTRDchi2 = 0;
1347   fTRDncls = 0;
1348   fTRDncls0 = 0;
1349   fTRDsignal = 0;
1350   fTRDNchamberdEdx = 0;
1351   fTRDNclusterdEdx = 0;
1352
1353   for (Int_t i=0;i<kTRDnPlanes;i++) {
1354     fTRDTimBin[i]  = 0;
1355   }
1356   if (fTRDr) for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0;
1357   fTRDLabel = 0;
1358   fTRDQuality  = 0;
1359   fTRDntracklets = 0;
1360   if(fTRDnSlices)
1361     delete[] fTRDslices;
1362   fTRDslices=0x0;
1363   fTRDnSlices=0;
1364   fTRDBudget  = 0;
1365
1366   // Reset TOF related track information
1367   fTOFchi2 = 0;
1368   fTOFindex = -1;
1369   fTOFsignal = 99999;
1370   fTOFCalChannel = -1;
1371   fTOFsignalToT = 99999;
1372   fTOFsignalRaw = 99999;
1373   fTOFsignalDz = 999;
1374   fTOFsignalDx = 999;
1375   fTOFdeltaBC = 999;
1376   fTOFl0l1 = 999;
1377   if (fTOFr) for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i] = 0;
1378   for (Int_t i=0;i<10;i++) fTOFInfo[i] = 0;
1379
1380   // Reset HMPID related track information
1381   fHMPIDchi2 = 0;
1382   fHMPIDqn = 0;
1383   fHMPIDcluIdx = -1;
1384   fHMPIDsignal = 0;
1385   if (fHMPIDr) for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i] = 0;
1386   fHMPIDtrkTheta = 0;
1387   fHMPIDtrkPhi = 0;
1388   fHMPIDtrkX = 0;
1389   fHMPIDtrkY = 0;
1390   fHMPIDmipX = 0;
1391   fHMPIDmipY = 0;
1392   fCaloIndex = kEMCALNoMatch;
1393
1394   // reset global track chi2
1395   fGlobalChi2 = 0;
1396
1397   fVertexID = -2; // an orphan track
1398   fPIDForTracking = AliPID::kPion;
1399   //
1400   delete fFriendTrack; fFriendTrack = 0;
1401 }
1402
1403 //_______________________________________________________________________
1404 Int_t AliESDtrack::GetPID(Bool_t tpcOnly) const
1405 {
1406   // Returns the particle most probable id. For backward compatibility first the prob. arrays
1407   // will be checked, but normally the GetPIDForTracking will be returned
1408   Int_t i;
1409   const Double32_t *prob = 0;
1410   if (tpcOnly) { // check if TPCpid is valid
1411     if (!fTPCr) return GetPIDForTracking();
1412     prob = fTPCr;
1413     for (i=0; i<AliPID::kSPECIES-1; i++) if (prob[i] != prob[i+1]) break;
1414     if (i == AliPID::kSPECIES-1) prob = 0; // not valid, try with combined pid
1415   }
1416   if (!prob) { // either requested TPCpid is not valid or comb.pid is requested
1417     if (!fR) return GetPIDForTracking();
1418     prob = fR;
1419     for (i=0; i<AliPID::kSPECIES-1; i++) if (prob[i] != prob[i+1]) break;
1420     if (i == AliPID::kSPECIES-1) return GetPIDForTracking();  // If all the probabilities are equal, return the pion mass
1421   }
1422   //
1423   Float_t max=0.;
1424   Int_t k=-1;
1425   for (i=0; i<AliPID::kSPECIES; i++) if (prob[i]>max) {k=i; max=prob[i];}
1426   //
1427   if (k==0) { // dE/dx "crossing points" in the TPC
1428     /*
1429     Double_t p=GetP();
1430     if ((p>0.38)&&(p<0.48))
1431       if (prob[0]<prob[3]*10.) return AliPID::kKaon;
1432     if ((p>0.75)&&(p<0.85))
1433       if (prob[0]<prob[4]*10.) return AliPID::kProton;
1434     */
1435     return AliPID::kElectron;
1436   }
1437   if (k==1) return AliPID::kMuon;
1438   if (k==2||k==-1) return AliPID::kPion;
1439   if (k==3) return AliPID::kKaon;
1440   if (k==4) return AliPID::kProton;
1441   //  AliWarning("Undefined PID !");
1442   return GetPIDForTracking();
1443 }
1444
1445 //_______________________________________________________________________
1446 Int_t AliESDtrack::GetTOFBunchCrossing(Double_t b, Bool_t pidTPConly) const
1447 {
1448   // Returns the number of bunch crossings after trigger (assuming 25ns spacing)
1449   const double kSpacing = 25e3; // min interbanch spacing
1450   const double kShift = 0;
1451   Int_t bcid = kTOFBCNA; // defualt one
1452   if (!IsOn(kTOFout)/* || !IsOn(kESDpid)*/) return bcid; // no info
1453   //
1454   double tdif = GetTOFsignal();
1455   if (IsOn(kTIME)) { // integrated time info is there
1456     int pid = GetPID(pidTPConly);
1457     Double_t times[AliPID::kSPECIESC];
1458     // old esd has only AliPID::kSPECIES times
1459     GetIntegratedTimes(times,pid>=AliPID::kSPECIES ? AliPID::kSPECIESC : AliPID::kSPECIES);
1460     tdif -= times[pid];
1461   }
1462   else { // assume integrated time info from TOF radius and momentum
1463     const double kRTOF = 385.;
1464     const double kCSpeed = 3.e-2; // cm/ps
1465     double p = GetP();
1466     if (p<0.01) return bcid;
1467     double m = GetMass(pidTPConly);
1468     double curv = GetC(b);
1469     double path = TMath::Abs(curv)>kAlmost0 ? // account for curvature
1470       2./curv*TMath::ASin(kRTOF*curv/2.)*TMath::Sqrt(1.+GetTgl()*GetTgl()) : kRTOF;
1471     tdif -= path/kCSpeed*TMath::Sqrt(1.+m*m/(p*p));
1472   }
1473   bcid = TMath::Nint((tdif - kShift)/kSpacing);
1474   return bcid;
1475 }
1476
1477 //______________________________________________________________________________
1478 Double_t AliESDtrack::M() const
1479 {
1480   // Returns the assumed mass
1481   // (the pion mass, if the particle can't be identified properly).
1482   static Bool_t printerr=kTRUE;
1483   if (printerr) {
1484      AliWarning("WARNING !!! ... THIS WILL BE PRINTED JUST ONCE !!!");
1485      printerr = kFALSE;
1486      AliWarning("This is the ESD mass. Use it with care !");
1487   }
1488   return GetMass();
1489 }
1490
1491 //______________________________________________________________________________
1492 Double_t AliESDtrack::E() const
1493 {
1494   // Returns the energy of the particle given its assumed mass.
1495   // Assumes the pion mass if the particle can't be identified properly.
1496
1497   Double_t m = M();
1498   Double_t p = P();
1499   return TMath::Sqrt(p*p + m*m);
1500 }
1501
1502 //______________________________________________________________________________
1503 Double_t AliESDtrack::Y() const
1504 {
1505   // Returns the rapidity of a particle given its assumed mass.
1506   // Assumes the pion mass if the particle can't be identified properly.
1507
1508   Double_t e = E();
1509   Double_t pz = Pz();
1510   if (e != TMath::Abs(pz)) { // energy was not equal to pz
1511     return 0.5*TMath::Log((e+pz)/(e-pz));
1512   } else { // energy was equal to pz
1513     return -999.;
1514   }
1515 }
1516
1517 //_______________________________________________________________________
1518 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
1519   //
1520   // This function updates track's running parameters
1521   //
1522   Bool_t rc=kTRUE;
1523
1524   SetStatus(flags);
1525   fLabel=t->GetLabel();
1526
1527   if (t->IsStartedTimeIntegral()) {
1528     SetStatus(kTIME);
1529     Double_t times[AliPID::kSPECIESC];
1530     t->GetIntegratedTimes(times);
1531     SetIntegratedTimes(times);
1532     SetIntegratedLength(t->GetIntegratedLength());
1533   }
1534
1535   Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1536   if (fFriendTrack) {
1537   if (flags==kITSout) fFriendTrack->SetITSOut(*t);
1538   if (flags==kTPCout) fFriendTrack->SetTPCOut(*t);
1539   if (flags==kTRDrefit) fFriendTrack->SetTRDIn(*t);
1540   }
1541
1542   switch (flags) {
1543
1544   case kITSin:
1545     fITSchi2Std[0] = t->GetChi2();
1546     //
1547   case kITSout:
1548     fITSchi2Std[1] = t->GetChi2();
1549   case kITSrefit:
1550     {
1551     fITSchi2Std[2] = t->GetChi2();
1552     fITSClusterMap=0;
1553     fITSncls=t->GetNumberOfClusters();
1554     if (fFriendTrack) {
1555     Int_t* indexITS = new Int_t[AliESDfriendTrack::kMaxITScluster];
1556     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
1557         indexITS[i]=t->GetClusterIndex(i);
1558
1559         if (i<fITSncls) {
1560           Int_t l=(indexITS[i] & 0xf0000000) >> 28;
1561            SETBIT(fITSClusterMap,l);
1562         }
1563     }
1564     fFriendTrack->SetITSIndices(indexITS,AliESDfriendTrack::kMaxITScluster);
1565     delete [] indexITS;
1566     }
1567
1568     fITSchi2=t->GetChi2();
1569     fITSsignal=t->GetPIDsignal();
1570     fITSLabel = t->GetLabel();
1571     // keep in fOp the parameters outside ITS for ITS stand-alone tracks
1572     if (flags==kITSout) {
1573       if (!fOp) fOp=new AliExternalTrackParam(*t);
1574       else
1575         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1576     }
1577     }
1578     break;
1579
1580   case kTPCin: case kTPCrefit:
1581     {
1582     fTPCLabel = t->GetLabel();
1583     if (flags==kTPCin)  {
1584       fTPCInner=new AliExternalTrackParam(*t);
1585       fTPCnclsIter1=t->GetNumberOfClusters();
1586       fTPCchi2Iter1=t->GetChi2();
1587     }
1588     if (!fIp) fIp=new AliExternalTrackParam(*t);
1589     else
1590       fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1591     }
1592     // Intentionally no break statement; need to set general TPC variables as well
1593   case kTPCout:
1594     {
1595     if (flags & kTPCout){
1596       if (!fOp) fOp=new AliExternalTrackParam(*t);
1597       else
1598         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1599     }
1600     fTPCncls=t->GetNumberOfClusters();
1601     fTPCchi2=t->GetChi2();
1602
1603     if (fFriendTrack) {  // Copy cluster indices
1604       Int_t* indexTPC = new Int_t[AliESDfriendTrack::kMaxTPCcluster];
1605       for (Int_t i=0;i<AliESDfriendTrack::kMaxTPCcluster;i++)
1606         indexTPC[i]=t->GetClusterIndex(i);
1607       fFriendTrack->SetTPCIndices(indexTPC,AliESDfriendTrack::kMaxTPCcluster);
1608       delete [] indexTPC;
1609     }
1610     fTPCsignal=t->GetPIDsignal();
1611     }
1612     break;
1613
1614   case kTRDin: case kTRDrefit:
1615     break;
1616   case kTRDout:
1617     {
1618     fTRDLabel = t->GetLabel();
1619     fTRDchi2  = t->GetChi2();
1620     fTRDncls  = t->GetNumberOfClusters();
1621     if (fFriendTrack) {
1622       Int_t* indexTRD = new Int_t[AliESDfriendTrack::kMaxTRDcluster];
1623       for (Int_t i=0;i<AliESDfriendTrack::kMaxTRDcluster;i++) indexTRD[i]=-2;
1624       for (Int_t i=0;i<6;i++) indexTRD[i]=t->GetTrackletIndex(i);
1625       fFriendTrack->SetTRDIndices(indexTRD,AliESDfriendTrack::kMaxTRDcluster);
1626       delete [] indexTRD;
1627     }
1628
1629     //commented out by Xianguo
1630     //fTRDsignal=t->GetPIDsignal();
1631     }
1632     break;
1633   case kTRDbackup:
1634     if (!fOp) fOp=new AliExternalTrackParam(*t);
1635     else
1636       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1637     fTRDncls0 = t->GetNumberOfClusters();
1638     break;
1639   case kTOFin:
1640     break;
1641   case kTOFout:
1642     break;
1643   case kTRDStop:
1644     break;
1645   case kHMPIDout:
1646   if (!fHMPIDp) fHMPIDp=new AliExternalTrackParam(*t);
1647     else
1648       fHMPIDp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
1649     break;
1650   default:
1651     AliError("Wrong flag !");
1652     return kFALSE;
1653   }
1654
1655   return rc;
1656 }
1657
1658 //_______________________________________________________________________
1659 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
1660   //---------------------------------------------------------------------
1661   // This function returns external representation of the track parameters
1662   //---------------------------------------------------------------------
1663   x=GetX();
1664   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
1665 }
1666
1667 //_______________________________________________________________________
1668 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
1669   //---------------------------------------------------------------------
1670   // This function returns external representation of the cov. matrix
1671   //---------------------------------------------------------------------
1672   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
1673 }
1674
1675 //_______________________________________________________________________
1676 Bool_t AliESDtrack::GetConstrainedExternalParameters
1677                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1678   //---------------------------------------------------------------------
1679   // This function returns the constrained external track parameters
1680   //---------------------------------------------------------------------
1681   if (!fCp) return kFALSE;
1682   alpha=fCp->GetAlpha();
1683   x=fCp->GetX();
1684   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
1685   return kTRUE;
1686 }
1687
1688 //_______________________________________________________________________
1689 Bool_t
1690 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
1691   //---------------------------------------------------------------------
1692   // This function returns the constrained external cov. matrix
1693   //---------------------------------------------------------------------
1694   if (!fCp) return kFALSE;
1695   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
1696   return kTRUE;
1697 }
1698
1699 Bool_t
1700 AliESDtrack::GetInnerExternalParameters
1701                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1702   //---------------------------------------------------------------------
1703   // This function returns external representation of the track parameters
1704   // at the inner layer of TPC
1705   //---------------------------------------------------------------------
1706   if (!fIp) return kFALSE;
1707   alpha=fIp->GetAlpha();
1708   x=fIp->GetX();
1709   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
1710   return kTRUE;
1711 }
1712
1713 Bool_t
1714 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
1715  //---------------------------------------------------------------------
1716  // This function returns external representation of the cov. matrix
1717  // at the inner layer of TPC
1718  //---------------------------------------------------------------------
1719   if (!fIp) return kFALSE;
1720   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
1721   return kTRUE;
1722 }
1723
1724 void
1725 AliESDtrack::SetOuterParam(const AliExternalTrackParam *p, ULong_t flags) {
1726   //
1727   // This is a direct setter for the outer track parameters
1728   //
1729   SetStatus(flags);
1730   if (fOp) delete fOp;
1731   fOp=new AliExternalTrackParam(*p);
1732 }
1733
1734 void
1735 AliESDtrack::SetOuterHmpParam(const AliExternalTrackParam *p, ULong_t flags) {
1736   //
1737   // This is a direct setter for the outer track parameters
1738   //
1739   SetStatus(flags);
1740   if (fHMPIDp) delete fHMPIDp;
1741   fHMPIDp=new AliExternalTrackParam(*p);
1742 }
1743
1744 Bool_t
1745 AliESDtrack::GetOuterExternalParameters
1746                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1747   //---------------------------------------------------------------------
1748   // This function returns external representation of the track parameters
1749   // at the inner layer of TRD
1750   //---------------------------------------------------------------------
1751   if (!fOp) return kFALSE;
1752   alpha=fOp->GetAlpha();
1753   x=fOp->GetX();
1754   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
1755   return kTRUE;
1756 }
1757
1758 Bool_t
1759 AliESDtrack::GetOuterHmpExternalParameters
1760                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
1761   //---------------------------------------------------------------------
1762   // This function returns external representation of the track parameters
1763   // at the inner layer of TRD
1764   //---------------------------------------------------------------------
1765   if (!fHMPIDp) return kFALSE;
1766   alpha=fHMPIDp->GetAlpha();
1767   x=fHMPIDp->GetX();
1768   for (Int_t i=0; i<5; i++) p[i]=fHMPIDp->GetParameter()[i];
1769   return kTRUE;
1770 }
1771
1772 Bool_t
1773 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
1774  //---------------------------------------------------------------------
1775  // This function returns external representation of the cov. matrix
1776  // at the inner layer of TRD
1777  //---------------------------------------------------------------------
1778   if (!fOp) return kFALSE;
1779   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
1780   return kTRUE;
1781 }
1782
1783 Bool_t
1784 AliESDtrack::GetOuterHmpExternalCovariance(Double_t cov[15]) const {
1785  //---------------------------------------------------------------------
1786  // This function returns external representation of the cov. matrix
1787  // at the inner layer of TRD
1788  //---------------------------------------------------------------------
1789   if (!fHMPIDp) return kFALSE;
1790   for (Int_t i=0; i<15; i++) cov[i]=fHMPIDp->GetCovariance()[i];
1791   return kTRUE;
1792 }
1793
1794 Int_t AliESDtrack::GetNcls(Int_t idet) const
1795 {
1796   // Get number of clusters by subdetector index
1797   //
1798   Int_t ncls = 0;
1799   switch(idet){
1800   case 0:
1801     ncls = fITSncls;
1802     break;
1803   case 1:
1804     ncls = fTPCncls;
1805     break;
1806   case 2:
1807     ncls = fTRDncls;
1808     break;
1809   case 3:
1810     if (fTOFindex != -1)
1811       ncls = 1;
1812     break;
1813   case 4: //PHOS
1814     break;
1815   case 5: //HMPID
1816     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1817       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1818         ncls = 1;
1819       }
1820     }
1821     break;
1822   default:
1823     break;
1824   }
1825   return ncls;
1826 }
1827
1828 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
1829 {
1830   // Get cluster index array by subdetector index
1831   //
1832   Int_t ncls = 0;
1833   switch(idet){
1834   case 0:
1835     ncls = GetITSclusters(idx);
1836     break;
1837   case 1:
1838     ncls = GetTPCclusters(idx);
1839     break;
1840   case 2:
1841     ncls = GetTRDclusters(idx);
1842     break;
1843   case 3:
1844     if (fTOFindex != -1) {
1845       idx[0] = fTOFindex;
1846       ncls = 1;
1847     }
1848     break;
1849   case 4: //PHOS
1850     break;
1851   case 5:
1852     if ((fHMPIDcluIdx >= 0) && (fHMPIDcluIdx < 7000000)) {
1853       if ((fHMPIDcluIdx%1000000 != 9999) && (fHMPIDcluIdx%1000000 != 99999)) {
1854         idx[0] = GetHMPIDcluIdx();
1855         ncls = 1;
1856       }
1857     }
1858     break;
1859   case 6: //EMCAL
1860     break;
1861   default:
1862     break;
1863   }
1864   return ncls;
1865 }
1866
1867 //_______________________________________________________________________
1868 void AliESDtrack::GetIntegratedTimes(Double_t *times, Int_t nspec) const
1869 {
1870   // get integrated time for requested N species
1871   if (nspec<1) return;
1872   if(fNtofClusters>0 && GetESDEvent()){
1873     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
1874     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
1875     //
1876     for(int i=tofcl->GetNMatchableTracks();i--;){
1877       if(tofcl->GetTrackIndex(i) == GetID()) {
1878         for (int j=nspec; j--;) times[j]=tofcl->GetIntegratedTime(j,i);
1879         return;
1880       }
1881     }
1882   }
1883   else if(fNtofClusters>0)
1884     AliInfo("No AliESDEvent available here!\n");
1885
1886   // Returns the array with integrated times for each particle hypothesis
1887   if(fTrackTime)
1888     for (int i=nspec; i--;) times[i]=fTrackTime[i];
1889   else
1890 // The line below is wrong since it does not honor the nspec value
1891 // The "times" array may have only AliPID::kSPECIES size, as called by:
1892 // AliESDpidCuts::AcceptTrack()
1893 //    for (int i=AliPID::kSPECIESC; i--;) times[i]=0.0;
1894     for (int i=nspec; i--;) times[i]=0.0;
1895 }
1896 //_______________________________________________________________________
1897 Double_t AliESDtrack::GetIntegratedLength() const{
1898   Int_t index = -1;
1899   if(fNtofClusters>0 && GetESDEvent()){
1900     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
1901     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
1902
1903     for(Int_t i=0;i < tofcl->GetNMatchableTracks();i++){
1904       if(tofcl->GetTrackIndex(i) == GetID()) index = i;
1905     }
1906
1907     if(fNtofClusters>0 && index > -1)
1908       return tofcl->GetLength(index);
1909   }
1910   else if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
1911
1912   return fTrackLength;
1913 }
1914
1915 //_______________________________________________________________________
1916 void AliESDtrack::SetIntegratedTimes(const Double_t *times)
1917 {
1918   // Sets the array with integrated times for each particle hypotesis
1919   if(!fTrackTime) fTrackTime = new Double32_t[AliPID::kSPECIESC];
1920   for (int i=AliPID::kSPECIESC; i--;) fTrackTime[i]=times[i];
1921 }
1922
1923 //_______________________________________________________________________
1924 void AliESDtrack::SetITSpid(const Double_t *p)
1925 {
1926   // Sets values for the probability of each particle type (in ITS)
1927   if (!fITSr) fITSr = new Double32_t[AliPID::kSPECIESC];
1928   SetPIDValues(fITSr,p,AliPID::kSPECIES);
1929   SetStatus(AliESDtrack::kITSpid);
1930 }
1931
1932 //_______________________________________________________________________
1933 void AliESDtrack::GetITSpid(Double_t *p) const {
1934   // Gets the probability of each particle type (in ITS)
1935   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i] = fITSr ? fITSr[i] : 0;
1936 }
1937
1938 //_______________________________________________________________________
1939 Char_t AliESDtrack::GetITSclusters(Int_t *idx) const {
1940   //---------------------------------------------------------------------
1941   // This function returns indices of the assgined ITS clusters
1942   //---------------------------------------------------------------------
1943   if (idx && fFriendTrack) {
1944     Int_t *index=fFriendTrack->GetITSindices();
1945     for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) {
1946       if ( (i>=fITSncls) && (i<6) ) idx[i]=-1;
1947       else {
1948         if (index) {
1949           idx[i]=index[i];
1950         }
1951         else idx[i]= -2;
1952       }
1953     }
1954   }
1955   return fITSncls;
1956 }
1957
1958 //_______________________________________________________________________
1959 Bool_t AliESDtrack::GetITSModuleIndexInfo(Int_t ilayer,Int_t &idet,Int_t &status,
1960                                          Float_t &xloc,Float_t &zloc) const {
1961   //----------------------------------------------------------------------
1962   // This function encodes in the module number also the status of cluster association
1963   // "status" can have the following values:
1964   // 1 "found" (cluster is associated),
1965   // 2 "dead" (module is dead from OCDB),
1966   // 3 "skipped" (module or layer forced to be skipped),
1967   // 4 "outinz" (track out of z acceptance),
1968   // 5 "nocls" (no clusters in the road),
1969   // 6 "norefit" (cluster rejected during refit),
1970   // 7 "deadzspd" (holes in z in SPD)
1971   // Also given are the coordinates of the crossing point of track and module
1972   // (in the local module ref. system)
1973   // WARNING: THIS METHOD HAS TO BE SYNCHRONIZED WITH AliITStrackV2::GetModuleIndexInfo()!
1974   //----------------------------------------------------------------------
1975
1976   if(fITSModule[ilayer]==-1) {
1977     idet = -1;
1978     status=0;
1979     xloc=-99.; zloc=-99.;
1980     return kFALSE;
1981   }
1982
1983   Int_t module = fITSModule[ilayer];
1984
1985   idet = Int_t(module/1000000);
1986
1987   module -= idet*1000000;
1988
1989   status = Int_t(module/100000);
1990
1991   module -= status*100000;
1992
1993   Int_t signs = Int_t(module/10000);
1994
1995   module-=signs*10000;
1996
1997   Int_t xInt = Int_t(module/100);
1998   module -= xInt*100;
1999
2000   Int_t zInt = module;
2001
2002   if(signs==1) { xInt*=1; zInt*=1; }
2003   if(signs==2) { xInt*=1; zInt*=-1; }
2004   if(signs==3) { xInt*=-1; zInt*=1; }
2005   if(signs==4) { xInt*=-1; zInt*=-1; }
2006
2007   xloc = 0.1*(Float_t)xInt;
2008   zloc = 0.1*(Float_t)zInt;
2009
2010   if(status==4) idet = -1;
2011
2012   return kTRUE;
2013 }
2014
2015 //_______________________________________________________________________
2016 UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
2017   //---------------------------------------------------------------------
2018   // This function returns indices of the assgined ITS clusters
2019   //---------------------------------------------------------------------
2020   if (idx && fFriendTrack) {
2021     Int_t *index=fFriendTrack->GetTPCindices();
2022
2023     if (index){
2024       for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
2025     }
2026     else {
2027       for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=-2;
2028     }
2029   }
2030   return fTPCncls;
2031 }
2032
2033 //_______________________________________________________________________
2034 Float_t AliESDtrack::GetTPCCrossedRows() const
2035 {
2036   // This function calls GetTPCClusterInfo with some default parameters which are used in the track selection and caches the outcome
2037   // because GetTPCClusterInfo is quite time-consuming
2038
2039   if (fCacheNCrossedRows > -1)
2040     return fCacheNCrossedRows;
2041
2042   fCacheNCrossedRows = GetTPCClusterInfo(2, 1);
2043   return fCacheNCrossedRows;
2044 }
2045
2046 //_______________________________________________________________________
2047 Float_t AliESDtrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1, Int_t bitType ) const
2048 {
2049   //
2050   // TPC cluster information
2051   // type 0: get fraction of found/findable clusters with neighbourhood definition
2052   //      1: findable clusters with neighbourhood definition
2053   //      2: found clusters
2054   // bitType:
2055   //      0 - all cluster used
2056   //      1 - clusters  used for the kalman update
2057   // definition of findable clusters:
2058   //            a cluster is defined as findable if there is another cluster
2059   //           within +- nNeighbours pad rows. The idea is to overcome threshold
2060   //           effects with a very simple algorithm.
2061   //
2062
2063
2064   Int_t found=0;
2065   Int_t findable=0;
2066   Int_t last=-nNeighbours;
2067   const TBits & clusterMap = (bitType%2==0) ? fTPCClusterMap : fTPCFitMap;
2068
2069   Int_t upperBound=clusterMap.GetNbits();
2070   if (upperBound>row1) upperBound=row1;
2071   for (Int_t i=row0; i<upperBound; ++i){
2072     //look to current row
2073     if (clusterMap[i]) {
2074       last=i;
2075       ++found;
2076       ++findable;
2077       continue;
2078     }
2079     //look to nNeighbours before
2080     if ((i-last)<=nNeighbours) {
2081       ++findable;
2082       continue;
2083     }
2084     //look to nNeighbours after
2085     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
2086       if (clusterMap[j]){
2087         ++findable;
2088         break;
2089       }
2090     }
2091   }
2092   if (type==2) return found;
2093   if (type==1) return findable;
2094
2095   if (type==0){
2096     Float_t fraction=0;
2097     if (findable>0)
2098       fraction=(Float_t)found/(Float_t)findable;
2099     else
2100       fraction=0;
2101     return fraction;
2102   }
2103   return 0;  // undefined type - default value
2104 }
2105
2106 //_______________________________________________________________________
2107 Float_t AliESDtrack::GetTPCClusterDensity(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1, Int_t bitType ) const
2108 {
2109   //
2110   // TPC cluster density -  only rows where signal before and after given row are used
2111   //                     -  slower function
2112   // type 0: get fraction of found/findable clusters with neighbourhood definition
2113   //      1: findable clusters with neighbourhood definition
2114   //      2: found clusters
2115   // bitType:
2116   //      0 - all cluster used
2117   //      1 - clusters  used for the kalman update
2118   // definition of findable clusters:
2119   //            a cluster is defined as findable if there is another cluster
2120   //           within +- nNeighbours pad rows. The idea is to overcome threshold
2121   //           effects with a very simple algorithm.
2122   //
2123   Int_t found=0;
2124   Int_t findable=0;
2125   //  Int_t last=-nNeighbours;
2126   const TBits & clusterMap = (bitType%2==0) ? fTPCClusterMap : fTPCFitMap;
2127   Int_t upperBound=clusterMap.GetNbits();
2128   if (upperBound>row1) upperBound=row1;
2129   for (Int_t i=row0; i<upperBound; ++i){
2130     Bool_t isUp=kFALSE;
2131     Bool_t isDown=kFALSE;
2132     for (Int_t idelta=1; idelta<=nNeighbours; idelta++){
2133       if (i-idelta>=0 && clusterMap[i-idelta]) isDown=kTRUE;
2134       if (i+idelta<upperBound && clusterMap[i+idelta]) isUp=kTRUE;
2135     }
2136     if (isUp&&isDown){
2137       ++findable;
2138       if (clusterMap[i]) ++found;
2139     }
2140   }
2141   if (type==2) return found;
2142   if (type==1) return findable;
2143
2144   if (type==0){
2145     Float_t fraction=0;
2146     if (findable>0)
2147       fraction=(Float_t)found/(Float_t)findable;
2148     else
2149       fraction=0;
2150     return fraction;
2151   }
2152   return 0;  // undefined type - default value
2153 }
2154
2155
2156
2157
2158 //_______________________________________________________________________
2159 Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
2160   //
2161   // GetDensity of the clusters on given region between row0 and row1
2162   // Dead zone effect takin into acoount
2163   //
2164   if (!fFriendTrack) return 0.0;
2165   Int_t good  = 0;
2166   Int_t found = 0;
2167   //
2168   Int_t *index=fFriendTrack->GetTPCindices();
2169   for (Int_t i=row0;i<=row1;i++){
2170     Int_t idx = index[i];
2171     if (idx!=-1)  good++;             // track outside of dead zone
2172     if (idx>0)    found++;
2173   }
2174   Float_t density=0.5;
2175   if (good>TMath::Max((row1-row0)*0.5,0.0)) density = Float_t(found)/Float_t(good);
2176   return density;
2177 }
2178
2179 //_______________________________________________________________________
2180 void AliESDtrack::SetTPCpid(const Double_t *p) {
2181   // Sets values for the probability of each particle type (in TPC)
2182   if (!fTPCr) fTPCr = new Double32_t[AliPID::kSPECIES];
2183   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
2184   SetStatus(AliESDtrack::kTPCpid);
2185 }
2186
2187 //_______________________________________________________________________
2188 void AliESDtrack::GetTPCpid(Double_t *p) const {
2189   // Gets the probability of each particle type (in TPC)
2190   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i] = fTPCr ? fTPCr[i] : 0;
2191 }
2192
2193 //_______________________________________________________________________
2194 UChar_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
2195   //---------------------------------------------------------------------
2196   // This function returns indices of the assgined TRD clusters
2197   //---------------------------------------------------------------------
2198   if (idx && fFriendTrack) {
2199     Int_t *index=fFriendTrack->GetTRDindices();
2200
2201     if (index) {
2202       for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
2203     }
2204     else {
2205       for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=-2;
2206     }
2207   }
2208   return fTRDncls;
2209 }
2210
2211 //_______________________________________________________________________
2212 UChar_t AliESDtrack::GetTRDtracklets(Int_t *idx) const {
2213 //
2214 // This function returns the number of TRD tracklets used in tracking
2215 // and it fills the indices of these tracklets in the array "idx" as they
2216 // are registered in the TRD track list.
2217 //
2218 // Caution :
2219 //   1. The idx array has to be allocated with a size >= AliESDtrack::kTRDnPlanes
2220 //   2. The idx array store not only the index but also the layer of the tracklet.
2221 //      Therefore tracks with TRD gaps contain default values for indices [-1]
2222
2223   if (!fFriendTrack) return 0;
2224   if (!idx) return GetTRDntracklets();
2225   Int_t *index=fFriendTrack->GetTRDindices();
2226   Int_t n = 0;
2227   for (Int_t i=0; i<kTRDnPlanes; i++){
2228     if (index){
2229       if(index[i]>=0) n++;
2230       idx[i]=index[i];
2231     }
2232     else idx[i] = -2;
2233   }
2234   return n;
2235 }
2236
2237 //_______________________________________________________________________
2238 void AliESDtrack::SetTRDpid(const Double_t *p) {
2239   // Sets values for the probability of each particle type (in TRD)
2240   if (!fTRDr) fTRDr = new Double32_t[AliPID::kSPECIES];
2241   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
2242   SetStatus(AliESDtrack::kTRDpid);
2243 }
2244
2245 //_______________________________________________________________________
2246 void AliESDtrack::GetTRDpid(Double_t *p) const {
2247   // Gets the probability of each particle type (in TRD)
2248   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i] = fTRDr ? fTRDr[i]:0;
2249 }
2250
2251 //_______________________________________________________________________
2252 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
2253 {
2254   // Sets the probability of particle type iSpecies to p (in TRD)
2255   if (!fTRDr) {
2256     fTRDr = new Double32_t[AliPID::kSPECIES];
2257     for (Int_t i=AliPID::kSPECIES; i--;) fTRDr[i] = 0;
2258   }
2259   fTRDr[iSpecies] = p;
2260 }
2261
2262 Double_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
2263 {
2264   // Returns the probability of particle type iSpecies (in TRD)
2265   return fTRDr ? fTRDr[iSpecies] : 0;
2266 }
2267
2268 //____________________________________________________
2269 Int_t AliESDtrack::GetNumberOfTRDslices() const
2270 {
2271   // built in backward compatibility
2272   Int_t idx = fTRDnSlices - (kTRDnPlanes<<1);
2273   return idx<18 ? fTRDnSlices/kTRDnPlanes : idx/kTRDnPlanes;
2274 }
2275
2276 //____________________________________________________
2277 Double_t AliESDtrack::GetTRDmomentum(Int_t plane, Double_t *sp) const
2278 {
2279 //Returns momentum estimation and optional its error (sp)
2280 // in TRD layer "plane".
2281
2282   if (!fTRDnSlices) {
2283     AliDebug(2, "No TRD info allocated for this track.");
2284     return -1.;
2285   }
2286   if ((plane<0) || (plane>=kTRDnPlanes)) {
2287     AliWarning(Form("Request for TRD plane[%d] outside range.", plane));
2288     return -1.;
2289   }
2290
2291   Int_t idx = fTRDnSlices-(kTRDnPlanes<<1)+plane;
2292   // Protection for backward compatibility
2293   if(idx<(GetNumberOfTRDslices()*kTRDnPlanes)) return -1.;
2294
2295   if(sp) (*sp) = fTRDslices[idx+kTRDnPlanes];
2296   return fTRDslices[idx];
2297 }
2298
2299 //____________________________________________________
2300 Double_t  AliESDtrack::GetTRDslice(Int_t plane, Int_t slice) const {
2301   //Gets the charge from the slice of the plane
2302
2303   if(!fTRDslices) {
2304     //AliError("No TRD slices allocated for this track !");
2305     return -1.;
2306   }
2307   if ((plane<0) || (plane>=kTRDnPlanes)) {
2308     AliError("Info for TRD plane not available !");
2309     return -1.;
2310   }
2311   Int_t ns=GetNumberOfTRDslices();
2312   if ((slice<-1) || (slice>=ns)) {
2313     //AliError("Wrong TRD slice !");
2314     return -1.;
2315   }
2316
2317   if(slice>=0) return fTRDslices[plane*ns + slice];
2318
2319   // return average of the dEdx measurements
2320   Double_t q=0.; Double32_t *s = &fTRDslices[plane*ns];
2321   for (Int_t i=0; i<ns; i++, s++) if((*s)>0.) q+=(*s);
2322   return q/ns;
2323 }
2324
2325 //____________________________________________________
2326 void  AliESDtrack::SetNumberOfTRDslices(Int_t n) {
2327   //Sets the number of slices used for PID
2328   if (fTRDnSlices) return;
2329
2330   fTRDnSlices=n;
2331   fTRDslices=new Double32_t[fTRDnSlices];
2332
2333   // set-up correctly the allocated memory
2334   memset(fTRDslices, 0, n*sizeof(Double32_t));
2335   for (Int_t i=GetNumberOfTRDslices(); i--;) fTRDslices[i]=-1.;
2336 }
2337
2338 //____________________________________________________
2339 void  AliESDtrack::SetTRDslice(Double_t q, Int_t plane, Int_t slice) {
2340   //Sets the charge q in the slice of the plane
2341   if(!fTRDslices) {
2342     AliError("No TRD slices allocated for this track !");
2343     return;
2344   }
2345   if ((plane<0) || (plane>=kTRDnPlanes)) {
2346     AliError("Info for TRD plane not allocated !");
2347     return;
2348   }
2349   Int_t ns=GetNumberOfTRDslices();
2350   if ((slice<0) || (slice>=ns)) {
2351     AliError(Form("Wrong TRD slice %d/%d, NSlices=%d",plane,slice,ns));
2352     return;
2353   }
2354   Int_t n=plane*ns + slice;
2355   fTRDslices[n]=q;
2356 }
2357
2358
2359 //____________________________________________________
2360 void AliESDtrack::SetTRDmomentum(Double_t p, Int_t plane, Double_t *sp)
2361 {
2362   if(!fTRDslices) {
2363     AliError("No TRD slices allocated for this track !");
2364     return;
2365   }
2366   if ((plane<0) || (plane>=kTRDnPlanes)) {
2367     AliError("Info for TRD plane not allocated !");
2368     return;
2369   }
2370
2371   Int_t idx = fTRDnSlices-(kTRDnPlanes<<1)+plane;
2372   // Protection for backward compatibility
2373   if(idx<GetNumberOfTRDslices()*kTRDnPlanes) return;
2374
2375   if(sp) fTRDslices[idx+kTRDnPlanes] = (*sp);
2376   fTRDslices[idx] = p;
2377 }
2378
2379
2380 //_______________________________________________________________________
2381 void AliESDtrack::SetTOFpid(const Double_t *p) {
2382   // Sets the probability of each particle type (in TOF)
2383   if (!fTOFr) fTOFr = new Double32_t[AliPID::kSPECIES];
2384   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
2385   SetStatus(AliESDtrack::kTOFpid);
2386 }
2387
2388 //_______________________________________________________________________
2389 void AliESDtrack::SetTOFLabel(const Int_t *p) {
2390
2391   if(fNtofClusters>0){
2392     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
2393     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
2394     AliESDTOFHit* hit = tofcl->GetTOFHit(0);
2395
2396     if(hit) hit->SetTOFLabel(p);
2397   }
2398   else{
2399     // Sets  (in TOF)
2400     if(!fTOFLabel) fTOFLabel = new Int_t[3];
2401     for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
2402   }
2403 }
2404
2405 //_______________________________________________________________________
2406 void AliESDtrack::GetTOFpid(Double_t *p) const {
2407   // Gets probabilities of each particle type (in TOF)
2408   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i] = fTOFr ? fTOFr[i]:0;
2409 }
2410
2411 //_______________________________________________________________________
2412 void AliESDtrack::GetTOFLabel(Int_t *p) const {
2413   // Gets (in TOF)
2414   if(fNtofClusters>0){
2415     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
2416     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
2417
2418     for (Int_t i=0; i<3; i++) p[i]=tofcl->GetLabel(i);
2419   }
2420   else{
2421     if(fTOFLabel) for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
2422     else for (int i=3;i--;) p[i] = -1;
2423   }
2424 }
2425
2426 //_______________________________________________________________________
2427 void AliESDtrack::GetTOFInfo(Float_t *info) const {
2428   // Gets (in TOF)
2429   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
2430 }
2431
2432 //_______________________________________________________________________
2433 void AliESDtrack::SetTOFInfo(Float_t*info) {
2434   // Gets (in TOF)
2435   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
2436 }
2437
2438
2439
2440 //_______________________________________________________________________
2441 void AliESDtrack::SetHMPIDpid(const Double_t *p) {
2442   // Sets the probability of each particle type (in HMPID)
2443   if (!fHMPIDr) fHMPIDr = new Double32_t[AliPID::kSPECIES];
2444   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
2445   SetStatus(AliESDtrack::kHMPIDpid);
2446 }
2447
2448 //_______________________________________________________________________
2449 void  AliESDtrack::SetTPCdEdxInfo(AliTPCdEdxInfo * dEdxInfo){
2450   if(fTPCdEdxInfo) delete fTPCdEdxInfo;
2451   fTPCdEdxInfo = dEdxInfo;
2452 }
2453
2454 //_______________________________________________________________________
2455 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
2456   // Gets probabilities of each particle type (in HMPID)
2457   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i] = fHMPIDr ? fHMPIDr[i]:0;
2458 }
2459
2460
2461
2462 //_______________________________________________________________________
2463 void AliESDtrack::SetESDpid(const Double_t *p) {
2464   // Sets the probability of each particle type for the ESD track
2465   if (!fR) fR = new Double32_t[AliPID::kSPECIES];
2466   SetPIDValues(fR,p,AliPID::kSPECIES);
2467   SetStatus(AliESDtrack::kESDpid);
2468 }
2469
2470 //_______________________________________________________________________
2471 void AliESDtrack::GetESDpid(Double_t *p) const {
2472   // Gets probability of each particle type for the ESD track
2473   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i] = fR ? fR[i]:0;
2474 }
2475
2476 //_______________________________________________________________________
2477 Bool_t AliESDtrack::RelateToVertexTPC(const AliESDVertex *vtx,
2478 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
2479   //
2480   // Try to relate the TPC-only track parameters to the vertex "vtx",
2481   // if the (rough) transverse impact parameter is not bigger then "maxd".
2482   //            Magnetic field is "b" (kG).
2483   //
2484   // a) The TPC-only paramters are extapolated to the DCA to the vertex.
2485   // b) The impact parameters and their covariance matrix are calculated.
2486   // c) An attempt to constrain the TPC-only params to the vertex is done.
2487   //    The constrained params are returned via "cParam".
2488   //
2489   // In the case of success, the returned value is kTRUE
2490   // otherwise, it's kFALSE)
2491   //
2492
2493   if (!fTPCInner) return kFALSE;
2494   if (!vtx) return kFALSE;
2495
2496   Double_t dz[2],cov[3];
2497   if (!fTPCInner->PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
2498
2499   fdTPC = dz[0];
2500   fzTPC = dz[1];
2501   fCddTPC = cov[0];
2502   fCdzTPC = cov[1];
2503   fCzzTPC = cov[2];
2504
2505   Double_t covar[6]; vtx->GetCovMatrix(covar);
2506   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2507   Double_t c[3]={covar[2],0.,covar[5]};
2508
2509   Double_t chi2=GetPredictedChi2(p,c);
2510   if (chi2>kVeryBig) return kFALSE;
2511
2512   fCchi2TPC=chi2;
2513
2514   if (!cParam) return kTRUE;
2515
2516   *cParam = *fTPCInner;
2517   if (!cParam->Update(p,c)) return kFALSE;
2518
2519   return kTRUE;
2520 }
2521
2522 //_______________________________________________________________________
2523 Bool_t AliESDtrack::RelateToVertexTPCBxByBz(const AliESDVertex *vtx,
2524 Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam) {
2525   //
2526   // Try to relate the TPC-only track parameters to the vertex "vtx",
2527   // if the (rough) transverse impact parameter is not bigger then "maxd".
2528   //
2529   // All three components of the magnetic field ,"b[3]" (kG),
2530   // are taken into account.
2531   //
2532   // a) The TPC-only paramters are extapolated to the DCA to the vertex.
2533   // b) The impact parameters and their covariance matrix are calculated.
2534   // c) An attempt to constrain the TPC-only params to the vertex is done.
2535   //    The constrained params are returned via "cParam".
2536   //
2537   // In the case of success, the returned value is kTRUE
2538   // otherwise, it's kFALSE)
2539   //
2540
2541   if (!fTPCInner) return kFALSE;
2542   if (!vtx) return kFALSE;
2543
2544   Double_t dz[2],cov[3];
2545   if (!fTPCInner->PropagateToDCABxByBz(vtx, b, maxd, dz, cov)) return kFALSE;
2546
2547   fdTPC = dz[0];
2548   fzTPC = dz[1];
2549   fCddTPC = cov[0];
2550   fCdzTPC = cov[1];
2551   fCzzTPC = cov[2];
2552
2553   Double_t covar[6]; vtx->GetCovMatrix(covar);
2554   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2555   Double_t c[3]={covar[2],0.,covar[5]};
2556
2557   Double_t chi2=GetPredictedChi2(p,c);
2558   if (chi2>kVeryBig) return kFALSE;
2559
2560   fCchi2TPC=chi2;
2561
2562   if (!cParam) return kTRUE;
2563
2564   *cParam = *fTPCInner;
2565   if (!cParam->Update(p,c)) return kFALSE;
2566
2567   return kTRUE;
2568 }
2569
2570 //_______________________________________________________________________
2571 Bool_t AliESDtrack::RelateToVertex(const AliESDVertex *vtx,
2572 Double_t b, Double_t maxd, AliExternalTrackParam *cParam) {
2573   //
2574   // Try to relate this track to the vertex "vtx",
2575   // if the (rough) transverse impact parameter is not bigger then "maxd".
2576   //            Magnetic field is "b" (kG).
2577   //
2578   // a) The track gets extapolated to the DCA to the vertex.
2579   // b) The impact parameters and their covariance matrix are calculated.
2580   // c) An attempt to constrain this track to the vertex is done.
2581   //    The constrained params are returned via "cParam".
2582   //
2583   // In the case of success, the returned value is kTRUE
2584   // (otherwise, it's kFALSE)
2585   //
2586
2587   if (!vtx) return kFALSE;
2588
2589   Double_t dz[2],cov[3];
2590   if (!PropagateToDCA(vtx, b, maxd, dz, cov)) return kFALSE;
2591
2592   fD = dz[0];
2593   fZ = dz[1];
2594   fCdd = cov[0];
2595   fCdz = cov[1];
2596   fCzz = cov[2];
2597
2598   Double_t covar[6]; vtx->GetCovMatrix(covar);
2599   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2600   Double_t c[3]={covar[2],0.,covar[5]};
2601
2602   Double_t chi2=GetPredictedChi2(p,c);
2603   if (chi2>kVeryBig) return kFALSE;
2604
2605   fCchi2=chi2;
2606
2607
2608   //--- Could now these lines be removed ? ---
2609   delete fCp;
2610   fCp=new AliExternalTrackParam(*this);
2611
2612   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
2613   //----------------------------------------
2614
2615   fVertexID = vtx->GetID();
2616
2617   if (!cParam) return kTRUE;
2618
2619   *cParam = *this;
2620   if (!cParam->Update(p,c)) return kFALSE;
2621
2622   return kTRUE;
2623 }
2624
2625 //_______________________________________________________________________
2626 Bool_t AliESDtrack::RelateToVertexBxByBz(const AliESDVertex *vtx,
2627 Double_t b[3], Double_t maxd, AliExternalTrackParam *cParam) {
2628   //
2629   // Try to relate this track to the vertex "vtx",
2630   // if the (rough) transverse impact parameter is not bigger then "maxd".
2631   //            Magnetic field is "b" (kG).
2632   //
2633   // a) The track gets extapolated to the DCA to the vertex.
2634   // b) The impact parameters and their covariance matrix are calculated.
2635   // c) An attempt to constrain this track to the vertex is done.
2636   //    The constrained params are returned via "cParam".
2637   //
2638   // In the case of success, the returned value is kTRUE
2639   // (otherwise, it's kFALSE)
2640   //
2641
2642   if (!vtx) return kFALSE;
2643
2644   Double_t dz[2],cov[3];
2645   if (!PropagateToDCABxByBz(vtx, b, maxd, dz, cov)) return kFALSE;
2646
2647   fD = dz[0];
2648   fZ = dz[1];
2649   fCdd = cov[0];
2650   fCdz = cov[1];
2651   fCzz = cov[2];
2652
2653   Double_t covar[6]; vtx->GetCovMatrix(covar);
2654   Double_t p[2]={GetParameter()[0]-dz[0],GetParameter()[1]-dz[1]};
2655   Double_t c[3]={covar[2],0.,covar[5]};
2656
2657   Double_t chi2=GetPredictedChi2(p,c);
2658   if (chi2>kVeryBig) return kFALSE;
2659
2660   fCchi2=chi2;
2661
2662
2663   //--- Could now these lines be removed ? ---
2664   delete fCp;
2665   fCp=new AliExternalTrackParam(*this);
2666
2667   if (!fCp->Update(p,c)) {delete fCp; fCp=0; return kFALSE;}
2668   //----------------------------------------
2669
2670   fVertexID = vtx->GetID();
2671
2672   if (!cParam) return kTRUE;
2673
2674   *cParam = *this;
2675   if (!cParam->Update(p,c)) return kFALSE;
2676
2677   return kTRUE;
2678 }
2679
2680 //_______________________________________________________________________
2681 void AliESDtrack::Print(Option_t *) const {
2682   // Prints info on the track
2683   AliExternalTrackParam::Print();
2684   printf("ESD track info\n") ;
2685   Double_t p[AliPID::kSPECIES] ;
2686   Int_t index = 0 ;
2687   if( IsOn(kITSpid) ){
2688     printf("From ITS: ") ;
2689     GetITSpid(p) ;
2690     for(index = 0 ; index < AliPID::kSPECIES; index++)
2691       printf("%f, ", p[index]) ;
2692     printf("\n           signal = %f\n", GetITSsignal()) ;
2693   }
2694   if( IsOn(kTPCpid) ){
2695     printf("From TPC: ") ;
2696     GetTPCpid(p) ;
2697     for(index = 0 ; index < AliPID::kSPECIES; index++)
2698       printf("%f, ", p[index]) ;
2699     printf("\n           signal = %f\n", GetTPCsignal()) ;
2700   }
2701   if( IsOn(kTRDpid) ){
2702     printf("From TRD: ") ;
2703     GetTRDpid(p) ;
2704     for(index = 0 ; index < AliPID::kSPECIES; index++)
2705       printf("%f, ", p[index]) ;
2706       printf("\n           signal = %f\n", GetTRDsignal()) ;
2707       printf("\n           NchamberdEdx = %d\n", GetTRDNchamberdEdx()) ;
2708       printf("\n           NclusterdEdx = %d\n", GetTRDNclusterdEdx()) ;
2709   }
2710   if( IsOn(kTOFpid) ){
2711     printf("From TOF: ") ;
2712     GetTOFpid(p) ;
2713     for(index = 0 ; index < AliPID::kSPECIES; index++)
2714       printf("%f, ", p[index]) ;
2715     printf("\n           signal = %f\n", GetTOFsignal()) ;
2716   }
2717   if( IsOn(kHMPIDpid) ){
2718     printf("From HMPID: ") ;
2719     GetHMPIDpid(p) ;
2720     for(index = 0 ; index < AliPID::kSPECIES; index++)
2721       printf("%f, ", p[index]) ;
2722     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
2723   }
2724 }
2725
2726
2727 //
2728 // Draw functionality
2729 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
2730 //
2731 void AliESDtrack::FillPolymarker(TPolyMarker3D *pol, Float_t magF, Float_t minR, Float_t maxR, Float_t stepR){
2732   //
2733   // Fill points in the polymarker
2734   //
2735   TObjArray arrayRef;
2736   arrayRef.AddLast(new AliExternalTrackParam(*this));
2737   if (fIp) arrayRef.AddLast(new AliExternalTrackParam(*fIp));
2738   if (fOp) arrayRef.AddLast(new AliExternalTrackParam(*fOp));
2739   if (fHMPIDp) arrayRef.AddLast(new AliExternalTrackParam(*fHMPIDp));
2740   //
2741   Double_t mpos[3]={0,0,0};
2742   Int_t entries=arrayRef.GetEntries();
2743   for (Int_t i=0;i<entries;i++){
2744     Double_t pos[3];
2745     ((AliExternalTrackParam*)arrayRef.At(i))->GetXYZ(pos);
2746     mpos[0]+=pos[0]/entries;
2747     mpos[1]+=pos[1]/entries;
2748     mpos[2]+=pos[2]/entries;
2749   }
2750   // Rotate to the mean position
2751   //
2752   Float_t fi= TMath::ATan2(mpos[1],mpos[0]);
2753   for (Int_t i=0;i<entries;i++){
2754     Bool_t res = ((AliExternalTrackParam*)arrayRef.At(i))->Rotate(fi);
2755     if (!res) delete arrayRef.RemoveAt(i);
2756   }
2757   Int_t counter=0;
2758   for (Double_t r=minR; r<maxR; r+=stepR){
2759     Double_t sweight=0;
2760     Double_t mlpos[3]={0,0,0};
2761     for (Int_t i=0;i<entries;i++){
2762       Double_t point[3]={0,0,0};
2763       AliExternalTrackParam *param = ((AliExternalTrackParam*)arrayRef.At(i));
2764       if (!param) continue;
2765       if (param->GetXYZAt(r,magF,point)){
2766         Double_t weight = 1./(10.+(r-param->GetX())*(r-param->GetX()));
2767         sweight+=weight;
2768         mlpos[0]+=point[0]*weight;
2769         mlpos[1]+=point[1]*weight;
2770         mlpos[2]+=point[2]*weight;
2771       }
2772     }
2773     if (sweight>0){
2774       mlpos[0]/=sweight;
2775       mlpos[1]/=sweight;
2776       mlpos[2]/=sweight;
2777       pol->SetPoint(counter,mlpos[0],mlpos[1], mlpos[2]);
2778       //      printf("xyz\t%f\t%f\t%f\n",mlpos[0], mlpos[1],mlpos[2]);
2779       counter++;
2780     }
2781   }
2782 }
2783
2784 //_______________________________________________________________________
2785 void AliESDtrack::SetITSdEdxSamples(const Double_t s[4]) {
2786   //
2787   // Store the dE/dx samples measured by the two SSD and two SDD layers.
2788   // These samples are corrected for the track segment length.
2789   //
2790   for (Int_t i=0; i<4; i++) fITSdEdxSamples[i]=s[i];
2791 }
2792
2793 //_______________________________________________________________________
2794 void AliESDtrack::GetITSdEdxSamples(Double_t s[4]) const {
2795   //
2796   // Get the dE/dx samples measured by the two SSD and two SDD layers.
2797   // These samples are corrected for the track segment length.
2798   //
2799   for (Int_t i=0; i<4; i++) s[i]=fITSdEdxSamples[i];
2800 }
2801
2802
2803 UShort_t   AliESDtrack::GetTPCnclsS(Int_t i0,Int_t i1) const{
2804   //
2805   // get number of shared TPC clusters
2806   //
2807   return  fTPCSharedMap.CountBits(i0)-fTPCSharedMap.CountBits(i1);
2808 }
2809
2810 UShort_t   AliESDtrack::GetTPCncls(Int_t i0,Int_t i1) const{
2811   //
2812   // get number of TPC clusters
2813   //
2814   return  fTPCClusterMap.CountBits(i0)-fTPCClusterMap.CountBits(i1);
2815 }
2816
2817 //____________________________________________________________________
2818 Double_t AliESDtrack::GetChi2TPCConstrainedVsGlobal(const AliESDVertex* vtx) const
2819 {
2820   // Calculates the chi2 between the TPC track (TPCinner) constrained to the primary vertex and the global track
2821   //
2822   // Returns -1 in case the calculation failed
2823   //
2824   // Value is cached as a non-persistent member.
2825   //
2826   // Code adapted from original code by GSI group (Jacek, Marian, Michael)
2827
2828   // cache, ignoring that a different vertex might be passed
2829   if (fCacheChi2TPCConstrainedVsGlobalVertex == vtx)
2830     return fCacheChi2TPCConstrainedVsGlobal;
2831
2832   fCacheChi2TPCConstrainedVsGlobal = -1;
2833   fCacheChi2TPCConstrainedVsGlobalVertex = vtx;
2834
2835   Double_t x[3];
2836   GetXYZ(x);
2837   Double_t b[3];
2838   AliTrackerBase::GetBxByBz(x,b);
2839
2840   if (!fTPCInner)  {
2841     AliWarning("Could not get TPC Inner Param.");
2842     return fCacheChi2TPCConstrainedVsGlobal;
2843   }
2844
2845   // clone for constraining
2846   AliExternalTrackParam* tpcInnerC = new AliExternalTrackParam(*fTPCInner);
2847   if (!tpcInnerC) {
2848     AliWarning("Clone of TPCInnerParam failed.");
2849     return fCacheChi2TPCConstrainedVsGlobal;
2850   }
2851
2852   // transform to the track reference frame
2853   Bool_t isOK = tpcInnerC->Rotate(GetAlpha());
2854   isOK &= tpcInnerC->PropagateTo(GetX(), b[2]);
2855   if (!isOK) {
2856     delete tpcInnerC;
2857     tpcInnerC = 0;
2858     AliWarning("Rotation/Propagation of track failed.") ;
2859     return fCacheChi2TPCConstrainedVsGlobal;
2860   }
2861
2862   // constrain TPCinner
2863   isOK = tpcInnerC->ConstrainToVertex(vtx, b);
2864
2865   // transform to the track reference frame
2866   isOK &= tpcInnerC->Rotate(GetAlpha());
2867   isOK &= tpcInnerC->PropagateTo(GetX(), b[2]);
2868
2869   if (!isOK) {
2870     AliWarning("ConstrainTPCInner failed.") ;
2871     delete tpcInnerC;
2872     tpcInnerC = 0;
2873     return fCacheChi2TPCConstrainedVsGlobal;
2874   }
2875
2876   // calculate chi2 between vi and vj vectors
2877   // with covi and covj covariance matrices
2878   // chi2ij = (vi-vj)^(T)*(covi+covj)^(-1)*(vi-vj)
2879   TMatrixD deltaT(5,1);
2880   TMatrixD delta(1,5);
2881   TMatrixD covarM(5,5);
2882
2883   for (Int_t ipar=0; ipar<5; ipar++) {
2884     deltaT(ipar,0) = tpcInnerC->GetParameter()[ipar] - GetParameter()[ipar];
2885     delta(0,ipar) = tpcInnerC->GetParameter()[ipar] - GetParameter()[ipar];
2886
2887     for (Int_t jpar=0; jpar<5; jpar++) {
2888       Int_t index = GetIndex(ipar,jpar);
2889       covarM(ipar,jpar) = GetCovariance()[index]+tpcInnerC->GetCovariance()[index];
2890     }
2891   }
2892   // chi2 distance TPC constrained and TPC+ITS
2893   TMatrixD covarMInv = covarM.Invert();
2894   TMatrixD mat2 = covarMInv*deltaT;
2895   TMatrixD chi2 = delta*mat2;
2896
2897   delete tpcInnerC;
2898   tpcInnerC = 0;
2899
2900   fCacheChi2TPCConstrainedVsGlobal = chi2(0,0);
2901   return fCacheChi2TPCConstrainedVsGlobal;
2902 }
2903
2904 void AliESDtrack::SetDetectorPID(const AliDetectorPID *pid)
2905 {
2906   //
2907   // Set the detector PID
2908   //
2909   if (fDetectorPID) delete fDetectorPID;
2910   fDetectorPID=pid;
2911
2912 }
2913
2914 Double_t AliESDtrack::GetLengthInActiveZone( Int_t mode, Double_t deltaY, Double_t deltaZ, Double_t bz, Double_t exbPhi , TTreeSRedirector * pcstream) const {
2915   //
2916   // Input parameters:
2917   //   mode  - type of external track parameters
2918   //   deltaY - user defined "dead region" in cm
2919   //   deltaZ - user defined "active region" in cm (250 cm drift lenght - 14 cm L1 delay
2920   //   bz     - magnetic field
2921   //   exbPhi - optional rotation due to the ExB effect
2922   // return value:
2923   //   the length of the track in cm in "active volume" of the TPC
2924   //
2925   if (mode==0) return GetLengthInActiveZone(this, deltaY,deltaZ,bz, exbPhi,pcstream);
2926   if (mode==1) return GetLengthInActiveZone(fIp, deltaY,deltaZ,bz, exbPhi,pcstream);
2927   if (mode==2) return GetLengthInActiveZone(fOp, deltaY,deltaZ,bz, exbPhi,pcstream);
2928   return 0;
2929 }
2930
2931 Double_t AliESDtrack::GetLengthInActiveZone(const AliExternalTrackParam  *paramT, Double_t deltaY, Double_t deltaZ, Double_t bz, Double_t exbPhi , TTreeSRedirector * pcstream) {
2932   //
2933   // Numerical code to calculate the length of the track in active region of the TPC
2934   // ( can be speed up if somebody wants to invest time - analysical version shoult be possible)
2935   //
2936   // Input parameters:
2937   //   paramT - external track parameters
2938   //   deltaY - user defined "dead region" in cm
2939   //   deltaZ - user defined "active region" in cm (250 cm drift lenght - 14 cm L1 delay
2940   //   bz     - magnetic field
2941   //   exbPhi - optional rotation due to the ExB effect
2942   // return value:
2943   //   the length of the track in cm in "active volume" of the TPC
2944   //
2945   const Double_t rIn=85;
2946   const Double_t rOut=245;
2947   Double_t xyz[3], pxyz[3];
2948   if (paramT->GetXYZAt(rIn,bz,xyz)){
2949     paramT->GetPxPyPzAt(rIn,bz,pxyz);
2950   }else{
2951     paramT->GetXYZ(xyz);
2952     paramT->GetPxPyPz(pxyz);
2953   }
2954   //
2955   Double_t dca   = -paramT->GetD(0,0,bz);  // get impact parameter distance to point (0,0)
2956   Double_t radius= TMath::Abs(1/paramT->GetC(bz));  //
2957   Double_t sign  = paramT->GetSign();
2958   Double_t R0    = TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]);   // radius at current point
2959   Double_t phiR0 = TMath::ATan2(xyz[1],xyz[0]);                // angle of given point
2960   Double_t dPhiR0= -TMath::ASin((dca*dca-2*dca*radius*sign+R0*R0)/(2*R0*(dca-radius*sign)));
2961   Double_t phi0  = phiR0-(dPhiR0);  // global phi offset to be added
2962   //
2963   //
2964   AliExternalTrackParam paramR=(*paramT);
2965   Double_t length=0;
2966   for (Double_t R=rIn; R<=rOut; R++){
2967     Double_t sinPhi=(dca*dca-2*dca*radius*sign+R*R)/(2*R*(dca-radius*sign));
2968     if (TMath::Abs(sinPhi)>=1) continue;
2969     Double_t dphi     = -TMath::ASin(sinPhi);
2970     Double_t phi      = phi0+dphi;                           // global phi
2971     Int_t    sector   = TMath::Nint(9*phi/(TMath::Pi()));
2972     Double_t dPhiEdge = phi-(sector*TMath::Pi()/9)+exbPhi;   // distance to sector boundary in rphi
2973     Double_t dX   = R*TMath::Cos(phi)-xyz[0];
2974     Double_t dY   = R*TMath::Sin(phi)-xyz[1];
2975     Double_t deltaPhi = 2*TMath::ASin(0.5*TMath::Sqrt(dX*dX+dY*dY)/radius);
2976     Double_t z = xyz[2]+deltaPhi*radius*paramT->GetTgl();
2977     if (TMath::Abs(dPhiEdge*R)>deltaY && TMath::Abs(z)<deltaZ){
2978       length++;
2979     }
2980     //    Double_t deltaZ= dphi*radius;
2981     if (pcstream){
2982       //should we keep debug possibility ?
2983       AliExternalTrackParam paramTcopy=(*paramT);
2984       paramR.Rotate(phi);
2985       paramR.PropagateTo(R,bz);
2986       (*pcstream)<<"debugEdge"<<
2987         "R="<<R<<                   // radius
2988         "dphiEdge="<<dPhiEdge<<     // distance to edge
2989         "phi0="<<phi0<<             // phi0 -phi at the track initial position
2990         "phi="<<phi<<               //
2991         "z="<<z<<
2992         "pT.="<<&paramTcopy<<
2993         "pR.="<<&paramR<<
2994         "\n";
2995     }
2996   }
2997   return length;
2998 }
2999
3000 Double_t AliESDtrack::GetMassForTracking() const
3001 {
3002   int pid = fPIDForTracking;
3003   if (pid<AliPID::kPion) pid = AliPID::kPion;
3004   double m = AliPID::ParticleMass(pid);
3005   return (fPIDForTracking==AliPID::kHe3 || fPIDForTracking==AliPID::kAlpha) ? -m : m;
3006 }
3007
3008
3009 void    AliESDtrack::SetTOFclusterArray(Int_t /*ncluster*/,Int_t */*TOFcluster*/){
3010   AliInfo("Method has to be implemented!");
3011 //   fNtofClusters=ncluster;
3012 //   if(TOFcluster == fTOFcluster) return;
3013 //   if(fTOFcluster){ // reset previous content
3014 //     delete[] fTOFcluster;
3015 //     fTOFcluster = NULL;
3016 //     fNtofClusters=0;
3017 //   }
3018
3019 //   if(ncluster){ // set new content
3020 //     fTOFcluster = new Int_t[fNtofClusters];
3021 //     for(Int_t i=0;i < fNtofClusters;i++) fTOFcluster[i] = TOFcluster[i];
3022 //   }
3023 //   else
3024 //     fTOFcluster = 0;
3025 }
3026
3027 //____________________________________________
3028 void  AliESDtrack::SuppressTOFMatches()
3029 {
3030   // remove reference to this track from TOF clusters
3031   if (!fNtofClusters || !GetESDEvent()) return;
3032   TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3033   for (;fNtofClusters--;) {
3034     AliESDTOFCluster* clTOF = (AliESDTOFCluster*)tofclArray->At(fTOFcluster[fNtofClusters]);
3035     clTOF->SuppressMatchedTrack(GetID());
3036     if (!clTOF->GetNMatchableTracks()) { // remove this cluster
3037       int last = tofclArray->GetEntriesFast()-1;
3038       AliESDTOFCluster* clTOFL = (AliESDTOFCluster*)tofclArray->At(last);
3039       if (last != fTOFcluster[fNtofClusters]) {
3040         *clTOF = *clTOFL; // move last cluster to the place of eliminated one
3041         // fix the references on this cluster
3042         clTOF->FixSelfReferences(last,fTOFcluster[fNtofClusters]);
3043       }
3044       tofclArray->RemoveAt(last);
3045     }
3046   }
3047 }
3048
3049 //____________________________________________
3050 void  AliESDtrack::ReplaceTOFTrackID(int oldID, int newID)
3051 {
3052   // replace the ID in TOF clusters references to this track
3053   if (!fNtofClusters || !GetESDEvent()) return;
3054   TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3055   for (int it=fNtofClusters;it--;) {
3056     AliESDTOFCluster* clTOF = (AliESDTOFCluster*)tofclArray->At(fTOFcluster[it]);
3057     clTOF->ReplaceMatchedTrackID(oldID,newID);
3058   }
3059 }
3060
3061 //____________________________________________
3062 void  AliESDtrack::ReplaceTOFClusterID(int oldID, int newID)
3063 {
3064   // replace the referenc on TOF cluster oldID by newID
3065   if (!fNtofClusters || !GetESDEvent()) return;
3066   for (int it=fNtofClusters;it--;) {
3067     if (fTOFcluster[it] == oldID) {
3068       fTOFcluster[it] = newID;
3069       return;
3070     }
3071   }
3072 }
3073
3074 //____________________________________________
3075 void  AliESDtrack::ReplaceTOFMatchID(int oldID, int newID)
3076 {
3077   // replace in the ESDTOFCluster associated with this track the id of the corresponding
3078   // ESDTOFMatch from oldID to newID
3079   if (!fNtofClusters || !GetESDEvent()) return;
3080   TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3081   for (int it=fNtofClusters;it--;) {
3082     AliESDTOFCluster* clTOF = (AliESDTOFCluster*)tofclArray->At(fTOFcluster[it]);
3083     clTOF->ReplaceMatchID(oldID,newID);
3084   }
3085 }
3086
3087 //____________________________________________
3088 void AliESDtrack::AddTOFcluster(Int_t icl)
3089 {
3090   fNtofClusters++;
3091
3092   Int_t *old = fTOFcluster;
3093   fTOFcluster = new Int_t[fNtofClusters];
3094
3095   for(Int_t i=0;i < fNtofClusters-1;i++) fTOFcluster[i] = old[i];
3096   fTOFcluster[fNtofClusters-1] = icl;
3097
3098   if(fNtofClusters-1)  delete[] old; // delete previous content
3099
3100 }
3101
3102 //____________________________________________
3103 void AliESDtrack::SetTOFsignal(Double_t tof)
3104 {
3105   if(fNtofClusters>0 && GetESDEvent()){
3106     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3107     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3108     AliESDTOFHit* hit = tofcl->GetTOFHit(0);
3109     if(hit) hit->SetTime(tof);
3110   }
3111   else{
3112     if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3113     fTOFsignal=tof;
3114   }
3115 }
3116 //____________________________________________
3117 void AliESDtrack::SetTOFCalChannel(Int_t index){
3118   if(fNtofClusters>0 && GetESDEvent()){
3119     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3120     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3121     AliESDTOFHit* hit = tofcl->GetTOFHit(0);
3122     if(hit) hit->SetTOFchannel(index);
3123   }
3124   else{
3125     if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3126     fTOFCalChannel=index;
3127   }
3128 }
3129 //____________________________________________
3130 void AliESDtrack::SetTOFsignalToT(Double_t ToT){
3131   if(fNtofClusters>0 && GetESDEvent()){
3132     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3133     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3134     AliESDTOFHit* hit = tofcl->GetTOFHit(0);
3135     if(hit) hit->SetTOT(ToT);
3136   }
3137   else{
3138     if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3139     fTOFsignalToT=ToT;
3140   }
3141 }
3142 //____________________________________________
3143 void AliESDtrack::SetTOFsignalRaw(Double_t tof){
3144   if(fNtofClusters>0 && GetESDEvent()){
3145     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3146     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3147     AliESDTOFHit* hit = tofcl->GetTOFHit(0);
3148     if(hit) hit->SetTimeRaw(tof);
3149   }
3150   else{
3151     if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3152     fTOFsignalRaw=tof;
3153   }
3154 }
3155 //____________________________________________
3156 void AliESDtrack::SetTOFsignalDz(Double_t dz){
3157   Int_t index = -1;
3158   AliESDTOFCluster *tofcl;
3159
3160   if(fNtofClusters>0 && GetESDEvent()){
3161     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3162     tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3163
3164     for(Int_t i=0;i < tofcl->GetNMatchableTracks();i++){
3165       if(tofcl->GetTrackIndex(i) == GetID()) index = i;
3166     }
3167
3168   }
3169   if(index > -1){
3170     AliESDTOFMatch* match = tofcl->GetTOFMatch(index);
3171     if(match){
3172       match->SetDz(dz);
3173     }
3174   }
3175   else{
3176     if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3177     fTOFsignalDz=dz;
3178   }
3179
3180
3181 }
3182 //____________________________________________
3183 void AliESDtrack::SetTOFsignalDx(Double_t dx){
3184   Int_t index = -1;
3185   AliESDTOFCluster *tofcl;
3186
3187   if(fNtofClusters>0 && GetESDEvent()){
3188     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3189     tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3190
3191     for(Int_t i=0;i < tofcl->GetNMatchableTracks();i++){
3192       if(tofcl->GetTrackIndex(i) == GetID()) index = i;
3193     }
3194
3195   }
3196   if(index > -1){
3197     AliESDTOFMatch* match = tofcl->GetTOFMatch(index);
3198     if(match){
3199       match->SetDx(dx);
3200     }
3201   }
3202   else{
3203     if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3204     fTOFsignalDx=dx;
3205   }
3206 }
3207 //____________________________________________
3208 void AliESDtrack::SetTOFDeltaBC(Short_t deltaBC){
3209   if(fNtofClusters>0 && GetESDEvent()){
3210     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3211     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3212     AliESDTOFHit* hit = tofcl->GetTOFHit(0);
3213     if(hit) hit->SetDeltaBC(deltaBC);
3214   }
3215   else{
3216     if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3217     fTOFdeltaBC=deltaBC;
3218   }
3219 }
3220 //____________________________________________
3221 void AliESDtrack::SetTOFL0L1(Short_t l0l1){
3222   if(fNtofClusters>0 && GetESDEvent()){
3223     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3224     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3225     AliESDTOFHit* hit = tofcl->GetTOFHit(0);
3226     if(hit) hit->SetL0L1Latency(l0l1);
3227   }
3228   else{
3229     if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3230     fTOFl0l1=l0l1;
3231   }
3232 }
3233 //____________________________________________
3234 Double_t AliESDtrack::GetTOFsignal() const
3235 {
3236   if(fNtofClusters>0 && GetESDEvent()){
3237     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3238     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3239
3240     return tofcl->GetTime();
3241   }
3242   else if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3243
3244   return fTOFsignal;
3245 }
3246
3247 //____________________________________________
3248 Double_t AliESDtrack::GetTOFsignalToT() const
3249 {
3250   if(fNtofClusters>0 && GetESDEvent()){
3251     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3252     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3253
3254     return tofcl->GetTOT();
3255   }
3256   else if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3257
3258   return fTOFsignalToT;
3259 }
3260
3261 //____________________________________________
3262 Double_t AliESDtrack::GetTOFsignalRaw() const
3263 {
3264   if(fNtofClusters>0 && GetESDEvent()){
3265     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3266     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3267
3268     return tofcl->GetTimeRaw();
3269   }
3270   else if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3271
3272   return fTOFsignalRaw;
3273 }
3274
3275 //____________________________________________
3276 Double_t AliESDtrack::GetTOFsignalDz() const
3277 {
3278
3279   AliESDTOFCluster *tofcl;
3280
3281   Int_t index = -1;
3282   if(fNtofClusters>0 && GetESDEvent()){
3283     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3284     tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3285
3286     for(Int_t i=0;i < tofcl->GetNMatchableTracks();i++){
3287       if(tofcl->GetTrackIndex(i) == GetID()) index = i;
3288     }
3289   }
3290   else if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3291
3292   if(fNtofClusters>0 && index > -1){
3293     return tofcl->GetDz(index);
3294   }
3295   return fTOFsignalDz;
3296 }
3297
3298 //____________________________________________
3299 Double_t AliESDtrack::GetTOFsignalDx() const
3300 {
3301   AliESDTOFCluster *tofcl;
3302
3303   Int_t index = -1;
3304   if(fNtofClusters>0 && GetESDEvent()){
3305     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3306     tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3307     for(Int_t i=0;i < tofcl->GetNMatchableTracks();i++){
3308       if(tofcl->GetTrackIndex(i) == GetID()) index = i;
3309     }
3310   }
3311   else if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3312   if(fNtofClusters>0 && index > -1){
3313     return tofcl->GetDx(index);
3314   }
3315   return fTOFsignalDx;
3316 }
3317
3318 //____________________________________________
3319 Short_t  AliESDtrack::GetTOFDeltaBC() const
3320 {
3321   if(fNtofClusters>0 && GetESDEvent()){
3322     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3323     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3324     return tofcl->GetDeltaBC();
3325   }
3326   else if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3327
3328   return fTOFdeltaBC;
3329 }
3330
3331 //____________________________________________
3332 Short_t  AliESDtrack::GetTOFL0L1() const
3333 {
3334   if(fNtofClusters>0 && GetESDEvent()){
3335     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3336     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3337
3338     return tofcl->GetL0L1Latency();
3339   }
3340   else if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3341
3342   return fTOFl0l1;
3343 }
3344
3345 //____________________________________________
3346 Int_t   AliESDtrack::GetTOFCalChannel() const
3347 {
3348   if(fNtofClusters>0 && GetESDEvent()){
3349     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3350     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3351
3352     return tofcl->GetTOFchannel();
3353   }
3354   else if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3355
3356   return fTOFCalChannel;
3357 }
3358
3359 //____________________________________________
3360 Int_t   AliESDtrack::GetTOFcluster() const
3361 {
3362   if(fNtofClusters>0 && GetESDEvent()){
3363     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3364     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3365
3366     return tofcl->GetClusterIndex();
3367   }
3368   else if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3369
3370   return fTOFindex;
3371 }
3372
3373 //____________________________________________
3374 Int_t   AliESDtrack::GetTOFclusterN() const
3375 {
3376   return fNtofClusters;
3377 }
3378
3379 //____________________________________________
3380 Bool_t  AliESDtrack::IsTOFHitAlreadyMatched() const{
3381   if(fNtofClusters>0 && GetESDEvent()){
3382     TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3383     AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[0]);
3384
3385     if (tofcl->GetNMatchableTracks() > 1)
3386       return kTRUE;
3387   }
3388   else if(fNtofClusters>0) AliInfo("No AliESDEvent available here!\n");
3389
3390   return kFALSE;
3391 }
3392
3393 //____________________________________________
3394 void AliESDtrack::ReMapTOFcluster(Int_t ncl,Int_t *mapping){
3395   for(Int_t i=0;i<fNtofClusters;i++){
3396     if(fTOFcluster[i]<ncl && fTOFcluster[i]>-1)
3397       fTOFcluster[i] = mapping[fTOFcluster[i]];
3398     else
3399       AliInfo(Form("TOF cluster re-mapping in AliESDtrack: out of range (%i > %i)\n",fTOFcluster[i],ncl));
3400   }
3401 }
3402
3403 //____________________________________________
3404 void AliESDtrack::SortTOFcluster(){
3405   TClonesArray *tofclArray = GetESDEvent()->GetESDTOFClusters();
3406
3407   for(Int_t i=0;i<fNtofClusters-1;i++){
3408     for(Int_t j=i+1;j<fNtofClusters;j++){
3409       AliESDTOFCluster *tofcl = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[i]);
3410       Int_t index1 = -1;
3411       for(Int_t it=0;it < tofcl->GetNMatchableTracks();it++){
3412          if(tofcl->GetTrackIndex(it) == GetID()) index1 = it;
3413       }
3414       Double_t timedist1 = 10000;
3415       for(Int_t isp=0; isp< AliPID::kSPECIESC;isp++){
3416         Double_t timec = TMath::Abs(tofcl->GetTime() - tofcl->GetIntegratedTime(isp));
3417         if(timec < timedist1) timedist1 = timec;
3418       }
3419       timedist1 *= 0.03; // in cm
3420       Double_t radius1 = tofcl->GetDx(index1)*tofcl->GetDx(index1) + tofcl->GetDz(index1)*tofcl->GetDz(index1) + timedist1*timedist1;
3421
3422       AliESDTOFCluster *tofcl2 = (AliESDTOFCluster *) tofclArray->At(fTOFcluster[j]);
3423       Int_t index2 = -1;
3424       for(Int_t it=0;it < tofcl2->GetNMatchableTracks();it++){
3425          if(tofcl2->GetTrackIndex(it) == GetID()) index2 = it;
3426       }
3427       if(index1 == -1 || index2 == -1){
3428       }
3429       Double_t timedist2 = 10000;
3430       for(Int_t isp=0; isp< AliPID::kSPECIESC;isp++){
3431         Double_t timec = TMath::Abs(tofcl2->GetTime() - tofcl2->GetIntegratedTime(isp));
3432         if(timec < timedist2) timedist2 = timec;
3433       }
3434       timedist2 *= 0.03; // in cm
3435       Double_t radius2 = tofcl2->GetDx(index2)*tofcl2->GetDx(index2) + tofcl2->GetDz(index2)*tofcl2->GetDz(index2) + timedist2*timedist2;
3436
3437       if(radius2 < radius1){
3438         Int_t change = fTOFcluster[i];
3439         fTOFcluster[i] = fTOFcluster[j];
3440         fTOFcluster[j] = change;
3441       }
3442     }
3443   }
3444 }
3445
3446 //____________________________________________
3447 const AliTOFHeader* AliESDtrack::GetTOFHeader() const {
3448   return fESDEvent->GetTOFHeader();
3449 }
3450
3451 //___________________________________________
3452 void AliESDtrack::SetID(Short_t id)
3453 {
3454   // set track ID taking care about dependencies
3455   if (fNtofClusters) ReplaceTOFTrackID(fID,id);
3456   fID=id;
3457 }
3458
3459
3460 Double_t  AliESDtrack::GetdEdxInfo(Int_t regionID, Int_t calibID, Int_t qID, Int_t valueID){
3461   //
3462   // Interface to get the calibrated dEdx information
3463   // For details of arguments and return values see
3464   //     AliTPCdEdxInfo::GetdEdxInfo(Int_t regionID, Int_t calibID, Int_t valueID)
3465   //
3466   if (!fTPCdEdxInfo) return 0;
3467   if (!fIp) return 0;
3468   return fTPCdEdxInfo->GetdEdxInfo(fIp, regionID, calibID, qID, valueID);
3469 }
3470
3471
3472 Double_t AliESDtrack::GetdEdxInfoTRD(Int_t method, Double_t p0, Double_t p1, Double_t p2){
3473   //
3474   // Methods
3475   // mean values:
3476   //     0.)   linear
3477   //     1.)   logarithmic
3478   //     2.)   1/sqrt
3479   //     3.)   power()
3480   // time COG:
3481   //     4.)   linear
3482   //     5.)   logarithmic
3483   //     6.)   square
3484   Int_t nSlicesPerLayer=GetNumberOfTRDslices();
3485   Int_t nSlicesAll=GetNumberOfTRDslices()*kTRDnPlanes;
3486
3487   if (method<=3){
3488     Double_t sumAmp=0;
3489     Int_t    sumW=0;
3490     for (Int_t ibin=0; ibin<nSlicesAll; ibin++){
3491       if (fTRDslices[ibin]<=0) continue;
3492       sumW++;
3493       if (method==0) sumAmp+=fTRDslices[ibin];
3494       if (method==1) sumAmp+=TMath::Log(TMath::Abs(fTRDslices[ibin])+p0);
3495       if (method==2) sumAmp+=1/TMath::Sqrt(TMath::Abs(fTRDslices[ibin])+p0);
3496       if (method==3) sumAmp+=TMath::Power(TMath::Abs(fTRDslices[ibin])+p0,p1);
3497     }
3498     if (sumW==0) return 0;
3499     Double_t dEdx=sumAmp/sumW;
3500     if (method==1) dEdx= TMath::Exp(dEdx);
3501     if (method==2) dEdx= 1/(dEdx*dEdx);
3502     if (method==3) dEdx= TMath::Power(dEdx,1/p1);
3503     return dEdx;
3504   }
3505   if (method>3){
3506     Double_t sumWT=0;
3507     Double_t sumW=0;
3508     for (Int_t ibin=0; ibin<nSlicesAll; ibin++){
3509       if (fTRDslices[ibin]<=0) continue;
3510       Double_t time=(ibin%nSlicesPerLayer);
3511       Double_t weight=fTRDslices[ibin];
3512       if (method==5) weight=TMath::Log((weight+p0)/p0);
3513       if (method==6) weight=TMath::Power(weight+p0,p1);
3514       sumWT+=time*weight;
3515       sumW+=weight;
3516     }
3517     if (sumW<=0) return 0;
3518     Double_t meanTime=sumWT/sumW;
3519     return meanTime;
3520   }
3521   return 0;
3522 }