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