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 <AliESDtrack.h>
  29 #include "AliHFEsecVtx.h"
  30 #include <AliKFParticle.h>
  31 #include <AliKFVertex.h>
  32 #include <AliLog.h>
  33 #include <AliStack.h>
  34
  35 ClassImp(AliHFEsecVtx)
  36
  37 //_______________________________________________________________________________________________
  38 AliHFEsecVtx::AliHFEsecVtx():
  39          fESD1(0x0)
  40         ,fStack(0x0)
  41         ,fNparents(0)
  42         ,fHistTagged()
  43         ,fPairTagged(0)
  44         ,fdistTwoSecVtx(-1)
  45         ,fcosPhi(-1)
  46         ,fsignedLxy(-1)
  47         ,finvmass(-1)
  48         ,finvmassSigma(-1)
  49         ,fBTagged(0)
  50         ,fBElec(0)
  51 {
  52         // Default constructor
  53
  54         Init();
  55
  56 }
  57
  58 //_______________________________________________________________________________________________
  59 AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
  60          TObject(p)
  61         ,fESD1(0x0)
  62         ,fStack(0x0)
  63         ,fNparents(p.fNparents)
  64         ,fHistTagged()
  65         ,fPairTagged(p.fPairTagged)
  66         ,fdistTwoSecVtx(p.fdistTwoSecVtx)
  67         ,fcosPhi(p.fcosPhi)
  68         ,fsignedLxy(p.fsignedLxy)
  69         ,finvmass(p.finvmass)
  70         ,finvmassSigma(p.finvmassSigma)
  71         ,fBTagged(p.fBTagged)
  72         ,fBElec(p.fBElec)
  73 {
  74         // Copy constructor
  75 }
  76
  77 //_______________________________________________________________________________________________
  78 AliHFEsecVtx&
  79 AliHFEsecVtx::operator=(const AliHFEsecVtx &)
  80 {
  81   // Assignment operator
  82
  83   AliInfo("Not yet implemented.");
  84   return *this;
  85 }
  86
  87 //_______________________________________________________________________________________________
  88 AliHFEsecVtx::~AliHFEsecVtx()
  89 {
  90         // Destructor
  91
  92         //cout << "Analysis Done." << endl;
  93 }
  94
  95 //__________________________________________
  96 void AliHFEsecVtx::Init()
  97 {
  98
  99         // set pdg code and index
 100
 101         fNparents = 7;
 102
 103         fParentSelect[0][0] =  411;
 104         fParentSelect[0][1] =  421;
 105         fParentSelect[0][2] =  431;
 106         fParentSelect[0][3] = 4122;
 107         fParentSelect[0][4] = 4132;
 108         fParentSelect[0][5] = 4232;
 109         fParentSelect[0][6] = 4332;
 110
 111         fParentSelect[1][0] =  511;
 112         fParentSelect[1][1] =  521;
 113         fParentSelect[1][2] =  531;
 114         fParentSelect[1][3] = 5122;
 115         fParentSelect[1][4] = 5132;
 116         fParentSelect[1][5] = 5232;
 117         fParentSelect[1][6] = 5332;
 118
 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 }
 166
 167 //__________________________________________
 168 void AliHFEsecVtx::ResetTagVar()
 169 {
 170         // reset tag variables
 171
 172         fdistTwoSecVtx = -1;
 173         fcosPhi = -1;
 174         fsignedLxy = -1;
 175         finvmass = -1;
 176         finvmassSigma = -1;
 177         fBTagged = kFALSE;
 178         fBElec = kFALSE;
 179 }
 180
 181 //__________________________________________
 182 void AliHFEsecVtx::InitAnaPair()
 183 {
 184         // initialize pair tagging variables
 185
 186         fPairTagged = 0;
 187         for (Int_t i=0; i<20; i++){
 188                 fpairedTrackID[i] = -1;
 189                 fpairedChi2[i] = -1;
 190                 fpairedInvMass[i] = -1;
 191                 fpairedSignedLxy[i] = -1;
 192         }
 193
 194         fid[0][0] = 0;
 195         fid[0][1] = 1;
 196         fid[0][2] = 2;
 197
 198         fid[1][0] = 0;
 199         fid[1][1] = 1;
 200         fid[1][2] = 3;
 201
 202         fid[2][0] = 0;
 203         fid[2][1] = 2;
 204         fid[2][2] = 3;
 205
 206         fid[3][0] = 1;
 207         fid[3][1] = 2;
 208         fid[3][2] = 3;
 209
 210         fia[0][0][0] = 0;
 211         fia[0][0][1] = 1;
 212         fia[0][1][0] = 0;
 213         fia[0][1][1] = 2;
 214         fia[0][2][0] = 1;
 215         fia[0][2][1] = 2;
 216
 217         fia[1][0][0] = 0;
 218         fia[1][0][1] = 1;
 219         fia[1][1][0] = 0;
 220         fia[1][1][1] = 3;
 221         fia[1][2][0] = 1;
 222         fia[1][2][1] = 3;
 223
 224         fia[2][0][0] = 0;
 225         fia[2][0][1] = 2;
 226         fia[2][1][0] = 0;
 227         fia[2][1][1] = 3;
 228         fia[2][2][0] = 2;
 229         fia[2][2][1] = 3;
 230
 231         fia[3][0][0] = 1;
 232         fia[3][0][1] = 2;
 233         fia[3][1][0] = 1;
 234         fia[3][1][1] = 3;
 235         fia[3][2][0] = 2;
 236         fia[3][2][1] = 3;
 237
 238 }
 239
 240 //_______________________________________________________________________________________________
 241 void AliHFEsecVtx::CreateHistograms(TString hnopt)
 242 {
 243         // create histograms
 244
 245         fkSourceLabel[kAll]="all";
 246         fkSourceLabel[kDirectCharm]="directCharm";
 247         fkSourceLabel[kDirectBeauty]="directBeauty";
 248         fkSourceLabel[kBeautyCharm]="beauty2charm";
 249         fkSourceLabel[kGamma]="gamma";
 250         fkSourceLabel[kPi0]="pi0";
 251         fkSourceLabel[kElse]="others";
 252         fkSourceLabel[kBeautyGamma]="beauty22gamma";
 253         fkSourceLabel[kBeautyPi0]="beauty22pi0";
 254         fkSourceLabel[kBeautyElse]="beauty22others";
 255
 256
 257         TString hname;
 258         for (Int_t isource = 0; isource < 10; isource++ ){
 259
 260            hname=hnopt+"InvMass_"+fkSourceLabel[isource];
 261            fHistPair[isource].fInvMass = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
 262            hname=hnopt+"InvMassCut1_"+fkSourceLabel[isource];
 263            fHistPair[isource].fInvMassCut1 = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
 264            hname=hnopt+"InvMassCut2_"+fkSourceLabel[isource];
 265            fHistPair[isource].fInvMassCut2 = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
 266            hname=hnopt+"KFChi2_"+fkSourceLabel[isource];
 267            fHistPair[isource].fKFChi2 = new TH1F(hname,hname,200,0,20);
 268            hname=hnopt+"OpenAngle_"+fkSourceLabel[isource];
 269            fHistPair[isource].fOpenAngle = new TH1F(hname,hname,100,0,3.14);
 270            hname=hnopt+"CosOpenAngle_"+fkSourceLabel[isource];
 271            fHistPair[isource].fCosOpenAngle = new TH1F(hname,hname,100,-1.1,1.1);
 272            hname=hnopt+"SignedLxy_"+fkSourceLabel[isource];
 273            fHistPair[isource].fSignedLxy = new TH2F(hname,hname,1000,-5,5,120,-2,10);
 274            hname=hnopt+"KFIP_"+fkSourceLabel[isource];
 275            fHistPair[isource].fKFIP = new TH1F(hname,hname,1000,-5,5);
 276            hname=hnopt+"IPMax_"+fkSourceLabel[isource];
 277            fHistPair[isource].fIPMax= new TH1F(hname,hname,500,0,5);
 278
 279         }
 280
 281         hname=hnopt+"pt_beautyelec";
 282         fHistTagged.fPtBeautyElec= new TH1F(hname,hname,250,0,50);
 283         hname=hnopt+"pt_taggedelec";
 284         fHistTagged.fPtTaggedElec= new TH1F(hname,hname,250,0,50);
 285         hname=hnopt+"pt_righttaggedelec";
 286         fHistTagged.fPtRightTaggedElec = new TH1F(hname,hname,250,0,50);
 287         hname=hnopt+"pt_wrongtaggedelec";
 288         fHistTagged.fPtWrongTaggedElec = new TH1F(hname,hname,250,0,50);
 289
 290         hname=hnopt+"InvmassBeautyElecSecVtx";
 291         fHistTagged.fInvmassBeautyElecSecVtx= new TH1F(hname,hname,120,-2,10);
 292         hname=hnopt+"InvmassNotBeautyElecSecVtx";
 293         fHistTagged.fInvmassNotBeautyElecSecVtx= new TH1F(hname,hname,120,-2,10);
 294
 295         hname=hnopt+"SignedLxyBeautyElecSecVtx";
 296         fHistTagged.fSignedLxyBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
 297         hname=hnopt+"SignedLxyNotBeautyElecSecVtx";
 298         fHistTagged.fSignedLxyNotBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
 299
 300         hname=hnopt+"DistTwoVtxBeautyElecSecVtx";
 301         fHistTagged.fDistTwoVtxBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
 302         hname=hnopt+"DistTwoVtxNotBeautyElecSecVtx";
 303         fHistTagged.fDistTwoVtxNotBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
 304
 305         hname=hnopt+"CosPhiBeautyElecSecVtx";
 306         fHistTagged.fCosPhiBeautyElecSecVtx= new TH1F(hname,hname,100,-1.1,1.1);
 307         hname=hnopt+"CosPhiNotBeautyElecSecVtx";
 308         fHistTagged.fCosPhiNotBeautyElecSecVtx= new TH1F(hname,hname,100,-1.1,1.1);
 309
 310         hname=hnopt+"Chi2BeautyElecSecVtx";
 311         fHistTagged.fChi2BeautyElecSecVtx= new TH1F(hname,hname,200,0,20);
 312         hname=hnopt+"Chi2NotBeautyElecSecVtx";
 313         fHistTagged.fChi2NotBeautyElecSecVtx= new TH1F(hname,hname,200,0,20);
 314
 315         hname=hnopt+"InvmassBeautyElec2trkVtx";
 316         fHistTagged.fInvmassBeautyElec2trkVtx= new TH1F(hname,hname,120,-2,10);
 317         hname=hnopt+"InvmassNotBeautyElec2trkVtx";
 318         fHistTagged.fInvmassNotBeautyElec2trkVtx= new TH1F(hname,hname,120,-2,10);
 319
 320         hname=hnopt+"SignedLxyBeautyElec2trkVtx";
 321         fHistTagged.fSignedLxyBeautyElec2trkVtx= new TH1F(hname,hname,1000,-5,5);
 322         hname=hnopt+"SignedLxyNotBeautyElec2trkVtx";
 323         fHistTagged.fSignedLxyNotBeautyElec2trkVtx= new TH1F(hname,hname,1000,-5,5);
 324
 325         hname=hnopt+"Chi2BeautyElec2trkVtx";
 326         fHistTagged.fChi2BeautyElec2trkVtx= new TH1F(hname,hname,200,0,20);
 327         hname=hnopt+"Chi2NotBeautyElec2trkVtx";
 328         fHistTagged.fChi2NotBeautyElec2trkVtx= new TH1F(hname,hname,200,0,20);
 329 }
 330
 331 //_______________________________________________________________________________________________
 332 void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t index2)
 333 {
 334         // calculate e-h pair characteristics and tag pair
 335
 336         Int_t sourcePart = PairCode(track1,track2);
 337
 338         // get KF particle input pid
 339         Int_t pdg1 = GetMCPID(track1);
 340         Int_t pdg2 = GetMCPID(track2);
 341
 342
 343         // create KF particle of pair
 344         AliKFParticle::SetField(fESD1->GetMagneticField());
 345         AliKFParticle kfTrack1(*track1, pdg1);
 346         AliKFParticle kfTrack2(*track2, pdg2);
 347
 348         AliKFParticle kfSecondary(kfTrack1,kfTrack2);
 349
 350         // copy primary vertex from ESD
 351         AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
 352         if( primVtxCopy.GetNDF() <1 ) return;
 353
 354         //primary vertex point
 355         Double_t pvx = primVtxCopy.GetX();
 356         Double_t pvy = primVtxCopy.GetY();
 357         //Double_t pvz = primVtxCopy.GetZ();
 358
 359         //secondary vertex point from kf particle
 360         Double_t kfx = kfSecondary.GetX();
 361         Double_t kfy = kfSecondary.GetY();
 362         //Double_t kfz = kfSecondary.GetZ();
 363
 364         //momentum at the decay point from kf particle
 365         Double_t kfpx = kfSecondary.GetPx();
 366         Double_t kfpy = kfSecondary.GetPy();
 367         //Double_t kfpz = kfSecondary.GetPz();
 368
 369
 370         Double_t dx = kfx-pvx;
 371         Double_t dy = kfy-pvy;
 372
 373
 374
 375         // discriminating variables ----------------------------------------------------------
 376
 377         // invariant mass of the KF particle
 378         Double_t invmass = -1;
 379         Double_t invmassSigma = -1;
 380         kfSecondary.GetMass(invmass,invmassSigma);
 381
 382         // chi2 of the KF particle
 383         Double_t kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 384
 385         // opening angle between two particles in XY plane
 386         Double_t phi = kfTrack1.GetAngleXY(kfTrack2);
 387         Double_t cosphi = TMath::Cos(phi);
 388
 389         // projection of kf vertex vector to the kf momentum direction
 390         Double_t costheta = ( dx*kfpx + dy*kfpy)/TMath::Sqrt(dx*dx+dy*dy)*TMath::Sqrt(kfpx*kfpx + kfpy*kfpy);
 391         Double_t signedLxy = TMath::Sqrt(dx*dx+dy*dy)*costheta;
 392
 393         // DCA from primary to e-h KF particle (impact parameter of KF particle)
 394         Double_t vtx[2]={pvx, pvy};
 395         Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
 396
 397
 398               Float_t dcaR=-1;
 399         Float_t dcaR1=-1, dcaR2=-1;
 400         Float_t dcaZ1=-1, dcaZ2=-1;
 401         track1->GetImpactParameters(dcaR1,dcaZ1);
 402         track2->GetImpactParameters(dcaR2,dcaZ2);
 403
 404               if (TMath::Abs(dcaR1) >= TMath::Abs(dcaR2)) dcaR=dcaR1;
 405               else dcaR=dcaR2;
 406
 407         // fill histograms
 408         fHistPair[sourcePart].fInvMass->Fill(invmass,invmassSigma);
 409         fHistPair[sourcePart].fKFChi2->Fill(kfchi2);
 410         fHistPair[sourcePart].fOpenAngle->Fill(phi);
 411         fHistPair[sourcePart].fCosOpenAngle->Fill(cosphi);
 412         fHistPair[sourcePart].fSignedLxy->Fill(signedLxy,invmass);
 413         fHistPair[sourcePart].fKFIP->Fill(kfip);
 414         fHistPair[sourcePart].fIPMax->Fill(TMath::Abs(dcaR));
 415
 416         // pair cuts
 417         if( kfchi2 >2. ) return;
 418
 419         if ( signedLxy > 0.05 && cosphi > 0.5) fHistPair[sourcePart].fInvMassCut1->Fill(invmass,invmassSigma);
 420
 421         // pair tagging if it passed the above cuts
 422         if(signedLxy > 0. && cosphi > 0.) fHistPair[sourcePart].fInvMassCut2->Fill(invmass,invmassSigma);
 423
 424
 425         // pair tagging condition
 426         if ( signedLxy > 0.0 && cosphi > 0) { // testing loose cut
 427         //if ( signedLxy > 0.06 && cosphi > 0) {
 428           fpairedTrackID[fPairTagged] = index2;
 429           fpairedChi2[fPairTagged] = kfchi2;
 430           fpairedInvMass[fPairTagged] = invmass;
 431           fpairedSignedLxy[fPairTagged] = signedLxy;
 432           fPairTagged++;
 433         }
 434
 435 }
 436
 437 //_______________________________________________________________________________________________
 438 void AliHFEsecVtx::RunSECVTX(AliESDtrack *track)
 439 {
 440         // run secondary vertexing algorithm and do tagging
 441
 442         ResetTagVar();
 443
 444         Int_t imclabel = TMath::Abs(track->GetLabel());
 445         if(imclabel<0) return;
 446         TParticle* mcpart = fStack->Particle(imclabel);
 447         Int_t esource = GetElectronSource(imclabel);
 448         if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){
 449                 fHistTagged.fPtBeautyElec->Fill(mcpart->Pt());
 450                 fBElec = kTRUE;
 451         }
 452
 453
 454         if (fPairTagged >= 4) {
 455           FindSECVTXCandid4Tracks(track);
 456         }
 457         else if (fPairTagged == 3) {
 458           FindSECVTXCandid3Tracks(track);
 459         }
 460         else if (fPairTagged == 2) {
 461           FindSECVTXCandid2Tracks(track);
 462         }
 463         else if (fPairTagged == 1) {
 464           ApplyPairTagCut();
 465         }
 466
 467
 468         if (fBTagged) {
 469                 fHistTagged.fPtTaggedElec->Fill(mcpart->Pt());
 470                 if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){
 471                         fHistTagged.fPtRightTaggedElec->Fill(mcpart->Pt());
 472                 }
 473                 else fHistTagged.fPtWrongTaggedElec->Fill(mcpart->Pt());
 474         }
 475
 476 }
 477
 478 //_______________________________________________________________________________________________
 479 void AliHFEsecVtx::ApplyPairTagCut()
 480 {
 481         // apply tagging cut for e-h pair
 482
 483         if (fBElec){
 484                 fHistTagged.fInvmassBeautyElec2trkVtx->Fill(fpairedInvMass[0]);
 485                 fHistTagged.fSignedLxyBeautyElec2trkVtx->Fill(fpairedSignedLxy[0]);
 486                 fHistTagged.fChi2BeautyElec2trkVtx->Fill(fpairedChi2[0]);
 487         }
 488         else {
 489                 fHistTagged.fInvmassNotBeautyElec2trkVtx->Fill(fpairedInvMass[0]);
 490                 fHistTagged.fSignedLxyNotBeautyElec2trkVtx->Fill(fpairedSignedLxy[0]);
 491                 fHistTagged.fChi2NotBeautyElec2trkVtx->Fill(fpairedChi2[0]);
 492         }
 493
 494         if (fpairedChi2[0] > 2.0) return;
 495         if (fpairedInvMass[0] < 1.5) return;
 496         if (fpairedSignedLxy[0] < 0.05) return;
 497
 498         fBTagged = kTRUE;
 499 }
 500
 501 //_______________________________________________________________________________________________
 502 Bool_t AliHFEsecVtx::ApplyPairTagCut(Int_t id)
 503 {
 504         // apply tagging cut for e-h pair of given indexed hadron
 505
 506         if (fpairedChi2[id] > 2.0) return kFALSE;
 507         if (fpairedInvMass[id] < 1.5) return kFALSE;
 508         if (fpairedSignedLxy[id] < 0.05) return kFALSE;
 509
 510         fBTagged = kTRUE;
 511         return kTRUE;
 512
 513 }
 514
 515 //_______________________________________________________________________________________________
 516 Bool_t AliHFEsecVtx::ApplyTagCut(Double_t chi2)
 517 {
 518         // apply tagging cut for secondary vertex
 519
 520         if (chi2 > 2.0) return kFALSE;
 521         if (finvmass < 1.5) return kFALSE;
 522         if (fsignedLxy < 0.05) return kFALSE;
 523         if (fcosPhi < 0.90) return kFALSE;
 524         if (fdistTwoSecVtx > 0.1) return kFALSE;
 525
 526         fBTagged = kTRUE;
 527         return kTRUE;
 528 }
 529
 530 //_______________________________________________________________________________________________
 531 void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2, Int_t id1, Int_t id2)
 532 {
 533         // apply tagging cut for e-h pair of given indexed hadron
 534
 535         if(!ApplyTagCut(chi2)){
 536                 if(!ApplyPairTagCut(id1)) ApplyPairTagCut(id2);
 537         }
 538
 539 }
 540
 541 //_______________________________________________________________________________________________
 542 void AliHFEsecVtx::FindSECVTXCandid4Tracks(AliESDtrack *track)
 543 {
 544         // find secondary vertex for >= 4 e-h pairs
 545
 546         // sort pair in increasing order (kFALSE-increasing order)
 547         Int_t index[20];
 548         Int_t indexA[4];
 549         Int_t indexB[3];
 550         Double_t sevchi2[4];
 551         AliESDtrack *htrack[4];
 552
 553         if(fPairTagged>20) return; // protection
 554
 555         TMath::Sort(fPairTagged,fpairedChi2,index,kFALSE);
 556
 557         // select 4 partner tracks retruning smallest pair chi2
 558
 559         for (Int_t i=0; i<4; i++){
 560                 htrack[i] = fESD1->GetTrack(fpairedTrackID[index[i]]);
 561         }
 562
 563         // calculated chi2 with 1 electron and 3 partner tracks
 564         for (Int_t i=0; i<4; i++){
 565                 sevchi2[i] = GetSecVtxChi2(track, htrack[fid[i][0]], htrack[fid[i][1]], htrack[fid[i][2]]);
 566         }
 567
 568         // select 3 partner tracks retruning smallest pair chi2
 569         // [think] if two smallest chi2 are similar, have to think about better handling of selection
 570         TMath::Sort(4,sevchi2,indexA,kFALSE);
 571
 572         // calculated chi2 with 1 electron and 2 partner tracks
 573         for (Int_t i=0; i<3; i++){
 574                 sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
 575         }
 576
 577         // select 2 partner tracks retruning smallest pair chi2
 578         TMath::Sort(3,sevchi2,indexB,kFALSE);
 579
 580         // calculate secondary vertex quality variables with 1 electron and 2 hadrons
 581         CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
 582
 583         if (fBElec){
 584                 fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
 585         }
 586         else {
 587                 fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
 588         }
 589
 590         ApplyVtxTagCut(sevchi2[indexB[0]],index[fia[indexA[0]][indexB[0]][0]],index[fia[indexA[0]][indexB[0]][1]]);
 591
 592 }
 593
 594 //_______________________________________________________________________________________________
 595 void AliHFEsecVtx::FindSECVTXCandid3Tracks(AliESDtrack *track)
 596 {
 597         // find secondary vertex for 3 e-h pairs
 598
 599         // sort pair in increasing order (kFALSE-increasing order)
 600         Int_t indexA[1] = { 0 };
 601         Int_t indexB[3];
 602         Double_t sevchi2[3];
 603         AliESDtrack *htrack[3];
 604
 605         // select 4 partner tracks retruning smallest pair chi2
 606
 607         for (Int_t i=0; i<3; i++){
 608                 htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
 609         }
 610
 611         // calculated chi2 with 1 electron and 2 partner tracks
 612         for (Int_t i=0; i<3; i++){
 613                 sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
 614         }
 615
 616         // select 2 partner tracks retruning smallest pair chi2
 617         TMath::Sort(3,sevchi2,indexB,kFALSE);
 618
 619         // calculate secondary vertex quality variables with 1 electron and 2 hadrons
 620         CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
 621
 622         if (fBElec){
 623                 fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
 624         }
 625         else {
 626                 fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
 627         }
 628
 629         ApplyVtxTagCut(sevchi2[indexB[0]],fia[indexA[0]][indexB[0]][0],fia[indexA[0]][indexB[0]][1]);
 630 }
 631
 632 //_______________________________________________________________________________________________
 633 void AliHFEsecVtx::FindSECVTXCandid2Tracks(AliESDtrack *track)
 634 {
 635         // find secondary vertex for 2 e-h pairs
 636
 637         Double_t sevchi2[1];
 638         AliESDtrack *htrack[2];
 639
 640         for (Int_t i=0; i<2; i++){
 641                 htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
 642         }
 643
 644         sevchi2[0] = GetSecVtxChi2(track, htrack[0], htrack[1]);
 645
 646         // calculate secondary vertex quality variables with 1 electron and 2 hadrons
 647         CalcSECVTXProperty(track,htrack[0],htrack[1]);
 648
 649         if (fBElec){
 650                 fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[0]);
 651         }
 652         else {
 653                 fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[0]);
 654         }
 655
 656         ApplyVtxTagCut(sevchi2[0],0,1);
 657 }
 658
 659 //_______________________________________________________________________________________________
 660 void AliHFEsecVtx::CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
 661 {
 662         // calculate secondary vertex properties
 663
 664         Int_t pdg1 = GetMCPID(track1);
 665         Int_t pdg2 = GetMCPID(track2);
 666         Int_t pdg3 = GetMCPID(track3);
 667
 668         AliKFParticle::SetField(fESD1->GetMagneticField());
 669         AliKFParticle kfTrack1(*track1, pdg1);
 670         AliKFParticle kfTrack2(*track2, pdg2);
 671         AliKFParticle kfTrack3(*track3, pdg3);
 672
 673         AliKFParticle kfSecondary12(kfTrack1,kfTrack2);
 674         AliKFParticle kfSecondary13(kfTrack1,kfTrack3);
 675         AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
 676
 677         // copy primary vertex from ESD
 678         AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
 679         //printf("primary ndf= %d\n",primVtxCopy.GetNDF());
 680         if( primVtxCopy.GetNDF() <1 ) return;
 681
 682         Double_t kdx12 = kfSecondary12.GetX()-primVtxCopy.GetX();
 683         Double_t kdy12 = kfSecondary12.GetY()-primVtxCopy.GetY();
 684         //Double_t kdz12 = kfSecondary12.GetZ()-primVtxCopy.GetZ();
 685
 686         Double_t kdx13 = kfSecondary13.GetX()-primVtxCopy.GetX();
 687         Double_t kdy13 = kfSecondary13.GetY()-primVtxCopy.GetY();
 688         //Double_t kdz13 = kfSecondary13.GetZ()-primVtxCopy.GetZ();
 689
 690         Double_t kdx = kfSecondary.GetX()-primVtxCopy.GetX();
 691         Double_t kdy = kfSecondary.GetY()-primVtxCopy.GetY();
 692         //Double_t kdz = kfSecondary.GetZ()-primVtxCopy.GetZ();
 693
 694         // calculate distance and angle between two secvtxes
 695         fdistTwoSecVtx = TMath::Sqrt((kdx12-kdx13)*(kdx12-kdx13) + (kdy12-kdy13)*(kdy12-kdy13));
 696         fcosPhi = ( kdx12*kdx13 + kdy12*kdy13 ) / ( TMath::Sqrt(kdx12*kdx12+kdy12*kdy12)*TMath::Sqrt(kdx13*kdx13+kdy13*kdy13) );
 697         //Double_t lengthdiff = TMath::Abs(TMath::Sqrt(kdx12*kdx12+kdy12*kdy12) - TMath::Sqrt(kdx13*kdx13+kdy13*kdy13));
 698
 699         // calculate angle between secondary vertex vector and secondary particle momentum vector in transverse plane
 700         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());
 701         // calculate signed Lxy
 702         fsignedLxy = TMath::Sqrt(kdx*kdx+kdy*kdy)*cosTheta;
 703
 704         // calculate invariant mass of the kf particle
 705         kfSecondary.GetMass(finvmass,finvmassSigma);
 706
 707         if (fBElec){
 708                 fHistTagged.fInvmassBeautyElecSecVtx->Fill(finvmass);
 709                 fHistTagged.fSignedLxyBeautyElecSecVtx->Fill(fsignedLxy);
 710                 fHistTagged.fDistTwoVtxBeautyElecSecVtx->Fill(fdistTwoSecVtx);
 711                 fHistTagged.fCosPhiBeautyElecSecVtx->Fill(fcosPhi);
 712         }
 713         else {
 714                 fHistTagged.fInvmassNotBeautyElecSecVtx->Fill(finvmass);
 715                 fHistTagged.fSignedLxyNotBeautyElecSecVtx->Fill(fsignedLxy);
 716                 fHistTagged.fDistTwoVtxNotBeautyElecSecVtx->Fill(fdistTwoSecVtx);
 717                 fHistTagged.fCosPhiNotBeautyElecSecVtx->Fill(fcosPhi);
 718         }
 719
 720 }
 721
 722 //_______________________________________________________________________________________________
 723 Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track)
 724 {
 725         // return mc pid
 726
 727         Int_t label = TMath::Abs(track->GetLabel());
 728         TParticle* mcpart = fStack->Particle(label);
 729         if ( !mcpart ) return 0;
 730         Int_t pdgCode = mcpart->GetPdgCode();
 731
 732         return pdgCode;
 733 }
 734
 735 //_______________________________________________________________________________________________
 736 Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3, AliESDtrack* track4)
 737 {
 738         // return 4 track secondary vertex chi2
 739
 740         Int_t pdg1 = GetMCPID(track1);
 741         Int_t pdg2 = GetMCPID(track2);
 742         Int_t pdg3 = GetMCPID(track3);
 743         Int_t pdg4 = GetMCPID(track4);
 744
 745         AliKFParticle::SetField(fESD1->GetMagneticField());
 746         AliKFParticle kfTrack1(*track1, pdg1);
 747         AliKFParticle kfTrack2(*track2, pdg2);
 748         AliKFParticle kfTrack3(*track3, pdg3);
 749         AliKFParticle kfTrack4(*track4, pdg4);
 750         AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3,kfTrack4);
 751
 752         return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 753
 754 }
 755
 756 //_______________________________________________________________________________________________
 757 Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
 758 {
 759         // return 3 track secondary vertex chi2
 760
 761         Int_t pdg1 = GetMCPID(track1);
 762         Int_t pdg2 = GetMCPID(track2);
 763         Int_t pdg3 = GetMCPID(track3);
 764
 765         AliKFParticle::SetField(fESD1->GetMagneticField());
 766         AliKFParticle kfTrack1(*track1, pdg1);
 767         AliKFParticle kfTrack2(*track2, pdg2);
 768         AliKFParticle kfTrack3(*track3, pdg3);
 769         AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
 770
 771         return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 772
 773 }
 774
 775 //_______________________________________________________________________________________________
 776 Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2)
 777 {
 778         // return 2 track secondary vertex chi2
 779
 780         Int_t pdg1 = GetMCPID(track1);
 781         Int_t pdg2 = GetMCPID(track2);
 782
 783         AliKFParticle::SetField(fESD1->GetMagneticField());
 784         AliKFParticle kfTrack1(*track1, pdg1);
 785         AliKFParticle kfTrack2(*track2, pdg2);
 786         AliKFParticle kfSecondary(kfTrack1,kfTrack2);
 787
 788         return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 789
 790 }
 791
 792 //_______________________________________________________________________________________________
 793 Int_t AliHFEsecVtx::PairOrigin(AliESDtrack* track1, AliESDtrack* track2)
 794 {
 795
 796         //
 797         // return pdg code of the origin(source) of the pair
 798         //
 799         //
 800         // ---*---*---*-----ancester A----- track1
 801         //                        |____*______
 802         //                             |______ track2
 803         // => if they originated from same ancester,
 804         //    then return "the absolute value of pdg code of ancester A"
 805         //
 806         // ---*---*---B hadron-----ancester A----- track1
 807         //                               |____*______
 808         //                                    |______ track2
 809         // => if they originated from same ancester, and this ancester originally comes from B hadrons
 810         //    then return -1*"the absolute value of pdg code of ancester A"
 811         //
 812         // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
 813         //
 814
 815         if (track1->GetLabel()<0 || track2->GetLabel()<0) return 0;
 816         TParticle* part1 = fStack->Particle(TMath::Abs(track1->GetLabel()));
 817         TParticle* part2 = fStack->Particle(TMath::Abs(track2->GetLabel()));
 818         if (!(part1)) return 0;
 819         if (!(part2)) return 0;
 820
 821         TParticle* part1Crtgen = part1; // copy track into current generation particle
 822         TParticle* part2Crtgen = part2; // copy track into current generation particle
 823
 824
 825         Int_t sourcePDG = 0;
 826
 827         // if the two tracks' mother's label is same, get the mother info
 828         // in case of charm, check if it originated from beauty
 829         for (Int_t i=0; i<100; i++){ // iterate 100
 830                 Int_t iLabel = part1Crtgen->GetFirstMother(); //label of mother of current generation for 1st partilce
 831                 if (iLabel < 0) return 0;
 832
 833                 for (Int_t j=0; j<100; j++){ // iterate 100
 834                         Int_t jLabel = part2Crtgen->GetFirstMother(); //label of mother of current generation for 2nd partilce
 835                         if (jLabel < 0) return 0; // if jLabel == -1
 836
 837                         if (iLabel == jLabel){ // check if two tracks are originated from same mother
 838                                 TParticle* thismother = fStack->Particle(jLabel); // if yes, get "thismother" info
 839                                 sourcePDG = abs(thismother->GetPdgCode()); // get the pdg code of "this mother"
 840
 841                                 // check ancester to see if it is originally from beauty
 842                                 for (Int_t k=0; k<10; k++){ // check up to 10 ancesters
 843                                         Int_t ancesterLabel = thismother->GetFirstMother();
 844                                         if (ancesterLabel < 0) return sourcePDG; // if "thismoter" doesn't have mother anymore, return thismother's pdg
 845
 846                                         TParticle* thisancester = fStack->Particle(ancesterLabel);
 847                                         Int_t ancesterPDG = abs(thisancester->GetPdgCode());
 848
 849                                         for (Int_t l=0; l<fNparents; l++){
 850                                                 if (abs(ancesterPDG)==fParentSelect[1][l]){
 851                                                         sourcePDG = -1*sourcePDG; // multiply -1 for charm from bottom
 852                                                         return sourcePDG;
 853                                                 }
 854                                         }
 855                                         thismother = thisancester;
 856                                 }
 857
 858                         }
 859                         part2Crtgen = fStack->Particle(jLabel); // if their mother is different, go up to previous generation of 2nd particle
 860                 }
 861                 part1Crtgen = fStack->Particle(iLabel); // if their mother is different, go up to previous generation of 2nd particle
 862
 863                 // if you don't find additionional(2nd particle) track originated from a given beauty hadron, break to save time
 864                 Int_t motherPDGtmp = abs(part1Crtgen->GetPdgCode());
 865                 for (Int_t l=0; l<fNparents; l++){
 866                         if (abs(motherPDGtmp)==fParentSelect[1][l]) return sourcePDG;
 867                 }
 868
 869         }
 870
 871         return sourcePDG;
 872
 873 }
 874
 875 //_______________________________________________________________________________________________
 876 Int_t AliHFEsecVtx::PairCode(AliESDtrack* track1,AliESDtrack* track2)
 877 {
 878
 879         //
 880         // return pair code which is predefinded as:
 881         //  kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
 882         //
 883
 884         Int_t pairOriginsPDG = PairOrigin(track1,track2);
 885
 886         Int_t sourcePart = kElse;
 887
 888         if (pairOriginsPDG < 0) {
 889                 sourcePart = kBeautyElse;
 890         }
 891         for (Int_t i=0; i<fNparents; i++){
 892                 if (abs(pairOriginsPDG)==fParentSelect[0][i]) {
 893                         if (pairOriginsPDG>0) sourcePart = kDirectCharm;
 894                         if (pairOriginsPDG<0) {
 895                                 sourcePart = kBeautyCharm;
 896                         }
 897                 }
 898                 if (abs(pairOriginsPDG)==fParentSelect[1][i]) {
 899                         if (pairOriginsPDG>0) {
 900                                 sourcePart = kDirectBeauty;
 901                         }
 902                         if (pairOriginsPDG<0)  return kElse;
 903                 }
 904         }
 905         if (pairOriginsPDG == 22) sourcePart = kGamma;
 906         if (pairOriginsPDG == -22) {
 907                 sourcePart = kBeautyGamma;
 908         }
 909         if (pairOriginsPDG == 111) sourcePart = kPi0;
 910         if (pairOriginsPDG == -111) {
 911                 sourcePart = kBeautyPi0;
 912         }
 913
 914         return sourcePart;
 915
 916 }
 917
 918 //_______________________________________________________________________________________________
 919 Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack)
 920 {
 921
 922     // return decay electron's origin
 923
 924     if (iTrack < 0) {
 925       AliDebug(1, "Stack label is negative, return\n");
 926       return -1;
 927     }
 928
 929     TParticle* mcpart = fStack->Particle(iTrack);
 930
 931     if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !!!!!!!!!!!!!!!!!
 932
 933     Int_t iLabel = mcpart->GetFirstMother();
 934     if (iLabel<0){
 935       AliDebug(1, "Stack label is negative, return\n");
 936       return -1;
 937     }
 938
 939     Int_t origin = -1;
 940     Bool_t isFinalOpenCharm = kFALSE;
 941
 942     TParticle *partMother = fStack->Particle(iLabel);
 943     Int_t maPdgcode = partMother->GetPdgCode();
 944
 945     // if the mother is charmed hadron
 946     if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 947
 948          for (Int_t i=0; i<fNparents; i++){
 949             if (abs(maPdgcode)==fParentSelect[0][i]){
 950               isFinalOpenCharm = kTRUE;
 951             }
 952          }
 953          if (!isFinalOpenCharm) return -1;
 954
 955
 956           // iterate until you find B hadron as a mother or become top ancester
 957           for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
 958
 959              Int_t jLabel = partMother->GetFirstMother();
 960              if (jLabel == -1){
 961                origin = kDirectCharm;
 962                return origin;
 963              }
 964              if (jLabel < 0){ // safety protection
 965                AliDebug(1, "Stack label is negative, return\n");
 966                return -1;
 967              }
 968
 969              // if there is an ancester
 970              TParticle* grandMa = fStack->Particle(jLabel);
 971              Int_t grandMaPDG = grandMa->GetPdgCode();
 972
 973              for (Int_t j=0; j<fNparents; j++){
 974                 if (abs(grandMaPDG)==fParentSelect[1][j]){
 975
 976                   origin = kBeautyCharm;
 977                   return origin;
 978                 }
 979              }
 980
 981              partMother = grandMa;
 982           } // end of iteration
 983     } // end of if
 984     else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
 985          for (Int_t i=0; i<fNparents; i++){
 986             if (abs(maPdgcode)==fParentSelect[1][i]){
 987               origin = kDirectBeauty;
 988               return origin;
 989             }
 990          }
 991     } // end of if
 992
 993
 994     //============ gamma ================
 995     else if ( abs(maPdgcode) == 22 ) {
 996          origin = kGamma;
 997
 998           // iterate until you find B hadron as a mother or become top ancester
 999           for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1000
1001              Int_t jLabel = partMother->GetFirstMother();
1002              if (jLabel == -1){
1003                origin = kGamma;
1004                return origin;
1005              }
1006              if (jLabel < 0){ // safety protection
1007                AliDebug(1, "Stack label is negative, return\n");
1008                return -1;
1009              }
1010
1011              // if there is an ancester
1012              TParticle* grandMa = fStack->Particle(jLabel);
1013              Int_t grandMaPDG = grandMa->GetPdgCode();
1014
1015              for (Int_t j=0; j<fNparents; j++){
1016                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1017                   origin = kBeautyGamma;
1018                   return origin;
1019                 }
1020              }
1021
1022              partMother = grandMa;
1023           } // end of iteration
1024
1025          return origin;
1026     } // end of if
1027
1028     //============ pi0 ================
1029     else if ( abs(maPdgcode) == 111 ) {
1030          origin = kPi0;
1031
1032           // iterate until you find B hadron as a mother or become top ancester
1033           for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1034
1035              Int_t jLabel = partMother->GetFirstMother();
1036              if (jLabel == -1){
1037                origin = kPi0;
1038                return origin;
1039              }
1040              if (jLabel < 0){ // safety protection
1041                AliDebug(1, "Stack label is negative, return\n");
1042                return -1;
1043              }
1044
1045              // if there is an ancester
1046              TParticle* grandMa = fStack->Particle(jLabel);
1047              Int_t grandMaPDG = grandMa->GetPdgCode();
1048
1049              for (Int_t j=0; j<fNparents; j++){
1050                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1051                   origin = kBeautyPi0;
1052                   return origin;
1053                 }
1054              }
1055
1056              partMother = grandMa;
1057           } // end of iteration
1058
1059          return origin;
1060     } // end of if
1061
1062
1063     else {
1064         origin = kElse;
1065
1066           // iterate until you find B hadron as a mother or become top ancester
1067           for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1068
1069              Int_t jLabel = partMother->GetFirstMother();
1070              if (jLabel == -1){
1071                origin = kElse;
1072                return origin;
1073              }
1074              if (jLabel < 0){ // safety protection
1075                AliDebug(1, "Stack label is negative, return\n");
1076                return -1;
1077              }
1078
1079              // if there is an ancester
1080              TParticle* grandMa = fStack->Particle(jLabel);
1081              Int_t grandMaPDG = grandMa->GetPdgCode();
1082
1083              for (Int_t j=0; j<fNparents; j++){
1084                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1085                   origin = kBeautyElse;
1086                   return origin;
1087                 }
1088              }
1089
1090              partMother = grandMa;
1091           } // end of iteration
1092
1093     }
1094
1095     return origin;
1096
1097 }
1098
1099 /*
1100 //_______________________________________________________________________________________________
1101 Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
1102 {
1103
1104         //
1105         // decay electron's origin
1106         //
1107
1108         if (iTrackLabel < 0) {
1109                 AliDebug(1, "Stack label is negative, return\n");
1110                 return -1;
1111         }
1112
1113         TParticle* mcpart = fStack->Particle(iTrackLabel);
1114         Int_t iLabel = mcpart->GetFirstMother();
1115         if (iLabel<0){
1116                 AliDebug(1, "Stack label is negative, return\n");
1117                 return -1;
1118         }
1119
1120         Int_t origin = -1;
1121         Bool_t isFinalOpenCharm = kFALSE;
1122
1123         TParticle *partMother = fStack->Particle(iLabel);
1124         Int_t maPdgcode = partMother->GetPdgCode(); // get mother's pdg code
1125
1126         //beauty --------------------------
1127         if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1128                 for (Int_t i=0; i<fNparents; i++){
1129                         if (abs(maPdgcode)==fParentSelect[1][i]){
1130                                 origin = kDirectBeauty;
1131                                 return origin;
1132                         }
1133                         else return -1; // this track is originated beauties not in the final B hadron list => excited beauties
1134                 }
1135         } // end of if
1136
1137         //charm --------------------------
1138         else if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1139
1140                 for (Int_t i=0; i<fNparents; i++){
1141                         if (abs(maPdgcode)==fParentSelect[0][i])
1142                                 isFinalOpenCharm = kTRUE;
1143                 }
1144                 if (!isFinalOpenCharm) return -1; // this track is originated charms not in the final D hadron list => excited charms
1145                                                   // to prevent any possible double counting
1146
1147                 for (Int_t i=0; i<100; i++){ // iterate 100 until you find B hadron as a mother or become top ancester
1148
1149                         Int_t jLabel = partMother->GetFirstMother();
1150                         if (jLabel == -1){
1151                                 origin = kDirectCharm;
1152                                 return origin;
1153                         }
1154                         if (jLabel < 0){ // safety protection even though not really necessary here
1155                                 AliDebug(1, "Stack label is negative, return\n");
1156                                 return -1;
1157                         }
1158
1159                         // if there is an ancester, check if it in the final B hadron list
1160                         TParticle* grandMa = fStack->Particle(jLabel);
1161                         Int_t grandMaPDG = grandMa->GetPdgCode();
1162
1163                         for (Int_t j=0; j<fNparents; j++){
1164                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1165                                 origin = kBeautyCharm;
1166                                 return origin;
1167                                 }
1168                         }
1169
1170                         partMother = grandMa;
1171                 } // end of iteration
1172         } // end of if
1173
1174         //gamma --------------------------
1175         else if ( abs(maPdgcode) == 22 ) {
1176                 origin = kGamma;
1177
1178                 // iterate until you find B hadron as a mother or become top ancester
1179                 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
1180
1181                         Int_t jLabel = partMother->GetFirstMother();
1182                         if (jLabel == -1){
1183                                 origin = kGamma;
1184                                 return origin;
1185                         }
1186                         if (jLabel < 0){ // safety protection
1187                                 AliDebug(1, "Stack label is negative, return\n");
1188                                 return -1;
1189                         }
1190
1191                         // if there is an ancester
1192                         TParticle* grandMa = fStack->Particle(jLabel);
1193                         Int_t grandMaPDG = grandMa->GetPdgCode();
1194
1195                         for (Int_t j=0; j<fNparents; j++){
1196                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1197                                         origin = kBeautyGamma;
1198                                         return origin;
1199                                 }
1200                         }
1201
1202                         partMother = grandMa;
1203                 } // end of iteration
1204
1205                 return origin;
1206         } // end of if
1207
1208         //pi0 --------------------------
1209         else if ( abs(maPdgcode) == 111 ) {
1210                 origin = kPi0;
1211
1212                 // iterate until you find B hadron as a mother or become top ancester
1213                 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
1214
1215                         Int_t jLabel = partMother->GetFirstMother();
1216                         if (jLabel == -1){
1217                                 origin = kPi0;
1218                                 return origin;
1219                         }
1220                         if (jLabel < 0){ // safety protection
1221                                 AliDebug(1, "Stack label is negative, return\n");
1222                                 return -1;
1223                         }
1224
1225                         // if there is an ancester
1226                         TParticle* grandMa = fStack->Particle(jLabel);
1227                         Int_t grandMaPDG = grandMa->GetPdgCode();
1228
1229                         for (Int_t j=0; j<fNparents; j++){
1230                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1231                                         origin = kBeautyPi0;
1232                                         return origin;
1233                                 }
1234                         }
1235
1236                         partMother = grandMa;
1237                 } // end of iteration
1238
1239                 return origin;
1240         } // end of if
1241
1242
1243         //else --------------------------
1244         else {
1245                 origin = kElse;
1246
1247                 // iterate until you find B hadron as a mother or become top ancester
1248                 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
1249
1250                         Int_t jLabel = partMother->GetFirstMother();
1251                         if (jLabel == -1){
1252                                 origin = kElse;
1253                                 return origin;
1254                         }
1255                         if (jLabel < 0){ // safety protection
1256                                 AliDebug(1, "Stack label is negative, return\n");
1257                                 return -1;
1258                         }
1259
1260                         // if there is an ancester
1261                         TParticle* grandMa = fStack->Particle(jLabel);
1262                         Int_t grandMaPDG = grandMa->GetPdgCode();
1263
1264                         for (Int_t j=0; j<fNparents; j++){
1265                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1266                                         origin = kBeautyElse;
1267                                         return origin;
1268                                 }
1269                         }
1270
1271                         partMother = grandMa;
1272                 } // end of iteration
1273
1274         }
1275
1276         return origin;
1277
1278 }
1279 */
1280
1281 //_______________________________________________________________________________________________
1282 Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track) const
1283 {
1284         // test cuts
1285
1286         //if (track->Pt() < 1.0) return kFALSE;
1287         //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
1288         //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
1289         //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1290         if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
1291         //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
1292
1293
1294 /*
1295         Float_t dcaR=-1;
1296         Float_t dcaZ=-1;
1297         track->GetImpactParameters(dcaR,dcaZ);
1298         if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
1299         if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
1300 */
1301         return kTRUE;
1302 }