STEER/AliESDtrack.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 //-----------------------------------------------------------------
  16 //           Implementation of the ESD track class
  17 //   ESD = Event Summary Data
  18 //   This is the class to deal with during the phisics analysis of data
  19 //      Origin: Iouri Belikov, CERN
  20 //      e-mail: Jouri.Belikov@cern.ch
  21 //-----------------------------------------------------------------
  22
  23 #include <TMath.h>
  24 #include <TParticle.h>
  25
  26 #include "AliESDVertex.h"
  27 #include "AliESDtrack.h"
  28 #include "AliKalmanTrack.h"
  29 #include "AliLog.h"
  30 #include "AliTrackPointArray.h"
  31
  32 ClassImp(AliESDtrack)
  33
  34 void SetPIDValues(Float_t * dest, const Double_t * src, Int_t n) {
  35   // This function copies "n" PID weights from "scr" to "dest"
  36   // and normalizes their sum to 1 thus producing conditional probabilities.
  37   // The negative weights are set to 0.
  38   // In case all the weights are non-positive they are replaced by
  39   // uniform probabilities
  40
  41   if (n<=0) return;
  42
  43   Float_t uniform = 1./(Float_t)n;
  44
  45   Float_t sum = 0;
  46   for (Int_t i=0; i<n; i++)
  47     if (src[i]>=0) {
  48       sum+=src[i];
  49       dest[i] = src[i];
  50     }
  51     else {
  52       dest[i] = 0;
  53     }
  54
  55   if(sum>0)
  56     for (Int_t i=0; i<n; i++) dest[i] /= sum;
  57   else
  58     for (Int_t i=0; i<n; i++) dest[i] = uniform;
  59 }
  60
  61 //_______________________________________________________________________
  62 AliESDtrack::AliESDtrack() :
  63   AliExternalTrackParam(),
  64   fFlags(0),
  65   fLabel(0),
  66   fID(0),
  67   fTrackLength(0),
  68   fD(0),fZ(0),
  69   fCdd(0),fCdz(0),fCzz(0),
  70   fStopVertex(0),
  71   fCp(0),
  72   fCchi2(1e10),
  73   fIp(0),
  74   fTPCInner(0),
  75   fOp(0),
  76   fITSchi2(0),
  77   fITSncls(0),
  78   fITSClusterMap(0),
  79   fITSsignal(0),
  80   fITSLabel(0),
  81   fTPCchi2(0),
  82   fTPCncls(0),
  83   fTPCnclsF(0),
  84   fTPCClusterMap(159),//number of padrows
  85   fTPCSharedMap(159),//number of padrows
  86   fTPCsignal(0),
  87   fTPCsignalN(0),
  88   fTPCsignalS(0),
  89   fTPCLabel(0),
  90   fTRDchi2(0),
  91   fTRDncls(0),
  92   fTRDncls0(0),
  93   fTRDsignal(0),
  94   fTRDLabel(0),
  95   fTRDQuality(0),
  96   fTRDBudget(0),
  97   fTOFchi2(0),
  98   fTOFindex(0),
  99   fTOFCalChannel(-1),
 100   fTOFsignal(-1),
 101   fTOFsignalToT(0),
 102   fTOFsignalRaw(0),
 103   fTOFsignalDz(0),
 104   fHMPIDchi2(1e10),
 105   fHMPIDqn(-1),
 106   fHMPIDcluIdx(-1),
 107   fHMPIDsignal(-1),
 108   fHMPIDtrkTheta(-1),
 109   fHMPIDtrkPhi(-1),
 110   fHMPIDtrkX(-1),
 111   fHMPIDtrkY(-1),
 112   fHMPIDmipX(-1),
 113   fHMPIDmipY(-1),
 114   fEMCALindex(kEMCALNoMatch),
 115   fFriendTrack(new AliESDfriendTrack())
 116 {
 117   //
 118   // The default ESD constructor
 119   //
 120   Int_t i, j;
 121   for (i=0; i<AliPID::kSPECIES; i++) {
 122     fTrackTime[i]=0.;
 123     fR[i]=1.;
 124     fITSr[i]=1.;
 125     fTPCr[i]=1.;
 126     fTRDr[i]=1.;
 127     fTOFr[i]=1.;
 128     fHMPIDr[i]=1.;
 129   }
 130
 131   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
 132   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
 133   for (i=0;i<kNPlane;i++) {
 134     for (j=0;j<kNSlice;j++) {
 135       fTRDsignals[i][j]=0.;
 136     }
 137     fTRDTimBin[i]=-1;
 138   }
 139   for (i=0;i<4;i++) {fTPCPoints[i]=-1;}
 140   for (i=0;i<3;i++) {fTOFLabel[i]=-1;}
 141   for (i=0;i<10;i++) {fTOFInfo[i]=-1;}
 142 }
 143
 144 //_______________________________________________________________________
 145 AliESDtrack::AliESDtrack(const AliESDtrack& track):
 146   AliExternalTrackParam(track),
 147   fFlags(track.fFlags),
 148   fLabel(track.fLabel),
 149   fID(track.fID),
 150   fTrackLength(track.fTrackLength),
 151   fD(track.fD),fZ(track.fZ),
 152   fCdd(track.fCdd),fCdz(track.fCdz),fCzz(track.fCzz),
 153   fStopVertex(track.fStopVertex),
 154   fCp(0),
 155   fCchi2(track.fCchi2),
 156   fIp(0),
 157   fTPCInner(0),
 158   fOp(0),
 159   fITSchi2(track.fITSchi2),
 160   fITSncls(track.fITSncls),
 161   fITSClusterMap(track.fITSClusterMap),
 162   fITSsignal(track.fITSsignal),
 163   fITSLabel(track.fITSLabel),
 164   fTPCchi2(track.fTPCchi2),
 165   fTPCncls(track.fTPCncls),
 166   fTPCnclsF(track.fTPCnclsF),
 167   fTPCClusterMap(track.fTPCClusterMap),
 168   fTPCSharedMap(track.fTPCSharedMap),
 169   fTPCsignal(track.fTPCsignal),
 170   fTPCsignalN(track.fTPCsignalN),
 171   fTPCsignalS(track.fTPCsignalS),
 172   fTPCLabel(track.fTPCLabel),
 173   fTRDchi2(track.fTRDchi2),
 174   fTRDncls(track.fTRDncls),
 175   fTRDncls0(track.fTRDncls0),
 176   fTRDsignal(track.fTRDsignal),
 177   fTRDLabel(track.fTRDLabel),
 178   fTRDQuality(track.fTRDQuality),
 179   fTRDBudget(track.fTRDBudget),
 180   fTOFchi2(track.fTOFchi2),
 181   fTOFindex(track.fTOFindex),
 182   fTOFCalChannel(track.fTOFCalChannel),
 183   fTOFsignal(track.fTOFsignal),
 184   fTOFsignalToT(track.fTOFsignalToT),
 185   fTOFsignalRaw(track.fTOFsignalRaw),
 186   fTOFsignalDz(track.fTOFsignalDz),
 187   fHMPIDchi2(track.fHMPIDchi2),
 188   fHMPIDqn(track.fHMPIDqn),
 189   fHMPIDcluIdx(track.fHMPIDcluIdx),
 190   fHMPIDsignal(track.fHMPIDsignal),
 191   fHMPIDtrkTheta(track.fHMPIDtrkTheta),
 192   fHMPIDtrkPhi(track.fHMPIDtrkPhi),
 193   fHMPIDtrkX(track.fHMPIDtrkX),
 194   fHMPIDtrkY(track.fHMPIDtrkY),
 195   fHMPIDmipX(track.fHMPIDmipX),
 196   fHMPIDmipY(track.fHMPIDmipY),
 197   fEMCALindex(track.fEMCALindex),
 198   fFriendTrack(0)
 199 {
 200   //
 201   //copy constructor
 202   //
 203   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i]=track.fTrackTime[i];
 204   for (Int_t i=0;i<AliPID::kSPECIES;i++)  fR[i]=track.fR[i];
 205   //
 206   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=track.fITSr[i];
 207   //
 208   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=track.fTPCr[i];
 209   for (Int_t i=0;i<4;i++) {fTPCPoints[i]=track.fTPCPoints[i];}
 210   for (Int_t i=0; i<3;i++)   { fKinkIndexes[i]=track.fKinkIndexes[i];}
 211   for (Int_t i=0; i<3;i++)   { fV0Indexes[i]=track.fV0Indexes[i];}
 212   //
 213   for (Int_t i=0;i<kNPlane;i++) {
 214     for (Int_t j=0;j<kNSlice;j++) {
 215       fTRDsignals[i][j]=track.fTRDsignals[i][j];
 216     }
 217     fTRDTimBin[i]=track.fTRDTimBin[i];
 218   }
 219   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i]=track.fTRDr[i];
 220   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i]=track.fTOFr[i];
 221   for (Int_t i=0;i<3;i++) fTOFLabel[i]=track.fTOFLabel[i];
 222   for (Int_t i=0;i<10;i++) fTOFInfo[i]=track.fTOFInfo[i];
 223   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i]=track.fHMPIDr[i];
 224
 225   if (track.fCp) fCp=new AliExternalTrackParam(*track.fCp);
 226   if (track.fIp) fIp=new AliExternalTrackParam(*track.fIp);
 227   if (track.fTPCInner) fTPCInner=new AliExternalTrackParam(*track.fTPCInner);
 228   if (track.fOp) fOp=new AliExternalTrackParam(*track.fOp);
 229
 230   if (track.fFriendTrack) fFriendTrack=new AliESDfriendTrack(*(track.fFriendTrack));
 231 }
 232
 233 //_______________________________________________________________________
 234 AliESDtrack::AliESDtrack(TParticle * part) :
 235   AliExternalTrackParam(),
 236   fFlags(0),
 237   fLabel(0),
 238   fID(0),
 239   fTrackLength(0),
 240   fD(0),fZ(0),
 241   fCdd(0),fCdz(0),fCzz(0),
 242   fStopVertex(0),
 243   fCp(0),
 244   fCchi2(1e10),
 245   fIp(0),
 246   fTPCInner(0),
 247   fOp(0),
 248   fITSchi2(0),
 249   fITSncls(0),
 250   fITSClusterMap(0),
 251   fITSsignal(0),
 252   fITSLabel(0),
 253   fTPCchi2(0),
 254   fTPCncls(0),
 255   fTPCnclsF(0),
 256   fTPCClusterMap(159),//number of padrows
 257   fTPCSharedMap(159),//number of padrows
 258   fTPCsignal(0),
 259   fTPCsignalN(0),
 260   fTPCsignalS(0),
 261   fTPCLabel(0),
 262   fTRDchi2(0),
 263   fTRDncls(0),
 264   fTRDncls0(0),
 265   fTRDsignal(0),
 266   fTRDLabel(0),
 267   fTRDQuality(0),
 268   fTRDBudget(0),
 269   fTOFchi2(0),
 270   fTOFindex(0),
 271   fTOFCalChannel(-1),
 272   fTOFsignal(-1),
 273   fTOFsignalToT(0),
 274   fTOFsignalRaw(0),
 275   fTOFsignalDz(0),
 276   fHMPIDchi2(1e10),
 277   fHMPIDqn(-1),
 278   fHMPIDcluIdx(-1),
 279   fHMPIDsignal(-1),
 280   fHMPIDtrkTheta(-1),
 281   fHMPIDtrkPhi(-1),
 282   fHMPIDtrkX(-1),
 283   fHMPIDtrkY(-1),
 284   fHMPIDmipX(-1),
 285   fHMPIDmipY(-1),
 286   fEMCALindex(kEMCALNoMatch),
 287   fFriendTrack(0)
 288 {
 289   //
 290   // ESD track from TParticle
 291   //
 292
 293   // Reset all the arrays
 294   Int_t i, j;
 295   for (i=0; i<AliPID::kSPECIES; i++) {
 296     fTrackTime[i]=0.;
 297     fR[i]=0.;
 298     fITSr[i]=0.;
 299     fTPCr[i]=0.;
 300     fTRDr[i]=0.;
 301     fTOFr[i]=0.;
 302     fHMPIDr[i]=0.;
 303   }
 304
 305   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
 306   for (i=0; i<3; i++)   { fV0Indexes[i]=-1;}
 307   for (i=0;i<kNPlane;i++) {
 308     for (j=0;j<kNSlice;j++) {
 309       fTRDsignals[i][j]=0.;
 310     }
 311     fTRDTimBin[i]=-1;
 312   }
 313   for (i=0;i<4;i++) {fTPCPoints[i]=-1;}
 314   for (i=0;i<3;i++) {fTOFLabel[i]=-1;}
 315   for (i=0;i<10;i++) {fTOFInfo[i]=-1;}
 316
 317   // Calculate the AliExternalTrackParam content
 318
 319   Double_t xref;
 320   Double_t alpha;
 321   Double_t param[5];
 322   Double_t covar[15];
 323
 324   // Calculate alpha: the rotation angle of the corresponding local system (TPC sector)
 325   alpha = part->Phi()*180./TMath::Pi();
 326   if (alpha<0) alpha+= 360.;
 327   if (alpha>360) alpha -= 360.;
 328
 329   Int_t sector = (Int_t)(alpha/20.);
 330   alpha = 10. + 20.*sector;
 331   alpha /= 180;
 332   alpha *= TMath::Pi();
 333
 334   // Covariance matrix: no errors, the parameters are exact
 335   for (i=0; i<15; i++) covar[i]=0.;
 336
 337   // Get the vertex of origin and the momentum
 338   TVector3 ver(part->Vx(),part->Vy(),part->Vz());
 339   TVector3 mom(part->Px(),part->Py(),part->Pz());
 340
 341   // Rotate to the local coordinate system (TPC sector)
 342   ver.RotateZ(-alpha);
 343   mom.RotateZ(-alpha);
 344
 345   // X of the referense plane
 346   xref = ver.X();
 347
 348   Int_t pdgCode = part->GetPdgCode();
 349
 350   Double_t charge =
 351     TDatabasePDG::Instance()->GetParticle(pdgCode)->Charge();
 352
 353   param[0] = ver.Y();
 354   param[1] = ver.Z();
 355   param[2] = TMath::Sin(mom.Phi());
 356   param[3] = mom.Pz()/mom.Pt();
 357   param[4] = TMath::Sign(1/mom.Pt(),charge);
 358
 359   // Set AliExternalTrackParam
 360   Set(xref, alpha, param, covar);
 361
 362   // Set the PID
 363   Int_t indexPID = 99;
 364
 365   switch (TMath::Abs(pdgCode)) {
 366
 367   case  11: // electron
 368     indexPID = 0;
 369     break;
 370
 371   case 13: // muon
 372     indexPID = 1;
 373     break;
 374
 375   case 211: // pion
 376     indexPID = 2;
 377     break;
 378
 379   case 321: // kaon
 380     indexPID = 3;
 381     break;
 382
 383   case 2212: // proton
 384     indexPID = 4;
 385     break;
 386
 387   default:
 388     break;
 389   }
 390
 391   // If the particle is not e,mu,pi,K or p the PID probabilities are set to 0
 392   if (indexPID < AliPID::kSPECIES) {
 393     fR[indexPID]=1.;
 394     fITSr[indexPID]=1.;
 395     fTPCr[indexPID]=1.;
 396     fTRDr[indexPID]=1.;
 397     fTOFr[indexPID]=1.;
 398     fHMPIDr[indexPID]=1.;
 399
 400   }
 401   // AliESD track label
 402   SetLabel(part->GetUniqueID());
 403
 404 }
 405
 406 //_______________________________________________________________________
 407 AliESDtrack::~AliESDtrack(){
 408   //
 409   // This is destructor according Coding Conventrions
 410   //
 411   //printf("Delete track\n");
 412   delete fIp;
 413   delete fTPCInner;
 414   delete fOp;
 415   delete fCp;
 416   delete fFriendTrack;
 417 }
 418
 419 void AliESDtrack::AddCalibObject(TObject * object){
 420   //
 421   // add calib object to the list
 422   //
 423   if (!fFriendTrack) fFriendTrack  = new AliESDfriendTrack;
 424   fFriendTrack->AddCalibObject(object);
 425 }
 426
 427 TObject *  AliESDtrack::GetCalibObject(Int_t index){
 428   //
 429   // return calib objct at given position
 430   //
 431   if (!fFriendTrack) return 0;
 432   return fFriendTrack->GetCalibObject(index);
 433 }
 434
 435
 436 //_______________________________________________________________________
 437 void AliESDtrack::MakeMiniESDtrack(){
 438   // Resets everything except
 439   // fFlags: Reconstruction status flags
 440   // fLabel: Track label
 441   // fID:  Unique ID of the track
 442   // fD: Impact parameter in XY-plane
 443   // fZ: Impact parameter in Z
 444   // fR[AliPID::kSPECIES]: combined "detector response probability"
 445   // Running track parameters in the base class (AliExternalTrackParam)
 446
 447   fTrackLength = 0;
 448   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i] = 0;
 449   fStopVertex = 0;
 450
 451   // Reset track parameters constrained to the primary vertex
 452   fCp = 0;
 453   fCchi2 = 0;
 454
 455   // Reset track parameters at the inner wall of TPC
 456   fIp = 0;
 457   fTPCInner=0;
 458   // Reset track parameters at the inner wall of the TRD
 459   fOp = 0;
 460
 461   // Reset ITS track related information
 462   fITSchi2 = 0;
 463   fITSncls = 0;
 464   fITSClusterMap=0;
 465   fITSsignal = 0;
 466   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=0;
 467   fITSLabel = 0;
 468
 469   // Reset TPC related track information
 470   fTPCchi2 = 0;
 471   fTPCncls = 0;
 472   fTPCnclsF = 0;
 473   fTPCClusterMap = 0;
 474   fTPCSharedMap = 0;
 475   fTPCsignal= 0;
 476   fTPCsignalS= 0;
 477   fTPCsignalN= 0;
 478   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0;
 479   fTPCLabel=0;
 480   for (Int_t i=0;i<4;i++) fTPCPoints[i] = 0;
 481   for (Int_t i=0; i<3;i++)   fKinkIndexes[i] = 0;
 482   for (Int_t i=0; i<3;i++)   fV0Indexes[i] = 0;
 483
 484   // Reset TRD related track information
 485   fTRDchi2 = 0;
 486   fTRDncls = 0;
 487   fTRDncls0 = 0;
 488   fTRDsignal = 0;
 489   for (Int_t i=0;i<kNPlane;i++) {
 490     for (Int_t j=0;j<kNSlice;j++) {
 491       fTRDsignals[i][j] = 0;
 492     }
 493     fTRDTimBin[i]  = 0;
 494   }
 495   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0;
 496   fTRDLabel = 0;
 497   fTRDQuality  = 0;
 498   fTRDBudget  = 0;
 499
 500   // Reset TOF related track information
 501   fTOFchi2 = 0;
 502   fTOFindex = 0;
 503   fTOFsignal = 0;
 504   fTOFCalChannel = -1;
 505   fTOFsignalToT = 0;
 506   fTOFsignalRaw = 0;
 507   fTOFsignalDz = 0;
 508   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i] = 0;
 509   for (Int_t i=0;i<3;i++) fTOFLabel[i] = 0;
 510   for (Int_t i=0;i<10;i++) fTOFInfo[i] = 0;
 511
 512   // Reset HMPID related track information
 513   fHMPIDchi2 = 0;
 514   fHMPIDqn = -1;
 515   fHMPIDcluIdx = -1;
 516   fHMPIDsignal = 0;
 517   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i] = 0;
 518   fHMPIDtrkTheta = -1;
 519   fHMPIDtrkPhi = -1;
 520   fHMPIDtrkX = -1;
 521   fHMPIDtrkY = -1;
 522   fHMPIDmipX = -1;
 523   fHMPIDmipY = -1;
 524   fEMCALindex = kEMCALNoMatch;
 525
 526   delete fFriendTrack; fFriendTrack = 0;
 527 }
 528 //_______________________________________________________________________
 529 Double_t AliESDtrack::GetMass() const {
 530   // Returns the mass of the most probable particle type
 531   Float_t max=0.;
 532   Int_t k=-1;
 533   for (Int_t i=0; i<AliPID::kSPECIES; i++) {
 534     if (fR[i]>max) {k=i; max=fR[i];}
 535   }
 536   if (k==0) { // dE/dx "crossing points" in the TPC
 537      Double_t p=GetP();
 538      if ((p>0.38)&&(p<0.48))
 539         if (fR[0]<fR[3]*10.) return AliPID::ParticleMass(AliPID::kKaon);
 540      if ((p>0.75)&&(p<0.85))
 541         if (fR[0]<fR[4]*10.) return AliPID::ParticleMass(AliPID::kProton);
 542      return 0.00051;
 543   }
 544   if (k==1) return AliPID::ParticleMass(AliPID::kMuon);
 545   if (k==2||k==-1) return AliPID::ParticleMass(AliPID::kPion);
 546   if (k==3) return AliPID::ParticleMass(AliPID::kKaon);
 547   if (k==4) return AliPID::ParticleMass(AliPID::kProton);
 548   AliWarning("Undefined mass !");
 549   return AliPID::ParticleMass(AliPID::kPion);
 550 }
 551
 552 //_______________________________________________________________________
 553 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
 554   //
 555   // This function updates track's running parameters
 556   //
 557   Int_t *index=0;
 558   Bool_t rc=kTRUE;
 559
 560   SetStatus(flags);
 561   fLabel=t->GetLabel();
 562
 563   if (t->IsStartedTimeIntegral()) {
 564     SetStatus(kTIME);
 565     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
 566     SetIntegratedLength(t->GetIntegratedLength());
 567   }
 568
 569   Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 570
 571   switch (flags) {
 572
 573   case kITSin: case kITSout: case kITSrefit:
 574     fITSncls=t->GetNumberOfClusters();
 575     index=fFriendTrack->GetITSindices();
 576     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
 577         index[i]=t->GetClusterIndex(i);
 578         if (i<fITSncls) {
 579            Int_t l=(index[i] & 0xf0000000) >> 28;
 580            SETBIT(fITSClusterMap,l);
 581         }
 582     }
 583     fITSchi2=t->GetChi2();
 584     fITSsignal=t->GetPIDsignal();
 585     fITSLabel = t->GetLabel();
 586     break;
 587
 588   case kTPCin: case kTPCrefit:
 589     fTPCLabel = t->GetLabel();
 590     if (flags==kTPCin)  fTPCInner=new AliExternalTrackParam(*t);
 591     if (!fIp) fIp=new AliExternalTrackParam(*t);
 592     else
 593       fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 594   case kTPCout:
 595     index=fFriendTrack->GetTPCindices();
 596     if (flags & kTPCout){
 597       if (!fOp) fOp=new AliExternalTrackParam(*t);
 598       else
 599         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 600     }
 601     fTPCncls=t->GetNumberOfClusters();
 602     fTPCchi2=t->GetChi2();
 603
 604      {//prevrow must be declared in separate namespace, otherwise compiler cries:
 605       //"jump to case label crosses initialization of `Int_t prevrow'"
 606        Int_t prevrow = -1;
 607        //       for (Int_t i=0;i<fTPCncls;i++)
 608        for (Int_t i=0;i<AliESDfriendTrack::kMaxTPCcluster;i++)
 609         {
 610           index[i]=t->GetClusterIndex(i);
 611           Int_t idx = index[i];
 612
 613           if (idx<0) continue;
 614
 615           // Piotr's Cluster Map for HBT
 616           // ### please change accordingly if cluster array is changing
 617           // to "New TPC Tracking" style (with gaps in array)
 618           Int_t sect = (idx&0xff000000)>>24;
 619           Int_t row = (idx&0x00ff0000)>>16;
 620           if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
 621
 622           fTPCClusterMap.SetBitNumber(row,kTRUE);
 623
 624           //Fill the gap between previous row and this row with 0 bits
 625           //In case  ###  pleas change it as well - just set bit 0 in case there
 626           //is no associated clusters for current "i"
 627           if (prevrow < 0)
 628            {
 629              prevrow = row;//if previous bit was not assigned yet == this is the first one
 630            }
 631           else
 632            { //we don't know the order (inner to outer or reverse)
 633              //just to be save in case it is going to change
 634              Int_t n = 0, m = 0;
 635              if (prevrow < row)
 636               {
 637                 n = prevrow;
 638                 m = row;
 639               }
 640              else
 641               {
 642                 n = row;
 643                 m = prevrow;
 644               }
 645
 646              for (Int_t j = n+1; j < m; j++)
 647               {
 648                 fTPCClusterMap.SetBitNumber(j,kFALSE);
 649               }
 650              prevrow = row;
 651            }
 652           // End Of Piotr's Cluster Map for HBT
 653         }
 654      }
 655     fTPCsignal=t->GetPIDsignal();
 656     break;
 657
 658   case kTRDout: case kTRDin: case kTRDrefit:
 659     index=fFriendTrack->GetTRDindices();
 660     fTRDLabel = t->GetLabel();
 661     fTRDncls=t->GetNumberOfClusters();
 662     fTRDchi2=t->GetChi2();
 663     for (Int_t i=0;i<fTRDncls;i++) index[i]=t->GetClusterIndex(i);
 664     fTRDsignal=t->GetPIDsignal();
 665     break;
 666   case kTRDbackup:
 667     if (!fOp) fOp=new AliExternalTrackParam(*t);
 668     else
 669       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 670     fTRDncls0 = t->GetNumberOfClusters();
 671     break;
 672   case kTOFin:
 673     break;
 674   case kTOFout:
 675     break;
 676   case kTRDStop:
 677     break;
 678   default:
 679     AliError("Wrong flag !");
 680     return kFALSE;
 681   }
 682
 683   return rc;
 684 }
 685
 686 //_______________________________________________________________________
 687 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
 688   //---------------------------------------------------------------------
 689   // This function returns external representation of the track parameters
 690   //---------------------------------------------------------------------
 691   x=GetX();
 692   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
 693 }
 694
 695 //_______________________________________________________________________
 696 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
 697   //---------------------------------------------------------------------
 698   // This function returns external representation of the cov. matrix
 699   //---------------------------------------------------------------------
 700   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
 701 }
 702
 703 //_______________________________________________________________________
 704 Bool_t AliESDtrack::GetConstrainedExternalParameters
 705                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 706   //---------------------------------------------------------------------
 707   // This function returns the constrained external track parameters
 708   //---------------------------------------------------------------------
 709   if (!fCp) return kFALSE;
 710   alpha=fCp->GetAlpha();
 711   x=fCp->GetX();
 712   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
 713   return kTRUE;
 714 }
 715
 716 //_______________________________________________________________________
 717 Bool_t
 718 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
 719   //---------------------------------------------------------------------
 720   // This function returns the constrained external cov. matrix
 721   //---------------------------------------------------------------------
 722   if (!fCp) return kFALSE;
 723   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
 724   return kTRUE;
 725 }
 726
 727 Bool_t
 728 AliESDtrack::GetInnerExternalParameters
 729                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 730   //---------------------------------------------------------------------
 731   // This function returns external representation of the track parameters
 732   // at the inner layer of TPC
 733   //---------------------------------------------------------------------
 734   if (!fIp) return kFALSE;
 735   alpha=fIp->GetAlpha();
 736   x=fIp->GetX();
 737   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
 738   return kTRUE;
 739 }
 740
 741 Bool_t
 742 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
 743  //---------------------------------------------------------------------
 744  // This function returns external representation of the cov. matrix
 745  // at the inner layer of TPC
 746  //---------------------------------------------------------------------
 747   if (!fIp) return kFALSE;
 748   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
 749   return kTRUE;
 750 }
 751
 752 Bool_t
 753 AliESDtrack::GetOuterExternalParameters
 754                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 755   //---------------------------------------------------------------------
 756   // This function returns external representation of the track parameters
 757   // at the inner layer of TRD
 758   //---------------------------------------------------------------------
 759   if (!fOp) return kFALSE;
 760   alpha=fOp->GetAlpha();
 761   x=fOp->GetX();
 762   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
 763   return kTRUE;
 764 }
 765
 766 Bool_t
 767 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
 768  //---------------------------------------------------------------------
 769  // This function returns external representation of the cov. matrix
 770  // at the inner layer of TRD
 771  //---------------------------------------------------------------------
 772   if (!fOp) return kFALSE;
 773   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
 774   return kTRUE;
 775 }
 776
 777 Int_t AliESDtrack::GetNcls(Int_t idet) const
 778 {
 779   // Get number of clusters by subdetector index
 780   //
 781   Int_t ncls = 0;
 782   switch(idet){
 783   case 0:
 784     ncls = fITSncls;
 785     break;
 786   case 1:
 787     ncls = fTPCncls;
 788     break;
 789   case 2:
 790     ncls = fTRDncls;
 791     break;
 792   case 3:
 793     if (fTOFindex != 0)
 794       ncls = 1;
 795     break;
 796   default:
 797     break;
 798   }
 799   return ncls;
 800 }
 801
 802 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
 803 {
 804   // Get cluster index array by subdetector index
 805   //
 806   Int_t ncls = 0;
 807   switch(idet){
 808   case 0:
 809     ncls = GetITSclusters(idx);
 810     break;
 811   case 1:
 812     ncls = GetTPCclusters(idx);
 813     break;
 814   case 2:
 815     ncls = GetTRDclusters(idx);
 816     break;
 817   case 3:
 818     if (fTOFindex != 0) {
 819       idx[0] = GetTOFcluster();
 820       ncls = 1;
 821     }
 822     break;
 823   default:
 824     break;
 825   }
 826   return ncls;
 827 }
 828
 829 //_______________________________________________________________________
 830 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
 831   // Returns the array with integrated times for each particle hypothesis
 832   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
 833 }
 834
 835 //_______________________________________________________________________
 836 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
 837   // Sets the array with integrated times for each particle hypotesis
 838   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
 839 }
 840
 841 //_______________________________________________________________________
 842 void AliESDtrack::SetITSpid(const Double_t *p) {
 843   // Sets values for the probability of each particle type (in ITS)
 844   SetPIDValues(fITSr,p,AliPID::kSPECIES);
 845   SetStatus(AliESDtrack::kITSpid);
 846 }
 847
 848 //_______________________________________________________________________
 849 void AliESDtrack::GetITSpid(Double_t *p) const {
 850   // Gets the probability of each particle type (in ITS)
 851   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
 852 }
 853
 854 //_______________________________________________________________________
 855 Int_t AliESDtrack::GetITSclusters(Int_t *idx) const {
 856   //---------------------------------------------------------------------
 857   // This function returns indices of the assgined ITS clusters
 858   //---------------------------------------------------------------------
 859   if (idx!=0) {
 860      Int_t *index=fFriendTrack->GetITSindices();
 861      for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) idx[i]=index[i];
 862   }
 863   return fITSncls;
 864 }
 865
 866 //_______________________________________________________________________
 867 Int_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
 868   //---------------------------------------------------------------------
 869   // This function returns indices of the assgined ITS clusters
 870   //---------------------------------------------------------------------
 871   if (idx!=0) {
 872     Int_t *index=fFriendTrack->GetTPCindices();
 873     for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
 874   }
 875   return fTPCncls;
 876 }
 877
 878 Float_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
 879   //
 880   // GetDensity of the clusters on given region between row0 and row1
 881   // Dead zone effect takin into acoount
 882   //
 883   Int_t good  = 0;
 884   Int_t found = 0;
 885   //
 886   Int_t *index=fFriendTrack->GetTPCindices();
 887   for (Int_t i=row0;i<=row1;i++){
 888     Int_t idx = index[i];
 889     if (idx!=-1)  good++;             // track outside of dead zone
 890     if (idx>0)    found++;
 891   }
 892   Float_t density=0.5;
 893   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
 894   return density;
 895 }
 896
 897 //_______________________________________________________________________
 898 void AliESDtrack::SetTPCpid(const Double_t *p) {
 899   // Sets values for the probability of each particle type (in TPC)
 900   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
 901   SetStatus(AliESDtrack::kTPCpid);
 902 }
 903
 904 //_______________________________________________________________________
 905 void AliESDtrack::GetTPCpid(Double_t *p) const {
 906   // Gets the probability of each particle type (in TPC)
 907   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
 908 }
 909
 910 //_______________________________________________________________________
 911 Int_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
 912   //---------------------------------------------------------------------
 913   // This function returns indices of the assgined TRD clusters
 914   //---------------------------------------------------------------------
 915   if (idx!=0) {
 916      Int_t *index=fFriendTrack->GetTRDindices();
 917      for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
 918   }
 919   return fTRDncls;
 920 }
 921
 922 //_______________________________________________________________________
 923 void AliESDtrack::SetTRDpid(const Double_t *p) {
 924   // Sets values for the probability of each particle type (in TRD)
 925   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
 926   SetStatus(AliESDtrack::kTRDpid);
 927 }
 928
 929 //_______________________________________________________________________
 930 void AliESDtrack::GetTRDpid(Double_t *p) const {
 931   // Gets the probability of each particle type (in TRD)
 932   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
 933 }
 934
 935 //_______________________________________________________________________
 936 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
 937 {
 938   // Sets the probability of particle type iSpecies to p (in TRD)
 939   fTRDr[iSpecies] = p;
 940 }
 941
 942 Float_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
 943 {
 944   // Returns the probability of particle type iSpecies (in TRD)
 945   return fTRDr[iSpecies];
 946 }
 947
 948 //_______________________________________________________________________
 949 void AliESDtrack::SetTOFpid(const Double_t *p) {
 950   // Sets the probability of each particle type (in TOF)
 951   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
 952   SetStatus(AliESDtrack::kTOFpid);
 953 }
 954
 955 //_______________________________________________________________________
 956 void AliESDtrack::SetTOFLabel(const Int_t *p) {
 957   // Sets  (in TOF)
 958   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
 959 }
 960
 961 //_______________________________________________________________________
 962 void AliESDtrack::GetTOFpid(Double_t *p) const {
 963   // Gets probabilities of each particle type (in TOF)
 964   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
 965 }
 966
 967 //_______________________________________________________________________
 968 void AliESDtrack::GetTOFLabel(Int_t *p) const {
 969   // Gets (in TOF)
 970   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
 971 }
 972
 973 //_______________________________________________________________________
 974 void AliESDtrack::GetTOFInfo(Float_t *info) const {
 975   // Gets (in TOF)
 976   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
 977 }
 978
 979 //_______________________________________________________________________
 980 void AliESDtrack::SetTOFInfo(Float_t*info) {
 981   // Gets (in TOF)
 982   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
 983 }
 984
 985
 986
 987 //_______________________________________________________________________
 988 void AliESDtrack::SetHMPIDpid(const Double_t *p) {
 989   // Sets the probability of each particle type (in HMPID)
 990   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
 991   SetStatus(AliESDtrack::kHMPIDpid);
 992 }
 993
 994 //_______________________________________________________________________
 995 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
 996   // Gets probabilities of each particle type (in HMPID)
 997   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
 998 }
 999
1000
1001
1002 //_______________________________________________________________________
1003 void AliESDtrack::SetESDpid(const Double_t *p) {
1004   // Sets the probability of each particle type for the ESD track
1005   SetPIDValues(fR,p,AliPID::kSPECIES);
1006   SetStatus(AliESDtrack::kESDpid);
1007 }
1008
1009 //_______________________________________________________________________
1010 void AliESDtrack::GetESDpid(Double_t *p) const {
1011   // Gets probability of each particle type for the ESD track
1012   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
1013 }
1014
1015 //_______________________________________________________________________
1016 Bool_t AliESDtrack::RelateToVertex
1017 (const AliESDVertex *vtx, Double_t b, Double_t maxd) {
1018   //
1019   // Try to relate this track to the vertex "vtx",
1020   // if the (rough) transverse impact parameter is not bigger then "maxd".
1021   //            Magnetic field is "b" (kG).
1022   //
1023   // a) The track gets extapolated to the DCA to the vertex.
1024   // b) The impact parameters and their covariance matrix are calculated.
1025   // c) An attempt to constrain this track to the vertex is done.
1026   //
1027   //    In the case of success, the returned value is kTRUE
1028   //    (otherwise, it's kFALSE)
1029   //
1030
1031   if (!vtx) return kFALSE;
1032
1033   Double_t alpha=GetAlpha();
1034   Double_t sn=TMath::Sin(alpha), cs=TMath::Cos(alpha);
1035   Double_t x=GetX(), y=GetParameter()[0], snp=GetParameter()[2];
1036   Double_t xv= vtx->GetXv()*cs + vtx->GetYv()*sn;
1037   Double_t yv=-vtx->GetXv()*sn + vtx->GetYv()*cs, zv=vtx->GetZv();
1038   x-=xv; y-=yv;
1039
1040   //Estimate the impact parameter neglecting the track curvature
1041   Double_t d=TMath::Abs(x*snp - y*TMath::Sqrt(1.- snp*snp));
1042   if (d > maxd) return kFALSE;
1043
1044   //Propagate to the DCA
1045   Double_t crv=kB2C*b*GetParameter()[4];
1046   if (TMath::Abs(b) < kAlmost0Field) crv=0.;
1047
1048   Double_t tgfv=-(crv*x - snp)/(crv*y + TMath::Sqrt(1.-snp*snp));
1049   sn=tgfv/TMath::Sqrt(1.+ tgfv*tgfv);
1050   if (TMath::Abs(tgfv)>0.) cs = sn/tgfv;
1051   else cs=1.;
1052
1053   x = xv*cs + yv*sn;
1054   yv=-xv*sn + yv*cs; xv=x;
1055
1056   if (!Propagate(alpha+TMath::ASin(sn),xv,b)) return kFALSE;
1057
1058   fD = GetParameter()[0] - yv;
1059   fZ = GetParameter()[1] - zv;
1060
1061   Double_t cov[6]; vtx->GetCovMatrix(cov);
1062
1063   //***** Improvements by A.Dainese
1064   alpha=GetAlpha(); sn=TMath::Sin(alpha); cs=TMath::Cos(alpha);
1065   Double_t s2ylocvtx = cov[0]*sn*sn + cov[2]*cs*cs - 2.*cov[1]*cs*sn;
1066   fCdd = GetCovariance()[0] + s2ylocvtx;   // neglecting correlations
1067   fCdz = GetCovariance()[1];               // between (x,y) and z
1068   fCzz = GetCovariance()[2] + cov[5];      // in vertex's covariance matrix
1069   //*****
1070
1071   {//Try to constrain
1072     Double_t p[2]={yv,zv}, c[3]={cov[2],0.,cov[5]};
1073     Double_t chi2=GetPredictedChi2(p,c);
1074
1075     if (chi2>77.) return kFALSE;
1076
1077     AliExternalTrackParam tmp(*this);
1078     if (!tmp.Update(p,c)) return kFALSE;
1079
1080     fCchi2=chi2;
1081     if (!fCp) fCp=new AliExternalTrackParam();
1082     new (fCp) AliExternalTrackParam(tmp);
1083   }
1084
1085   return kTRUE;
1086 }
1087
1088 //_______________________________________________________________________
1089 void AliESDtrack::Print(Option_t *) const {
1090   // Prints info on the track
1091
1092   printf("ESD track info\n") ;
1093   Double_t p[AliPID::kSPECIESN] ;
1094   Int_t index = 0 ;
1095   if( IsOn(kITSpid) ){
1096     printf("From ITS: ") ;
1097     GetITSpid(p) ;
1098     for(index = 0 ; index < AliPID::kSPECIES; index++)
1099       printf("%f, ", p[index]) ;
1100     printf("\n           signal = %f\n", GetITSsignal()) ;
1101   }
1102   if( IsOn(kTPCpid) ){
1103     printf("From TPC: ") ;
1104     GetTPCpid(p) ;
1105     for(index = 0 ; index < AliPID::kSPECIES; index++)
1106       printf("%f, ", p[index]) ;
1107     printf("\n           signal = %f\n", GetTPCsignal()) ;
1108   }
1109   if( IsOn(kTRDpid) ){
1110     printf("From TRD: ") ;
1111     GetTRDpid(p) ;
1112     for(index = 0 ; index < AliPID::kSPECIES; index++)
1113       printf("%f, ", p[index]) ;
1114     printf("\n           signal = %f\n", GetTRDsignal()) ;
1115   }
1116   if( IsOn(kTOFpid) ){
1117     printf("From TOF: ") ;
1118     GetTOFpid(p) ;
1119     for(index = 0 ; index < AliPID::kSPECIES; index++)
1120       printf("%f, ", p[index]) ;
1121     printf("\n           signal = %f\n", GetTOFsignal()) ;
1122   }
1123   if( IsOn(kHMPIDpid) ){
1124     printf("From HMPID: ") ;
1125     GetHMPIDpid(p) ;
1126     for(index = 0 ; index < AliPID::kSPECIES; index++)
1127       printf("%f, ", p[index]) ;
1128     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
1129   }
1130 }
1131