]> 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
Fixes for the analysis modes (Marek)
[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 *
18 * *
19 * Authors: *
20 * MinJung Kweon <minjung@physi.uni-heidelberg.de> *
21 * *
22 **************************************************************************/
23
24#include <iostream>
25#include <TH2F.h>
26#include <TH3F.h>
27#include <TLorentzVector.h>
28#include <TNtuple.h>
29#include <TParticle.h>
30
31#include <AliESDEvent.h>
32#include <AliESDtrack.h>
33#include "AliHFEsecVtx.h"
34#include <AliKFParticle.h>
35#include <AliKFVertex.h>
36#include <AliLog.h>
37#include <AliPID.h>
38#include <AliStack.h>
39
40ClassImp(AliHFEsecVtx)
41
42//_______________________________________________________________________________________________
43AliHFEsecVtx::AliHFEsecVtx():
44 fESD1(0x0)
45 ,fStack(0x0)
46 ,fNparents(0)
47 ,fHistTagged()
48 ,fPairTagged(0)
49 ,fdistTwoSecVtx(-1)
50 ,fcosPhi(-1)
51 ,fsignedLxy(-1)
52 ,finvmass(-1)
53 ,finvmassSigma(-1)
54 ,fBTagged(0)
55{
56 // Default constructor
57
58 Init();
59
60}
61
62//_______________________________________________________________________________________________
63AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
64 TObject(p)
65 ,fESD1(0x0)
66 ,fStack(0x0)
67 ,fNparents(p.fNparents)
68 ,fHistTagged()
69 ,fPairTagged(p.fPairTagged)
70 ,fdistTwoSecVtx(p.fdistTwoSecVtx)
71 ,fcosPhi(p.fcosPhi)
72 ,fsignedLxy(p.fsignedLxy)
73 ,finvmass(p.finvmass)
74 ,finvmassSigma(p.finvmassSigma)
75 ,fBTagged(p.fBTagged)
76{
77 // Copy constructor
78}
79
80//_______________________________________________________________________________________________
81AliHFEsecVtx&
82AliHFEsecVtx::operator=(const AliHFEsecVtx &)
83{
84 // Assignment operator
85
86 AliInfo("Not yet implemented.");
87 return *this;
88}
89
90//_______________________________________________________________________________________________
91AliHFEsecVtx::~AliHFEsecVtx()
92{
93 // Destructor
94
95 //cout << "Analysis Done." << endl;
96}
97
98//__________________________________________
99void AliHFEsecVtx::Init()
100{
101
102 fNparents = 7;
103
104 fParentSelect[0][0] = 411;
105 fParentSelect[0][1] = 421;
106 fParentSelect[0][2] = 431;
107 fParentSelect[0][3] = 4122;
108 fParentSelect[0][4] = 4132;
109 fParentSelect[0][5] = 4232;
110 fParentSelect[0][6] = 4332;
111
112 fParentSelect[1][0] = 511;
113 fParentSelect[1][1] = 521;
114 fParentSelect[1][2] = 531;
115 fParentSelect[1][3] = 5122;
116 fParentSelect[1][4] = 5132;
117 fParentSelect[1][5] = 5232;
118 fParentSelect[1][6] = 5332;
119
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;
163
164
165 //fid[4][3] = {0,1,2, 0,1,3, 0,2,3, 1,2,3};
166 //fia[4][3][2] = {0,1, 0,2, 1,2, 0,1, 0,3, 1,3, 0,2, 0,3, 2,3, 1,2, 1,3, 2,3};
167
168}
169
170//__________________________________________
171void AliHFEsecVtx::ResetTagVar()
172{
173 fdistTwoSecVtx = -1;
174 fcosPhi = -1;
175 fsignedLxy = -1;
176 finvmass = -1;
177 finvmassSigma = -1;
178 fBTagged = kFALSE;
179}
180
181//__________________________________________
182void AliHFEsecVtx::InitAnaPair()
183{
184 fPairTagged = 0;
185 for (Int_t i=0; i<20; i++){
186 fpairedTrackID[i] = -1;
187 fpairedChi2[i] = -1;
188 fpairedInvMass[i] = -1;
189 fpairedSignedLxy[i] = -1;
190 }
191
192}
193
194//_______________________________________________________________________________________________
195void AliHFEsecVtx::CreateHistograms(TString hnopt)
196{
197 // create histograms
198
199 fkSourceLabel[fkAll]="all";
200 fkSourceLabel[fkDirectCharm]="directCharm";
201 fkSourceLabel[fkDirectBeauty]="directBeauty";
202 fkSourceLabel[fkBeautyCharm]="beauty2charm";
203 fkSourceLabel[fkGamma]="gamma";
204 fkSourceLabel[fkPi0]="pi0";
205 fkSourceLabel[fkElse]="others";
206 fkSourceLabel[fkBeautyGamma]="beauty22gamma";
207 fkSourceLabel[fkBeautyPi0]="beauty22pi0";
208 fkSourceLabel[fkBeautyElse]="beauty22others";
209
210
211 TString hname;
212 for (Int_t isource = 0; isource < 10; isource++ ){
213
214 hname=hnopt+"InvMass_"+fkSourceLabel[isource];
215 fHistPair[isource].fInvMass = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
216 hname=hnopt+"InvMassCut1_"+fkSourceLabel[isource];
217 fHistPair[isource].fInvMassCut1 = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
218 hname=hnopt+"InvMassCut2_"+fkSourceLabel[isource];
219 fHistPair[isource].fInvMassCut2 = new TH2F(hname,hname+";invMass;invMassSigma",120,-2,10,100,0,10);
220 hname=hnopt+"KFChi2_"+fkSourceLabel[isource];
221 fHistPair[isource].fKFChi2 = new TH1F(hname,hname,200,0,20);
222 hname=hnopt+"CosOpenAngle_"+fkSourceLabel[isource];
223 fHistPair[isource].fCosOpenAngle = new TH1F(hname,hname,100,-1.1,1.1);
224 hname=hnopt+"SignedLxy_"+fkSourceLabel[isource];
225 fHistPair[isource].fSignedLxy = new TH2F(hname,hname,1000,-5,5,120,-2,10);
226 hname=hnopt+"KFIP_"+fkSourceLabel[isource];
227 fHistPair[isource].fKFIP = new TH1F(hname,hname,1000,-5,5);
228
229 }
230
231 hname=hnopt+"pt_beautyelec";
232 fHistTagged.fPtBeautyElec= new TH1F(hname,hname,150,0,30);
233 hname=hnopt+"pt_taggedelec";
234 fHistTagged.fPtTaggedElec= new TH1F(hname,hname,150,0,30);
235 hname=hnopt+"pt_righttaggedelec";
236 fHistTagged.fPtRightTaggedElec = new TH1F(hname,hname,150,0,30);
237 hname=hnopt+"pt_wrongtaggedelec";
238 fHistTagged.fPtWrongTaggedElec = new TH1F(hname,hname,150,0,30);
239
240}
241
242//_______________________________________________________________________________________________
243void AliHFEsecVtx::AnaPair(AliESDtrack* track1, AliESDtrack* track2, Int_t sourcePart, Int_t index2)
244{
245
246 // get KF particle input pid
247 Int_t pdg1 = GetMCPID(track1);
248 Int_t pdg2 = GetMCPID(track2);
249
250
251 // create KF particle of pair
252 AliKFParticle::SetField(fESD1->GetMagneticField());
253 AliKFParticle kfTrack1(*track1, pdg1);
254 AliKFParticle kfTrack2(*track2, pdg2);
255
256 AliKFParticle kfSecondary(kfTrack1,kfTrack2);
257
258 // copy primary vertex from ESD
259 AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
260 if( primVtxCopy.GetNDF() <1 ) return;
261
262 //primary vertex point
263 Double_t pvx = primVtxCopy.GetX();
264 Double_t pvy = primVtxCopy.GetY();
265 //Double_t pvz = primVtxCopy.GetZ();
266
267 //secondary vertex point from kf particle
268 Double_t kfx = kfSecondary.GetX();
269 Double_t kfy = kfSecondary.GetY();
270 //Double_t kfz = kfSecondary.GetZ();
271
272 //momentum at the decay point from kf particle
273 Double_t kfpx = kfSecondary.GetPx();
274 Double_t kfpy = kfSecondary.GetPy();
275 //Double_t kfpz = kfSecondary.GetPz();
276
277
278 Double_t dx = kfx-pvx;
279 Double_t dy = kfy-pvy;
280
281
282
283 // discriminating variables ----------------------------------------------------------
284
285 // invariant mass of the KF particle
286 Double_t invmass = -1;
287 Double_t invmassSigma = -1;
288 kfSecondary.GetMass(invmass,invmassSigma);
289
290 // chi2 of the KF particle
291 Double_t kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
292
293 // opening angle between two particles in XY plane
294 Double_t cosphi = TMath::Cos(kfTrack1.GetAngleXY(kfTrack2));
295
296 // projection of kf vertex vector to the kf momentum direction
297 Double_t costheta = ( dx*kfpx + dy*kfpy)/TMath::Sqrt(dx*dx+dy*dy)*TMath::Sqrt(kfpx*kfpx + kfpy*kfpy);
298 Double_t signedLxy = TMath::Sqrt(dx*dx+dy*dy)*costheta;
299
300 // DCA from primary to e-h KF particle (impact parameter of KF particle)
301 Double_t vtx[2]={pvx, pvy};
302 Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
303
304
305
306
307 // pair cuts
308 if( kfchi2 >2. ) return;
309
310 // fill histograms
311 fHistPair[sourcePart].fInvMass->Fill(invmass,invmassSigma);
312 fHistPair[sourcePart].fKFChi2->Fill(kfchi2);
313 fHistPair[sourcePart].fCosOpenAngle->Fill(cosphi);
314 fHistPair[sourcePart].fSignedLxy->Fill(signedLxy,invmass);
315 fHistPair[sourcePart].fKFIP->Fill(kfip);
316
317 if ( signedLxy > 0.05 && cosphi > 0.5) {
318 fHistPair[sourcePart].fInvMassCut1->Fill(invmass,invmassSigma);
319 }
320
321 // pair cuts
322 if (signedLxy < 0.05) return;
323 if (cosphi < 0.0) return;
324
325 // pair tagging if it passed the above cuts
326 fHistPair[sourcePart].fInvMassCut2->Fill(invmass,invmassSigma);
327
328
329 // pair tagging condition
330 if ( signedLxy > 0.0 && cosphi > 0) {
331 //if ( signedLxy > 0.06 && cosphi > 0) {
332 fpairedTrackID[fPairTagged] = index2;
333 fpairedChi2[fPairTagged] = kfchi2;
334 fpairedInvMass[fPairTagged] = invmass;
335 fpairedSignedLxy[fPairTagged] = signedLxy;
336 fPairTagged++;
337 }
338
339}
340
341//_______________________________________________________________________________________________
342void AliHFEsecVtx::RunSECVTX(AliESDtrack *track)
343{
344
345 ResetTagVar();
346
347 if (fPairTagged >= 4) {
348 FindSECVTXCandid4Tracks(track);
349 }
350 else if (fPairTagged == 3) {
351 FindSECVTXCandid3Tracks(track);
352 }
353 else if (fPairTagged == 2) {
354 FindSECVTXCandid2Tracks(track);
355 }
356 else if (fPairTagged == 1) {
357 ApplyPairTagCut();
358 }
359
360
361 Int_t imclabel = TMath::Abs(track->GetLabel());
362 if(imclabel<0) return;
363 TParticle* mcpart = fStack->Particle(imclabel);
364 Int_t esource = GetElectronSource(imclabel);
365 if (esource == fkDirectBeauty || esource == fkBeautyCharm || esource == fkBeautyGamma || esource == fkBeautyPi0 || esource == fkBeautyElse){
366 fHistTagged.fPtBeautyElec->Fill(mcpart->Pt());
367 }
368
369
370 if (fBTagged) {
371 fHistTagged.fPtTaggedElec->Fill(mcpart->Pt());
372 if (esource == fkDirectBeauty || esource == fkBeautyCharm || esource == fkBeautyGamma || esource == fkBeautyPi0 || esource == fkBeautyElse){
373 fHistTagged.fPtRightTaggedElec->Fill(mcpart->Pt());
374 }
375 else fHistTagged.fPtWrongTaggedElec->Fill(mcpart->Pt());
376 }
377
378}
379
380//_______________________________________________________________________________________________
381void AliHFEsecVtx::ApplyPairTagCut()
382{
383
384 if (fpairedChi2[0] > 2.0) return;
385 if (fpairedInvMass[0] < 1.5) return;
386 if (fpairedSignedLxy[0] < 0.5) return;
387
388 fBTagged = kTRUE;
389
390}
391
392//_______________________________________________________________________________________________
393Bool_t AliHFEsecVtx::ApplyPairTagCut(Int_t id)
394{
395
396 if (fpairedChi2[id] > 2.0) return kFALSE;
397 if (fpairedInvMass[id] < 1.5) return kFALSE;
398 if (fpairedSignedLxy[id] < 0.5) return kFALSE;
399
400 fBTagged = kTRUE;
401 return kTRUE;
402
403}
404
405//_______________________________________________________________________________________________
406Bool_t AliHFEsecVtx::ApplyTagCut(Double_t chi2)
407{
408
409 if (chi2 > 2.0) return kFALSE;
410 if (finvmass < 1.5) return kFALSE;
411 if (fsignedLxy < 0.5) return kFALSE;
412 if (fcosPhi < 0.90) return kFALSE;
413 if (fdistTwoSecVtx > 0.1) return kFALSE;
414
415 fBTagged = kTRUE;
416 return kTRUE;
417}
418
419//_______________________________________________________________________________________________
420//void AliHFEsecVtx::ApplyTagCut(Double_t chi2, AliESDtrack *track, AliESDtrack *htrack1, AliESDtrack *htrack2)
421void AliHFEsecVtx::ApplyVtxTagCut(Double_t chi2)
422{
423
424 if(!ApplyTagCut(chi2)){
425 if(!ApplyPairTagCut(0)) ApplyPairTagCut(1);
426 }
427
428}
429
430//_______________________________________________________________________________________________
431void AliHFEsecVtx::FindSECVTXCandid4Tracks(AliESDtrack *track)
432{
433
434 // sort pair in increasing order (kFALSE-increasing order)
435 Int_t index[20];
436 Int_t indexA[4];
437 Int_t indexB[3];
438 Double_t sevchi2[4];
439 AliESDtrack *htrack[4];
440
441 TMath::Sort(fPairTagged,fpairedChi2,index,kFALSE);
442
443 // select 4 partner tracks retruning smallest pair chi2
444
445 for (Int_t i=0; i<4; i++){
446 htrack[i] = fESD1->GetTrack(fpairedTrackID[index[i]]);
447 }
448
449 // calculated chi2 with 1 electron and 3 partner tracks
450 for (Int_t i=0; i<4; i++){
451 sevchi2[i] = GetSecVtxChi2(track, htrack[fid[i][0]], htrack[fid[i][1]], htrack[fid[i][2]]);
452 }
453
454 // select 3 partner tracks retruning smallest pair chi2
455 // [think] if two smallest chi2 are similar, have to think about better handling of selection
456 TMath::Sort(4,sevchi2,indexA,kFALSE);
457
458 // calculated chi2 with 1 electron and 2 partner tracks
459 for (Int_t i=0; i<3; i++){
460 sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
461 }
462
463 // select 2 partner tracks retruning smallest pair chi2
464 TMath::Sort(3,sevchi2,indexB,kFALSE);
465
466 // calculate secondary vertex quality variables with 1 electron and 2 hadrons
467 CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
468
469 ApplyVtxTagCut(sevchi2[indexB[0]]);
470 //ApplyTagCut(sevchi2[indexB[0]],track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
471
472}
473
474//_______________________________________________________________________________________________
475void AliHFEsecVtx::FindSECVTXCandid3Tracks(AliESDtrack *track)
476{
477
478 // sort pair in increasing order (kFALSE-increasing order)
479 Int_t indexA[1] = { 0 };
480 Int_t indexB[3];
481 Double_t sevchi2[3];
482 AliESDtrack *htrack[3];
483
484 // select 4 partner tracks retruning smallest pair chi2
485
486 for (Int_t i=0; i<3; i++){
487 htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
488 }
489
490 // calculated chi2 with 1 electron and 2 partner tracks
491 for (Int_t i=0; i<3; i++){
492 sevchi2[i] = GetSecVtxChi2(track, htrack[fia[indexA[0]][i][0]], htrack[fia[indexA[0]][i][1]]);
493 }
494
495 // select 2 partner tracks retruning smallest pair chi2
496 TMath::Sort(3,sevchi2,indexB,kFALSE);
497
498 // calculate secondary vertex quality variables with 1 electron and 2 hadrons
499 CalcSECVTXProperty(track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
500
501 ApplyVtxTagCut(sevchi2[indexB[0]]);
502 //ApplyTagCut(sevchi2[indexB[0]],track,htrack[fia[indexA[0]][indexB[0]][0]],htrack[fia[indexA[0]][indexB[0]][1]]);
503}
504
505//_______________________________________________________________________________________________
506void AliHFEsecVtx::FindSECVTXCandid2Tracks(AliESDtrack *track)
507{
508
509 Double_t sevchi2[1];
510 AliESDtrack *htrack[2];
511
512 for (Int_t i=0; i<2; i++){
513 htrack[i] = fESD1->GetTrack(fpairedTrackID[i]);
514 }
515
516 sevchi2[0] = GetSecVtxChi2(track, htrack[0], htrack[1]);
517
518 // calculate secondary vertex quality variables with 1 electron and 2 hadrons
519 CalcSECVTXProperty(track,htrack[0],htrack[1]);
520
521 ApplyVtxTagCut(sevchi2[0]);
522 //ApplyTagCut(sevchi2[0],track,htrack[0],htrack[1]);
523}
524
525//_______________________________________________________________________________________________
526void AliHFEsecVtx::CalcSECVTXProperty(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
527{
528
529 Int_t pdg1 = GetMCPID(track1);
530 Int_t pdg2 = GetMCPID(track2);
531 Int_t pdg3 = GetMCPID(track3);
532
533 AliKFParticle::SetField(fESD1->GetMagneticField());
534 AliKFParticle kfTrack1(*track1, pdg1);
535 AliKFParticle kfTrack2(*track2, pdg2);
536 AliKFParticle kfTrack3(*track3, pdg3);
537
538 AliKFParticle kfSecondary12(kfTrack1,kfTrack2);
539 AliKFParticle kfSecondary13(kfTrack1,kfTrack3);
540 AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
541
542 // copy primary vertex from ESD
543 AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
544 //printf("primary ndf= %d\n",primVtxCopy.GetNDF());
545 if( primVtxCopy.GetNDF() <1 ) return;
546
547 Double_t kdx12 = kfSecondary12.GetX()-primVtxCopy.GetX();
548 Double_t kdy12 = kfSecondary12.GetY()-primVtxCopy.GetY();
549 //Double_t kdz12 = kfSecondary12.GetZ()-primVtxCopy.GetZ();
550
551 Double_t kdx13 = kfSecondary13.GetX()-primVtxCopy.GetX();
552 Double_t kdy13 = kfSecondary13.GetY()-primVtxCopy.GetY();
553 //Double_t kdz13 = kfSecondary13.GetZ()-primVtxCopy.GetZ();
554
555 Double_t kdx = kfSecondary.GetX()-primVtxCopy.GetX();
556 Double_t kdy = kfSecondary.GetY()-primVtxCopy.GetY();
557 //Double_t kdz = kfSecondary.GetZ()-primVtxCopy.GetZ();
558
559 // calculate distance and angle between two secvtxes
560 fdistTwoSecVtx = TMath::Sqrt((kdx12-kdx13)*(kdx12-kdx13) + (kdy12-kdy13)*(kdy12-kdy13));
561 fcosPhi = ( kdx12*kdx13 + kdy12*kdy13 ) / ( TMath::Sqrt(kdx12*kdx12+kdy12*kdy12)*TMath::Sqrt(kdx13*kdx13+kdy13*kdy13) );
562 //Double_t lengthdiff = TMath::Abs(TMath::Sqrt(kdx12*kdx12+kdy12*kdy12) - TMath::Sqrt(kdx13*kdx13+kdy13*kdy13));
563
564 // calculate angle between secondary vertex vector and secondary particle momentum vector in transverse plane
565 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());
566 // calculate signed Lxy
567 fsignedLxy = TMath::Sqrt(kdx*kdx+kdy*kdy)*cosTheta;
568
569 // calculate invariant mass of the kf particle
570 kfSecondary.GetMass(finvmass,finvmassSigma);
571}
572
573//_______________________________________________________________________________________________
574Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track)
575{
576
577 Int_t label = TMath::Abs(track->GetLabel());
578 TParticle* mcpart = fStack->Particle(label);
579 if ( !mcpart ) return 0;
580 Int_t pdgCode = mcpart->GetPdgCode();
581
582 return pdgCode;
583}
584
585//_______________________________________________________________________________________________
586Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3, AliESDtrack* track4)
587{
588
589 Int_t pdg1 = GetMCPID(track1);
590 Int_t pdg2 = GetMCPID(track2);
591 Int_t pdg3 = GetMCPID(track3);
592 Int_t pdg4 = GetMCPID(track4);
593
594 AliKFParticle::SetField(fESD1->GetMagneticField());
595 AliKFParticle kfTrack1(*track1, pdg1);
596 AliKFParticle kfTrack2(*track2, pdg2);
597 AliKFParticle kfTrack3(*track3, pdg3);
598 AliKFParticle kfTrack4(*track4, pdg4);
599 AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3,kfTrack4);
600
601 return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
602
603}
604
605//_______________________________________________________________________________________________
606Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2, AliESDtrack* track3)
607{
608
609 Int_t pdg1 = GetMCPID(track1);
610 Int_t pdg2 = GetMCPID(track2);
611 Int_t pdg3 = GetMCPID(track3);
612
613 AliKFParticle::SetField(fESD1->GetMagneticField());
614 AliKFParticle kfTrack1(*track1, pdg1);
615 AliKFParticle kfTrack2(*track2, pdg2);
616 AliKFParticle kfTrack3(*track3, pdg3);
617 AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
618
619 return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
620
621}
622
623//_______________________________________________________________________________________________
624Double_t AliHFEsecVtx::GetSecVtxChi2(AliESDtrack* track1, AliESDtrack* track2)
625{
626
627 Int_t pdg1 = GetMCPID(track1);
628 Int_t pdg2 = GetMCPID(track2);
629
630 AliKFParticle::SetField(fESD1->GetMagneticField());
631 AliKFParticle kfTrack1(*track1, pdg1);
632 AliKFParticle kfTrack2(*track2, pdg2);
633 AliKFParticle kfSecondary(kfTrack1,kfTrack2);
634
635 return TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
636
637}
638
639//_______________________________________________________________________________________________
640Int_t AliHFEsecVtx::PairOrigin(AliESDtrack* track1, AliESDtrack* track2)
641{
642
643 //
644 // return pdg code of the origin(source) of the pair
645 //
646 //
647 // ---*---*---*-----ancester A----- track1
648 // |____*______
649 // |______ track2
650 // => if they originated from same ancester,
651 // then return "the absolute value of pdg code of ancester A"
652 //
653 // ---*---*---B hadron-----ancester A----- track1
654 // |____*______
655 // |______ track2
656 // => if they originated from same ancester, and this ancester originally comes from B hadrons
657 // then return -1*"the absolute value of pdg code of ancester A"
658 //
659 // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
660 //
661
662 TParticle* part1 = fStack->Particle(TMath::Abs(track1->GetLabel()));
663 TParticle* part2 = fStack->Particle(TMath::Abs(track2->GetLabel()));
664 if (!(part1)) return 0;
665 if (!(part2)) return 0;
666
667 TParticle* part1Crtgen = part1; // copy track into current generation particle
668 TParticle* part2Crtgen = part2; // copy track into current generation particle
669
670
671 Int_t sourcePDG = 0;
672
673 // if the two tracks' mother's label is same, get the mother info
674 // in case of charm, check if it originated from beauty
675 for (Int_t i=0; i<100; i++){ // iterate 100
676 Int_t iLabel = part1Crtgen->GetFirstMother(); //label of mother of current generation for 1st partilce
677 if (iLabel < 0) return 0;
678
679 for (Int_t j=0; j<100; j++){ // iterate 100
680 Int_t jLabel = part2Crtgen->GetFirstMother(); //label of mother of current generation for 2nd partilce
681 if (jLabel < 0) return 0; // if jLabel == -1
682
683 if (iLabel == jLabel){ // check if two tracks are originated from same mother
684 TParticle* thismother = fStack->Particle(jLabel); // if yes, get "thismother" info
685 sourcePDG = abs(thismother->GetPdgCode()); // get the pdg code of "this mother"
686
687 // check ancester to see if it is originally from beauty
688 for (Int_t k=0; k<10; k++){ // check up to 10 ancesters
689 Int_t ancesterLabel = thismother->GetFirstMother();
690 if (ancesterLabel < 0) return sourcePDG; // if "thismoter" doesn't have mother anymore, return thismother's pdg
691
692 TParticle* thisancester = fStack->Particle(ancesterLabel);
693 Int_t ancesterPDG = abs(thisancester->GetPdgCode());
694
695 for (Int_t l=0; l<fNparents; l++){
696 if (abs(ancesterPDG)==fParentSelect[1][l]){
697 sourcePDG = -1*sourcePDG; // multiply -1 for charm from bottom
698 return sourcePDG;
699 }
700 }
701 thismother = thisancester;
702 }
703
704 }
705 part2Crtgen = fStack->Particle(jLabel); // if their mother is different, go up to previous generation of 2nd particle
706 }
707 part1Crtgen = fStack->Particle(iLabel); // if their mother is different, go up to previous generation of 2nd particle
708
709 // if you don't find additionional(2nd particle) track originated from a given beauty hadron, break to save time
710 Int_t motherPDGtmp = abs(part1Crtgen->GetPdgCode());
711 for (Int_t l=0; l<fNparents; l++){
712 if (abs(motherPDGtmp)==fParentSelect[1][l]) return sourcePDG;
713 }
714
715 }
716
717 return sourcePDG;
718
719}
720
721//_______________________________________________________________________________________________
722Int_t AliHFEsecVtx::PairCode(AliESDtrack* track1,AliESDtrack* track2)
723{
724
725 //
726 // return pair code which is predefind as:
727 // fkDirectCharm, fkDirectBeauty, fkBeautyCharm, fkGamma, fkPi0, fkElse, fkBeautyGamma, fkBeautyPi0, fkBeautyElse
728 //
729
730 Int_t pairOriginsPDG = PairOrigin(track1,track2);
731
732 Int_t sourcePart = fkElse;
733
734 if (pairOriginsPDG < 0) {
735 sourcePart = fkBeautyElse;
736 }
737 for (Int_t i=0; i<fNparents; i++){
738 if (abs(pairOriginsPDG)==fParentSelect[0][i]) {
739 if (pairOriginsPDG>0) sourcePart = fkDirectCharm;
740 if (pairOriginsPDG<0) {
741 sourcePart = fkBeautyCharm;
742 }
743 }
744 if (abs(pairOriginsPDG)==fParentSelect[1][i]) {
745 if (pairOriginsPDG>0) {
746 sourcePart = fkDirectBeauty;
747 }
748 if (pairOriginsPDG<0) return fkElse;
749 }
750 }
751 if (pairOriginsPDG == 22) sourcePart = fkGamma;
752 if (pairOriginsPDG == -22) {
753 sourcePart = fkBeautyGamma;
754 }
755 if (pairOriginsPDG == 111) sourcePart = fkPi0;
756 if (pairOriginsPDG == -111) {
757 sourcePart = fkBeautyPi0;
758 }
759
760 return sourcePart;
761
762}
763
764//_______________________________________________________________________________________________
765Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrackLabel)
766{
767
768 //
769 // decay electron's origin
770 //
771
772 if (iTrackLabel < 0) {
773 AliDebug(1, "Stack label is negative, return\n");
774 return -1;
775 }
776
777 TParticle* mcpart = fStack->Particle(iTrackLabel);
778 Int_t iLabel = mcpart->GetFirstMother();
779 if (iLabel<0){
780 AliDebug(1, "Stack label is negative, return\n");
781 return -1;
782 }
783
784 Int_t origin = -1;
785 Bool_t IsFinalOpenCharm = kFALSE;
786
787 TParticle *partMother = fStack->Particle(iLabel);
788 Int_t maPdgcode = partMother->GetPdgCode(); // get mother's pdg code
789
790 //beauty --------------------------
791 if ( int(abs(maPdgcode)/100.) == fkBeauty || int(abs(maPdgcode)/1000.) == fkBeauty ) {
792 for (Int_t i=0; i<fNparents; i++){
793 if (abs(maPdgcode)==fParentSelect[1][i]){
794 origin = fkDirectBeauty;
795 return origin;
796 }
797 else return -1; // this track is originated beauties not in the final B hadron list => excited beauties
798 }
799 } // end of if
800
801 //charm --------------------------
802 else if ( int(abs(maPdgcode)/100.) == fkCharm || int(abs(maPdgcode)/1000.) == fkCharm ) {
803
804 for (Int_t i=0; i<fNparents; i++){
805 if (abs(maPdgcode)==fParentSelect[0][i])
806 IsFinalOpenCharm = kTRUE;
807 }
808 if (!IsFinalOpenCharm) return -1; // this track is originated charms not in the final D hadron list => excited charms
809 // to prevent any possible double counting
810
811 for (Int_t i=0; i<100; i++){ // iterate 100 until you find B hadron as a mother or become top ancester
812
813 Int_t jLabel = partMother->GetFirstMother();
814 if (jLabel == -1){
815 origin = fkDirectCharm;
816 return origin;
817 }
818 if (jLabel < 0){ // safety protection even though not really necessary here
819 AliDebug(1, "Stack label is negative, return\n");
820 return -1;
821 }
822
823 // if there is an ancester, check if it in the final B hadron list
824 TParticle* grandMa = fStack->Particle(jLabel);
825 Int_t grandMaPDG = grandMa->GetPdgCode();
826
827 for (Int_t j=0; j<fNparents; j++){
828 if (abs(grandMaPDG)==fParentSelect[1][j]){
829 origin = fkBeautyCharm;
830 return origin;
831 }
832 }
833
834 partMother = grandMa;
835 } // end of iteration
836 } // end of if
837
838 //gamma --------------------------
839 else if ( abs(maPdgcode) == 22 ) {
840 origin = fkGamma;
841
842 // iterate until you find B hadron as a mother or become top ancester
843 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
844
845 Int_t jLabel = partMother->GetFirstMother();
846 if (jLabel == -1){
847 origin = fkGamma;
848 return origin;
849 }
850 if (jLabel < 0){ // safety protection
851 AliDebug(1, "Stack label is negative, return\n");
852 return -1;
853 }
854
855 // if there is an ancester
856 TParticle* grandMa = fStack->Particle(jLabel);
857 Int_t grandMaPDG = grandMa->GetPdgCode();
858
859 for (Int_t j=0; j<fNparents; j++){
860 if (abs(grandMaPDG)==fParentSelect[1][j]){
861 origin = fkBeautyGamma;
862 return origin;
863 }
864 }
865
866 partMother = grandMa;
867 } // end of iteration
868
869 return origin;
870 } // end of if
871
872 //pi0 --------------------------
873 else if ( abs(maPdgcode) == 111 ) {
874 origin = fkPi0;
875
876 // iterate until you find B hadron as a mother or become top ancester
877 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
878
879 Int_t jLabel = partMother->GetFirstMother();
880 if (jLabel == -1){
881 origin = fkPi0;
882 return origin;
883 }
884 if (jLabel < 0){ // safety protection
885 AliDebug(1, "Stack label is negative, return\n");
886 return -1;
887 }
888
889 // if there is an ancester
890 TParticle* grandMa = fStack->Particle(jLabel);
891 Int_t grandMaPDG = grandMa->GetPdgCode();
892
893 for (Int_t j=0; j<fNparents; j++){
894 if (abs(grandMaPDG)==fParentSelect[1][j]){
895 origin = fkBeautyPi0;
896 return origin;
897 }
898 }
899
900 partMother = grandMa;
901 } // end of iteration
902
903 return origin;
904 } // end of if
905
906
907 //else --------------------------
908 else {
909 origin = fkElse;
910
911 // iterate until you find B hadron as a mother or become top ancester
912 for (Int_t i=0; i<100; i++){ // check back to the 100 generation older
913
914 Int_t jLabel = partMother->GetFirstMother();
915 if (jLabel == -1){
916 origin = fkElse;
917 return origin;
918 }
919 if (jLabel < 0){ // safety protection
920 AliDebug(1, "Stack label is negative, return\n");
921 return -1;
922 }
923
924 // if there is an ancester
925 TParticle* grandMa = fStack->Particle(jLabel);
926 Int_t grandMaPDG = grandMa->GetPdgCode();
927
928 for (Int_t j=0; j<fNparents; j++){
929 if (abs(grandMaPDG)==fParentSelect[1][j]){
930 origin = fkBeautyElse;
931 return origin;
932 }
933 }
934
935 partMother = grandMa;
936 } // end of iteration
937
938 }
939
940 return origin;
941
942}
943
944//_______________________________________________________________________________________________
945Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track)
946{
947
948 //if (track->Pt() < 1.0) return kFALSE;
949 //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
950 //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
951 //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
952 if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
953 //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
954
955
956/*
957 Float_t dcaR=-1;
958 Float_t dcaZ=-1;
959 track->GetImpactParameters(dcaR,dcaZ);
960 if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
961 if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
962*/
963 return kTRUE;
964}