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