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