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