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  *                                                                        *
  17  * Secondary vertexing construction Class                                 *
  18  *                                                                        *
  19  * Authors:                                                               *
  20  *   MinJung Kweon <minjung@physi.uni-heidelberg.de>                      *
  21  *                                                                        *
  22  **************************************************************************/
  23
  24 #include <iostream>
  25 #include <TH2F.h>
  26 #include <TH3F.h>
  27 #include <TLorentzVector.h>
  28 #include <TNtuple.h>
  29 #include <TParticle.h>
  30
  31 #include <AliESDEvent.h>
  32 #include <AliESDtrack.h>
  33 #include "AliHFEsecVtx.h"
  34 #include <AliKFParticle.h>
  35 #include <AliKFVertex.h>
  36 #include <AliLog.h>
  37 #include <AliPID.h>
  38 #include <AliStack.h>
  39
  40 ClassImp(AliHFEsecVtx)
  41
  42 //_______________________________________________________________________________________________
  43 AliHFEsecVtx::AliHFEsecVtx():
  44          fESD1(0x0)
  45         ,fStack(0x0)
  46         ,fNparents(0)
  47         ,fHistTagged()
  48         ,fPairTagged(0)
  49         ,fdistTwoSecVtx(-1)
  50         ,fcosPhi(-1)
  51         ,fsignedLxy(-1)
  52         ,finvmass(-1)
  53         ,finvmassSigma(-1)
  54         ,fBTagged(0)
  55 {
  56         // Default constructor
  57
  58         Init();
  59
  60 }
  61
  62 //_______________________________________________________________________________________________
  63 AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
  64          TObject(p)
  65         ,fESD1(0x0)
  66         ,fStack(0x0)
  67         ,fNparents(p.fNparents)
  68         ,fHistTagged()
  69         ,fPairTagged(p.fPairTagged)
  70         ,fdistTwoSecVtx(p.fdistTwoSecVtx)
  71         ,fcosPhi(p.fcosPhi)
  72         ,fsignedLxy(p.fsignedLxy)
  73         ,finvmass(p.finvmass)
  74         ,finvmassSigma(p.finvmassSigma)
  75         ,fBTagged(p.fBTagged)
  76 {
  77         // Copy constructor
  78 }
  79
  80 //_______________________________________________________________________________________________
  81 AliHFEsecVtx&
  82 AliHFEsecVtx::operator=(const AliHFEsecVtx &)
  83 {
  84   // Assignment operator
  85
  86   AliInfo("Not yet implemented.");
  87   return *this;
  88 }
  89
  90 //_______________________________________________________________________________________________
  91 AliHFEsecVtx::~AliHFEsecVtx()
  92 {
  93         // Destructor
  94
  95         //cout << "Analysis Done." << endl;
  96 }
  97
  98 //__________________________________________
  99 void AliHFEsecVtx::Init()
 100 {
 101
 102         fNparents = 7;
 103
 104         fParentSelect[0][0] =  411;
 105         fParentSelect[0][1] =  421;
 106         fParentSelect[0][2] =  431;
 107         fParentSelect[0][3] = 4122;
 108         fParentSelect[0][4] = 4132;
 109         fParentSelect[0][5] = 4232;
 110         fParentSelect[0][6] = 4332;
 111
 112         fParentSelect[1][0] =  511;
 113         fParentSelect[1][1] =  521;
 114         fParentSelect[1][2] =  531;
 115         fParentSelect[1][3] = 5122;
 116         fParentSelect[1][4] = 5132;
 117         fParentSelect[1][5] = 5232;
 118         fParentSelect[1][6] = 5332;
 119
 120         fid[0][0] = 0;
 121         fid[0][1] = 1;
 122         fid[0][2] = 2;
 123
 124         fid[1][0] = 0;
 125         fid[1][1] = 1;
 126         fid[1][2] = 3;
 127
 128         fid[2][0] = 0;
 129         fid[2][1] = 2;
 130         fid[2][2] = 3;
 131
 132         fid[3][0] = 1;
 133         fid[3][1] = 2;
 134         fid[3][2] = 3;
 135
 136         fia[0][0][0] = 0;
 137         fia[0][0][1] = 1;
 138         fia[0][1][0] = 0;
 139         fia[0][1][1] = 2;
 140         fia[0][2][0] = 1;
 141         fia[0][2][1] = 2;
 142
 143         fia[1][0][0] = 0;
 144         fia[1][0][1] = 1;
 145         fia[1][1][0] = 0;
 146         fia[1][1][1] = 3;
 147         fia[1][2][0] = 1;
 148         fia[1][2][1] = 3;
 149
 150         fia[2][0][0] = 0;
 151         fia[2][0][1] = 2;
 152         fia[2][1][0] = 0;
 153         fia[2][1][1] = 3;
 154         fia[2][2][0] = 2;
 155         fia[2][2][1] = 3;
 156
 157         fia[3][0][0] = 1;
 158         fia[3][0][1] = 2;
 159         fia[3][1][0] = 1;
 160         fia[3][1][1] = 3;
 161         fia[3][2][0] = 2;
 162         fia[3][2][1] = 3;
 163
 164
 165   //fid[4][3] = {0,1,2, 0,1,3, 0,2,3, 1,2,3};
 166   //fia[4][3][2] = {0,1, 0,2, 1,2,  0,1, 0,3, 1,3,  0,2, 0,3, 2,3,  1,2, 1,3, 2,3};
 167
 168 }
 169
 170 //__________________________________________
 171 void AliHFEsecVtx::ResetTagVar()
 172 {
 173         fdistTwoSecVtx = -1;
 174         fcosPhi = -1;
 175         fsignedLxy = -1;
 176         finvmass = -1;
 177         finvmassSigma = -1;
 178         fBTagged = kFALSE;
 179 }
 180
 181 //__________________________________________
 182 void AliHFEsecVtx::InitAnaPair()
 183 {
 184         fPairTagged = 0;
 185         for (Int_t i=0; i<20; i++){
 186                 fpairedTrackID[i] = -1;
 187                 fpairedChi2[i] = -1;
 188                 fpairedInvMass[i] = -1;
 189                 fpairedSignedLxy[i] = -1;
 190         }
 191
 192 }
 193
 194 //_______________________________________________________________________________________________
 195 void AliHFEsecVtx::CreateHistograms(TString hnopt)
 196 {
 197         // create histograms
 198
 199         fkSourceLabel[fkAll]="all";
 200         fkSourceLabel[fkDirectCharm]="directCharm";
 201         fkSourceLabel[fkDirectBeauty]="directBeauty";
 202         fkSourceLabel[fkBeautyCharm]="beauty2charm";
 203         fkSourceLabel[fkGamma]="gamma";
 204         fkSourceLabel[fkPi0]="pi0";
 205         fkSourceLabel[fkElse]="others";
 206         fkSourceLabel[fkBeautyGamma]="beauty22gamma";
 207         fkSourceLabel[fkBeautyPi0]="beauty22pi0";
 208         fkSourceLabel[fkBeautyElse]="beauty22others";
 209
 210
 211         TString hname;
 212         for (Int_t isource = 0; isource < 10; isource++ ){
 213
 214            hname=hnopt+"InvMass_"+fkSourceLabel[isource];
 215            fHistPair[isource].fInvMass = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
 216            hname=hnopt+"InvMassCut1_"+fkSourceLabel[isource];
 217            fHistPair[isource].fInvMassCut1 = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
 218            hname=hnopt+"InvMassCut2_"+fkSourceLabel[isource];
 219            fHistPair[isource].fInvMassCut2 = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
 220            hname=hnopt+"KFChi2_"+fkSourceLabel[isource];
 221            fHistPair[isource].fKFChi2 = new TH1F(hname,hname,200,0,20);
 222            hname=hnopt+"CosOpenAngle_"+fkSourceLabel[isource];
 223            fHistPair[isource].fCosOpenAngle = new TH1F(hname,hname,100,-1.1,1.1);
 224            hname=hnopt+"SignedLxy_"+fkSourceLabel[isource];
 225            fHistPair[isource].fSignedLxy = new TH2F(hname,hname,1000,-5,5,120,-2,10);
 226            hname=hnopt+"KFIP_"+fkSourceLabel[isource];
 227            fHistPair[isource].fKFIP = new TH1F(hname,hname,1000,-5,5);
 228
 229         }
 230
 231         hname=hnopt+"pt_beautyelec";
 232         fHistTagged.fPtBeautyElec= new TH1F(hname,hname,150,0,30);
 233         hname=hnopt+"pt_taggedelec";
 234         fHistTagged.fPtTaggedElec= new TH1F(hname,hname,150,0,30);
 235         hname=hnopt+"pt_righttaggedelec";
 236         fHistTagged.fPtRightTaggedElec = new TH1F(hname,hname,150,0,30);
 237         hname=hnopt+"pt_wrongtaggedelec";
 238         fHistTagged.fPtWrongTaggedElec = new TH1F(hname,hname,150,0,30);
 239
 240 }
 241
 242 //_______________________________________________________________________________________________
 243 void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t sourcePart, Int_t index2)
 244 {
 245
 246         // get KF particle input pid
 247         Int_t pdg1 = GetMCPID(track1);
 248         Int_t pdg2 = GetMCPID(track2);
 249
 250
 251         // create KF particle of pair
 252         AliKFParticle::SetField(fESD1->GetMagneticField());
 253         AliKFParticle kfTrack1(*track1, pdg1);
 254         AliKFParticle kfTrack2(*track2, pdg2);
 255
 256         AliKFParticle kfSecondary(kfTrack1,kfTrack2);
 257
 258         // copy primary vertex from ESD
 259         AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
 260         if( primVtxCopy.GetNDF() <1 ) return;
 261
 262         //primary vertex point
 263         Double_t pvx = primVtxCopy.GetX();
 264         Double_t pvy = primVtxCopy.GetY();
 265         //Double_t pvz = primVtxCopy.GetZ();
 266
 267         //secondary vertex point from kf particle
 268         Double_t kfx = kfSecondary.GetX();
 269         Double_t kfy = kfSecondary.GetY();
 270         //Double_t kfz = kfSecondary.GetZ();
 271
 272         //momentum at the decay point from kf particle
 273         Double_t kfpx = kfSecondary.GetPx();
 274         Double_t kfpy = kfSecondary.GetPy();
 275         //Double_t kfpz = kfSecondary.GetPz();
 276
 277
 278         Double_t dx = kfx-pvx;
 279         Double_t dy = kfy-pvy;
 280
 281
 282
 283         // discriminating variables ----------------------------------------------------------
 284
 285         // invariant mass of the KF particle
 286         Double_t invmass = -1;
 287         Double_t invmassSigma = -1;
 288         kfSecondary.GetMass(invmass,invmassSigma);
 289
 290         // chi2 of the KF particle
 291         Double_t kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 292
 293         // opening angle between two particles in XY plane
 294         Double_t cosphi = TMath::Cos(kfTrack1.GetAngleXY(kfTrack2));
 295
 296         // projection of kf vertex vector to the kf momentum direction
 297         Double_t costheta = ( dx*kfpx + dy*kfpy)/TMath::Sqrt(dx*dx+dy*dy)*TMath::Sqrt(kfpx*kfpx + kfpy*kfpy);
 298         Double_t signedLxy = TMath::Sqrt(dx*dx+dy*dy)*costheta;
 299
 300         // DCA from primary to e-h KF particle (impact parameter of KF particle)
 301         Double_t vtx[2]={pvx, pvy};
 302         Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
 303
 304
 305
 306
 307         // pair cuts
 308         if( kfchi2 >2. ) return;
 309
 310         // fill histograms
 311         fHistPair[sourcePart].fInvMass->Fill(invmass,invmassSigma);
 312         fHistPair[sourcePart].fKFChi2->Fill(kfchi2);
 313         fHistPair[sourcePart].fCosOpenAngle->Fill(cosphi);
 314         fHistPair[sourcePart].fSignedLxy->Fill(signedLxy,invmass);
 315         fHistPair[sourcePart].fKFIP->Fill(kfip);
 316
 317         if ( signedLxy > 0.05 && cosphi > 0.5) {
 318           fHistPair[sourcePart].fInvMassCut1->Fill(invmass,invmassSigma);
 319         }
 320
 321         // pair cuts
 322         if (signedLxy < 0.05) return;
 323         if (cosphi < 0.0) return;
 324
 325         // pair tagging if it passed the above cuts
 326         fHistPair[sourcePart].fInvMassCut2->Fill(invmass,invmassSigma);
 327
 328
 329         // pair tagging condition
 330         if ( signedLxy > 0.0 && cosphi > 0) {
 331         //if ( signedLxy > 0.06 && cosphi > 0) {
 332           fpairedTrackID[fPairTagged] = index2;
 333           fpairedChi2[fPairTagged] = kfchi2;
 334           fpairedInvMass[fPairTagged] = invmass;
 335           fpairedSignedLxy[fPairTagged] = signedLxy;
 336           fPairTagged++;
 337         }
 338
 339 }
 340
 341 //_______________________________________________________________________________________________
 342 void AliHFEsecVtx::RunSECVTX(AliESDtrack *track)
 343 {
 344
 345         ResetTagVar();
 346
 347         if (fPairTagged >= 4) {
 348           FindSECVTXCandid4Tracks(track);
 349         }
 350         else if (fPairTagged == 3) {
 351           FindSECVTXCandid3Tracks(track);
 352         }
 353         else if (fPairTagged == 2) {
 354           FindSECVTXCandid2Tracks(track);
 355         }
 356         else if (fPairTagged == 1) {
 357           ApplyPairTagCut();
 358         }
 359
 360
 361         Int_t imclabel = TMath::Abs(track->GetLabel());
 362         if(imclabel<0) return;
 363         TParticle* mcpart = fStack->Particle(imclabel);
 364         Int_t esource = GetElectronSource(imclabel);
 365         if (esource == fkDirectBeauty || esource == fkBeautyCharm || esource == fkBeautyGamma || esource == fkBeautyPi0 || esource == fkBeautyElse){
 366                 fHistTagged.fPtBeautyElec->Fill(mcpart->Pt());
 367         }
 368
 369
 370         if (fBTagged) {
 371                 fHistTagged.fPtTaggedElec->Fill(mcpart->Pt());
 372                 if (esource == fkDirectBeauty || esource == fkBeautyCharm || esource == fkBeautyGamma || esource == fkBeautyPi0 || esource == fkBeautyElse){
 373                         fHistTagged.fPtRightTaggedElec->Fill(mcpart->Pt());
 374                 }
 375                 else fHistTagged.fPtWrongTaggedElec->Fill(mcpart->Pt());
 376         }
 377
 378 }
 379
 380 //_______________________________________________________________________________________________
 381 void AliHFEsecVtx::ApplyPairTagCut()
 382 {
 383
 384         if (fpairedChi2[0] > 2.0) return;
 385         if (fpairedInvMass[0] < 1.5) return;
 386         if (fpairedSignedLxy[0] < 0.5) return;
 387
 388         fBTagged = kTRUE;
 389
 390 }
 391
 392 //_______________________________________________________________________________________________
 393 Bool_t AliHFEsecVtx::ApplyPairTagCut(Int_t id)
 394 {
 395
 396         if (fpairedChi2[id] > 2.0) return kFALSE;
 397         if (fpairedInvMass[id] < 1.5) return kFALSE;
 398         if (fpairedSignedLxy[id] < 0.5) return kFALSE;
 399
 400         fBTagged = kTRUE;
 401         return kTRUE;
 402
 403 }
 404
 405 //_______________________________________________________________________________________________
 406 Bool_t AliHFEsecVtx::ApplyTagCut(Double_t chi2)
 407 {
 408
 409         if (chi2 > 2.0) return kFALSE;
 410         if (finvmass < 1.5) return kFALSE;
 411         if (fsignedLxy < 0.5) return kFALSE;
 412         if (fcosPhi < 0.90) return kFALSE;
 413         if (fdistTwoSecVtx > 0.1) return kFALSE;
 414
 415         fBTagged = kTRUE;
 416         return kTRUE;
 417 }
 418
 419 //_______________________________________________________________________________________________
 420 //void AliHFEsecVtx::ApplyTagCut(Double_t chi2, AliESDtrack *track, AliESDtrack *htrack1, AliESDtrack *htrack2)
 421 void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2)
 422 {
 423
 424         if(!ApplyTagCut(chi2)){
 425                 if(!ApplyPairTagCut(0)) ApplyPairTagCut(1);
 426         }
 427
 428 }
 429
 430 //_______________________________________________________________________________________________
 431 void AliHFEsecVtx::FindSECVTXCandid4Tracks(AliESDtrack *track)
 432 {
 433
 434         // sort pair in increasing order (kFALSE-increasing order)
 435         Int_t index[20];
 436         Int_t indexA[4];
 437         Int_t indexB[3];
 438         Double_t sevchi2[4];
 439         AliESDtrack *htrack[4];
 440
 441         TMath::Sort(fPairTagged,fpairedChi2,index,kFALSE);
 442
 443         // select 4 partner tracks retruning smallest pair chi2
 444
 445         for (Int_t i=0; i<4; i++){
 446                 htrack[i] = fESD1->GetTrack(fpairedTrackID[index[i]]);
 447         }
 448
 449         // calculated chi2 with 1 electron and 3 partner tracks
 450         for (Int_t i=0; i<4; i++){
 451                 sevchi2[i] = GetSecVtxChi2(track, htrack[fid[i][0]], htrack[fid[i][1]], htrack[fid[i][2]]);
 452         }
 453
 454         // select 3 partner tracks retruning smallest pair chi2
 455         // [think] if two smallest chi2 are similar, have to think about better handling of selection
 456         TMath::Sort(4,sevchi2,indexA,kFALSE);
 457
 458         // calculated chi2 with 1 electron and 2 partner tracks
 459         for (Int_t i=0; i<3; i++){
 460                 sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
 461         }
 462
 463         // select 2 partner tracks retruning smallest pair chi2
 464         TMath::Sort(3,sevchi2,indexB,kFALSE);
 465
 466         // calculate secondary vertex quality variables with 1 electron and 2 hadrons
 467         CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
 468
 469         ApplyVtxTagCut(sevchi2[indexB[0]]);
 470         //ApplyTagCut(sevchi2[indexB[0]],track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
 471
 472 }
 473
 474 //_______________________________________________________________________________________________
 475 void AliHFEsecVtx::FindSECVTXCandid3Tracks(AliESDtrack *track)
 476 {
 477
 478         // sort pair in increasing order (kFALSE-increasing order)
 479         Int_t indexA[1] = { 0 };
 480         Int_t indexB[3];
 481         Double_t sevchi2[3];
 482         AliESDtrack *htrack[3];
 483
 484         // select 4 partner tracks retruning smallest pair chi2
 485
 486         for (Int_t i=0; i<3; i++){
 487                 htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
 488         }
 489
 490         // calculated chi2 with 1 electron and 2 partner tracks
 491         for (Int_t i=0; i<3; i++){
 492                 sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
 493         }
 494
 495         // select 2 partner tracks retruning smallest pair chi2
 496         TMath::Sort(3,sevchi2,indexB,kFALSE);
 497
 498         // calculate secondary vertex quality variables with 1 electron and 2 hadrons
 499         CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
 500
 501         ApplyVtxTagCut(sevchi2[indexB[0]]);
 502         //ApplyTagCut(sevchi2[indexB[0]],track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
 503 }
 504
 505 //_______________________________________________________________________________________________
 506 void AliHFEsecVtx::FindSECVTXCandid2Tracks(AliESDtrack *track)
 507 {
 508
 509         Double_t sevchi2[1];
 510         AliESDtrack *htrack[2];
 511
 512         for (Int_t i=0; i<2; i++){
 513                 htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
 514         }
 515
 516         sevchi2[0] = GetSecVtxChi2(track, htrack[0], htrack[1]);
 517
 518         // calculate secondary vertex quality variables with 1 electron and 2 hadrons
 519         CalcSECVTXProperty(track,htrack[0],htrack[1]);
 520
 521         ApplyVtxTagCut(sevchi2[0]);
 522         //ApplyTagCut(sevchi2[0],track,htrack[0],htrack[1]);
 523 }
 524
 525 //_______________________________________________________________________________________________
 526 void AliHFEsecVtx::CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
 527 {
 528
 529         Int_t pdg1 = GetMCPID(track1);
 530         Int_t pdg2 = GetMCPID(track2);
 531         Int_t pdg3 = GetMCPID(track3);
 532
 533         AliKFParticle::SetField(fESD1->GetMagneticField());
 534         AliKFParticle kfTrack1(*track1, pdg1);
 535         AliKFParticle kfTrack2(*track2, pdg2);
 536         AliKFParticle kfTrack3(*track3, pdg3);
 537
 538         AliKFParticle kfSecondary12(kfTrack1,kfTrack2);
 539         AliKFParticle kfSecondary13(kfTrack1,kfTrack3);
 540         AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
 541
 542         // copy primary vertex from ESD
 543         AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
 544         //printf("primary ndf= %d\n",primVtxCopy.GetNDF());
 545         if( primVtxCopy.GetNDF() <1 ) return;
 546
 547         Double_t kdx12 = kfSecondary12.GetX()-primVtxCopy.GetX();
 548         Double_t kdy12 = kfSecondary12.GetY()-primVtxCopy.GetY();
 549         //Double_t kdz12 = kfSecondary12.GetZ()-primVtxCopy.GetZ();
 550
 551         Double_t kdx13 = kfSecondary13.GetX()-primVtxCopy.GetX();
 552         Double_t kdy13 = kfSecondary13.GetY()-primVtxCopy.GetY();
 553         //Double_t kdz13 = kfSecondary13.GetZ()-primVtxCopy.GetZ();
 554
 555         Double_t kdx = kfSecondary.GetX()-primVtxCopy.GetX();
 556         Double_t kdy = kfSecondary.GetY()-primVtxCopy.GetY();
 557         //Double_t kdz = kfSecondary.GetZ()-primVtxCopy.GetZ();
 558
 559         // calculate distance and angle between two secvtxes
 560         fdistTwoSecVtx = TMath::Sqrt((kdx12-kdx13)*(kdx12-kdx13) + (kdy12-kdy13)*(kdy12-kdy13));
 561         fcosPhi = ( kdx12*kdx13 + kdy12*kdy13 ) / ( TMath::Sqrt(kdx12*kdx12+kdy12*kdy12)*TMath::Sqrt(kdx13*kdx13+kdy13*kdy13) );
 562         //Double_t lengthdiff = TMath::Abs(TMath::Sqrt(kdx12*kdx12+kdy12*kdy12) - TMath::Sqrt(kdx13*kdx13+kdy13*kdy13));
 563
 564         // calculate angle between secondary vertex vector and secondary particle momentum vector in transverse plane
 565         Double_t cosTheta = ( kdx*kfSecondary.GetPx() + kdy*kfSecondary.GetPy()) / TMath::Sqrt(kdx*kdx+kdy*kdy)*TMath::Sqrt(kfSecondary.GetPx()*kfSecondary.GetPx()+kfSecondary.GetPy()*kfSecondary.GetPy());
 566         // calculate signed Lxy
 567         fsignedLxy = TMath::Sqrt(kdx*kdx+kdy*kdy)*cosTheta;
 568
 569         // calculate invariant mass of the kf particle
 570         kfSecondary.GetMass(finvmass,finvmassSigma);
 571 }
 572
 573 //_______________________________________________________________________________________________
 574 Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track)
 575 {
 576
 577         Int_t label = TMath::Abs(track->GetLabel());
 578         TParticle* mcpart = fStack->Particle(label);
 579         if ( !mcpart ) return 0;
 580         Int_t pdgCode = mcpart->GetPdgCode();
 581
 582         return pdgCode;
 583 }
 584
 585 //_______________________________________________________________________________________________
 586 Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3, AliESDtrack* track4)
 587 {
 588
 589         Int_t pdg1 = GetMCPID(track1);
 590         Int_t pdg2 = GetMCPID(track2);
 591         Int_t pdg3 = GetMCPID(track3);
 592         Int_t pdg4 = GetMCPID(track4);
 593
 594         AliKFParticle::SetField(fESD1->GetMagneticField());
 595         AliKFParticle kfTrack1(*track1, pdg1);
 596         AliKFParticle kfTrack2(*track2, pdg2);
 597         AliKFParticle kfTrack3(*track3, pdg3);
 598         AliKFParticle kfTrack4(*track4, pdg4);
 599         AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3,kfTrack4);
 600
 601         return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 602
 603 }
 604
 605 //_______________________________________________________________________________________________
 606 Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
 607 {
 608
 609         Int_t pdg1 = GetMCPID(track1);
 610         Int_t pdg2 = GetMCPID(track2);
 611         Int_t pdg3 = GetMCPID(track3);
 612
 613         AliKFParticle::SetField(fESD1->GetMagneticField());
 614         AliKFParticle kfTrack1(*track1, pdg1);
 615         AliKFParticle kfTrack2(*track2, pdg2);
 616         AliKFParticle kfTrack3(*track3, pdg3);
 617         AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
 618
 619         return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 620
 621 }
 622
 623 //_______________________________________________________________________________________________
 624 Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2)
 625 {
 626
 627         Int_t pdg1 = GetMCPID(track1);
 628         Int_t pdg2 = GetMCPID(track2);
 629
 630         AliKFParticle::SetField(fESD1->GetMagneticField());
 631         AliKFParticle kfTrack1(*track1, pdg1);
 632         AliKFParticle kfTrack2(*track2, pdg2);
 633         AliKFParticle kfSecondary(kfTrack1,kfTrack2);
 634
 635         return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 636
 637 }
 638
 639 //_______________________________________________________________________________________________
 640 Int_t AliHFEsecVtx::PairOrigin(AliESDtrack* track1, AliESDtrack* track2)
 641 {
 642
 643         //
 644         // return pdg code of the origin(source) of the pair
 645         //
 646         //
 647         // ---*---*---*-----ancester A----- track1
 648         //                        |____*______
 649         //                             |______ track2
 650         // => if they originated from same ancester,
 651         //    then return "the absolute value of pdg code of ancester A"
 652         //
 653         // ---*---*---B hadron-----ancester A----- track1
 654         //                               |____*______
 655         //                                    |______ track2
 656         // => if they originated from same ancester, and this ancester originally comes from B hadrons
 657         //    then return -1*"the absolute value of pdg code of ancester A"
 658         //
 659         // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
 660         //
 661
 662         TParticle* part1 = fStack->Particle(TMath::Abs(track1->GetLabel()));
 663         TParticle* part2 = fStack->Particle(TMath::Abs(track2->GetLabel()));
 664         if (!(part1)) return 0;
 665         if (!(part2)) return 0;
 666
 667         TParticle* part1Crtgen = part1; // copy track into current generation particle
 668         TParticle* part2Crtgen = part2; // copy track into current generation particle
 669
 670
 671         Int_t sourcePDG = 0;
 672
 673         // if the two tracks' mother's label is same, get the mother info
 674         // in case of charm, check if it originated from beauty
 675         for (Int_t i=0; i<100; i++){ // iterate 100
 676                 Int_t iLabel = part1Crtgen->GetFirstMother(); //label of mother of current generation for 1st partilce
 677                 if (iLabel < 0) return 0;
 678
 679                 for (Int_t j=0; j<100; j++){ // iterate 100
 680                         Int_t jLabel = part2Crtgen->GetFirstMother(); //label of mother of current generation for 2nd partilce
 681                         if (jLabel < 0) return 0; // if jLabel == -1
 682
 683                         if (iLabel == jLabel){ // check if two tracks are originated from same mother
 684                                 TParticle* thismother = fStack->Particle(jLabel); // if yes, get "thismother" info
 685                                 sourcePDG = abs(thismother->GetPdgCode()); // get the pdg code of "this mother"
 686
 687                                 // check ancester to see if it is originally from beauty
 688                                 for (Int_t k=0; k<10; k++){ // check up to 10 ancesters
 689                                         Int_t ancesterLabel = thismother->GetFirstMother();
 690                                         if (ancesterLabel < 0) return sourcePDG; // if "thismoter" doesn't have mother anymore, return thismother's pdg
 691
 692                                         TParticle* thisancester = fStack->Particle(ancesterLabel);
 693                                         Int_t ancesterPDG = abs(thisancester->GetPdgCode());
 694
 695                                         for (Int_t l=0; l<fNparents; l++){
 696                                                 if (abs(ancesterPDG)==fParentSelect[1][l]){
 697                                                         sourcePDG = -1*sourcePDG; // multiply -1 for charm from bottom
 698                                                         return sourcePDG;
 699                                                 }
 700                                         }
 701                                         thismother = thisancester;
 702                                 }
 703
 704                         }
 705                         part2Crtgen = fStack->Particle(jLabel); // if their mother is different, go up to previous generation of 2nd particle
 706                 }
 707                 part1Crtgen = fStack->Particle(iLabel); // if their mother is different, go up to previous generation of 2nd particle
 708
 709                 // if you don't find additionional(2nd particle) track originated from a given beauty hadron, break to save time
 710                 Int_t motherPDGtmp = abs(part1Crtgen->GetPdgCode());
 711                 for (Int_t l=0; l<fNparents; l++){
 712                         if (abs(motherPDGtmp)==fParentSelect[1][l]) return sourcePDG;
 713                 }
 714
 715         }
 716
 717         return sourcePDG;
 718
 719 }
 720
 721 //_______________________________________________________________________________________________
 722 Int_t AliHFEsecVtx::PairCode(AliESDtrack* track1,AliESDtrack* track2)
 723 {
 724
 725         //
 726         // return pair code which is predefind as:
 727         //  fkDirectCharm, fkDirectBeauty, fkBeautyCharm, fkGamma, fkPi0, fkElse, fkBeautyGamma, fkBeautyPi0, fkBeautyElse
 728         //
 729
 730         Int_t pairOriginsPDG = PairOrigin(track1,track2);
 731
 732         Int_t sourcePart = fkElse;
 733
 734         if (pairOriginsPDG < 0) {
 735                 sourcePart = fkBeautyElse;
 736         }
 737         for (Int_t i=0; i<fNparents; i++){
 738                 if (abs(pairOriginsPDG)==fParentSelect[0][i]) {
 739                         if (pairOriginsPDG>0) sourcePart = fkDirectCharm;
 740                         if (pairOriginsPDG<0) {
 741                                 sourcePart = fkBeautyCharm;
 742                         }
 743                 }
 744                 if (abs(pairOriginsPDG)==fParentSelect[1][i]) {
 745                         if (pairOriginsPDG>0) {
 746                                 sourcePart = fkDirectBeauty;
 747                         }
 748                         if (pairOriginsPDG<0)  return fkElse;
 749                 }
 750         }
 751         if (pairOriginsPDG == 22) sourcePart = fkGamma;
 752         if (pairOriginsPDG == -22) {
 753                 sourcePart = fkBeautyGamma;
 754         }
 755         if (pairOriginsPDG == 111) sourcePart = fkPi0;
 756         if (pairOriginsPDG == -111) {
 757                 sourcePart = fkBeautyPi0;
 758         }
 759
 760         return sourcePart;
 761
 762 }
 763
 764 //_______________________________________________________________________________________________
 765 Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
 766 {
 767
 768         //
 769         // decay electron's origin
 770         //
 771
 772         if (iTrackLabel < 0) {
 773                 AliDebug(1, "Stack label is negative, return\n");
 774                 return -1;
 775         }
 776
 777         TParticle* mcpart = fStack->Particle(iTrackLabel);
 778         Int_t iLabel = mcpart->GetFirstMother();
 779         if (iLabel<0){
 780                 AliDebug(1, "Stack label is negative, return\n");
 781                 return -1;
 782         }
 783
 784         Int_t origin = -1;
 785         Bool_t IsFinalOpenCharm = kFALSE;
 786
 787         TParticle *partMother = fStack->Particle(iLabel);
 788         Int_t maPdgcode = partMother->GetPdgCode(); // get mother's pdg code
 789
 790         //beauty --------------------------
 791         if ( int(abs(maPdgcode)/100.) == fkBeauty || int(abs(maPdgcode)/1000.) == fkBeauty ) {
 792                 for (Int_t i=0; i<fNparents; i++){
 793                         if (abs(maPdgcode)==fParentSelect[1][i]){
 794                                 origin = fkDirectBeauty;
 795                                 return origin;
 796                         }
 797                         else return -1; // this track is originated beauties not in the final B hadron list => excited beauties
 798                 }
 799         } // end of if
 800
 801         //charm --------------------------
 802         else if ( int(abs(maPdgcode)/100.) == fkCharm || int(abs(maPdgcode)/1000.) == fkCharm ) {
 803
 804                 for (Int_t i=0; i<fNparents; i++){
 805                         if (abs(maPdgcode)==fParentSelect[0][i])
 806                                 IsFinalOpenCharm = kTRUE;
 807                 }
 808                 if (!IsFinalOpenCharm) return -1; // this track is originated charms not in the final D hadron list => excited charms
 809                                                   // to prevent any possible double counting
 810
 811                 for (Int_t i=0; i<100; i++){ // iterate 100 until you find B hadron as a mother or become top ancester
 812
 813                         Int_t jLabel = partMother->GetFirstMother();
 814                         if (jLabel == -1){
 815                                 origin = fkDirectCharm;
 816                                 return origin;
 817                         }
 818                         if (jLabel < 0){ // safety protection even though not really necessary here
 819                                 AliDebug(1, "Stack label is negative, return\n");
 820                                 return -1;
 821                         }
 822
 823                         // if there is an ancester, check if it in the final B hadron list
 824                         TParticle* grandMa = fStack->Particle(jLabel);
 825                         Int_t grandMaPDG = grandMa->GetPdgCode();
 826
 827                         for (Int_t j=0; j<fNparents; j++){
 828                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
 829                                 origin = fkBeautyCharm;
 830                                 return origin;
 831                                 }
 832                         }
 833
 834                         partMother = grandMa;
 835                 } // end of iteration
 836         } // end of if
 837
 838         //gamma --------------------------
 839         else if ( abs(maPdgcode) == 22 ) {
 840                 origin = fkGamma;
 841
 842                 // iterate until you find B hadron as a mother or become top ancester
 843                 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
 844
 845                         Int_t jLabel = partMother->GetFirstMother();
 846                         if (jLabel == -1){
 847                                 origin = fkGamma;
 848                                 return origin;
 849                         }
 850                         if (jLabel < 0){ // safety protection
 851                                 AliDebug(1, "Stack label is negative, return\n");
 852                                 return -1;
 853                         }
 854
 855                         // if there is an ancester
 856                         TParticle* grandMa = fStack->Particle(jLabel);
 857                         Int_t grandMaPDG = grandMa->GetPdgCode();
 858
 859                         for (Int_t j=0; j<fNparents; j++){
 860                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
 861                                         origin = fkBeautyGamma;
 862                                         return origin;
 863                                 }
 864                         }
 865
 866                         partMother = grandMa;
 867                 } // end of iteration
 868
 869                 return origin;
 870         } // end of if
 871
 872         //pi0 --------------------------
 873         else if ( abs(maPdgcode) == 111 ) {
 874                 origin = fkPi0;
 875
 876                 // iterate until you find B hadron as a mother or become top ancester
 877                 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
 878
 879                         Int_t jLabel = partMother->GetFirstMother();
 880                         if (jLabel == -1){
 881                                 origin = fkPi0;
 882                                 return origin;
 883                         }
 884                         if (jLabel < 0){ // safety protection
 885                                 AliDebug(1, "Stack label is negative, return\n");
 886                                 return -1;
 887                         }
 888
 889                         // if there is an ancester
 890                         TParticle* grandMa = fStack->Particle(jLabel);
 891                         Int_t grandMaPDG = grandMa->GetPdgCode();
 892
 893                         for (Int_t j=0; j<fNparents; j++){
 894                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
 895                                         origin = fkBeautyPi0;
 896                                         return origin;
 897                                 }
 898                         }
 899
 900                         partMother = grandMa;
 901                 } // end of iteration
 902
 903                 return origin;
 904         } // end of if
 905
 906
 907         //else --------------------------
 908         else {
 909                 origin = fkElse;
 910
 911                 // iterate until you find B hadron as a mother or become top ancester
 912                 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
 913
 914                         Int_t jLabel = partMother->GetFirstMother();
 915                         if (jLabel == -1){
 916                                 origin = fkElse;
 917                                 return origin;
 918                         }
 919                         if (jLabel < 0){ // safety protection
 920                                 AliDebug(1, "Stack label is negative, return\n");
 921                                 return -1;
 922                         }
 923
 924                         // if there is an ancester
 925                         TParticle* grandMa = fStack->Particle(jLabel);
 926                         Int_t grandMaPDG = grandMa->GetPdgCode();
 927
 928                         for (Int_t j=0; j<fNparents; j++){
 929                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
 930                                         origin = fkBeautyElse;
 931                                         return origin;
 932                                 }
 933                         }
 934
 935                         partMother = grandMa;
 936                 } // end of iteration
 937
 938         }
 939
 940         return origin;
 941
 942 }
 943
 944 //_______________________________________________________________________________________________
 945 Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track)
 946 {
 947
 948         //if (track->Pt() < 1.0) return kFALSE;
 949         //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
 950         //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
 951         //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
 952         if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
 953         //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
 954
 955
 956 /*
 957         Float_t dcaR=-1;
 958         Float_t dcaZ=-1;
 959         track->GetImpactParameters(dcaR,dcaZ);
 960         if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
 961         if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
 962 */
 963         return kTRUE;
 964 }