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 <TIterator.h>
  26 #include <TParticle.h>
  27
  28 #include <AliESDVertex.h>
  29 #include <AliESDEvent.h>
  30 #include <AliAODEvent.h>
  31 #include <AliVTrack.h>
  32 #include <AliESDtrack.h>
  33 #include <AliAODTrack.h>
  34 #include "AliHFEsecVtx.h"
  35 #include <AliKFParticle.h>
  36 #include <AliKFVertex.h>
  37 #include <AliLog.h>
  38 #include <AliMCEvent.h>
  39 #include <AliAODMCParticle.h>
  40 #include "AliHFEpairs.h"
  41 #include "AliHFEsecVtxs.h"
  42 #include "AliHFEtrackFilter.h"
  43 #include "AliHFEmcQA.h"
  44 #include "AliHFEtools.h"
  45
  46 ClassImp(AliHFEsecVtx)
  47
  48 //_______________________________________________________________________________________________
  49 AliHFEsecVtx::AliHFEsecVtx():
  50   fFilter(0x0)
  51   ,fESD1(0x0)
  52   ,fAOD1(0x0)
  53   ,fMCEvent(0x0)
  54   ,fMCQA(0x0)
  55   ,fUseMCPID(kFALSE)
  56   ,fkSourceLabel()
  57   ,fNparents(0)
  58   ,fParentSelect()
  59   ,fPtRng()
  60   ,fDcaCut()
  61   ,fNoOfHFEpairs(0)
  62   ,fNoOfHFEsecvtxs(0)
  63   ,fArethereSecVtx(0)
  64   ,fHFEpairs(0x0)
  65   ,fHFEsecvtxs(0x0)
  66   ,fMCArray(0x0)
  67   ,fPVx(-999)
  68   ,fPVy(-999)
  69   ,fPVx2(999)
  70   ,fPVy2(-999)
  71   ,fCosPhi(-1)
  72   ,fSignedLxy(-1)
  73   ,fSignedLxy2(-1)
  74   ,fKFchi2(-1)
  75   ,fInvmass(-1)
  76   ,fInvmassSigma(-1)
  77   ,fKFip(0)
  78   ,fKFip2(0)
  79   ,fNsectrk2prim(0)
  80   ,fVtxchi2Tightcut(3.)
  81   ,fVtxchi2Loosecut(5.)
  82   ,fPairQA(0x0)
  83   ,fSecvtxQA(0x0)
  84   ,fSecVtxList(0x0)
  85 {
  86   //
  87   // Default constructor
  88   //
  89
  90   Init();
  91 }
  92
  93 //_______________________________________________________________________________________________
  94 AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
  95   TObject(p)
  96   ,fFilter(0x0)
  97   ,fESD1(0x0)
  98   ,fAOD1(0x0)
  99   ,fMCEvent(0x0)
 100   ,fMCQA(0x0)
 101   ,fUseMCPID(p.fUseMCPID)
 102   ,fkSourceLabel()
 103   ,fNparents(p.fNparents)
 104   ,fParentSelect()
 105   ,fPtRng()
 106   ,fDcaCut()
 107   ,fNoOfHFEpairs(p.fNoOfHFEpairs)
 108   ,fNoOfHFEsecvtxs(p.fNoOfHFEsecvtxs)
 109   ,fArethereSecVtx(p.fArethereSecVtx)
 110   ,fHFEpairs(0x0)
 111   ,fHFEsecvtxs(0x0)
 112   ,fMCArray(0x0)
 113   ,fPVx(p.fPVx)
 114   ,fPVy(p.fPVy)
 115   ,fPVx2(p.fPVx2)
 116   ,fPVy2(p.fPVy2)
 117   ,fCosPhi(p.fCosPhi)
 118   ,fSignedLxy(p.fSignedLxy)
 119   ,fSignedLxy2(p.fSignedLxy2)
 120   ,fKFchi2(p.fKFchi2)
 121   ,fInvmass(p.fInvmass)
 122   ,fInvmassSigma(p.fInvmassSigma)
 123   ,fKFip(p.fKFip)
 124   ,fKFip2(p.fKFip2)
 125   ,fNsectrk2prim(p.fNsectrk2prim)
 126   ,fVtxchi2Tightcut(p.fVtxchi2Tightcut)
 127   ,fVtxchi2Loosecut(p.fVtxchi2Loosecut)
 128   ,fPairQA(0x0)
 129   ,fSecvtxQA(0x0)
 130   ,fSecVtxList(0x0)
 131 {
 132   //
 133   // Copy constructor
 134   //
 135   fFilter = new AliHFEtrackFilter(*p.fFilter);
 136 }
 137
 138 //_______________________________________________________________________________________________
 139 AliHFEsecVtx&
 140 AliHFEsecVtx::operator=(const AliHFEsecVtx &)
 141 {
 142   //
 143   // Assignment operator
 144   //
 145
 146   AliInfo("Not yet implemented.");
 147   return *this;
 148 }
 149
 150 //_______________________________________________________________________________________________
 151 AliHFEsecVtx::~AliHFEsecVtx()
 152 {
 153   //
 154   // Destructor
 155   //
 156
 157   //cout << "Analysis Done." << endl;
 158   delete fFilter;
 159 }
 160
 161 //__________________________________________
 162 void AliHFEsecVtx::Init()
 163 {
 164   //
 165   // set pdg code and index
 166   //
 167
 168   fNparents = 7;
 169
 170   fParentSelect[0][0] =  411;
 171   fParentSelect[0][1] =  421;
 172   fParentSelect[0][2] =  431;
 173   fParentSelect[0][3] = 4122;
 174   fParentSelect[0][4] = 4132;
 175   fParentSelect[0][5] = 4232;
 176   fParentSelect[0][6] = 4332;
 177
 178   fParentSelect[1][0] =  511;
 179   fParentSelect[1][1] =  521;
 180   fParentSelect[1][2] =  531;
 181   fParentSelect[1][3] = 5122;
 182   fParentSelect[1][4] = 5132;
 183   fParentSelect[1][5] = 5232;
 184   fParentSelect[1][6] = 5332;
 185
 186   // momentum ranges to apply pt dependent cuts
 187   fPtRng[0] = 1.0;
 188   fPtRng[1] = 2.0;
 189   fPtRng[2] = 2.5;
 190   fPtRng[3] = 3.0;
 191   fPtRng[4] = 5.0;
 192   fPtRng[5] = 12.0;
 193   fPtRng[6] = 20.0;
 194
 195   // momentum dependent DCA cut values (preliminary)
 196   fDcaCut[0] = 0.04;
 197   fDcaCut[1] = 0.03;
 198   fDcaCut[2] = 0.02;
 199   fDcaCut[3] = 0.015;
 200   fDcaCut[4] = 0.01;
 201   fDcaCut[5] = 0.005;
 202
 203   fFilter = new AliHFEtrackFilter("SecVtxFilter");
 204   fFilter->MakeCutStepRecKineITSTPC();
 205   fFilter->MakeCutStepPrimary();
 206 }
 207
 208 void AliHFEsecVtx::Process(AliVTrack *signalTrack){
 209   //
 210   // Run Process
 211   //
 212   //if(signalTrack->Pt() < 1.0) return;
 213
 214
 215   AliESDtrack *track = dynamic_cast<AliESDtrack *>(signalTrack);
 216
 217   FillHistos(0,track); // wo any cuts
 218
 219   InitHFEpairs();
 220   InitHFEsecvtxs();
 221   AliESDtrack *htrack = 0x0;
 222   fFilter->Flush();
 223   fFilter->FilterTracks(fESD1);
 224   TObjArray *candidates = fFilter->GetFilteredTracks();
 225   TIterator *trackIter = candidates->MakeIterator();
 226   while((htrack = dynamic_cast<AliESDtrack *>(trackIter->Next()))){
 227     if(track->GetID() == htrack->GetID()) continue; // since it is for tagging single electron, don't need additional condition
 228     if (htrack->Pt()<1.0) continue;
 229     PairAnalysis(track, htrack, htrack->GetID()); // e-h pairing
 230   }
 231   delete trackIter;
 232   for(int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
 233     //if(HasMCData()){
 234     AliHFEpairs *pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
 235     //if(!(pair->GetPairCode()>1. && pair->GetPairCode()<4.))  // apply various cuts
 236     if(pair->GetKFChi2()>5.) // only apply vertex chi2 cut for the moment
 237       HFEpairs()->RemoveAt(ip);
 238     //}
 239   }
 240   HFEpairs()->Compress();
 241   if(HFEpairs()->GetEntriesFast()) FillHistos(1,track); // after paired
 242   if(HFEpairs()->GetEntriesFast()) RunSECVTX(track); // secondary vertexing with e,h1,h2,.. tracks
 243   if(HFEsecvtxs()->GetEntriesFast()) FillHistos(2,track); // after secvtxing
 244   for(int ip=0; ip<HFEsecvtxs()->GetEntriesFast(); ip++){
 245     AliHFEsecVtxs *secvtx=0x0;
 246     secvtx = (AliHFEsecVtxs*) (HFEsecvtxs()->UncheckedAt(ip));
 247     // here you apply cuts, then if it doesn't pass the cut, remove it from the HFEsecvtxs()
 248     if(!(secvtx->GetInvmass()>2.0 && secvtx->GetInvmass()<5.2) || !(secvtx->GetSignedLxy2()>0.08 && secvtx->GetSignedLxy2()<1.5) || !(secvtx->GetKFIP2()>-0.1 && secvtx->GetKFIP2()<0.1))
 249       HFEsecvtxs()->RemoveAt(ip);
 250   }
 251
 252   // fill histos for raw spectra
 253   if(HFEsecvtxs()->GetEntriesFast()) FillHistos(3,track); //after secvtx cut
 254
 255   DeleteHFEpairs();
 256   DeleteHFEsecvtxs();
 257 }
 258
 259 //_______________________________________________________________________________________________
 260 void AliHFEsecVtx::CreateHistograms(TList * const qaList)
 261 {
 262   //
 263   // create histograms
 264   //
 265
 266   if(!qaList) return;
 267
 268   fSecVtxList = qaList;
 269   fSecVtxList->SetName("SecVtx");
 270
 271   MakeContainer();
 272   MakeHistos(0);
 273   MakeHistos(1);
 274   MakeHistos(2);
 275   MakeHistos(3);
 276
 277   /*
 278   fkSourceLabel[kAll]="all";
 279   fkSourceLabel[kDirectCharm]="directCharm";
 280   fkSourceLabel[kDirectBeauty]="directBeauty";
 281   fkSourceLabel[kBeautyCharm]="beauty2charm";
 282   fkSourceLabel[kGamma]="gamma";
 283   fkSourceLabel[kPi0]="pi0";
 284   fkSourceLabel[kElse]="others";
 285   fkSourceLabel[kBeautyGamma]="beauty22gamma";
 286   fkSourceLabel[kBeautyPi0]="beauty22pi0";
 287   fkSourceLabel[kBeautyElse]="beauty22others";
 288
 289   TString hname;
 290   TString hnopt="secvtx_";
 291   for (Int_t isource = 0; isource < 10; isource++ ){
 292   }
 293   */
 294
 295   //qaList->AddLast(fSecVtxList);
 296 }
 297
 298 //_______________________________________________________________________________________________
 299 void AliHFEsecVtx::GetPrimaryCondition()
 300 {
 301   //
 302   // get primary characteristics and set
 303   //
 304
 305   if (fESD1) {
 306     AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
 307     if( primVtxCopy.GetNDF() <1 ) return;
 308     fPVx = primVtxCopy.GetX();
 309     fPVy = primVtxCopy.GetY();
 310   }
 311   else if(fAOD1) {
 312     AliKFVertex primVtxCopy(*(fAOD1->GetPrimaryVertex()));
 313     if( primVtxCopy.GetNDF() <1 ) return;
 314     fPVx = primVtxCopy.GetX();
 315     fPVy = primVtxCopy.GetY();
 316   }
 317 }
 318
 319 //_______________________________________________________________________________________________
 320 void AliHFEsecVtx::PairAnalysis(AliVTrack* track1, AliVTrack* track2, Int_t index2)
 321 {
 322   //
 323   // calculate e-h pair characteristics and tag pair
 324   //
 325
 326   //later consider the below
 327   Float_t dca1[2]={-999.,-999.}, dca2[2]={-999.,-999.};
 328   Float_t cov1[3]={-999.,-999.,-999.}, cov2[3]={-999.,-999.,-999.};
 329
 330   Double_t dca1aod[2]={-999.,-999.}, dca2aod[2]={-999.,-999.};
 331   Double_t cov1aod[3]={-999.,-999.,-999.}, cov2aod[3]={-999.,-999.,-999.};
 332
 333   if (IsAODanalysis()){
 334     const AliAODVertex *primVtx = fAOD1->GetPrimaryVertex();
 335     AliESDtrack esdTrk1(track1);
 336     AliESDtrack esdTrk2(track2);
 337     esdTrk1.PropagateToDCA(primVtx,0.,10000.,dca1aod,cov1aod);
 338     esdTrk2.PropagateToDCA(primVtx,0.,10000.,dca2aod,cov2aod);
 339   }
 340   else {
 341     ((AliESDtrack*)track1)->GetImpactParameters(dca1,cov1);
 342     ((AliESDtrack*)track2)->GetImpactParameters(dca2,cov2);
 343   }
 344
 345   // apply pt dependent dca cut on hadrons
 346   /*for(int ibin=0; ibin<6; ibin++){
 347     if((track2->Pt()>fPtRng[ibin] && track2->Pt()<fPtRng[ibin+1]) && TMath::Abs(dca2[0])<fDcaCut[ibin]) return;
 348   }*/
 349
 350   // get KF particle input pid
 351   Int_t pdg1 = GetPDG(track1);
 352   Int_t pdg2 = GetPDG(track2);
 353
 354   if(pdg1==-1 || pdg2==-1) {
 355     //printf("out if considered pid range \n");
 356     return;
 357   }
 358
 359   // create KF particle of pair
 360   if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
 361   else AliKFParticle::SetField(fESD1->GetMagneticField());
 362
 363   AliKFParticle kfTrack[2];
 364   kfTrack[0] = AliKFParticle(*track1, pdg1);
 365   kfTrack[1] = AliKFParticle(*track2, pdg2);
 366
 367   AliKFParticle kfSecondary(kfTrack[0],kfTrack[1]);
 368
 369   //secondary vertex point from kf particle
 370   Double_t kfx = kfSecondary.GetX();
 371   Double_t kfy = kfSecondary.GetY();
 372   //Double_t kfz = kfSecondary.GetZ();
 373
 374   //momentum at the decay point from kf particle
 375   Double_t kfpx = kfSecondary.GetPx();
 376   Double_t kfpy = kfSecondary.GetPy();
 377   //Double_t kfpz = kfSecondary.GetPz();
 378
 379
 380 /* //directly use of ESD vertex
 381   const AliESDVertex *pvertex = fESD1->GetPrimaryVertex();
 382   Double_t xyzVtx[3];
 383   pvertex->GetXYZ(xyzVtx);
 384
 385   Double_t dx = kfx-xyzVtx[0];
 386   Double_t dy = kfy-xyzVtx[1];*/
 387
 388   AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
 389   if( primVtxCopy.GetNDF() <1 ) return;
 390   fPVx = primVtxCopy.GetX();
 391   fPVy = primVtxCopy.GetY();
 392
 393  // printf("esdx= %lf kfx= %lf esdy= %lf kfy= %lf\n",xyzVtx[0],fPVx,xyzVtx[1],fPVy);
 394
 395   Double_t dx = kfx-fPVx;
 396   Double_t dy = kfy-fPVy;
 397
 398   // discriminating variables
 399   // invariant mass of the KF particle
 400   Double_t invmass = -1;
 401   Double_t invmassSigma = -1;
 402   kfSecondary.GetMass(invmass,invmassSigma);
 403
 404   // chi2 of the KF particle
 405   Double_t kfchi2 = -1;
 406   if(kfSecondary.GetNDF()>0) kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 407
 408   // opening angle between two particles in XY plane
 409   Double_t phi = kfTrack[0].GetAngleXY(kfTrack[1]);
 410   Double_t cosphi = TMath::Cos(phi);
 411
 412   // DCA from primary to e-h KF particle (impact parameter of KF particle)
 413   Double_t vtx[2]={fPVx, fPVy};
 414   Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
 415
 416   // projection of kf vertex vector to the kf momentum direction
 417   Double_t signedLxy=-999.;
 418   if((dx*kfpx+dy*kfpy)>0) signedLxy = TMath::Sqrt(dx*dx+dy*dy);
 419   if((dx*kfpx+dy*kfpy)<0) signedLxy = -1*TMath::Sqrt(dx*dx+dy*dy);
 420   //[the other way to think about]
 421   //Double_t psqr = kfpx*kfpx+kfpy*kfpy;
 422   //if(psqr>0) signedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
 423
 424   //recalculating primary vertex after removing secvtx tracks --------------------------
 425   Int_t trkid[2];
 426   trkid[0] = track1->GetID();
 427   trkid[1] = track2->GetID();
 428
 429   RecalcPrimvtx(2, trkid, kfTrack);
 430   Double_t dx2 = kfx-fPVx2;
 431   Double_t dy2 = kfy-fPVy2;
 432
 433   // IP of sec particle recalculated based on recalculated primary vertex
 434   Double_t vtx2[2]={fPVx2, fPVy2};
 435   Double_t kfip2 = kfSecondary.GetDistanceFromVertexXY(vtx2);
 436   // signed Lxy recalculated based on recalculated primary vertex
 437   Double_t signedLxy2=-999.;
 438   if((dx2*kfpx+dy2*kfpy)>0) signedLxy2 = TMath::Sqrt(dx2*dx2+dy2*dy2);
 439   if((dx2*kfpx+dy2*kfpy)<0) signedLxy2 = -1*TMath::Sqrt(dx2*dx2+dy2*dy2);
 440   //------------------------------------------------------------------------------------
 441
 442
 443   Int_t paircode = -1;
 444   if (HasMCData()) paircode = GetPairCode(track1,track2);
 445
 446   AliHFEpairs hfepair;
 447   hfepair.SetTrkLabel(index2);
 448   hfepair.SetInvmass(invmass);
 449   hfepair.SetKFChi2(kfchi2);
 450   hfepair.SetOpenangle(phi);
 451   hfepair.SetCosOpenangle(cosphi);
 452   hfepair.SetSignedLxy(signedLxy);
 453   hfepair.SetSignedLxy2(signedLxy2);
 454   hfepair.SetKFIP(kfip);
 455   hfepair.SetKFIP2(kfip2);
 456   hfepair.SetPairCode(paircode);
 457   AddHFEpairToArray(&hfepair);
 458   fNoOfHFEpairs++;
 459
 460   // fill into container for later QA
 461   Double_t dataE[6];
 462   dataE[0]=invmass;
 463   dataE[1]=kfchi2;
 464   dataE[2]=phi;
 465   dataE[3]=signedLxy2;
 466   dataE[4]=kfip2;
 467   dataE[5]=paircode;
 468   //if(cov1[0]>0) dataE[6]=Double_t(dca1[0]/cov1[0]);
 469   //if(cov2[0]>0) dataE[7]=Double_t(dca2[0]/cov2[0]);
 470   //dataE[6]=track1->Pt();
 471   //dataE[7]=track2->Pt();
 472   //dataE[6]=dca1[0]; //mjtmp
 473   //dataE[7]=dca2[0]; //mjtmp
 474   //dataE[8]=TMath::Abs(dca1[0]);
 475   //dataE[9]=TMath::Abs(dca2[0]);
 476
 477   fPairQA->Fill(dataE);
 478
 479 }
 480
 481 //_______________________________________________________________________________________________
 482 void AliHFEsecVtx::RunSECVTX(AliVTrack *track)
 483 {
 484   //
 485   // run secondary vertexing algorithm and do tagging
 486   //
 487
 488   //printf("number of considered pairs= %d\n",HFEpairs()->GetEntriesFast());
 489   FindSECVTXCandid(track);
 490 }
 491
 492 //_______________________________________________________________________________________________
 493 void AliHFEsecVtx::FindSECVTXCandid(AliVTrack *track)
 494 {
 495   //
 496   // find secondary vertex candidate and store them into container
 497   //
 498
 499   AliVTrack *htrack[20];
 500   //Int_t htracklabel[20];
 501   //Int_t paircode[20];
 502   //Double_t vtxchi2[20];
 503   //Double_t dataE[7]={-999.,-999.,-999.,-999.,-1.,0,0};
 504
 505   fVtxchi2Tightcut=3.; // tight cut for pair
 506   fVtxchi2Loosecut=5.; // loose cut for secvtx
 507
 508   if (HFEpairs()->GetEntriesFast()>20){
 509     AliDebug(3, "number of paired hadron is over maximum(20)");
 510     return;
 511   }
 512
 513   // get paired track objects
 514   AliHFEpairs *pair=0x0;
 515   if (IsAODanalysis()){
 516     for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
 517        pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
 518        //htracklabel[ip] = pair->GetTrkLabel();
 519        htrack[ip] = fAOD1->GetTrack(pair->GetTrkLabel());
 520        //if(pair->GetPairCode()==2 || pair->GetPairCode()==3) paircode[ip]=1;
 521        //else paircode[ip]=0;
 522        //vtxchi2[ip] = pair->GetKFChi2();
 523     }
 524   }
 525   else{
 526     for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
 527        pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
 528        //htracklabel[ip] = pair->GetTrkLabel();
 529        htrack[ip] = fESD1->GetTrack(pair->GetTrkLabel());
 530        //if(pair->GetPairCode()==2 || pair->GetPairCode()==3) paircode[ip]=1;
 531        //else paircode[ip]=0;
 532        //vtxchi2[ip] = pair->GetKFChi2();
 533     }
 534   }
 535
 536   Int_t nPairs = HFEpairs()->GetEntriesFast();
 537
 538
 539   // 1 electron candidate + 1 track
 540   if (nPairs == 1){
 541     if (pair->GetKFChi2() < fVtxchi2Tightcut) { // you can also put single track cut -> here you apply very tight cut for the pair
 542       Fill2TrkSECVTX(track, pair);
 543     }
 544     return;
 545   }
 546   //--------------------------------------------------------------
 547
 548   // 1 electron candidate + 2 tracks
 549   if (nPairs == 2){
 550     CalcSECVTXProperty(track, htrack[0], htrack[1]); // calculate secondary vertex property
 551
 552     if (fKFchi2 < fVtxchi2Loosecut) { // -> here you apply rather loose cut
 553       Fill3TrkSECVTX(track, 0, 1);
 554     }
 555     else{ // if doesn't pass the sec vtx chi2 cut
 556       for(int jp=0; jp<2; jp++){
 557          pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(jp));
 558          if (pair->GetKFChi2() < fVtxchi2Tightcut){
 559            Fill2TrkSECVTX(track, pair);
 560          }
 561       }
 562     }
 563     return;
 564   }
 565   //--------------------------------------------------------------
 566
 567   // 1 electron candidate + 3 tracks
 568   if (nPairs == 3){
 569     CalcSECVTXProperty(track, htrack[0], htrack[1], htrack[2]); // calculate secondary vertex property
 570
 571     if (fKFchi2 < fVtxchi2Loosecut) {
 572       Fill4TrkSECVTX(track, 0, 1, 2);
 573     }
 574     else {
 575       fArethereSecVtx=0;
 576       for (int i=0; i<nPairs-1; i++){
 577          for (int j=i+1; j<nPairs; j++){
 578             CalcSECVTXProperty(track, htrack[i], htrack[j]);
 579             if (fKFchi2 < fVtxchi2Loosecut) {
 580               fArethereSecVtx++;
 581               Fill3TrkSECVTX(track, i, j);
 582             }
 583          }
 584       }
 585       if(!fArethereSecVtx){
 586         for(int jp=0; jp<nPairs; jp++){
 587            pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(jp));
 588            if (pair->GetKFChi2() < fVtxchi2Tightcut){
 589              Fill2TrkSECVTX(track, pair);
 590            }
 591         }
 592       }
 593     }
 594     return;
 595   }
 596   //--------------------------------------------------------------
 597
 598   // 1 electron candidate + more than 3 tracks
 599   if (nPairs > 3){
 600     fArethereSecVtx=0;
 601     for (int ih1=0; ih1<nPairs-2; ih1++){
 602       for (int ih2=ih1+1; ih2<nPairs-1; ih2++){
 603         for (int ih3=ih2+1; ih3<nPairs; ih3++){
 604           CalcSECVTXProperty(track, htrack[ih1], htrack[ih2], htrack[ih3]); // calculate secondary vertex property
 605           if (fKFchi2 < fVtxchi2Loosecut) {
 606             fArethereSecVtx++;
 607             Fill4TrkSECVTX(track, ih1, ih2, ih3);
 608           }
 609         }
 610       }
 611     }
 612     if (!fArethereSecVtx){
 613       fArethereSecVtx=0;
 614       for (int i=0; i<nPairs-1; i++){
 615         for (int j=i+1; j<nPairs; j++){
 616           CalcSECVTXProperty(track, htrack[i], htrack[j]);
 617           if (fKFchi2 < fVtxchi2Loosecut) {
 618             fArethereSecVtx++;
 619             Fill3TrkSECVTX(track, i, j);
 620           }
 621         }
 622       }
 623     }
 624     if (!fArethereSecVtx){
 625       for(int jp=0; jp<nPairs; jp++){
 626         pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(jp));
 627         if (pair->GetKFChi2() < fVtxchi2Tightcut){
 628           Fill2TrkSECVTX(track, pair);
 629         }
 630       }
 631     }
 632     return;
 633   }
 634   //--------------------------------------------------------------
 635
 636 }
 637
 638 //_______________________________________________________________________________________________
 639 void AliHFEsecVtx::Fill4TrkSECVTX(AliVTrack* track, Int_t ipair, Int_t jpair, Int_t kpair)
 640 {
 641   //
 642   // fill 3 tracks' secondary vertex properties
 643   //
 644
 645   Double_t dataE[9]={-999.,-999.,-999.,-999.,-1.,0,0,-999.,-999.};
 646
 647   Int_t paircode1 = 0, paircode2 = 0, paircode3 = 0;
 648   Int_t htracklabel1 = 0, htracklabel2= 0;
 649
 650   if (HasMCData()){
 651     AliHFEpairs *pair1=0x0;
 652     AliHFEpairs *pair2=0x0;
 653     AliHFEpairs *pair3=0x0;
 654     pair1 = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ipair));
 655     pair2 = (AliHFEpairs*) (HFEpairs()->UncheckedAt(jpair));
 656     pair3 = (AliHFEpairs*) (HFEpairs()->UncheckedAt(kpair));
 657
 658     htracklabel1 = pair1->GetTrkLabel();
 659     htracklabel2 = pair2->GetTrkLabel();
 660
 661     if (pair1->GetPairCode()==2 || pair1->GetPairCode()==3) paircode1=1;
 662     else paircode1=0;
 663     if (pair2->GetPairCode()==2 || pair2->GetPairCode()==3) paircode2=1;
 664     else paircode2=0;
 665     if (pair3->GetPairCode()==2 || pair3->GetPairCode()==3) paircode3=1;
 666     else paircode3=0;
 667   }
 668
 669   AliHFEsecVtxs hfesecvtx;
 670   hfesecvtx.SetTrkLabel1(htracklabel1); // mj: not much meaningful for the moment
 671   hfesecvtx.SetTrkLabel2(htracklabel2); // mj: not much meaningful for the moment
 672   if(HasMCData()) hfesecvtx.SetMCCode(GetElectronSource(TMath::Abs(track->GetLabel())));
 673   hfesecvtx.SetKFChi2(fKFchi2);
 674   hfesecvtx.SetInvmass(fInvmass);
 675   hfesecvtx.SetSignedLxy(fSignedLxy);
 676   hfesecvtx.SetSignedLxy2(fSignedLxy2);
 677   hfesecvtx.SetKFIP(fKFip);
 678   hfesecvtx.SetKFIP2(fKFip2);
 679   AddHFEsecvtxToArray(&hfesecvtx);
 680   fNoOfHFEsecvtxs++;
 681
 682   dataE[0]=fInvmass;
 683   dataE[1]=fKFchi2;
 684   dataE[2]=fSignedLxy;
 685   dataE[3]=fKFip;
 686   if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
 687   dataE[5]=4; //# of associated tracks
 688   if(paircode1 & paircode2 & paircode3) dataE[6]=1;
 689   else if(!paircode1 & !paircode2 & !paircode3) dataE[6]=0;
 690   else dataE[6]=3;
 691   dataE[7]=fSignedLxy2;
 692   dataE[8]=track->Pt();
 693   fSecvtxQA->Fill(dataE);
 694 }
 695
 696 //_______________________________________________________________________________________________
 697 void AliHFEsecVtx::Fill3TrkSECVTX(AliVTrack* track, Int_t ipair, Int_t jpair)
 698 {
 699   //
 700   // fill 3 tracks' secondary vertex properties
 701   //
 702
 703   Double_t dataE[9]={-999.,-999.,-999.,-999.,-1.,0,0,-999.,-999.};
 704
 705   Int_t paircode1 = 0, paircode2 = 0;
 706   Int_t htracklabel1 = 0, htracklabel2 = 0;
 707
 708   if (HasMCData()){
 709     AliHFEpairs *pair1=0x0;
 710     AliHFEpairs *pair2=0x0;
 711     pair1 = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ipair));
 712     pair2 = (AliHFEpairs*) (HFEpairs()->UncheckedAt(jpair));
 713
 714     htracklabel1 = pair1->GetTrkLabel();
 715     htracklabel2 = pair2->GetTrkLabel();
 716
 717     if (pair1->GetPairCode()==2 || pair1->GetPairCode()==3) paircode1=1;
 718     else paircode1=0;
 719     if (pair2->GetPairCode()==2 || pair2->GetPairCode()==3) paircode2=1;
 720     else paircode2=0;
 721   }
 722
 723   // fill secondary vertex container
 724   AliHFEsecVtxs hfesecvtx;
 725   hfesecvtx.SetTrkLabel1(htracklabel1);
 726   hfesecvtx.SetTrkLabel2(htracklabel2);
 727   if(HasMCData()) hfesecvtx.SetMCCode(GetElectronSource(TMath::Abs(track->GetLabel())));
 728   hfesecvtx.SetKFChi2(fKFchi2);
 729   hfesecvtx.SetInvmass(fInvmass);
 730   hfesecvtx.SetSignedLxy(fSignedLxy);
 731   hfesecvtx.SetSignedLxy2(fSignedLxy2);
 732   hfesecvtx.SetKFIP(fKFip);
 733   hfesecvtx.SetKFIP2(fKFip2);
 734   AddHFEsecvtxToArray(&hfesecvtx);
 735   fNoOfHFEsecvtxs++;
 736
 737   // fill debugging THnSparse
 738   dataE[0]=fInvmass;
 739   dataE[1]=fKFchi2;
 740   dataE[2]=fSignedLxy;
 741   dataE[3]=fKFip;
 742   if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
 743   dataE[5]=3;
 744   if(paircode1 & paircode2) dataE[6]=1;
 745   else if(!paircode1 & !paircode2) dataE[6]=0;
 746   else dataE[6]=3;
 747   dataE[7]=fSignedLxy2;
 748   dataE[8]=track->Pt();
 749   fSecvtxQA->Fill(dataE);
 750
 751 }
 752
 753 //_______________________________________________________________________________________________
 754 void AliHFEsecVtx::Fill2TrkSECVTX(AliVTrack* track, AliHFEpairs *pair)
 755 {
 756   //
 757   // fill 2 tracks' secondary vertex properties
 758   //
 759
 760   Double_t dataE[9]={-999.,-999.,-999.,-999.,-1.,0,0,-999.,-999.};
 761
 762   Int_t paircode;
 763   if (pair->GetPairCode()==2 || pair->GetPairCode()==3) paircode=1;
 764   else paircode=0;
 765
 766   // fill secondary vertex container
 767   AliHFEsecVtxs hfesecvtx;
 768   hfesecvtx.SetTrkLabel1(pair->GetTrkLabel());
 769   hfesecvtx.SetTrkLabel2(-999);
 770   if(HasMCData()) hfesecvtx.SetMCCode(GetElectronSource(TMath::Abs(track->GetLabel())));
 771   hfesecvtx.SetInvmass(pair->GetInvmass());
 772   hfesecvtx.SetKFChi2(pair->GetKFChi2());
 773   hfesecvtx.SetSignedLxy(pair->GetSignedLxy());
 774   hfesecvtx.SetSignedLxy2(pair->GetSignedLxy2());
 775   hfesecvtx.SetKFIP(pair->GetKFIP());
 776   hfesecvtx.SetKFIP2(pair->GetKFIP2());
 777   AddHFEsecvtxToArray(&hfesecvtx);
 778   fNoOfHFEsecvtxs++;
 779
 780   // fill debugging THnSparse
 781   dataE[0]=pair->GetInvmass();
 782   dataE[1]=pair->GetKFChi2();
 783   dataE[2]=pair->GetSignedLxy();
 784   dataE[3]=pair->GetKFIP();
 785   if (HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
 786   dataE[5]=2;              //# of associated tracks
 787   dataE[6]=paircode;
 788   dataE[7]=pair->GetSignedLxy2();
 789   dataE[8]=track->Pt();
 790   fSecvtxQA->Fill(dataE);
 791
 792 }
 793
 794 //_______________________________________________________________________________________________
 795 void AliHFEsecVtx::CalcSECVTXProperty(AliVTrack* track1, AliVTrack* track2, AliVTrack* track3)
 796 {
 797   //
 798   // calculate secondary vertex properties
 799   //
 800
 801   // get KF particle input pid
 802   Int_t pdg1 = GetPDG(track1);
 803   Int_t pdg2 = GetPDG(track2);
 804   Int_t pdg3 = GetPDG(track3);
 805
 806   if(pdg1==-1 || pdg2==-1 || pdg3==-1) {
 807     //printf("out if considered pid range \n");
 808     return;
 809   }
 810
 811   // create KF particle of pair
 812   if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
 813   else AliKFParticle::SetField(fESD1->GetMagneticField());
 814   AliKFParticle kfTrack[3];
 815   kfTrack[0] = AliKFParticle(*track1, pdg1);
 816   kfTrack[1] = AliKFParticle(*track2, pdg2);
 817   kfTrack[2] = AliKFParticle(*track3, pdg3);
 818
 819   AliKFParticle kfSecondary(kfTrack[0],kfTrack[1],kfTrack[2]);
 820   //AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
 821
 822   //secondary vertex point from kf particle
 823   Double_t kfx = kfSecondary.GetX();
 824   Double_t kfy = kfSecondary.GetY();
 825   //Double_t kfz = kfSecondary.GetZ();
 826
 827   //momentum at the decay point from kf particle
 828   Double_t kfpx = kfSecondary.GetPx();
 829   Double_t kfpy = kfSecondary.GetPy();
 830   //Double_t kfpz = kfSecondary.GetPz();
 831
 832   Double_t dx = kfx-fPVx;
 833   Double_t dy = kfy-fPVy;
 834
 835   // discriminating variables ----------------------------------------------------------
 836
 837   if(kfSecondary.GetNDF()>0) fKFchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 838
 839   // invariant mass of the KF particle
 840   kfSecondary.GetMass(fInvmass,fInvmassSigma);
 841
 842   // DCA from primary to e-h KF particle (impact parameter of KF particle)
 843   Double_t vtx[2]={fPVx, fPVy};
 844   fKFip = kfSecondary.GetDistanceFromVertexXY(vtx);
 845
 846   if((dx*kfpx+dy*kfpy)>0) fSignedLxy= TMath::Sqrt(dx*dx+dy*dy);
 847   if((dx*kfpx+dy*kfpy)<0) fSignedLxy= -1*TMath::Sqrt(dx*dx+dy*dy);
 848   //[the other way to think about] - projection of kf vertex vector to the kf momentum direction
 849   //Double_t psqr = kfpx*kfpx+kfpy*kfpy;
 850   //if(psqr>0) fSignedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
 851
 852
 853   //recalculating primary vertex after removing secvtx tracks --------------------------
 854   Int_t trkid[3];
 855   trkid[0] = track1->GetID();
 856   trkid[1] = track2->GetID();
 857   trkid[2] = track3->GetID();
 858
 859   RecalcPrimvtx(3, trkid, kfTrack);
 860   Double_t dx2 = kfx-fPVx2;
 861   Double_t dy2 = kfy-fPVy2;
 862
 863   // IP of sec particle recalculated based on recalculated primary vertex
 864   Double_t vtx2[2]={fPVx2, fPVy2};
 865   fKFip2 = kfSecondary.GetDistanceFromVertexXY(vtx2);
 866   // signed Lxy recalculated based on recalculated primary vertex
 867   if((dx2*kfpx+dy2*kfpy)>0) fSignedLxy2= TMath::Sqrt(dx2*dx2+dy2*dy2);
 868   if((dx2*kfpx+dy2*kfpy)<0) fSignedLxy2= -1*TMath::Sqrt(dx2*dx2+dy2*dy2);
 869   //------------------------------------------------------------------------------------
 870
 871 }
 872
 873 //_______________________________________________________________________________________________
 874 void AliHFEsecVtx::CalcSECVTXProperty(AliVTrack* track1, AliVTrack* track2, AliVTrack* track3, AliVTrack* track4)
 875 {
 876   //
 877   // calculate secondary vertex properties
 878   //
 879
 880   // get KF particle input pid
 881   Int_t pdg1 = GetPDG(track1);
 882   Int_t pdg2 = GetPDG(track2);
 883   Int_t pdg3 = GetPDG(track3);
 884   Int_t pdg4 = GetPDG(track4);
 885
 886   if(pdg1==-1 || pdg2==-1 || pdg3==-1 || pdg4==-1) {
 887     //printf("out if considered pid range \n");
 888     return;
 889   }
 890
 891   // create KF particle of pair
 892   if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
 893   else AliKFParticle::SetField(fESD1->GetMagneticField());
 894
 895   AliKFParticle kfTrack[4];
 896   kfTrack[0] = AliKFParticle(*track1, pdg1);
 897   kfTrack[1] = AliKFParticle(*track2, pdg2);
 898   kfTrack[2] = AliKFParticle(*track3, pdg3);
 899   kfTrack[3] = AliKFParticle(*track4, pdg4);
 900
 901   AliKFParticle kfSecondary(kfTrack[0],kfTrack[1],kfTrack[2],kfTrack[3]);
 902
 903   //secondary vertex point from kf particle
 904   Double_t kfx = kfSecondary.GetX();
 905   Double_t kfy = kfSecondary.GetY();
 906   //Double_t kfz = kfSecondary.GetZ();
 907
 908   //momentum at the decay point from kf particle
 909   Double_t kfpx = kfSecondary.GetPx();
 910   Double_t kfpy = kfSecondary.GetPy();
 911   //Double_t kfpz = kfSecondary.GetPz();
 912
 913   Double_t dx = kfx-fPVx;
 914   Double_t dy = kfy-fPVy;
 915
 916   // discriminating variables ----------------------------------------------------------
 917
 918   if(kfSecondary.GetNDF()>0) fKFchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 919
 920   // invariant mass of the KF particle
 921   kfSecondary.GetMass(fInvmass,fInvmassSigma);
 922
 923   // DCA from primary to e-h KF particle (impact parameter of KF particle)
 924   Double_t vtx[2]={fPVx, fPVy};
 925   fKFip = kfSecondary.GetDistanceFromVertexXY(vtx);
 926
 927   if((dx*kfpx+dy*kfpy)>0) fSignedLxy= TMath::Sqrt(dx*dx+dy*dy);
 928   if((dx*kfpx+dy*kfpy)<0) fSignedLxy= -1*TMath::Sqrt(dx*dx+dy*dy);
 929   //[the other way to think about] - projection of kf vertex vector to the kf momentum direction
 930   //Double_t psqr = kfpx*kfpx+kfpy*kfpy;
 931   //if(psqr>0) fSignedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
 932
 933   //recalculating primary vertex after removing secvtx tracks --------------------------
 934   Int_t trkid[4];
 935   trkid[0] = track1->GetID();
 936   trkid[1] = track2->GetID();
 937   trkid[2] = track3->GetID();
 938   trkid[3] = track4->GetID();
 939
 940   RecalcPrimvtx(4, trkid, kfTrack);
 941   Double_t dx2 = kfx-fPVx2;
 942   Double_t dy2 = kfy-fPVy2;
 943
 944   // IP of sec particle recalculated based on recalculated primary vertex
 945   Double_t vtx2[2]={fPVx2, fPVy2};
 946   fKFip2 = kfSecondary.GetDistanceFromVertexXY(vtx2);
 947   // signed Lxy recalculated based on recalculated primary vertex
 948   if((dx2*kfpx+dy2*kfpy)>0) fSignedLxy2= TMath::Sqrt(dx2*dx2+dy2*dy2);
 949   if((dx2*kfpx+dy2*kfpy)<0) fSignedLxy2= -1*TMath::Sqrt(dx2*dx2+dy2*dy2);
 950   //------------------------------------------------------------------------------------
 951
 952 }
 953
 954 //_______________________________________________________________________________________________
 955 void AliHFEsecVtx::RecalcPrimvtx(Int_t nkftrk, const Int_t * const trkid, const AliKFParticle * const kftrk){
 956   //
 957   // reccalculate primary vertex after removing considering track in the calculation
 958   //
 959
 960   const AliESDVertex *primvtx = fESD1->GetPrimaryVertex();
 961
 962   AliKFVertex kfESDprimary;
 963   Int_t n = primvtx->GetNIndices();
 964   fNsectrk2prim = 0;
 965   fPVx2 = -999.;
 966   fPVy2 = -999.;
 967
 968   if (n>0 && primvtx->GetStatus()){
 969     kfESDprimary = AliKFVertex(*primvtx);
 970     UShort_t *priIndex = primvtx->GetIndices();
 971     for(Int_t j=0; j<nkftrk; j++){
 972       for (Int_t i=0;i<n;i++){
 973         Int_t idx = Int_t(priIndex[i]);
 974         if (idx == trkid[j]){
 975           kfESDprimary -= kftrk[j];
 976           fNsectrk2prim++;
 977         }
 978       }
 979     }
 980   }
 981
 982   fPVx2 = kfESDprimary.GetX();
 983   fPVy2 = kfESDprimary.GetY();
 984
 985 }
 986
 987 //_______________________________________________________________________________________________
 988 Int_t AliHFEsecVtx::GetMCPID(const AliESDtrack *track)
 989 {
 990   //
 991   // return mc pid
 992   //
 993
 994   AliMCParticle *mctrack = NULL;
 995   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) return 0;
 996   TParticle *mcpart = mctrack->Particle();
 997
 998   if ( !mcpart ) return 0;
 999   Int_t pdgCode = mcpart->GetPdgCode();
1000
1001   return pdgCode;
1002 }
1003
1004 //_______________________________________________________________________________________________
1005 Int_t AliHFEsecVtx::GetPairOriginESD(AliESDtrack* trk1, AliESDtrack* trk2)
1006 {
1007   //
1008   // return pdg code of the origin(source) of the pair
1009   //
1010   //
1011   // ---*---*---*-----ancester A----- track1
1012   //                        |____*______
1013   //                             |______ track2
1014   // => if they originated from same ancester,
1015   //    then return "the absolute value of pdg code of ancester A"
1016   //
1017   // ---*---*---B hadron-----ancester A----- track1
1018   //                               |____*______
1019   //                                    |______ track2
1020   // => if they originated from same ancester, and this ancester originally comes from B hadrons
1021   //    then return -1*"the absolute value of pdg code of ancester A"
1022   //
1023   // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
1024   //
1025
1026   if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
1027
1028   AliMCParticle *mctrack = NULL;
1029   AliMCParticle *mctrack1 = NULL;
1030   AliMCParticle *mctrack2 = NULL;
1031   if(!(mctrack1 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(trk1->GetLabel()))))) return 0;
1032   if(!(mctrack2 = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(trk2->GetLabel()))))) return 0;
1033   TParticle *part1 = mctrack1->Particle();
1034   TParticle *part2 = mctrack2->Particle();
1035
1036   TParticle* part2cp = part2;
1037   if (!(part1) || !(part2)) return 0;
1038
1039   Int_t srcpdg = 0;
1040
1041   //if the two tracks' mother's label is same, get the mother info
1042   //in case of charm, check if it originated from beauty
1043   for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
1044      Int_t label1 = part1->GetFirstMother();
1045      if (label1 < 0) return 0;
1046
1047      for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
1048         Int_t label2 = part2->GetFirstMother();
1049         if (label2 < 0) break;
1050
1051         if (label1 == label2){ //check if two tracks are originated from same mother
1052           if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(label2))))) return 0;
1053           TParticle* commonmom = mctrack->Particle();
1054
1055           srcpdg = abs(commonmom->GetPdgCode());
1056
1057           //check ancester to see if it is originally from beauty
1058           for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
1059              Int_t ancesterlabel = commonmom->GetFirstMother();
1060              if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
1061
1062              if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(ancesterlabel))))) return 0;
1063              TParticle* commonancester = mctrack->Particle();
1064
1065              Int_t ancesterpdg = abs(commonancester->GetPdgCode());
1066
1067              for (Int_t l=0; l<fNparents; l++){
1068                 if (abs(ancesterpdg)==fParentSelect[1][l]){
1069                   srcpdg = -1*srcpdg; //multiply -1 for hadron from bottom
1070                   return srcpdg;
1071                 }
1072              }
1073              commonmom = commonancester;
1074           }
1075         }
1076         if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(label2))))) return 0;
1077         part2 = mctrack->Particle(); //if their mother is different, go to earlier generation of 2nd particle
1078
1079         if (!(part2)) break;
1080      }
1081      if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(label1))))) return 0;
1082      part1 = mctrack->Particle(); //if their mother is different, go to earlier generation of 1st particle
1083      part2 = part2cp;
1084      if (!(part1)) return 0;
1085   }
1086
1087   return srcpdg;
1088 }
1089
1090 //_______________________________________________________________________________________________
1091 Int_t AliHFEsecVtx::GetPairOriginAOD(AliAODTrack* trk1, AliAODTrack* trk2)
1092 {
1093
1094   //
1095   // return pdg code of the origin(source) of the pair
1096   //
1097   //
1098   // ---*---*---*-----ancester A----- track1
1099   //                        |____*______
1100   //                             |______ track2
1101   // => if they originated from same ancester,
1102   //    then return "the absolute value of pdg code of ancester A"
1103   //
1104   // ---*---*---B hadron-----ancester A----- track1
1105   //                               |____*______
1106   //                                    |______ track2
1107   // => if they originated from same ancester, and this ancester originally comes from B hadrons
1108   //    then return -1*"the absolute value of pdg code of ancester A"
1109   //
1110   // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
1111   //
1112
1113   if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
1114   AliAODMCParticle *part1 = (AliAODMCParticle*)fMCArray->At(trk1->GetLabel());
1115   AliAODMCParticle *part2 = (AliAODMCParticle*)fMCArray->At(trk2->GetLabel());
1116   AliAODMCParticle *part2cp = part2;
1117   if (!(part1) || !(part2)) return 0;
1118
1119   Int_t srcpdg = 0;
1120
1121   //if the two tracks' mother's label is same, get the mother info
1122   //in case of charm, check if it originated from beauty
1123   for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
1124      Int_t label1 = part1->GetMother();
1125      if (label1 < 0) return 0;
1126
1127      for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
1128         Int_t label2 = part2->GetMother();
1129         if (label2 < 0) return 0;
1130
1131         if (label1 == label2){ //check if two tracks are originated from same mother
1132           AliAODMCParticle *commonmom = (AliAODMCParticle*)fMCArray->At(label1);
1133           srcpdg = abs(commonmom->GetPdgCode());
1134
1135           //check ancester to see if it is originally from beauty
1136           for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
1137              Int_t ancesterlabel = commonmom->GetMother();
1138              if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
1139
1140              AliAODMCParticle *commonancester = (AliAODMCParticle*)fMCArray->At(ancesterlabel);
1141              Int_t ancesterpdg = abs(commonancester->GetPdgCode());
1142
1143              for (Int_t l=0; l<fNparents; l++){
1144                 if (abs(ancesterpdg)==fParentSelect[1][l]){
1145                   srcpdg = -1*srcpdg; //multiply -1 for charm from bottom
1146                   return srcpdg;
1147                 }
1148              }
1149              commonmom = commonancester;
1150           }
1151         }
1152         part2 = (AliAODMCParticle*)fMCArray->At(label2); //if their mother is different, go to earlier generation of 2nd particle
1153         if (!(part2)) break;
1154      }
1155      part1 = (AliAODMCParticle*)fMCArray->At(label1); //if their mother is different, go to earlier generation of 2nd particle
1156      part2 = part2cp;
1157      if (!(part1)) return 0;
1158   }
1159
1160   return srcpdg;
1161 }
1162
1163 //_______________________________________________________________________________________________
1164 Int_t AliHFEsecVtx::GetPairCode(const AliVTrack* const trk1, const AliVTrack* const trk2)
1165 {
1166   //
1167   // return pair code which is predefinded as:
1168   //  kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0,
1169   //  kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
1170   //
1171
1172   Int_t srcpdg = -1;
1173   Int_t srccode = kElse;
1174
1175   if(IsAODanalysis()) srcpdg = GetPairOriginAOD((AliAODTrack*)trk1,(AliAODTrack*)trk2);
1176   else srcpdg = GetPairOriginESD((AliESDtrack*)trk1,(AliESDtrack*)trk2);
1177
1178   if (srcpdg < 0) srccode = kBeautyElse;
1179   for (Int_t i=0; i<fNparents; i++){
1180     if (abs(srcpdg)==fParentSelect[0][i]) {
1181       if (srcpdg>0) srccode = kDirectCharm;
1182       if (srcpdg<0) srccode = kBeautyCharm;
1183     }
1184     if (abs(srcpdg)==fParentSelect[1][i]) {
1185       if (srcpdg>0) srccode = kDirectBeauty;
1186       if (srcpdg<0)  return kElse;
1187     }
1188   }
1189   if (srcpdg == 22) srccode = kGamma;
1190   if (srcpdg == -22) srccode = kBeautyGamma;
1191   if (srcpdg == 111) srccode = kPi0;
1192   if (srcpdg == -111) srccode = kBeautyPi0;
1193
1194   return srccode;
1195 }
1196
1197 //_______________________________________________________________________________________________
1198 Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack)
1199 {
1200   //
1201   // return decay electron's origin
1202   //
1203
1204   if (iTrack < 0) {
1205     AliDebug(1, "Stack label is negative, return\n");
1206     return -1;
1207   }
1208
1209   AliMCParticle *mctrack = NULL;
1210   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iTrack))))) return -1;
1211   TParticle *mcpart = mctrack->Particle();
1212
1213   if(!mcpart){
1214     AliDebug(1, "no mc particle, return\n");
1215     return -1;
1216   }
1217
1218   if ( abs(mcpart->GetPdgCode()) != 11 ) return kMisID;
1219
1220 //  if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !
1221
1222   Int_t iLabel = mcpart->GetFirstMother();
1223   if (iLabel<0){
1224     AliDebug(1, "Stack label is negative, return\n");
1225     return -1;
1226   }
1227
1228   Int_t origin = -1;
1229   Bool_t isFinalOpenCharm = kFALSE;
1230
1231   if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1;
1232   TParticle *partMother = mctrack->Particle();
1233
1234   Int_t maPdgcode = partMother->GetPdgCode();
1235
1236   // if the mother is charmed hadron
1237   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1238
1239     for (Int_t i=0; i<fNparents; i++){
1240       if (abs(maPdgcode)==fParentSelect[0][i]){
1241         isFinalOpenCharm = kTRUE;
1242       }
1243     }
1244     if (!isFinalOpenCharm) return -1;
1245
1246     // iterate until you find B hadron as a mother or become top ancester
1247     for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1248
1249       Int_t jLabel = partMother->GetFirstMother();
1250       if (jLabel == -1){
1251         origin = kDirectCharm;
1252         return origin;
1253       }
1254       if (jLabel < 0){ // safety protection
1255         AliDebug(1, "Stack label is negative, return\n");
1256         return -1;
1257       }
1258
1259       // if there is an ancester
1260       if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1261       TParticle *grandMa = mctrack->Particle();
1262
1263       Int_t grandMaPDG = grandMa->GetPdgCode();
1264
1265       for (Int_t j=0; j<fNparents; j++){
1266         if (abs(grandMaPDG)==fParentSelect[1][j]){
1267           origin = kBeautyCharm;
1268           return origin;
1269         }
1270       }
1271
1272       partMother = grandMa;
1273     } // end of iteration
1274   } // end of if
1275   else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1276     for (Int_t i=0; i<fNparents; i++){
1277       if (abs(maPdgcode)==fParentSelect[1][i]){
1278         origin = kDirectBeauty;
1279         return origin;
1280       }
1281     }
1282   } // end of if
1283
1284   //============ gamma ================
1285   else if ( abs(maPdgcode) == 22 ) {
1286     origin = kGamma;
1287
1288     // iterate until you find B hadron as a mother or become top ancester
1289     for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1290
1291       Int_t jLabel = partMother->GetFirstMother();
1292       if (jLabel == -1){
1293         origin = kGamma;
1294         return origin;
1295       }
1296       if (jLabel < 0){ // safety protection
1297         AliDebug(1, "Stack label is negative, return\n");
1298         return -1;
1299       }
1300
1301       // if there is an ancester
1302       if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1303       TParticle *grandMa = mctrack->Particle();
1304
1305       Int_t grandMaPDG = grandMa->GetPdgCode();
1306
1307       for (Int_t j=0; j<fNparents; j++){
1308         if (abs(grandMaPDG)==fParentSelect[1][j]){
1309           origin = kBeautyGamma;
1310           return origin;
1311         }
1312       }
1313
1314       partMother = grandMa;
1315     } // end of iteration
1316
1317     return origin;
1318   } // end of if
1319
1320   //============ pi0 ================
1321   else if ( abs(maPdgcode) == 111 ) {
1322     origin = kPi0;
1323
1324     // iterate until you find B hadron as a mother or become top ancester
1325     for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1326
1327       Int_t jLabel = partMother->GetFirstMother();
1328       if (jLabel == -1){
1329         origin = kPi0;
1330         return origin;
1331       }
1332       if (jLabel < 0){ // safety protection
1333         AliDebug(1, "Stack label is negative, return\n");
1334         return -1;
1335       }
1336
1337       // if there is an ancester
1338       if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1339       TParticle *grandMa = mctrack->Particle();
1340
1341       Int_t grandMaPDG = grandMa->GetPdgCode();
1342
1343       for (Int_t j=0; j<fNparents; j++){
1344         if (abs(grandMaPDG)==fParentSelect[1][j]){
1345           origin = kBeautyPi0;
1346           return origin;
1347         }
1348       }
1349
1350       partMother = grandMa;
1351     } // end of iteration
1352
1353     return origin;
1354   } // end of if
1355
1356   else {
1357     origin = kElse;
1358     // iterate until you find B hadron as a mother or become top ancester
1359     for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1360
1361       Int_t jLabel = partMother->GetFirstMother();
1362       if (jLabel == -1){
1363         origin = kElse;
1364         return origin;
1365       }
1366       if (jLabel < 0){ // safety protection
1367         AliDebug(1, "Stack label is negative, return\n");
1368         return -1;
1369       }
1370
1371       // if there is an ancester
1372       if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1;
1373       TParticle *grandMa = mctrack->Particle();
1374
1375       Int_t grandMaPDG = grandMa->GetPdgCode();
1376
1377       for (Int_t j=0; j<fNparents; j++){
1378         if (abs(grandMaPDG)==fParentSelect[1][j]){
1379           origin = kBeautyElse;
1380           return origin;
1381         }
1382       }
1383
1384       partMother = grandMa;
1385     } // end of iteration
1386   }
1387
1388   return origin;
1389 }
1390
1391 //_______________________________________________________________________________________________
1392 Int_t AliHFEsecVtx::GetPDG(AliVTrack *track)
1393 {
1394   //
1395   // get KF particle input pdg for mass hypothesis
1396   //
1397
1398   Int_t pdgCode=-1;
1399
1400   if (fUseMCPID && HasMCData()){
1401     pdgCode = GetMCPDG(track);
1402   }
1403   else if(fESD1){
1404     Int_t pid=0;
1405     Double_t prob;
1406     GetESDPID((AliESDtrack*)track, pid, prob);
1407     switch(pid){
1408       case 0:  pdgCode = 11; break;
1409       case 1:  pdgCode = 13; break;
1410       case 2:  pdgCode = 211; break;
1411       case 3:  pdgCode = 321; break;
1412       case 4:  pdgCode = 2212; break;
1413       default: pdgCode = -1;
1414     }
1415   }
1416   else if(fAOD1){
1417     Int_t pid = ((AliAODTrack*)track)->GetMostProbablePID();
1418     switch(pid){
1419       case 0:  pdgCode = 11; break;
1420       case 1:  pdgCode = 13; break;
1421       case 2:  pdgCode = 211; break;
1422       case 3:  pdgCode = 321; break;
1423       case 4:  pdgCode = 2212; break;
1424       default: pdgCode = -1;
1425     }
1426   }
1427
1428   return pdgCode;
1429 }
1430
1431 //_______________________________________________________________________________________________
1432 Int_t AliHFEsecVtx::GetMCPDG(const AliVTrack *track)
1433 {
1434   //
1435   // return mc pdg code
1436   //
1437
1438   Int_t label = TMath::Abs(track->GetLabel());
1439   Int_t pdgCode;
1440   AliMCParticle *mctrack = NULL;
1441
1442   if (IsAODanalysis()) {
1443     AliAODMCParticle *mcpart = (AliAODMCParticle*)fMCArray->At(label);
1444     if ( !mcpart ) return 0;
1445       pdgCode = mcpart->GetPdgCode();
1446   }
1447   else {
1448     if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(label))))) return 0;
1449     TParticle *mcpart = mctrack->Particle();
1450
1451     if ( !mcpart ) return 0;
1452       pdgCode = mcpart->GetPdgCode();
1453   }
1454
1455     return pdgCode;
1456 }
1457
1458 //______________________________________________________________________________
1459 void AliHFEsecVtx::GetESDPID(AliESDtrack *track, Int_t &recpid, Double_t &recprob)
1460 {
1461   //
1462   // calculate likehood for esd pid
1463   //
1464
1465   recpid = -1;
1466   recprob = -1;
1467
1468   Int_t ipid=-1;
1469   Double_t max=0.;
1470
1471   Double_t probability[5];
1472
1473   // get probability of the diffenrent particle types
1474   track->GetESDpid(probability);
1475
1476   // find most probable particle in ESD pid
1477   // 0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
1478   ipid = TMath::LocMax(5,probability);
1479   max = TMath::MaxElement(5,probability);
1480
1481   recpid = ipid;
1482   recprob = max;
1483 }
1484
1485 //_____________________________________________________________________________
1486 void AliHFEsecVtx::AddHFEpairToArray(const AliHFEpairs* const pair)
1487 {
1488   //
1489   // Add a HFE pair to the array
1490   //
1491
1492   Int_t n = HFEpairs()->GetEntriesFast();
1493   if(n!=fNoOfHFEpairs)AliError(Form("fNoOfHFEpairs != HFEpairs()->GetEntriesFast %i != %i \n", fNoOfHFEpairs, n));
1494   new((*HFEpairs())[n]) AliHFEpairs(*pair);
1495 }
1496
1497 //_____________________________________________________________________________
1498 TClonesArray *AliHFEsecVtx::HFEpairs()
1499 {
1500   //
1501   // Returns the list of HFE pairs
1502   //
1503
1504   if (!fHFEpairs) {
1505       fHFEpairs = new TClonesArray("AliHFEpairs", 200);
1506   }
1507   return fHFEpairs;
1508 }
1509
1510 //_____________________________________________________________________________
1511 void AliHFEsecVtx::DeleteHFEpairs()
1512 {
1513   //
1514   // Delete the list of HFE pairs
1515   //
1516
1517   if (fHFEpairs) {
1518     fHFEpairs->Delete();
1519     //delete fHFEpairs;
1520   }
1521 }
1522
1523 //_____________________________________________________________________________
1524 void AliHFEsecVtx::InitHFEpairs()
1525 {
1526   //
1527   // Initialization should be done before make all possible pairs for a given electron candidate
1528   //
1529
1530   fNoOfHFEpairs = 0;
1531 }
1532
1533 //_____________________________________________________________________________
1534 void AliHFEsecVtx::AddHFEsecvtxToArray(const AliHFEsecVtxs* const secvtx)
1535 {
1536   //
1537   // Add a HFE secondary vertex to the array
1538   //
1539
1540   Int_t n = HFEsecvtxs()->GetEntriesFast();
1541   if(n!=fNoOfHFEsecvtxs)AliError(Form("fNoOfHFEsecvtxs != HFEsecvtxs()->GetEntriesFast %i != %i \n", fNoOfHFEsecvtxs, n));
1542   new((*HFEsecvtxs())[n]) AliHFEsecVtxs(*secvtx);
1543 }
1544
1545 //_____________________________________________________________________________
1546 TClonesArray *AliHFEsecVtx::HFEsecvtxs()
1547 {
1548   //
1549   // Returns the list of HFE secvtx
1550   //
1551
1552   if (!fHFEsecvtxs) {
1553       fHFEsecvtxs = new TClonesArray("AliHFEsecVtxs", 200);
1554   }
1555   return fHFEsecvtxs;
1556 }
1557
1558 //_____________________________________________________________________________
1559 void AliHFEsecVtx::DeleteHFEsecvtxs()
1560 {
1561   //
1562   // Delete the list of HFE pairs
1563   //
1564
1565   if (fHFEsecvtxs) {
1566     fHFEsecvtxs->Delete();
1567     //delete fHFEsecvtx;
1568   }
1569 }
1570
1571 //_____________________________________________________________________________
1572 void AliHFEsecVtx::InitHFEsecvtxs()
1573 {
1574   //
1575   // Initialization should be done
1576   //
1577
1578   fNoOfHFEsecvtxs = 0;
1579 }
1580
1581 //____________________________________________________________
1582 void AliHFEsecVtx::MakeContainer(){
1583
1584   //
1585   // make container
1586   //
1587
1588   const Int_t nDimPair=6;
1589   Int_t nBinPair[nDimPair] = {200, 500, 314, 2000, 2000, 13};
1590   //Int_t nBinPair[nDimPair] = {200, 500, 314, 2000, 2000, 13, 60, 60, 2000, 2000};
1591   const Double_t kInvmassmin = 0., kInvmassmax = 20.;
1592   const Double_t kKFChi2min = 0, kKFChi2max= 50;
1593   const Double_t kOpenanglemin = 0, kOpenanglemax = 3.14;
1594   const Double_t kSignedLxymin = -10, kSignedLxymax= 10;
1595   const Double_t kKFIPmin = -10, kKFIPmax= 10;
1596   const Double_t kPairCodemin = -1, kPairCodemax= 12;
1597   //const Double_t kPtmin = 0, kPtmax= 30;
1598   //const Double_t kDCAsigmin = -5, kDCAsigmax= 5;
1599
1600   Double_t* binEdgesPair[nDimPair];
1601   for(Int_t ivar = 0; ivar < nDimPair; ivar++)
1602     binEdgesPair[ivar] = new Double_t[nBinPair[ivar] + 1];
1603
1604   for(Int_t i=0; i<=nBinPair[0]; i++) binEdgesPair[0][i]=(Double_t)kInvmassmin + (kInvmassmax - kInvmassmin)/nBinPair[0]*(Double_t)i;
1605   for(Int_t i=0; i<=nBinPair[1]; i++) binEdgesPair[1][i]=(Double_t)kKFChi2min + (kKFChi2max - kKFChi2min)/nBinPair[1]*(Double_t)i;
1606   for(Int_t i=0; i<=nBinPair[2]; i++) binEdgesPair[2][i]=(Double_t)kOpenanglemin + (kOpenanglemax - kOpenanglemin)/nBinPair[2]*(Double_t)i;
1607   for(Int_t i=0; i<=nBinPair[3]; i++) binEdgesPair[3][i]=(Double_t)kSignedLxymin + (kSignedLxymax - kSignedLxymin)/nBinPair[3]*(Double_t)i;
1608   for(Int_t i=0; i<=nBinPair[4]; i++) binEdgesPair[4][i]=(Double_t)kKFIPmin + (kKFIPmax - kKFIPmin)/nBinPair[4]*(Double_t)i;
1609   for(Int_t i=0; i<=nBinPair[5]; i++) binEdgesPair[5][i]=(Double_t)kPairCodemin + (kPairCodemax - kPairCodemin)/nBinPair[5]*(Double_t)i;
1610   //for(Int_t i=0; i<=nBinPair[6]; i++) binEdgesPair[6][i]=(Double_t)kPtmin + (kPtmax - kPtmin)/nBinPair[6]*(Double_t)i;
1611   //for(Int_t i=0; i<=nBinPair[7]; i++) binEdgesPair[7][i]=binEdgesPair[6][i];
1612   //for(Int_t i=0; i<=nBinPair[6]; i++) binEdgesPair[6][i]=(Double_t)kDCAsigmin + (kDCAsigmax - kDCAsigmin)/nBinPair[6]*(Double_t)i;
1613   //for(Int_t i=0; i<=nBinPair[7]; i++) binEdgesPair[7][i]=binEdgesPair[6][i];
1614
1615   fPairQA = new THnSparseF("pairQA", "QA for Pair; invmass[GeV/c^2]; KF chi2; opening angle; signed Lxy; KF ip; pair code; dca1; dca2", nDimPair, nBinPair);
1616   //fPairQA = new THnSparseF("pairQA", "QA for Pair; invmass[GeV/c^2]; KF chi2; opening angle; signed Lxy; KF ip; pair code; pt1; pt2; dca1; dca2", nDimPair, nBinPair);
1617   for(Int_t idim = 0; idim < nDimPair; idim++){
1618     fPairQA->SetBinEdges(idim, binEdgesPair[idim]);
1619   }
1620
1621   fSecVtxList->AddAt(fPairQA,0);
1622
1623   const Int_t nDimSecvtx=9;
1624   Double_t* binEdgesSecvtx[nDimSecvtx];
1625   Int_t nBinSecvtx[nDimSecvtx] = {200, 500, 2000, 2000, 13, 10, 4, 2000, 500};
1626   const Double_t kNtrksmin = 0, kNtrksmax= 10;
1627   const Double_t kTrueBmin = 0, kTrueBmax= 4;
1628   const Double_t kPtmin = 0, kPtmax= 50;
1629   for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
1630     binEdgesSecvtx[ivar] = new Double_t[nBinSecvtx[ivar] + 1];
1631
1632   for(Int_t i=0; i<=nBinSecvtx[0]; i++) binEdgesSecvtx[0][i]=binEdgesPair[0][i];
1633   for(Int_t i=0; i<=nBinSecvtx[1]; i++) binEdgesSecvtx[1][i]=binEdgesPair[1][i];
1634   for(Int_t i=0; i<=nBinSecvtx[2]; i++) binEdgesSecvtx[2][i]=binEdgesPair[3][i];
1635   for(Int_t i=0; i<=nBinSecvtx[3]; i++) binEdgesSecvtx[3][i]=binEdgesPair[4][i];
1636   for(Int_t i=0; i<=nBinSecvtx[4]; i++) binEdgesSecvtx[4][i]=binEdgesPair[5][i];
1637   for(Int_t i=0; i<=nBinSecvtx[5]; i++) binEdgesSecvtx[5][i]=(Double_t)kNtrksmin + (kNtrksmax - kNtrksmin)/nBinSecvtx[5]*(Double_t)i;
1638   for(Int_t i=0; i<=nBinSecvtx[6]; i++) binEdgesSecvtx[6][i]=(Double_t)kTrueBmin + (kTrueBmax - kTrueBmin)/nBinSecvtx[6]*(Double_t)i;
1639   for(Int_t i=0; i<=nBinSecvtx[7]; i++) binEdgesSecvtx[7][i]=binEdgesPair[3][i];
1640   for(Int_t i=0; i<=nBinSecvtx[8]; i++) binEdgesSecvtx[8][i]=(Double_t)kPtmin + (kPtmax - kPtmin)/nBinSecvtx[8]*(Double_t)i;
1641
1642   fSecvtxQA = new THnSparseF("secvtxQA", "QA for Secvtx; invmass[GeV/c^2]; KF chi2; signed Lxy; KF ip; pair code; n tracks ", nDimSecvtx, nBinSecvtx);
1643   for(Int_t idim = 0; idim < nDimSecvtx; idim++){
1644     fSecvtxQA->SetBinEdges(idim, binEdgesSecvtx[idim]);
1645   }
1646
1647   fSecVtxList->AddAt(fSecvtxQA,1);
1648
1649   for(Int_t ivar = 0; ivar < nDimPair; ivar++)
1650     delete binEdgesPair[ivar];
1651   for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
1652     delete binEdgesSecvtx[ivar];
1653 }
1654
1655 //____________________________________________________________
1656 void AliHFEsecVtx::MakeHistos(Int_t step){
1657
1658   //
1659   // make container
1660   //
1661
1662   TString hname=Form("step%d",step);
1663   step = step*7;
1664
1665   const Double_t kPtbound[2] = {0.1, 20.};
1666   Int_t iBin[1];
1667   iBin[0] = 44; // bins in pt
1668   Double_t* binEdges[1];
1669   binEdges[0] =  AliHFEtools::MakeLogarithmicBinning(iBin[0], kPtbound[0], kPtbound[1]);
1670
1671   fSecVtxList->AddAt(new TH1F(hname+"taggedElec", "pT of e", iBin[0],binEdges[0]), step);
1672   fSecVtxList->AddAt(new TH1F(hname+"charmElec", "pT of e", iBin[0],binEdges[0]), step+1);
1673   fSecVtxList->AddAt(new TH1F(hname+"beautyElec", "pT of e", iBin[0],binEdges[0]), step+2);
1674   fSecVtxList->AddAt(new TH1F(hname+"conversionElec", "pT of e", iBin[0],binEdges[0]), step+3);
1675   fSecVtxList->AddAt(new TH1F(hname+"ebgElec", "pT of e", iBin[0],binEdges[0]), step+4);
1676   fSecVtxList->AddAt(new TH1F(hname+"hcontaminElec", "pT of e", iBin[0],binEdges[0]), step+5);
1677   fSecVtxList->AddAt(new TH1F(hname+"elseElec", "pT of e", iBin[0],binEdges[0]), step+6);
1678 }
1679
1680 //____________________________________________________________
1681 void AliHFEsecVtx::FillHistos(Int_t step, const AliESDtrack *track){
1682
1683   //
1684   // make container
1685   //
1686
1687   step = step*7;
1688
1689   AliMCParticle *mctrack = NULL;
1690   TParticle* mcpart = NULL;
1691
1692   (dynamic_cast<TH1F *>(fSecVtxList->At(step)))->Fill(track->Pt()); // electrons tagged
1693
1694   if(HasMCData()){
1695     if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) return;
1696     mcpart = mctrack->Particle();
1697
1698     Int_t esource=fMCQA->GetElecSource(mcpart);
1699     if(esource==1) {
1700       (dynamic_cast<TH1F *>(fSecVtxList->At(step+1)))->Fill(mcpart->Pt()); //charm
1701     }
1702     else if(esource==2 || esource==3) {
1703       (dynamic_cast<TH1F *>(fSecVtxList->At(step+2)))->Fill(mcpart->Pt()); //beauty
1704     }
1705     else if(esource==4) {
1706       (dynamic_cast<TH1F *>(fSecVtxList->At(step+3)))->Fill(mcpart->Pt()); //conversion
1707     }
1708     else if(esource==7) {
1709       (dynamic_cast<TH1F *>(fSecVtxList->At(step+5)))->Fill(mcpart->Pt()); //contamination
1710     }
1711     else if(!(esource<0)) {
1712       (dynamic_cast<TH1F *>(fSecVtxList->At(step+4)))->Fill(mcpart->Pt()); //e backgrounds
1713     }
1714     else {
1715       (dynamic_cast<TH1F *>(fSecVtxList->At(step+6)))->Fill(mcpart->Pt()); //something else?
1716     }
1717   }
1718
1719 }