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