PWG3/hfe/AliHFEsecVtx.cxx

   1  /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 // Secondary vertexing construction Class
  17 //  Construct secondary vertex from Beauty hadron with electron and
  18 //  hadrons, then apply selection criteria
  19 //
  20 // Authors:
  21 //   MinJung Kweon <minjung@physi.uni-heidelberg.de>
  22 //
  23
  24 #include <TH2F.h>
  25 #include <TParticle.h>
  26
  27 #include <AliESDEvent.h>
  28 #include <AliAODEvent.h>
  29 #include <AliVTrack.h>
  30 #include <AliESDtrack.h>
  31 #include <AliAODTrack.h>
  32 #include "AliHFEsecVtx.h"
  33 #include <AliKFParticle.h>
  34 #include <AliKFVertex.h>
  35 #include <AliLog.h>
  36 #include <AliStack.h>
  37 #include <AliAODMCParticle.h>
  38 #include "AliHFEpairs.h"
  39 #include "AliHFEsecVtxs.h"
  40
  41 ClassImp(AliHFEsecVtx)
  42
  43 //_______________________________________________________________________________________________
  44 AliHFEsecVtx::AliHFEsecVtx():
  45   fESD1(0x0)
  46   ,fAOD1(0x0)
  47   ,fStack(0x0)
  48   ,fUseMCPID(kFALSE)
  49   ,fkSourceLabel()
  50   ,fNparents(0)
  51   ,fParentSelect()
  52   ,fNoOfHFEpairs(0)
  53   ,fNoOfHFEsecvtxs(0)
  54   ,fHFEpairs(0x0)
  55   ,fHFEsecvtxs(0x0)
  56   ,fMCArray(0x0)
  57   ,fPVx(0)
  58   ,fPVy(0)
  59   ,fCosPhi(-1)
  60   ,fSignedLxy(-1)
  61   ,fKFchi2(-1)
  62   ,fInvmass(-1)
  63   ,fInvmassSigma(-1)
  64   ,fKFip(0)
  65   ,fPairQA(0x0)
  66   ,fSecvtxQA(0x0)
  67   ,fSecVtxList(0x0)
  68 {
  69   //
  70   // Default constructor
  71   //
  72
  73   Init();
  74 }
  75
  76 //_______________________________________________________________________________________________
  77 AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
  78   TObject(p)
  79   ,fESD1(0x0)
  80   ,fAOD1(0x0)
  81   ,fStack(0x0)
  82   ,fUseMCPID(p.fUseMCPID)
  83   ,fkSourceLabel()
  84   ,fNparents(p.fNparents)
  85   ,fParentSelect()
  86   ,fNoOfHFEpairs(p.fNoOfHFEpairs)
  87   ,fNoOfHFEsecvtxs(p.fNoOfHFEsecvtxs)
  88   ,fHFEpairs(0x0)
  89   ,fHFEsecvtxs(0x0)
  90   ,fMCArray(0x0)
  91   ,fPVx(p.fPVx)
  92   ,fPVy(p.fPVy)
  93   ,fCosPhi(p.fCosPhi)
  94   ,fSignedLxy(p.fSignedLxy)
  95   ,fKFchi2(p.fKFchi2)
  96   ,fInvmass(p.fInvmass)
  97   ,fInvmassSigma(p.fInvmassSigma)
  98   ,fKFip(p.fKFip)
  99   ,fPairQA(0x0)
 100   ,fSecvtxQA(0x0)
 101   ,fSecVtxList(0x0)
 102 {
 103   //
 104   // Copy constructor
 105   //
 106 }
 107
 108 //_______________________________________________________________________________________________
 109 AliHFEsecVtx&
 110 AliHFEsecVtx::operator=(const AliHFEsecVtx &)
 111 {
 112   //
 113   // Assignment operator
 114   //
 115
 116   AliInfo("Not yet implemented.");
 117   return *this;
 118 }
 119
 120 //_______________________________________________________________________________________________
 121 AliHFEsecVtx::~AliHFEsecVtx()
 122 {
 123   //
 124   // Destructor
 125   //
 126
 127   //cout << "Analysis Done." << endl;
 128 }
 129
 130 //__________________________________________
 131 void AliHFEsecVtx::Init()
 132 {
 133   //
 134   // set pdg code and index
 135   //
 136
 137   fNparents = 7;
 138
 139   fParentSelect[0][0] =  411;
 140   fParentSelect[0][1] =  421;
 141   fParentSelect[0][2] =  431;
 142   fParentSelect[0][3] = 4122;
 143   fParentSelect[0][4] = 4132;
 144   fParentSelect[0][5] = 4232;
 145   fParentSelect[0][6] = 4332;
 146
 147   fParentSelect[1][0] =  511;
 148   fParentSelect[1][1] =  521;
 149   fParentSelect[1][2] =  531;
 150   fParentSelect[1][3] = 5122;
 151   fParentSelect[1][4] = 5132;
 152   fParentSelect[1][5] = 5232;
 153   fParentSelect[1][6] = 5332;
 154 }
 155
 156 //_______________________________________________________________________________________________
 157 void AliHFEsecVtx::CreateHistograms(TList *qaList)
 158 {
 159   //
 160   // create histograms
 161   //
 162
 163   fSecVtxList = new TList;
 164   fSecVtxList->SetName("SecVtx");
 165
 166   MakeContainer();
 167   /*
 168   fkSourceLabel[kAll]="all";
 169   fkSourceLabel[kDirectCharm]="directCharm";
 170   fkSourceLabel[kDirectBeauty]="directBeauty";
 171   fkSourceLabel[kBeautyCharm]="beauty2charm";
 172   fkSourceLabel[kGamma]="gamma";
 173   fkSourceLabel[kPi0]="pi0";
 174   fkSourceLabel[kElse]="others";
 175   fkSourceLabel[kBeautyGamma]="beauty22gamma";
 176   fkSourceLabel[kBeautyPi0]="beauty22pi0";
 177   fkSourceLabel[kBeautyElse]="beauty22others";
 178
 179   TString hname;
 180   TString hnopt="secvtx_";
 181   for (Int_t isource = 0; isource < 10; isource++ ){
 182   }
 183   */
 184
 185   qaList->AddLast(fSecVtxList);
 186 }
 187
 188 //_______________________________________________________________________________________________
 189 void AliHFEsecVtx::GetPrimaryCondition()
 190 {
 191   //
 192   // get primary characteristics and set
 193   //
 194
 195   if (fESD1) {
 196     AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
 197     if( primVtxCopy.GetNDF() <1 ) return;
 198     fPVx = primVtxCopy.GetX();
 199     fPVx = primVtxCopy.GetY();
 200   }
 201   else if(fAOD1) {
 202     AliKFVertex primVtxCopy(*(fAOD1->GetPrimaryVertex()));
 203     if( primVtxCopy.GetNDF() <1 ) return;
 204     fPVx = primVtxCopy.GetX();
 205     fPVx = primVtxCopy.GetY();
 206   }
 207 }
 208
 209 //_______________________________________________________________________________________________
 210 void AliHFEsecVtx::PairAnalysis(AliVTrack* track1, AliVTrack* track2, Int_t index2)
 211 {
 212   //
 213   // calculate e-h pair characteristics and tag pair
 214   //
 215
 216   // get KF particle input pid
 217   Int_t pdg1 = GetPDG(track1);
 218   Int_t pdg2 = GetPDG(track2);
 219
 220   if(pdg1==-1 || pdg2==-1) {
 221     //printf("out if considered pid range \n");
 222     return;
 223   }
 224
 225   // create KF particle of pair
 226   if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
 227   else AliKFParticle::SetField(fESD1->GetMagneticField());
 228   AliKFParticle kfTrack1(*track1, pdg1);
 229   AliKFParticle kfTrack2(*track2, pdg2);
 230
 231   AliKFParticle kfSecondary(kfTrack1,kfTrack2);
 232
 233   //secondary vertex point from kf particle
 234   Double_t kfx = kfSecondary.GetX();
 235   Double_t kfy = kfSecondary.GetY();
 236   //Double_t kfz = kfSecondary.GetZ();
 237
 238   //momentum at the decay point from kf particle
 239   Double_t kfpx = kfSecondary.GetPx();
 240   Double_t kfpy = kfSecondary.GetPy();
 241   //Double_t kfpz = kfSecondary.GetPz();
 242
 243   Double_t dx = kfx-fPVx;
 244   Double_t dy = kfy-fPVy;
 245
 246   // discriminating variables
 247   // invariant mass of the KF particle
 248   Double_t invmass = -1;
 249   Double_t invmassSigma = -1;
 250   kfSecondary.GetMass(invmass,invmassSigma);
 251
 252   // chi2 of the KF particle
 253   Double_t kfchi2 = -1;
 254   if(kfSecondary.GetNDF()>0) kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 255
 256   // opening angle between two particles in XY plane
 257   Double_t phi = kfTrack1.GetAngleXY(kfTrack2);
 258   Double_t cosphi = TMath::Cos(phi);
 259
 260   // projection of kf vertex vector to the kf momentum direction
 261   Double_t signedLxy=-999.;
 262   Double_t psqr = kfpx*kfpx+kfpy*kfpy;
 263   /*
 264   //[the other way to think about]
 265   if(psqr>0) {
 266     if((dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr)>0) signedLxy = TMath::Sqrt(dx*dx+dy*dy);
 267     if((dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr)<0) signedLxy = -1*TMath::Sqrt(dx*dx+dy*dy);
 268   }*/
 269   if(psqr>0) signedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
 270
 271   // DCA from primary to e-h KF particle (impact parameter of KF particle)
 272   Double_t vtx[2]={fPVx, fPVy};
 273   Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
 274
 275   /* //later consider the below
 276   Float_t dcaR1=-999., dcaR2=-999.;
 277   Float_t dcaZ1=-999., dcaZ2=-999.;
 278
 279   if(IsAODanalysis()){
 280     Double_t trkIpPar1[2];
 281     Double_t trkIpCov1[3];
 282     Double_t trkIpPar2[2];
 283     Double_t trkIpCov2[3];
 284
 285     AliESDtrack esdTrk1(track1);
 286     AliESDtrack esdTrk2(track2);
 287     esdTrk1.PropagateToDCA(fAOD1->GetPrimaryVertex(),0,10000,trkIpPar1,trkIpCov1);
 288     esdTrk2.PropagateToDCA(fAOD1->GetPrimaryVertex(),0,10000,trkIpPar2,trkIpCov2);
 289
 290     dcaR1=trkIpPar1[0];
 291     dcaZ1=trkIpPar1[1];
 292     dcaR2=trkIpPar2[0];
 293     dcaZ2=trkIpPar2[1];
 294   }
 295   else {
 296     ((AliESDtrack*)track1)->GetImpactParameters(dcaR1,dcaZ1);
 297     ((AliESDtrack*)track2)->GetImpactParameters(dcaR2,dcaZ2);
 298   }*/
 299
 300   Int_t paircode = -1;
 301   if (HasMCData()) paircode = GetPairCode(track1,track2);
 302
 303   AliHFEpairs hfepair;
 304   hfepair.SetTrkLabel(index2);
 305   hfepair.SetInvmass(invmass);
 306   hfepair.SetKFChi2(kfchi2);
 307   hfepair.SetOpenangle(phi);
 308   hfepair.SetCosOpenangle(cosphi);
 309   hfepair.SetSignedLxy(signedLxy);
 310   hfepair.SetKFIP(kfip);
 311   hfepair.SetPairCode(paircode);
 312   AddHFEpairToArray(&hfepair);
 313   fNoOfHFEpairs++;
 314
 315   // fill into container for later QA
 316   Double_t dataE[6];
 317   dataE[0]=invmass;
 318   dataE[1]=kfchi2;
 319   dataE[2]=phi;
 320   dataE[3]=signedLxy;
 321   dataE[4]=kfip;
 322   dataE[5]=paircode;
 323   fPairQA->Fill(dataE);
 324 }
 325
 326 //_______________________________________________________________________________________________
 327 void AliHFEsecVtx::RunSECVTX(AliVTrack *track)
 328 {
 329   //
 330   // run secondary vertexing algorithm and do tagging
 331   //
 332
 333   //printf("number of considered pairs= %d\n",HFEpairs()->GetEntriesFast());
 334   FindSECVTXCandid(track);
 335 }
 336
 337 //_______________________________________________________________________________________________
 338 void AliHFEsecVtx::FindSECVTXCandid(AliVTrack *track)
 339 {
 340   //
 341   // find secondary vertex candidate and store them into container
 342   //
 343
 344   AliVTrack *htrack[20];
 345   Int_t htracklabel[20];
 346   Double_t vtxchi2cut=3.; // testing cut
 347   Double_t dataE[6]={-999.,-999.,-999.,-999.,-1.,0};
 348   if (HFEpairs()->GetEntriesFast()>20){
 349     AliDebug(3, "number of paired hadron is over maximum(20)");
 350     return;
 351   }
 352
 353   // get paired track objects
 354   AliHFEpairs *pair=0x0;
 355   if (IsAODanalysis()){
 356     for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
 357        pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
 358        htracklabel[ip] = pair->GetTrkLabel();
 359        htrack[ip] = fAOD1->GetTrack(pair->GetTrkLabel());
 360     }
 361   }
 362   else{
 363     for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
 364        pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
 365        htracklabel[ip] = pair->GetTrkLabel();
 366        htrack[ip] = fESD1->GetTrack(pair->GetTrkLabel());
 367     }
 368   }
 369   // in case there is only one paired track with the electron, put pair characteristics into secvtx container
 370   // for the moment, I only apply pair vertex chi2 cut
 371   if (HFEpairs()->GetEntriesFast() == 1){
 372     if (pair->GetKFChi2()<vtxchi2cut) { // you can also put single track cut
 373       AliHFEsecVtxs hfesecvtx;
 374       hfesecvtx.SetTrkLabel1(pair->GetTrkLabel());
 375       hfesecvtx.SetTrkLabel2(-999);
 376       hfesecvtx.SetInvmass(pair->GetInvmass());
 377       hfesecvtx.SetKFChi2(pair->GetKFChi2());
 378       hfesecvtx.SetSignedLxy(pair->GetSignedLxy());
 379       hfesecvtx.SetKFIP(pair->GetKFIP());
 380       AddHFEsecvtxToArray(&hfesecvtx);
 381       fNoOfHFEsecvtxs++;
 382
 383       dataE[0]=pair->GetInvmass();
 384       dataE[1]=pair->GetKFChi2();
 385       dataE[2]=pair->GetSignedLxy();
 386       dataE[3]=pair->GetKFIP();
 387       if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
 388       dataE[5]=2;
 389       fSecvtxQA->Fill(dataE);
 390     }
 391     return;
 392   }
 393
 394   // in case there are multiple paired track with the electron, calculate secvtx characteristics
 395   // put the secvtx characteristics into container if it passes cuts
 396   for (int i=0; i<HFEpairs()->GetEntriesFast()-1; i++){
 397      for (int j=i+1; j<HFEpairs()->GetEntriesFast(); j++){
 398         CalcSECVTXProperty(track, htrack[i], htrack[j]);
 399         if (fKFchi2<vtxchi2cut) {
 400           AliHFEsecVtxs hfesecvtx;
 401           hfesecvtx.SetTrkLabel1(htracklabel[i]);
 402           hfesecvtx.SetTrkLabel2(htracklabel[j]);
 403           hfesecvtx.SetKFChi2(fKFchi2);
 404           hfesecvtx.SetInvmass(fInvmass);
 405           hfesecvtx.SetSignedLxy(fSignedLxy);
 406           hfesecvtx.SetKFIP(fKFip);
 407           AddHFEsecvtxToArray(&hfesecvtx);
 408           fNoOfHFEsecvtxs++;
 409
 410           dataE[0]=pair->GetInvmass();
 411           dataE[1]=pair->GetKFChi2();
 412           dataE[2]=pair->GetSignedLxy();
 413           dataE[3]=pair->GetKFIP();
 414           if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
 415           dataE[5]=3;
 416           fSecvtxQA->Fill(dataE);
 417         }
 418      }
 419   }
 420 }
 421
 422 //_______________________________________________________________________________________________
 423 void AliHFEsecVtx::CalcSECVTXProperty(AliVTrack* track1, AliVTrack* track2, AliVTrack* track3)
 424 {
 425   //
 426   // calculate secondary vertex properties
 427   //
 428
 429   // get KF particle input pid
 430   Int_t pdg1 = GetPDG(track1);
 431   Int_t pdg2 = GetPDG(track2);
 432   Int_t pdg3 = GetPDG(track3);
 433
 434   if(pdg1==-1 || pdg2==-1 || pdg3==-1) {
 435     //printf("out if considered pid range \n");
 436     return;
 437   }
 438
 439   // create KF particle of pair
 440   if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
 441   else AliKFParticle::SetField(fESD1->GetMagneticField());
 442   AliKFParticle kfTrack1(*track1, pdg1);
 443   AliKFParticle kfTrack2(*track2, pdg2);
 444   AliKFParticle kfTrack3(*track3, pdg3);
 445
 446   AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
 447
 448   //secondary vertex point from kf particle
 449   Double_t kfx = kfSecondary.GetX();
 450   Double_t kfy = kfSecondary.GetY();
 451   //Double_t kfz = kfSecondary.GetZ();
 452
 453   //momentum at the decay point from kf particle
 454   Double_t kfpx = kfSecondary.GetPx();
 455   Double_t kfpy = kfSecondary.GetPy();
 456   //Double_t kfpz = kfSecondary.GetPz();
 457
 458   Double_t dx = kfx-fPVx;
 459   Double_t dy = kfy-fPVy;
 460
 461   // discriminating variables ----------------------------------------------------------
 462
 463   if(kfSecondary.GetNDF()>0) fKFchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 464
 465   // invariant mass of the KF particle
 466   kfSecondary.GetMass(fInvmass,fInvmassSigma);
 467
 468   // DCA from primary to e-h KF particle (impact parameter of KF particle)
 469   Double_t vtx[2]={fPVx, fPVy};
 470   fKFip = kfSecondary.GetDistanceFromVertexXY(vtx);
 471
 472   // projection of kf vertex vector to the kf momentum direction
 473   Double_t psqr = kfpx*kfpx+kfpy*kfpy;
 474   if(psqr>0) fSignedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
 475 }
 476
 477 //_______________________________________________________________________________________________
 478 Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track)
 479 {
 480   //
 481   // return mc pid
 482   //
 483
 484   Int_t label = TMath::Abs(track->GetLabel());
 485   TParticle* mcpart = fStack->Particle(label);
 486   if ( !mcpart ) return 0;
 487   Int_t pdgCode = mcpart->GetPdgCode();
 488
 489   return pdgCode;
 490 }
 491
 492 //_______________________________________________________________________________________________
 493 Int_t AliHFEsecVtx::GetPairOriginESD(AliESDtrack* trk1, AliESDtrack* trk2)
 494 {
 495   //
 496   // return pdg code of the origin(source) of the pair
 497   //
 498   //
 499   // ---*---*---*-----ancester A----- track1
 500   //                        |____*______
 501   //                             |______ track2
 502   // => if they originated from same ancester,
 503   //    then return "the absolute value of pdg code of ancester A"
 504   //
 505   // ---*---*---B hadron-----ancester A----- track1
 506   //                               |____*______
 507   //                                    |______ track2
 508   // => if they originated from same ancester, and this ancester originally comes from B hadrons
 509   //    then return -1*"the absolute value of pdg code of ancester A"
 510   //
 511   // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
 512   //
 513
 514   if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
 515   TParticle* part1 = fStack->Particle(trk1->GetLabel());
 516   TParticle* part2 = fStack->Particle(trk2->GetLabel());
 517   TParticle* part2cp = part2;
 518   if (!(part1) || !(part2)) return 0;
 519
 520   Int_t srcpdg = 0;
 521
 522   //if the two tracks' mother's label is same, get the mother info
 523   //in case of charm, check if it originated from beauty
 524   for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
 525      Int_t label1 = part1->GetFirstMother();
 526      if (label1 < 0) return 0;
 527
 528      for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
 529         Int_t label2 = part2->GetFirstMother();
 530         if (label2 < 0) break;
 531
 532         if (label1 == label2){ //check if two tracks are originated from same mother
 533           TParticle* commonmom = fStack->Particle(label2);
 534           srcpdg = abs(commonmom->GetPdgCode());
 535
 536           //check ancester to see if it is originally from beauty
 537           for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
 538              Int_t ancesterlabel = commonmom->GetFirstMother();
 539              if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
 540
 541              TParticle* commonancester = fStack->Particle(ancesterlabel);
 542              Int_t ancesterpdg = abs(commonancester->GetPdgCode());
 543
 544              for (Int_t l=0; l<fNparents; l++){
 545                 if (abs(ancesterpdg)==fParentSelect[1][l]){
 546                   srcpdg = -1*srcpdg; //multiply -1 for hadron from bottom
 547                   return srcpdg;
 548                 }
 549              }
 550              commonmom = commonancester;
 551           }
 552         }
 553         part2 = fStack->Particle(label2); //if their mother is different, go to earlier generation of 2nd particle
 554         if (!(part2)) break;
 555      }
 556      part1 = fStack->Particle(label1); //if their mother is different, go to earlier generation of 1st particle
 557      part2 = part2cp;
 558      if (!(part1)) return 0;
 559   }
 560
 561   return srcpdg;
 562 }
 563
 564 //_______________________________________________________________________________________________
 565 Int_t AliHFEsecVtx::GetPairOriginAOD(AliAODTrack* trk1, AliAODTrack* trk2)
 566 {
 567
 568   //
 569   // return pdg code of the origin(source) of the pair
 570   //
 571   //
 572   // ---*---*---*-----ancester A----- track1
 573   //                        |____*______
 574   //                             |______ track2
 575   // => if they originated from same ancester,
 576   //    then return "the absolute value of pdg code of ancester A"
 577   //
 578   // ---*---*---B hadron-----ancester A----- track1
 579   //                               |____*______
 580   //                                    |______ track2
 581   // => if they originated from same ancester, and this ancester originally comes from B hadrons
 582   //    then return -1*"the absolute value of pdg code of ancester A"
 583   //
 584   // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
 585   //
 586
 587   if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
 588   AliAODMCParticle *part1 = (AliAODMCParticle*)fMCArray->At(trk1->GetLabel());
 589   AliAODMCParticle *part2 = (AliAODMCParticle*)fMCArray->At(trk2->GetLabel());
 590   AliAODMCParticle *part2cp = part2;
 591   if (!(part1) || !(part2)) return 0;
 592
 593   Int_t srcpdg = 0;
 594
 595   //if the two tracks' mother's label is same, get the mother info
 596   //in case of charm, check if it originated from beauty
 597   for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
 598      Int_t label1 = part1->GetMother();
 599      if (label1 < 0) return 0;
 600
 601      for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
 602         Int_t label2 = part2->GetMother();
 603         if (label2 < 0) return 0;
 604
 605         if (label1 == label2){ //check if two tracks are originated from same mother
 606           AliAODMCParticle *commonmom = (AliAODMCParticle*)fMCArray->At(label1);
 607           srcpdg = abs(commonmom->GetPdgCode());
 608
 609           //check ancester to see if it is originally from beauty
 610           for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
 611              Int_t ancesterlabel = commonmom->GetMother();
 612              if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
 613
 614              AliAODMCParticle *commonancester = (AliAODMCParticle*)fMCArray->At(ancesterlabel);
 615              Int_t ancesterpdg = abs(commonancester->GetPdgCode());
 616
 617              for (Int_t l=0; l<fNparents; l++){
 618                 if (abs(ancesterpdg)==fParentSelect[1][l]){
 619                   srcpdg = -1*srcpdg; //multiply -1 for charm from bottom
 620                   return srcpdg;
 621                 }
 622              }
 623              commonmom = commonancester;
 624           }
 625         }
 626         part2 = (AliAODMCParticle*)fMCArray->At(label2); //if their mother is different, go to earlier generation of 2nd particle
 627         if (!(part2)) break;
 628      }
 629      part1 = (AliAODMCParticle*)fMCArray->At(label1); //if their mother is different, go to earlier generation of 2nd particle
 630      part2 = part2cp;
 631      if (!(part1)) return 0;
 632   }
 633
 634   return srcpdg;
 635 }
 636
 637 //_______________________________________________________________________________________________
 638 Int_t AliHFEsecVtx::GetPairCode(const AliVTrack* const trk1, const AliVTrack* const trk2)
 639 {
 640   //
 641   // return pair code which is predefinded as:
 642   //  kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0,
 643   //  kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
 644   //
 645
 646   Int_t srcpdg = -1;
 647   Int_t srccode = kElse;
 648
 649   if(IsAODanalysis()) srcpdg = GetPairOriginAOD((AliAODTrack*)trk1,(AliAODTrack*)trk2);
 650   else srcpdg = GetPairOriginESD((AliESDtrack*)trk1,(AliESDtrack*)trk2);
 651
 652   if (srcpdg < 0) srccode = kBeautyElse;
 653   for (Int_t i=0; i<fNparents; i++){
 654     if (abs(srcpdg)==fParentSelect[0][i]) {
 655       if (srcpdg>0) srccode = kDirectCharm;
 656       if (srcpdg<0) srccode = kBeautyCharm;
 657     }
 658     if (abs(srcpdg)==fParentSelect[1][i]) {
 659       if (srcpdg>0) srccode = kDirectBeauty;
 660       if (srcpdg<0)  return kElse;
 661     }
 662   }
 663   if (srcpdg == 22) srccode = kGamma;
 664   if (srcpdg == -22) srccode = kBeautyGamma;
 665   if (srcpdg == 111) srccode = kPi0;
 666   if (srcpdg == -111) srccode = kBeautyPi0;
 667
 668   return srccode;
 669 }
 670
 671 //_______________________________________________________________________________________________
 672 Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack)
 673 {
 674   //
 675   // return decay electron's origin
 676   //
 677
 678   if (iTrack < 0) {
 679     AliDebug(1, "Stack label is negative, return\n");
 680     return -1;
 681   }
 682
 683   TParticle* mcpart = fStack->Particle(iTrack);
 684
 685 //  if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !
 686
 687   Int_t iLabel = mcpart->GetFirstMother();
 688   if (iLabel<0){
 689     AliDebug(1, "Stack label is negative, return\n");
 690     return -1;
 691   }
 692
 693   Int_t origin = -1;
 694   Bool_t isFinalOpenCharm = kFALSE;
 695
 696   TParticle *partMother = fStack->Particle(iLabel);
 697   Int_t maPdgcode = partMother->GetPdgCode();
 698
 699   // if the mother is charmed hadron
 700   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 701
 702     for (Int_t i=0; i<fNparents; i++){
 703       if (abs(maPdgcode)==fParentSelect[0][i]){
 704         isFinalOpenCharm = kTRUE;
 705       }
 706     }
 707     if (!isFinalOpenCharm) return -1;
 708
 709     // iterate until you find B hadron as a mother or become top ancester
 710     for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
 711
 712       Int_t jLabel = partMother->GetFirstMother();
 713       if (jLabel == -1){
 714         origin = kDirectCharm;
 715         return origin;
 716       }
 717       if (jLabel < 0){ // safety protection
 718         AliDebug(1, "Stack label is negative, return\n");
 719         return -1;
 720       }
 721
 722       // if there is an ancester
 723       TParticle* grandMa = fStack->Particle(jLabel);
 724       Int_t grandMaPDG = grandMa->GetPdgCode();
 725
 726       for (Int_t j=0; j<fNparents; j++){
 727         if (abs(grandMaPDG)==fParentSelect[1][j]){
 728           origin = kBeautyCharm;
 729           return origin;
 730         }
 731       }
 732
 733       partMother = grandMa;
 734     } // end of iteration
 735   } // end of if
 736   else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
 737     for (Int_t i=0; i<fNparents; i++){
 738       if (abs(maPdgcode)==fParentSelect[1][i]){
 739         origin = kDirectBeauty;
 740         return origin;
 741       }
 742     }
 743   } // end of if
 744
 745   //============ gamma ================
 746   else if ( abs(maPdgcode) == 22 ) {
 747     origin = kGamma;
 748
 749     // iterate until you find B hadron as a mother or become top ancester
 750     for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
 751
 752       Int_t jLabel = partMother->GetFirstMother();
 753       if (jLabel == -1){
 754         origin = kGamma;
 755         return origin;
 756       }
 757       if (jLabel < 0){ // safety protection
 758         AliDebug(1, "Stack label is negative, return\n");
 759         return -1;
 760       }
 761
 762       // if there is an ancester
 763       TParticle* grandMa = fStack->Particle(jLabel);
 764       Int_t grandMaPDG = grandMa->GetPdgCode();
 765
 766       for (Int_t j=0; j<fNparents; j++){
 767         if (abs(grandMaPDG)==fParentSelect[1][j]){
 768           origin = kBeautyGamma;
 769           return origin;
 770         }
 771       }
 772
 773       partMother = grandMa;
 774     } // end of iteration
 775
 776     return origin;
 777   } // end of if
 778
 779   //============ pi0 ================
 780   else if ( abs(maPdgcode) == 111 ) {
 781     origin = kPi0;
 782
 783     // iterate until you find B hadron as a mother or become top ancester
 784     for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
 785
 786       Int_t jLabel = partMother->GetFirstMother();
 787       if (jLabel == -1){
 788         origin = kPi0;
 789         return origin;
 790       }
 791       if (jLabel < 0){ // safety protection
 792         AliDebug(1, "Stack label is negative, return\n");
 793         return -1;
 794       }
 795
 796       // if there is an ancester
 797       TParticle* grandMa = fStack->Particle(jLabel);
 798       Int_t grandMaPDG = grandMa->GetPdgCode();
 799
 800       for (Int_t j=0; j<fNparents; j++){
 801         if (abs(grandMaPDG)==fParentSelect[1][j]){
 802           origin = kBeautyPi0;
 803           return origin;
 804         }
 805       }
 806
 807       partMother = grandMa;
 808     } // end of iteration
 809
 810     return origin;
 811   } // end of if
 812
 813   else {
 814     origin = kElse;
 815     // iterate until you find B hadron as a mother or become top ancester
 816     for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
 817
 818       Int_t jLabel = partMother->GetFirstMother();
 819       if (jLabel == -1){
 820         origin = kElse;
 821         return origin;
 822       }
 823       if (jLabel < 0){ // safety protection
 824         AliDebug(1, "Stack label is negative, return\n");
 825         return -1;
 826       }
 827
 828       // if there is an ancester
 829       TParticle* grandMa = fStack->Particle(jLabel);
 830       Int_t grandMaPDG = grandMa->GetPdgCode();
 831
 832       for (Int_t j=0; j<fNparents; j++){
 833         if (abs(grandMaPDG)==fParentSelect[1][j]){
 834           origin = kBeautyElse;
 835           return origin;
 836         }
 837       }
 838
 839       partMother = grandMa;
 840     } // end of iteration
 841   }
 842
 843   return origin;
 844 }
 845
 846 //_______________________________________________________________________________________________
 847 Int_t AliHFEsecVtx::GetPDG(AliVTrack *track)
 848 {
 849   //
 850   // get KF particle input pdg for mass hypothesis
 851   //
 852
 853   Int_t pdgCode=-1;
 854
 855   if (fUseMCPID && HasMCData()){
 856     pdgCode = GetMCPDG(track);
 857   }
 858   else if(fESD1){
 859     Int_t pid=0;
 860     Double_t prob;
 861     GetESDPID((AliESDtrack*)track, pid, prob);
 862     switch(pid){
 863       case 0:  pdgCode = 11; break;
 864       case 1:  pdgCode = 13; break;
 865       case 2:  pdgCode = 211; break;
 866       case 3:  pdgCode = 321; break;
 867       case 4:  pdgCode = 2212; break;
 868       default: pdgCode = -1;
 869     }
 870   }
 871   else if(fAOD1){
 872     Int_t pid = ((AliAODTrack*)track)->GetMostProbablePID();
 873     switch(pid){
 874       case 0:  pdgCode = 11; break;
 875       case 1:  pdgCode = 13; break;
 876       case 2:  pdgCode = 211; break;
 877       case 3:  pdgCode = 321; break;
 878       case 4:  pdgCode = 2212; break;
 879       default: pdgCode = -1;
 880     }
 881   }
 882
 883   return pdgCode;
 884 }
 885
 886 //_______________________________________________________________________________________________
 887 Int_t AliHFEsecVtx::GetMCPDG(AliVTrack *track)
 888 {
 889   //
 890   // return mc pdg code
 891   //
 892
 893   Int_t label = TMath::Abs(track->GetLabel());
 894   Int_t pdgCode;
 895
 896   if (IsAODanalysis()) {
 897     AliAODMCParticle *mcpart = (AliAODMCParticle*)fMCArray->At(label);
 898     if ( !mcpart ) return 0;
 899       pdgCode = mcpart->GetPdgCode();
 900   }
 901   else {
 902     TParticle* mcpart = fStack->Particle(label);
 903     if ( !mcpart ) return 0;
 904       pdgCode = mcpart->GetPdgCode();
 905   }
 906
 907     return pdgCode;
 908 }
 909
 910 //______________________________________________________________________________
 911 void AliHFEsecVtx::GetESDPID(AliESDtrack *track, Int_t &recpid, Double_t &recprob)
 912 {
 913   //
 914   // calculate likehood for esd pid
 915   //
 916
 917   recpid = -1;
 918   recprob = -1;
 919
 920   Int_t ipid=-1;
 921   Double_t max=0.;
 922
 923   Double_t probability[5];
 924
 925   // get probability of the diffenrent particle types
 926   track->GetESDpid(probability);
 927
 928   // find most probable particle in ESD pid
 929   // 0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
 930   ipid = TMath::LocMax(5,probability);
 931   max = TMath::MaxElement(5,probability);
 932
 933   recpid = ipid;
 934   recprob = max;
 935 }
 936
 937 //_____________________________________________________________________________
 938 void AliHFEsecVtx::AddHFEpairToArray(const AliHFEpairs* const pair)
 939 {
 940   //
 941   // Add a HFE pair to the array
 942   //
 943
 944   Int_t n = HFEpairs()->GetEntriesFast();
 945   if(n!=fNoOfHFEpairs)AliError(Form("fNoOfHFEpairs != HFEpairs()->GetEntriesFast %i != %i \n", fNoOfHFEpairs, n));
 946   new((*HFEpairs())[n]) AliHFEpairs(*pair);
 947 }
 948
 949 //_____________________________________________________________________________
 950 TClonesArray *AliHFEsecVtx::HFEpairs()
 951 {
 952   //
 953   // Returns the list of HFE pairs
 954   //
 955
 956   if (!fHFEpairs) {
 957       fHFEpairs = new TClonesArray("AliHFEpairs", 200);
 958   }
 959   return fHFEpairs;
 960 }
 961
 962 //_____________________________________________________________________________
 963 void AliHFEsecVtx::DeleteHFEpairs()
 964 {
 965   //
 966   // Delete the list of HFE pairs
 967   //
 968
 969   if (fHFEpairs) {
 970     fHFEpairs->Delete();
 971     //delete fHFEpairs;
 972   }
 973 }
 974
 975 //_____________________________________________________________________________
 976 void AliHFEsecVtx::InitHFEpairs()
 977 {
 978   //
 979   // Initialization should be done before make all possible pairs for a given electron candidate
 980   //
 981
 982   fNoOfHFEpairs = 0;
 983 }
 984
 985 //_____________________________________________________________________________
 986 void AliHFEsecVtx::AddHFEsecvtxToArray(const AliHFEsecVtxs* const secvtx)
 987 {
 988   //
 989   // Add a HFE secondary vertex to the array
 990   //
 991
 992   Int_t n = HFEsecvtxs()->GetEntriesFast();
 993   if(n!=fNoOfHFEsecvtxs)AliError(Form("fNoOfHFEsecvtxs != HFEsecvtxs()->GetEntriesFast %i != %i \n", fNoOfHFEsecvtxs, n));
 994   new((*HFEsecvtxs())[n]) AliHFEsecVtxs(*secvtx);
 995 }
 996
 997 //_____________________________________________________________________________
 998 TClonesArray *AliHFEsecVtx::HFEsecvtxs()
 999 {
1000   //
1001   // Returns the list of HFE secvtx
1002   //
1003
1004   if (!fHFEsecvtxs) {
1005       fHFEsecvtxs = new TClonesArray("AliHFEsecVtxs", 200);
1006   }
1007   return fHFEsecvtxs;
1008 }
1009
1010 //_____________________________________________________________________________
1011 void AliHFEsecVtx::DeleteHFEsecvtxs()
1012 {
1013   //
1014   // Delete the list of HFE pairs
1015   //
1016
1017   if (fHFEsecvtxs) {
1018     fHFEsecvtxs->Delete();
1019     //delete fHFEsecvtx;
1020   }
1021 }
1022
1023 //_____________________________________________________________________________
1024 void AliHFEsecVtx::InitHFEsecvtxs()
1025 {
1026   //
1027   // Initialization should be done
1028   //
1029
1030   fNoOfHFEsecvtxs = 0;
1031 }
1032
1033 //_______________________________________________________________________________________________
1034 Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track) const
1035 {
1036   //if (track->Pt() < 1.0) return kFALSE;
1037   //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
1038   //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
1039   //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1040   if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
1041   //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
1042
1043 /*
1044   Float_t dcaR=-1;
1045   Float_t dcaZ=-1;
1046   track->GetImpactParameters(dcaR,dcaZ);
1047   if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
1048   if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
1049 */
1050   return kTRUE;
1051 }
1052 //____________________________________________________________
1053 void AliHFEsecVtx::MakeContainer(){
1054
1055   //
1056   // make container
1057   //
1058
1059   const Int_t nDimPair=6;
1060   Int_t nBinPair[nDimPair] = {200, 500, 314, 2000, 2000, 11};
1061   const Double_t kInvmassmin = 0., kInvmassmax = 20.;
1062   const Double_t kKFChi2min = 0, kKFChi2max= 50;
1063   const Double_t kOpenanglemin = 0, kOpenanglemax = 3.14;
1064   const Double_t kSignedLxymin = -10, kSignedLxymax= 10;
1065   const Double_t kKFIPmin = -10, kKFIPmax= 10;
1066   const Double_t kPairCodemin = -1, kPairCodemax= 10;
1067
1068   Double_t* binEdgesPair[nDimPair];
1069   for(Int_t ivar = 0; ivar < nDimPair; ivar++)
1070     binEdgesPair[ivar] = new Double_t[nBinPair[ivar] + 1];
1071
1072   for(Int_t i=0; i<=nBinPair[0]; i++) binEdgesPair[0][i]=(Double_t)kInvmassmin + (kInvmassmax - kInvmassmin)/nBinPair[0]*(Double_t)i;
1073   for(Int_t i=0; i<=nBinPair[1]; i++) binEdgesPair[1][i]=(Double_t)kKFChi2min + (kKFChi2max - kKFChi2min)/nBinPair[1]*(Double_t)i;
1074   for(Int_t i=0; i<=nBinPair[2]; i++) binEdgesPair[2][i]=(Double_t)kOpenanglemin + (kOpenanglemax - kOpenanglemin)/nBinPair[2]*(Double_t)i;
1075   for(Int_t i=0; i<=nBinPair[3]; i++) binEdgesPair[3][i]=(Double_t)kSignedLxymin + (kSignedLxymax - kSignedLxymin)/nBinPair[3]*(Double_t)i;
1076   for(Int_t i=0; i<=nBinPair[4]; i++) binEdgesPair[4][i]=(Double_t)kKFIPmin + (kKFIPmax - kKFIPmin)/nBinPair[4]*(Double_t)i;
1077   for(Int_t i=0; i<=nBinPair[5]; i++) binEdgesPair[5][i]=(Double_t)kPairCodemin + (kPairCodemax - kPairCodemin)/nBinPair[5]*(Double_t)i;
1078
1079   fPairQA = new THnSparseF("pairQA", "QA for Pair; invmass[GeV/c^2]; KF chi2; opening angle; signed Lxy; KF ip; pair code ", nDimPair, nBinPair);
1080   for(Int_t idim = 0; idim < nDimPair; idim++){
1081     fPairQA->SetBinEdges(idim, binEdgesPair[idim]);
1082   }
1083
1084   fSecVtxList->AddAt(fPairQA,0);
1085
1086   const Int_t nDimSecvtx=6;
1087   Double_t* binEdgesSecvtx[nDimSecvtx];
1088   Int_t nBinSecvtx[nDimSecvtx] = {200, 500, 2000, 2000, 11, 3};
1089   const Double_t kNtrksmin = 0, kNtrksmax= 3;
1090   for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
1091     binEdgesSecvtx[ivar] = new Double_t[nBinSecvtx[ivar] + 1];
1092
1093   for(Int_t i=0; i<=nBinSecvtx[0]; i++) binEdgesSecvtx[0][i]=binEdgesPair[0][i];
1094   for(Int_t i=0; i<=nBinSecvtx[1]; i++) binEdgesSecvtx[1][i]=binEdgesPair[1][i];
1095   for(Int_t i=0; i<=nBinSecvtx[2]; i++) binEdgesSecvtx[2][i]=binEdgesPair[3][i];
1096   for(Int_t i=0; i<=nBinSecvtx[3]; i++) binEdgesSecvtx[3][i]=binEdgesPair[4][i];
1097   for(Int_t i=0; i<=nBinSecvtx[4]; i++) binEdgesSecvtx[4][i]=binEdgesPair[5][i];
1098   for(Int_t i=0; i<=nBinSecvtx[5]; i++) binEdgesSecvtx[5][i]=(Double_t)kNtrksmin + (kNtrksmax - kNtrksmin)/nBinSecvtx[5]*(Double_t)i;
1099
1100   fSecvtxQA = new THnSparseF("secvtxQA", "QA for Secvtx; invmass[GeV/c^2]; KF chi2; signed Lxy; KF ip; pair code; n tracks ", nDimSecvtx, nBinSecvtx);
1101   for(Int_t idim = 0; idim < nDimSecvtx; idim++){
1102     fSecvtxQA->SetBinEdges(idim, binEdgesSecvtx[idim]);
1103   }
1104
1105   fSecVtxList->AddAt(fSecvtxQA,1);
1106 }