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