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 Double_t AliESDtrack::E() const
 554 {
 555   // Returns the energy of the particle given its assumed mass.
 556   // Assumes the pion mass if the particle can't be identified properly.
 557
 558   Double_t m = M();
 559   Double_t p = P();
 560   return TMath::Sqrt(p*p + m*m);
 561 }
 562
 563 //______________________________________________________________________________
 564 Double_t AliESDtrack::Y() const
 565 {
 566   // Returns the rapidity of a particle given its assumed mass.
 567   // Assumes the pion mass if the particle can't be identified properly.
 568
 569   Double_t e = E();
 570   Double_t pz = Pz();
 571   if (e!=pz) { // energy was not equal to pz
 572     return 0.5*TMath::Log((e+pz)/(e-pz));
 573   } else { // energy was equal to pz
 574     return -999.;
 575   }
 576 }
 577
 578 //_______________________________________________________________________
 579 Bool_t AliESDtrack::UpdateTrackParams(const AliKalmanTrack *t, ULong_t flags){
 580   //
 581   // This function updates track's running parameters
 582   //
 583   Int_t *index=0;
 584   Bool_t rc=kTRUE;
 585
 586   SetStatus(flags);
 587   fLabel=t->GetLabel();
 588
 589   if (t->IsStartedTimeIntegral()) {
 590     SetStatus(kTIME);
 591     Double_t times[10];t->GetIntegratedTimes(times); SetIntegratedTimes(times);
 592     SetIntegratedLength(t->GetIntegratedLength());
 593   }
 594
 595   Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 596
 597   switch (flags) {
 598
 599   case kITSin: case kITSout: case kITSrefit:
 600     fITSncls=t->GetNumberOfClusters();
 601     index=fFriendTrack->GetITSindices();
 602     for (Int_t i=0;i<AliESDfriendTrack::kMaxITScluster;i++) {
 603         index[i]=t->GetClusterIndex(i);
 604         if (i<fITSncls) {
 605            Int_t l=(index[i] & 0xf0000000) >> 28;
 606            SETBIT(fITSClusterMap,l);
 607         }
 608     }
 609     fITSchi2=t->GetChi2();
 610     fITSsignal=t->GetPIDsignal();
 611     fITSLabel = t->GetLabel();
 612     break;
 613
 614   case kTPCin: case kTPCrefit:
 615     fTPCLabel = t->GetLabel();
 616     if (flags==kTPCin)  fTPCInner=new AliExternalTrackParam(*t);
 617     if (!fIp) fIp=new AliExternalTrackParam(*t);
 618     else
 619       fIp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 620   case kTPCout:
 621     index=fFriendTrack->GetTPCindices();
 622     if (flags & kTPCout){
 623       if (!fOp) fOp=new AliExternalTrackParam(*t);
 624       else
 625         fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 626     }
 627     fTPCncls=t->GetNumberOfClusters();
 628     fTPCchi2=t->GetChi2();
 629
 630      {//prevrow must be declared in separate namespace, otherwise compiler cries:
 631       //"jump to case label crosses initialization of `Int_t prevrow'"
 632        Int_t prevrow = -1;
 633        //       for (Int_t i=0;i<fTPCncls;i++)
 634        for (Int_t i=0;i<AliESDfriendTrack::kMaxTPCcluster;i++)
 635         {
 636           index[i]=t->GetClusterIndex(i);
 637           Int_t idx = index[i];
 638
 639           if (idx<0) continue;
 640
 641           // Piotr's Cluster Map for HBT
 642           // ### please change accordingly if cluster array is changing
 643           // to "New TPC Tracking" style (with gaps in array)
 644           Int_t sect = (idx&0xff000000)>>24;
 645           Int_t row = (idx&0x00ff0000)>>16;
 646           if (sect > 18) row +=63; //if it is outer sector, add number of inner sectors
 647
 648           fTPCClusterMap.SetBitNumber(row,kTRUE);
 649
 650           //Fill the gap between previous row and this row with 0 bits
 651           //In case  ###  pleas change it as well - just set bit 0 in case there
 652           //is no associated clusters for current "i"
 653           if (prevrow < 0)
 654            {
 655              prevrow = row;//if previous bit was not assigned yet == this is the first one
 656            }
 657           else
 658            { //we don't know the order (inner to outer or reverse)
 659              //just to be save in case it is going to change
 660              Int_t n = 0, m = 0;
 661              if (prevrow < row)
 662               {
 663                 n = prevrow;
 664                 m = row;
 665               }
 666              else
 667               {
 668                 n = row;
 669                 m = prevrow;
 670               }
 671
 672              for (Int_t j = n+1; j < m; j++)
 673               {
 674                 fTPCClusterMap.SetBitNumber(j,kFALSE);
 675               }
 676              prevrow = row;
 677            }
 678           // End Of Piotr's Cluster Map for HBT
 679         }
 680      }
 681     fTPCsignal=t->GetPIDsignal();
 682     break;
 683
 684   case kTRDout: case kTRDin: case kTRDrefit:
 685     index=fFriendTrack->GetTRDindices();
 686     fTRDLabel = t->GetLabel();
 687     fTRDncls=t->GetNumberOfClusters();
 688     fTRDchi2=t->GetChi2();
 689     for (Int_t i=0;i<fTRDncls;i++) index[i]=t->GetClusterIndex(i);
 690     fTRDsignal=t->GetPIDsignal();
 691     break;
 692   case kTRDbackup:
 693     if (!fOp) fOp=new AliExternalTrackParam(*t);
 694     else
 695       fOp->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
 696     fTRDncls0 = t->GetNumberOfClusters();
 697     break;
 698   case kTOFin:
 699     break;
 700   case kTOFout:
 701     break;
 702   case kTRDStop:
 703     break;
 704   default:
 705     AliError("Wrong flag !");
 706     return kFALSE;
 707   }
 708
 709   return rc;
 710 }
 711
 712 //_______________________________________________________________________
 713 void AliESDtrack::GetExternalParameters(Double_t &x, Double_t p[5]) const {
 714   //---------------------------------------------------------------------
 715   // This function returns external representation of the track parameters
 716   //---------------------------------------------------------------------
 717   x=GetX();
 718   for (Int_t i=0; i<5; i++) p[i]=GetParameter()[i];
 719 }
 720
 721 //_______________________________________________________________________
 722 void AliESDtrack::GetExternalCovariance(Double_t cov[15]) const {
 723   //---------------------------------------------------------------------
 724   // This function returns external representation of the cov. matrix
 725   //---------------------------------------------------------------------
 726   for (Int_t i=0; i<15; i++) cov[i]=AliExternalTrackParam::GetCovariance()[i];
 727 }
 728
 729 //_______________________________________________________________________
 730 Bool_t AliESDtrack::GetConstrainedExternalParameters
 731                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 732   //---------------------------------------------------------------------
 733   // This function returns the constrained external track parameters
 734   //---------------------------------------------------------------------
 735   if (!fCp) return kFALSE;
 736   alpha=fCp->GetAlpha();
 737   x=fCp->GetX();
 738   for (Int_t i=0; i<5; i++) p[i]=fCp->GetParameter()[i];
 739   return kTRUE;
 740 }
 741
 742 //_______________________________________________________________________
 743 Bool_t
 744 AliESDtrack::GetConstrainedExternalCovariance(Double_t c[15]) const {
 745   //---------------------------------------------------------------------
 746   // This function returns the constrained external cov. matrix
 747   //---------------------------------------------------------------------
 748   if (!fCp) return kFALSE;
 749   for (Int_t i=0; i<15; i++) c[i]=fCp->GetCovariance()[i];
 750   return kTRUE;
 751 }
 752
 753 Bool_t
 754 AliESDtrack::GetInnerExternalParameters
 755                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 756   //---------------------------------------------------------------------
 757   // This function returns external representation of the track parameters
 758   // at the inner layer of TPC
 759   //---------------------------------------------------------------------
 760   if (!fIp) return kFALSE;
 761   alpha=fIp->GetAlpha();
 762   x=fIp->GetX();
 763   for (Int_t i=0; i<5; i++) p[i]=fIp->GetParameter()[i];
 764   return kTRUE;
 765 }
 766
 767 Bool_t
 768 AliESDtrack::GetInnerExternalCovariance(Double_t cov[15]) const {
 769  //---------------------------------------------------------------------
 770  // This function returns external representation of the cov. matrix
 771  // at the inner layer of TPC
 772  //---------------------------------------------------------------------
 773   if (!fIp) return kFALSE;
 774   for (Int_t i=0; i<15; i++) cov[i]=fIp->GetCovariance()[i];
 775   return kTRUE;
 776 }
 777
 778 Bool_t
 779 AliESDtrack::GetOuterExternalParameters
 780                  (Double_t &alpha, Double_t &x, Double_t p[5]) const {
 781   //---------------------------------------------------------------------
 782   // This function returns external representation of the track parameters
 783   // at the inner layer of TRD
 784   //---------------------------------------------------------------------
 785   if (!fOp) return kFALSE;
 786   alpha=fOp->GetAlpha();
 787   x=fOp->GetX();
 788   for (Int_t i=0; i<5; i++) p[i]=fOp->GetParameter()[i];
 789   return kTRUE;
 790 }
 791
 792 Bool_t
 793 AliESDtrack::GetOuterExternalCovariance(Double_t cov[15]) const {
 794  //---------------------------------------------------------------------
 795  // This function returns external representation of the cov. matrix
 796  // at the inner layer of TRD
 797  //---------------------------------------------------------------------
 798   if (!fOp) return kFALSE;
 799   for (Int_t i=0; i<15; i++) cov[i]=fOp->GetCovariance()[i];
 800   return kTRUE;
 801 }
 802
 803 Int_t AliESDtrack::GetNcls(Int_t idet) const
 804 {
 805   // Get number of clusters by subdetector index
 806   //
 807   Int_t ncls = 0;
 808   switch(idet){
 809   case 0:
 810     ncls = fITSncls;
 811     break;
 812   case 1:
 813     ncls = fTPCncls;
 814     break;
 815   case 2:
 816     ncls = fTRDncls;
 817     break;
 818   case 3:
 819     if (fTOFindex != 0)
 820       ncls = 1;
 821     break;
 822   default:
 823     break;
 824   }
 825   return ncls;
 826 }
 827
 828 Int_t AliESDtrack::GetClusters(Int_t idet, Int_t *idx) const
 829 {
 830   // Get cluster index array by subdetector index
 831   //
 832   Int_t ncls = 0;
 833   switch(idet){
 834   case 0:
 835     ncls = GetITSclusters(idx);
 836     break;
 837   case 1:
 838     ncls = GetTPCclusters(idx);
 839     break;
 840   case 2:
 841     ncls = GetTRDclusters(idx);
 842     break;
 843   case 3:
 844     if (fTOFindex != 0) {
 845       idx[0] = GetTOFcluster();
 846       ncls = 1;
 847     }
 848     break;
 849   default:
 850     break;
 851   }
 852   return ncls;
 853 }
 854
 855 //_______________________________________________________________________
 856 void AliESDtrack::GetIntegratedTimes(Double_t *times) const {
 857   // Returns the array with integrated times for each particle hypothesis
 858   for (Int_t i=0; i<AliPID::kSPECIES; i++) times[i]=fTrackTime[i];
 859 }
 860
 861 //_______________________________________________________________________
 862 void AliESDtrack::SetIntegratedTimes(const Double_t *times) {
 863   // Sets the array with integrated times for each particle hypotesis
 864   for (Int_t i=0; i<AliPID::kSPECIES; i++) fTrackTime[i]=times[i];
 865 }
 866
 867 //_______________________________________________________________________
 868 void AliESDtrack::SetITSpid(const Double_t *p) {
 869   // Sets values for the probability of each particle type (in ITS)
 870   SetPIDValues(fITSr,p,AliPID::kSPECIES);
 871   SetStatus(AliESDtrack::kITSpid);
 872 }
 873
 874 //_______________________________________________________________________
 875 void AliESDtrack::GetITSpid(Double_t *p) const {
 876   // Gets the probability of each particle type (in ITS)
 877   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fITSr[i];
 878 }
 879
 880 //_______________________________________________________________________
 881 Int_t AliESDtrack::GetITSclusters(Int_t *idx) const {
 882   //---------------------------------------------------------------------
 883   // This function returns indices of the assgined ITS clusters
 884   //---------------------------------------------------------------------
 885   if (idx!=0) {
 886      Int_t *index=fFriendTrack->GetITSindices();
 887      for (Int_t i=0; i<AliESDfriendTrack::kMaxITScluster; i++) idx[i]=index[i];
 888   }
 889   return fITSncls;
 890 }
 891
 892 //_______________________________________________________________________
 893 Int_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
 894   //---------------------------------------------------------------------
 895   // This function returns indices of the assgined ITS clusters
 896   //---------------------------------------------------------------------
 897   if (idx!=0) {
 898     Int_t *index=fFriendTrack->GetTPCindices();
 899     for (Int_t i=0; i<AliESDfriendTrack::kMaxTPCcluster; i++) idx[i]=index[i];
 900   }
 901   return fTPCncls;
 902 }
 903
 904 Float_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
 905   //
 906   // GetDensity of the clusters on given region between row0 and row1
 907   // Dead zone effect takin into acoount
 908   //
 909   Int_t good  = 0;
 910   Int_t found = 0;
 911   //
 912   Int_t *index=fFriendTrack->GetTPCindices();
 913   for (Int_t i=row0;i<=row1;i++){
 914     Int_t idx = index[i];
 915     if (idx!=-1)  good++;             // track outside of dead zone
 916     if (idx>0)    found++;
 917   }
 918   Float_t density=0.5;
 919   if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
 920   return density;
 921 }
 922
 923 //_______________________________________________________________________
 924 void AliESDtrack::SetTPCpid(const Double_t *p) {
 925   // Sets values for the probability of each particle type (in TPC)
 926   SetPIDValues(fTPCr,p,AliPID::kSPECIES);
 927   SetStatus(AliESDtrack::kTPCpid);
 928 }
 929
 930 //_______________________________________________________________________
 931 void AliESDtrack::GetTPCpid(Double_t *p) const {
 932   // Gets the probability of each particle type (in TPC)
 933   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTPCr[i];
 934 }
 935
 936 //_______________________________________________________________________
 937 Int_t AliESDtrack::GetTRDclusters(Int_t *idx) const {
 938   //---------------------------------------------------------------------
 939   // This function returns indices of the assgined TRD clusters
 940   //---------------------------------------------------------------------
 941   if (idx!=0) {
 942      Int_t *index=fFriendTrack->GetTRDindices();
 943      for (Int_t i=0; i<AliESDfriendTrack::kMaxTRDcluster; i++) idx[i]=index[i];
 944   }
 945   return fTRDncls;
 946 }
 947
 948 //_______________________________________________________________________
 949 void AliESDtrack::SetTRDpid(const Double_t *p) {
 950   // Sets values for the probability of each particle type (in TRD)
 951   SetPIDValues(fTRDr,p,AliPID::kSPECIES);
 952   SetStatus(AliESDtrack::kTRDpid);
 953 }
 954
 955 //_______________________________________________________________________
 956 void AliESDtrack::GetTRDpid(Double_t *p) const {
 957   // Gets the probability of each particle type (in TRD)
 958   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTRDr[i];
 959 }
 960
 961 //_______________________________________________________________________
 962 void    AliESDtrack::SetTRDpid(Int_t iSpecies, Float_t p)
 963 {
 964   // Sets the probability of particle type iSpecies to p (in TRD)
 965   fTRDr[iSpecies] = p;
 966 }
 967
 968 Float_t AliESDtrack::GetTRDpid(Int_t iSpecies) const
 969 {
 970   // Returns the probability of particle type iSpecies (in TRD)
 971   return fTRDr[iSpecies];
 972 }
 973
 974 //_______________________________________________________________________
 975 void AliESDtrack::SetTOFpid(const Double_t *p) {
 976   // Sets the probability of each particle type (in TOF)
 977   SetPIDValues(fTOFr,p,AliPID::kSPECIES);
 978   SetStatus(AliESDtrack::kTOFpid);
 979 }
 980
 981 //_______________________________________________________________________
 982 void AliESDtrack::SetTOFLabel(const Int_t *p) {
 983   // Sets  (in TOF)
 984   for (Int_t i=0; i<3; i++) fTOFLabel[i]=p[i];
 985 }
 986
 987 //_______________________________________________________________________
 988 void AliESDtrack::GetTOFpid(Double_t *p) const {
 989   // Gets probabilities of each particle type (in TOF)
 990   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fTOFr[i];
 991 }
 992
 993 //_______________________________________________________________________
 994 void AliESDtrack::GetTOFLabel(Int_t *p) const {
 995   // Gets (in TOF)
 996   for (Int_t i=0; i<3; i++) p[i]=fTOFLabel[i];
 997 }
 998
 999 //_______________________________________________________________________
1000 void AliESDtrack::GetTOFInfo(Float_t *info) const {
1001   // Gets (in TOF)
1002   for (Int_t i=0; i<10; i++) info[i]=fTOFInfo[i];
1003 }
1004
1005 //_______________________________________________________________________
1006 void AliESDtrack::SetTOFInfo(Float_t*info) {
1007   // Gets (in TOF)
1008   for (Int_t i=0; i<10; i++) fTOFInfo[i]=info[i];
1009 }
1010
1011
1012
1013 //_______________________________________________________________________
1014 void AliESDtrack::SetHMPIDpid(const Double_t *p) {
1015   // Sets the probability of each particle type (in HMPID)
1016   SetPIDValues(fHMPIDr,p,AliPID::kSPECIES);
1017   SetStatus(AliESDtrack::kHMPIDpid);
1018 }
1019
1020 //_______________________________________________________________________
1021 void AliESDtrack::GetHMPIDpid(Double_t *p) const {
1022   // Gets probabilities of each particle type (in HMPID)
1023   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fHMPIDr[i];
1024 }
1025
1026
1027
1028 //_______________________________________________________________________
1029 void AliESDtrack::SetESDpid(const Double_t *p) {
1030   // Sets the probability of each particle type for the ESD track
1031   SetPIDValues(fR,p,AliPID::kSPECIES);
1032   SetStatus(AliESDtrack::kESDpid);
1033 }
1034
1035 //_______________________________________________________________________
1036 void AliESDtrack::GetESDpid(Double_t *p) const {
1037   // Gets probability of each particle type for the ESD track
1038   for (Int_t i=0; i<AliPID::kSPECIES; i++) p[i]=fR[i];
1039 }
1040
1041 //_______________________________________________________________________
1042 Bool_t AliESDtrack::RelateToVertex
1043 (const AliESDVertex *vtx, Double_t b, Double_t maxd) {
1044   //
1045   // Try to relate this track to the vertex "vtx",
1046   // if the (rough) transverse impact parameter is not bigger then "maxd".
1047   //            Magnetic field is "b" (kG).
1048   //
1049   // a) The track gets extapolated to the DCA to the vertex.
1050   // b) The impact parameters and their covariance matrix are calculated.
1051   // c) An attempt to constrain this track to the vertex is done.
1052   //
1053   //    In the case of success, the returned value is kTRUE
1054   //    (otherwise, it's kFALSE)
1055   //
1056
1057   if (!vtx) return kFALSE;
1058
1059   Double_t alpha=GetAlpha();
1060   Double_t sn=TMath::Sin(alpha), cs=TMath::Cos(alpha);
1061   Double_t x=GetX(), y=GetParameter()[0], snp=GetParameter()[2];
1062   Double_t xv= vtx->GetXv()*cs + vtx->GetYv()*sn;
1063   Double_t yv=-vtx->GetXv()*sn + vtx->GetYv()*cs, zv=vtx->GetZv();
1064   x-=xv; y-=yv;
1065
1066   //Estimate the impact parameter neglecting the track curvature
1067   Double_t d=TMath::Abs(x*snp - y*TMath::Sqrt(1.- snp*snp));
1068   if (d > maxd) return kFALSE;
1069
1070   //Propagate to the DCA
1071   Double_t crv=kB2C*b*GetParameter()[4];
1072   if (TMath::Abs(b) < kAlmost0Field) crv=0.;
1073
1074   Double_t tgfv=-(crv*x - snp)/(crv*y + TMath::Sqrt(1.-snp*snp));
1075   sn=tgfv/TMath::Sqrt(1.+ tgfv*tgfv);
1076   if (TMath::Abs(tgfv)>0.) cs = sn/tgfv;
1077   else cs=1.;
1078
1079   x = xv*cs + yv*sn;
1080   yv=-xv*sn + yv*cs; xv=x;
1081
1082   if (!Propagate(alpha+TMath::ASin(sn),xv,b)) return kFALSE;
1083
1084   fD = GetParameter()[0] - yv;
1085   fZ = GetParameter()[1] - zv;
1086
1087   Double_t cov[6]; vtx->GetCovMatrix(cov);
1088
1089   //***** Improvements by A.Dainese
1090   alpha=GetAlpha(); sn=TMath::Sin(alpha); cs=TMath::Cos(alpha);
1091   Double_t s2ylocvtx = cov[0]*sn*sn + cov[2]*cs*cs - 2.*cov[1]*cs*sn;
1092   fCdd = GetCovariance()[0] + s2ylocvtx;   // neglecting correlations
1093   fCdz = GetCovariance()[1];               // between (x,y) and z
1094   fCzz = GetCovariance()[2] + cov[5];      // in vertex's covariance matrix
1095   //*****
1096
1097   {//Try to constrain
1098     Double_t p[2]={yv,zv}, c[3]={cov[2],0.,cov[5]};
1099     Double_t chi2=GetPredictedChi2(p,c);
1100
1101     if (chi2>77.) return kFALSE;
1102
1103     AliExternalTrackParam tmp(*this);
1104     if (!tmp.Update(p,c)) return kFALSE;
1105
1106     fCchi2=chi2;
1107     if (!fCp) fCp=new AliExternalTrackParam();
1108     new (fCp) AliExternalTrackParam(tmp);
1109   }
1110
1111   return kTRUE;
1112 }
1113
1114 //_______________________________________________________________________
1115 void AliESDtrack::Print(Option_t *) const {
1116   // Prints info on the track
1117
1118   printf("ESD track info\n") ;
1119   Double_t p[AliPID::kSPECIESN] ;
1120   Int_t index = 0 ;
1121   if( IsOn(kITSpid) ){
1122     printf("From ITS: ") ;
1123     GetITSpid(p) ;
1124     for(index = 0 ; index < AliPID::kSPECIES; index++)
1125       printf("%f, ", p[index]) ;
1126     printf("\n           signal = %f\n", GetITSsignal()) ;
1127   }
1128   if( IsOn(kTPCpid) ){
1129     printf("From TPC: ") ;
1130     GetTPCpid(p) ;
1131     for(index = 0 ; index < AliPID::kSPECIES; index++)
1132       printf("%f, ", p[index]) ;
1133     printf("\n           signal = %f\n", GetTPCsignal()) ;
1134   }
1135   if( IsOn(kTRDpid) ){
1136     printf("From TRD: ") ;
1137     GetTRDpid(p) ;
1138     for(index = 0 ; index < AliPID::kSPECIES; index++)
1139       printf("%f, ", p[index]) ;
1140     printf("\n           signal = %f\n", GetTRDsignal()) ;
1141   }
1142   if( IsOn(kTOFpid) ){
1143     printf("From TOF: ") ;
1144     GetTOFpid(p) ;
1145     for(index = 0 ; index < AliPID::kSPECIES; index++)
1146       printf("%f, ", p[index]) ;
1147     printf("\n           signal = %f\n", GetTOFsignal()) ;
1148   }
1149   if( IsOn(kHMPIDpid) ){
1150     printf("From HMPID: ") ;
1151     GetHMPIDpid(p) ;
1152     for(index = 0 ; index < AliPID::kSPECIES; index++)
1153       printf("%f, ", p[index]) ;
1154     printf("\n           signal = %f\n", GetHMPIDsignal()) ;
1155   }
1156 }
1157