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 base class
  20 //     Author: Markus Oldenburg, CERN
  21 //-------------------------------------------------------------------------
  22
  23 #include "AliAODTrack.h"
  24
  25 //#include <TPDGCode.h>
  26 //#include <TDatabasePDG.h>
  27
  28 ClassImp(AliAODTrack)
  29
  30 //______________________________________________________________________________
  31 AliAODTrack::AliAODTrack() :
  32   AliVParticle(),
  33   fChi2perNDF(-999.),
  34   fChi2MatchTrigger(0.),
  35   fLabel(-999),
  36   fITSMuonClusterMap(0),
  37   fFilterMap(0),
  38   fID(-999),
  39   fCharge(-99),
  40   fType(kUndef),
  41   fCovMatrix(NULL),
  42   fProdVertex(0x0)
  43 {
  44   // default constructor
  45
  46   SetP();
  47   SetPosition((Float_t*)NULL);
  48   SetPID((Float_t*)NULL);
  49 }
  50
  51 //______________________________________________________________________________
  52 AliAODTrack::AliAODTrack(Short_t id,
  53                          Int_t label,
  54                          Double_t p[3],
  55                          Bool_t cartesian,
  56                          Double_t x[3],
  57                          Bool_t isDCA,
  58                          Double_t covMatrix[21],
  59                          Short_t charge,
  60                          UChar_t itsClusMap,
  61                          Double_t pid[10],
  62                          AliAODVertex *prodVertex,
  63                          Bool_t usedForVtxFit,
  64                          Bool_t usedForPrimVtxFit,
  65                          AODTrk_t ttype,
  66                          UInt_t selectInfo) :
  67   AliVParticle(),
  68   fChi2perNDF(-999.),
  69   fChi2MatchTrigger(0.),
  70   fLabel(label),
  71   fITSMuonClusterMap(itsClusMap),
  72   fFilterMap(selectInfo),
  73   fID(id),
  74   fCharge(charge),
  75   fType(ttype),
  76   fCovMatrix(NULL),
  77   fProdVertex(prodVertex)
  78 {
  79   // constructor
  80
  81   SetP(p, cartesian);
  82   SetPosition(x, isDCA);
  83   SetUsedForVtxFit(usedForVtxFit);
  84   SetUsedForPrimVtxFit(usedForPrimVtxFit);
  85   if(covMatrix) SetCovMatrix(covMatrix);
  86   SetPID(pid);
  87
  88 }
  89
  90 //______________________________________________________________________________
  91 AliAODTrack::AliAODTrack(Short_t id,
  92                          Int_t label,
  93                          Float_t p[3],
  94                          Bool_t cartesian,
  95                          Float_t x[3],
  96                          Bool_t isDCA,
  97                          Float_t covMatrix[21],
  98                          Short_t charge,
  99                          UChar_t itsClusMap,
 100                          Float_t pid[10],
 101                          AliAODVertex *prodVertex,
 102                          Bool_t usedForVtxFit,
 103                          Bool_t usedForPrimVtxFit,
 104                          AODTrk_t ttype,
 105                          UInt_t selectInfo) :
 106   AliVParticle(),
 107   fChi2perNDF(-999.),
 108   fChi2MatchTrigger(0.),
 109   fLabel(label),
 110   fITSMuonClusterMap(itsClusMap),
 111   fFilterMap(selectInfo),
 112   fID(id),
 113   fCharge(charge),
 114   fType(ttype),
 115   fCovMatrix(NULL),
 116   fProdVertex(prodVertex)
 117 {
 118   // constructor
 119
 120   SetP(p, cartesian);
 121   SetPosition(x, isDCA);
 122   SetUsedForVtxFit(usedForVtxFit);
 123   SetUsedForPrimVtxFit(usedForPrimVtxFit);
 124   if(covMatrix) SetCovMatrix(covMatrix);
 125   SetPID(pid);
 126 }
 127
 128 //______________________________________________________________________________
 129 AliAODTrack::~AliAODTrack()
 130 {
 131   // destructor
 132   delete fCovMatrix;
 133 }
 134
 135
 136 //______________________________________________________________________________
 137 AliAODTrack::AliAODTrack(const AliAODTrack& trk) :
 138   AliVParticle(trk),
 139   fChi2perNDF(trk.fChi2perNDF),
 140   fChi2MatchTrigger(trk.fChi2MatchTrigger),
 141   fLabel(trk.fLabel),
 142   fITSMuonClusterMap(trk.fITSMuonClusterMap),
 143   fFilterMap(trk.fFilterMap),
 144   fID(trk.fID),
 145   fCharge(trk.fCharge),
 146   fType(trk.fType),
 147   fCovMatrix(NULL),
 148   fProdVertex(trk.fProdVertex)
 149 {
 150   // Copy constructor
 151
 152   trk.GetP(fMomentum);
 153   trk.GetPosition(fPosition);
 154   SetUsedForVtxFit(trk.GetUsedForVtxFit());
 155   SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
 156   if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
 157   SetPID(trk.fPID);
 158
 159 }
 160
 161 //______________________________________________________________________________
 162 AliAODTrack& AliAODTrack::operator=(const AliAODTrack& trk)
 163 {
 164   // Assignment operator
 165   if(this!=&trk) {
 166
 167     AliVParticle::operator=(trk);
 168
 169     trk.GetP(fMomentum);
 170     trk.GetPosition(fPosition);
 171     trk.GetPID(fPID);
 172
 173     fChi2perNDF = trk.fChi2perNDF;
 174     fChi2MatchTrigger = trk.fChi2MatchTrigger;
 175
 176     fID = trk.fID;
 177     fLabel = trk.fLabel;
 178
 179     fITSMuonClusterMap = trk.fITSMuonClusterMap;
 180     fFilterMap = trk.fFilterMap;
 181
 182     fCharge = trk.fCharge;
 183     fType = trk.fType;
 184
 185     delete fCovMatrix;
 186     if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
 187     else fCovMatrix=NULL;
 188     fProdVertex = trk.fProdVertex;
 189
 190     SetUsedForVtxFit(trk.GetUsedForVtxFit());
 191     SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
 192   }
 193
 194   return *this;
 195 }
 196
 197 //______________________________________________________________________________
 198 Double_t AliAODTrack::M(AODTrkPID_t pid) const
 199 {
 200   // Returns the mass.
 201   // Masses for nuclei don't exist in the PDG tables, therefore they were put by hand.
 202
 203   switch (pid) {
 204
 205   case kElectron :
 206     return 0.000510999; //TDatabasePDG::Instance()->GetParticle(11/*::kElectron*/)->Mass();
 207     break;
 208
 209   case kMuon :
 210     return 0.1056584; //TDatabasePDG::Instance()->GetParticle(13/*::kMuonMinus*/)->Mass();
 211     break;
 212
 213   case kPion :
 214     return 0.13957; //TDatabasePDG::Instance()->GetParticle(211/*::kPiPlus*/)->Mass();
 215     break;
 216
 217   case kKaon :
 218     return 0.4937; //TDatabasePDG::Instance()->GetParticle(321/*::kKPlus*/)->Mass();
 219     break;
 220
 221   case kProton :
 222     return 0.9382720; //TDatabasePDG::Instance()->GetParticle(2212/*::kProton*/)->Mass();
 223     break;
 224
 225   case kDeuteron :
 226     return 1.8756; //TDatabasePDG::Instance()->GetParticle(1000010020)->Mass();
 227     break;
 228
 229   case kTriton :
 230     return 2.8089; //TDatabasePDG::Instance()->GetParticle(1000010030)->Mass();
 231     break;
 232
 233   case kHelium3 :
 234     return 2.8084; //TDatabasePDG::Instance()->GetParticle(1000020030)->Mass();
 235     break;
 236
 237   case kAlpha :
 238     return 3.7274; //TDatabasePDG::Instance()->GetParticle(1000020040)->Mass();
 239     break;
 240
 241   case kUnknown :
 242     return -999.;
 243     break;
 244
 245   default :
 246     return -999.;
 247   }
 248 }
 249
 250 //______________________________________________________________________________
 251 Double_t AliAODTrack::E(AODTrkPID_t pid) const
 252 {
 253   // Returns the energy of the particle of a given pid.
 254
 255   if (pid != kUnknown) { // particle was identified
 256     Double_t m = M(pid);
 257     return TMath::Sqrt(P()*P() + m*m);
 258   } else { // pid unknown
 259     return -999.;
 260   }
 261 }
 262
 263 //______________________________________________________________________________
 264 Double_t AliAODTrack::Y(AODTrkPID_t pid) const
 265 {
 266   // Returns the rapidity of a particle of a given pid.
 267
 268   if (pid != kUnknown) { // particle was identified
 269     Double_t e = E(pid);
 270     Double_t pz = Pz();
 271     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
 272       return 0.5*TMath::Log((e+pz)/(e-pz));
 273     } else { // energy not known or equal to pz
 274       return -999.;
 275     }
 276   } else { // pid unknown
 277     return -999.;
 278   }
 279 }
 280
 281 //______________________________________________________________________________
 282 Double_t AliAODTrack::Y(Double_t m) const
 283 {
 284   // Returns the rapidity of a particle of a given mass.
 285
 286   if (m >= 0.) { // mass makes sense
 287     Double_t e = E(m);
 288     Double_t pz = Pz();
 289     if (e>=0 && e!=pz) { // energy was positive (e.g. not -999.) and not equal to pz
 290       return 0.5*TMath::Log((e+pz)/(e-pz));
 291     } else { // energy not known or equal to pz
 292       return -999.;
 293     }
 294   } else { // pid unknown
 295     return -999.;
 296   }
 297 }
 298
 299 //______________________________________________________________________________
 300 AliAODTrack::AODTrkPID_t AliAODTrack::GetMostProbablePID() const
 301 {
 302   // Returns the most probable PID array element.
 303
 304   Int_t nPID = 10;
 305   if (fPID) {
 306     AODTrkPID_t loc = kUnknown;
 307     Double_t max = 0.;
 308     Bool_t allTheSame = kTRUE;
 309
 310     for (Int_t iPID = 0; iPID < nPID; iPID++) {
 311       if (fPID[iPID] >= max) {
 312         if (fPID[iPID] > max) {
 313           allTheSame = kFALSE;
 314           max = fPID[iPID];
 315           loc = (AODTrkPID_t)iPID;
 316         } else {
 317           allTheSame = kTRUE;
 318         }
 319       }
 320     }
 321
 322     return allTheSame ? kUnknown : loc;
 323   } else {
 324     return kUnknown;
 325   }
 326 }
 327
 328 //______________________________________________________________________________
 329 void AliAODTrack::ConvertAliPIDtoAODPID()
 330 {
 331   // Converts AliPID array.
 332   // The numbering scheme is the same for electrons, muons, pions, kaons, and protons.
 333   // Everything else has to be set to zero.
 334
 335   fPID[kDeuteron] = 0.;
 336   fPID[kTriton]   = 0.;
 337   fPID[kHelium3]  = 0.;
 338   fPID[kAlpha]    = 0.;
 339   fPID[kUnknown]  = 0.;
 340
 341   return;
 342 }
 343
 344
 345 //______________________________________________________________________________
 346 template <class T> void AliAODTrack::SetP(const T *p, const Bool_t cartesian)
 347 {
 348   // Set the momentum
 349
 350   if (p) {
 351     if (cartesian) {
 352       Double_t pt2 = p[0]*p[0] + p[1]*p[1];
 353       Double_t P = TMath::Sqrt(pt2 + p[2]*p[2]);
 354
 355       fMomentum[0] = TMath::Sqrt(pt2); // pt
 356       fMomentum[1] = (pt2 != 0.) ? TMath::Pi()+TMath::ATan2(-p[1], -p[0]) : -999; // phi
 357       fMomentum[2] = (P != 0.) ? TMath::ACos(p[2]/P) : -999.; // theta
 358     } else {
 359       fMomentum[0] = p[0];  // pt
 360       fMomentum[1] = p[1];  // phi
 361       fMomentum[2] = p[2];  // theta
 362     }
 363   } else {
 364     fMomentum[0] = -999.;
 365     fMomentum[1] = -999.;
 366     fMomentum[2] = -999.;
 367   }
 368 }
 369
 370 //______________________________________________________________________________
 371 template <class T> void AliAODTrack::SetPosition(const T *x, const Bool_t dca)
 372 {
 373   // set the position
 374
 375   if (x) {
 376     if (!dca) {
 377       ResetBit(kIsDCA);
 378
 379       fPosition[0] = x[0];
 380       fPosition[1] = x[1];
 381       fPosition[2] = x[2];
 382     } else {
 383       SetBit(kIsDCA);
 384       // don't know any better yet
 385       fPosition[0] = -999.;
 386       fPosition[1] = -999.;
 387       fPosition[2] = -999.;
 388     }
 389   } else {
 390     ResetBit(kIsDCA);
 391
 392     fPosition[0] = -999.;
 393     fPosition[1] = -999.;
 394     fPosition[2] = -999.;
 395   }
 396 }
 397
 398 //______________________________________________________________________________
 399 void AliAODTrack::SetDCA(Double_t d, Double_t z)
 400 {
 401   // set the dca
 402   fPosition[0] = d;
 403   fPosition[1] = z;
 404   fPosition[2] = 0.;
 405   SetBit(kIsDCA);
 406 }
 407
 408 //______________________________________________________________________________
 409 void AliAODTrack::Print(Option_t* /* option */) const
 410 {
 411   // prints information about AliAODTrack
 412
 413   printf("Object name: %s   Track type: %s\n", GetName(), GetTitle());
 414   printf("        px = %f\n", Px());
 415   printf("        py = %f\n", Py());
 416   printf("        pz = %f\n", Pz());
 417   printf("        pt = %f\n", Pt());
 418   printf("      1/pt = %f\n", OneOverPt());
 419   printf("     theta = %f\n", Theta());
 420   printf("       phi = %f\n", Phi());
 421   printf("  chi2/NDF = %f\n", Chi2perNDF());
 422   printf("    charge = %d\n", Charge());
 423   printf(" PID object: %p\n", PID());
 424 }
 425
 426 void AliAODTrack::SetMatchTrigger(Int_t MatchTrigger){
 427 //
 428 // Set the MUON trigger information
 429   switch(MatchTrigger){
 430     case 0: // 0 track does not match trigger
 431       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
 432       break;
 433     case 1: // 1 track match but does not pass pt cut
 434       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x40000000;
 435       break;
 436     case 2: // 2 track match Low pt cut
 437       fITSMuonClusterMap=(fITSMuonClusterMap&0x3fffffff)|0x80000000;
 438       break;
 439     case 3: // 3 track match High pt cut
 440       fITSMuonClusterMap=fITSMuonClusterMap|0xc0000000;
 441       break;
 442     default:
 443       fITSMuonClusterMap=fITSMuonClusterMap&0x3fffffff;
 444       printf("AliAODTrack::SetMatchTrigger unknown case for MatchTrigger: %d\n",MatchTrigger);
 445   }
 446 }
 447
 448 void AliAODTrack::SetHitsPatternInTrigCh(UShort_t hitsPatternInTrigCh){
 449 //
 450 // Set the MUON hit pattern (1 bit per chamber)
 451   fITSMuonClusterMap=(fITSMuonClusterMap&0xffff00ff)|(hitsPatternInTrigCh<<8);
 452 }
 453
 454 Int_t AliAODTrack::HitsMT(Int_t istation, Int_t iplane, Char_t *cathode){
 455 //
 456 // Retrieve hit information for MUON identified by  (station, plane, cathode)
 457   if(cathode){
 458     if(cathode[0]=='x'||cathode[0]=='X'){
 459       if(istation==1){
 460         if(iplane==1)
 461           return (fITSMuonClusterMap&0x8000)?1:0;
 462         else if(iplane==2)
 463           return (fITSMuonClusterMap&0x4000)?1:0;
 464         else
 465           return 0;
 466       }else if(istation==2){
 467         if(iplane==1)
 468           return (fITSMuonClusterMap&0x2000)?1:0;
 469         else if(iplane==2)
 470           return (fITSMuonClusterMap&0x1000)?1:0;
 471         else
 472           return 0;
 473       }else{
 474         return 0;
 475       }
 476     }else if(cathode[0]=='y'||cathode[0]=='Y'){
 477       if(istation==1){
 478         if(iplane==1)
 479           return (fITSMuonClusterMap&0x0800)?1:0;
 480         else if(iplane==2)
 481           return (fITSMuonClusterMap&0x0400)?1:0;
 482         else
 483           return 0;
 484       }else if(istation==2){
 485         if(iplane==1)
 486           return (fITSMuonClusterMap&0x0200)?1:0;
 487         else if(iplane==2)
 488           return (fITSMuonClusterMap&0x0100)?1:0;
 489         else
 490           return 0;
 491       }else{
 492         return 0;
 493       }
 494     }else{
 495       return 0;
 496     }
 497   }else{
 498     if(istation==1){
 499       if(iplane==1)
 500         return (HitsMT(1,1,"X")||HitsMT(1,1,"Y"))?1:0;
 501       else if(iplane==2)
 502         return (HitsMT(1,2,"X")||HitsMT(1,2,"Y"))?1:0;
 503       else
 504         return 0;
 505     }else if(istation==2){
 506       if(iplane==1)
 507         return (HitsMT(2,1,"X")||HitsMT(2,1,"Y"))?1:0;
 508       else if(iplane==2)
 509         return (HitsMT(2,2,"X")||HitsMT(2,2,"Y"))?1:0;
 510       else
 511         return 0;
 512     }else{
 513       return 0;
 514     }
 515   }
 516 }
 517
 518 Int_t AliAODTrack::HitsMuonChamber(Int_t MuonChamber){
 519 // Retrieve hit information for MUON Chamber
 520   switch(MuonChamber){
 521     case 11:
 522       return HitsMT(1,1);
 523     case 12:
 524       return HitsMT(1,2);
 525     case 13:
 526       return HitsMT(2,1);
 527     case 14:
 528       return HitsMT(2,2);
 529     default:
 530       printf("Unknown MUON chamber: %d\n",MuonChamber);
 531       return 0;
 532   }
 533 }
 534
 535
 536
 537 Bool_t AliAODTrack::PropagateTo(Double_t xk, Double_t b) {
 538   //----------------------------------------------------------------
 539   // Propagate this track to the plane X=xk (cm) in the field "b" (kG)
 540   // This is in local coordinates!!!
 541   //----------------------------------------------------------------
 542
 543   Double_t alpha = 0.;
 544   Double_t localP[3] = {Px(), Py(), Pz()}; // set global (sic!) p
 545   Global2LocalMomentum(localP, Charge(), alpha); // convert global to local momentum
 546
 547   AliAODVertex *origin = (AliAODVertex*)fProdVertex.GetObject();
 548   Double_t localX[3] = {origin->GetX(), origin->GetY(), origin->GetZ()}; // set global (sic!) location of first track point
 549   Global2LocalPosition(localX, alpha); // convert global to local position
 550
 551   Double_t &fX = localX[0];
 552
 553   Double_t dx=xk-fX;
 554   if (TMath::Abs(dx)<=kAlmost0)  return kTRUE;
 555
 556   Double_t crv=localP[0]*b*kB2C;
 557   if (TMath::Abs(b) < kAlmost0Field) crv=0.;
 558
 559   Double_t f1=localP[1], f2=f1 + crv*dx;
 560   if (TMath::Abs(f1) >= kAlmost1) return kFALSE;
 561   if (TMath::Abs(f2) >= kAlmost1) return kFALSE;
 562
 563   Double_t &fP0=localX[1], &fP1=localX[2], &fP2=localP[0], &fP3=localP[1], &fP4=localP[2];
 564   /* covariance matrix to be fixed!
 565   Double_t
 566   &fC00=fC[0],
 567   &fC10=fC[1],   &fC11=fC[2],
 568   &fC20=fC[3],   &fC21=fC[4],   &fC22=fC[5],
 569   &fC30=fC[6],   &fC31=fC[7],   &fC32=fC[8],   &fC33=fC[9],
 570   &fC40=fC[10],  &fC41=fC[11],  &fC42=fC[12],  &fC43=fC[13], &fC44=fC[14];
 571   */
 572   Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
 573
 574   fX=xk;
 575   fP0 += dx*(f1+f2)/(r1+r2);
 576   fP1 += dx*(r2 + f2*(f1+f2)/(r1+r2))*fP3;
 577   fP2 += dx*crv;
 578
 579   //f = F - 1
 580
 581   Double_t f02=    dx/(r1*r1*r1);            Double_t cc=crv/fP4;
 582   Double_t f04=0.5*dx*dx/(r1*r1*r1);         f04*=cc;
 583   Double_t f12=    dx*fP3*f1/(r1*r1*r1);
 584   Double_t f14=0.5*dx*dx*fP3*f1/(r1*r1*r1);  f14*=cc;
 585   Double_t f13=    dx/r1;
 586   Double_t f24=    dx;                       f24*=cc;
 587
 588   /* covariance matrix to be fixed!
 589   //b = C*ft
 590   Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
 591   Double_t b02=f24*fC40;
 592   Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
 593   Double_t b12=f24*fC41;
 594   Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
 595   Double_t b22=f24*fC42;
 596   Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
 597   Double_t b42=f24*fC44;
 598   Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
 599   Double_t b32=f24*fC43;
 600
 601   //a = f*b = f*C*ft
 602   Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a02=f02*b22+f04*b42;
 603   Double_t a11=f12*b21+f14*b41+f13*b31,a12=f12*b22+f14*b42+f13*b32;
 604   Double_t a22=f24*b42;
 605
 606   //F*C*Ft = C + (b + bt + a)
 607   fC00 += b00 + b00 + a00;
 608   fC10 += b10 + b01 + a01;
 609   fC20 += b20 + b02 + a02;
 610   fC30 += b30;
 611   fC40 += b40;
 612   fC11 += b11 + b11 + a11;
 613   fC21 += b21 + b12 + a12;
 614   fC31 += b31;
 615   fC41 += b41;
 616   fC22 += b22 + b22 + a22;
 617   fC32 += b32;
 618   fC42 += b42;
 619   */
 620
 621   Local2GlobalMomentum(localP, alpha); // convert local to global momentum
 622   SetP(localP);
 623
 624   return kTRUE;
 625 }