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