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