]> 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 classes (Markus)
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEsecVtx.cxx
CommitLineData
9bcfd1ab 1 /**************************************************************************
259c3296 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>
9bcfd1ab 28#include <AliAODEvent.h>
29#include <AliVTrack.h>
259c3296 30#include <AliESDtrack.h>
9bcfd1ab 31#include <AliAODTrack.h>
259c3296 32#include "AliHFEsecVtx.h"
33#include <AliKFParticle.h>
34#include <AliKFVertex.h>
35#include <AliLog.h>
259c3296 36#include <AliStack.h>
9bcfd1ab 37#include <AliAODMCParticle.h>
38#include "AliHFEpairs.h"
39#include "AliHFEsecVtxs.h"
259c3296 40
41ClassImp(AliHFEsecVtx)
42
43//_______________________________________________________________________________________________
44AliHFEsecVtx::AliHFEsecVtx():
9bcfd1ab 45 fESD1(0x0)
46 ,fAOD1(0x0)
47 ,fStack(0x0)
48 ,fUseMCPID(kFALSE)
49 ,fkSourceLabel()
50 ,fNparents(0)
51 ,fParentSelect()
52 ,fNoOfHFEpairs(0)
53 ,fNoOfHFEsecvtxs(0)
54 ,fHFEpairs(0x0)
55 ,fHFEsecvtxs(0x0)
56 ,fMCArray(0x0)
57 ,fPVx(0)
58 ,fPVy(0)
59 ,fCosPhi(-1)
60 ,fSignedLxy(-1)
61 ,fKFchi2(-1)
62 ,fInvmass(-1)
63 ,fInvmassSigma(-1)
64 ,fKFip(0)
65 ,fPairQA(0x0)
66 ,fSecvtxQA(0x0)
67 ,fSecVtxList(0x0)
259c3296 68{
9bcfd1ab 69 //
70 // Default constructor
71 //
259c3296 72
9bcfd1ab 73 Init();
259c3296 74}
75
76//_______________________________________________________________________________________________
77AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
9bcfd1ab 78 TObject(p)
79 ,fESD1(0x0)
80 ,fAOD1(0x0)
81 ,fStack(0x0)
82 ,fUseMCPID(p.fUseMCPID)
83 ,fkSourceLabel()
84 ,fNparents(p.fNparents)
85 ,fParentSelect()
86 ,fNoOfHFEpairs(p.fNoOfHFEpairs)
87 ,fNoOfHFEsecvtxs(p.fNoOfHFEsecvtxs)
88 ,fHFEpairs(0x0)
89 ,fHFEsecvtxs(0x0)
90 ,fMCArray(0x0)
91 ,fPVx(p.fPVx)
92 ,fPVy(p.fPVy)
93 ,fCosPhi(p.fCosPhi)
94 ,fSignedLxy(p.fSignedLxy)
95 ,fKFchi2(p.fKFchi2)
96 ,fInvmass(p.fInvmass)
97 ,fInvmassSigma(p.fInvmassSigma)
98 ,fKFip(p.fKFip)
99 ,fPairQA(0x0)
100 ,fSecvtxQA(0x0)
101 ,fSecVtxList(0x0)
259c3296 102{
9bcfd1ab 103 //
104 // Copy constructor
105 //
259c3296 106}
107
108//_______________________________________________________________________________________________
109AliHFEsecVtx&
110AliHFEsecVtx::operator=(const AliHFEsecVtx &)
111{
9bcfd1ab 112 //
259c3296 113 // Assignment operator
9bcfd1ab 114 //
259c3296 115
116 AliInfo("Not yet implemented.");
117 return *this;
118}
119
120//_______________________________________________________________________________________________
121AliHFEsecVtx::~AliHFEsecVtx()
122{
9bcfd1ab 123 //
124 // Destructor
125 //
259c3296 126
9bcfd1ab 127 //cout << "Analysis Done." << endl;
259c3296 128}
129
130//__________________________________________
131void AliHFEsecVtx::Init()
132{
9bcfd1ab 133 //
134 // set pdg code and index
135 //
136
137 fNparents = 7;
138
139 fParentSelect[0][0] = 411;
140 fParentSelect[0][1] = 421;
141 fParentSelect[0][2] = 431;
142 fParentSelect[0][3] = 4122;
143 fParentSelect[0][4] = 4132;
144 fParentSelect[0][5] = 4232;
145 fParentSelect[0][6] = 4332;
146
147 fParentSelect[1][0] = 511;
148 fParentSelect[1][1] = 521;
149 fParentSelect[1][2] = 531;
150 fParentSelect[1][3] = 5122;
151 fParentSelect[1][4] = 5132;
152 fParentSelect[1][5] = 5232;
153 fParentSelect[1][6] = 5332;
259c3296 154}
155
9bcfd1ab 156//_______________________________________________________________________________________________
157void AliHFEsecVtx::CreateHistograms(TList *qaList)
158{
159 //
160 // create histograms
161 //
162
163 fSecVtxList = new TList;
164 fSecVtxList->SetName("SecVtx");
165
166 MakeContainer();
167 /*
168 fkSourceLabel[kAll]="all";
169 fkSourceLabel[kDirectCharm]="directCharm";
170 fkSourceLabel[kDirectBeauty]="directBeauty";
171 fkSourceLabel[kBeautyCharm]="beauty2charm";
172 fkSourceLabel[kGamma]="gamma";
173 fkSourceLabel[kPi0]="pi0";
174 fkSourceLabel[kElse]="others";
175 fkSourceLabel[kBeautyGamma]="beauty22gamma";
176 fkSourceLabel[kBeautyPi0]="beauty22pi0";
177 fkSourceLabel[kBeautyElse]="beauty22others";
178
179 TString hname;
180 TString hnopt="secvtx_";
181 for (Int_t isource = 0; isource < 10; isource++ ){
182 }
183 */
184
185 qaList->AddLast(fSecVtxList);
259c3296 186}
187
9bcfd1ab 188//_______________________________________________________________________________________________
189void AliHFEsecVtx::GetPrimaryCondition()
259c3296 190{
9bcfd1ab 191 //
192 // get primary characteristics and set
193 //
194
195 if (fESD1) {
196 AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
197 if( primVtxCopy.GetNDF() <1 ) return;
198 fPVx = primVtxCopy.GetX();
199 fPVx = primVtxCopy.GetY();
200 }
201 else if(fAOD1) {
202 AliKFVertex primVtxCopy(*(fAOD1->GetPrimaryVertex()));
203 if( primVtxCopy.GetNDF() <1 ) return;
204 fPVx = primVtxCopy.GetX();
205 fPVx = primVtxCopy.GetY();
206 }
259c3296 207}
208
209//_______________________________________________________________________________________________
9bcfd1ab 210void AliHFEsecVtx::PairAnalysis(AliVTrack* track1, AliVTrack* track2, Int_t index2)
211{
212 //
213 // calculate e-h pair characteristics and tag pair
214 //
259c3296 215
9bcfd1ab 216 // get KF particle input pid
217 Int_t pdg1 = GetPDG(track1);
218 Int_t pdg2 = GetPDG(track2);
219
220 if(pdg1==-1 || pdg2==-1) {
221 //printf("out if considered pid range \n");
222 return;
223 }
224
225 // create KF particle of pair
226 if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
227 else AliKFParticle::SetField(fESD1->GetMagneticField());
228 AliKFParticle kfTrack1(*track1, pdg1);
229 AliKFParticle kfTrack2(*track2, pdg2);
230
231 AliKFParticle kfSecondary(kfTrack1,kfTrack2);
232
233 //secondary vertex point from kf particle
234 Double_t kfx = kfSecondary.GetX();
235 Double_t kfy = kfSecondary.GetY();
236 //Double_t kfz = kfSecondary.GetZ();
237
238 //momentum at the decay point from kf particle
239 Double_t kfpx = kfSecondary.GetPx();
240 Double_t kfpy = kfSecondary.GetPy();
241 //Double_t kfpz = kfSecondary.GetPz();
242
243 Double_t dx = kfx-fPVx;
244 Double_t dy = kfy-fPVy;
245
246 // discriminating variables
247 // invariant mass of the KF particle
248 Double_t invmass = -1;
249 Double_t invmassSigma = -1;
250 kfSecondary.GetMass(invmass,invmassSigma);
251
252 // chi2 of the KF particle
253 Double_t kfchi2 = -1;
254 if(kfSecondary.GetNDF()>0) kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
255
256 // opening angle between two particles in XY plane
257 Double_t phi = kfTrack1.GetAngleXY(kfTrack2);
258 Double_t cosphi = TMath::Cos(phi);
259
260 // projection of kf vertex vector to the kf momentum direction
261 Double_t signedLxy=-999.;
262 Double_t psqr = kfpx*kfpx+kfpy*kfpy;
263 /*
264 //[the other way to think about]
265 if(psqr>0) {
266 if((dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr)>0) signedLxy = TMath::Sqrt(dx*dx+dy*dy);
267 if((dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr)<0) signedLxy = -1*TMath::Sqrt(dx*dx+dy*dy);
268 }*/
269 if(psqr>0) signedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
270
271 // DCA from primary to e-h KF particle (impact parameter of KF particle)
272 Double_t vtx[2]={fPVx, fPVy};
273 Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
274
275 /* //later consider the below
276 Float_t dcaR1=-999., dcaR2=-999.;
277 Float_t dcaZ1=-999., dcaZ2=-999.;
278
279 if(IsAODanalysis()){
280 Double_t trkIpPar1[2];
281 Double_t trkIpCov1[3];
282 Double_t trkIpPar2[2];
283 Double_t trkIpCov2[3];
284
285 AliESDtrack esdTrk1(track1);
286 AliESDtrack esdTrk2(track2);
287 esdTrk1.PropagateToDCA(fAOD1->GetPrimaryVertex(),0,10000,trkIpPar1,trkIpCov1);
288 esdTrk2.PropagateToDCA(fAOD1->GetPrimaryVertex(),0,10000,trkIpPar2,trkIpCov2);
289
290 dcaR1=trkIpPar1[0];
291 dcaZ1=trkIpPar1[1];
292 dcaR2=trkIpPar2[0];
293 dcaZ2=trkIpPar2[1];
294 }
295 else {
296 ((AliESDtrack*)track1)->GetImpactParameters(dcaR1,dcaZ1);
297 ((AliESDtrack*)track2)->GetImpactParameters(dcaR2,dcaZ2);
298 }*/
299
300 Int_t paircode = -1;
301 if (HasMCData()) paircode = GetPairCode(track1,track2);
302
303 AliHFEpairs hfepair;
304 hfepair.SetTrkLabel(index2);
305 hfepair.SetInvmass(invmass);
306 hfepair.SetKFChi2(kfchi2);
307 hfepair.SetOpenangle(phi);
308 hfepair.SetCosOpenangle(cosphi);
309 hfepair.SetSignedLxy(signedLxy);
310 hfepair.SetKFIP(kfip);
311 hfepair.SetPairCode(paircode);
312 AddHFEpairToArray(&hfepair);
313 fNoOfHFEpairs++;
314
315 // fill into container for later QA
316 Double_t dataE[6];
317 dataE[0]=invmass;
318 dataE[1]=kfchi2;
319 dataE[2]=phi;
320 dataE[3]=signedLxy;
321 dataE[4]=kfip;
322 dataE[5]=paircode;
323 fPairQA->Fill(dataE);
259c3296 324}
325
326//_______________________________________________________________________________________________
9bcfd1ab 327void AliHFEsecVtx::RunSECVTX(AliVTrack *track)
259c3296 328{
9bcfd1ab 329 //
330 // run secondary vertexing algorithm and do tagging
331 //
259c3296 332
9bcfd1ab 333 //printf("number of considered pairs= %d\n",HFEpairs()->GetEntriesFast());
334 FindSECVTXCandid(track);
335}
259c3296 336
9bcfd1ab 337//_______________________________________________________________________________________________
338void AliHFEsecVtx::FindSECVTXCandid(AliVTrack *track)
339{
340 //
341 // find secondary vertex candidate and store them into container
342 //
343
344 AliVTrack *htrack[20];
345 Int_t htracklabel[20];
346 Double_t vtxchi2cut=3.; // testing cut
347 Double_t dataE[6]={-999.,-999.,-999.,-999.,-1.,0};
348 if (HFEpairs()->GetEntriesFast()>20){
349 AliDebug(3, "number of paired hadron is over maximum(20)");
350 return;
351 }
259c3296 352
9bcfd1ab 353 // get paired track objects
354 AliHFEpairs *pair=0x0;
355 if (IsAODanalysis()){
356 for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
357 pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
358 htracklabel[ip] = pair->GetTrkLabel();
359 htrack[ip] = fAOD1->GetTrack(pair->GetTrkLabel());
360 }
361 }
362 else{
363 for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
364 pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
365 htracklabel[ip] = pair->GetTrkLabel();
366 htrack[ip] = fESD1->GetTrack(pair->GetTrkLabel());
367 }
368 }
369 // in case there is only one paired track with the electron, put pair characteristics into secvtx container
370 // for the moment, I only apply pair vertex chi2 cut
371 if (HFEpairs()->GetEntriesFast() == 1){
372 if (pair->GetKFChi2()<vtxchi2cut) { // you can also put single track cut
373 AliHFEsecVtxs hfesecvtx;
374 hfesecvtx.SetTrkLabel1(pair->GetTrkLabel());
375 hfesecvtx.SetTrkLabel2(-999);
376 hfesecvtx.SetInvmass(pair->GetInvmass());
377 hfesecvtx.SetKFChi2(pair->GetKFChi2());
378 hfesecvtx.SetSignedLxy(pair->GetSignedLxy());
379 hfesecvtx.SetKFIP(pair->GetKFIP());
380 AddHFEsecvtxToArray(&hfesecvtx);
381 fNoOfHFEsecvtxs++;
382
383 dataE[0]=pair->GetInvmass();
384 dataE[1]=pair->GetKFChi2();
385 dataE[2]=pair->GetSignedLxy();
386 dataE[3]=pair->GetKFIP();
387 if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
388 dataE[5]=2;
389 fSecvtxQA->Fill(dataE);
390 }
391 return;
392 }
393
394 // in case there are multiple paired track with the electron, calculate secvtx characteristics
395 // put the secvtx characteristics into container if it passes cuts
396 for (int i=0; i<HFEpairs()->GetEntriesFast()-1; i++){
397 for (int j=i+1; j<HFEpairs()->GetEntriesFast(); j++){
398 CalcSECVTXProperty(track, htrack[i], htrack[j]);
399 if (fKFchi2<vtxchi2cut) {
400 AliHFEsecVtxs hfesecvtx;
401 hfesecvtx.SetTrkLabel1(htracklabel[i]);
402 hfesecvtx.SetTrkLabel2(htracklabel[j]);
403 hfesecvtx.SetKFChi2(fKFchi2);
404 hfesecvtx.SetInvmass(fInvmass);
405 hfesecvtx.SetSignedLxy(fSignedLxy);
406 hfesecvtx.SetKFIP(fKFip);
407 AddHFEsecvtxToArray(&hfesecvtx);
408 fNoOfHFEsecvtxs++;
409
410 dataE[0]=pair->GetInvmass();
411 dataE[1]=pair->GetKFChi2();
412 dataE[2]=pair->GetSignedLxy();
413 dataE[3]=pair->GetKFIP();
414 if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
415 dataE[5]=3;
416 fSecvtxQA->Fill(dataE);
259c3296 417 }
9bcfd1ab 418 }
419 }
259c3296 420}
421
422//_______________________________________________________________________________________________
9bcfd1ab 423void AliHFEsecVtx::CalcSECVTXProperty(AliVTrack* track1, AliVTrack* track2, AliVTrack* track3)
259c3296 424{
9bcfd1ab 425 //
426 // calculate secondary vertex properties
427 //
428
429 // get KF particle input pid
430 Int_t pdg1 = GetPDG(track1);
431 Int_t pdg2 = GetPDG(track2);
432 Int_t pdg3 = GetPDG(track3);
433
434 if(pdg1==-1 || pdg2==-1 || pdg3==-1) {
435 //printf("out if considered pid range \n");
436 return;
437 }
438
439 // create KF particle of pair
440 if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
441 else AliKFParticle::SetField(fESD1->GetMagneticField());
442 AliKFParticle kfTrack1(*track1, pdg1);
443 AliKFParticle kfTrack2(*track2, pdg2);
444 AliKFParticle kfTrack3(*track3, pdg3);
445
446 AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
447
448 //secondary vertex point from kf particle
449 Double_t kfx = kfSecondary.GetX();
450 Double_t kfy = kfSecondary.GetY();
451 //Double_t kfz = kfSecondary.GetZ();
259c3296 452
9bcfd1ab 453 //momentum at the decay point from kf particle
454 Double_t kfpx = kfSecondary.GetPx();
455 Double_t kfpy = kfSecondary.GetPy();
456 //Double_t kfpz = kfSecondary.GetPz();
259c3296 457
9bcfd1ab 458 Double_t dx = kfx-fPVx;
459 Double_t dy = kfy-fPVy;
dbe3abbe 460
9bcfd1ab 461 // discriminating variables ----------------------------------------------------------
dbe3abbe 462
9bcfd1ab 463 if(kfSecondary.GetNDF()>0) fKFchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
259c3296 464
9bcfd1ab 465 // invariant mass of the KF particle
466 kfSecondary.GetMass(fInvmass,fInvmassSigma);
259c3296 467
9bcfd1ab 468 // DCA from primary to e-h KF particle (impact parameter of KF particle)
469 Double_t vtx[2]={fPVx, fPVy};
470 fKFip = kfSecondary.GetDistanceFromVertexXY(vtx);
259c3296 471
9bcfd1ab 472 // projection of kf vertex vector to the kf momentum direction
473 Double_t psqr = kfpx*kfpx+kfpy*kfpy;
474 if(psqr>0) fSignedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
259c3296 475}
476
477//_______________________________________________________________________________________________
9bcfd1ab 478Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track)
259c3296 479{
9bcfd1ab 480 //
481 // return mc pid
482 //
259c3296 483
9bcfd1ab 484 Int_t label = TMath::Abs(track->GetLabel());
485 TParticle* mcpart = fStack->Particle(label);
486 if ( !mcpart ) return 0;
487 Int_t pdgCode = mcpart->GetPdgCode();
259c3296 488
9bcfd1ab 489 return pdgCode;
259c3296 490}
491
492//_______________________________________________________________________________________________
9bcfd1ab 493Int_t AliHFEsecVtx::GetPairOriginESD(AliESDtrack* trk1, AliESDtrack* trk2)
259c3296 494{
9bcfd1ab 495 //
496 // return pdg code of the origin(source) of the pair
497 //
498 //
499 // ---*---*---*-----ancester A----- track1
500 // |____*______
501 // |______ track2
502 // => if they originated from same ancester,
503 // then return "the absolute value of pdg code of ancester A"
504 //
505 // ---*---*---B hadron-----ancester A----- track1
506 // |____*______
507 // |______ track2
508 // => if they originated from same ancester, and this ancester originally comes from B hadrons
509 // then return -1*"the absolute value of pdg code of ancester A"
510 //
511 // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
512 //
513
514 if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
515 TParticle* part1 = fStack->Particle(trk1->GetLabel());
516 TParticle* part2 = fStack->Particle(trk2->GetLabel());
517 TParticle* part2cp = part2;
518 if (!(part1) || !(part2)) return 0;
519
520 Int_t srcpdg = 0;
521
522 //if the two tracks' mother's label is same, get the mother info
523 //in case of charm, check if it originated from beauty
524 for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
525 Int_t label1 = part1->GetFirstMother();
526 if (label1 < 0) return 0;
527
528 for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
529 Int_t label2 = part2->GetFirstMother();
530 if (label2 < 0) break;
531
532 if (label1 == label2){ //check if two tracks are originated from same mother
533 TParticle* commonmom = fStack->Particle(label2);
534 srcpdg = abs(commonmom->GetPdgCode());
535
536 //check ancester to see if it is originally from beauty
537 for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
538 Int_t ancesterlabel = commonmom->GetFirstMother();
539 if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
540
541 TParticle* commonancester = fStack->Particle(ancesterlabel);
542 Int_t ancesterpdg = abs(commonancester->GetPdgCode());
543
544 for (Int_t l=0; l<fNparents; l++){
545 if (abs(ancesterpdg)==fParentSelect[1][l]){
546 srcpdg = -1*srcpdg; //multiply -1 for hadron from bottom
547 return srcpdg;
548 }
549 }
550 commonmom = commonancester;
551 }
552 }
553 part2 = fStack->Particle(label2); //if their mother is different, go to earlier generation of 2nd particle
554 if (!(part2)) break;
555 }
556 part1 = fStack->Particle(label1); //if their mother is different, go to earlier generation of 1st particle
557 part2 = part2cp;
558 if (!(part1)) return 0;
559 }
560
561 return srcpdg;
259c3296 562}
563
564//_______________________________________________________________________________________________
9bcfd1ab 565Int_t AliHFEsecVtx::GetPairOriginAOD(AliAODTrack* trk1, AliAODTrack* trk2)
259c3296 566{
567
9bcfd1ab 568 //
569 // return pdg code of the origin(source) of the pair
570 //
571 //
572 // ---*---*---*-----ancester A----- track1
573 // |____*______
574 // |______ track2
575 // => if they originated from same ancester,
576 // then return "the absolute value of pdg code of ancester A"
577 //
578 // ---*---*---B hadron-----ancester A----- track1
579 // |____*______
580 // |______ track2
581 // => if they originated from same ancester, and this ancester originally comes from B hadrons
582 // then return -1*"the absolute value of pdg code of ancester A"
583 //
584 // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
585 //
586
587 if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
588 AliAODMCParticle *part1 = (AliAODMCParticle*)fMCArray->At(trk1->GetLabel());
589 AliAODMCParticle *part2 = (AliAODMCParticle*)fMCArray->At(trk2->GetLabel());
590 AliAODMCParticle *part2cp = part2;
591 if (!(part1) || !(part2)) return 0;
592
593 Int_t srcpdg = 0;
594
595 //if the two tracks' mother's label is same, get the mother info
596 //in case of charm, check if it originated from beauty
597 for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
598 Int_t label1 = part1->GetMother();
599 if (label1 < 0) return 0;
600
601 for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
602 Int_t label2 = part2->GetMother();
603 if (label2 < 0) return 0;
604
605 if (label1 == label2){ //check if two tracks are originated from same mother
606 AliAODMCParticle *commonmom = (AliAODMCParticle*)fMCArray->At(label1);
607 srcpdg = abs(commonmom->GetPdgCode());
608
609 //check ancester to see if it is originally from beauty
610 for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
611 Int_t ancesterlabel = commonmom->GetMother();
612 if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
613
614 AliAODMCParticle *commonancester = (AliAODMCParticle*)fMCArray->At(ancesterlabel);
615 Int_t ancesterpdg = abs(commonancester->GetPdgCode());
616
617 for (Int_t l=0; l<fNparents; l++){
618 if (abs(ancesterpdg)==fParentSelect[1][l]){
619 srcpdg = -1*srcpdg; //multiply -1 for charm from bottom
620 return srcpdg;
621 }
622 }
623 commonmom = commonancester;
624 }
625 }
626 part2 = (AliAODMCParticle*)fMCArray->At(label2); //if their mother is different, go to earlier generation of 2nd particle
627 if (!(part2)) break;
628 }
629 part1 = (AliAODMCParticle*)fMCArray->At(label1); //if their mother is different, go to earlier generation of 2nd particle
630 part2 = part2cp;
631 if (!(part1)) return 0;
632 }
633
634 return srcpdg;
259c3296 635}
636
637//_______________________________________________________________________________________________
9bcfd1ab 638Int_t AliHFEsecVtx::GetPairCode(const AliVTrack* const trk1, const AliVTrack* const trk2)
259c3296 639{
9bcfd1ab 640 //
641 // return pair code which is predefinded as:
642 // kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0,
643 // kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
644 //
645
646 Int_t srcpdg = -1;
647 Int_t srccode = kElse;
648
649 if(IsAODanalysis()) srcpdg = GetPairOriginAOD((AliAODTrack*)trk1,(AliAODTrack*)trk2);
650 else srcpdg = GetPairOriginESD((AliESDtrack*)trk1,(AliESDtrack*)trk2);
651
652 if (srcpdg < 0) srccode = kBeautyElse;
653 for (Int_t i=0; i<fNparents; i++){
654 if (abs(srcpdg)==fParentSelect[0][i]) {
655 if (srcpdg>0) srccode = kDirectCharm;
656 if (srcpdg<0) srccode = kBeautyCharm;
657 }
658 if (abs(srcpdg)==fParentSelect[1][i]) {
659 if (srcpdg>0) srccode = kDirectBeauty;
660 if (srcpdg<0) return kElse;
661 }
662 }
663 if (srcpdg == 22) srccode = kGamma;
664 if (srcpdg == -22) srccode = kBeautyGamma;
665 if (srcpdg == 111) srccode = kPi0;
666 if (srcpdg == -111) srccode = kBeautyPi0;
259c3296 667
9bcfd1ab 668 return srccode;
259c3296 669}
670
671//_______________________________________________________________________________________________
9bcfd1ab 672Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack)
259c3296 673{
9bcfd1ab 674 //
675 // return decay electron's origin
676 //
259c3296 677
9bcfd1ab 678 if (iTrack < 0) {
679 AliDebug(1, "Stack label is negative, return\n");
680 return -1;
681 }
259c3296 682
9bcfd1ab 683 TParticle* mcpart = fStack->Particle(iTrack);
41d0c3b3 684
9bcfd1ab 685// if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !
259c3296 686
9bcfd1ab 687 Int_t iLabel = mcpart->GetFirstMother();
688 if (iLabel<0){
689 AliDebug(1, "Stack label is negative, return\n");
690 return -1;
691 }
259c3296 692
9bcfd1ab 693 Int_t origin = -1;
694 Bool_t isFinalOpenCharm = kFALSE;
259c3296 695
9bcfd1ab 696 TParticle *partMother = fStack->Particle(iLabel);
697 Int_t maPdgcode = partMother->GetPdgCode();
259c3296 698
9bcfd1ab 699 // if the mother is charmed hadron
700 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
259c3296 701
9bcfd1ab 702 for (Int_t i=0; i<fNparents; i++){
703 if (abs(maPdgcode)==fParentSelect[0][i]){
704 isFinalOpenCharm = kTRUE;
705 }
706 }
707 if (!isFinalOpenCharm) return -1;
259c3296 708
9bcfd1ab 709 // iterate until you find B hadron as a mother or become top ancester
710 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
259c3296 711
9bcfd1ab 712 Int_t jLabel = partMother->GetFirstMother();
713 if (jLabel == -1){
714 origin = kDirectCharm;
715 return origin;
716 }
717 if (jLabel < 0){ // safety protection
718 AliDebug(1, "Stack label is negative, return\n");
719 return -1;
720 }
721
722 // if there is an ancester
723 TParticle* grandMa = fStack->Particle(jLabel);
724 Int_t grandMaPDG = grandMa->GetPdgCode();
725
726 for (Int_t j=0; j<fNparents; j++){
727 if (abs(grandMaPDG)==fParentSelect[1][j]){
728 origin = kBeautyCharm;
729 return origin;
dbe3abbe 730 }
9bcfd1ab 731 }
732
733 partMother = grandMa;
734 } // end of iteration
735 } // end of if
736 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
737 for (Int_t i=0; i<fNparents; i++){
738 if (abs(maPdgcode)==fParentSelect[1][i]){
739 origin = kDirectBeauty;
740 return origin;
741 }
742 }
743 } // end of if
dbe3abbe 744
9bcfd1ab 745 //============ gamma ================
746 else if ( abs(maPdgcode) == 22 ) {
747 origin = kGamma;
259c3296 748
9bcfd1ab 749 // iterate until you find B hadron as a mother or become top ancester
750 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
259c3296 751
9bcfd1ab 752 Int_t jLabel = partMother->GetFirstMother();
753 if (jLabel == -1){
754 origin = kGamma;
755 return origin;
756 }
757 if (jLabel < 0){ // safety protection
758 AliDebug(1, "Stack label is negative, return\n");
759 return -1;
760 }
761
762 // if there is an ancester
763 TParticle* grandMa = fStack->Particle(jLabel);
764 Int_t grandMaPDG = grandMa->GetPdgCode();
765
766 for (Int_t j=0; j<fNparents; j++){
767 if (abs(grandMaPDG)==fParentSelect[1][j]){
768 origin = kBeautyGamma;
769 return origin;
770 }
771 }
259c3296 772
9bcfd1ab 773 partMother = grandMa;
774 } // end of iteration
259c3296 775
9bcfd1ab 776 return origin;
777 } // end of if
259c3296 778
9bcfd1ab 779 //============ pi0 ================
780 else if ( abs(maPdgcode) == 111 ) {
781 origin = kPi0;
259c3296 782
9bcfd1ab 783 // iterate until you find B hadron as a mother or become top ancester
784 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
259c3296 785
9bcfd1ab 786 Int_t jLabel = partMother->GetFirstMother();
787 if (jLabel == -1){
788 origin = kPi0;
789 return origin;
790 }
791 if (jLabel < 0){ // safety protection
792 AliDebug(1, "Stack label is negative, return\n");
793 return -1;
794 }
795
796 // if there is an ancester
797 TParticle* grandMa = fStack->Particle(jLabel);
798 Int_t grandMaPDG = grandMa->GetPdgCode();
799
800 for (Int_t j=0; j<fNparents; j++){
801 if (abs(grandMaPDG)==fParentSelect[1][j]){
802 origin = kBeautyPi0;
803 return origin;
dbe3abbe 804 }
9bcfd1ab 805 }
dbe3abbe 806
9bcfd1ab 807 partMother = grandMa;
808 } // end of iteration
259c3296 809
9bcfd1ab 810 return origin;
811 } // end of if
259c3296 812
9bcfd1ab 813 else {
814 origin = kElse;
815 // iterate until you find B hadron as a mother or become top ancester
816 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
259c3296 817
9bcfd1ab 818 Int_t jLabel = partMother->GetFirstMother();
819 if (jLabel == -1){
820 origin = kElse;
821 return origin;
822 }
823 if (jLabel < 0){ // safety protection
824 AliDebug(1, "Stack label is negative, return\n");
825 return -1;
826 }
827
828 // if there is an ancester
829 TParticle* grandMa = fStack->Particle(jLabel);
830 Int_t grandMaPDG = grandMa->GetPdgCode();
831
832 for (Int_t j=0; j<fNparents; j++){
833 if (abs(grandMaPDG)==fParentSelect[1][j]){
834 origin = kBeautyElse;
835 return origin;
259c3296 836 }
9bcfd1ab 837 }
259c3296 838
9bcfd1ab 839 partMother = grandMa;
840 } // end of iteration
841 }
dbe3abbe 842
9bcfd1ab 843 return origin;
259c3296 844}
845
846//_______________________________________________________________________________________________
9bcfd1ab 847Int_t AliHFEsecVtx::GetPDG(AliVTrack *track)
259c3296 848{
9bcfd1ab 849 //
850 // get KF particle input pdg for mass hypothesis
851 //
852
853 Int_t pdgCode=-1;
854
855 if (fUseMCPID && HasMCData()){
856 pdgCode = GetMCPDG(track);
857 }
858 else if(fESD1){
859 Int_t pid=0;
860 Double_t prob;
861 GetESDPID((AliESDtrack*)track, pid, prob);
862 switch(pid){
863 case 0: pdgCode = 11; break;
864 case 1: pdgCode = 13; break;
865 case 2: pdgCode = 211; break;
866 case 3: pdgCode = 321; break;
867 case 4: pdgCode = 2212; break;
868 default: pdgCode = -1;
869 }
870 }
871 else if(fAOD1){
872 Int_t pid = ((AliAODTrack*)track)->GetMostProbablePID();
873 switch(pid){
874 case 0: pdgCode = 11; break;
875 case 1: pdgCode = 13; break;
876 case 2: pdgCode = 211; break;
877 case 3: pdgCode = 321; break;
878 case 4: pdgCode = 2212; break;
879 default: pdgCode = -1;
880 }
881 }
dbe3abbe 882
9bcfd1ab 883 return pdgCode;
259c3296 884}
885
886//_______________________________________________________________________________________________
9bcfd1ab 887Int_t AliHFEsecVtx::GetMCPDG(AliVTrack *track)
259c3296 888{
9bcfd1ab 889 //
890 // return mc pdg code
891 //
892
893 Int_t label = TMath::Abs(track->GetLabel());
894 Int_t pdgCode;
895
896 if (IsAODanalysis()) {
897 AliAODMCParticle *mcpart = (AliAODMCParticle*)fMCArray->At(label);
898 if ( !mcpart ) return 0;
899 pdgCode = mcpart->GetPdgCode();
900 }
901 else {
902 TParticle* mcpart = fStack->Particle(label);
903 if ( !mcpart ) return 0;
904 pdgCode = mcpart->GetPdgCode();
905 }
906
907 return pdgCode;
259c3296 908}
909
9bcfd1ab 910//______________________________________________________________________________
911void AliHFEsecVtx::GetESDPID(AliESDtrack *track, Int_t &recpid, Double_t &recprob)
259c3296 912{
9bcfd1ab 913 //
914 // calculate likehood for esd pid
915 //
259c3296 916
9bcfd1ab 917 recpid = -1;
918 recprob = -1;
259c3296 919
9bcfd1ab 920 Int_t ipid=-1;
921 Double_t max=0.;
259c3296 922
9bcfd1ab 923 Double_t probability[5];
259c3296 924
9bcfd1ab 925 // get probability of the diffenrent particle types
926 track->GetESDpid(probability);
259c3296 927
9bcfd1ab 928 // find most probable particle in ESD pid
929 // 0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
930 ipid = TMath::LocMax(5,probability);
931 max = TMath::MaxElement(5,probability);
259c3296 932
9bcfd1ab 933 recpid = ipid;
934 recprob = max;
259c3296 935}
936
9bcfd1ab 937//_____________________________________________________________________________
938void AliHFEsecVtx::AddHFEpairToArray(const AliHFEpairs* const pair)
259c3296 939{
9bcfd1ab 940 //
941 // Add a HFE pair to the array
942 //
259c3296 943
9bcfd1ab 944 Int_t n = HFEpairs()->GetEntriesFast();
945 if(n!=fNoOfHFEpairs)AliError(Form("fNoOfHFEpairs != HFEpairs()->GetEntriesFast %i != %i \n", fNoOfHFEpairs, n));
946 new((*HFEpairs())[n]) AliHFEpairs(*pair);
259c3296 947}
948
9bcfd1ab 949//_____________________________________________________________________________
950TClonesArray *AliHFEsecVtx::HFEpairs()
259c3296 951{
9bcfd1ab 952 //
953 // Returns the list of HFE pairs
954 //
955
956 if (!fHFEpairs) {
957 fHFEpairs = new TClonesArray("AliHFEpairs", 200);
958 }
959 return fHFEpairs;
259c3296 960}
961
9bcfd1ab 962//_____________________________________________________________________________
963void AliHFEsecVtx::DeleteHFEpairs()
259c3296 964{
9bcfd1ab 965 //
966 // Delete the list of HFE pairs
967 //
968
969 if (fHFEpairs) {
970 fHFEpairs->Delete();
971 //delete fHFEpairs;
972 }
259c3296 973}
974
9bcfd1ab 975//_____________________________________________________________________________
976void AliHFEsecVtx::InitHFEpairs()
dbe3abbe 977{
9bcfd1ab 978 //
979 // Initialization should be done before make all possible pairs for a given electron candidate
980 //
dbe3abbe 981
9bcfd1ab 982 fNoOfHFEpairs = 0;
dbe3abbe 983}
984
9bcfd1ab 985//_____________________________________________________________________________
986void AliHFEsecVtx::AddHFEsecvtxToArray(const AliHFEsecVtxs* const secvtx)
259c3296 987{
9bcfd1ab 988 //
989 // Add a HFE secondary vertex to the array
990 //
259c3296 991
9bcfd1ab 992 Int_t n = HFEsecvtxs()->GetEntriesFast();
993 if(n!=fNoOfHFEsecvtxs)AliError(Form("fNoOfHFEsecvtxs != HFEsecvtxs()->GetEntriesFast %i != %i \n", fNoOfHFEsecvtxs, n));
994 new((*HFEsecvtxs())[n]) AliHFEsecVtxs(*secvtx);
995}
259c3296 996
9bcfd1ab 997//_____________________________________________________________________________
998TClonesArray *AliHFEsecVtx::HFEsecvtxs()
999{
1000 //
1001 // Returns the list of HFE secvtx
1002 //
1003
1004 if (!fHFEsecvtxs) {
1005 fHFEsecvtxs = new TClonesArray("AliHFEsecVtxs", 200);
1006 }
1007 return fHFEsecvtxs;
1008}
259c3296 1009
9bcfd1ab 1010//_____________________________________________________________________________
1011void AliHFEsecVtx::DeleteHFEsecvtxs()
1012{
1013 //
1014 // Delete the list of HFE pairs
1015 //
1016
1017 if (fHFEsecvtxs) {
1018 fHFEsecvtxs->Delete();
1019 //delete fHFEsecvtx;
1020 }
1021}
259c3296 1022
9bcfd1ab 1023//_____________________________________________________________________________
1024void AliHFEsecVtx::InitHFEsecvtxs()
1025{
1026 //
1027 // Initialization should be done
1028 //
259c3296 1029
9bcfd1ab 1030 fNoOfHFEsecvtxs = 0;
259c3296 1031}
1032
1033//_______________________________________________________________________________________________
75d81601 1034Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track) const
259c3296 1035{
9bcfd1ab 1036 //if (track->Pt() < 1.0) return kFALSE;
1037 //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
1038 //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
1039 //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1040 if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
1041 //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
259c3296 1042
1043/*
9bcfd1ab 1044 Float_t dcaR=-1;
1045 Float_t dcaZ=-1;
1046 track->GetImpactParameters(dcaR,dcaZ);
1047 if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
1048 if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
259c3296 1049*/
9bcfd1ab 1050 return kTRUE;
1051}
1052//____________________________________________________________
1053void AliHFEsecVtx::MakeContainer(){
1054
1055 //
1056 // make container
1057 //
1058
1059 const Int_t nDimPair=6;
1060 Int_t nBinPair[nDimPair] = {200, 500, 314, 2000, 2000, 11};
1061 const Double_t kInvmassmin = 0., kInvmassmax = 20.;
1062 const Double_t kKFChi2min = 0, kKFChi2max= 50;
1063 const Double_t kOpenanglemin = 0, kOpenanglemax = 3.14;
1064 const Double_t kSignedLxymin = -10, kSignedLxymax= 10;
1065 const Double_t kKFIPmin = -10, kKFIPmax= 10;
1066 const Double_t kPairCodemin = -1, kPairCodemax= 10;
1067
1068 Double_t* binEdgesPair[nDimPair];
1069 for(Int_t ivar = 0; ivar < nDimPair; ivar++)
1070 binEdgesPair[ivar] = new Double_t[nBinPair[ivar] + 1];
1071
1072 for(Int_t i=0; i<=nBinPair[0]; i++) binEdgesPair[0][i]=(Double_t)kInvmassmin + (kInvmassmax - kInvmassmin)/nBinPair[0]*(Double_t)i;
1073 for(Int_t i=0; i<=nBinPair[1]; i++) binEdgesPair[1][i]=(Double_t)kKFChi2min + (kKFChi2max - kKFChi2min)/nBinPair[1]*(Double_t)i;
1074 for(Int_t i=0; i<=nBinPair[2]; i++) binEdgesPair[2][i]=(Double_t)kOpenanglemin + (kOpenanglemax - kOpenanglemin)/nBinPair[2]*(Double_t)i;
1075 for(Int_t i=0; i<=nBinPair[3]; i++) binEdgesPair[3][i]=(Double_t)kSignedLxymin + (kSignedLxymax - kSignedLxymin)/nBinPair[3]*(Double_t)i;
1076 for(Int_t i=0; i<=nBinPair[4]; i++) binEdgesPair[4][i]=(Double_t)kKFIPmin + (kKFIPmax - kKFIPmin)/nBinPair[4]*(Double_t)i;
1077 for(Int_t i=0; i<=nBinPair[5]; i++) binEdgesPair[5][i]=(Double_t)kPairCodemin + (kPairCodemax - kPairCodemin)/nBinPair[5]*(Double_t)i;
1078
1079 fPairQA = new THnSparseF("pairQA", "QA for Pair; invmass[GeV/c^2]; KF chi2; opening angle; signed Lxy; KF ip; pair code ", nDimPair, nBinPair);
1080 for(Int_t idim = 0; idim < nDimPair; idim++){
1081 fPairQA->SetBinEdges(idim, binEdgesPair[idim]);
1082 }
1083
1084 fSecVtxList->AddAt(fPairQA,0);
1085
1086 const Int_t nDimSecvtx=6;
1087 Double_t* binEdgesSecvtx[nDimSecvtx];
1088 Int_t nBinSecvtx[nDimSecvtx] = {200, 500, 2000, 2000, 11, 3};
1089 const Double_t kNtrksmin = 0, kNtrksmax= 3;
1090 for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
1091 binEdgesSecvtx[ivar] = new Double_t[nBinSecvtx[ivar] + 1];
1092
1093 for(Int_t i=0; i<=nBinSecvtx[0]; i++) binEdgesSecvtx[0][i]=binEdgesPair[0][i];
1094 for(Int_t i=0; i<=nBinSecvtx[1]; i++) binEdgesSecvtx[1][i]=binEdgesPair[1][i];
1095 for(Int_t i=0; i<=nBinSecvtx[2]; i++) binEdgesSecvtx[2][i]=binEdgesPair[3][i];
1096 for(Int_t i=0; i<=nBinSecvtx[3]; i++) binEdgesSecvtx[3][i]=binEdgesPair[4][i];
1097 for(Int_t i=0; i<=nBinSecvtx[4]; i++) binEdgesSecvtx[4][i]=binEdgesPair[5][i];
1098 for(Int_t i=0; i<=nBinSecvtx[5]; i++) binEdgesSecvtx[5][i]=(Double_t)kNtrksmin + (kNtrksmax - kNtrksmin)/nBinSecvtx[5]*(Double_t)i;
1099
1100 fSecvtxQA = new THnSparseF("secvtxQA", "QA for Secvtx; invmass[GeV/c^2]; KF chi2; signed Lxy; KF ip; pair code; n tracks ", nDimSecvtx, nBinSecvtx);
1101 for(Int_t idim = 0; idim < nDimSecvtx; idim++){
1102 fSecvtxQA->SetBinEdges(idim, binEdgesSecvtx[idim]);
1103 }
1104
1105 fSecVtxList->AddAt(fSecvtxQA,1);
259c3296 1106}