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