Since it contains fixes of coding rule violations, all classes are involved. Further...
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEsecVtx.cxx
CommitLineData
259c3296 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 *
dbe3abbe 18 * Construct secondary vertex from Beauty hadron with electron and *
19 * hadrons, then apply selection criteria *
259c3296 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
42ClassImp(AliHFEsecVtx)
43
44//_______________________________________________________________________________________________
45AliHFEsecVtx::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)
dbe3abbe 57 ,fBElec(0)
259c3296 58{
59 // Default constructor
60
61 Init();
62
63}
64
65//_______________________________________________________________________________________________
66AliHFEsecVtx::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)
dbe3abbe 79 ,fBElec(p.fBElec)
259c3296 80{
81 // Copy constructor
82}
83
84//_______________________________________________________________________________________________
85AliHFEsecVtx&
86AliHFEsecVtx::operator=(const AliHFEsecVtx &)
87{
88 // Assignment operator
89
90 AliInfo("Not yet implemented.");
91 return *this;
92}
93
94//_______________________________________________________________________________________________
95AliHFEsecVtx::~AliHFEsecVtx()
96{
97 // Destructor
98
99 //cout << "Analysis Done." << endl;
100}
101
102//__________________________________________
103void AliHFEsecVtx::Init()
104{
105
dbe3abbe 106 // set pdg code and index
107
259c3296 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
41d0c3b3 126/*
259c3296 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;
41d0c3b3 170*/
259c3296 171
259c3296 172}
173
174//__________________________________________
175void AliHFEsecVtx::ResetTagVar()
176{
dbe3abbe 177 // reset tag variables
178
259c3296 179 fdistTwoSecVtx = -1;
180 fcosPhi = -1;
181 fsignedLxy = -1;
182 finvmass = -1;
183 finvmassSigma = -1;
184 fBTagged = kFALSE;
dbe3abbe 185 fBElec = kFALSE;
259c3296 186}
187
188//__________________________________________
189void AliHFEsecVtx::InitAnaPair()
190{
dbe3abbe 191 // initialize pair tagging variables
192
259c3296 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
41d0c3b3 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
259c3296 245}
246
247//_______________________________________________________________________________________________
248void AliHFEsecVtx::CreateHistograms(TString hnopt)
249{
250 // create histograms
251
dbe3abbe 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";
259c3296 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);
dbe3abbe 275 hname=hnopt+"OpenAngle_"+fkSourceLabel[isource];
276 fHistPair[isource].fOpenAngle = new TH1F(hname,hname,100,0,3.14);
259c3296 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);
dbe3abbe 283 hname=hnopt+"IPMax_"+fkSourceLabel[isource];
284 fHistPair[isource].fIPMax= new TH1F(hname,hname,500,0,5);
259c3296 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
dbe3abbe 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);
259c3296 336}
337
338//_______________________________________________________________________________________________
dbe3abbe 339void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t index2)
259c3296 340{
dbe3abbe 341 // calculate e-h pair characteristics and tag pair
342
343 Int_t sourcePart = PairCode(track1,track2);
259c3296 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);
75d81601 354
259c3296 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
dbe3abbe 393 Double_t phi = kfTrack1.GetAngleXY(kfTrack2);
394 Double_t cosphi = TMath::Cos(phi);
259c3296 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
dbe3abbe 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);
259c3296 410
dbe3abbe 411 if (TMath::Abs(dcaR1) >= TMath::Abs(dcaR2)) dcaR=dcaR1;
412 else dcaR=dcaR2;
259c3296 413
414 // fill histograms
415 fHistPair[sourcePart].fInvMass->Fill(invmass,invmassSigma);
416 fHistPair[sourcePart].fKFChi2->Fill(kfchi2);
dbe3abbe 417 fHistPair[sourcePart].fOpenAngle->Fill(phi);
259c3296 418 fHistPair[sourcePart].fCosOpenAngle->Fill(cosphi);
419 fHistPair[sourcePart].fSignedLxy->Fill(signedLxy,invmass);
420 fHistPair[sourcePart].fKFIP->Fill(kfip);
dbe3abbe 421 fHistPair[sourcePart].fIPMax->Fill(TMath::Abs(dcaR));
259c3296 422
423 // pair cuts
dbe3abbe 424 if( kfchi2 >2. ) return;
425
426 if ( signedLxy > 0.05 && cosphi > 0.5) fHistPair[sourcePart].fInvMassCut1->Fill(invmass,invmassSigma);
259c3296 427
428 // pair tagging if it passed the above cuts
dbe3abbe 429 if(signedLxy > 0. && cosphi > 0.) fHistPair[sourcePart].fInvMassCut2->Fill(invmass,invmassSigma);
259c3296 430
431
432 // pair tagging condition
dbe3abbe 433 if ( signedLxy > 0.0 && cosphi > 0) { // testing loose cut
259c3296 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//_______________________________________________________________________________________________
445void AliHFEsecVtx::RunSECVTX(AliESDtrack *track)
446{
dbe3abbe 447 // run secondary vertexing algorithm and do tagging
259c3296 448
449 ResetTagVar();
450
dbe3abbe 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
259c3296 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
259c3296 474
475 if (fBTagged) {
476 fHistTagged.fPtTaggedElec->Fill(mcpart->Pt());
dbe3abbe 477 if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){
259c3296 478 fHistTagged.fPtRightTaggedElec->Fill(mcpart->Pt());
479 }
480 else fHistTagged.fPtWrongTaggedElec->Fill(mcpart->Pt());
481 }
482
483}
484
485//_______________________________________________________________________________________________
486void AliHFEsecVtx::ApplyPairTagCut()
487{
dbe3abbe 488 // apply tagging cut for e-h pair
259c3296 489
dbe3abbe 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
259c3296 501 if (fpairedChi2[0] > 2.0) return;
502 if (fpairedInvMass[0] < 1.5) return;
dbe3abbe 503 if (fpairedSignedLxy[0] < 0.05) return;
259c3296 504
505 fBTagged = kTRUE;
259c3296 506}
507
508//_______________________________________________________________________________________________
509Bool_t AliHFEsecVtx::ApplyPairTagCut(Int_t id)
510{
dbe3abbe 511 // apply tagging cut for e-h pair of given indexed hadron
259c3296 512
513 if (fpairedChi2[id] > 2.0) return kFALSE;
514 if (fpairedInvMass[id] < 1.5) return kFALSE;
dbe3abbe 515 if (fpairedSignedLxy[id] < 0.05) return kFALSE;
259c3296 516
517 fBTagged = kTRUE;
518 return kTRUE;
519
520}
521
522//_______________________________________________________________________________________________
523Bool_t AliHFEsecVtx::ApplyTagCut(Double_t chi2)
524{
dbe3abbe 525 // apply tagging cut for secondary vertex
259c3296 526
527 if (chi2 > 2.0) return kFALSE;
528 if (finvmass < 1.5) return kFALSE;
dbe3abbe 529 if (fsignedLxy < 0.05) return kFALSE;
259c3296 530 if (fcosPhi < 0.90) return kFALSE;
531 if (fdistTwoSecVtx > 0.1) return kFALSE;
532
533 fBTagged = kTRUE;
534 return kTRUE;
535}
536
537//_______________________________________________________________________________________________
dbe3abbe 538void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2, Int_t id1, Int_t id2)
259c3296 539{
dbe3abbe 540 // apply tagging cut for e-h pair of given indexed hadron
259c3296 541
542 if(!ApplyTagCut(chi2)){
dbe3abbe 543 if(!ApplyPairTagCut(id1)) ApplyPairTagCut(id2);
259c3296 544 }
545
546}
547
548//_______________________________________________________________________________________________
549void AliHFEsecVtx::FindSECVTXCandid4Tracks(AliESDtrack *track)
550{
dbe3abbe 551 // find secondary vertex for >= 4 e-h pairs
259c3296 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
41d0c3b3 560 if(fPairTagged>20) return; // protection
561
259c3296 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
dbe3abbe 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]]);
259c3296 598
599}
600
601//_______________________________________________________________________________________________
602void AliHFEsecVtx::FindSECVTXCandid3Tracks(AliESDtrack *track)
603{
dbe3abbe 604 // find secondary vertex for 3 e-h pairs
259c3296 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
dbe3abbe 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]);
259c3296 637}
638
639//_______________________________________________________________________________________________
640void AliHFEsecVtx::FindSECVTXCandid2Tracks(AliESDtrack *track)
641{
dbe3abbe 642 // find secondary vertex for 2 e-h pairs
259c3296 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
dbe3abbe 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);
259c3296 664}
665
666//_______________________________________________________________________________________________
667void AliHFEsecVtx::CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
668{
dbe3abbe 669 // calculate secondary vertex properties
259c3296 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);
dbe3abbe 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
259c3296 727}
728
729//_______________________________________________________________________________________________
730Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track)
731{
dbe3abbe 732 // return mc pid
259c3296 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//_______________________________________________________________________________________________
743Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3, AliESDtrack* track4)
744{
dbe3abbe 745 // return 4 track secondary vertex chi2
259c3296 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//_______________________________________________________________________________________________
764Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
765{
dbe3abbe 766 // return 3 track secondary vertex chi2
259c3296 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//_______________________________________________________________________________________________
783Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2)
784{
dbe3abbe 785 // return 2 track secondary vertex chi2
259c3296 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//_______________________________________________________________________________________________
800Int_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//_______________________________________________________________________________________________
882Int_t AliHFEsecVtx::PairCode(AliESDtrack* track1,AliESDtrack* track2)
883{
884
885 //
dbe3abbe 886 // return pair code which is predefinded as:
887 // kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
259c3296 888 //
889
890 Int_t pairOriginsPDG = PairOrigin(track1,track2);
891
dbe3abbe 892 Int_t sourcePart = kElse;
259c3296 893
894 if (pairOriginsPDG < 0) {
dbe3abbe 895 sourcePart = kBeautyElse;
259c3296 896 }
897 for (Int_t i=0; i<fNparents; i++){
898 if (abs(pairOriginsPDG)==fParentSelect[0][i]) {
dbe3abbe 899 if (pairOriginsPDG>0) sourcePart = kDirectCharm;
259c3296 900 if (pairOriginsPDG<0) {
dbe3abbe 901 sourcePart = kBeautyCharm;
259c3296 902 }
903 }
904 if (abs(pairOriginsPDG)==fParentSelect[1][i]) {
905 if (pairOriginsPDG>0) {
dbe3abbe 906 sourcePart = kDirectBeauty;
259c3296 907 }
dbe3abbe 908 if (pairOriginsPDG<0) return kElse;
259c3296 909 }
910 }
dbe3abbe 911 if (pairOriginsPDG == 22) sourcePart = kGamma;
259c3296 912 if (pairOriginsPDG == -22) {
dbe3abbe 913 sourcePart = kBeautyGamma;
259c3296 914 }
dbe3abbe 915 if (pairOriginsPDG == 111) sourcePart = kPi0;
259c3296 916 if (pairOriginsPDG == -111) {
dbe3abbe 917 sourcePart = kBeautyPi0;
259c3296 918 }
919
920 return sourcePart;
921
922}
923
dbe3abbe 924//_______________________________________________________________________________________________
925Int_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
41d0c3b3 979 for (Int_t j=0; j<fNparents; j++){
980 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 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
41d0c3b3 1021 for (Int_t j=0; j<fNparents; j++){
1022 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 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
41d0c3b3 1055 for (Int_t j=0; j<fNparents; j++){
1056 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 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
41d0c3b3 1089 for (Int_t j=0; j<fNparents; j++){
1090 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 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/*
259c3296 1106//_______________________________________________________________________________________________
1107Int_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;
dbe3abbe 1127 Bool_t isFinalOpenCharm = kFALSE;
259c3296 1128
1129 TParticle *partMother = fStack->Particle(iLabel);
1130 Int_t maPdgcode = partMother->GetPdgCode(); // get mother's pdg code
1131
1132 //beauty --------------------------
dbe3abbe 1133 if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
259c3296 1134 for (Int_t i=0; i<fNparents; i++){
1135 if (abs(maPdgcode)==fParentSelect[1][i]){
dbe3abbe 1136 origin = kDirectBeauty;
259c3296 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 --------------------------
dbe3abbe 1144 else if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
259c3296 1145
1146 for (Int_t i=0; i<fNparents; i++){
1147 if (abs(maPdgcode)==fParentSelect[0][i])
dbe3abbe 1148 isFinalOpenCharm = kTRUE;
259c3296 1149 }
dbe3abbe 1150 if (!isFinalOpenCharm) return -1; // this track is originated charms not in the final D hadron list => excited charms
259c3296 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){
dbe3abbe 1157 origin = kDirectCharm;
259c3296 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]){
dbe3abbe 1171 origin = kBeautyCharm;
259c3296 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 ) {
dbe3abbe 1182 origin = kGamma;
259c3296 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){
dbe3abbe 1189 origin = kGamma;
259c3296 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]){
dbe3abbe 1203 origin = kBeautyGamma;
259c3296 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 ) {
dbe3abbe 1216 origin = kPi0;
259c3296 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){
dbe3abbe 1223 origin = kPi0;
259c3296 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]){
dbe3abbe 1237 origin = kBeautyPi0;
259c3296 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 {
dbe3abbe 1251 origin = kElse;
259c3296 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){
dbe3abbe 1258 origin = kElse;
259c3296 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]){
dbe3abbe 1272 origin = kBeautyElse;
259c3296 1273 return origin;
1274 }
1275 }
1276
1277 partMother = grandMa;
1278 } // end of iteration
1279
1280 }
1281
1282 return origin;
1283
1284}
dbe3abbe 1285*/
259c3296 1286
1287//_______________________________________________________________________________________________
75d81601 1288Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track) const
259c3296 1289{
dbe3abbe 1290 // test cuts
259c3296 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}