Since it contains fixes of coding rule violations, all classes are involved. Further...
[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  *                                                                        *
17  * Secondary vertexing construction Class                                 *
18  *  Construct secondary vertex from Beauty hadron with electron and       *
19  *  hadrons, then apply selection criteria                                *
20  *                                                                        *
21  * Authors:                                                               *
22  *   MinJung Kweon <minjung@physi.uni-heidelberg.de>                      *
23  *                                                                        *
24  **************************************************************************/
25
26 #include <iostream>
27 #include <TH2F.h>
28 #include <TH3F.h>
29 #include <TLorentzVector.h>
30 #include <TNtuple.h>
31 #include <TParticle.h>
32
33 #include <AliESDEvent.h>
34 #include <AliESDtrack.h>
35 #include "AliHFEsecVtx.h"
36 #include <AliKFParticle.h>
37 #include <AliKFVertex.h>
38 #include <AliLog.h>
39 #include <AliPID.h>
40 #include <AliStack.h>
41
42 ClassImp(AliHFEsecVtx)
43
44 //_______________________________________________________________________________________________
45 AliHFEsecVtx::AliHFEsecVtx():
46          fESD1(0x0)
47         ,fStack(0x0)
48         ,fNparents(0)
49         ,fHistTagged()
50         ,fPairTagged(0)
51         ,fdistTwoSecVtx(-1)
52         ,fcosPhi(-1)
53         ,fsignedLxy(-1)
54         ,finvmass(-1)
55         ,finvmassSigma(-1)
56         ,fBTagged(0)
57         ,fBElec(0)
58
59         // Default constructor
60
61         Init();
62
63 }
64
65 //_______________________________________________________________________________________________
66 AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
67          TObject(p)
68         ,fESD1(0x0)
69         ,fStack(0x0)
70         ,fNparents(p.fNparents)
71         ,fHistTagged()
72         ,fPairTagged(p.fPairTagged)
73         ,fdistTwoSecVtx(p.fdistTwoSecVtx)
74         ,fcosPhi(p.fcosPhi)
75         ,fsignedLxy(p.fsignedLxy)
76         ,finvmass(p.finvmass)
77         ,finvmassSigma(p.finvmassSigma)
78         ,fBTagged(p.fBTagged)
79         ,fBElec(p.fBElec)
80
81         // Copy constructor
82 }
83
84 //_______________________________________________________________________________________________
85 AliHFEsecVtx&
86 AliHFEsecVtx::operator=(const AliHFEsecVtx &)
87 {
88   // Assignment operator
89
90   AliInfo("Not yet implemented.");
91   return *this;
92 }
93
94 //_______________________________________________________________________________________________
95 AliHFEsecVtx::~AliHFEsecVtx()
96 {
97         // Destructor
98
99         //cout << "Analysis Done." << endl;
100 }
101
102 //__________________________________________
103 void AliHFEsecVtx::Init()
104 {
105
106         // set pdg code and index
107
108         fNparents = 7;
109
110         fParentSelect[0][0] =  411;
111         fParentSelect[0][1] =  421;
112         fParentSelect[0][2] =  431;
113         fParentSelect[0][3] = 4122;
114         fParentSelect[0][4] = 4132;
115         fParentSelect[0][5] = 4232;
116         fParentSelect[0][6] = 4332;
117
118         fParentSelect[1][0] =  511;
119         fParentSelect[1][1] =  521;
120         fParentSelect[1][2] =  531;
121         fParentSelect[1][3] = 5122;
122         fParentSelect[1][4] = 5132;
123         fParentSelect[1][5] = 5232;
124         fParentSelect[1][6] = 5332;
125
126 /*
127         fid[0][0] = 0;
128         fid[0][1] = 1;
129         fid[0][2] = 2;
130
131         fid[1][0] = 0;
132         fid[1][1] = 1;
133         fid[1][2] = 3;
134
135         fid[2][0] = 0;
136         fid[2][1] = 2;
137         fid[2][2] = 3;
138
139         fid[3][0] = 1;
140         fid[3][1] = 2;
141         fid[3][2] = 3;
142
143         fia[0][0][0] = 0;
144         fia[0][0][1] = 1;
145         fia[0][1][0] = 0;
146         fia[0][1][1] = 2;
147         fia[0][2][0] = 1;
148         fia[0][2][1] = 2;
149
150         fia[1][0][0] = 0;
151         fia[1][0][1] = 1;
152         fia[1][1][0] = 0;
153         fia[1][1][1] = 3;
154         fia[1][2][0] = 1;
155         fia[1][2][1] = 3;
156
157         fia[2][0][0] = 0;
158         fia[2][0][1] = 2;
159         fia[2][1][0] = 0;
160         fia[2][1][1] = 3;
161         fia[2][2][0] = 2;
162         fia[2][2][1] = 3;
163
164         fia[3][0][0] = 1;
165         fia[3][0][1] = 2;
166         fia[3][1][0] = 1;
167         fia[3][1][1] = 3;
168         fia[3][2][0] = 2;
169         fia[3][2][1] = 3;
170 */
171
172
173
174 //__________________________________________
175 void AliHFEsecVtx::ResetTagVar()
176 {
177         // reset tag variables
178
179         fdistTwoSecVtx = -1;
180         fcosPhi = -1;
181         fsignedLxy = -1;
182         finvmass = -1;
183         finvmassSigma = -1;
184         fBTagged = kFALSE;
185         fBElec = kFALSE;
186 }
187
188 //__________________________________________
189 void AliHFEsecVtx::InitAnaPair()
190 {
191         // initialize pair tagging variables
192
193         fPairTagged = 0;
194         for (Int_t i=0; i<20; i++){
195                 fpairedTrackID[i] = -1;
196                 fpairedChi2[i] = -1;
197                 fpairedInvMass[i] = -1;
198                 fpairedSignedLxy[i] = -1;
199         }
200
201         fid[0][0] = 0;
202         fid[0][1] = 1;
203         fid[0][2] = 2;
204
205         fid[1][0] = 0;
206         fid[1][1] = 1;
207         fid[1][2] = 3;
208
209         fid[2][0] = 0;
210         fid[2][1] = 2;
211         fid[2][2] = 3;
212
213         fid[3][0] = 1;
214         fid[3][1] = 2;
215         fid[3][2] = 3;
216
217         fia[0][0][0] = 0;
218         fia[0][0][1] = 1;
219         fia[0][1][0] = 0;
220         fia[0][1][1] = 2;
221         fia[0][2][0] = 1;
222         fia[0][2][1] = 2;
223
224         fia[1][0][0] = 0;
225         fia[1][0][1] = 1;
226         fia[1][1][0] = 0;
227         fia[1][1][1] = 3;
228         fia[1][2][0] = 1;
229         fia[1][2][1] = 3;
230
231         fia[2][0][0] = 0;
232         fia[2][0][1] = 2;
233         fia[2][1][0] = 0;
234         fia[2][1][1] = 3;
235         fia[2][2][0] = 2;
236         fia[2][2][1] = 3;
237
238         fia[3][0][0] = 1;
239         fia[3][0][1] = 2;
240         fia[3][1][0] = 1;
241         fia[3][1][1] = 3;
242         fia[3][2][0] = 2;
243         fia[3][2][1] = 3;
244
245 }
246
247 //_______________________________________________________________________________________________
248 void AliHFEsecVtx::CreateHistograms(TString hnopt)
249
250         // create histograms
251
252         fkSourceLabel[kAll]="all";
253         fkSourceLabel[kDirectCharm]="directCharm";
254         fkSourceLabel[kDirectBeauty]="directBeauty";
255         fkSourceLabel[kBeautyCharm]="beauty2charm";
256         fkSourceLabel[kGamma]="gamma";
257         fkSourceLabel[kPi0]="pi0";
258         fkSourceLabel[kElse]="others";
259         fkSourceLabel[kBeautyGamma]="beauty22gamma";
260         fkSourceLabel[kBeautyPi0]="beauty22pi0";
261         fkSourceLabel[kBeautyElse]="beauty22others";
262
263
264         TString hname;
265         for (Int_t isource = 0; isource < 10; isource++ ){
266
267            hname=hnopt+"InvMass_"+fkSourceLabel[isource];
268            fHistPair[isource].fInvMass = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
269            hname=hnopt+"InvMassCut1_"+fkSourceLabel[isource];
270            fHistPair[isource].fInvMassCut1 = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
271            hname=hnopt+"InvMassCut2_"+fkSourceLabel[isource];
272            fHistPair[isource].fInvMassCut2 = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
273            hname=hnopt+"KFChi2_"+fkSourceLabel[isource];
274            fHistPair[isource].fKFChi2 = new TH1F(hname,hname,200,0,20);
275            hname=hnopt+"OpenAngle_"+fkSourceLabel[isource];
276            fHistPair[isource].fOpenAngle = new TH1F(hname,hname,100,0,3.14);
277            hname=hnopt+"CosOpenAngle_"+fkSourceLabel[isource];
278            fHistPair[isource].fCosOpenAngle = new TH1F(hname,hname,100,-1.1,1.1);
279            hname=hnopt+"SignedLxy_"+fkSourceLabel[isource];
280            fHistPair[isource].fSignedLxy = new TH2F(hname,hname,1000,-5,5,120,-2,10);
281            hname=hnopt+"KFIP_"+fkSourceLabel[isource];
282            fHistPair[isource].fKFIP = new TH1F(hname,hname,1000,-5,5);
283            hname=hnopt+"IPMax_"+fkSourceLabel[isource];
284            fHistPair[isource].fIPMax= new TH1F(hname,hname,500,0,5);
285
286         }
287
288         hname=hnopt+"pt_beautyelec";
289         fHistTagged.fPtBeautyElec= new TH1F(hname,hname,150,0,30);
290         hname=hnopt+"pt_taggedelec";
291         fHistTagged.fPtTaggedElec= new TH1F(hname,hname,150,0,30);
292         hname=hnopt+"pt_righttaggedelec";
293         fHistTagged.fPtRightTaggedElec = new TH1F(hname,hname,150,0,30);
294         hname=hnopt+"pt_wrongtaggedelec";
295         fHistTagged.fPtWrongTaggedElec = new TH1F(hname,hname,150,0,30);
296
297         hname=hnopt+"InvmassBeautyElecSecVtx";
298         fHistTagged.fInvmassBeautyElecSecVtx= new TH1F(hname,hname,120,-2,10);
299         hname=hnopt+"InvmassNotBeautyElecSecVtx";
300         fHistTagged.fInvmassNotBeautyElecSecVtx= new TH1F(hname,hname,120,-2,10);
301
302         hname=hnopt+"SignedLxyBeautyElecSecVtx";
303         fHistTagged.fSignedLxyBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
304         hname=hnopt+"SignedLxyNotBeautyElecSecVtx";
305         fHistTagged.fSignedLxyNotBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
306
307         hname=hnopt+"DistTwoVtxBeautyElecSecVtx";
308         fHistTagged.fDistTwoVtxBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
309         hname=hnopt+"DistTwoVtxNotBeautyElecSecVtx";
310         fHistTagged.fDistTwoVtxNotBeautyElecSecVtx= new TH1F(hname,hname,1000,-5,5);
311
312         hname=hnopt+"CosPhiBeautyElecSecVtx";
313         fHistTagged.fCosPhiBeautyElecSecVtx= new TH1F(hname,hname,100,-1.1,1.1);
314         hname=hnopt+"CosPhiNotBeautyElecSecVtx";
315         fHistTagged.fCosPhiNotBeautyElecSecVtx= new TH1F(hname,hname,100,-1.1,1.1);
316
317         hname=hnopt+"Chi2BeautyElecSecVtx";
318         fHistTagged.fChi2BeautyElecSecVtx= new TH1F(hname,hname,200,0,20);
319         hname=hnopt+"Chi2NotBeautyElecSecVtx";
320         fHistTagged.fChi2NotBeautyElecSecVtx= new TH1F(hname,hname,200,0,20);
321
322         hname=hnopt+"InvmassBeautyElec2trkVtx";
323         fHistTagged.fInvmassBeautyElec2trkVtx= new TH1F(hname,hname,120,-2,10);
324         hname=hnopt+"InvmassNotBeautyElec2trkVtx";
325         fHistTagged.fInvmassNotBeautyElec2trkVtx= new TH1F(hname,hname,120,-2,10);
326
327         hname=hnopt+"SignedLxyBeautyElec2trkVtx";
328         fHistTagged.fSignedLxyBeautyElec2trkVtx= new TH1F(hname,hname,1000,-5,5);
329         hname=hnopt+"SignedLxyNotBeautyElec2trkVtx";
330         fHistTagged.fSignedLxyNotBeautyElec2trkVtx= new TH1F(hname,hname,1000,-5,5);
331
332         hname=hnopt+"Chi2BeautyElec2trkVtx";
333         fHistTagged.fChi2BeautyElec2trkVtx= new TH1F(hname,hname,200,0,20);
334         hname=hnopt+"Chi2NotBeautyElec2trkVtx";
335         fHistTagged.fChi2NotBeautyElec2trkVtx= new TH1F(hname,hname,200,0,20);
336 }
337
338 //_______________________________________________________________________________________________
339 void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t index2)
340 {
341         // calculate e-h pair characteristics and tag pair 
342
343         Int_t sourcePart = PairCode(track1,track2);
344
345         // get KF particle input pid
346         Int_t pdg1 = GetMCPID(track1);
347         Int_t pdg2 = GetMCPID(track2);
348         
349
350         // create KF particle of pair
351         AliKFParticle::SetField(fESD1->GetMagneticField());
352         AliKFParticle kfTrack1(*track1, pdg1);
353         AliKFParticle kfTrack2(*track2, pdg2);
354
355         AliKFParticle kfSecondary(kfTrack1,kfTrack2);
356
357         // copy primary vertex from ESD
358         AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
359         if( primVtxCopy.GetNDF() <1 ) return;
360
361         //primary vertex point
362         Double_t pvx = primVtxCopy.GetX();
363         Double_t pvy = primVtxCopy.GetY();
364         //Double_t pvz = primVtxCopy.GetZ();
365
366         //secondary vertex point from kf particle
367         Double_t kfx = kfSecondary.GetX();
368         Double_t kfy = kfSecondary.GetY();
369         //Double_t kfz = kfSecondary.GetZ();
370         
371         //momentum at the decay point from kf particle
372         Double_t kfpx = kfSecondary.GetPx();
373         Double_t kfpy = kfSecondary.GetPy();
374         //Double_t kfpz = kfSecondary.GetPz();
375         
376
377         Double_t dx = kfx-pvx;
378         Double_t dy = kfy-pvy;
379
380
381
382         // discriminating variables ----------------------------------------------------------
383
384         // invariant mass of the KF particle
385         Double_t invmass = -1;
386         Double_t invmassSigma = -1;
387         kfSecondary.GetMass(invmass,invmassSigma);
388
389         // chi2 of the KF particle
390         Double_t kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
391
392         // opening angle between two particles in XY plane
393         Double_t phi = kfTrack1.GetAngleXY(kfTrack2);
394         Double_t cosphi = TMath::Cos(phi);
395
396         // projection of kf vertex vector to the kf momentum direction 
397         Double_t costheta = ( dx*kfpx + dy*kfpy)/TMath::Sqrt(dx*dx+dy*dy)*TMath::Sqrt(kfpx*kfpx + kfpy*kfpy);
398         Double_t signedLxy = TMath::Sqrt(dx*dx+dy*dy)*costheta;
399
400         // DCA from primary to e-h KF particle (impact parameter of KF particle)
401         Double_t vtx[2]={pvx, pvy}; 
402         Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
403
404
405               Float_t dcaR=-1; 
406         Float_t dcaR1=-1, dcaR2=-1;
407         Float_t dcaZ1=-1, dcaZ2=-1;
408         track1->GetImpactParameters(dcaR1,dcaZ1);
409         track2->GetImpactParameters(dcaR2,dcaZ2);
410
411               if (TMath::Abs(dcaR1) >= TMath::Abs(dcaR2)) dcaR=dcaR1;
412               else dcaR=dcaR2;
413
414         // fill histograms 
415         fHistPair[sourcePart].fInvMass->Fill(invmass,invmassSigma);
416         fHistPair[sourcePart].fKFChi2->Fill(kfchi2);
417         fHistPair[sourcePart].fOpenAngle->Fill(phi);
418         fHistPair[sourcePart].fCosOpenAngle->Fill(cosphi);
419         fHistPair[sourcePart].fSignedLxy->Fill(signedLxy,invmass);
420         fHistPair[sourcePart].fKFIP->Fill(kfip);
421         fHistPair[sourcePart].fIPMax->Fill(TMath::Abs(dcaR));
422
423         // pair cuts 
424         if( kfchi2 >2. ) return;
425
426         if ( signedLxy > 0.05 && cosphi > 0.5) fHistPair[sourcePart].fInvMassCut1->Fill(invmass,invmassSigma);
427
428         // pair tagging if it passed the above cuts
429         if(signedLxy > 0. && cosphi > 0.) fHistPair[sourcePart].fInvMassCut2->Fill(invmass,invmassSigma);
430
431
432         // pair tagging condition
433         if ( signedLxy > 0.0 && cosphi > 0) { // testing loose cut
434         //if ( signedLxy > 0.06 && cosphi > 0) {
435           fpairedTrackID[fPairTagged] = index2;
436           fpairedChi2[fPairTagged] = kfchi2;
437           fpairedInvMass[fPairTagged] = invmass;
438           fpairedSignedLxy[fPairTagged] = signedLxy;
439           fPairTagged++;
440         }
441
442 }
443
444 //_______________________________________________________________________________________________
445 void AliHFEsecVtx::RunSECVTX(AliESDtrack *track)
446 {
447         // run secondary vertexing algorithm and do tagging
448
449         ResetTagVar();
450
451         Int_t imclabel = TMath::Abs(track->GetLabel());
452         if(imclabel<0) return;
453         TParticle* mcpart = fStack->Particle(imclabel);
454         Int_t esource = GetElectronSource(imclabel);
455         if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){
456                 fHistTagged.fPtBeautyElec->Fill(mcpart->Pt());
457                 fBElec = kTRUE;
458         }
459
460
461         if (fPairTagged >= 4) {
462           FindSECVTXCandid4Tracks(track);         
463         }
464         else if (fPairTagged == 3) {
465           FindSECVTXCandid3Tracks(track);         
466         }
467         else if (fPairTagged == 2) {
468           FindSECVTXCandid2Tracks(track);         
469         }
470         else if (fPairTagged == 1) {
471           ApplyPairTagCut();      
472         }
473
474
475         if (fBTagged) {
476                 fHistTagged.fPtTaggedElec->Fill(mcpart->Pt());
477                 if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){
478                         fHistTagged.fPtRightTaggedElec->Fill(mcpart->Pt());
479                 }
480                 else fHistTagged.fPtWrongTaggedElec->Fill(mcpart->Pt());
481         }
482
483 }
484
485 //_______________________________________________________________________________________________
486 void AliHFEsecVtx::ApplyPairTagCut()
487 {
488         // apply tagging cut for e-h pair
489
490         if (fBElec){
491                 fHistTagged.fInvmassBeautyElec2trkVtx->Fill(fpairedInvMass[0]);
492                 fHistTagged.fSignedLxyBeautyElec2trkVtx->Fill(fpairedSignedLxy[0]);
493                 fHistTagged.fChi2BeautyElec2trkVtx->Fill(fpairedChi2[0]);
494         }
495         else {
496                 fHistTagged.fInvmassNotBeautyElec2trkVtx->Fill(fpairedInvMass[0]);
497                 fHistTagged.fSignedLxyNotBeautyElec2trkVtx->Fill(fpairedSignedLxy[0]);
498                 fHistTagged.fChi2NotBeautyElec2trkVtx->Fill(fpairedChi2[0]);
499         }
500                 
501         if (fpairedChi2[0] > 2.0) return;
502         if (fpairedInvMass[0] < 1.5) return;
503         if (fpairedSignedLxy[0] < 0.05) return;
504
505         fBTagged = kTRUE;
506 }
507
508 //_______________________________________________________________________________________________
509 Bool_t AliHFEsecVtx::ApplyPairTagCut(Int_t id) 
510 {
511         // apply tagging cut for e-h pair of given indexed hadron
512
513         if (fpairedChi2[id] > 2.0) return kFALSE;
514         if (fpairedInvMass[id] < 1.5) return kFALSE;
515         if (fpairedSignedLxy[id] < 0.05) return kFALSE;
516
517         fBTagged = kTRUE;
518         return kTRUE;
519                 
520 }
521
522 //_______________________________________________________________________________________________
523 Bool_t AliHFEsecVtx::ApplyTagCut(Double_t chi2) 
524 {
525         // apply tagging cut for secondary vertex
526
527         if (chi2 > 2.0) return kFALSE;
528         if (finvmass < 1.5) return kFALSE;
529         if (fsignedLxy < 0.05) return kFALSE;
530         if (fcosPhi < 0.90) return kFALSE;
531         if (fdistTwoSecVtx > 0.1) return kFALSE;
532         
533         fBTagged = kTRUE;
534         return kTRUE;
535 }
536
537 //_______________________________________________________________________________________________
538 void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2, Int_t id1, Int_t id2)
539 {
540         // apply tagging cut for e-h pair of given indexed hadron
541
542         if(!ApplyTagCut(chi2)){
543                 if(!ApplyPairTagCut(id1)) ApplyPairTagCut(id2);
544         }
545         
546 }
547
548 //_______________________________________________________________________________________________
549 void AliHFEsecVtx::FindSECVTXCandid4Tracks(AliESDtrack *track) 
550 {
551         // find secondary vertex for >= 4 e-h pairs
552
553         // sort pair in increasing order (kFALSE-increasing order)
554         Int_t index[20];
555         Int_t indexA[4];
556         Int_t indexB[3];
557         Double_t sevchi2[4];
558         AliESDtrack *htrack[4];
559
560         if(fPairTagged>20) return; // protection
561
562         TMath::Sort(fPairTagged,fpairedChi2,index,kFALSE);
563
564         // select 4 partner tracks retruning smallest pair chi2 
565
566         for (Int_t i=0; i<4; i++){
567                 htrack[i] = fESD1->GetTrack(fpairedTrackID[index[i]]);
568         }
569         
570         // calculated chi2 with 1 electron and 3 partner tracks 
571         for (Int_t i=0; i<4; i++){
572                 sevchi2[i] = GetSecVtxChi2(track, htrack[fid[i][0]], htrack[fid[i][1]], htrack[fid[i][2]]);
573         }
574
575         // select 3 partner tracks retruning smallest pair chi2
576         // [think] if two smallest chi2 are similar, have to think about better handling of selection
577         TMath::Sort(4,sevchi2,indexA,kFALSE);
578
579         // calculated chi2 with 1 electron and 2 partner tracks 
580         for (Int_t i=0; i<3; i++){
581                 sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
582         }
583
584         // select 2 partner tracks retruning smallest pair chi2
585         TMath::Sort(3,sevchi2,indexB,kFALSE);
586
587         // calculate secondary vertex quality variables with 1 electron and 2 hadrons
588         CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
589
590         if (fBElec){
591                 fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
592         }
593         else {
594                 fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
595         }
596
597         ApplyVtxTagCut(sevchi2[indexB[0]],index[fia[indexA[0]][indexB[0]][0]],index[fia[indexA[0]][indexB[0]][1]]);
598
599 }
600
601 //_______________________________________________________________________________________________
602 void AliHFEsecVtx::FindSECVTXCandid3Tracks(AliESDtrack *track) 
603 {
604         // find secondary vertex for 3 e-h pairs
605
606         // sort pair in increasing order (kFALSE-increasing order)
607         Int_t indexA[1] = { 0 };
608         Int_t indexB[3];
609         Double_t sevchi2[3];
610         AliESDtrack *htrack[3];
611
612         // select 4 partner tracks retruning smallest pair chi2 
613
614         for (Int_t i=0; i<3; i++){
615                 htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
616         }
617         
618         // calculated chi2 with 1 electron and 2 partner tracks 
619         for (Int_t i=0; i<3; i++){
620                 sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
621         }
622
623         // select 2 partner tracks retruning smallest pair chi2
624         TMath::Sort(3,sevchi2,indexB,kFALSE);
625
626         // calculate secondary vertex quality variables with 1 electron and 2 hadrons
627         CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
628
629         if (fBElec){
630                 fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
631         }
632         else {
633                 fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
634         }
635
636         ApplyVtxTagCut(sevchi2[indexB[0]],fia[indexA[0]][indexB[0]][0],fia[indexA[0]][indexB[0]][1]);
637 }
638
639 //_______________________________________________________________________________________________
640 void AliHFEsecVtx::FindSECVTXCandid2Tracks(AliESDtrack *track) 
641 {
642         // find secondary vertex for 2 e-h pairs
643
644         Double_t sevchi2[1];
645         AliESDtrack *htrack[2];
646
647         for (Int_t i=0; i<2; i++){
648                 htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
649         }
650         
651         sevchi2[0] = GetSecVtxChi2(track, htrack[0], htrack[1]);
652
653         // calculate secondary vertex quality variables with 1 electron and 2 hadrons
654         CalcSECVTXProperty(track,htrack[0],htrack[1]);
655
656         if (fBElec){
657                 fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[0]);
658         }
659         else {
660                 fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[0]);
661         }
662
663         ApplyVtxTagCut(sevchi2[0],0,1);
664 }
665
666 //_______________________________________________________________________________________________
667 void AliHFEsecVtx::CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
668 {
669         // calculate secondary vertex properties
670
671         Int_t pdg1 = GetMCPID(track1);
672         Int_t pdg2 = GetMCPID(track2);
673         Int_t pdg3 = GetMCPID(track3);
674
675         AliKFParticle::SetField(fESD1->GetMagneticField());
676         AliKFParticle kfTrack1(*track1, pdg1);
677         AliKFParticle kfTrack2(*track2, pdg2);
678         AliKFParticle kfTrack3(*track3, pdg3);
679
680         AliKFParticle kfSecondary12(kfTrack1,kfTrack2);
681         AliKFParticle kfSecondary13(kfTrack1,kfTrack3);
682         AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
683
684         // copy primary vertex from ESD
685         AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));        
686         //printf("primary ndf= %d\n",primVtxCopy.GetNDF());
687         if( primVtxCopy.GetNDF() <1 ) return;
688
689         Double_t kdx12 = kfSecondary12.GetX()-primVtxCopy.GetX();
690         Double_t kdy12 = kfSecondary12.GetY()-primVtxCopy.GetY();
691         //Double_t kdz12 = kfSecondary12.GetZ()-primVtxCopy.GetZ();
692
693         Double_t kdx13 = kfSecondary13.GetX()-primVtxCopy.GetX();
694         Double_t kdy13 = kfSecondary13.GetY()-primVtxCopy.GetY();
695         //Double_t kdz13 = kfSecondary13.GetZ()-primVtxCopy.GetZ();
696
697         Double_t kdx = kfSecondary.GetX()-primVtxCopy.GetX();
698         Double_t kdy = kfSecondary.GetY()-primVtxCopy.GetY();
699         //Double_t kdz = kfSecondary.GetZ()-primVtxCopy.GetZ();
700
701         // calculate distance and angle between two secvtxes 
702         fdistTwoSecVtx = TMath::Sqrt((kdx12-kdx13)*(kdx12-kdx13) + (kdy12-kdy13)*(kdy12-kdy13));
703         fcosPhi = ( kdx12*kdx13 + kdy12*kdy13 ) / ( TMath::Sqrt(kdx12*kdx12+kdy12*kdy12)*TMath::Sqrt(kdx13*kdx13+kdy13*kdy13) );
704         //Double_t lengthdiff = TMath::Abs(TMath::Sqrt(kdx12*kdx12+kdy12*kdy12) - TMath::Sqrt(kdx13*kdx13+kdy13*kdy13));
705
706         // calculate angle between secondary vertex vector and secondary particle momentum vector in transverse plane
707         Double_t cosTheta = ( kdx*kfSecondary.GetPx() + kdy*kfSecondary.GetPy()) / TMath::Sqrt(kdx*kdx+kdy*kdy)*TMath::Sqrt(kfSecondary.GetPx()*kfSecondary.GetPx()+kfSecondary.GetPy()*kfSecondary.GetPy());
708         // calculate signed Lxy
709         fsignedLxy = TMath::Sqrt(kdx*kdx+kdy*kdy)*cosTheta;
710
711         // calculate invariant mass of the kf particle
712         kfSecondary.GetMass(finvmass,finvmassSigma);
713
714         if (fBElec){
715                 fHistTagged.fInvmassBeautyElecSecVtx->Fill(finvmass);
716                 fHistTagged.fSignedLxyBeautyElecSecVtx->Fill(fsignedLxy);
717                 fHistTagged.fDistTwoVtxBeautyElecSecVtx->Fill(fdistTwoSecVtx);
718                 fHistTagged.fCosPhiBeautyElecSecVtx->Fill(fcosPhi);
719         }
720         else {
721                 fHistTagged.fInvmassNotBeautyElecSecVtx->Fill(finvmass);
722                 fHistTagged.fSignedLxyNotBeautyElecSecVtx->Fill(fsignedLxy);
723                 fHistTagged.fDistTwoVtxNotBeautyElecSecVtx->Fill(fdistTwoSecVtx);
724                 fHistTagged.fCosPhiNotBeautyElecSecVtx->Fill(fcosPhi);
725         }
726
727 }
728
729 //_______________________________________________________________________________________________
730 Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track) 
731 {
732         // return mc pid
733
734         Int_t label = TMath::Abs(track->GetLabel());
735         TParticle* mcpart = fStack->Particle(label);
736         if ( !mcpart ) return 0;
737         Int_t pdgCode = mcpart->GetPdgCode();
738
739         return pdgCode;
740 }
741
742 //_______________________________________________________________________________________________
743 Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3, AliESDtrack* track4)
744 {
745         // return 4 track secondary vertex chi2
746
747         Int_t pdg1 = GetMCPID(track1);
748         Int_t pdg2 = GetMCPID(track2);
749         Int_t pdg3 = GetMCPID(track3);
750         Int_t pdg4 = GetMCPID(track4);
751
752         AliKFParticle::SetField(fESD1->GetMagneticField());
753         AliKFParticle kfTrack1(*track1, pdg1);
754         AliKFParticle kfTrack2(*track2, pdg2);
755         AliKFParticle kfTrack3(*track3, pdg3);
756         AliKFParticle kfTrack4(*track4, pdg4);
757         AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3,kfTrack4);
758
759         return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
760
761 }
762
763 //_______________________________________________________________________________________________
764 Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
765 {
766         // return 3 track secondary vertex chi2
767
768         Int_t pdg1 = GetMCPID(track1);
769         Int_t pdg2 = GetMCPID(track2);
770         Int_t pdg3 = GetMCPID(track3);
771
772         AliKFParticle::SetField(fESD1->GetMagneticField());
773         AliKFParticle kfTrack1(*track1, pdg1);
774         AliKFParticle kfTrack2(*track2, pdg2);
775         AliKFParticle kfTrack3(*track3, pdg3);
776         AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
777
778         return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
779
780 }
781
782 //_______________________________________________________________________________________________
783 Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2)
784 {
785         // return 2 track secondary vertex chi2
786
787         Int_t pdg1 = GetMCPID(track1);
788         Int_t pdg2 = GetMCPID(track2);
789
790         AliKFParticle::SetField(fESD1->GetMagneticField());
791         AliKFParticle kfTrack1(*track1, pdg1);
792         AliKFParticle kfTrack2(*track2, pdg2);
793         AliKFParticle kfSecondary(kfTrack1,kfTrack2);
794
795         return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
796
797 }
798
799 //_______________________________________________________________________________________________
800 Int_t AliHFEsecVtx::PairOrigin(AliESDtrack* track1, AliESDtrack* track2)
801 {
802
803         //
804         // return pdg code of the origin(source) of the pair 
805         // 
806         //
807         // ---*---*---*-----ancester A----- track1
808         //                        |____*______ 
809         //                             |______ track2
810         // => if they originated from same ancester, 
811         //    then return "the absolute value of pdg code of ancester A"
812         //
813         // ---*---*---B hadron-----ancester A----- track1
814         //                               |____*______ 
815         //                                    |______ track2
816         // => if they originated from same ancester, and this ancester originally comes from B hadrons
817         //    then return -1*"the absolute value of pdg code of ancester A"
818         //
819         // caution : it can also return parton pdg code if it originated from same string or gluon spliting. 
820         //           
821
822         TParticle* part1 = fStack->Particle(TMath::Abs(track1->GetLabel()));
823         TParticle* part2 = fStack->Particle(TMath::Abs(track2->GetLabel()));
824         if (!(part1)) return 0;
825         if (!(part2)) return 0;
826
827         TParticle* part1Crtgen = part1; // copy track into current generation particle
828         TParticle* part2Crtgen = part2; // copy track into current generation particle
829
830
831         Int_t sourcePDG = 0;
832
833         // if the two tracks' mother's label is same, get the mother info
834         // in case of charm, check if it originated from beauty
835         for (Int_t i=0; i<100; i++){ // iterate 100
836                 Int_t iLabel = part1Crtgen->GetFirstMother(); //label of mother of current generation for 1st partilce
837                 if (iLabel < 0) return 0;
838
839                 for (Int_t j=0; j<100; j++){ // iterate 100
840                         Int_t jLabel = part2Crtgen->GetFirstMother(); //label of mother of current generation for 2nd partilce
841                         if (jLabel < 0) return 0; // if jLabel == -1
842
843                         if (iLabel == jLabel){ // check if two tracks are originated from same mother
844                                 TParticle* thismother = fStack->Particle(jLabel); // if yes, get "thismother" info 
845                                 sourcePDG = abs(thismother->GetPdgCode()); // get the pdg code of "this mother"
846
847                                 // check ancester to see if it is originally from beauty 
848                                 for (Int_t k=0; k<10; k++){ // check up to 10 ancesters
849                                         Int_t ancesterLabel = thismother->GetFirstMother();
850                                         if (ancesterLabel < 0) return sourcePDG; // if "thismoter" doesn't have mother anymore, return thismother's pdg  
851
852                                         TParticle* thisancester = fStack->Particle(ancesterLabel);
853                                         Int_t ancesterPDG = abs(thisancester->GetPdgCode());
854
855                                         for (Int_t l=0; l<fNparents; l++){
856                                                 if (abs(ancesterPDG)==fParentSelect[1][l]){
857                                                         sourcePDG = -1*sourcePDG; // multiply -1 for charm from bottom
858                                                         return sourcePDG;
859                                                 }
860                                         }
861                                         thismother = thisancester;
862                                 }
863
864                         }
865                         part2Crtgen = fStack->Particle(jLabel); // if their mother is different, go up to previous generation of 2nd particle
866                 }
867                 part1Crtgen = fStack->Particle(iLabel); // if their mother is different, go up to previous generation of 2nd particle
868
869                 // if you don't find additionional(2nd particle) track originated from a given beauty hadron, break to save time
870                 Int_t motherPDGtmp = abs(part1Crtgen->GetPdgCode());
871                 for (Int_t l=0; l<fNparents; l++){
872                         if (abs(motherPDGtmp)==fParentSelect[1][l]) return sourcePDG;
873                 }
874
875         }
876
877         return sourcePDG;
878
879 }
880
881 //_______________________________________________________________________________________________
882 Int_t AliHFEsecVtx::PairCode(AliESDtrack* track1,AliESDtrack* track2)
883 {
884
885         //           
886         // return pair code which is predefinded as:
887         //  kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
888         //           
889
890         Int_t pairOriginsPDG = PairOrigin(track1,track2);
891
892         Int_t sourcePart = kElse;
893
894         if (pairOriginsPDG < 0) {
895                 sourcePart = kBeautyElse;
896         }
897         for (Int_t i=0; i<fNparents; i++){
898                 if (abs(pairOriginsPDG)==fParentSelect[0][i]) {
899                         if (pairOriginsPDG>0) sourcePart = kDirectCharm;
900                         if (pairOriginsPDG<0) {
901                                 sourcePart = kBeautyCharm;
902                         }
903                 }
904                 if (abs(pairOriginsPDG)==fParentSelect[1][i]) {
905                         if (pairOriginsPDG>0) {
906                                 sourcePart = kDirectBeauty;
907                         }
908                         if (pairOriginsPDG<0)  return kElse;
909                 }
910         }
911         if (pairOriginsPDG == 22) sourcePart = kGamma;
912         if (pairOriginsPDG == -22) {
913                 sourcePart = kBeautyGamma;
914         }
915         if (pairOriginsPDG == 111) sourcePart = kPi0;
916         if (pairOriginsPDG == -111) {
917                 sourcePart = kBeautyPi0;
918         }
919
920         return sourcePart;
921
922 }
923
924 //_______________________________________________________________________________________________
925 Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack) 
926 {
927
928     // return decay electron's origin 
929
930     if (iTrack < 0) {
931       AliDebug(1, "Stack label is negative, return\n");
932       return -1;
933     }
934
935     TParticle* mcpart = fStack->Particle(iTrack);
936
937     if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !!!!!!!!!!!!!!!!!
938
939     Int_t iLabel = mcpart->GetFirstMother();
940     if (iLabel<0){
941       AliDebug(1, "Stack label is negative, return\n");
942       return -1;
943     }
944
945     Int_t origin = -1;
946     Bool_t isFinalOpenCharm = kFALSE;
947
948     TParticle *partMother = fStack->Particle(iLabel);
949     Int_t maPdgcode = partMother->GetPdgCode();
950
951     // if the mother is charmed hadron  
952     if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
953
954          for (Int_t i=0; i<fNparents; i++){
955             if (abs(maPdgcode)==fParentSelect[0][i]){
956               isFinalOpenCharm = kTRUE;
957             }
958          }
959          if (!isFinalOpenCharm) return -1;
960
961
962           // iterate until you find B hadron as a mother or become top ancester 
963           for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
964
965              Int_t jLabel = partMother->GetFirstMother();
966              if (jLabel == -1){
967                origin = kDirectCharm;
968                return origin;
969              }
970              if (jLabel < 0){ // safety protection
971                AliDebug(1, "Stack label is negative, return\n");
972                return -1;
973              }
974
975              // if there is an ancester
976              TParticle* grandMa = fStack->Particle(jLabel);
977              Int_t grandMaPDG = grandMa->GetPdgCode();
978
979              for (Int_t j=0; j<fNparents; j++){
980                 if (abs(grandMaPDG)==fParentSelect[1][j]){
981
982                   origin = kBeautyCharm;
983                   return origin;
984                 }
985              }
986
987              partMother = grandMa;
988           } // end of iteration 
989     } // end of if
990     else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
991          for (Int_t i=0; i<fNparents; i++){
992             if (abs(maPdgcode)==fParentSelect[1][i]){
993               origin = kDirectBeauty;
994               return origin;
995             }
996          }
997     } // end of if
998
999
1000     //============ gamma ================
1001     else if ( abs(maPdgcode) == 22 ) {
1002          origin = kGamma;
1003
1004           // iterate until you find B hadron as a mother or become top ancester 
1005           for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1006
1007              Int_t jLabel = partMother->GetFirstMother();
1008              if (jLabel == -1){
1009                origin = kGamma;
1010                return origin;
1011              }
1012              if (jLabel < 0){ // safety protection
1013                AliDebug(1, "Stack label is negative, return\n");
1014                return -1;
1015              }
1016
1017              // if there is an ancester
1018              TParticle* grandMa = fStack->Particle(jLabel);
1019              Int_t grandMaPDG = grandMa->GetPdgCode();
1020
1021              for (Int_t j=0; j<fNparents; j++){
1022                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1023                   origin = kBeautyGamma;
1024                   return origin;
1025                 }
1026              }
1027
1028              partMother = grandMa;
1029           } // end of iteration 
1030
1031          return origin;
1032     } // end of if
1033
1034     //============ pi0 ================
1035     else if ( abs(maPdgcode) == 111 ) {
1036          origin = kPi0;
1037
1038           // iterate until you find B hadron as a mother or become top ancester 
1039           for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1040
1041              Int_t jLabel = partMother->GetFirstMother();
1042              if (jLabel == -1){
1043                origin = kPi0;
1044                return origin;
1045              }
1046              if (jLabel < 0){ // safety protection
1047                AliDebug(1, "Stack label is negative, return\n");
1048                return -1;
1049              }
1050
1051              // if there is an ancester
1052              TParticle* grandMa = fStack->Particle(jLabel);
1053              Int_t grandMaPDG = grandMa->GetPdgCode();
1054
1055              for (Int_t j=0; j<fNparents; j++){
1056                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1057                   origin = kBeautyPi0;
1058                   return origin;
1059                 }
1060              }
1061
1062              partMother = grandMa;
1063           } // end of iteration 
1064
1065          return origin;
1066     } // end of if
1067
1068
1069     else {
1070         origin = kElse;
1071
1072           // iterate until you find B hadron as a mother or become top ancester 
1073           for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1074
1075              Int_t jLabel = partMother->GetFirstMother();
1076              if (jLabel == -1){
1077                origin = kElse;
1078                return origin;
1079              }
1080              if (jLabel < 0){ // safety protection
1081                AliDebug(1, "Stack label is negative, return\n");
1082                return -1;
1083              }
1084
1085              // if there is an ancester
1086              TParticle* grandMa = fStack->Particle(jLabel);
1087              Int_t grandMaPDG = grandMa->GetPdgCode();
1088
1089              for (Int_t j=0; j<fNparents; j++){
1090                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1091                   origin = kBeautyElse;
1092                   return origin;
1093                 }
1094              }
1095
1096              partMother = grandMa;
1097           } // end of iteration 
1098
1099     }
1100
1101     return origin;
1102
1103 }
1104
1105 /*
1106 //_______________________________________________________________________________________________
1107 Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel) 
1108 {
1109
1110         //           
1111         // decay electron's origin 
1112         //           
1113              
1114         if (iTrackLabel < 0) {
1115                 AliDebug(1, "Stack label is negative, return\n");
1116                 return -1; 
1117         }           
1118
1119         TParticle* mcpart = fStack->Particle(iTrackLabel);
1120         Int_t iLabel = mcpart->GetFirstMother();
1121         if (iLabel<0){
1122                 AliDebug(1, "Stack label is negative, return\n");
1123                 return -1;
1124         }
1125
1126         Int_t origin = -1;
1127         Bool_t isFinalOpenCharm = kFALSE;
1128
1129         TParticle *partMother = fStack->Particle(iLabel);
1130         Int_t maPdgcode = partMother->GetPdgCode(); // get mother's pdg code
1131
1132         //beauty --------------------------
1133         if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
1134                 for (Int_t i=0; i<fNparents; i++){
1135                         if (abs(maPdgcode)==fParentSelect[1][i]){
1136                                 origin = kDirectBeauty;
1137                                 return origin;
1138                         }
1139                         else return -1; // this track is originated beauties not in the final B hadron list => excited beauties
1140                 }
1141         } // end of if
1142
1143         //charm --------------------------
1144         else if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
1145         
1146                 for (Int_t i=0; i<fNparents; i++){
1147                         if (abs(maPdgcode)==fParentSelect[0][i])
1148                                 isFinalOpenCharm = kTRUE;
1149                 }
1150                 if (!isFinalOpenCharm) return -1; // this track is originated charms not in the final D hadron list => excited charms   
1151                                                   // to prevent any possible double counting  
1152
1153                 for (Int_t i=0; i<100; i++){ // iterate 100 until you find B hadron as a mother or become top ancester
1154
1155                         Int_t jLabel = partMother->GetFirstMother();
1156                         if (jLabel == -1){
1157                                 origin = kDirectCharm;
1158                                 return origin;
1159                         }
1160                         if (jLabel < 0){ // safety protection even though not really necessary here
1161                                 AliDebug(1, "Stack label is negative, return\n");
1162                                 return -1;
1163                         }
1164
1165                         // if there is an ancester, check if it in the final B hadron list 
1166                         TParticle* grandMa = fStack->Particle(jLabel);
1167                         Int_t grandMaPDG = grandMa->GetPdgCode();
1168
1169                         for (Int_t j=0; j<fNparents; j++){
1170                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1171                                 origin = kBeautyCharm;
1172                                 return origin;
1173                                 }
1174                         } 
1175
1176                         partMother = grandMa;
1177                 } // end of iteration 
1178         } // end of if
1179
1180         //gamma --------------------------
1181         else if ( abs(maPdgcode) == 22 ) {
1182                 origin = kGamma;
1183
1184                 // iterate until you find B hadron as a mother or become top ancester 
1185                 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
1186
1187                         Int_t jLabel = partMother->GetFirstMother();
1188                         if (jLabel == -1){
1189                                 origin = kGamma;
1190                                 return origin;
1191                         }
1192                         if (jLabel < 0){ // safety protection
1193                                 AliDebug(1, "Stack label is negative, return\n");
1194                                 return -1;
1195                         }
1196
1197                         // if there is an ancester
1198                         TParticle* grandMa = fStack->Particle(jLabel);
1199                         Int_t grandMaPDG = grandMa->GetPdgCode();
1200
1201                         for (Int_t j=0; j<fNparents; j++){
1202                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1203                                         origin = kBeautyGamma;
1204                                         return origin;
1205                                 }
1206                         }
1207
1208                         partMother = grandMa;
1209                 } // end of iteration 
1210
1211                 return origin;
1212         } // end of if
1213
1214         //pi0 --------------------------
1215         else if ( abs(maPdgcode) == 111 ) {
1216                 origin = kPi0;
1217
1218                 // iterate until you find B hadron as a mother or become top ancester 
1219                 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
1220
1221                         Int_t jLabel = partMother->GetFirstMother();
1222                         if (jLabel == -1){
1223                                 origin = kPi0;
1224                                 return origin;
1225                         }
1226                         if (jLabel < 0){ // safety protection
1227                                 AliDebug(1, "Stack label is negative, return\n");
1228                                 return -1;
1229                         }
1230
1231                         // if there is an ancester
1232                         TParticle* grandMa = fStack->Particle(jLabel);
1233                         Int_t grandMaPDG = grandMa->GetPdgCode();
1234
1235                         for (Int_t j=0; j<fNparents; j++){
1236                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1237                                         origin = kBeautyPi0;
1238                                         return origin;
1239                                 }
1240                         }
1241
1242                         partMother = grandMa;
1243                 } // end of iteration 
1244
1245                 return origin;
1246         } // end of if
1247
1248
1249         //else --------------------------
1250         else {
1251                 origin = kElse;
1252
1253                 // iterate until you find B hadron as a mother or become top ancester 
1254                 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
1255
1256                         Int_t jLabel = partMother->GetFirstMother();
1257                         if (jLabel == -1){
1258                                 origin = kElse;
1259                                 return origin;
1260                         }
1261                         if (jLabel < 0){ // safety protection
1262                                 AliDebug(1, "Stack label is negative, return\n");
1263                                 return -1;
1264                         }
1265
1266                         // if there is an ancester
1267                         TParticle* grandMa = fStack->Particle(jLabel);
1268                         Int_t grandMaPDG = grandMa->GetPdgCode();
1269
1270                         for (Int_t j=0; j<fNparents; j++){
1271                                 if (abs(grandMaPDG)==fParentSelect[1][j]){
1272                                         origin = kBeautyElse;
1273                                         return origin;
1274                                 }
1275                         }
1276
1277                         partMother = grandMa;
1278                 } // end of iteration 
1279
1280         }
1281
1282         return origin;
1283
1284 }
1285 */
1286
1287 //_______________________________________________________________________________________________
1288 Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track) const
1289 {
1290         // test cuts 
1291
1292         //if (track->Pt() < 1.0) return kFALSE;
1293         //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
1294         //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
1295         //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1296         if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
1297         //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
1298
1299
1300 /*
1301         Float_t dcaR=-1;
1302         Float_t dcaZ=-1;
1303         track->GetImpactParameters(dcaR,dcaZ);
1304         if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
1305         if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
1306 */
1307         return kTRUE;
1308 }