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 <AliESDEvent.h>
  29 #include <AliAODEvent.h>
  30 #include <AliVTrack.h>
  31 #include <AliESDtrack.h>
  32 #include <AliAODTrack.h>
  33 #include "AliHFEsecVtx.h"
  34 #include <AliKFParticle.h>
  35 #include <AliKFVertex.h>
  36 #include <AliLog.h>
  37 #include <AliStack.h>
  38 #include <AliAODMCParticle.h>
  39 #include "AliHFEpairs.h"
  40 #include "AliHFEsecVtxs.h"
  41 #include "AliHFEtrackFilter.h"
  42
  43 ClassImp(AliHFEsecVtx)
  44
  45 //_______________________________________________________________________________________________
  46 AliHFEsecVtx::AliHFEsecVtx():
  47   fFilter(0x0)
  48   ,fESD1(0x0)
  49   ,fAOD1(0x0)
  50   ,fStack(0x0)
  51   ,fUseMCPID(kFALSE)
  52   ,fkSourceLabel()
  53   ,fNparents(0)
  54   ,fParentSelect()
  55   ,fPtRng()
  56   ,fDcaCut()
  57   ,fNoOfHFEpairs(0)
  58   ,fNoOfHFEsecvtxs(0)
  59   ,fHFEpairs(0x0)
  60   ,fHFEsecvtxs(0x0)
  61   ,fMCArray(0x0)
  62   ,fPVx(0)
  63   ,fPVy(0)
  64   ,fCosPhi(-1)
  65   ,fSignedLxy(-1)
  66   ,fKFchi2(-1)
  67   ,fInvmass(-1)
  68   ,fInvmassSigma(-1)
  69   ,fKFip(0)
  70   ,fPairQA(0x0)
  71   ,fSecvtxQA(0x0)
  72   ,fSecVtxList(0x0)
  73 {
  74   //
  75   // Default constructor
  76   //
  77
  78   Init();
  79 }
  80
  81 //_______________________________________________________________________________________________
  82 AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
  83   TObject(p)
  84   ,fFilter(0x0)
  85   ,fESD1(0x0)
  86   ,fAOD1(0x0)
  87   ,fStack(0x0)
  88   ,fUseMCPID(p.fUseMCPID)
  89   ,fkSourceLabel()
  90   ,fNparents(p.fNparents)
  91   ,fParentSelect()
  92   ,fPtRng()
  93   ,fDcaCut()
  94   ,fNoOfHFEpairs(p.fNoOfHFEpairs)
  95   ,fNoOfHFEsecvtxs(p.fNoOfHFEsecvtxs)
  96   ,fHFEpairs(0x0)
  97   ,fHFEsecvtxs(0x0)
  98   ,fMCArray(0x0)
  99   ,fPVx(p.fPVx)
 100   ,fPVy(p.fPVy)
 101   ,fCosPhi(p.fCosPhi)
 102   ,fSignedLxy(p.fSignedLxy)
 103   ,fKFchi2(p.fKFchi2)
 104   ,fInvmass(p.fInvmass)
 105   ,fInvmassSigma(p.fInvmassSigma)
 106   ,fKFip(p.fKFip)
 107   ,fPairQA(0x0)
 108   ,fSecvtxQA(0x0)
 109   ,fSecVtxList(0x0)
 110 {
 111   //
 112   // Copy constructor
 113   //
 114   fFilter = new AliHFEtrackFilter(*p.fFilter);
 115 }
 116
 117 //_______________________________________________________________________________________________
 118 AliHFEsecVtx&
 119 AliHFEsecVtx::operator=(const AliHFEsecVtx &)
 120 {
 121   //
 122   // Assignment operator
 123   //
 124
 125   AliInfo("Not yet implemented.");
 126   return *this;
 127 }
 128
 129 //_______________________________________________________________________________________________
 130 AliHFEsecVtx::~AliHFEsecVtx()
 131 {
 132   //
 133   // Destructor
 134   //
 135
 136   //cout << "Analysis Done." << endl;
 137   delete fFilter;
 138 }
 139
 140 //__________________________________________
 141 void AliHFEsecVtx::Init()
 142 {
 143   //
 144   // set pdg code and index
 145   //
 146
 147   fNparents = 7;
 148
 149   fParentSelect[0][0] =  411;
 150   fParentSelect[0][1] =  421;
 151   fParentSelect[0][2] =  431;
 152   fParentSelect[0][3] = 4122;
 153   fParentSelect[0][4] = 4132;
 154   fParentSelect[0][5] = 4232;
 155   fParentSelect[0][6] = 4332;
 156
 157   fParentSelect[1][0] =  511;
 158   fParentSelect[1][1] =  521;
 159   fParentSelect[1][2] =  531;
 160   fParentSelect[1][3] = 5122;
 161   fParentSelect[1][4] = 5132;
 162   fParentSelect[1][5] = 5232;
 163   fParentSelect[1][6] = 5332;
 164
 165   // momentum ranges to apply pt dependent cuts
 166   fPtRng[0] = 1.0;
 167   fPtRng[1] = 2.0;
 168   fPtRng[2] = 2.5;
 169   fPtRng[3] = 3.0;
 170   fPtRng[4] = 5.0;
 171   fPtRng[5] = 12.0;
 172   fPtRng[6] = 20.0;
 173
 174   // momentum dependent DCA cut values (preliminary)
 175   fDcaCut[0] = 0.04;
 176   fDcaCut[1] = 0.03;
 177   fDcaCut[2] = 0.02;
 178   fDcaCut[3] = 0.015;
 179   fDcaCut[4] = 0.01;
 180   fDcaCut[5] = 0.005;
 181
 182   fFilter = new AliHFEtrackFilter("SecVtxFilter");
 183   fFilter->MakeCutStepRecKineITSTPC();
 184   fFilter->MakeCutStepPrimary();
 185 }
 186
 187 void AliHFEsecVtx::Process(AliVTrack *signalTrack){
 188   if(signalTrack->Pt() < 1.0) return;
 189   AliESDtrack *track = dynamic_cast<AliESDtrack *>(signalTrack);
 190   InitHFEpairs();
 191   InitHFEsecvtxs();
 192   AliESDtrack *htrack = 0x0;
 193   fFilter->Flush();
 194   fFilter->FilterTracks(fESD1);
 195   TObjArray *candidates = fFilter->GetFilteredTracks();
 196   TIterator *trackIter = candidates->MakeIterator();
 197   while((htrack = dynamic_cast<AliESDtrack *>(trackIter->Next()))){
 198     if(track->GetID() == htrack->GetID()) continue; // since it is for tagging single electron, don't need additional condition
 199     if (htrack->Pt()<1.0) continue;
 200     PairAnalysis(track, htrack, htrack->GetID()); // e-h pairing
 201   }
 202   delete trackIter;
 203   /*for(int ip=0; ip<fSecVtx->HFEpairs()->GetEntriesFast(); ip++){
 204       if(HasMCData()){
 205         AliHFEpairs *pair = (AliHFEpairs*) (fSecVtx->HFEpairs()->UncheckedAt(ip));
 206         if(!(pair->GetPairCode()>1. && pair->GetPairCode()<4.))  // apply various cuts
 207         fSecVtx->HFEpairs()->RemoveAt(ip);
 208       }
 209     }*/
 210   HFEpairs()->Compress();
 211   RunSECVTX(track); // secondary vertexing with e,h1,h2,.. tracks
 212   for(int ip=0; ip<HFEsecvtxs()->GetEntriesFast(); ip++){
 213     AliHFEsecVtxs *secvtx=0x0;
 214     secvtx = (AliHFEsecVtxs*) (HFEsecvtxs()->UncheckedAt(ip));
 215     // here you apply cuts, then if it doesn't pass the cut, remove it from the fSecVtx->HFEsecvtxs()
 216   }
 217   DeleteHFEpairs();
 218   DeleteHFEsecvtxs();
 219 }
 220
 221 //_______________________________________________________________________________________________
 222 void AliHFEsecVtx::CreateHistograms(TList *qaList)
 223 {
 224   //
 225   // create histograms
 226   //
 227
 228   fSecVtxList = new TList;
 229   fSecVtxList->SetName("SecVtx");
 230
 231   MakeContainer();
 232   /*
 233   fkSourceLabel[kAll]="all";
 234   fkSourceLabel[kDirectCharm]="directCharm";
 235   fkSourceLabel[kDirectBeauty]="directBeauty";
 236   fkSourceLabel[kBeautyCharm]="beauty2charm";
 237   fkSourceLabel[kGamma]="gamma";
 238   fkSourceLabel[kPi0]="pi0";
 239   fkSourceLabel[kElse]="others";
 240   fkSourceLabel[kBeautyGamma]="beauty22gamma";
 241   fkSourceLabel[kBeautyPi0]="beauty22pi0";
 242   fkSourceLabel[kBeautyElse]="beauty22others";
 243
 244   TString hname;
 245   TString hnopt="secvtx_";
 246   for (Int_t isource = 0; isource < 10; isource++ ){
 247   }
 248   */
 249
 250   qaList->AddLast(fSecVtxList);
 251 }
 252
 253 //_______________________________________________________________________________________________
 254 void AliHFEsecVtx::GetPrimaryCondition()
 255 {
 256   //
 257   // get primary characteristics and set
 258   //
 259
 260   if (fESD1) {
 261     AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
 262     if( primVtxCopy.GetNDF() <1 ) return;
 263     fPVx = primVtxCopy.GetX();
 264     fPVx = primVtxCopy.GetY();
 265   }
 266   else if(fAOD1) {
 267     AliKFVertex primVtxCopy(*(fAOD1->GetPrimaryVertex()));
 268     if( primVtxCopy.GetNDF() <1 ) return;
 269     fPVx = primVtxCopy.GetX();
 270     fPVx = primVtxCopy.GetY();
 271   }
 272 }
 273
 274 //_______________________________________________________________________________________________
 275 void AliHFEsecVtx::PairAnalysis(AliVTrack* track1, AliVTrack* track2, Int_t index2)
 276 {
 277   //
 278   // calculate e-h pair characteristics and tag pair
 279   //
 280
 281   //later consider the below
 282   Float_t dca1[2]={-999.,-999.}, dca2[2]={-999.,-999.};
 283   Float_t cov1[3]={-999.,-999.,-999.}, cov2[3]={-999.,-999.,-999.};
 284
 285   if (IsAODanalysis()){
 286     AliESDtrack esdTrk1(track1);
 287     AliESDtrack esdTrk2(track2);
 288     esdTrk1.PropagateToDCA(fAOD1->GetPrimaryVertex(),0,10000,(Double_t*)dca1,(Double_t*)cov1);
 289     esdTrk2.PropagateToDCA(fAOD1->GetPrimaryVertex(),0,10000,(Double_t*)dca2,(Double_t*)cov2);
 290   }
 291   else {
 292     ((AliESDtrack*)track1)->GetImpactParameters(dca1,cov1);
 293     ((AliESDtrack*)track2)->GetImpactParameters(dca2,cov2);
 294   }
 295
 296   // apply pt dependent dca cut on hadrons
 297   for(int ibin=0; ibin<6; ibin++){
 298     if((track2->Pt()>fPtRng[ibin] && track2->Pt()<fPtRng[ibin+1]) && TMath::Abs(dca2[0])<fDcaCut[ibin]) return;
 299   }
 300
 301   // get KF particle input pid
 302   Int_t pdg1 = GetPDG(track1);
 303   Int_t pdg2 = GetPDG(track2);
 304
 305   if(pdg1==-1 || pdg2==-1) {
 306     //printf("out if considered pid range \n");
 307     return;
 308   }
 309
 310   // create KF particle of pair
 311   if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
 312   else AliKFParticle::SetField(fESD1->GetMagneticField());
 313   AliKFParticle kfTrack1(*track1, pdg1);
 314   AliKFParticle kfTrack2(*track2, pdg2);
 315
 316   AliKFParticle kfSecondary(kfTrack1,kfTrack2);
 317
 318   //secondary vertex point from kf particle
 319   Double_t kfx = kfSecondary.GetX();
 320   Double_t kfy = kfSecondary.GetY();
 321   //Double_t kfz = kfSecondary.GetZ();
 322
 323   //momentum at the decay point from kf particle
 324   Double_t kfpx = kfSecondary.GetPx();
 325   Double_t kfpy = kfSecondary.GetPy();
 326   //Double_t kfpz = kfSecondary.GetPz();
 327
 328   Double_t dx = kfx-fPVx;
 329   Double_t dy = kfy-fPVy;
 330
 331   // discriminating variables
 332   // invariant mass of the KF particle
 333   Double_t invmass = -1;
 334   Double_t invmassSigma = -1;
 335   kfSecondary.GetMass(invmass,invmassSigma);
 336
 337   // chi2 of the KF particle
 338   Double_t kfchi2 = -1;
 339   if(kfSecondary.GetNDF()>0) kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 340
 341   // opening angle between two particles in XY plane
 342   Double_t phi = kfTrack1.GetAngleXY(kfTrack2);
 343   Double_t cosphi = TMath::Cos(phi);
 344
 345   // projection of kf vertex vector to the kf momentum direction
 346   Double_t signedLxy=-999.;
 347   if((dx*kfpx+dy*kfpy)>0) signedLxy = TMath::Sqrt(dx*dx+dy*dy);
 348   if((dx*kfpx+dy*kfpy)<0) signedLxy = -1*TMath::Sqrt(dx*dx+dy*dy);
 349   //[the other way to think about]
 350   //Double_t psqr = kfpx*kfpx+kfpy*kfpy;
 351   //if(psqr>0) signedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
 352
 353   // DCA from primary to e-h KF particle (impact parameter of KF particle)
 354   Double_t vtx[2]={fPVx, fPVy};
 355   Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
 356
 357   Int_t paircode = -1;
 358   if (HasMCData()) paircode = GetPairCode(track1,track2);
 359
 360   AliHFEpairs hfepair;
 361   hfepair.SetTrkLabel(index2);
 362   hfepair.SetInvmass(invmass);
 363   hfepair.SetKFChi2(kfchi2);
 364   hfepair.SetOpenangle(phi);
 365   hfepair.SetCosOpenangle(cosphi);
 366   hfepair.SetSignedLxy(signedLxy);
 367   hfepair.SetKFIP(kfip);
 368   hfepair.SetPairCode(paircode);
 369   AddHFEpairToArray(&hfepair);
 370   fNoOfHFEpairs++;
 371
 372   // fill into container for later QA
 373   Double_t dataE[6];
 374   dataE[0]=invmass;
 375   dataE[1]=kfchi2;
 376   dataE[2]=phi;
 377   dataE[3]=signedLxy;
 378   dataE[4]=kfip;
 379   dataE[5]=paircode;
 380   /*
 381   dataE[6]=TMath::Abs(dca1[0]);
 382   dataE[7]=TMath::Abs(dca2[0]);
 383   //if(cov1[0]>0) dataE[6]=Double_t(dca1[0]/cov1[0]);
 384   //if(cov2[0]>0) dataE[7]=Double_t(dca2[0]/cov2[0]);
 385   dataE[8]=track1->Pt();
 386   dataE[9]=track2->Pt();
 387   */
 388   fPairQA->Fill(dataE);
 389
 390 }
 391
 392 //_______________________________________________________________________________________________
 393 void AliHFEsecVtx::RunSECVTX(AliVTrack *track)
 394 {
 395   //
 396   // run secondary vertexing algorithm and do tagging
 397   //
 398
 399   //printf("number of considered pairs= %d\n",HFEpairs()->GetEntriesFast());
 400   FindSECVTXCandid(track);
 401 }
 402
 403 //_______________________________________________________________________________________________
 404 void AliHFEsecVtx::FindSECVTXCandid(AliVTrack *track)
 405 {
 406   //
 407   // find secondary vertex candidate and store them into container
 408   //
 409
 410   AliVTrack *htrack[20];
 411   Int_t htracklabel[20];
 412   Double_t vtxchi2cut=3.; // testing cut
 413   Double_t dataE[6]={-999.,-999.,-999.,-999.,-1.,0};
 414   if (HFEpairs()->GetEntriesFast()>20){
 415     AliDebug(3, "number of paired hadron is over maximum(20)");
 416     return;
 417   }
 418
 419   // get paired track objects
 420   AliHFEpairs *pair=0x0;
 421   if (IsAODanalysis()){
 422     for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
 423        pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
 424        htracklabel[ip] = pair->GetTrkLabel();
 425        htrack[ip] = fAOD1->GetTrack(pair->GetTrkLabel());
 426     }
 427   }
 428   else{
 429     for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
 430        pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
 431        htracklabel[ip] = pair->GetTrkLabel();
 432        htrack[ip] = fESD1->GetTrack(pair->GetTrkLabel());
 433     }
 434   }
 435   // in case there is only one paired track with the electron, put pair characteristics into secvtx container
 436   // for the moment, I only apply pair vertex chi2 cut
 437   if (HFEpairs()->GetEntriesFast() == 1){
 438     if (pair->GetKFChi2()<vtxchi2cut) { // you can also put single track cut
 439       AliHFEsecVtxs hfesecvtx;
 440       hfesecvtx.SetTrkLabel1(pair->GetTrkLabel());
 441       hfesecvtx.SetTrkLabel2(-999);
 442       hfesecvtx.SetInvmass(pair->GetInvmass());
 443       hfesecvtx.SetKFChi2(pair->GetKFChi2());
 444       hfesecvtx.SetSignedLxy(pair->GetSignedLxy());
 445       hfesecvtx.SetKFIP(pair->GetKFIP());
 446       AddHFEsecvtxToArray(&hfesecvtx);
 447       fNoOfHFEsecvtxs++;
 448
 449       dataE[0]=pair->GetInvmass();
 450       dataE[1]=pair->GetKFChi2();
 451       dataE[2]=pair->GetSignedLxy();
 452       dataE[3]=pair->GetKFIP();
 453       if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
 454       dataE[5]=2;
 455       fSecvtxQA->Fill(dataE);
 456     }
 457     return;
 458   }
 459
 460   // in case there are multiple paired track with the electron, calculate secvtx characteristics
 461   // put the secvtx characteristics into container if it passes cuts
 462   for (int i=0; i<HFEpairs()->GetEntriesFast()-1; i++){
 463      for (int j=i+1; j<HFEpairs()->GetEntriesFast(); j++){
 464         CalcSECVTXProperty(track, htrack[i], htrack[j]);
 465         if (fKFchi2<vtxchi2cut) {
 466           AliHFEsecVtxs hfesecvtx;
 467           hfesecvtx.SetTrkLabel1(htracklabel[i]);
 468           hfesecvtx.SetTrkLabel2(htracklabel[j]);
 469           hfesecvtx.SetKFChi2(fKFchi2);
 470           hfesecvtx.SetInvmass(fInvmass);
 471           hfesecvtx.SetSignedLxy(fSignedLxy);
 472           hfesecvtx.SetKFIP(fKFip);
 473           AddHFEsecvtxToArray(&hfesecvtx);
 474           fNoOfHFEsecvtxs++;
 475
 476           dataE[0]=fInvmass;
 477           dataE[1]=fKFchi2;
 478           dataE[2]=fSignedLxy;
 479           dataE[3]=fKFip;
 480           if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
 481           dataE[5]=3;
 482           fSecvtxQA->Fill(dataE);
 483         }
 484      }
 485   }
 486 }
 487
 488 //_______________________________________________________________________________________________
 489 void AliHFEsecVtx::CalcSECVTXProperty(AliVTrack* track1, AliVTrack* track2, AliVTrack* track3)
 490 {
 491   //
 492   // calculate secondary vertex properties
 493   //
 494
 495   // get KF particle input pid
 496   Int_t pdg1 = GetPDG(track1);
 497   Int_t pdg2 = GetPDG(track2);
 498   Int_t pdg3 = GetPDG(track3);
 499
 500   if(pdg1==-1 || pdg2==-1 || pdg3==-1) {
 501     //printf("out if considered pid range \n");
 502     return;
 503   }
 504
 505   // create KF particle of pair
 506   if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
 507   else AliKFParticle::SetField(fESD1->GetMagneticField());
 508   AliKFParticle kfTrack1(*track1, pdg1);
 509   AliKFParticle kfTrack2(*track2, pdg2);
 510   AliKFParticle kfTrack3(*track3, pdg3);
 511
 512   AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
 513
 514   //secondary vertex point from kf particle
 515   Double_t kfx = kfSecondary.GetX();
 516   Double_t kfy = kfSecondary.GetY();
 517   //Double_t kfz = kfSecondary.GetZ();
 518
 519   //momentum at the decay point from kf particle
 520   Double_t kfpx = kfSecondary.GetPx();
 521   Double_t kfpy = kfSecondary.GetPy();
 522   //Double_t kfpz = kfSecondary.GetPz();
 523
 524   Double_t dx = kfx-fPVx;
 525   Double_t dy = kfy-fPVy;
 526
 527   // discriminating variables ----------------------------------------------------------
 528
 529   if(kfSecondary.GetNDF()>0) fKFchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
 530
 531   // invariant mass of the KF particle
 532   kfSecondary.GetMass(fInvmass,fInvmassSigma);
 533
 534   // DCA from primary to e-h KF particle (impact parameter of KF particle)
 535   Double_t vtx[2]={fPVx, fPVy};
 536   fKFip = kfSecondary.GetDistanceFromVertexXY(vtx);
 537
 538   if((dx*kfpx+dy*kfpy)>0) fSignedLxy= TMath::Sqrt(dx*dx+dy*dy);
 539   if((dx*kfpx+dy*kfpy)<0) fSignedLxy= -1*TMath::Sqrt(dx*dx+dy*dy);
 540   //[the other way to think about] - projection of kf vertex vector to the kf momentum direction
 541   //Double_t psqr = kfpx*kfpx+kfpy*kfpy;
 542   //if(psqr>0) fSignedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
 543 }
 544
 545 //_______________________________________________________________________________________________
 546 Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track)
 547 {
 548   //
 549   // return mc pid
 550   //
 551
 552   Int_t label = TMath::Abs(track->GetLabel());
 553   TParticle* mcpart = fStack->Particle(label);
 554   if ( !mcpart ) return 0;
 555   Int_t pdgCode = mcpart->GetPdgCode();
 556
 557   return pdgCode;
 558 }
 559
 560 //_______________________________________________________________________________________________
 561 Int_t AliHFEsecVtx::GetPairOriginESD(AliESDtrack* trk1, AliESDtrack* trk2)
 562 {
 563   //
 564   // return pdg code of the origin(source) of the pair
 565   //
 566   //
 567   // ---*---*---*-----ancester A----- track1
 568   //                        |____*______
 569   //                             |______ track2
 570   // => if they originated from same ancester,
 571   //    then return "the absolute value of pdg code of ancester A"
 572   //
 573   // ---*---*---B hadron-----ancester A----- track1
 574   //                               |____*______
 575   //                                    |______ track2
 576   // => if they originated from same ancester, and this ancester originally comes from B hadrons
 577   //    then return -1*"the absolute value of pdg code of ancester A"
 578   //
 579   // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
 580   //
 581
 582   if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
 583   TParticle* part1 = fStack->Particle(trk1->GetLabel());
 584   TParticle* part2 = fStack->Particle(trk2->GetLabel());
 585   TParticle* part2cp = part2;
 586   if (!(part1) || !(part2)) return 0;
 587
 588   Int_t srcpdg = 0;
 589
 590   //if the two tracks' mother's label is same, get the mother info
 591   //in case of charm, check if it originated from beauty
 592   for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
 593      Int_t label1 = part1->GetFirstMother();
 594      if (label1 < 0) return 0;
 595
 596      for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
 597         Int_t label2 = part2->GetFirstMother();
 598         if (label2 < 0) break;
 599
 600         if (label1 == label2){ //check if two tracks are originated from same mother
 601           TParticle* commonmom = fStack->Particle(label2);
 602           srcpdg = abs(commonmom->GetPdgCode());
 603
 604           //check ancester to see if it is originally from beauty
 605           for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
 606              Int_t ancesterlabel = commonmom->GetFirstMother();
 607              if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
 608
 609              TParticle* commonancester = fStack->Particle(ancesterlabel);
 610              Int_t ancesterpdg = abs(commonancester->GetPdgCode());
 611
 612              for (Int_t l=0; l<fNparents; l++){
 613                 if (abs(ancesterpdg)==fParentSelect[1][l]){
 614                   srcpdg = -1*srcpdg; //multiply -1 for hadron from bottom
 615                   return srcpdg;
 616                 }
 617              }
 618              commonmom = commonancester;
 619           }
 620         }
 621         part2 = fStack->Particle(label2); //if their mother is different, go to earlier generation of 2nd particle
 622         if (!(part2)) break;
 623      }
 624      part1 = fStack->Particle(label1); //if their mother is different, go to earlier generation of 1st particle
 625      part2 = part2cp;
 626      if (!(part1)) return 0;
 627   }
 628
 629   return srcpdg;
 630 }
 631
 632 //_______________________________________________________________________________________________
 633 Int_t AliHFEsecVtx::GetPairOriginAOD(AliAODTrack* trk1, AliAODTrack* trk2)
 634 {
 635
 636   //
 637   // return pdg code of the origin(source) of the pair
 638   //
 639   //
 640   // ---*---*---*-----ancester A----- track1
 641   //                        |____*______
 642   //                             |______ track2
 643   // => if they originated from same ancester,
 644   //    then return "the absolute value of pdg code of ancester A"
 645   //
 646   // ---*---*---B hadron-----ancester A----- track1
 647   //                               |____*______
 648   //                                    |______ track2
 649   // => if they originated from same ancester, and this ancester originally comes from B hadrons
 650   //    then return -1*"the absolute value of pdg code of ancester A"
 651   //
 652   // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
 653   //
 654
 655   if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
 656   AliAODMCParticle *part1 = (AliAODMCParticle*)fMCArray->At(trk1->GetLabel());
 657   AliAODMCParticle *part2 = (AliAODMCParticle*)fMCArray->At(trk2->GetLabel());
 658   AliAODMCParticle *part2cp = part2;
 659   if (!(part1) || !(part2)) return 0;
 660
 661   Int_t srcpdg = 0;
 662
 663   //if the two tracks' mother's label is same, get the mother info
 664   //in case of charm, check if it originated from beauty
 665   for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
 666      Int_t label1 = part1->GetMother();
 667      if (label1 < 0) return 0;
 668
 669      for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
 670         Int_t label2 = part2->GetMother();
 671         if (label2 < 0) return 0;
 672
 673         if (label1 == label2){ //check if two tracks are originated from same mother
 674           AliAODMCParticle *commonmom = (AliAODMCParticle*)fMCArray->At(label1);
 675           srcpdg = abs(commonmom->GetPdgCode());
 676
 677           //check ancester to see if it is originally from beauty
 678           for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
 679              Int_t ancesterlabel = commonmom->GetMother();
 680              if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
 681
 682              AliAODMCParticle *commonancester = (AliAODMCParticle*)fMCArray->At(ancesterlabel);
 683              Int_t ancesterpdg = abs(commonancester->GetPdgCode());
 684
 685              for (Int_t l=0; l<fNparents; l++){
 686                 if (abs(ancesterpdg)==fParentSelect[1][l]){
 687                   srcpdg = -1*srcpdg; //multiply -1 for charm from bottom
 688                   return srcpdg;
 689                 }
 690              }
 691              commonmom = commonancester;
 692           }
 693         }
 694         part2 = (AliAODMCParticle*)fMCArray->At(label2); //if their mother is different, go to earlier generation of 2nd particle
 695         if (!(part2)) break;
 696      }
 697      part1 = (AliAODMCParticle*)fMCArray->At(label1); //if their mother is different, go to earlier generation of 2nd particle
 698      part2 = part2cp;
 699      if (!(part1)) return 0;
 700   }
 701
 702   return srcpdg;
 703 }
 704
 705 //_______________________________________________________________________________________________
 706 Int_t AliHFEsecVtx::GetPairCode(const AliVTrack* const trk1, const AliVTrack* const trk2)
 707 {
 708   //
 709   // return pair code which is predefinded as:
 710   //  kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0,
 711   //  kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
 712   //
 713
 714   Int_t srcpdg = -1;
 715   Int_t srccode = kElse;
 716
 717   if(IsAODanalysis()) srcpdg = GetPairOriginAOD((AliAODTrack*)trk1,(AliAODTrack*)trk2);
 718   else srcpdg = GetPairOriginESD((AliESDtrack*)trk1,(AliESDtrack*)trk2);
 719
 720   if (srcpdg < 0) srccode = kBeautyElse;
 721   for (Int_t i=0; i<fNparents; i++){
 722     if (abs(srcpdg)==fParentSelect[0][i]) {
 723       if (srcpdg>0) srccode = kDirectCharm;
 724       if (srcpdg<0) srccode = kBeautyCharm;
 725     }
 726     if (abs(srcpdg)==fParentSelect[1][i]) {
 727       if (srcpdg>0) srccode = kDirectBeauty;
 728       if (srcpdg<0)  return kElse;
 729     }
 730   }
 731   if (srcpdg == 22) srccode = kGamma;
 732   if (srcpdg == -22) srccode = kBeautyGamma;
 733   if (srcpdg == 111) srccode = kPi0;
 734   if (srcpdg == -111) srccode = kBeautyPi0;
 735
 736   return srccode;
 737 }
 738
 739 //_______________________________________________________________________________________________
 740 Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack)
 741 {
 742   //
 743   // return decay electron's origin
 744   //
 745
 746   if (iTrack < 0) {
 747     AliDebug(1, "Stack label is negative, return\n");
 748     return -1;
 749   }
 750
 751   TParticle* mcpart = fStack->Particle(iTrack);
 752
 753 //  if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !
 754
 755   Int_t iLabel = mcpart->GetFirstMother();
 756   if (iLabel<0){
 757     AliDebug(1, "Stack label is negative, return\n");
 758     return -1;
 759   }
 760
 761   Int_t origin = -1;
 762   Bool_t isFinalOpenCharm = kFALSE;
 763
 764   TParticle *partMother = fStack->Particle(iLabel);
 765   Int_t maPdgcode = partMother->GetPdgCode();
 766
 767   // if the mother is charmed hadron
 768   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
 769
 770     for (Int_t i=0; i<fNparents; i++){
 771       if (abs(maPdgcode)==fParentSelect[0][i]){
 772         isFinalOpenCharm = kTRUE;
 773       }
 774     }
 775     if (!isFinalOpenCharm) return -1;
 776
 777     // iterate until you find B hadron as a mother or become top ancester
 778     for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
 779
 780       Int_t jLabel = partMother->GetFirstMother();
 781       if (jLabel == -1){
 782         origin = kDirectCharm;
 783         return origin;
 784       }
 785       if (jLabel < 0){ // safety protection
 786         AliDebug(1, "Stack label is negative, return\n");
 787         return -1;
 788       }
 789
 790       // if there is an ancester
 791       TParticle* grandMa = fStack->Particle(jLabel);
 792       Int_t grandMaPDG = grandMa->GetPdgCode();
 793
 794       for (Int_t j=0; j<fNparents; j++){
 795         if (abs(grandMaPDG)==fParentSelect[1][j]){
 796           origin = kBeautyCharm;
 797           return origin;
 798         }
 799       }
 800
 801       partMother = grandMa;
 802     } // end of iteration
 803   } // end of if
 804   else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
 805     for (Int_t i=0; i<fNparents; i++){
 806       if (abs(maPdgcode)==fParentSelect[1][i]){
 807         origin = kDirectBeauty;
 808         return origin;
 809       }
 810     }
 811   } // end of if
 812
 813   //============ gamma ================
 814   else if ( abs(maPdgcode) == 22 ) {
 815     origin = kGamma;
 816
 817     // iterate until you find B hadron as a mother or become top ancester
 818     for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
 819
 820       Int_t jLabel = partMother->GetFirstMother();
 821       if (jLabel == -1){
 822         origin = kGamma;
 823         return origin;
 824       }
 825       if (jLabel < 0){ // safety protection
 826         AliDebug(1, "Stack label is negative, return\n");
 827         return -1;
 828       }
 829
 830       // if there is an ancester
 831       TParticle* grandMa = fStack->Particle(jLabel);
 832       Int_t grandMaPDG = grandMa->GetPdgCode();
 833
 834       for (Int_t j=0; j<fNparents; j++){
 835         if (abs(grandMaPDG)==fParentSelect[1][j]){
 836           origin = kBeautyGamma;
 837           return origin;
 838         }
 839       }
 840
 841       partMother = grandMa;
 842     } // end of iteration
 843
 844     return origin;
 845   } // end of if
 846
 847   //============ pi0 ================
 848   else if ( abs(maPdgcode) == 111 ) {
 849     origin = kPi0;
 850
 851     // iterate until you find B hadron as a mother or become top ancester
 852     for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
 853
 854       Int_t jLabel = partMother->GetFirstMother();
 855       if (jLabel == -1){
 856         origin = kPi0;
 857         return origin;
 858       }
 859       if (jLabel < 0){ // safety protection
 860         AliDebug(1, "Stack label is negative, return\n");
 861         return -1;
 862       }
 863
 864       // if there is an ancester
 865       TParticle* grandMa = fStack->Particle(jLabel);
 866       Int_t grandMaPDG = grandMa->GetPdgCode();
 867
 868       for (Int_t j=0; j<fNparents; j++){
 869         if (abs(grandMaPDG)==fParentSelect[1][j]){
 870           origin = kBeautyPi0;
 871           return origin;
 872         }
 873       }
 874
 875       partMother = grandMa;
 876     } // end of iteration
 877
 878     return origin;
 879   } // end of if
 880
 881   else {
 882     origin = kElse;
 883     // iterate until you find B hadron as a mother or become top ancester
 884     for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
 885
 886       Int_t jLabel = partMother->GetFirstMother();
 887       if (jLabel == -1){
 888         origin = kElse;
 889         return origin;
 890       }
 891       if (jLabel < 0){ // safety protection
 892         AliDebug(1, "Stack label is negative, return\n");
 893         return -1;
 894       }
 895
 896       // if there is an ancester
 897       TParticle* grandMa = fStack->Particle(jLabel);
 898       Int_t grandMaPDG = grandMa->GetPdgCode();
 899
 900       for (Int_t j=0; j<fNparents; j++){
 901         if (abs(grandMaPDG)==fParentSelect[1][j]){
 902           origin = kBeautyElse;
 903           return origin;
 904         }
 905       }
 906
 907       partMother = grandMa;
 908     } // end of iteration
 909   }
 910
 911   return origin;
 912 }
 913
 914 //_______________________________________________________________________________________________
 915 Int_t AliHFEsecVtx::GetPDG(AliVTrack *track)
 916 {
 917   //
 918   // get KF particle input pdg for mass hypothesis
 919   //
 920
 921   Int_t pdgCode=-1;
 922
 923   if (fUseMCPID && HasMCData()){
 924     pdgCode = GetMCPDG(track);
 925   }
 926   else if(fESD1){
 927     Int_t pid=0;
 928     Double_t prob;
 929     GetESDPID((AliESDtrack*)track, pid, prob);
 930     switch(pid){
 931       case 0:  pdgCode = 11; break;
 932       case 1:  pdgCode = 13; break;
 933       case 2:  pdgCode = 211; break;
 934       case 3:  pdgCode = 321; break;
 935       case 4:  pdgCode = 2212; break;
 936       default: pdgCode = -1;
 937     }
 938   }
 939   else if(fAOD1){
 940     Int_t pid = ((AliAODTrack*)track)->GetMostProbablePID();
 941     switch(pid){
 942       case 0:  pdgCode = 11; break;
 943       case 1:  pdgCode = 13; break;
 944       case 2:  pdgCode = 211; break;
 945       case 3:  pdgCode = 321; break;
 946       case 4:  pdgCode = 2212; break;
 947       default: pdgCode = -1;
 948     }
 949   }
 950
 951   return pdgCode;
 952 }
 953
 954 //_______________________________________________________________________________________________
 955 Int_t AliHFEsecVtx::GetMCPDG(AliVTrack *track)
 956 {
 957   //
 958   // return mc pdg code
 959   //
 960
 961   Int_t label = TMath::Abs(track->GetLabel());
 962   Int_t pdgCode;
 963
 964   if (IsAODanalysis()) {
 965     AliAODMCParticle *mcpart = (AliAODMCParticle*)fMCArray->At(label);
 966     if ( !mcpart ) return 0;
 967       pdgCode = mcpart->GetPdgCode();
 968   }
 969   else {
 970     TParticle* mcpart = fStack->Particle(label);
 971     if ( !mcpart ) return 0;
 972       pdgCode = mcpart->GetPdgCode();
 973   }
 974
 975     return pdgCode;
 976 }
 977
 978 //______________________________________________________________________________
 979 void AliHFEsecVtx::GetESDPID(AliESDtrack *track, Int_t &recpid, Double_t &recprob)
 980 {
 981   //
 982   // calculate likehood for esd pid
 983   //
 984
 985   recpid = -1;
 986   recprob = -1;
 987
 988   Int_t ipid=-1;
 989   Double_t max=0.;
 990
 991   Double_t probability[5];
 992
 993   // get probability of the diffenrent particle types
 994   track->GetESDpid(probability);
 995
 996   // find most probable particle in ESD pid
 997   // 0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
 998   ipid = TMath::LocMax(5,probability);
 999   max = TMath::MaxElement(5,probability);
1000
1001   recpid = ipid;
1002   recprob = max;
1003 }
1004
1005 //_____________________________________________________________________________
1006 void AliHFEsecVtx::AddHFEpairToArray(const AliHFEpairs* const pair)
1007 {
1008   //
1009   // Add a HFE pair to the array
1010   //
1011
1012   Int_t n = HFEpairs()->GetEntriesFast();
1013   if(n!=fNoOfHFEpairs)AliError(Form("fNoOfHFEpairs != HFEpairs()->GetEntriesFast %i != %i \n", fNoOfHFEpairs, n));
1014   new((*HFEpairs())[n]) AliHFEpairs(*pair);
1015 }
1016
1017 //_____________________________________________________________________________
1018 TClonesArray *AliHFEsecVtx::HFEpairs()
1019 {
1020   //
1021   // Returns the list of HFE pairs
1022   //
1023
1024   if (!fHFEpairs) {
1025       fHFEpairs = new TClonesArray("AliHFEpairs", 200);
1026   }
1027   return fHFEpairs;
1028 }
1029
1030 //_____________________________________________________________________________
1031 void AliHFEsecVtx::DeleteHFEpairs()
1032 {
1033   //
1034   // Delete the list of HFE pairs
1035   //
1036
1037   if (fHFEpairs) {
1038     fHFEpairs->Delete();
1039     //delete fHFEpairs;
1040   }
1041 }
1042
1043 //_____________________________________________________________________________
1044 void AliHFEsecVtx::InitHFEpairs()
1045 {
1046   //
1047   // Initialization should be done before make all possible pairs for a given electron candidate
1048   //
1049
1050   fNoOfHFEpairs = 0;
1051 }
1052
1053 //_____________________________________________________________________________
1054 void AliHFEsecVtx::AddHFEsecvtxToArray(const AliHFEsecVtxs* const secvtx)
1055 {
1056   //
1057   // Add a HFE secondary vertex to the array
1058   //
1059
1060   Int_t n = HFEsecvtxs()->GetEntriesFast();
1061   if(n!=fNoOfHFEsecvtxs)AliError(Form("fNoOfHFEsecvtxs != HFEsecvtxs()->GetEntriesFast %i != %i \n", fNoOfHFEsecvtxs, n));
1062   new((*HFEsecvtxs())[n]) AliHFEsecVtxs(*secvtx);
1063 }
1064
1065 //_____________________________________________________________________________
1066 TClonesArray *AliHFEsecVtx::HFEsecvtxs()
1067 {
1068   //
1069   // Returns the list of HFE secvtx
1070   //
1071
1072   if (!fHFEsecvtxs) {
1073       fHFEsecvtxs = new TClonesArray("AliHFEsecVtxs", 200);
1074   }
1075   return fHFEsecvtxs;
1076 }
1077
1078 //_____________________________________________________________________________
1079 void AliHFEsecVtx::DeleteHFEsecvtxs()
1080 {
1081   //
1082   // Delete the list of HFE pairs
1083   //
1084
1085   if (fHFEsecvtxs) {
1086     fHFEsecvtxs->Delete();
1087     //delete fHFEsecvtx;
1088   }
1089 }
1090
1091 //_____________________________________________________________________________
1092 void AliHFEsecVtx::InitHFEsecvtxs()
1093 {
1094   //
1095   // Initialization should be done
1096   //
1097
1098   fNoOfHFEsecvtxs = 0;
1099 }
1100
1101 //_______________________________________________________________________________________________
1102 Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track) const
1103 {
1104   //if (track->Pt() < 1.0) return kFALSE;
1105   //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
1106   //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
1107   //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1108   if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
1109   //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
1110
1111 /*
1112   Float_t dcaR=-1;
1113   Float_t dcaZ=-1;
1114   track->GetImpactParameters(dcaR,dcaZ);
1115   if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
1116   if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
1117 */
1118   return kTRUE;
1119 }
1120 //____________________________________________________________
1121 void AliHFEsecVtx::MakeContainer(){
1122
1123   //
1124   // make container
1125   //
1126
1127   const Int_t nDimPair=6;
1128   Int_t nBinPair[nDimPair] = {200, 500, 314, 2000, 2000, 11};
1129   //Int_t nBinPair[nDimPair] = {200, 500, 314, 2000, 2000, 11, 1000, 1000, 60, 60};
1130   const Double_t kInvmassmin = 0., kInvmassmax = 20.;
1131   const Double_t kKFChi2min = 0, kKFChi2max= 50;
1132   const Double_t kOpenanglemin = 0, kOpenanglemax = 3.14;
1133   const Double_t kSignedLxymin = -10, kSignedLxymax= 10;
1134   const Double_t kKFIPmin = -10, kKFIPmax= 10;
1135   const Double_t kPairCodemin = -1, kPairCodemax= 10;
1136   //const Double_t kDCAsigmin = 0, kDCAsigmax= 5;
1137   //const Double_t kPtmin = 0, kPtmax= 30;
1138
1139   Double_t* binEdgesPair[nDimPair];
1140   for(Int_t ivar = 0; ivar < nDimPair; ivar++)
1141     binEdgesPair[ivar] = new Double_t[nBinPair[ivar] + 1];
1142
1143   for(Int_t i=0; i<=nBinPair[0]; i++) binEdgesPair[0][i]=(Double_t)kInvmassmin + (kInvmassmax - kInvmassmin)/nBinPair[0]*(Double_t)i;
1144   for(Int_t i=0; i<=nBinPair[1]; i++) binEdgesPair[1][i]=(Double_t)kKFChi2min + (kKFChi2max - kKFChi2min)/nBinPair[1]*(Double_t)i;
1145   for(Int_t i=0; i<=nBinPair[2]; i++) binEdgesPair[2][i]=(Double_t)kOpenanglemin + (kOpenanglemax - kOpenanglemin)/nBinPair[2]*(Double_t)i;
1146   for(Int_t i=0; i<=nBinPair[3]; i++) binEdgesPair[3][i]=(Double_t)kSignedLxymin + (kSignedLxymax - kSignedLxymin)/nBinPair[3]*(Double_t)i;
1147   for(Int_t i=0; i<=nBinPair[4]; i++) binEdgesPair[4][i]=(Double_t)kKFIPmin + (kKFIPmax - kKFIPmin)/nBinPair[4]*(Double_t)i;
1148   for(Int_t i=0; i<=nBinPair[5]; i++) binEdgesPair[5][i]=(Double_t)kPairCodemin + (kPairCodemax - kPairCodemin)/nBinPair[5]*(Double_t)i;
1149   /*for(Int_t i=0; i<=nBinPair[6]; i++) binEdgesPair[6][i]=(Double_t)kDCAsigmin + (kDCAsigmax - kDCAsigmin)/nBinPair[6]*(Double_t)i;
1150   for(Int_t i=0; i<=nBinPair[7]; i++) binEdgesPair[7][i]=binEdgesPair[6][i];
1151   for(Int_t i=0; i<=nBinPair[8]; i++) binEdgesPair[8][i]=(Double_t)kPtmin + (kPtmax - kPtmin)/nBinPair[8]*(Double_t)i;
1152   for(Int_t i=0; i<=nBinPair[9]; i++) binEdgesPair[9][i]=binEdgesPair[8][i];*/
1153
1154   fPairQA = new THnSparseF("pairQA", "QA for Pair; invmass[GeV/c^2]; KF chi2; opening angle; signed Lxy; KF ip; pair code", nDimPair, nBinPair);
1155   //fPairQA = new THnSparseF("pairQA", "QA for Pair; invmass[GeV/c^2]; KF chi2; opening angle; signed Lxy; KF ip; pair code; dca sig trk1; dca sig trk2; pt trk1; pt trk2 ", nDimPair, nBinPair);
1156   for(Int_t idim = 0; idim < nDimPair; idim++){
1157     fPairQA->SetBinEdges(idim, binEdgesPair[idim]);
1158   }
1159
1160   fSecVtxList->AddAt(fPairQA,0);
1161
1162   const Int_t nDimSecvtx=6;
1163   Double_t* binEdgesSecvtx[nDimSecvtx];
1164   Int_t nBinSecvtx[nDimSecvtx] = {200, 500, 2000, 2000, 11, 3};
1165   const Double_t kNtrksmin = 0, kNtrksmax= 3;
1166   for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
1167     binEdgesSecvtx[ivar] = new Double_t[nBinSecvtx[ivar] + 1];
1168
1169   for(Int_t i=0; i<=nBinSecvtx[0]; i++) binEdgesSecvtx[0][i]=binEdgesPair[0][i];
1170   for(Int_t i=0; i<=nBinSecvtx[1]; i++) binEdgesSecvtx[1][i]=binEdgesPair[1][i];
1171   for(Int_t i=0; i<=nBinSecvtx[2]; i++) binEdgesSecvtx[2][i]=binEdgesPair[3][i];
1172   for(Int_t i=0; i<=nBinSecvtx[3]; i++) binEdgesSecvtx[3][i]=binEdgesPair[4][i];
1173   for(Int_t i=0; i<=nBinSecvtx[4]; i++) binEdgesSecvtx[4][i]=binEdgesPair[5][i];
1174   for(Int_t i=0; i<=nBinSecvtx[5]; i++) binEdgesSecvtx[5][i]=(Double_t)kNtrksmin + (kNtrksmax - kNtrksmin)/nBinSecvtx[5]*(Double_t)i;
1175
1176   fSecvtxQA = new THnSparseF("secvtxQA", "QA for Secvtx; invmass[GeV/c^2]; KF chi2; signed Lxy; KF ip; pair code; n tracks ", nDimSecvtx, nBinSecvtx);
1177   for(Int_t idim = 0; idim < nDimSecvtx; idim++){
1178     fSecvtxQA->SetBinEdges(idim, binEdgesSecvtx[idim]);
1179   }
1180
1181   fSecVtxList->AddAt(fSecvtxQA,1);
1182   for(Int_t ivar = 0; ivar < nDimPair; ivar++)
1183     delete binEdgesPair[ivar];
1184   for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
1185     delete binEdgesSecvtx[ivar];
1186 }