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