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(Double_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   fCp(0),
  65   fIp(0),
  66   fTPCInner(0),
  67   fOp(0),
  68   fFriendTrack(new AliESDfriendTrack()),
  69   fTPCClusterMap(159),//number of padrows
  70   fTPCSharedMap(159),//number of padrows
  71   fFlags(0),
  72   fID(0),
  73   fLabel(0),
  74   fITSLabel(0),
  75   fTPCLabel(0),
  76   fTRDLabel(0),
  77   fTOFCalChannel(0),
  78   fTOFindex(0),
  79   fHMPIDqn(0),
  80   fHMPIDcluIdx(0),
  81   fEMCALindex(kEMCALNoMatch),
  82   fHMPIDtrkTheta(0),
  83   fHMPIDtrkPhi(0),
  84   fHMPIDsignal(0),
  85   fTrackLength(0),
  86   fD(0),fZ(0),
  87   fCdd(0),fCdz(0),fCzz(0),
  88   fCchi2(0),
  89   fITSchi2(0),
  90   fTPCchi2(0),
  91   fTRDchi2(0),
  92   fTOFchi2(0),
  93   fHMPIDchi2(0),
  94   fITSsignal(0),
  95   fTPCsignal(0),
  96   fTPCsignalS(0),
  97   fTRDsignal(0),
  98   fTRDQuality(0),
  99   fTRDBudget(0),
 100   fTOFsignal(0),
 101   fTOFsignalToT(0),
 102   fTOFsignalRaw(0),
 103   fTOFsignalDz(0),
 104   fHMPIDtrkX(0),
 105   fHMPIDtrkY(0),
 106   fHMPIDmipX(0),
 107   fHMPIDmipY(0),
 108   fTPCncls(0),
 109   fTPCnclsF(0),
 110   fTPCsignalN(0),
 111   fITSncls(0),
 112   fITSClusterMap(0),
 113   fTRDncls(0),
 114   fTRDncls0(0),
 115   fTRDpidQuality(0)
 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]=0.;
 124     fITSr[i]=0.;
 125     fTPCr[i]=0.;
 126     fTRDr[i]=0.;
 127     fTOFr[i]=0.;
 128     fHMPIDr[i]=0.;
 129   }
 130
 131   for (i=0; i<3; i++)   { fKinkIndexes[i]=0;}
 132   for (i=0; i<3; i++)   { fV0Indexes[i]=0;}
 133   for (i=0;i<kNPlane;i++) {
 134     for (j=0;j<kNSlice;j++) {
 135       fTRDsignals[i][j]=0.;
 136     }
 137     fTRDTimBin[i]=0;
 138   }
 139   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
 140   for (i=0;i<3;i++) {fTOFLabel[i]=0;}
 141   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
 142 }
 143
 144 //_______________________________________________________________________
 145 AliESDtrack::AliESDtrack(const AliESDtrack& track):
 146   AliExternalTrackParam(track),
 147   fCp(0),
 148   fIp(0),
 149   fTPCInner(0),
 150   fOp(0),
 151   fFriendTrack(0),
 152   fTPCClusterMap(track.fTPCClusterMap),
 153   fTPCSharedMap(track.fTPCSharedMap),
 154   fFlags(track.fFlags),
 155   fID(track.fID),
 156   fLabel(track.fLabel),
 157   fITSLabel(track.fITSLabel),
 158   fTPCLabel(track.fTPCLabel),
 159   fTRDLabel(track.fTRDLabel),
 160   fTOFCalChannel(track.fTOFCalChannel),
 161   fTOFindex(track.fTOFindex),
 162   fHMPIDqn(track.fHMPIDqn),
 163   fHMPIDcluIdx(track.fHMPIDcluIdx),
 164   fEMCALindex(track.fEMCALindex),
 165   fHMPIDtrkTheta(track.fHMPIDtrkTheta),
 166   fHMPIDtrkPhi(track.fHMPIDtrkPhi),
 167   fHMPIDsignal(track.fHMPIDsignal),
 168   fTrackLength(track.fTrackLength),
 169   fD(track.fD),fZ(track.fZ),
 170   fCdd(track.fCdd),fCdz(track.fCdz),fCzz(track.fCzz),
 171   fCchi2(track.fCchi2),
 172   fITSchi2(track.fITSchi2),
 173   fTPCchi2(track.fTPCchi2),
 174   fTRDchi2(track.fTRDchi2),
 175   fTOFchi2(track.fTOFchi2),
 176   fHMPIDchi2(track.fHMPIDchi2),
 177   fITSsignal(track.fITSsignal),
 178   fTPCsignal(track.fTPCsignal),
 179   fTPCsignalS(track.fTPCsignalS),
 180   fTRDsignal(track.fTRDsignal),
 181   fTRDQuality(track.fTRDQuality),
 182   fTRDBudget(track.fTRDBudget),
 183   fTOFsignal(track.fTOFsignal),
 184   fTOFsignalToT(track.fTOFsignalToT),
 185   fTOFsignalRaw(track.fTOFsignalRaw),
 186   fTOFsignalDz(track.fTOFsignalDz),
 187   fHMPIDtrkX(track.fHMPIDtrkX),
 188   fHMPIDtrkY(track.fHMPIDtrkY),
 189   fHMPIDmipX(track.fHMPIDmipX),
 190   fHMPIDmipY(track.fHMPIDmipY),
 191   fTPCncls(track.fTPCncls),
 192   fTPCnclsF(track.fTPCnclsF),
 193   fTPCsignalN(track.fTPCsignalN),
 194   fITSncls(track.fITSncls),
 195   fITSClusterMap(track.fITSClusterMap),
 196   fTRDncls(track.fTRDncls),
 197   fTRDncls0(track.fTRDncls0),
 198   fTRDpidQuality(track.fTRDpidQuality)
 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   fCp(0),
 237   fIp(0),
 238   fTPCInner(0),
 239   fOp(0),
 240   fFriendTrack(0),
 241   fTPCClusterMap(159),//number of padrows
 242   fTPCSharedMap(159),//number of padrows
 243   fFlags(0),
 244   fID(0),
 245   fLabel(0),
 246   fITSLabel(0),
 247   fTPCLabel(0),
 248   fTRDLabel(0),
 249   fTOFCalChannel(0),
 250   fTOFindex(0),
 251   fHMPIDqn(0),
 252   fHMPIDcluIdx(0),
 253   fEMCALindex(kEMCALNoMatch),
 254   fHMPIDtrkTheta(0),
 255   fHMPIDtrkPhi(0),
 256   fHMPIDsignal(0),
 257   fTrackLength(0),
 258   fD(0),fZ(0),
 259   fCdd(0),fCdz(0),fCzz(0),
 260   fCchi2(0),
 261   fITSchi2(0),
 262   fTPCchi2(0),
 263   fTRDchi2(0),
 264   fTOFchi2(0),
 265   fHMPIDchi2(0),
 266   fITSsignal(0),
 267   fTPCsignal(0),
 268   fTPCsignalS(0),
 269   fTRDsignal(0),
 270   fTRDQuality(0),
 271   fTRDBudget(0),
 272   fTOFsignal(0),
 273   fTOFsignalToT(0),
 274   fTOFsignalRaw(0),
 275   fTOFsignalDz(0),
 276   fHMPIDtrkX(0),
 277   fHMPIDtrkY(0),
 278   fHMPIDmipX(0),
 279   fHMPIDmipY(0),
 280   fTPCncls(0),
 281   fTPCnclsF(0),
 282   fTPCsignalN(0),
 283   fITSncls(0),
 284   fITSClusterMap(0),
 285   fTRDncls(0),
 286   fTRDncls0(0),
 287   fTRDpidQuality(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]=0;
 312   }
 313   for (i=0;i<4;i++) {fTPCPoints[i]=0;}
 314   for (i=0;i<3;i++) {fTOFLabel[i]=0;}
 315   for (i=0;i<10;i++) {fTOFInfo[i]=0;}
 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
 449   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTrackTime[i] = 0;
 450
 451   // Reset track parameters constrained to the primary vertex
 452   delete fCp;fCp = 0;
 453   fCchi2 = 0;
 454
 455   // Reset track parameters at the inner wall of TPC
 456   delete fIp;fIp = 0;
 457   delete fTPCInner;fTPCInner=0;
 458   // Reset track parameters at the inner wall of the TRD
 459   delete fOp;fOp = 0;
 460
 461
 462   // Reset ITS track related information
 463   fITSchi2 = 0;
 464   fITSncls = 0;
 465   fITSClusterMap=0;
 466   fITSsignal = 0;
 467   for (Int_t i=0;i<AliPID::kSPECIES;i++) fITSr[i]=0;
 468   fITSLabel = 0;
 469
 470   // Reset TPC related track information
 471   fTPCchi2 = 0;
 472   fTPCncls = 0;
 473   fTPCnclsF = 0;
 474   fTPCClusterMap = 0;
 475   fTPCSharedMap = 0;
 476   fTPCsignal= 0;
 477   fTPCsignalS= 0;
 478   fTPCsignalN= 0;
 479   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTPCr[i]=0;
 480   fTPCLabel=0;
 481   for (Int_t i=0;i<4;i++) fTPCPoints[i] = 0;
 482   for (Int_t i=0; i<3;i++)   fKinkIndexes[i] = 0;
 483   for (Int_t i=0; i<3;i++)   fV0Indexes[i] = 0;
 484
 485   // Reset TRD related track information
 486   fTRDchi2 = 0;
 487   fTRDncls = 0;
 488   fTRDncls0 = 0;
 489   fTRDsignal = 0;
 490   for (Int_t i=0;i<kNPlane;i++) {
 491     for (Int_t j=0;j<kNSlice;j++) {
 492       fTRDsignals[i][j] = 0;
 493     }
 494     fTRDTimBin[i]  = 0;
 495   }
 496   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i] = 0;
 497   fTRDLabel = 0;
 498   fTRDQuality  = 0;
 499   fTRDpidQuality = 0;
 500   fTRDBudget  = 0;
 501
 502   // Reset TOF related track information
 503   fTOFchi2 = 0;
 504   fTOFindex = 0;
 505   fTOFsignal = 0;
 506   fTOFCalChannel = 0;
 507   fTOFsignalToT = 0;
 508   fTOFsignalRaw = 0;
 509   fTOFsignalDz = 0;
 510   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i] = 0;
 511   for (Int_t i=0;i<3;i++) fTOFLabel[i] = 0;
 512   for (Int_t i=0;i<10;i++) fTOFInfo[i] = 0;
 513
 514   // Reset HMPID related track information
 515   fHMPIDchi2 = 0;
 516   fHMPIDqn = 0;
 517   fHMPIDcluIdx = 0;
 518   fHMPIDsignal = 0;
 519   for (Int_t i=0;i<AliPID::kSPECIES;i++) fHMPIDr[i] = 0;
 520   fHMPIDtrkTheta = 0;
 521   fHMPIDtrkPhi = 0;
 522   fHMPIDtrkX = 0;
 523   fHMPIDtrkY = 0;
 524   fHMPIDmipX = 0;
 525   fHMPIDmipY = 0;
 526   fEMCALindex = kEMCALNoMatch;
 527
 528   delete fFriendTrack; fFriendTrack = 0;
 529 }
 530 //_______________________________________________________________________
 531 Double_t AliESDtrack::GetMass() const {
 532   // Returns the mass of the most probable particle type
 533   Float_t max=0.;
 534   Int_t k=-1;
 535   for (Int_t i=0; i<AliPID::kSPECIES; i++) {
 536     if (fR[i]>max) {k=i; max=fR[i];}
 537   }
 538   if (k==0) { // dE/dx "crossing points" in the TPC
 539      Double_t p=GetP();
 540      if ((p>0.38)&&(p<0.48))
 541         if (fR[0]<fR[3]*10.) return AliPID::ParticleMass(AliPID::kKaon);
 542      if ((p>0.75)&&(p<0.85))
 543         if (fR[0]<fR[4]*10.) return AliPID::ParticleMass(AliPID::kProton);
 544      return 0.00051;
 545   }
 546   if (k==1) return AliPID::ParticleMass(AliPID::kMuon);
 547   if (k==2||k==-1) return AliPID::ParticleMass(AliPID::kPion);
 548   if (k==3) return AliPID::ParticleMass(AliPID::kKaon);
 549   if (k==4) return AliPID::ParticleMass(AliPID::kProton);
 550   AliWarning("Undefined mass !");
 551   return AliPID::ParticleMass(AliPID::kPion);
 552 }
 553
 554 //______________________________________________________________________________
 555 Double_t AliESDtrack::E() const
 556 {
 557   // Returns the energy of the particle given its assumed mass.
 558   // Assumes the pion mass if the particle can't be identified properly.
 559
 560   Double_t m = M();
 561   Double_t p = P();
 562   return TMath::Sqrt(p*p + m*m);
 563 }
 564
 565 //______________________________________________________________________________
 566 Double_t AliESDtrack::Y() const
 567 {
 568   // Returns the rapidity of a particle given its assumed mass.
 569   // Assumes the pion mass if the particle can't be identified properly.
 570
 571   Double_t e = E();
 572   Double_t pz = Pz();
 573   if (e != TMath::Abs(pz)) { // energy was not equal to pz
 574     return 0.5*TMath::Log((e+pz)/(e-pz));
 575   } else { // energy was equal to pz
 576     return -999.;
 577   }
 578 }
 579
 580 //_______________________________________________________________________
 581 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
 582   //
 583   // This function updates track's running parameters
 584   //
 585   Int_t *index=0;
 586   Bool_t rc=kTRUE;
 587
 588   SetStatus(flags);
 589   fLabel=t->GetLabel();
 590
 591   if (t->IsStartedTimeIntegral()) {
 592     SetStatus(kTIME);
 593     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
 594     SetIntegratedLength(t->GetIntegratedLength());
 595   }
 596
 597   Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 598
 599   switch (flags) {
 600
 601   case kITSin: case kITSout: case kITSrefit:
 602     fITSClusterMap=0;
 603     fITSncls=t->GetNumberOfClusters();
 604     index=fFriendTrack->GetITSindices();
 605     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
 606         index[i]=t->GetClusterIndex(i);
 607         if (i<fITSncls) {
 608            Int_t l=(index[i] & 0xf0000000) >> 28;
 609            SETBIT(fITSClusterMap,l);
 610         }
 611     }
 612     fITSchi2=t->GetChi2();
 613     fITSsignal=t->GetPIDsignal();
 614     fITSLabel = t->GetLabel();
 615     break;
 616
 617   case kTPCin: case kTPCrefit:
 618     fTPCLabel = t->GetLabel();
 619     if (flags==kTPCin)  fTPCInner=new AliExternalTrackParam(*t);
 620     if (!fIp) fIp=new AliExternalTrackParam(*t);
 621     else
 622       fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 623   case kTPCout:
 624     index=fFriendTrack->GetTPCindices();
 625     if (flags & kTPCout){
 626       if (!fOp) fOp=new AliExternalTrackParam(*t);
 627       else
 628         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 629     }
 630     fTPCncls=t->GetNumberOfClusters();
 631     fTPCchi2=t->GetChi2();
 632
 633      {//prevrow must be declared in separate namespace, otherwise compiler cries:
 634       //"jump to case label crosses initialization of `Int_t prevrow'"
 635        Int_t prevrow = -1;
 636        //       for (Int_t i=0;i<fTPCncls;i++)
 637        for (Int_t i=0;i<AliESDfriendTrack::kMaxTPCcluster;i++)
 638         {
 639           index[i]=t->GetClusterIndex(i);
 640           Int_t idx = index[i];
 641
 642           if (idx<0) continue;
 643
 644           // Piotr's Cluster Map for HBT
 645           // ### please change accordingly if cluster array is changing
 646           // to "New TPC Tracking" style (with gaps in array)
 647           Int_t sect = (idx&0xff000000)>>24;
 648           Int_t row = (idx&0x00ff0000)>>16;
 649           if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
 650
 651           fTPCClusterMap.SetBitNumber(row,kTRUE);
 652
 653           //Fill the gap between previous row and this row with 0 bits
 654           //In case  ###  pleas change it as well - just set bit 0 in case there
 655           //is no associated clusters for current "i"
 656           if (prevrow < 0)
 657            {
 658              prevrow = row;//if previous bit was not assigned yet == this is the first one
 659            }
 660           else
 661            { //we don't know the order (inner to outer or reverse)
 662              //just to be save in case it is going to change
 663              Int_t n = 0, m = 0;
 664              if (prevrow < row)
 665               {
 666                 n = prevrow;
 667                 m = row;
 668               }
 669              else
 670               {
 671                 n = row;
 672                 m = prevrow;
 673               }
 674
 675              for (Int_t j = n+1; j < m; j++)
 676               {
 677                 fTPCClusterMap.SetBitNumber(j,kFALSE);
 678               }
 679              prevrow = row;
 680            }
 681           // End Of Piotr's Cluster Map for HBT
 682         }
 683      }
 684     fTPCsignal=t->GetPIDsignal();
 685     break;
 686
 687   case kTRDout: case kTRDin: case kTRDrefit:
 688     index=fFriendTrack->GetTRDindices();
 689     fTRDLabel = t->GetLabel();
 690     fTRDncls=t->GetNumberOfClusters();
 691     fTRDchi2=t->GetChi2();
 692     for (Int_t i=0;i<fTRDncls;i++) index[i]=t->GetClusterIndex(i);
 693     fTRDsignal=t->GetPIDsignal();
 694     break;
 695   case kTRDbackup:
 696     if (!fOp) fOp=new AliExternalTrackParam(*t);
 697     else
 698       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 699     fTRDncls0 = t->GetNumberOfClusters();
 700     break;
 701   case kTOFin:
 702     break;
 703   case kTOFout:
 704     break;
 705   case kTRDStop:
 706     break;
 707   default:
 708     AliError("Wrong flag !");
 709     return kFALSE;
 710   }
 711
 712   return rc;
 713 }
 714
 715 //_______________________________________________________________________
 716 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
 717   //---------------------------------------------------------------------
 718   // This function returns external representation of the track parameters
 719   //---------------------------------------------------------------------
 720   x=GetX();
 721   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
 722 }
 723
 724 //_______________________________________________________________________
 725 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
 726   //---------------------------------------------------------------------
 727   // This function returns external representation of the cov. matrix
 728   //---------------------------------------------------------------------
 729   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
 730 }
 731
 732 //_______________________________________________________________________
 733 Bool_t AliESDtrack::GetConstrainedExternalParameters
 734                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 735   //---------------------------------------------------------------------
 736   // This function returns the constrained external track parameters
 737   //---------------------------------------------------------------------
 738   if (!fCp) return kFALSE;
 739   alpha=fCp->GetAlpha();
 740   x=fCp->GetX();
 741   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
 742   return kTRUE;
 743 }
 744
 745 //_______________________________________________________________________
 746 Bool_t
 747 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
 748   //---------------------------------------------------------------------
 749   // This function returns the constrained external cov. matrix
 750   //---------------------------------------------------------------------
 751   if (!fCp) return kFALSE;
 752   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
 753   return kTRUE;
 754 }
 755
 756 Bool_t
 757 AliESDtrack::GetInnerExternalParameters
 758                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 759   //---------------------------------------------------------------------
 760   // This function returns external representation of the track parameters
 761   // at the inner layer of TPC
 762   //---------------------------------------------------------------------
 763   if (!fIp) return kFALSE;
 764   alpha=fIp->GetAlpha();
 765   x=fIp->GetX();
 766   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
 767   return kTRUE;
 768 }
 769
 770 Bool_t
 771 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
 772  //---------------------------------------------------------------------
 773  // This function returns external representation of the cov. matrix
 774  // at the inner layer of TPC
 775  //---------------------------------------------------------------------
 776   if (!fIp) return kFALSE;
 777   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
 778   return kTRUE;
 779 }
 780
 781 Bool_t
 782 AliESDtrack::GetOuterExternalParameters
 783                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 784   //---------------------------------------------------------------------
 785   // This function returns external representation of the track parameters
 786   // at the inner layer of TRD
 787   //---------------------------------------------------------------------
 788   if (!fOp) return kFALSE;
 789   alpha=fOp->GetAlpha();
 790   x=fOp->GetX();
 791   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
 792   return kTRUE;
 793 }
 794
 795 Bool_t
 796 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
 797  //---------------------------------------------------------------------
 798  // This function returns external representation of the cov. matrix
 799  // at the inner layer of TRD
 800  //---------------------------------------------------------------------
 801   if (!fOp) return kFALSE;
 802   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
 803   return kTRUE;
 804 }
 805
 806 Int_t AliESDtrack::GetNcls(Int_t idet) const
 807 {
 808   // Get number of clusters by subdetector index
 809   //
 810   Int_t ncls = 0;
 811   switch(idet){
 812   case 0:
 813     ncls = fITSncls;
 814     break;
 815   case 1:
 816     ncls = fTPCncls;
 817     break;
 818   case 2:
 819     ncls = fTRDncls;
 820     break;
 821   case 3:
 822     if (fTOFindex != 0)
 823       ncls = 1;
 824     break;
 825   default:
 826     break;
 827   }
 828   return ncls;
 829 }
 830
 831 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
 832 {
 833   // Get cluster index array by subdetector index
 834   //
 835   Int_t ncls = 0;
 836   switch(idet){
 837   case 0:
 838     ncls = GetITSclusters(idx);
 839     break;
 840   case 1:
 841     ncls = GetTPCclusters(idx);
 842     break;
 843   case 2:
 844     ncls = GetTRDclusters(idx);
 845     break;
 846   case 3:
 847     if (fTOFindex != 0) {
 848       idx[0] = GetTOFcluster();
 849       ncls = 1;
 850     }
 851     break;
 852   default:
 853     break;
 854   }
 855   return ncls;
 856 }
 857
 858 //_______________________________________________________________________
 859 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
 860   // Returns the array with integrated times for each particle hypothesis
 861   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
 862 }
 863
 864 //_______________________________________________________________________
 865 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
 866   // Sets the array with integrated times for each particle hypotesis
 867   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
 868 }
 869
 870 //_______________________________________________________________________
 871 void AliESDtrack::SetITSpid(const Double_t *p) {
 872   // Sets values for the probability of each particle type (in ITS)
 873   SetPIDValues(fITSr,p,AliPID::kSPECIES);
 874   SetStatus(AliESDtrack::kITSpid);
 875 }
 876
 877 //_______________________________________________________________________
 878 void AliESDtrack::GetITSpid(Double_t *p) const {
 879   // Gets the probability of each particle type (in ITS)
 880   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
 881 }
 882
 883 //_______________________________________________________________________
 884 Char_t AliESDtrack::GetITSclusters(Int_t *idx) const {
 885   //---------------------------------------------------------------------
 886   // This function returns indices of the assgined ITS clusters
 887   //---------------------------------------------------------------------
 888   if (idx!=0) {
 889      Int_t *index=fFriendTrack->GetITSindices();
 890      for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) idx[i]=index[i];
 891   }
 892   return fITSncls;
 893 }
 894
 895 //_______________________________________________________________________
 896 UShort_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
 897   //---------------------------------------------------------------------
 898   // This function returns indices of the assgined ITS clusters
 899   //---------------------------------------------------------------------
 900   if (idx!=0) {
 901     Int_t *index=fFriendTrack->GetTPCindices();
 902     for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
 903   }
 904   return fTPCncls;
 905 }
 906
 907 Double_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
 908   //
 909   // GetDensity of the clusters on given region between row0 and row1
 910   // Dead zone effect takin into acoount
 911   //
 912   Int_t good  = 0;
 913   Int_t found = 0;
 914   //
 915   Int_t *index=fFriendTrack->GetTPCindices();
 916   for (Int_t i=row0;i<=row1;i++){
 917     Int_t idx = index[i];
 918     if (idx!=-1)  good++;             // track outside of dead zone
 919     if (idx>0)    found++;
 920   }
 921   Float_t density=0.5;
 922   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
 923   return density;
 924 }
 925
 926 //_______________________________________________________________________
 927 void AliESDtrack::SetTPCpid(const Double_t *p) {
 928   // Sets values for the probability of each particle type (in TPC)
 929   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
 930   SetStatus(AliESDtrack::kTPCpid);
 931 }
 932
 933 //_______________________________________________________________________
 934 void AliESDtrack::GetTPCpid(Double_t *p) const {
 935   // Gets the probability of each particle type (in TPC)
 936   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
 937 }
 938
 939 //_______________________________________________________________________
 940 UChar_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
 941   //---------------------------------------------------------------------
 942   // This function returns indices of the assgined TRD clusters
 943   //---------------------------------------------------------------------
 944   if (idx!=0) {
 945      Int_t *index=fFriendTrack->GetTRDindices();
 946      for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
 947   }
 948   return fTRDncls;
 949 }
 950
 951 //_______________________________________________________________________
 952 void AliESDtrack::SetTRDpid(const Double_t *p) {
 953   // Sets values for the probability of each particle type (in TRD)
 954   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
 955   SetStatus(AliESDtrack::kTRDpid);
 956 }
 957
 958 //_______________________________________________________________________
 959 void AliESDtrack::GetTRDpid(Double_t *p) const {
 960   // Gets the probability of each particle type (in TRD)
 961   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
 962 }
 963
 964 //_______________________________________________________________________
 965 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
 966 {
 967   // Sets the probability of particle type iSpecies to p (in TRD)
 968   fTRDr[iSpecies] = p;
 969 }
 970
 971 Double_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
 972 {
 973   // Returns the probability of particle type iSpecies (in TRD)
 974   return fTRDr[iSpecies];
 975 }
 976
 977 //_______________________________________________________________________
 978 void AliESDtrack::SetTOFpid(const Double_t *p) {
 979   // Sets the probability of each particle type (in TOF)
 980   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
 981   SetStatus(AliESDtrack::kTOFpid);
 982 }
 983
 984 //_______________________________________________________________________
 985 void AliESDtrack::SetTOFLabel(const Int_t *p) {
 986   // Sets  (in TOF)
 987   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
 988 }
 989
 990 //_______________________________________________________________________
 991 void AliESDtrack::GetTOFpid(Double_t *p) const {
 992   // Gets probabilities of each particle type (in TOF)
 993   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
 994 }
 995
 996 //_______________________________________________________________________
 997 void AliESDtrack::GetTOFLabel(Int_t *p) const {
 998   // Gets (in TOF)
 999   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
1000 }
1001
1002 //_______________________________________________________________________
1003 void AliESDtrack::GetTOFInfo(Float_t *info) const {
1004   // Gets (in TOF)
1005   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
1006 }
1007
1008 //_______________________________________________________________________
1009 void AliESDtrack::SetTOFInfo(Float_t*info) {
1010   // Gets (in TOF)
1011   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
1012 }
1013
1014
1015
1016 //_______________________________________________________________________
1017 void AliESDtrack::SetHMPIDpid(const Double_t *p) {
1018   // Sets the probability of each particle type (in HMPID)
1019   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
1020   SetStatus(AliESDtrack::kHMPIDpid);
1021 }
1022
1023 //_______________________________________________________________________
1024 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
1025   // Gets probabilities of each particle type (in HMPID)
1026   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
1027 }
1028
1029
1030
1031 //_______________________________________________________________________
1032 void AliESDtrack::SetESDpid(const Double_t *p) {
1033   // Sets the probability of each particle type for the ESD track
1034   SetPIDValues(fR,p,AliPID::kSPECIES);
1035   SetStatus(AliESDtrack::kESDpid);
1036 }
1037
1038 //_______________________________________________________________________
1039 void AliESDtrack::GetESDpid(Double_t *p) const {
1040   // Gets probability of each particle type for the ESD track
1041   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
1042 }
1043
1044 //_______________________________________________________________________
1045 Bool_t AliESDtrack::RelateToVertex
1046 (const AliESDVertex *vtx, Double_t b, Double_t maxd) {
1047   //
1048   // Try to relate this track to the vertex "vtx",
1049   // if the (rough) transverse impact parameter is not bigger then "maxd".
1050   //            Magnetic field is "b" (kG).
1051   //
1052   // a) The track gets extapolated to the DCA to the vertex.
1053   // b) The impact parameters and their covariance matrix are calculated.
1054   // c) An attempt to constrain this track to the vertex is done.
1055   //
1056   //    In the case of success, the returned value is kTRUE
1057   //    (otherwise, it's kFALSE)
1058   //
1059
1060   if (!vtx) return kFALSE;
1061
1062   Double_t alpha=GetAlpha();
1063   Double_t sn=TMath::Sin(alpha), cs=TMath::Cos(alpha);
1064   Double_t x=GetX(), y=GetParameter()[0], snp=GetParameter()[2];
1065   Double_t xv= vtx->GetXv()*cs + vtx->GetYv()*sn;
1066   Double_t yv=-vtx->GetXv()*sn + vtx->GetYv()*cs, zv=vtx->GetZv();
1067   x-=xv; y-=yv;
1068
1069   //Estimate the impact parameter neglecting the track curvature
1070   Double_t d=TMath::Abs(x*snp - y*TMath::Sqrt(1.- snp*snp));
1071   if (d > maxd) return kFALSE;
1072
1073   //Propagate to the DCA
1074   Double_t crv=kB2C*b*GetParameter()[4];
1075   if (TMath::Abs(b) < kAlmost0Field) crv=0.;
1076
1077   Double_t tgfv=-(crv*x - snp)/(crv*y + TMath::Sqrt(1.-snp*snp));
1078   sn=tgfv/TMath::Sqrt(1.+ tgfv*tgfv);
1079   if (TMath::Abs(tgfv)>0.) cs = sn/tgfv;
1080   else cs=1.;
1081
1082   x = xv*cs + yv*sn;
1083   yv=-xv*sn + yv*cs; xv=x;
1084
1085   if (!Propagate(alpha+TMath::ASin(sn),xv,b)) return kFALSE;
1086
1087   fD = GetParameter()[0] - yv;
1088   fZ = GetParameter()[1] - zv;
1089
1090   Double_t cov[6]; vtx->GetCovMatrix(cov);
1091
1092   //***** Improvements by A.Dainese
1093   alpha=GetAlpha(); sn=TMath::Sin(alpha); cs=TMath::Cos(alpha);
1094   Double_t s2ylocvtx = cov[0]*sn*sn + cov[2]*cs*cs - 2.*cov[1]*cs*sn;
1095   fCdd = GetCovariance()[0] + s2ylocvtx;   // neglecting correlations
1096   fCdz = GetCovariance()[1];               // between (x,y) and z
1097   fCzz = GetCovariance()[2] + cov[5];      // in vertex's covariance matrix
1098   //*****
1099
1100   {//Try to constrain
1101     Double_t p[2]={yv,zv}, c[3]={cov[2],0.,cov[5]};
1102     Double_t chi2=GetPredictedChi2(p,c);
1103
1104     if (chi2>77.) return kFALSE;
1105
1106     AliExternalTrackParam tmp(*this);
1107     if (!tmp.Update(p,c)) return kFALSE;
1108
1109     fCchi2=chi2;
1110     if (!fCp) fCp=new AliExternalTrackParam();
1111     new (fCp) AliExternalTrackParam(tmp);
1112   }
1113
1114   return kTRUE;
1115 }
1116
1117 //_______________________________________________________________________
1118 void AliESDtrack::Print(Option_t *) const {
1119   // Prints info on the track
1120
1121   printf("ESD track info\n") ;
1122   Double_t p[AliPID::kSPECIESN] ;
1123   Int_t index = 0 ;
1124   if( IsOn(kITSpid) ){
1125     printf("From ITS: ") ;
1126     GetITSpid(p) ;
1127     for(index = 0 ; index < AliPID::kSPECIES; index++)
1128       printf("%f, ", p[index]) ;
1129     printf("\n           signal = %f\n", GetITSsignal()) ;
1130   }
1131   if( IsOn(kTPCpid) ){
1132     printf("From TPC: ") ;
1133     GetTPCpid(p) ;
1134     for(index = 0 ; index < AliPID::kSPECIES; index++)
1135       printf("%f, ", p[index]) ;
1136     printf("\n           signal = %f\n", GetTPCsignal()) ;
1137   }
1138   if( IsOn(kTRDpid) ){
1139     printf("From TRD: ") ;
1140     GetTRDpid(p) ;
1141     for(index = 0 ; index < AliPID::kSPECIES; index++)
1142       printf("%f, ", p[index]) ;
1143     printf("\n           signal = %f\n", GetTRDsignal()) ;
1144   }
1145   if( IsOn(kTOFpid) ){
1146     printf("From TOF: ") ;
1147     GetTOFpid(p) ;
1148     for(index = 0 ; index < AliPID::kSPECIES; index++)
1149       printf("%f, ", p[index]) ;
1150     printf("\n           signal = %f\n", GetTOFsignal()) ;
1151   }
1152   if( IsOn(kHMPIDpid) ){
1153     printf("From HMPID: ") ;
1154     GetHMPIDpid(p) ;
1155     for(index = 0 ; index < AliPID::kSPECIES; index++)
1156       printf("%f, ", p[index]) ;
1157     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
1158   }
1159 }
1160