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