STEER/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 "AliLog.h"
  25 #include "AliExternalTrackParam.h"
  26 #include "AliVVertex.h"
  27 #include "AliAODTrack.h"
  28
  29 ClassImp(AliAODTrack)
  30
  31 //______________________________________________________________________________
  32 AliAODTrack::AliAODTrack() :
  33   AliVTrack(),
  34   fRAtAbsorberEnd(0.),
  35   fChi2perNDF(-999.),
  36   fChi2MatchTrigger(0.),
  37   fFlags(0),
  38   fLabel(-999),
  39   fITSMuonClusterMap(0),
  40   fFilterMap(0),
  41   fTPCClusterMap(),
  42   fTPCSharedMap(),
  43   fTPCnclsF(0),
  44   fID(-999),
  45   fCharge(-99),
  46   fType(kUndef),
  47   fCovMatrix(NULL),
  48   fDetPid(NULL),
  49   fProdVertex(NULL)
  50 {
  51   // default constructor
  52
  53   SetP();
  54   SetPosition((Float_t*)NULL);
  55   SetXYAtDCA(-999., -999.);
  56   SetPxPyPzAtDCA(-999., -999., -999.);
  57   SetPID((Float_t*)NULL);
  58 }
  59
  60 //______________________________________________________________________________
  61 AliAODTrack::AliAODTrack(Short_t id,
  62                          Int_t label,
  63                          Double_t p[3],
  64                          Bool_t cartesian,
  65                          Double_t x[3],
  66                          Bool_t isDCA,
  67                          Double_t covMatrix[21],
  68                          Short_t charge,
  69                          UChar_t itsClusMap,
  70                          Double_t pid[10],
  71                          AliAODVertex *prodVertex,
  72                          Bool_t usedForVtxFit,
  73                          Bool_t usedForPrimVtxFit,
  74                          AODTrk_t ttype,
  75                          UInt_t selectInfo,
  76                          Float_t chi2perNDF) :
  77   AliVTrack(),
  78   fRAtAbsorberEnd(0.),
  79   fChi2perNDF(chi2perNDF),
  80   fChi2MatchTrigger(0.),
  81   fFlags(0),
  82   fLabel(label),
  83   fITSMuonClusterMap(0),
  84   fFilterMap(selectInfo),
  85   fTPCClusterMap(),
  86   fTPCSharedMap(),
  87   fTPCnclsF(0),
  88   fID(id),
  89   fCharge(charge),
  90   fType(ttype),
  91   fCovMatrix(NULL),
  92   fDetPid(NULL),
  93   fProdVertex(prodVertex)
  94 {
  95   // constructor
  96
  97   SetP(p, cartesian);
  98   SetPosition(x, isDCA);
  99   SetXYAtDCA(-999., -999.);
 100   SetPxPyPzAtDCA(-999., -999., -999.);
 101   SetUsedForVtxFit(usedForVtxFit);
 102   SetUsedForPrimVtxFit(usedForPrimVtxFit);
 103   if(covMatrix) SetCovMatrix(covMatrix);
 104   SetPID(pid);
 105   SetITSClusterMap(itsClusMap);
 106 }
 107
 108 //______________________________________________________________________________
 109 AliAODTrack::AliAODTrack(Short_t id,
 110                          Int_t label,
 111                          Float_t p[3],
 112                          Bool_t cartesian,
 113                          Float_t x[3],
 114                          Bool_t isDCA,
 115                          Float_t covMatrix[21],
 116                          Short_t charge,
 117                          UChar_t itsClusMap,
 118                          Float_t pid[10],
 119                          AliAODVertex *prodVertex,
 120                          Bool_t usedForVtxFit,
 121                          Bool_t usedForPrimVtxFit,
 122                          AODTrk_t ttype,
 123                          UInt_t selectInfo,
 124                          Float_t chi2perNDF) :
 125   AliVTrack(),
 126   fRAtAbsorberEnd(0.),
 127   fChi2perNDF(chi2perNDF),
 128   fChi2MatchTrigger(0.),
 129   fFlags(0),
 130   fLabel(label),
 131   fITSMuonClusterMap(0),
 132   fFilterMap(selectInfo),
 133   fTPCClusterMap(),
 134   fTPCSharedMap(),
 135   fTPCnclsF(0),
 136   fID(id),
 137   fCharge(charge),
 138   fType(ttype),
 139   fCovMatrix(NULL),
 140   fDetPid(NULL),
 141   fProdVertex(prodVertex)
 142 {
 143   // constructor
 144
 145   SetP(p, cartesian);
 146   SetPosition(x, isDCA);
 147   SetXYAtDCA(-999., -999.);
 148   SetPxPyPzAtDCA(-999., -999., -999.);
 149   SetUsedForVtxFit(usedForVtxFit);
 150   SetUsedForPrimVtxFit(usedForPrimVtxFit);
 151   if(covMatrix) SetCovMatrix(covMatrix);
 152   SetPID(pid);
 153   SetITSClusterMap(itsClusMap);
 154 }
 155
 156 //______________________________________________________________________________
 157 AliAODTrack::~AliAODTrack()
 158 {
 159   // destructor
 160   delete fCovMatrix;
 161   delete fDetPid;
 162 }
 163
 164
 165 //______________________________________________________________________________
 166 AliAODTrack::AliAODTrack(const AliAODTrack& trk) :
 167   AliVTrack(trk),
 168   fRAtAbsorberEnd(trk.fRAtAbsorberEnd),
 169   fChi2perNDF(trk.fChi2perNDF),
 170   fChi2MatchTrigger(trk.fChi2MatchTrigger),
 171   fFlags(trk.fFlags),
 172   fLabel(trk.fLabel),
 173   fITSMuonClusterMap(trk.fITSMuonClusterMap),
 174   fFilterMap(trk.fFilterMap),
 175   fTPCClusterMap(trk.fTPCClusterMap),
 176   fTPCSharedMap(trk.fTPCSharedMap),
 177   fTPCnclsF(trk.fTPCnclsF),
 178   fID(trk.fID),
 179   fCharge(trk.fCharge),
 180   fType(trk.fType),
 181   fCovMatrix(NULL),
 182   fDetPid(NULL),
 183   fProdVertex(trk.fProdVertex)
 184 {
 185   // Copy constructor
 186
 187   trk.GetP(fMomentum);
 188   trk.GetPosition(fPosition);
 189   SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
 190   SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
 191   SetUsedForVtxFit(trk.GetUsedForVtxFit());
 192   SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
 193   if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
 194   if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
 195   SetPID(trk.fPID);
 196 }
 197
 198 //______________________________________________________________________________
 199 AliAODTrack& AliAODTrack::operator=(const AliAODTrack& trk)
 200 {
 201   // Assignment operator
 202   if(this!=&trk) {
 203
 204     AliVTrack::operator=(trk);
 205
 206     trk.GetP(fMomentum);
 207     trk.GetPosition(fPosition);
 208     trk.GetPID(fPID);
 209
 210     SetXYAtDCA(trk.XAtDCA(), trk.YAtDCA());
 211     SetPxPyPzAtDCA(trk.PxAtDCA(), trk.PyAtDCA(), trk.PzAtDCA());
 212
 213     fRAtAbsorberEnd = trk.fRAtAbsorberEnd;
 214
 215     fChi2perNDF = trk.fChi2perNDF;
 216     fChi2MatchTrigger = trk.fChi2MatchTrigger;
 217
 218     fFlags = trk.fFlags;
 219     fLabel = trk.fLabel;
 220
 221     fITSMuonClusterMap = trk.fITSMuonClusterMap;
 222     fFilterMap = trk.fFilterMap;
 223
 224     fID = trk.fID;
 225
 226     fCharge = trk.fCharge;
 227     fType = trk.fType;
 228
 229     delete fCovMatrix;
 230     if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
 231     else fCovMatrix=NULL;
 232     fProdVertex = trk.fProdVertex;
 233
 234     SetUsedForVtxFit(trk.GetUsedForVtxFit());
 235     SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
 236
 237     delete fDetPid;
 238     if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
 239     else fDetPid=NULL;
 240   }
 241
 242   return *this;
 243 }
 244
 245 //______________________________________________________________________________
 246 Double_t AliAODTrack::M(AODTrkPID_t pid) const
 247 {
 248   // Returns the mass.
 249   // Masses for nuclei don't exist in the PDG tables, therefore they were put by hand.
 250
 251   switch (pid) {
 252
 253   case kElectron :
 254     return 0.000510999; //TDatabasePDG::Instance()->GetParticle(11/*::kElectron*/)->Mass();
 255     break;
 256
 257   case kMuon :
 258     return 0.1056584; //TDatabasePDG::Instance()->GetParticle(13/*::kMuonMinus*/)->Mass();
 259     break;
 260
 261   case kPion :
 262     return 0.13957; //TDatabasePDG::Instance()->GetParticle(211/*::kPiPlus*/)->Mass();
 263     break;
 264
 265   case kKaon :
 266     return 0.4937; //TDatabasePDG::Instance()->GetParticle(321/*::kKPlus*/)->Mass();
 267     break;
 268
 269   case kProton :
 270     return 0.9382720; //TDatabasePDG::Instance()->GetParticle(2212/*::kProton*/)->Mass();
 271     break;
 272
 273   case kDeuteron :
 274     return 1.8756; //TDatabasePDG::Instance()->GetParticle(1000010020)->Mass();
 275     break;
 276
 277   case kTriton :
 278     return 2.8089; //TDatabasePDG::Instance()->GetParticle(1000010030)->Mass();
 279     break;
 280
 281   case kHelium3 :
 282     return 2.8084; //TDatabasePDG::Instance()->GetParticle(1000020030)->Mass();
 283     break;
 284
 285   case kAlpha :
 286     return 3.7274; //TDatabasePDG::Instance()->GetParticle(1000020040)->Mass();
 287     break;
 288
 289   case kUnknown :
 290     return -999.;
 291     break;
 292
 293   default :
 294     return -999.;
 295   }
 296 }
 297
 298 //______________________________________________________________________________
 299 Double_t AliAODTrack::E(AODTrkPID_t pid) const
 300 {
 301   // Returns the energy of the particle of a given pid.
 302
 303   if (pid != kUnknown) { // particle was identified
 304     Double_t m = M(pid);
 305     return TMath::Sqrt(P()*P() + m*m);
 306   } else { // pid unknown
 307     return -999.;
 308   }
 309 }
 310
 311 //______________________________________________________________________________
 312 Double_t AliAODTrack::Y(AODTrkPID_t pid) const
 313 {
 314   // Returns the rapidity of a particle of a given pid.
 315
 316   if (pid != kUnknown) { // particle was identified
 317     Double_t e = E(pid);
 318     Double_t pz = Pz();
 319     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
 320       return 0.5*TMath::Log((e+pz)/(e-pz));
 321     } else { // energy not known or equal to pz
 322       return -999.;
 323     }
 324   } else { // pid unknown
 325     return -999.;
 326   }
 327 }
 328
 329 //______________________________________________________________________________
 330 Double_t AliAODTrack::Y(Double_t m) const
 331 {
 332   // Returns the rapidity of a particle of a given mass.
 333
 334   if (m >= 0.) { // mass makes sense
 335     Double_t e = E(m);
 336     Double_t pz = Pz();
 337     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
 338       return 0.5*TMath::Log((e+pz)/(e-pz));
 339     } else { // energy not known or equal to pz
 340       return -999.;
 341     }
 342   } else { // pid unknown
 343     return -999.;
 344   }
 345 }
 346
 347 //______________________________________________________________________________
 348 AliAODTrack::AODTrkPID_t AliAODTrack::GetMostProbablePID() const
 349 {
 350   // Returns the most probable PID array element.
 351
 352   Int_t nPID = 10;
 353   AODTrkPID_t loc = kUnknown;
 354   Double_t max = 0.;
 355   Bool_t allTheSame = kTRUE;
 356
 357   for (Int_t iPID = 0; iPID < nPID; iPID++) {
 358     if (fPID[iPID] >= max) {
 359       if (fPID[iPID] > max) {
 360         allTheSame = kFALSE;
 361         max = fPID[iPID];
 362         loc = (AODTrkPID_t)iPID;
 363       } else {
 364         allTheSame = kTRUE;
 365       }
 366     }
 367   }
 368   return allTheSame ? kUnknown : loc;
 369 }
 370
 371 //______________________________________________________________________________
 372 void AliAODTrack::ConvertAliPIDtoAODPID()
 373 {
 374   // Converts AliPID array.
 375   // The numbering scheme is the same for electrons, muons, pions, kaons, and protons.
 376   // Everything else has to be set to zero.
 377
 378   fPID[kDeuteron] = 0.;
 379   fPID[kTriton]   = 0.;
 380   fPID[kHelium3]  = 0.;
 381   fPID[kAlpha]    = 0.;
 382   fPID[kUnknown]  = 0.;
 383
 384   return;
 385 }
 386
 387
 388 //______________________________________________________________________________
 389 template <class T> void AliAODTrack::SetP(const T *p, const Bool_t cartesian)
 390 {
 391   // Set the momentum
 392
 393   if (p) {
 394     if (cartesian) {
 395       Double_t pt2 = p[0]*p[0] + p[1]*p[1];
 396       Double_t pp  = TMath::Sqrt(pt2 + p[2]*p[2]);
 397
 398       fMomentum[0] = TMath::Sqrt(pt2); // pt
 399       fMomentum[1] = (pt2 != 0.) ? TMath::Pi()+TMath::ATan2(-p[1], -p[0]) : -999; // phi
 400       fMomentum[2] = (pp != 0.) ? TMath::ACos(p[2] / pp) : -999.; // theta
 401     } else {
 402       fMomentum[0] = p[0];  // pt
 403       fMomentum[1] = p[1];  // phi
 404       fMomentum[2] = p[2];  // theta
 405     }
 406   } else {
 407     fMomentum[0] = -999.;
 408     fMomentum[1] = -999.;
 409     fMomentum[2] = -999.;
 410   }
 411 }
 412
 413 //______________________________________________________________________________
 414 template <class T> void AliAODTrack::SetPosition(const T *x, const Bool_t dca)
 415 {
 416   // set the position
 417
 418   if (x) {
 419     if (!dca) {
 420       ResetBit(kIsDCA);
 421
 422       fPosition[0] = x[0];
 423       fPosition[1] = x[1];
 424       fPosition[2] = x[2];
 425     } else {
 426       SetBit(kIsDCA);
 427       // don't know any better yet
 428       fPosition[0] = -999.;
 429       fPosition[1] = -999.;
 430       fPosition[2] = -999.;
 431     }
 432   } else {
 433     ResetBit(kIsDCA);
 434
 435     fPosition[0] = -999.;
 436     fPosition[1] = -999.;
 437     fPosition[2] = -999.;
 438   }
 439 }
 440
 441 //______________________________________________________________________________
 442 void AliAODTrack::SetDCA(Double_t d, Double_t z)
 443 {
 444   // set the dca
 445   fPosition[0] = d;
 446   fPosition[1] = z;
 447   fPosition[2] = 0.;
 448   SetBit(kIsDCA);
 449 }
 450
 451 //______________________________________________________________________________
 452 void AliAODTrack::Print(Option_t* /* option */) const
 453 {
 454   // prints information about AliAODTrack
 455
 456   printf("Object name: %s   Track type: %s\n", GetName(), GetTitle());
 457   printf("        px = %f\n", Px());
 458   printf("        py = %f\n", Py());
 459   printf("        pz = %f\n", Pz());
 460   printf("        pt = %f\n", Pt());
 461   printf("      1/pt = %f\n", OneOverPt());
 462   printf("     theta = %f\n", Theta());
 463   printf("       phi = %f\n", Phi());
 464   printf("  chi2/NDF = %f\n", Chi2perNDF());
 465   printf("    charge = %d\n", Charge());
 466 }
 467
 468 //______________________________________________________________________________
 469 void AliAODTrack::SetMatchTrigger(Int_t matchTrig)
 470 {
 471   // Set the MUON trigger information
 472   switch(matchTrig){
 473     case 0: // 0 track does not match trigger
 474       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
 475       break;
 476     case 1: // 1 track match but does not pass pt cut
 477       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x40000000;
 478       break;
 479     case 2: // 2 track match Low pt cut
 480       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x80000000;
 481       break;
 482     case 3: // 3 track match High pt cut
 483       fITSMuonClusterMap=fITSMuonClusterMap|0xc0000000;
 484       break;
 485     default:
 486       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
 487       AliWarning(Form("unknown case for matchTrig: %d\n",matchTrig));
 488   }
 489 }
 490
 491 //______________________________________________________________________________
 492 Bool_t AliAODTrack::HitsMuonChamber(Int_t MuonChamber, Int_t cathode) const
 493 {
 494   // return kTRUE if the track fires the given tracking or trigger chamber.
 495   // If the chamber is a trigger one:
 496   // - if cathode = 0 or 1, the track matches the corresponding cathode
 497   // - if cathode = -1, the track matches both cathodes
 498
 499   if (MuonChamber < 0) return kFALSE;
 500
 501   if (MuonChamber < 10) return TESTBIT(GetMUONClusterMap(), MuonChamber);
 502
 503   if (MuonChamber < 14) {
 504
 505     if (cathode < 0) return TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber) &&
 506                             TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber+4);
 507
 508     if (cathode < 2) return TESTBIT(GetHitsPatternInTrigCh(), 13-MuonChamber+(1-cathode)*4);
 509
 510   }
 511
 512   return kFALSE;
 513 }
 514
 515 //______________________________________________________________________________
 516 Bool_t AliAODTrack::MatchTriggerDigits() const
 517 {
 518   // return kTRUE if the track matches a digit on both planes of at least 2 trigger chambers
 519
 520   Int_t nMatchedChambers = 0;
 521   for (Int_t ich=10; ich<14; ich++) if (HitsMuonChamber(ich)) nMatchedChambers++;
 522
 523   return (nMatchedChambers >= 2);
 524 }
 525
 526 //______________________________________________________________________________
 527 Bool_t AliAODTrack::PropagateToDCA(const AliVVertex *vtx,
 528     Double_t b, Double_t maxd, Double_t dz[2], Double_t covar[3])
 529 {
 530   // compute impact parameters to the vertex vtx and their covariance matrix
 531   // b is the Bz, needed to propagate correctly the track to vertex
 532   // only the track parameters are update after the propagation (pos and mom),
 533   // not the covariance matrix. This is OK for propagation over short distance
 534   // inside the beam pipe.
 535   // return kFALSE is something went wrong
 536
 537   // convert to AliExternalTrackParam
 538   AliExternalTrackParam etp(this);
 539
 540   Float_t xstart = etp.GetX();
 541   if(xstart>3.) {
 542     AliError("This method can be used only for propagation inside the beam pipe");
 543     return kFALSE;
 544   }
 545
 546   if(!etp.PropagateToDCA(vtx,b,maxd,dz,covar)) return kFALSE;
 547
 548   // update track position and momentum
 549   Double_t mom[3];
 550   etp.GetPxPyPz(mom);
 551   SetP(mom,kTRUE);
 552   etp.GetXYZ(mom);
 553   SetPosition(mom,kFALSE);
 554
 555
 556   return kTRUE;
 557 }
 558
 559 //______________________________________________________________________________
 560 Bool_t AliAODTrack::GetPxPyPz(Double_t p[3]) const
 561 {
 562     //---------------------------------------------------------------------
 563     // This function returns the global track momentum components
 564     //---------------------------------------------------------------------
 565   p[0]=Px(); p[1]=Py(); p[2]=Pz();
 566   return kTRUE;
 567 }
 568
 569 //______________________________________________________________________________
 570 Float_t AliAODTrack::GetTPCClusterInfo(Int_t nNeighbours/*=3*/, Int_t type/*=0*/, Int_t row0, Int_t row1) const
 571 {
 572   //
 573   // TPC cluster information
 574   // type 0: get fraction of found/findable clusters with neighbourhood definition
 575   //      1: findable clusters with neighbourhood definition
 576   //      2: found clusters
 577   //
 578   // definition of findable clusters:
 579   //            a cluster is defined as findable if there is another cluster
 580   //           within +- nNeighbours pad rows. The idea is to overcome threshold
 581   //           effects with a very simple algorithm.
 582   //
 583
 584   if (type==2) return fTPCClusterMap.CountBits();
 585
 586   Int_t found=0;
 587   Int_t findable=0;
 588   Int_t last=-nNeighbours;
 589
 590   for (Int_t i=row0; i<row1; ++i){
 591     //look to current row
 592     if (fTPCClusterMap[i]) {
 593       last=i;
 594       ++found;
 595       ++findable;
 596       continue;
 597     }
 598     //look to nNeighbours before
 599     if ((i-last)<=nNeighbours) {
 600       ++findable;
 601       continue;
 602     }
 603     //look to nNeighbours after
 604     for (Int_t j=i+1; j<i+1+nNeighbours; ++j){
 605       if (fTPCClusterMap[j]){
 606         ++findable;
 607         break;
 608       }
 609     }
 610   }
 611   if (type==1) return findable;
 612
 613   if (type==0){
 614     Float_t fraction=0;
 615     if (findable>0)
 616       fraction=(Float_t)found/(Float_t)findable;
 617     else
 618       fraction=0;
 619     return fraction;
 620   }
 621   return 0;  // undefined type - default value
 622 }