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