STEER/AOD/AliAODTrack.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2007, 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 /* $Id$ */
  17
  18 //-------------------------------------------------------------------------
  19 //     AOD track implementation of AliVTrack
  20 //     Author: Markus Oldenburg, CERN
  21 //     Markus.Oldenburg@cern.ch
  22 //-------------------------------------------------------------------------
  23
  24 #include <TVector3.h>
  25 #include "AliLog.h"
  26 #include "AliExternalTrackParam.h"
  27 #include "AliVVertex.h"
  28 #include "AliAODTrack.h"
  29
  30 ClassImp(AliAODTrack)
  31
  32 //______________________________________________________________________________
  33 AliAODTrack::AliAODTrack() :
  34   AliVTrack(),
  35   fRAtAbsorberEnd(0.),
  36   fChi2perNDF(-999.),
  37   fChi2MatchTrigger(0.),
  38   fFlags(0),
  39   fLabel(-999),
  40   fITSMuonClusterMap(0),
  41   fMUONtrigHitsMapTrg(0),
  42   fMUONtrigHitsMapTrk(0),
  43   fFilterMap(0),
  44   fTPCFitMap(),
  45   fTPCClusterMap(),
  46   fTPCSharedMap(),
  47   fTPCnclsF(0),
  48   fID(-999),
  49   fCharge(-99),
  50   fType(kUndef),
  51   fCaloIndex(kEMCALNoMatch),
  52   fCovMatrix(NULL),
  53   fDetPid(NULL),
  54   fProdVertex(NULL),
  55   fTrackPhiOnEMCal(-999),
  56   fTrackEtaOnEMCal(-999)
  57 {
  58   // default constructor
  59
  60   SetP();
  61   SetPosition((Float_t*)NULL);
  62   SetXYAtDCA(-999., -999.);
  63   SetPxPyPzAtDCA(-999., -999., -999.);
  64   SetPID((Float_t*)NULL);
  65 }
  66
  67 //______________________________________________________________________________
  68 AliAODTrack::AliAODTrack(Short_t id,
  69                          Int_t label,
  70                          Double_t p[3],
  71                          Bool_t cartesian,
  72                          Double_t x[3],
  73                          Bool_t isDCA,
  74                          Double_t covMatrix[21],
  75                          Short_t charge,
  76                          UChar_t itsClusMap,
  77                          Double_t pid[10],
  78                          AliAODVertex *prodVertex,
  79                          Bool_t usedForVtxFit,
  80                          Bool_t usedForPrimVtxFit,
  81                          AODTrk_t ttype,
  82                          UInt_t selectInfo,
  83                          Float_t chi2perNDF) :
  84   AliVTrack(),
  85   fRAtAbsorberEnd(0.),
  86   fChi2perNDF(chi2perNDF),
  87   fChi2MatchTrigger(0.),
  88   fFlags(0),
  89   fLabel(label),
  90   fITSMuonClusterMap(0),
  91   fMUONtrigHitsMapTrg(0),
  92   fMUONtrigHitsMapTrk(0),
  93   fFilterMap(selectInfo),
  94   fTPCFitMap(),
  95   fTPCClusterMap(),
  96   fTPCSharedMap(),
  97   fTPCnclsF(0),
  98   fID(id),
  99   fCharge(charge),
 100   fType(ttype),
 101   fCaloIndex(kEMCALNoMatch),
 102   fCovMatrix(NULL),
 103   fDetPid(NULL),
 104   fProdVertex(prodVertex),
 105   fTrackPhiOnEMCal(-999),
 106   fTrackEtaOnEMCal(-999)
 107 {
 108   // constructor
 109
 110   SetP(p, cartesian);
 111   SetPosition(x, isDCA);
 112   SetXYAtDCA(-999., -999.);
 113   SetPxPyPzAtDCA(-999., -999., -999.);
 114   SetUsedForVtxFit(usedForVtxFit);
 115   SetUsedForPrimVtxFit(usedForPrimVtxFit);
 116   if(covMatrix) SetCovMatrix(covMatrix);
 117   SetPID(pid);
 118   SetITSClusterMap(itsClusMap);
 119 }
 120
 121 //______________________________________________________________________________
 122 AliAODTrack::AliAODTrack(Short_t id,
 123                          Int_t label,
 124                          Float_t p[3],
 125                          Bool_t cartesian,
 126                          Float_t x[3],
 127                          Bool_t isDCA,
 128                          Float_t covMatrix[21],
 129                          Short_t charge,
 130                          UChar_t itsClusMap,
 131                          Float_t pid[10],
 132                          AliAODVertex *prodVertex,
 133                          Bool_t usedForVtxFit,
 134                          Bool_t usedForPrimVtxFit,
 135                          AODTrk_t ttype,
 136                          UInt_t selectInfo,
 137                          Float_t chi2perNDF) :
 138   AliVTrack(),
 139   fRAtAbsorberEnd(0.),
 140   fChi2perNDF(chi2perNDF),
 141   fChi2MatchTrigger(0.),
 142   fFlags(0),
 143   fLabel(label),
 144   fITSMuonClusterMap(0),
 145   fMUONtrigHitsMapTrg(0),
 146   fMUONtrigHitsMapTrk(0),
 147   fFilterMap(selectInfo),
 148   fTPCFitMap(),
 149   fTPCClusterMap(),
 150   fTPCSharedMap(),
 151   fTPCnclsF(0),
 152   fID(id),
 153   fCharge(charge),
 154   fType(ttype),
 155   fCaloIndex(kEMCALNoMatch),
 156   fCovMatrix(NULL),
 157   fDetPid(NULL),
 158   fProdVertex(prodVertex),
 159   fTrackPhiOnEMCal(-999),
 160   fTrackEtaOnEMCal(-999)
 161 {
 162   // constructor
 163
 164   SetP(p, cartesian);
 165   SetPosition(x, isDCA);
 166   SetXYAtDCA(-999., -999.);
 167   SetPxPyPzAtDCA(-999., -999., -999.);
 168   SetUsedForVtxFit(usedForVtxFit);
 169   SetUsedForPrimVtxFit(usedForPrimVtxFit);
 170   if(covMatrix) SetCovMatrix(covMatrix);
 171   SetPID(pid);
 172   SetITSClusterMap(itsClusMap);
 173 }
 174
 175 //______________________________________________________________________________
 176 AliAODTrack::~AliAODTrack()
 177 {
 178   // destructor
 179   delete fCovMatrix;
 180   delete fDetPid;
 181 }
 182
 183
 184 //______________________________________________________________________________
 185 AliAODTrack::AliAODTrack(const AliAODTrack& trk) :
 186   AliVTrack(trk),
 187   fRAtAbsorberEnd(trk.fRAtAbsorberEnd),
 188   fChi2perNDF(trk.fChi2perNDF),
 189   fChi2MatchTrigger(trk.fChi2MatchTrigger),
 190   fFlags(trk.fFlags),
 191   fLabel(trk.fLabel),
 192   fITSMuonClusterMap(trk.fITSMuonClusterMap),
 193   fMUONtrigHitsMapTrg(trk.fMUONtrigHitsMapTrg),
 194   fMUONtrigHitsMapTrk(trk.fMUONtrigHitsMapTrk),
 195   fFilterMap(trk.fFilterMap),
 196   fTPCFitMap(trk.fTPCFitMap),
 197   fTPCClusterMap(trk.fTPCClusterMap),
 198   fTPCSharedMap(trk.fTPCSharedMap),
 199   fTPCnclsF(trk.fTPCnclsF),
 200   fID(trk.fID),
 201   fCharge(trk.fCharge),
 202   fType(trk.fType),
 203   fCaloIndex(trk.fCaloIndex),
 204   fCovMatrix(NULL),
 205   fDetPid(NULL),
 206   fProdVertex(trk.fProdVertex),
 207   fTrackPhiOnEMCal(trk.fTrackPhiOnEMCal),
 208   fTrackEtaOnEMCal(trk.fTrackEtaOnEMCal)
 209 {
 210   // Copy constructor
 211
 212   trk.GetP(fMomentum);
 213   trk.GetPosition(fPosition);
 214   SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
 215   SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
 216   SetUsedForVtxFit(trk.GetUsedForVtxFit());
 217   SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
 218   if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
 219   if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
 220   SetPID(trk.fPID);
 221 }
 222
 223 //______________________________________________________________________________
 224 AliAODTrack& AliAODTrack::operator=(const AliAODTrack& trk)
 225 {
 226   // Assignment operator
 227   if(this!=&trk) {
 228
 229     AliVTrack::operator=(trk);
 230
 231     trk.GetP(fMomentum);
 232     trk.GetPosition(fPosition);
 233     SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
 234     SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
 235     fRAtAbsorberEnd    = trk.fRAtAbsorberEnd;
 236     fChi2perNDF        = trk.fChi2perNDF;
 237     fChi2MatchTrigger  = trk.fChi2MatchTrigger;
 238     trk.GetPID(fPID);
 239     fFlags             = trk.fFlags;
 240     fLabel             = trk.fLabel;
 241     fITSMuonClusterMap = trk.fITSMuonClusterMap;
 242     fMUONtrigHitsMapTrg = trk.fMUONtrigHitsMapTrg;
 243     fMUONtrigHitsMapTrk = trk.fMUONtrigHitsMapTrk;
 244     fFilterMap         = trk.fFilterMap;
 245     fTPCFitMap         = trk.fTPCFitMap;
 246     fTPCClusterMap     = trk.fTPCClusterMap;
 247     fTPCSharedMap      = trk.fTPCSharedMap;
 248     fTPCnclsF          = trk.fTPCnclsF;
 249     fID                = trk.fID;
 250     fCharge            = trk.fCharge;
 251     fType              = trk.fType;
 252     fCaloIndex         = trk.fCaloIndex;
 253     fTrackPhiOnEMCal   = trk.fTrackPhiOnEMCal;
 254     fTrackEtaOnEMCal   = trk.fTrackEtaOnEMCal;
 255
 256     delete fCovMatrix;
 257     if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
 258     else fCovMatrix=NULL;
 259
 260
 261     fProdVertex        = trk.fProdVertex;
 262     SetUsedForVtxFit(trk.GetUsedForVtxFit());
 263     SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
 264
 265     delete fDetPid;
 266     if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
 267     else fDetPid=NULL;
 268
 269
 270   }
 271
 272   return *this;
 273 }
 274
 275 //______________________________________________________________________________
 276 Double_t AliAODTrack::M(AODTrkPID_t pid) const
 277 {
 278   // Returns the mass.
 279   // Masses for nuclei don't exist in the PDG tables, therefore they were put by hand.
 280
 281   switch (pid) {
 282
 283   case kElectron :
 284     return 0.000510999; //TDatabasePDG::Instance()->GetParticle(11/*::kElectron*/)->Mass();
 285     break;
 286
 287   case kMuon :
 288     return 0.1056584; //TDatabasePDG::Instance()->GetParticle(13/*::kMuonMinus*/)->Mass();
 289     break;
 290
 291   case kPion :
 292     return 0.13957; //TDatabasePDG::Instance()->GetParticle(211/*::kPiPlus*/)->Mass();
 293     break;
 294
 295   case kKaon :
 296     return 0.4937; //TDatabasePDG::Instance()->GetParticle(321/*::kKPlus*/)->Mass();
 297     break;
 298
 299   case kProton :
 300     return 0.9382720; //TDatabasePDG::Instance()->GetParticle(2212/*::kProton*/)->Mass();
 301     break;
 302
 303   case kDeuteron :
 304     return 1.8756; //TDatabasePDG::Instance()->GetParticle(1000010020)->Mass();
 305     break;
 306
 307   case kTriton :
 308     return 2.8089; //TDatabasePDG::Instance()->GetParticle(1000010030)->Mass();
 309     break;
 310
 311   case kHelium3 :
 312     return 2.8084; //TDatabasePDG::Instance()->GetParticle(1000020030)->Mass();
 313     break;
 314
 315   case kAlpha :
 316     return 3.7274; //TDatabasePDG::Instance()->GetParticle(1000020040)->Mass();
 317     break;
 318
 319   case kUnknown :
 320     return -999.;
 321     break;
 322
 323   default :
 324     return -999.;
 325   }
 326 }
 327
 328 //______________________________________________________________________________
 329 Double_t AliAODTrack::E(AODTrkPID_t pid) const
 330 {
 331   // Returns the energy of the particle of a given pid.
 332
 333   if (pid != kUnknown) { // particle was identified
 334     Double_t m = M(pid);
 335     return TMath::Sqrt(P()*P() + m*m);
 336   } else { // pid unknown
 337     return -999.;
 338   }
 339 }
 340
 341 //______________________________________________________________________________
 342 Double_t AliAODTrack::Y(AODTrkPID_t pid) const
 343 {
 344   // Returns the rapidity of a particle of a given pid.
 345
 346   if (pid != kUnknown) { // particle was identified
 347     Double_t e = E(pid);
 348     Double_t pz = Pz();
 349     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
 350       return 0.5*TMath::Log((e+pz)/(e-pz));
 351     } else { // energy not known or equal to pz
 352       return -999.;
 353     }
 354   } else { // pid unknown
 355     return -999.;
 356   }
 357 }
 358
 359 //______________________________________________________________________________
 360 Double_t AliAODTrack::Y(Double_t m) const
 361 {
 362   // Returns the rapidity of a particle of a given mass.
 363
 364   if (m >= 0.) { // mass makes sense
 365     Double_t e = E(m);
 366     Double_t pz = Pz();
 367     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
 368       return 0.5*TMath::Log((e+pz)/(e-pz));
 369     } else { // energy not known or equal to pz
 370       return -999.;
 371     }
 372   } else { // pid unknown
 373     return -999.;
 374   }
 375 }
 376
 377 //______________________________________________________________________________
 378 AliAODTrack::AODTrkPID_t AliAODTrack::GetMostProbablePID() const
 379 {
 380   // Returns the most probable PID array element.
 381
 382   Int_t nPID = 10;
 383   AODTrkPID_t loc = kUnknown;
 384   Double_t max = 0.;
 385   Bool_t allTheSame = kTRUE;
 386
 387   for (Int_t iPID = 0; iPID < nPID; iPID++) {
 388     if (fPID[iPID] >= max) {
 389       if (fPID[iPID] > max) {
 390         allTheSame = kFALSE;
 391         max = fPID[iPID];
 392         loc = (AODTrkPID_t)iPID;
 393       } else {
 394         allTheSame = kTRUE;
 395       }
 396     }
 397   }
 398   return allTheSame ? kUnknown : loc;
 399 }
 400
 401 //______________________________________________________________________________
 402 void AliAODTrack::ConvertAliPIDtoAODPID()
 403 {
 404   // Converts AliPID array.
 405   // The numbering scheme is the same for electrons, muons, pions, kaons, and protons.
 406   // Everything else has to be set to zero.
 407
 408   fPID[kDeuteron] = 0.;
 409   fPID[kTriton]   = 0.;
 410   fPID[kHelium3]  = 0.;
 411   fPID[kAlpha]    = 0.;
 412   fPID[kUnknown]  = 0.;
 413
 414   return;
 415 }
 416
 417
 418 //______________________________________________________________________________
 419 template <class T> void AliAODTrack::SetP(const T *p, const Bool_t cartesian)
 420 {
 421   // Set the momentum
 422
 423   if (p) {
 424     if (cartesian) {
 425       Double_t pt2 = p[0]*p[0] + p[1]*p[1];
 426       Double_t pp  = TMath::Sqrt(pt2 + p[2]*p[2]);
 427
 428       fMomentum[0] = TMath::Sqrt(pt2); // pt
 429       fMomentum[1] = (pt2 != 0.) ? TMath::Pi()+TMath::ATan2(-p[1], -p[0]) : -999; // phi
 430       fMomentum[2] = (pp != 0.) ? TMath::ACos(p[2] / pp) : -999.; // theta
 431     } else {
 432       fMomentum[0] = p[0];  // pt
 433       fMomentum[1] = p[1];  // phi
 434       fMomentum[2] = p[2];  // theta
 435     }
 436   } else {
 437     fMomentum[0] = -999.;
 438     fMomentum[1] = -999.;
 439     fMomentum[2] = -999.;
 440   }
 441 }
 442
 443 //______________________________________________________________________________
 444 template <class T> void AliAODTrack::SetPosition(const T *x, const Bool_t dca)
 445 {
 446   // set the position
 447
 448   if (x) {
 449     if (!dca) {
 450       ResetBit(kIsDCA);
 451
 452       fPosition[0] = x[0];
 453       fPosition[1] = x[1];
 454       fPosition[2] = x[2];
 455     } else {
 456       SetBit(kIsDCA);
 457       // don't know any better yet
 458       fPosition[0] = -999.;
 459       fPosition[1] = -999.;
 460       fPosition[2] = -999.;
 461     }
 462   } else {
 463     ResetBit(kIsDCA);
 464
 465     fPosition[0] = -999.;
 466     fPosition[1] = -999.;
 467     fPosition[2] = -999.;
 468   }
 469 }
 470
 471 //______________________________________________________________________________
 472 void AliAODTrack::SetDCA(Double_t d, Double_t z)
 473 {
 474   // set the dca
 475   fPosition[0] = d;
 476   fPosition[1] = z;
 477   fPosition[2] = 0.;
 478   SetBit(kIsDCA);
 479 }
 480
 481 //______________________________________________________________________________
 482 void AliAODTrack::Print(Option_t* /* option */) const
 483 {
 484   // prints information about AliAODTrack
 485
 486   printf("Object name: %s   Track type: %s\n", GetName(), GetTitle());
 487   printf("        px = %f\n", Px());
 488   printf("        py = %f\n", Py());
 489   printf("        pz = %f\n", Pz());
 490   printf("        pt = %f\n", Pt());
 491   printf("      1/pt = %f\n", OneOverPt());
 492   printf("     theta = %f\n", Theta());
 493   printf("       phi = %f\n", Phi());
 494   printf("  chi2/NDF = %f\n", Chi2perNDF());
 495   printf("    charge = %d\n", Charge());
 496 }
 497
 498 //______________________________________________________________________________
 499 void AliAODTrack::SetMatchTrigger(Int_t matchTrig)
 500 {
 501   // Set the MUON trigger information
 502   switch(matchTrig){
 503     case 0: // 0 track does not match trigger
 504       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
 505       break;
 506     case 1: // 1 track match but does not pass pt cut
 507       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x40000000;
 508       break;
 509     case 2: // 2 track match Low pt cut
 510       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x80000000;
 511       break;
 512     case 3: // 3 track match High pt cut
 513       fITSMuonClusterMap=fITSMuonClusterMap|0xc0000000;
 514       break;
 515     default:
 516       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
 517       AliWarning(Form("unknown case for matchTrig: %d\n",matchTrig));
 518   }
 519 }
 520
 521 //______________________________________________________________________________
 522 Bool_t AliAODTrack::HitsMuonChamber(Int_t MuonChamber, Int_t cathode) const
 523 {
 524   // return kTRUE if the track fires the given tracking or trigger chamber.
 525   // If the chamber is a trigger one:
 526   // - if cathode = 0 or 1, the track matches the corresponding cathode
 527   // - if cathode = -1, the track matches both cathodes
 528
 529   if (MuonChamber < 0) return kFALSE;
 530
 531   if (MuonChamber < 10) return TESTBIT(GetMUONClusterMap(), MuonChamber);
 532
 533   if (MuonChamber < 14) {
 534
 535     if (cathode < 0) return TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber) &&
 536                             TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber+4);
 537
 538     if (cathode < 2) return TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber+(1-cathode)*4);
 539
 540   }
 541
 542   return kFALSE;
 543 }
 544
 545 //______________________________________________________________________________
 546 Bool_t AliAODTrack::MatchTriggerDigits() const
 547 {
 548   // return kTRUE if the track matches a digit on both planes of at least 2 trigger chambers
 549
 550   Int_t nMatchedChambers = 0;
 551   for (Int_t ich=10; ich<14; ich++) if (HitsMuonChamber(ich)) nMatchedChambers++;
 552
 553   return (nMatchedChambers >= 2);
 554 }
 555
 556 //______________________________________________________________________________
 557 Bool_t AliAODTrack::PropagateToDCA(const AliVVertex *vtx,
 558     Double_t b, Double_t maxd, Double_t dz[2], Double_t covar[3])
 559 {
 560   // compute impact parameters to the vertex vtx and their covariance matrix
 561   // b is the Bz, needed to propagate correctly the track to vertex
 562   // only the track parameters are update after the propagation (pos and mom),
 563   // not the covariance matrix. This is OK for propagation over short distance
 564   // inside the beam pipe.
 565   // return kFALSE is something went wrong
 566
 567   // convert to AliExternalTrackParam
 568   AliExternalTrackParam etp; etp.CopyFromVTrack(this);
 569
 570   Float_t xstart = etp.GetX();
 571   if(xstart>3.) {
 572     AliError("This method can be used only for propagation inside the beam pipe");
 573     return kFALSE;
 574   }
 575
 576   if(!etp.PropagateToDCA(vtx,b,maxd,dz,covar)) return kFALSE;
 577
 578   // update track position and momentum
 579   Double_t mom[3];
 580   etp.GetPxPyPz(mom);
 581   SetP(mom,kTRUE);
 582   etp.GetXYZ(mom);
 583   SetPosition(mom,kFALSE);
 584
 585
 586   return kTRUE;
 587 }
 588
 589 //______________________________________________________________________________
 590 Bool_t AliAODTrack::GetPxPyPz(Double_t p[3]) const
 591 {
 592     //---------------------------------------------------------------------
 593     // This function returns the global track momentum components
 594     //---------------------------------------------------------------------
 595   p[0]=Px(); p[1]=Py(); p[2]=Pz();
 596   return kTRUE;
 597 }
 598
 599
 600 //_______________________________________________________________________
 601 Float_t AliAODTrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1, Int_t bitType ) const
 602 {
 603   //
 604   // TPC cluster information
 605   // type 0: get fraction of found/findable clusters with neighbourhood definition
 606   //      1: findable clusters with neighbourhood definition
 607   //      2: found clusters
 608   // bitType:
 609   //      0 - all cluster used
 610   //      1 - clusters  used for the kalman update
 611   // definition of findable clusters:
 612   //            a cluster is defined as findable if there is another cluster
 613   //           within +- nNeighbours pad rows. The idea is to overcome threshold
 614   //           effects with a very simple algorithm.
 615   //
 616
 617
 618   Int_t found=0;
 619   Int_t findable=0;
 620   Int_t last=-nNeighbours;
 621   const TBits & clusterMap = (bitType%2==0) ? fTPCClusterMap : fTPCFitMap;
 622
 623   Int_t upperBound=clusterMap.GetNbits();
 624   if (upperBound>row1) upperBound=row1;
 625   for (Int_t i=row0; i<upperBound; ++i){
 626     //look to current row
 627     if (clusterMap[i]) {
 628       last=i;
 629       ++found;
 630       ++findable;
 631       continue;
 632     }
 633     //look to nNeighbours before
 634     if ((i-last)<=nNeighbours) {
 635       ++findable;
 636       continue;
 637     }
 638     //look to nNeighbours after
 639     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
 640       if (clusterMap[j]){
 641         ++findable;
 642         break;
 643       }
 644     }
 645   }
 646   if (type==2) return found;
 647   if (type==1) return findable;
 648
 649   if (type==0){
 650     Float_t fraction=0;
 651     if (findable>0)
 652       fraction=(Float_t)found/(Float_t)findable;
 653     else
 654       fraction=0;
 655     return fraction;
 656   }
 657   return 0;  // undefined type - default value
 658 }
 659
 660
 661 //______________________________________________________________________________
 662 Double_t  AliAODTrack::GetTRDslice(Int_t plane, Int_t slice) const {
 663   //
 664   // return TRD Pid information
 665   //
 666   if (!fDetPid) return -1;
 667   Double32_t *trdSlices=fDetPid->GetTRDsignal();
 668   if (!trdSlices) return -1;
 669   if ((plane<0) || (plane>=kTRDnPlanes)) {
 670     return -1.;
 671   }
 672
 673   Int_t ns=fDetPid->GetTRDnSlices()/kTRDnPlanes;
 674   if ((slice<-1) || (slice>=ns)) {
 675     return -1.;
 676   }
 677
 678   if(slice>=0) return trdSlices[plane*ns + slice];
 679
 680   // return average of the dEdx measurements
 681   Double_t q=0.; Double32_t *s = &trdSlices[plane*ns];
 682   for (Int_t i=0; i<ns; i++, s++) if((*s)>0.) q+=(*s);
 683   return q/ns;
 684 }
 685
 686 //______________________________________________________________________________
 687 UChar_t AliAODTrack::GetTRDntrackletsPID() const{
 688   //
 689   // return number of tracklets calculated from the slices
 690   //
 691   if(!fDetPid) return -1;
 692   return fDetPid->GetTRDntrackletsPID();
 693 }
 694
 695 //______________________________________________________________________________
 696 UChar_t AliAODTrack::GetTRDncls(Int_t layer) const {
 697   //
 698   // return number of TRD clusters
 699   //
 700   if(!fDetPid || layer > 5) return -1;
 701   if(layer < 0) return fDetPid->GetTRDncls();
 702   else return fDetPid->GetTRDncls(layer);
 703 }
 704
 705 //______________________________________________________________________________
 706 Double_t AliAODTrack::GetTRDmomentum(Int_t plane, Double_t */*sp*/) const
 707 {
 708   //Returns momentum estimation
 709   // in TRD layer "plane".
 710
 711   if (!fDetPid) return -1;
 712   const Float_t *trdMomentum=fDetPid->GetTRDmomentum();
 713
 714   if (!trdMomentum) {
 715     return -1.;
 716   }
 717   if ((plane<0) || (plane>=kTRDnPlanes)) {
 718     return -1.;
 719   }
 720
 721   return trdMomentum[plane];
 722 }
 723
 724 //_______________________________________________________________________
 725 Int_t AliAODTrack::GetTOFBunchCrossing(Double_t b, Bool_t) const
 726 {
 727   // Returns the number of bunch crossings after trigger (assuming 25ns spacing)
 728   const double kSpacing = 25e3; // min interbanch spacing
 729   const double kShift = 0;
 730   Int_t bcid = kTOFBCNA; // defualt one
 731   if (!IsOn(kTOFout) || !IsOn(kESDpid)) return bcid; // no info
 732   //
 733   double tdif = GetTOFsignal();
 734   if (IsOn(kTIME)) { // integrated time info is there
 735     int pid = (int)GetMostProbablePID();
 736     double ttimes[10];
 737     GetIntegratedTimes(ttimes);
 738     tdif -= ttimes[pid];
 739   }
 740   else { // assume integrated time info from TOF radius and momentum
 741     const double kRTOF = 385.;
 742     const double kCSpeed = 3.e-2; // cm/ps
 743     double p = P();
 744     if (p<0.001) p = 1.0;
 745     double m = M();
 746     double path =  kRTOF;     // mean TOF radius
 747     if (TMath::Abs(b)>kAlmost0) {  // account for curvature
 748       double curv = Pt()/(b*kB2C);
 749       if (curv>kAlmost0) {
 750         double tgl = Pz()/Pt();
 751         path = 2./curv*TMath::ASin(kRTOF*curv/2.)*TMath::Sqrt(1.+tgl*tgl);
 752       }
 753     }
 754     tdif -= path/kCSpeed*TMath::Sqrt(1.+m*m/(p*p));
 755   }
 756   bcid = TMath::Nint((tdif - kShift)/kSpacing);
 757   return bcid;
 758 }
 759
 760 //_____________________________________________________________________________
 761 Bool_t AliAODTrack::GetXYZAt(Double_t x, Double_t b, Double_t *r) const
 762 {
 763   //---------------------------------------------------------------------
 764   // This function returns the global track position extrapolated to
 765   // the radial position "x" (cm) in the magnetic field "b" (kG)
 766   //---------------------------------------------------------------------
 767
 768   //conversion of track parameter representation is
 769   //based on the implementation of AliExternalTrackParam::Set(...)
 770   //maybe some of this code can be moved to AliVTrack to avoid code duplication
 771   const double kSafe = 1e-5;
 772   Double_t alpha=0.0;
 773   Double_t radPos2 = fPosition[0]*fPosition[0]+fPosition[1]*fPosition[1];
 774   Double_t radMax  = 45.; // approximately ITS outer radius
 775   if (radPos2 < radMax*radMax) { // inside the ITS
 776      alpha = TMath::ATan2(fMomentum[1],fMomentum[0]);
 777   } else { // outside the ITS
 778      Float_t phiPos = TMath::Pi()+TMath::ATan2(-fPosition[1], -fPosition[0]);
 779      alpha =
 780      TMath::DegToRad()*(20*((((Int_t)(phiPos*TMath::RadToDeg()))/20))+10);
 781   }
 782   //
 783   Double_t cs=TMath::Cos(alpha), sn=TMath::Sin(alpha);
 784   // protection:  avoid alpha being too close to 0 or +-pi/2
 785   if (TMath::Abs(sn)<kSafe) {
 786     alpha = kSafe;
 787     cs=TMath::Cos(alpha);
 788     sn=TMath::Sin(alpha);
 789   }
 790   else if (cs<kSafe) {
 791     alpha -= TMath::Sign(kSafe, alpha);
 792     cs=TMath::Cos(alpha);
 793     sn=TMath::Sin(alpha);
 794   }
 795
 796   // Get the vertex of origin and the momentum
 797   TVector3 ver(fPosition[0],fPosition[1],fPosition[2]);
 798   TVector3 mom(fMomentum[0],fMomentum[1],fMomentum[2]);
 799   //
 800   // avoid momenta along axis
 801   if (TMath::Abs(mom[0])<kSafe) mom[0] = TMath::Sign(kSafe*TMath::Abs(mom[1]), mom[0]);
 802   if (TMath::Abs(mom[1])<kSafe) mom[1] = TMath::Sign(kSafe*TMath::Abs(mom[0]), mom[1]);
 803
 804   // Rotate to the local coordinate system
 805   ver.RotateZ(-alpha);
 806   mom.RotateZ(-alpha);
 807
 808   Double_t param0 = ver.Y();
 809   Double_t param1 = ver.Z();
 810   Double_t param2 = TMath::Sin(mom.Phi());
 811   Double_t param3 = mom.Pz()/mom.Pt();
 812   Double_t param4 = TMath::Sign(1/mom.Pt(),(Double_t)fCharge);
 813
 814   //calculate the propagated coordinates
 815   //this is based on AliExternalTrackParam::GetXYZAt(Double_t x, Double_t b, Double_t *r)
 816   Double_t dx=x-ver.X();
 817   if(TMath::Abs(dx)<=kAlmost0) return GetXYZ(r);
 818
 819   Double_t f1=param2;
 820   Double_t f2=f1 + dx*param4*b*kB2C;
 821
 822   if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
 823   if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
 824
 825   Double_t r1=TMath::Sqrt((1.-f1)*(1.+f1)), r2=TMath::Sqrt((1.-f2)*(1.+f2));
 826   r[0] = x;
 827   r[1] = param0 + dx*(f1+f2)/(r1+r2);
 828   r[2] = param1 + dx*(r2 + f2*(f1+f2)/(r1+r2))*param3;//Thanks to Andrea & Peter
 829
 830   return Local2GlobalPosition(r,alpha);
 831 }
 832
 833
 834