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