]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliHFEsecVtx.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Classes for efficiency corrections (Xaver, Roberta)
[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);
259c3296 354 AliKFParticle kfSecondary(kfTrack1,kfTrack2);
355
356 // copy primary vertex from ESD
357 AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
358 if( primVtxCopy.GetNDF() <1 ) return;
359
360 //primary vertex point
361 Double_t pvx = primVtxCopy.GetX();
362 Double_t pvy = primVtxCopy.GetY();
363 //Double_t pvz = primVtxCopy.GetZ();
364
365 //secondary vertex point from kf particle
366 Double_t kfx = kfSecondary.GetX();
367 Double_t kfy = kfSecondary.GetY();
368 //Double_t kfz = kfSecondary.GetZ();
369
370 //momentum at the decay point from kf particle
371 Double_t kfpx = kfSecondary.GetPx();
372 Double_t kfpy = kfSecondary.GetPy();
373 //Double_t kfpz = kfSecondary.GetPz();
374
375
376 Double_t dx = kfx-pvx;
377 Double_t dy = kfy-pvy;
378
379
380
381 // discriminating variables ----------------------------------------------------------
382
383 // invariant mass of the KF particle
384 Double_t invmass = -1;
385 Double_t invmassSigma = -1;
386 kfSecondary.GetMass(invmass,invmassSigma);
387
388 // chi2 of the KF particle
389 Double_t kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
390
391 // opening angle between two particles in XY plane
dbe3abbe 392 Double_t phi = kfTrack1.GetAngleXY(kfTrack2);
393 Double_t cosphi = TMath::Cos(phi);
259c3296 394
395 // projection of kf vertex vector to the kf momentum direction
396 Double_t costheta = ( dx*kfpx + dy*kfpy)/TMath::Sqrt(dx*dx+dy*dy)*TMath::Sqrt(kfpx*kfpx + kfpy*kfpy);
397 Double_t signedLxy = TMath::Sqrt(dx*dx+dy*dy)*costheta;
398
399 // DCA from primary to e-h KF particle (impact parameter of KF particle)
400 Double_t vtx[2]={pvx, pvy};
401 Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
402
403
dbe3abbe 404 Float_t dcaR=-1;
405 Float_t dcaR1=-1, dcaR2=-1;
406 Float_t dcaZ1=-1, dcaZ2=-1;
407 track1->GetImpactParameters(dcaR1,dcaZ1);
408 track2->GetImpactParameters(dcaR2,dcaZ2);
259c3296 409
dbe3abbe 410 if (TMath::Abs(dcaR1) >= TMath::Abs(dcaR2)) dcaR=dcaR1;
411 else dcaR=dcaR2;
259c3296 412
413 // fill histograms
414 fHistPair[sourcePart].fInvMass->Fill(invmass,invmassSigma);
415 fHistPair[sourcePart].fKFChi2->Fill(kfchi2);
dbe3abbe 416 fHistPair[sourcePart].fOpenAngle->Fill(phi);
259c3296 417 fHistPair[sourcePart].fCosOpenAngle->Fill(cosphi);
418 fHistPair[sourcePart].fSignedLxy->Fill(signedLxy,invmass);
419 fHistPair[sourcePart].fKFIP->Fill(kfip);
dbe3abbe 420 fHistPair[sourcePart].fIPMax->Fill(TMath::Abs(dcaR));
259c3296 421
422 // pair cuts
dbe3abbe 423 if( kfchi2 >2. ) return;
424
425 if ( signedLxy > 0.05 && cosphi > 0.5) fHistPair[sourcePart].fInvMassCut1->Fill(invmass,invmassSigma);
259c3296 426
427 // pair tagging if it passed the above cuts
dbe3abbe 428 if(signedLxy > 0. && cosphi > 0.) fHistPair[sourcePart].fInvMassCut2->Fill(invmass,invmassSigma);
259c3296 429
430
431 // pair tagging condition
dbe3abbe 432 if ( signedLxy > 0.0 && cosphi > 0) { // testing loose cut
259c3296 433 //if ( signedLxy > 0.06 && cosphi > 0) {
434 fpairedTrackID[fPairTagged] = index2;
435 fpairedChi2[fPairTagged] = kfchi2;
436 fpairedInvMass[fPairTagged] = invmass;
437 fpairedSignedLxy[fPairTagged] = signedLxy;
438 fPairTagged++;
439 }
440
441}
442
443//_______________________________________________________________________________________________
444void AliHFEsecVtx::RunSECVTX(AliESDtrack *track)
445{
dbe3abbe 446 // run secondary vertexing algorithm and do tagging
259c3296 447
448 ResetTagVar();
449
dbe3abbe 450 Int_t imclabel = TMath::Abs(track->GetLabel());
451 if(imclabel<0) return;
452 TParticle* mcpart = fStack->Particle(imclabel);
453 Int_t esource = GetElectronSource(imclabel);
454 if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){
455 fHistTagged.fPtBeautyElec->Fill(mcpart->Pt());
456 fBElec = kTRUE;
457 }
458
459
259c3296 460 if (fPairTagged >= 4) {
461 FindSECVTXCandid4Tracks(track);
462 }
463 else if (fPairTagged == 3) {
464 FindSECVTXCandid3Tracks(track);
465 }
466 else if (fPairTagged == 2) {
467 FindSECVTXCandid2Tracks(track);
468 }
469 else if (fPairTagged == 1) {
470 ApplyPairTagCut();
471 }
472
259c3296 473
474 if (fBTagged) {
475 fHistTagged.fPtTaggedElec->Fill(mcpart->Pt());
dbe3abbe 476 if (esource == kDirectBeauty || esource == kBeautyCharm || esource == kBeautyGamma || esource == kBeautyPi0 || esource == kBeautyElse){
259c3296 477 fHistTagged.fPtRightTaggedElec->Fill(mcpart->Pt());
478 }
479 else fHistTagged.fPtWrongTaggedElec->Fill(mcpart->Pt());
480 }
481
482}
483
484//_______________________________________________________________________________________________
485void AliHFEsecVtx::ApplyPairTagCut()
486{
dbe3abbe 487 // apply tagging cut for e-h pair
259c3296 488
dbe3abbe 489 if (fBElec){
490 fHistTagged.fInvmassBeautyElec2trkVtx->Fill(fpairedInvMass[0]);
491 fHistTagged.fSignedLxyBeautyElec2trkVtx->Fill(fpairedSignedLxy[0]);
492 fHistTagged.fChi2BeautyElec2trkVtx->Fill(fpairedChi2[0]);
493 }
494 else {
495 fHistTagged.fInvmassNotBeautyElec2trkVtx->Fill(fpairedInvMass[0]);
496 fHistTagged.fSignedLxyNotBeautyElec2trkVtx->Fill(fpairedSignedLxy[0]);
497 fHistTagged.fChi2NotBeautyElec2trkVtx->Fill(fpairedChi2[0]);
498 }
499
259c3296 500 if (fpairedChi2[0] > 2.0) return;
501 if (fpairedInvMass[0] < 1.5) return;
dbe3abbe 502 if (fpairedSignedLxy[0] < 0.05) return;
259c3296 503
504 fBTagged = kTRUE;
259c3296 505}
506
507//_______________________________________________________________________________________________
508Bool_t AliHFEsecVtx::ApplyPairTagCut(Int_t id)
509{
dbe3abbe 510 // apply tagging cut for e-h pair of given indexed hadron
259c3296 511
512 if (fpairedChi2[id] > 2.0) return kFALSE;
513 if (fpairedInvMass[id] < 1.5) return kFALSE;
dbe3abbe 514 if (fpairedSignedLxy[id] < 0.05) return kFALSE;
259c3296 515
516 fBTagged = kTRUE;
517 return kTRUE;
518
519}
520
521//_______________________________________________________________________________________________
522Bool_t AliHFEsecVtx::ApplyTagCut(Double_t chi2)
523{
dbe3abbe 524 // apply tagging cut for secondary vertex
259c3296 525
526 if (chi2 > 2.0) return kFALSE;
527 if (finvmass < 1.5) return kFALSE;
dbe3abbe 528 if (fsignedLxy < 0.05) return kFALSE;
259c3296 529 if (fcosPhi < 0.90) return kFALSE;
530 if (fdistTwoSecVtx > 0.1) return kFALSE;
531
532 fBTagged = kTRUE;
533 return kTRUE;
534}
535
536//_______________________________________________________________________________________________
dbe3abbe 537void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2, Int_t id1, Int_t id2)
259c3296 538{
dbe3abbe 539 // apply tagging cut for e-h pair of given indexed hadron
259c3296 540
541 if(!ApplyTagCut(chi2)){
dbe3abbe 542 if(!ApplyPairTagCut(id1)) ApplyPairTagCut(id2);
259c3296 543 }
544
545}
546
547//_______________________________________________________________________________________________
548void AliHFEsecVtx::FindSECVTXCandid4Tracks(AliESDtrack *track)
549{
dbe3abbe 550 // find secondary vertex for >= 4 e-h pairs
259c3296 551
552 // sort pair in increasing order (kFALSE-increasing order)
553 Int_t index[20];
554 Int_t indexA[4];
555 Int_t indexB[3];
556 Double_t sevchi2[4];
557 AliESDtrack *htrack[4];
558
41d0c3b3 559 if(fPairTagged>20) return; // protection
560
259c3296 561 TMath::Sort(fPairTagged,fpairedChi2,index,kFALSE);
562
563 // select 4 partner tracks retruning smallest pair chi2
564
565 for (Int_t i=0; i<4; i++){
566 htrack[i] = fESD1->GetTrack(fpairedTrackID[index[i]]);
567 }
568
569 // calculated chi2 with 1 electron and 3 partner tracks
570 for (Int_t i=0; i<4; i++){
571 sevchi2[i] = GetSecVtxChi2(track, htrack[fid[i][0]], htrack[fid[i][1]], htrack[fid[i][2]]);
572 }
573
574 // select 3 partner tracks retruning smallest pair chi2
575 // [think] if two smallest chi2 are similar, have to think about better handling of selection
576 TMath::Sort(4,sevchi2,indexA,kFALSE);
577
578 // calculated chi2 with 1 electron and 2 partner tracks
579 for (Int_t i=0; i<3; i++){
580 sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
581 }
582
583 // select 2 partner tracks retruning smallest pair chi2
584 TMath::Sort(3,sevchi2,indexB,kFALSE);
585
586 // calculate secondary vertex quality variables with 1 electron and 2 hadrons
587 CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
588
dbe3abbe 589 if (fBElec){
590 fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
591 }
592 else {
593 fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
594 }
595
596 ApplyVtxTagCut(sevchi2[indexB[0]],index[fia[indexA[0]][indexB[0]][0]],index[fia[indexA[0]][indexB[0]][1]]);
259c3296 597
598}
599
600//_______________________________________________________________________________________________
601void AliHFEsecVtx::FindSECVTXCandid3Tracks(AliESDtrack *track)
602{
dbe3abbe 603 // find secondary vertex for 3 e-h pairs
259c3296 604
605 // sort pair in increasing order (kFALSE-increasing order)
606 Int_t indexA[1] = { 0 };
607 Int_t indexB[3];
608 Double_t sevchi2[3];
609 AliESDtrack *htrack[3];
610
611 // select 4 partner tracks retruning smallest pair chi2
612
613 for (Int_t i=0; i<3; i++){
614 htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
615 }
616
617 // calculated chi2 with 1 electron and 2 partner tracks
618 for (Int_t i=0; i<3; i++){
619 sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
620 }
621
622 // select 2 partner tracks retruning smallest pair chi2
623 TMath::Sort(3,sevchi2,indexB,kFALSE);
624
625 // calculate secondary vertex quality variables with 1 electron and 2 hadrons
626 CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
627
dbe3abbe 628 if (fBElec){
629 fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
630 }
631 else {
632 fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[indexB[0]]);
633 }
634
635 ApplyVtxTagCut(sevchi2[indexB[0]],fia[indexA[0]][indexB[0]][0],fia[indexA[0]][indexB[0]][1]);
259c3296 636}
637
638//_______________________________________________________________________________________________
639void AliHFEsecVtx::FindSECVTXCandid2Tracks(AliESDtrack *track)
640{
dbe3abbe 641 // find secondary vertex for 2 e-h pairs
259c3296 642
643 Double_t sevchi2[1];
644 AliESDtrack *htrack[2];
645
646 for (Int_t i=0; i<2; i++){
647 htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
648 }
649
650 sevchi2[0] = GetSecVtxChi2(track, htrack[0], htrack[1]);
651
652 // calculate secondary vertex quality variables with 1 electron and 2 hadrons
653 CalcSECVTXProperty(track,htrack[0],htrack[1]);
654
dbe3abbe 655 if (fBElec){
656 fHistTagged.fChi2BeautyElecSecVtx->Fill(sevchi2[0]);
657 }
658 else {
659 fHistTagged.fChi2NotBeautyElecSecVtx->Fill(sevchi2[0]);
660 }
661
662 ApplyVtxTagCut(sevchi2[0],0,1);
259c3296 663}
664
665//_______________________________________________________________________________________________
666void AliHFEsecVtx::CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
667{
dbe3abbe 668 // calculate secondary vertex properties
259c3296 669
670 Int_t pdg1 = GetMCPID(track1);
671 Int_t pdg2 = GetMCPID(track2);
672 Int_t pdg3 = GetMCPID(track3);
673
674 AliKFParticle::SetField(fESD1->GetMagneticField());
675 AliKFParticle kfTrack1(*track1, pdg1);
676 AliKFParticle kfTrack2(*track2, pdg2);
677 AliKFParticle kfTrack3(*track3, pdg3);
678
679 AliKFParticle kfSecondary12(kfTrack1,kfTrack2);
680 AliKFParticle kfSecondary13(kfTrack1,kfTrack3);
681 AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
682
683 // copy primary vertex from ESD
684 AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
685 //printf("primary ndf= %d\n",primVtxCopy.GetNDF());
686 if( primVtxCopy.GetNDF() <1 ) return;
687
688 Double_t kdx12 = kfSecondary12.GetX()-primVtxCopy.GetX();
689 Double_t kdy12 = kfSecondary12.GetY()-primVtxCopy.GetY();
690 //Double_t kdz12 = kfSecondary12.GetZ()-primVtxCopy.GetZ();
691
692 Double_t kdx13 = kfSecondary13.GetX()-primVtxCopy.GetX();
693 Double_t kdy13 = kfSecondary13.GetY()-primVtxCopy.GetY();
694 //Double_t kdz13 = kfSecondary13.GetZ()-primVtxCopy.GetZ();
695
696 Double_t kdx = kfSecondary.GetX()-primVtxCopy.GetX();
697 Double_t kdy = kfSecondary.GetY()-primVtxCopy.GetY();
698 //Double_t kdz = kfSecondary.GetZ()-primVtxCopy.GetZ();
699
700 // calculate distance and angle between two secvtxes
701 fdistTwoSecVtx = TMath::Sqrt((kdx12-kdx13)*(kdx12-kdx13) + (kdy12-kdy13)*(kdy12-kdy13));
702 fcosPhi = ( kdx12*kdx13 + kdy12*kdy13 ) / ( TMath::Sqrt(kdx12*kdx12+kdy12*kdy12)*TMath::Sqrt(kdx13*kdx13+kdy13*kdy13) );
703 //Double_t lengthdiff = TMath::Abs(TMath::Sqrt(kdx12*kdx12+kdy12*kdy12) - TMath::Sqrt(kdx13*kdx13+kdy13*kdy13));
704
705 // calculate angle between secondary vertex vector and secondary particle momentum vector in transverse plane
706 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());
707 // calculate signed Lxy
708 fsignedLxy = TMath::Sqrt(kdx*kdx+kdy*kdy)*cosTheta;
709
710 // calculate invariant mass of the kf particle
711 kfSecondary.GetMass(finvmass,finvmassSigma);
dbe3abbe 712
713 if (fBElec){
714 fHistTagged.fInvmassBeautyElecSecVtx->Fill(finvmass);
715 fHistTagged.fSignedLxyBeautyElecSecVtx->Fill(fsignedLxy);
716 fHistTagged.fDistTwoVtxBeautyElecSecVtx->Fill(fdistTwoSecVtx);
717 fHistTagged.fCosPhiBeautyElecSecVtx->Fill(fcosPhi);
718 }
719 else {
720 fHistTagged.fInvmassNotBeautyElecSecVtx->Fill(finvmass);
721 fHistTagged.fSignedLxyNotBeautyElecSecVtx->Fill(fsignedLxy);
722 fHistTagged.fDistTwoVtxNotBeautyElecSecVtx->Fill(fdistTwoSecVtx);
723 fHistTagged.fCosPhiNotBeautyElecSecVtx->Fill(fcosPhi);
724 }
725
259c3296 726}
727
728//_______________________________________________________________________________________________
729Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track)
730{
dbe3abbe 731 // return mc pid
259c3296 732
733 Int_t label = TMath::Abs(track->GetLabel());
734 TParticle* mcpart = fStack->Particle(label);
735 if ( !mcpart ) return 0;
736 Int_t pdgCode = mcpart->GetPdgCode();
737
738 return pdgCode;
739}
740
741//_______________________________________________________________________________________________
742Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3, AliESDtrack* track4)
743{
dbe3abbe 744 // return 4 track secondary vertex chi2
259c3296 745
746 Int_t pdg1 = GetMCPID(track1);
747 Int_t pdg2 = GetMCPID(track2);
748 Int_t pdg3 = GetMCPID(track3);
749 Int_t pdg4 = GetMCPID(track4);
750
751 AliKFParticle::SetField(fESD1->GetMagneticField());
752 AliKFParticle kfTrack1(*track1, pdg1);
753 AliKFParticle kfTrack2(*track2, pdg2);
754 AliKFParticle kfTrack3(*track3, pdg3);
755 AliKFParticle kfTrack4(*track4, pdg4);
756 AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3,kfTrack4);
757
758 return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
759
760}
761
762//_______________________________________________________________________________________________
763Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
764{
dbe3abbe 765 // return 3 track secondary vertex chi2
259c3296 766
767 Int_t pdg1 = GetMCPID(track1);
768 Int_t pdg2 = GetMCPID(track2);
769 Int_t pdg3 = GetMCPID(track3);
770
771 AliKFParticle::SetField(fESD1->GetMagneticField());
772 AliKFParticle kfTrack1(*track1, pdg1);
773 AliKFParticle kfTrack2(*track2, pdg2);
774 AliKFParticle kfTrack3(*track3, pdg3);
775 AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
776
777 return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
778
779}
780
781//_______________________________________________________________________________________________
782Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2)
783{
dbe3abbe 784 // return 2 track secondary vertex chi2
259c3296 785
786 Int_t pdg1 = GetMCPID(track1);
787 Int_t pdg2 = GetMCPID(track2);
788
789 AliKFParticle::SetField(fESD1->GetMagneticField());
790 AliKFParticle kfTrack1(*track1, pdg1);
791 AliKFParticle kfTrack2(*track2, pdg2);
792 AliKFParticle kfSecondary(kfTrack1,kfTrack2);
793
794 return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
795
796}
797
798//_______________________________________________________________________________________________
799Int_t AliHFEsecVtx::PairOrigin(AliESDtrack* track1, AliESDtrack* track2)
800{
801
802 //
803 // return pdg code of the origin(source) of the pair
804 //
805 //
806 // ---*---*---*-----ancester A----- track1
807 // |____*______
808 // |______ track2
809 // => if they originated from same ancester,
810 // then return "the absolute value of pdg code of ancester A"
811 //
812 // ---*---*---B hadron-----ancester A----- track1
813 // |____*______
814 // |______ track2
815 // => if they originated from same ancester, and this ancester originally comes from B hadrons
816 // then return -1*"the absolute value of pdg code of ancester A"
817 //
818 // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
819 //
820
821 TParticle* part1 = fStack->Particle(TMath::Abs(track1->GetLabel()));
822 TParticle* part2 = fStack->Particle(TMath::Abs(track2->GetLabel()));
823 if (!(part1)) return 0;
824 if (!(part2)) return 0;
825
826 TParticle* part1Crtgen = part1; // copy track into current generation particle
827 TParticle* part2Crtgen = part2; // copy track into current generation particle
828
829
830 Int_t sourcePDG = 0;
831
832 // if the two tracks' mother's label is same, get the mother info
833 // in case of charm, check if it originated from beauty
834 for (Int_t i=0; i<100; i++){ // iterate 100
835 Int_t iLabel = part1Crtgen->GetFirstMother(); //label of mother of current generation for 1st partilce
836 if (iLabel < 0) return 0;
837
838 for (Int_t j=0; j<100; j++){ // iterate 100
839 Int_t jLabel = part2Crtgen->GetFirstMother(); //label of mother of current generation for 2nd partilce
840 if (jLabel < 0) return 0; // if jLabel == -1
841
842 if (iLabel == jLabel){ // check if two tracks are originated from same mother
843 TParticle* thismother = fStack->Particle(jLabel); // if yes, get "thismother" info
844 sourcePDG = abs(thismother->GetPdgCode()); // get the pdg code of "this mother"
845
846 // check ancester to see if it is originally from beauty
847 for (Int_t k=0; k<10; k++){ // check up to 10 ancesters
848 Int_t ancesterLabel = thismother->GetFirstMother();
849 if (ancesterLabel < 0) return sourcePDG; // if "thismoter" doesn't have mother anymore, return thismother's pdg
850
851 TParticle* thisancester = fStack->Particle(ancesterLabel);
852 Int_t ancesterPDG = abs(thisancester->GetPdgCode());
853
854 for (Int_t l=0; l<fNparents; l++){
855 if (abs(ancesterPDG)==fParentSelect[1][l]){
856 sourcePDG = -1*sourcePDG; // multiply -1 for charm from bottom
857 return sourcePDG;
858 }
859 }
860 thismother = thisancester;
861 }
862
863 }
864 part2Crtgen = fStack->Particle(jLabel); // if their mother is different, go up to previous generation of 2nd particle
865 }
866 part1Crtgen = fStack->Particle(iLabel); // if their mother is different, go up to previous generation of 2nd particle
867
868 // if you don't find additionional(2nd particle) track originated from a given beauty hadron, break to save time
869 Int_t motherPDGtmp = abs(part1Crtgen->GetPdgCode());
870 for (Int_t l=0; l<fNparents; l++){
871 if (abs(motherPDGtmp)==fParentSelect[1][l]) return sourcePDG;
872 }
873
874 }
875
876 return sourcePDG;
877
878}
879
880//_______________________________________________________________________________________________
881Int_t AliHFEsecVtx::PairCode(AliESDtrack* track1,AliESDtrack* track2)
882{
883
884 //
dbe3abbe 885 // return pair code which is predefinded as:
886 // kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
259c3296 887 //
888
889 Int_t pairOriginsPDG = PairOrigin(track1,track2);
890
dbe3abbe 891 Int_t sourcePart = kElse;
259c3296 892
893 if (pairOriginsPDG < 0) {
dbe3abbe 894 sourcePart = kBeautyElse;
259c3296 895 }
896 for (Int_t i=0; i<fNparents; i++){
897 if (abs(pairOriginsPDG)==fParentSelect[0][i]) {
dbe3abbe 898 if (pairOriginsPDG>0) sourcePart = kDirectCharm;
259c3296 899 if (pairOriginsPDG<0) {
dbe3abbe 900 sourcePart = kBeautyCharm;
259c3296 901 }
902 }
903 if (abs(pairOriginsPDG)==fParentSelect[1][i]) {
904 if (pairOriginsPDG>0) {
dbe3abbe 905 sourcePart = kDirectBeauty;
259c3296 906 }
dbe3abbe 907 if (pairOriginsPDG<0) return kElse;
259c3296 908 }
909 }
dbe3abbe 910 if (pairOriginsPDG == 22) sourcePart = kGamma;
259c3296 911 if (pairOriginsPDG == -22) {
dbe3abbe 912 sourcePart = kBeautyGamma;
259c3296 913 }
dbe3abbe 914 if (pairOriginsPDG == 111) sourcePart = kPi0;
259c3296 915 if (pairOriginsPDG == -111) {
dbe3abbe 916 sourcePart = kBeautyPi0;
259c3296 917 }
918
919 return sourcePart;
920
921}
922
dbe3abbe 923//_______________________________________________________________________________________________
924Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack)
925{
926
927 // return decay electron's origin
928
929 if (iTrack < 0) {
930 AliDebug(1, "Stack label is negative, return\n");
931 return -1;
932 }
933
934 TParticle* mcpart = fStack->Particle(iTrack);
935
936 if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !!!!!!!!!!!!!!!!!
937
938 Int_t iLabel = mcpart->GetFirstMother();
939 if (iLabel<0){
940 AliDebug(1, "Stack label is negative, return\n");
941 return -1;
942 }
943
944 Int_t origin = -1;
945 Bool_t isFinalOpenCharm = kFALSE;
946
947 TParticle *partMother = fStack->Particle(iLabel);
948 Int_t maPdgcode = partMother->GetPdgCode();
949
950 // if the mother is charmed hadron
951 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
952
953 for (Int_t i=0; i<fNparents; i++){
954 if (abs(maPdgcode)==fParentSelect[0][i]){
955 isFinalOpenCharm = kTRUE;
956 }
957 }
958 if (!isFinalOpenCharm) return -1;
959
960
961 // iterate until you find B hadron as a mother or become top ancester
962 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
963
964 Int_t jLabel = partMother->GetFirstMother();
965 if (jLabel == -1){
966 origin = kDirectCharm;
967 return origin;
968 }
969 if (jLabel < 0){ // safety protection
970 AliDebug(1, "Stack label is negative, return\n");
971 return -1;
972 }
973
974 // if there is an ancester
975 TParticle* grandMa = fStack->Particle(jLabel);
976 Int_t grandMaPDG = grandMa->GetPdgCode();
977
41d0c3b3 978 for (Int_t j=0; j<fNparents; j++){
979 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 980
981 origin = kBeautyCharm;
982 return origin;
983 }
984 }
985
986 partMother = grandMa;
987 } // end of iteration
988 } // end of if
989 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
990 for (Int_t i=0; i<fNparents; i++){
991 if (abs(maPdgcode)==fParentSelect[1][i]){
992 origin = kDirectBeauty;
993 return origin;
994 }
995 }
996 } // end of if
997
998
999 //============ gamma ================
1000 else if ( abs(maPdgcode) == 22 ) {
1001 origin = kGamma;
1002
1003 // iterate until you find B hadron as a mother or become top ancester
1004 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1005
1006 Int_t jLabel = partMother->GetFirstMother();
1007 if (jLabel == -1){
1008 origin = kGamma;
1009 return origin;
1010 }
1011 if (jLabel < 0){ // safety protection
1012 AliDebug(1, "Stack label is negative, return\n");
1013 return -1;
1014 }
1015
1016 // if there is an ancester
1017 TParticle* grandMa = fStack->Particle(jLabel);
1018 Int_t grandMaPDG = grandMa->GetPdgCode();
1019
41d0c3b3 1020 for (Int_t j=0; j<fNparents; j++){
1021 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1022 origin = kBeautyGamma;
1023 return origin;
1024 }
1025 }
1026
1027 partMother = grandMa;
1028 } // end of iteration
1029
1030 return origin;
1031 } // end of if
1032
1033 //============ pi0 ================
1034 else if ( abs(maPdgcode) == 111 ) {
1035 origin = kPi0;
1036
1037 // iterate until you find B hadron as a mother or become top ancester
1038 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1039
1040 Int_t jLabel = partMother->GetFirstMother();
1041 if (jLabel == -1){
1042 origin = kPi0;
1043 return origin;
1044 }
1045 if (jLabel < 0){ // safety protection
1046 AliDebug(1, "Stack label is negative, return\n");
1047 return -1;
1048 }
1049
1050 // if there is an ancester
1051 TParticle* grandMa = fStack->Particle(jLabel);
1052 Int_t grandMaPDG = grandMa->GetPdgCode();
1053
41d0c3b3 1054 for (Int_t j=0; j<fNparents; j++){
1055 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1056 origin = kBeautyPi0;
1057 return origin;
1058 }
1059 }
1060
1061 partMother = grandMa;
1062 } // end of iteration
1063
1064 return origin;
1065 } // end of if
1066
1067
1068 else {
1069 origin = kElse;
1070
1071 // iterate until you find B hadron as a mother or become top ancester
1072 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1073
1074 Int_t jLabel = partMother->GetFirstMother();
1075 if (jLabel == -1){
1076 origin = kElse;
1077 return origin;
1078 }
1079 if (jLabel < 0){ // safety protection
1080 AliDebug(1, "Stack label is negative, return\n");
1081 return -1;
1082 }
1083
1084 // if there is an ancester
1085 TParticle* grandMa = fStack->Particle(jLabel);
1086 Int_t grandMaPDG = grandMa->GetPdgCode();
1087
41d0c3b3 1088 for (Int_t j=0; j<fNparents; j++){
1089 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1090 origin = kBeautyElse;
1091 return origin;
1092 }
1093 }
1094
1095 partMother = grandMa;
1096 } // end of iteration
1097
1098 }
1099
1100 return origin;
1101
1102}
1103
1104/*
259c3296 1105//_______________________________________________________________________________________________
1106Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
1107{
1108
1109 //
1110 // decay electron's origin
1111 //
1112
1113 if (iTrackLabel < 0) {
1114 AliDebug(1, "Stack label is negative, return\n");
1115 return -1;
1116 }
1117
1118 TParticle* mcpart = fStack->Particle(iTrackLabel);
1119 Int_t iLabel = mcpart->GetFirstMother();
1120 if (iLabel<0){
1121 AliDebug(1, "Stack label is negative, return\n");
1122 return -1;
1123 }
1124
1125 Int_t origin = -1;
dbe3abbe 1126 Bool_t isFinalOpenCharm = kFALSE;
259c3296 1127
1128 TParticle *partMother = fStack->Particle(iLabel);
1129 Int_t maPdgcode = partMother->GetPdgCode(); // get mother's pdg code
1130
1131 //beauty --------------------------
dbe3abbe 1132 if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
259c3296 1133 for (Int_t i=0; i<fNparents; i++){
1134 if (abs(maPdgcode)==fParentSelect[1][i]){
dbe3abbe 1135 origin = kDirectBeauty;
259c3296 1136 return origin;
1137 }
1138 else return -1; // this track is originated beauties not in the final B hadron list => excited beauties
1139 }
1140 } // end of if
1141
1142 //charm --------------------------
dbe3abbe 1143 else if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
259c3296 1144
1145 for (Int_t i=0; i<fNparents; i++){
1146 if (abs(maPdgcode)==fParentSelect[0][i])
dbe3abbe 1147 isFinalOpenCharm = kTRUE;
259c3296 1148 }
dbe3abbe 1149 if (!isFinalOpenCharm) return -1; // this track is originated charms not in the final D hadron list => excited charms
259c3296 1150 // to prevent any possible double counting
1151
1152 for (Int_t i=0; i<100; i++){ // iterate 100 until you find B hadron as a mother or become top ancester
1153
1154 Int_t jLabel = partMother->GetFirstMother();
1155 if (jLabel == -1){
dbe3abbe 1156 origin = kDirectCharm;
259c3296 1157 return origin;
1158 }
1159 if (jLabel < 0){ // safety protection even though not really necessary here
1160 AliDebug(1, "Stack label is negative, return\n");
1161 return -1;
1162 }
1163
1164 // if there is an ancester, check if it in the final B hadron list
1165 TParticle* grandMa = fStack->Particle(jLabel);
1166 Int_t grandMaPDG = grandMa->GetPdgCode();
1167
1168 for (Int_t j=0; j<fNparents; j++){
1169 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1170 origin = kBeautyCharm;
259c3296 1171 return origin;
1172 }
1173 }
1174
1175 partMother = grandMa;
1176 } // end of iteration
1177 } // end of if
1178
1179 //gamma --------------------------
1180 else if ( abs(maPdgcode) == 22 ) {
dbe3abbe 1181 origin = kGamma;
259c3296 1182
1183 // iterate until you find B hadron as a mother or become top ancester
1184 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
1185
1186 Int_t jLabel = partMother->GetFirstMother();
1187 if (jLabel == -1){
dbe3abbe 1188 origin = kGamma;
259c3296 1189 return origin;
1190 }
1191 if (jLabel < 0){ // safety protection
1192 AliDebug(1, "Stack label is negative, return\n");
1193 return -1;
1194 }
1195
1196 // if there is an ancester
1197 TParticle* grandMa = fStack->Particle(jLabel);
1198 Int_t grandMaPDG = grandMa->GetPdgCode();
1199
1200 for (Int_t j=0; j<fNparents; j++){
1201 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1202 origin = kBeautyGamma;
259c3296 1203 return origin;
1204 }
1205 }
1206
1207 partMother = grandMa;
1208 } // end of iteration
1209
1210 return origin;
1211 } // end of if
1212
1213 //pi0 --------------------------
1214 else if ( abs(maPdgcode) == 111 ) {
dbe3abbe 1215 origin = kPi0;
259c3296 1216
1217 // iterate until you find B hadron as a mother or become top ancester
1218 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
1219
1220 Int_t jLabel = partMother->GetFirstMother();
1221 if (jLabel == -1){
dbe3abbe 1222 origin = kPi0;
259c3296 1223 return origin;
1224 }
1225 if (jLabel < 0){ // safety protection
1226 AliDebug(1, "Stack label is negative, return\n");
1227 return -1;
1228 }
1229
1230 // if there is an ancester
1231 TParticle* grandMa = fStack->Particle(jLabel);
1232 Int_t grandMaPDG = grandMa->GetPdgCode();
1233
1234 for (Int_t j=0; j<fNparents; j++){
1235 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1236 origin = kBeautyPi0;
259c3296 1237 return origin;
1238 }
1239 }
1240
1241 partMother = grandMa;
1242 } // end of iteration
1243
1244 return origin;
1245 } // end of if
1246
1247
1248 //else --------------------------
1249 else {
dbe3abbe 1250 origin = kElse;
259c3296 1251
1252 // iterate until you find B hadron as a mother or become top ancester
1253 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
1254
1255 Int_t jLabel = partMother->GetFirstMother();
1256 if (jLabel == -1){
dbe3abbe 1257 origin = kElse;
259c3296 1258 return origin;
1259 }
1260 if (jLabel < 0){ // safety protection
1261 AliDebug(1, "Stack label is negative, return\n");
1262 return -1;
1263 }
1264
1265 // if there is an ancester
1266 TParticle* grandMa = fStack->Particle(jLabel);
1267 Int_t grandMaPDG = grandMa->GetPdgCode();
1268
1269 for (Int_t j=0; j<fNparents; j++){
1270 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1271 origin = kBeautyElse;
259c3296 1272 return origin;
1273 }
1274 }
1275
1276 partMother = grandMa;
1277 } // end of iteration
1278
1279 }
1280
1281 return origin;
1282
1283}
dbe3abbe 1284*/
259c3296 1285
1286//_______________________________________________________________________________________________
1287Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track)
1288{
dbe3abbe 1289 // test cuts
259c3296 1290
1291 //if (track->Pt() < 1.0) return kFALSE;
1292 //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
1293 //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
1294 //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1295 if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
1296 //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
1297
1298
1299/*
1300 Float_t dcaR=-1;
1301 Float_t dcaZ=-1;
1302 track->GetImpactParameters(dcaR,dcaZ);
1303 if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
1304 if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
1305*/
1306 return kTRUE;
1307}
41d0c3b3 1308