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