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