]> 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
correcting compilation warning
[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";
78ea5ef4 289 fHistTagged.fPtBeautyElec= new TH1F(hname,hname,250,0,50);
259c3296 290 hname=hnopt+"pt_taggedelec";
78ea5ef4 291 fHistTagged.fPtTaggedElec= new TH1F(hname,hname,250,0,50);
259c3296 292 hname=hnopt+"pt_righttaggedelec";
78ea5ef4 293 fHistTagged.fPtRightTaggedElec = new TH1F(hname,hname,250,0,50);
259c3296 294 hname=hnopt+"pt_wrongtaggedelec";
78ea5ef4 295 fHistTagged.fPtWrongTaggedElec = new TH1F(hname,hname,250,0,50);
259c3296 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
0792aa82 822 if (track1->GetLabel()<0 || track2->GetLabel()<0) return 0;
259c3296 823 TParticle* part1 = fStack->Particle(TMath::Abs(track1->GetLabel()));
824 TParticle* part2 = fStack->Particle(TMath::Abs(track2->GetLabel()));
825 if (!(part1)) return 0;
826 if (!(part2)) return 0;
827
828 TParticle* part1Crtgen = part1; // copy track into current generation particle
829 TParticle* part2Crtgen = part2; // copy track into current generation particle
830
831
832 Int_t sourcePDG = 0;
833
834 // if the two tracks' mother's label is same, get the mother info
835 // in case of charm, check if it originated from beauty
836 for (Int_t i=0; i<100; i++){ // iterate 100
837 Int_t iLabel = part1Crtgen->GetFirstMother(); //label of mother of current generation for 1st partilce
838 if (iLabel < 0) return 0;
839
840 for (Int_t j=0; j<100; j++){ // iterate 100
841 Int_t jLabel = part2Crtgen->GetFirstMother(); //label of mother of current generation for 2nd partilce
842 if (jLabel < 0) return 0; // if jLabel == -1
843
844 if (iLabel == jLabel){ // check if two tracks are originated from same mother
845 TParticle* thismother = fStack->Particle(jLabel); // if yes, get "thismother" info
846 sourcePDG = abs(thismother->GetPdgCode()); // get the pdg code of "this mother"
847
848 // check ancester to see if it is originally from beauty
849 for (Int_t k=0; k<10; k++){ // check up to 10 ancesters
850 Int_t ancesterLabel = thismother->GetFirstMother();
851 if (ancesterLabel < 0) return sourcePDG; // if "thismoter" doesn't have mother anymore, return thismother's pdg
852
853 TParticle* thisancester = fStack->Particle(ancesterLabel);
854 Int_t ancesterPDG = abs(thisancester->GetPdgCode());
855
856 for (Int_t l=0; l<fNparents; l++){
857 if (abs(ancesterPDG)==fParentSelect[1][l]){
858 sourcePDG = -1*sourcePDG; // multiply -1 for charm from bottom
859 return sourcePDG;
860 }
861 }
862 thismother = thisancester;
863 }
864
865 }
866 part2Crtgen = fStack->Particle(jLabel); // if their mother is different, go up to previous generation of 2nd particle
867 }
868 part1Crtgen = fStack->Particle(iLabel); // if their mother is different, go up to previous generation of 2nd particle
869
870 // if you don't find additionional(2nd particle) track originated from a given beauty hadron, break to save time
871 Int_t motherPDGtmp = abs(part1Crtgen->GetPdgCode());
872 for (Int_t l=0; l<fNparents; l++){
873 if (abs(motherPDGtmp)==fParentSelect[1][l]) return sourcePDG;
874 }
875
876 }
877
878 return sourcePDG;
879
880}
881
882//_______________________________________________________________________________________________
883Int_t AliHFEsecVtx::PairCode(AliESDtrack* track1,AliESDtrack* track2)
884{
885
886 //
dbe3abbe 887 // return pair code which is predefinded as:
888 // kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0, kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
259c3296 889 //
890
891 Int_t pairOriginsPDG = PairOrigin(track1,track2);
892
dbe3abbe 893 Int_t sourcePart = kElse;
259c3296 894
895 if (pairOriginsPDG < 0) {
dbe3abbe 896 sourcePart = kBeautyElse;
259c3296 897 }
898 for (Int_t i=0; i<fNparents; i++){
899 if (abs(pairOriginsPDG)==fParentSelect[0][i]) {
dbe3abbe 900 if (pairOriginsPDG>0) sourcePart = kDirectCharm;
259c3296 901 if (pairOriginsPDG<0) {
dbe3abbe 902 sourcePart = kBeautyCharm;
259c3296 903 }
904 }
905 if (abs(pairOriginsPDG)==fParentSelect[1][i]) {
906 if (pairOriginsPDG>0) {
dbe3abbe 907 sourcePart = kDirectBeauty;
259c3296 908 }
dbe3abbe 909 if (pairOriginsPDG<0) return kElse;
259c3296 910 }
911 }
dbe3abbe 912 if (pairOriginsPDG == 22) sourcePart = kGamma;
259c3296 913 if (pairOriginsPDG == -22) {
dbe3abbe 914 sourcePart = kBeautyGamma;
259c3296 915 }
dbe3abbe 916 if (pairOriginsPDG == 111) sourcePart = kPi0;
259c3296 917 if (pairOriginsPDG == -111) {
dbe3abbe 918 sourcePart = kBeautyPi0;
259c3296 919 }
920
921 return sourcePart;
922
923}
924
dbe3abbe 925//_______________________________________________________________________________________________
926Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack)
927{
928
929 // return decay electron's origin
930
931 if (iTrack < 0) {
932 AliDebug(1, "Stack label is negative, return\n");
933 return -1;
934 }
935
936 TParticle* mcpart = fStack->Particle(iTrack);
937
938 if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !!!!!!!!!!!!!!!!!
939
940 Int_t iLabel = mcpart->GetFirstMother();
941 if (iLabel<0){
942 AliDebug(1, "Stack label is negative, return\n");
943 return -1;
944 }
945
946 Int_t origin = -1;
947 Bool_t isFinalOpenCharm = kFALSE;
948
949 TParticle *partMother = fStack->Particle(iLabel);
950 Int_t maPdgcode = partMother->GetPdgCode();
951
952 // if the mother is charmed hadron
953 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
954
955 for (Int_t i=0; i<fNparents; i++){
956 if (abs(maPdgcode)==fParentSelect[0][i]){
957 isFinalOpenCharm = kTRUE;
958 }
959 }
960 if (!isFinalOpenCharm) return -1;
961
962
963 // iterate until you find B hadron as a mother or become top ancester
964 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
965
966 Int_t jLabel = partMother->GetFirstMother();
967 if (jLabel == -1){
968 origin = kDirectCharm;
969 return origin;
970 }
971 if (jLabel < 0){ // safety protection
972 AliDebug(1, "Stack label is negative, return\n");
973 return -1;
974 }
975
976 // if there is an ancester
977 TParticle* grandMa = fStack->Particle(jLabel);
978 Int_t grandMaPDG = grandMa->GetPdgCode();
979
41d0c3b3 980 for (Int_t j=0; j<fNparents; j++){
981 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 982
983 origin = kBeautyCharm;
984 return origin;
985 }
986 }
987
988 partMother = grandMa;
989 } // end of iteration
990 } // end of if
991 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
992 for (Int_t i=0; i<fNparents; i++){
993 if (abs(maPdgcode)==fParentSelect[1][i]){
994 origin = kDirectBeauty;
995 return origin;
996 }
997 }
998 } // end of if
999
1000
1001 //============ gamma ================
1002 else if ( abs(maPdgcode) == 22 ) {
1003 origin = kGamma;
1004
1005 // iterate until you find B hadron as a mother or become top ancester
1006 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1007
1008 Int_t jLabel = partMother->GetFirstMother();
1009 if (jLabel == -1){
1010 origin = kGamma;
1011 return origin;
1012 }
1013 if (jLabel < 0){ // safety protection
1014 AliDebug(1, "Stack label is negative, return\n");
1015 return -1;
1016 }
1017
1018 // if there is an ancester
1019 TParticle* grandMa = fStack->Particle(jLabel);
1020 Int_t grandMaPDG = grandMa->GetPdgCode();
1021
41d0c3b3 1022 for (Int_t j=0; j<fNparents; j++){
1023 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1024 origin = kBeautyGamma;
1025 return origin;
1026 }
1027 }
1028
1029 partMother = grandMa;
1030 } // end of iteration
1031
1032 return origin;
1033 } // end of if
1034
1035 //============ pi0 ================
1036 else if ( abs(maPdgcode) == 111 ) {
1037 origin = kPi0;
1038
1039 // iterate until you find B hadron as a mother or become top ancester
1040 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1041
1042 Int_t jLabel = partMother->GetFirstMother();
1043 if (jLabel == -1){
1044 origin = kPi0;
1045 return origin;
1046 }
1047 if (jLabel < 0){ // safety protection
1048 AliDebug(1, "Stack label is negative, return\n");
1049 return -1;
1050 }
1051
1052 // if there is an ancester
1053 TParticle* grandMa = fStack->Particle(jLabel);
1054 Int_t grandMaPDG = grandMa->GetPdgCode();
1055
41d0c3b3 1056 for (Int_t j=0; j<fNparents; j++){
1057 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1058 origin = kBeautyPi0;
1059 return origin;
1060 }
1061 }
1062
1063 partMother = grandMa;
1064 } // end of iteration
1065
1066 return origin;
1067 } // end of if
1068
1069
1070 else {
1071 origin = kElse;
1072
1073 // iterate until you find B hadron as a mother or become top ancester
1074 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
1075
1076 Int_t jLabel = partMother->GetFirstMother();
1077 if (jLabel == -1){
1078 origin = kElse;
1079 return origin;
1080 }
1081 if (jLabel < 0){ // safety protection
1082 AliDebug(1, "Stack label is negative, return\n");
1083 return -1;
1084 }
1085
1086 // if there is an ancester
1087 TParticle* grandMa = fStack->Particle(jLabel);
1088 Int_t grandMaPDG = grandMa->GetPdgCode();
1089
41d0c3b3 1090 for (Int_t j=0; j<fNparents; j++){
1091 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1092 origin = kBeautyElse;
1093 return origin;
1094 }
1095 }
1096
1097 partMother = grandMa;
1098 } // end of iteration
1099
1100 }
1101
1102 return origin;
1103
1104}
1105
1106/*
259c3296 1107//_______________________________________________________________________________________________
1108Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
1109{
1110
1111 //
1112 // decay electron's origin
1113 //
1114
1115 if (iTrackLabel < 0) {
1116 AliDebug(1, "Stack label is negative, return\n");
1117 return -1;
1118 }
1119
1120 TParticle* mcpart = fStack->Particle(iTrackLabel);
1121 Int_t iLabel = mcpart->GetFirstMother();
1122 if (iLabel<0){
1123 AliDebug(1, "Stack label is negative, return\n");
1124 return -1;
1125 }
1126
1127 Int_t origin = -1;
dbe3abbe 1128 Bool_t isFinalOpenCharm = kFALSE;
259c3296 1129
1130 TParticle *partMother = fStack->Particle(iLabel);
1131 Int_t maPdgcode = partMother->GetPdgCode(); // get mother's pdg code
1132
1133 //beauty --------------------------
dbe3abbe 1134 if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
259c3296 1135 for (Int_t i=0; i<fNparents; i++){
1136 if (abs(maPdgcode)==fParentSelect[1][i]){
dbe3abbe 1137 origin = kDirectBeauty;
259c3296 1138 return origin;
1139 }
1140 else return -1; // this track is originated beauties not in the final B hadron list => excited beauties
1141 }
1142 } // end of if
1143
1144 //charm --------------------------
dbe3abbe 1145 else if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
259c3296 1146
1147 for (Int_t i=0; i<fNparents; i++){
1148 if (abs(maPdgcode)==fParentSelect[0][i])
dbe3abbe 1149 isFinalOpenCharm = kTRUE;
259c3296 1150 }
dbe3abbe 1151 if (!isFinalOpenCharm) return -1; // this track is originated charms not in the final D hadron list => excited charms
259c3296 1152 // to prevent any possible double counting
1153
1154 for (Int_t i=0; i<100; i++){ // iterate 100 until you find B hadron as a mother or become top ancester
1155
1156 Int_t jLabel = partMother->GetFirstMother();
1157 if (jLabel == -1){
dbe3abbe 1158 origin = kDirectCharm;
259c3296 1159 return origin;
1160 }
1161 if (jLabel < 0){ // safety protection even though not really necessary here
1162 AliDebug(1, "Stack label is negative, return\n");
1163 return -1;
1164 }
1165
1166 // if there is an ancester, check if it in the final B hadron list
1167 TParticle* grandMa = fStack->Particle(jLabel);
1168 Int_t grandMaPDG = grandMa->GetPdgCode();
1169
1170 for (Int_t j=0; j<fNparents; j++){
1171 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1172 origin = kBeautyCharm;
259c3296 1173 return origin;
1174 }
1175 }
1176
1177 partMother = grandMa;
1178 } // end of iteration
1179 } // end of if
1180
1181 //gamma --------------------------
1182 else if ( abs(maPdgcode) == 22 ) {
dbe3abbe 1183 origin = kGamma;
259c3296 1184
1185 // iterate until you find B hadron as a mother or become top ancester
1186 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
1187
1188 Int_t jLabel = partMother->GetFirstMother();
1189 if (jLabel == -1){
dbe3abbe 1190 origin = kGamma;
259c3296 1191 return origin;
1192 }
1193 if (jLabel < 0){ // safety protection
1194 AliDebug(1, "Stack label is negative, return\n");
1195 return -1;
1196 }
1197
1198 // if there is an ancester
1199 TParticle* grandMa = fStack->Particle(jLabel);
1200 Int_t grandMaPDG = grandMa->GetPdgCode();
1201
1202 for (Int_t j=0; j<fNparents; j++){
1203 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1204 origin = kBeautyGamma;
259c3296 1205 return origin;
1206 }
1207 }
1208
1209 partMother = grandMa;
1210 } // end of iteration
1211
1212 return origin;
1213 } // end of if
1214
1215 //pi0 --------------------------
1216 else if ( abs(maPdgcode) == 111 ) {
dbe3abbe 1217 origin = kPi0;
259c3296 1218
1219 // iterate until you find B hadron as a mother or become top ancester
1220 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
1221
1222 Int_t jLabel = partMother->GetFirstMother();
1223 if (jLabel == -1){
dbe3abbe 1224 origin = kPi0;
259c3296 1225 return origin;
1226 }
1227 if (jLabel < 0){ // safety protection
1228 AliDebug(1, "Stack label is negative, return\n");
1229 return -1;
1230 }
1231
1232 // if there is an ancester
1233 TParticle* grandMa = fStack->Particle(jLabel);
1234 Int_t grandMaPDG = grandMa->GetPdgCode();
1235
1236 for (Int_t j=0; j<fNparents; j++){
1237 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1238 origin = kBeautyPi0;
259c3296 1239 return origin;
1240 }
1241 }
1242
1243 partMother = grandMa;
1244 } // end of iteration
1245
1246 return origin;
1247 } // end of if
1248
1249
1250 //else --------------------------
1251 else {
dbe3abbe 1252 origin = kElse;
259c3296 1253
1254 // iterate until you find B hadron as a mother or become top ancester
1255 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
1256
1257 Int_t jLabel = partMother->GetFirstMother();
1258 if (jLabel == -1){
dbe3abbe 1259 origin = kElse;
259c3296 1260 return origin;
1261 }
1262 if (jLabel < 0){ // safety protection
1263 AliDebug(1, "Stack label is negative, return\n");
1264 return -1;
1265 }
1266
1267 // if there is an ancester
1268 TParticle* grandMa = fStack->Particle(jLabel);
1269 Int_t grandMaPDG = grandMa->GetPdgCode();
1270
1271 for (Int_t j=0; j<fNparents; j++){
1272 if (abs(grandMaPDG)==fParentSelect[1][j]){
dbe3abbe 1273 origin = kBeautyElse;
259c3296 1274 return origin;
1275 }
1276 }
1277
1278 partMother = grandMa;
1279 } // end of iteration
1280
1281 }
1282
1283 return origin;
1284
1285}
dbe3abbe 1286*/
259c3296 1287
1288//_______________________________________________________________________________________________
75d81601 1289Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track) const
259c3296 1290{
dbe3abbe 1291 // test cuts
259c3296 1292
1293 //if (track->Pt() < 1.0) return kFALSE;
1294 //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
1295 //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
1296 //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1297 if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
1298 //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
1299
1300
1301/*
1302 Float_t dcaR=-1;
1303 Float_t dcaZ=-1;
1304 track->GetImpactParameters(dcaR,dcaZ);
1305 if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
1306 if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
1307*/
1308 return kTRUE;
1309}