1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 // Secondary vertexing construction Class
17 // Construct secondary vertex from Beauty hadron with electron and
18 // hadrons, then apply selection criteria
21 // MinJung Kweon <minjung@physi.uni-heidelberg.de>
25 #include <TIterator.h>
26 #include <TParticle.h>
28 #include <AliESDEvent.h>
29 #include <AliAODEvent.h>
30 #include <AliVTrack.h>
31 #include <AliESDtrack.h>
32 #include <AliAODTrack.h>
33 #include "AliHFEsecVtx.h"
34 #include <AliKFParticle.h>
35 #include <AliKFVertex.h>
38 #include <AliAODMCParticle.h>
39 #include "AliHFEpairs.h"
40 #include "AliHFEsecVtxs.h"
41 #include "AliHFEtrackFilter.h"
43 ClassImp(AliHFEsecVtx)
45 //_______________________________________________________________________________________________
46 AliHFEsecVtx::AliHFEsecVtx():
75 // Default constructor
81 //_______________________________________________________________________________________________
82 AliHFEsecVtx::AliHFEsecVtx(const AliHFEsecVtx &p):
88 ,fUseMCPID(p.fUseMCPID)
90 ,fNparents(p.fNparents)
94 ,fNoOfHFEpairs(p.fNoOfHFEpairs)
95 ,fNoOfHFEsecvtxs(p.fNoOfHFEsecvtxs)
102 ,fSignedLxy(p.fSignedLxy)
104 ,fInvmass(p.fInvmass)
105 ,fInvmassSigma(p.fInvmassSigma)
114 fFilter = new AliHFEtrackFilter(*p.fFilter);
117 //_______________________________________________________________________________________________
119 AliHFEsecVtx::operator=(const AliHFEsecVtx &)
122 // Assignment operator
125 AliInfo("Not yet implemented.");
129 //_______________________________________________________________________________________________
130 AliHFEsecVtx::~AliHFEsecVtx()
136 //cout << "Analysis Done." << endl;
140 //__________________________________________
141 void AliHFEsecVtx::Init()
144 // set pdg code and index
149 fParentSelect[0][0] = 411;
150 fParentSelect[0][1] = 421;
151 fParentSelect[0][2] = 431;
152 fParentSelect[0][3] = 4122;
153 fParentSelect[0][4] = 4132;
154 fParentSelect[0][5] = 4232;
155 fParentSelect[0][6] = 4332;
157 fParentSelect[1][0] = 511;
158 fParentSelect[1][1] = 521;
159 fParentSelect[1][2] = 531;
160 fParentSelect[1][3] = 5122;
161 fParentSelect[1][4] = 5132;
162 fParentSelect[1][5] = 5232;
163 fParentSelect[1][6] = 5332;
165 // momentum ranges to apply pt dependent cuts
174 // momentum dependent DCA cut values (preliminary)
182 fFilter = new AliHFEtrackFilter("SecVtxFilter");
183 fFilter->MakeCutStepRecKineITSTPC();
184 fFilter->MakeCutStepPrimary();
187 void AliHFEsecVtx::Process(AliVTrack *signalTrack){
188 if(signalTrack->Pt() < 1.0) return;
189 AliESDtrack *track = dynamic_cast<AliESDtrack *>(signalTrack);
192 AliESDtrack *htrack = 0x0;
194 fFilter->FilterTracks(fESD1);
195 TObjArray *candidates = fFilter->GetFilteredTracks();
196 TIterator *trackIter = candidates->MakeIterator();
197 while((htrack = dynamic_cast<AliESDtrack *>(trackIter->Next()))){
198 if(track->GetID() == htrack->GetID()) continue; // since it is for tagging single electron, don't need additional condition
199 if (htrack->Pt()<1.0) continue;
200 PairAnalysis(track, htrack, htrack->GetID()); // e-h pairing
203 /*for(int ip=0; ip<fSecVtx->HFEpairs()->GetEntriesFast(); ip++){
205 AliHFEpairs *pair = (AliHFEpairs*) (fSecVtx->HFEpairs()->UncheckedAt(ip));
206 if(!(pair->GetPairCode()>1. && pair->GetPairCode()<4.)) // apply various cuts
207 fSecVtx->HFEpairs()->RemoveAt(ip);
210 HFEpairs()->Compress();
211 RunSECVTX(track); // secondary vertexing with e,h1,h2,.. tracks
212 for(int ip=0; ip<HFEsecvtxs()->GetEntriesFast(); ip++){
213 AliHFEsecVtxs *secvtx=0x0;
214 secvtx = (AliHFEsecVtxs*) (HFEsecvtxs()->UncheckedAt(ip));
215 // here you apply cuts, then if it doesn't pass the cut, remove it from the fSecVtx->HFEsecvtxs()
221 //_______________________________________________________________________________________________
222 void AliHFEsecVtx::CreateHistograms(TList *qaList)
228 fSecVtxList = new TList;
229 fSecVtxList->SetName("SecVtx");
233 fkSourceLabel[kAll]="all";
234 fkSourceLabel[kDirectCharm]="directCharm";
235 fkSourceLabel[kDirectBeauty]="directBeauty";
236 fkSourceLabel[kBeautyCharm]="beauty2charm";
237 fkSourceLabel[kGamma]="gamma";
238 fkSourceLabel[kPi0]="pi0";
239 fkSourceLabel[kElse]="others";
240 fkSourceLabel[kBeautyGamma]="beauty22gamma";
241 fkSourceLabel[kBeautyPi0]="beauty22pi0";
242 fkSourceLabel[kBeautyElse]="beauty22others";
245 TString hnopt="secvtx_";
246 for (Int_t isource = 0; isource < 10; isource++ ){
250 qaList->AddLast(fSecVtxList);
253 //_______________________________________________________________________________________________
254 void AliHFEsecVtx::GetPrimaryCondition()
257 // get primary characteristics and set
261 AliKFVertex primVtxCopy(*(fESD1->GetPrimaryVertex()));
262 if( primVtxCopy.GetNDF() <1 ) return;
263 fPVx = primVtxCopy.GetX();
264 fPVx = primVtxCopy.GetY();
267 AliKFVertex primVtxCopy(*(fAOD1->GetPrimaryVertex()));
268 if( primVtxCopy.GetNDF() <1 ) return;
269 fPVx = primVtxCopy.GetX();
270 fPVx = primVtxCopy.GetY();
274 //_______________________________________________________________________________________________
275 void AliHFEsecVtx::PairAnalysis(AliVTrack* track1, AliVTrack* track2, Int_t index2)
278 // calculate e-h pair characteristics and tag pair
281 //later consider the below
282 Float_t dca1[2]={-999.,-999.}, dca2[2]={-999.,-999.};
283 Float_t cov1[3]={-999.,-999.,-999.}, cov2[3]={-999.,-999.,-999.};
285 if (IsAODanalysis()){
286 AliESDtrack esdTrk1(track1);
287 AliESDtrack esdTrk2(track2);
288 esdTrk1.PropagateToDCA(fAOD1->GetPrimaryVertex(),0,10000,(Double_t*)dca1,(Double_t*)cov1);
289 esdTrk2.PropagateToDCA(fAOD1->GetPrimaryVertex(),0,10000,(Double_t*)dca2,(Double_t*)cov2);
292 ((AliESDtrack*)track1)->GetImpactParameters(dca1,cov1);
293 ((AliESDtrack*)track2)->GetImpactParameters(dca2,cov2);
296 // apply pt dependent dca cut on hadrons
297 for(int ibin=0; ibin<6; ibin++){
298 if((track2->Pt()>fPtRng[ibin] && track2->Pt()<fPtRng[ibin+1]) && TMath::Abs(dca2[0])<fDcaCut[ibin]) return;
301 // get KF particle input pid
302 Int_t pdg1 = GetPDG(track1);
303 Int_t pdg2 = GetPDG(track2);
305 if(pdg1==-1 || pdg2==-1) {
306 //printf("out if considered pid range \n");
310 // create KF particle of pair
311 if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
312 else AliKFParticle::SetField(fESD1->GetMagneticField());
313 AliKFParticle kfTrack1(*track1, pdg1);
314 AliKFParticle kfTrack2(*track2, pdg2);
316 AliKFParticle kfSecondary(kfTrack1,kfTrack2);
318 //secondary vertex point from kf particle
319 Double_t kfx = kfSecondary.GetX();
320 Double_t kfy = kfSecondary.GetY();
321 //Double_t kfz = kfSecondary.GetZ();
323 //momentum at the decay point from kf particle
324 Double_t kfpx = kfSecondary.GetPx();
325 Double_t kfpy = kfSecondary.GetPy();
326 //Double_t kfpz = kfSecondary.GetPz();
328 Double_t dx = kfx-fPVx;
329 Double_t dy = kfy-fPVy;
331 // discriminating variables
332 // invariant mass of the KF particle
333 Double_t invmass = -1;
334 Double_t invmassSigma = -1;
335 kfSecondary.GetMass(invmass,invmassSigma);
337 // chi2 of the KF particle
338 Double_t kfchi2 = -1;
339 if(kfSecondary.GetNDF()>0) kfchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
341 // opening angle between two particles in XY plane
342 Double_t phi = kfTrack1.GetAngleXY(kfTrack2);
343 Double_t cosphi = TMath::Cos(phi);
345 // projection of kf vertex vector to the kf momentum direction
346 Double_t signedLxy=-999.;
347 if((dx*kfpx+dy*kfpy)>0) signedLxy = TMath::Sqrt(dx*dx+dy*dy);
348 if((dx*kfpx+dy*kfpy)<0) signedLxy = -1*TMath::Sqrt(dx*dx+dy*dy);
349 //[the other way to think about]
350 //Double_t psqr = kfpx*kfpx+kfpy*kfpy;
351 //if(psqr>0) signedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
353 // DCA from primary to e-h KF particle (impact parameter of KF particle)
354 Double_t vtx[2]={fPVx, fPVy};
355 Double_t kfip = kfSecondary.GetDistanceFromVertexXY(vtx);
358 if (HasMCData()) paircode = GetPairCode(track1,track2);
361 hfepair.SetTrkLabel(index2);
362 hfepair.SetInvmass(invmass);
363 hfepair.SetKFChi2(kfchi2);
364 hfepair.SetOpenangle(phi);
365 hfepair.SetCosOpenangle(cosphi);
366 hfepair.SetSignedLxy(signedLxy);
367 hfepair.SetKFIP(kfip);
368 hfepair.SetPairCode(paircode);
369 AddHFEpairToArray(&hfepair);
372 // fill into container for later QA
381 dataE[6]=TMath::Abs(dca1[0]);
382 dataE[7]=TMath::Abs(dca2[0]);
383 //if(cov1[0]>0) dataE[6]=Double_t(dca1[0]/cov1[0]);
384 //if(cov2[0]>0) dataE[7]=Double_t(dca2[0]/cov2[0]);
385 dataE[8]=track1->Pt();
386 dataE[9]=track2->Pt();
388 fPairQA->Fill(dataE);
392 //_______________________________________________________________________________________________
393 void AliHFEsecVtx::RunSECVTX(AliVTrack *track)
396 // run secondary vertexing algorithm and do tagging
399 //printf("number of considered pairs= %d\n",HFEpairs()->GetEntriesFast());
400 FindSECVTXCandid(track);
403 //_______________________________________________________________________________________________
404 void AliHFEsecVtx::FindSECVTXCandid(AliVTrack *track)
407 // find secondary vertex candidate and store them into container
410 AliVTrack *htrack[20];
411 Int_t htracklabel[20];
412 Double_t vtxchi2cut=3.; // testing cut
413 Double_t dataE[6]={-999.,-999.,-999.,-999.,-1.,0};
414 if (HFEpairs()->GetEntriesFast()>20){
415 AliDebug(3, "number of paired hadron is over maximum(20)");
419 // get paired track objects
420 AliHFEpairs *pair=0x0;
421 if (IsAODanalysis()){
422 for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
423 pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
424 htracklabel[ip] = pair->GetTrkLabel();
425 htrack[ip] = fAOD1->GetTrack(pair->GetTrkLabel());
429 for (int ip=0; ip<HFEpairs()->GetEntriesFast(); ip++){
430 pair = (AliHFEpairs*) (HFEpairs()->UncheckedAt(ip));
431 htracklabel[ip] = pair->GetTrkLabel();
432 htrack[ip] = fESD1->GetTrack(pair->GetTrkLabel());
435 // in case there is only one paired track with the electron, put pair characteristics into secvtx container
436 // for the moment, I only apply pair vertex chi2 cut
437 if (HFEpairs()->GetEntriesFast() == 1){
438 if (pair->GetKFChi2()<vtxchi2cut) { // you can also put single track cut
439 AliHFEsecVtxs hfesecvtx;
440 hfesecvtx.SetTrkLabel1(pair->GetTrkLabel());
441 hfesecvtx.SetTrkLabel2(-999);
442 hfesecvtx.SetInvmass(pair->GetInvmass());
443 hfesecvtx.SetKFChi2(pair->GetKFChi2());
444 hfesecvtx.SetSignedLxy(pair->GetSignedLxy());
445 hfesecvtx.SetKFIP(pair->GetKFIP());
446 AddHFEsecvtxToArray(&hfesecvtx);
449 dataE[0]=pair->GetInvmass();
450 dataE[1]=pair->GetKFChi2();
451 dataE[2]=pair->GetSignedLxy();
452 dataE[3]=pair->GetKFIP();
453 if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
455 fSecvtxQA->Fill(dataE);
460 // in case there are multiple paired track with the electron, calculate secvtx characteristics
461 // put the secvtx characteristics into container if it passes cuts
462 for (int i=0; i<HFEpairs()->GetEntriesFast()-1; i++){
463 for (int j=i+1; j<HFEpairs()->GetEntriesFast(); j++){
464 CalcSECVTXProperty(track, htrack[i], htrack[j]);
465 if (fKFchi2<vtxchi2cut) {
466 AliHFEsecVtxs hfesecvtx;
467 hfesecvtx.SetTrkLabel1(htracklabel[i]);
468 hfesecvtx.SetTrkLabel2(htracklabel[j]);
469 hfesecvtx.SetKFChi2(fKFchi2);
470 hfesecvtx.SetInvmass(fInvmass);
471 hfesecvtx.SetSignedLxy(fSignedLxy);
472 hfesecvtx.SetKFIP(fKFip);
473 AddHFEsecvtxToArray(&hfesecvtx);
480 if(HasMCData()) dataE[4]=GetElectronSource(TMath::Abs(track->GetLabel()));
482 fSecvtxQA->Fill(dataE);
488 //_______________________________________________________________________________________________
489 void AliHFEsecVtx::CalcSECVTXProperty(AliVTrack* track1, AliVTrack* track2, AliVTrack* track3)
492 // calculate secondary vertex properties
495 // get KF particle input pid
496 Int_t pdg1 = GetPDG(track1);
497 Int_t pdg2 = GetPDG(track2);
498 Int_t pdg3 = GetPDG(track3);
500 if(pdg1==-1 || pdg2==-1 || pdg3==-1) {
501 //printf("out if considered pid range \n");
505 // create KF particle of pair
506 if(IsAODanalysis()) AliKFParticle::SetField(fAOD1->GetMagneticField());
507 else AliKFParticle::SetField(fESD1->GetMagneticField());
508 AliKFParticle kfTrack1(*track1, pdg1);
509 AliKFParticle kfTrack2(*track2, pdg2);
510 AliKFParticle kfTrack3(*track3, pdg3);
512 AliKFParticle kfSecondary(kfTrack1,kfTrack2,kfTrack3);
514 //secondary vertex point from kf particle
515 Double_t kfx = kfSecondary.GetX();
516 Double_t kfy = kfSecondary.GetY();
517 //Double_t kfz = kfSecondary.GetZ();
519 //momentum at the decay point from kf particle
520 Double_t kfpx = kfSecondary.GetPx();
521 Double_t kfpy = kfSecondary.GetPy();
522 //Double_t kfpz = kfSecondary.GetPz();
524 Double_t dx = kfx-fPVx;
525 Double_t dy = kfy-fPVy;
527 // discriminating variables ----------------------------------------------------------
529 if(kfSecondary.GetNDF()>0) fKFchi2 = TMath::Sqrt(TMath::Abs(kfSecondary.GetChi2()/kfSecondary.GetNDF()));
531 // invariant mass of the KF particle
532 kfSecondary.GetMass(fInvmass,fInvmassSigma);
534 // DCA from primary to e-h KF particle (impact parameter of KF particle)
535 Double_t vtx[2]={fPVx, fPVy};
536 fKFip = kfSecondary.GetDistanceFromVertexXY(vtx);
538 if((dx*kfpx+dy*kfpy)>0) fSignedLxy= TMath::Sqrt(dx*dx+dy*dy);
539 if((dx*kfpx+dy*kfpy)<0) fSignedLxy= -1*TMath::Sqrt(dx*dx+dy*dy);
540 //[the other way to think about] - projection of kf vertex vector to the kf momentum direction
541 //Double_t psqr = kfpx*kfpx+kfpy*kfpy;
542 //if(psqr>0) fSignedLxy=(dx*kfpx+dy*kfpy)/TMath::Sqrt(psqr);
545 //_______________________________________________________________________________________________
546 Int_t AliHFEsecVtx::GetMCPID(AliESDtrack *track)
552 Int_t label = TMath::Abs(track->GetLabel());
553 TParticle* mcpart = fStack->Particle(label);
554 if ( !mcpart ) return 0;
555 Int_t pdgCode = mcpart->GetPdgCode();
560 //_______________________________________________________________________________________________
561 Int_t AliHFEsecVtx::GetPairOriginESD(AliESDtrack* trk1, AliESDtrack* trk2)
564 // return pdg code of the origin(source) of the pair
567 // ---*---*---*-----ancester A----- track1
570 // => if they originated from same ancester,
571 // then return "the absolute value of pdg code of ancester A"
573 // ---*---*---B hadron-----ancester A----- track1
576 // => if they originated from same ancester, and this ancester originally comes from B hadrons
577 // then return -1*"the absolute value of pdg code of ancester A"
579 // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
582 if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
583 TParticle* part1 = fStack->Particle(trk1->GetLabel());
584 TParticle* part2 = fStack->Particle(trk2->GetLabel());
585 TParticle* part2cp = part2;
586 if (!(part1) || !(part2)) return 0;
590 //if the two tracks' mother's label is same, get the mother info
591 //in case of charm, check if it originated from beauty
592 for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
593 Int_t label1 = part1->GetFirstMother();
594 if (label1 < 0) return 0;
596 for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
597 Int_t label2 = part2->GetFirstMother();
598 if (label2 < 0) break;
600 if (label1 == label2){ //check if two tracks are originated from same mother
601 TParticle* commonmom = fStack->Particle(label2);
602 srcpdg = abs(commonmom->GetPdgCode());
604 //check ancester to see if it is originally from beauty
605 for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
606 Int_t ancesterlabel = commonmom->GetFirstMother();
607 if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
609 TParticle* commonancester = fStack->Particle(ancesterlabel);
610 Int_t ancesterpdg = abs(commonancester->GetPdgCode());
612 for (Int_t l=0; l<fNparents; l++){
613 if (abs(ancesterpdg)==fParentSelect[1][l]){
614 srcpdg = -1*srcpdg; //multiply -1 for hadron from bottom
618 commonmom = commonancester;
621 part2 = fStack->Particle(label2); //if their mother is different, go to earlier generation of 2nd particle
624 part1 = fStack->Particle(label1); //if their mother is different, go to earlier generation of 1st particle
626 if (!(part1)) return 0;
632 //_______________________________________________________________________________________________
633 Int_t AliHFEsecVtx::GetPairOriginAOD(AliAODTrack* trk1, AliAODTrack* trk2)
637 // return pdg code of the origin(source) of the pair
640 // ---*---*---*-----ancester A----- track1
643 // => if they originated from same ancester,
644 // then return "the absolute value of pdg code of ancester A"
646 // ---*---*---B hadron-----ancester A----- track1
649 // => if they originated from same ancester, and this ancester originally comes from B hadrons
650 // then return -1*"the absolute value of pdg code of ancester A"
652 // caution : it can also return parton pdg code if it originated from same string or gluon spliting.
655 if (trk1->GetLabel()<0 || trk2->GetLabel()<0) return 0;
656 AliAODMCParticle *part1 = (AliAODMCParticle*)fMCArray->At(trk1->GetLabel());
657 AliAODMCParticle *part2 = (AliAODMCParticle*)fMCArray->At(trk2->GetLabel());
658 AliAODMCParticle *part2cp = part2;
659 if (!(part1) || !(part2)) return 0;
663 //if the two tracks' mother's label is same, get the mother info
664 //in case of charm, check if it originated from beauty
665 for (Int_t i=0; i<10; i++){ //check up to 10 ancesters
666 Int_t label1 = part1->GetMother();
667 if (label1 < 0) return 0;
669 for (Int_t j=0; j<10; j++){ //check up to 10 ancesters
670 Int_t label2 = part2->GetMother();
671 if (label2 < 0) return 0;
673 if (label1 == label2){ //check if two tracks are originated from same mother
674 AliAODMCParticle *commonmom = (AliAODMCParticle*)fMCArray->At(label1);
675 srcpdg = abs(commonmom->GetPdgCode());
677 //check ancester to see if it is originally from beauty
678 for (Int_t k=0; k<10; k++){ //check up to 10 ancesters
679 Int_t ancesterlabel = commonmom->GetMother();
680 if (ancesterlabel < 0) return srcpdg; // if there is no more commonancester, return commonmom's pdg
682 AliAODMCParticle *commonancester = (AliAODMCParticle*)fMCArray->At(ancesterlabel);
683 Int_t ancesterpdg = abs(commonancester->GetPdgCode());
685 for (Int_t l=0; l<fNparents; l++){
686 if (abs(ancesterpdg)==fParentSelect[1][l]){
687 srcpdg = -1*srcpdg; //multiply -1 for charm from bottom
691 commonmom = commonancester;
694 part2 = (AliAODMCParticle*)fMCArray->At(label2); //if their mother is different, go to earlier generation of 2nd particle
697 part1 = (AliAODMCParticle*)fMCArray->At(label1); //if their mother is different, go to earlier generation of 2nd particle
699 if (!(part1)) return 0;
705 //_______________________________________________________________________________________________
706 Int_t AliHFEsecVtx::GetPairCode(const AliVTrack* const trk1, const AliVTrack* const trk2)
709 // return pair code which is predefinded as:
710 // kDirectCharm, kDirectBeauty, kBeautyCharm, kGamma, kPi0,
711 // kElse, kBeautyGamma, kBeautyPi0, kBeautyElse
715 Int_t srccode = kElse;
717 if(IsAODanalysis()) srcpdg = GetPairOriginAOD((AliAODTrack*)trk1,(AliAODTrack*)trk2);
718 else srcpdg = GetPairOriginESD((AliESDtrack*)trk1,(AliESDtrack*)trk2);
720 if (srcpdg < 0) srccode = kBeautyElse;
721 for (Int_t i=0; i<fNparents; i++){
722 if (abs(srcpdg)==fParentSelect[0][i]) {
723 if (srcpdg>0) srccode = kDirectCharm;
724 if (srcpdg<0) srccode = kBeautyCharm;
726 if (abs(srcpdg)==fParentSelect[1][i]) {
727 if (srcpdg>0) srccode = kDirectBeauty;
728 if (srcpdg<0) return kElse;
731 if (srcpdg == 22) srccode = kGamma;
732 if (srcpdg == -22) srccode = kBeautyGamma;
733 if (srcpdg == 111) srccode = kPi0;
734 if (srcpdg == -111) srccode = kBeautyPi0;
739 //_______________________________________________________________________________________________
740 Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack)
743 // return decay electron's origin
747 AliDebug(1, "Stack label is negative, return\n");
751 TParticle* mcpart = fStack->Particle(iTrack);
753 // if ( abs(mcpart->GetPdgCode()) != 11 ) return -1; // check if it is electron !
755 Int_t iLabel = mcpart->GetFirstMother();
757 AliDebug(1, "Stack label is negative, return\n");
762 Bool_t isFinalOpenCharm = kFALSE;
764 TParticle *partMother = fStack->Particle(iLabel);
765 Int_t maPdgcode = partMother->GetPdgCode();
767 // if the mother is charmed hadron
768 if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
770 for (Int_t i=0; i<fNparents; i++){
771 if (abs(maPdgcode)==fParentSelect[0][i]){
772 isFinalOpenCharm = kTRUE;
775 if (!isFinalOpenCharm) return -1;
777 // iterate until you find B hadron as a mother or become top ancester
778 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
780 Int_t jLabel = partMother->GetFirstMother();
782 origin = kDirectCharm;
785 if (jLabel < 0){ // safety protection
786 AliDebug(1, "Stack label is negative, return\n");
790 // if there is an ancester
791 TParticle* grandMa = fStack->Particle(jLabel);
792 Int_t grandMaPDG = grandMa->GetPdgCode();
794 for (Int_t j=0; j<fNparents; j++){
795 if (abs(grandMaPDG)==fParentSelect[1][j]){
796 origin = kBeautyCharm;
801 partMother = grandMa;
802 } // end of iteration
804 else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
805 for (Int_t i=0; i<fNparents; i++){
806 if (abs(maPdgcode)==fParentSelect[1][i]){
807 origin = kDirectBeauty;
813 //============ gamma ================
814 else if ( abs(maPdgcode) == 22 ) {
817 // iterate until you find B hadron as a mother or become top ancester
818 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
820 Int_t jLabel = partMother->GetFirstMother();
825 if (jLabel < 0){ // safety protection
826 AliDebug(1, "Stack label is negative, return\n");
830 // if there is an ancester
831 TParticle* grandMa = fStack->Particle(jLabel);
832 Int_t grandMaPDG = grandMa->GetPdgCode();
834 for (Int_t j=0; j<fNparents; j++){
835 if (abs(grandMaPDG)==fParentSelect[1][j]){
836 origin = kBeautyGamma;
841 partMother = grandMa;
842 } // end of iteration
847 //============ pi0 ================
848 else if ( abs(maPdgcode) == 111 ) {
851 // iterate until you find B hadron as a mother or become top ancester
852 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
854 Int_t jLabel = partMother->GetFirstMother();
859 if (jLabel < 0){ // safety protection
860 AliDebug(1, "Stack label is negative, return\n");
864 // if there is an ancester
865 TParticle* grandMa = fStack->Particle(jLabel);
866 Int_t grandMaPDG = grandMa->GetPdgCode();
868 for (Int_t j=0; j<fNparents; j++){
869 if (abs(grandMaPDG)==fParentSelect[1][j]){
875 partMother = grandMa;
876 } // end of iteration
883 // iterate until you find B hadron as a mother or become top ancester
884 for (Int_t i=1; i<100; i++){ // check back to the 100 generation older
886 Int_t jLabel = partMother->GetFirstMother();
891 if (jLabel < 0){ // safety protection
892 AliDebug(1, "Stack label is negative, return\n");
896 // if there is an ancester
897 TParticle* grandMa = fStack->Particle(jLabel);
898 Int_t grandMaPDG = grandMa->GetPdgCode();
900 for (Int_t j=0; j<fNparents; j++){
901 if (abs(grandMaPDG)==fParentSelect[1][j]){
902 origin = kBeautyElse;
907 partMother = grandMa;
908 } // end of iteration
914 //_______________________________________________________________________________________________
915 Int_t AliHFEsecVtx::GetPDG(AliVTrack *track)
918 // get KF particle input pdg for mass hypothesis
923 if (fUseMCPID && HasMCData()){
924 pdgCode = GetMCPDG(track);
929 GetESDPID((AliESDtrack*)track, pid, prob);
931 case 0: pdgCode = 11; break;
932 case 1: pdgCode = 13; break;
933 case 2: pdgCode = 211; break;
934 case 3: pdgCode = 321; break;
935 case 4: pdgCode = 2212; break;
936 default: pdgCode = -1;
940 Int_t pid = ((AliAODTrack*)track)->GetMostProbablePID();
942 case 0: pdgCode = 11; break;
943 case 1: pdgCode = 13; break;
944 case 2: pdgCode = 211; break;
945 case 3: pdgCode = 321; break;
946 case 4: pdgCode = 2212; break;
947 default: pdgCode = -1;
954 //_______________________________________________________________________________________________
955 Int_t AliHFEsecVtx::GetMCPDG(AliVTrack *track)
958 // return mc pdg code
961 Int_t label = TMath::Abs(track->GetLabel());
964 if (IsAODanalysis()) {
965 AliAODMCParticle *mcpart = (AliAODMCParticle*)fMCArray->At(label);
966 if ( !mcpart ) return 0;
967 pdgCode = mcpart->GetPdgCode();
970 TParticle* mcpart = fStack->Particle(label);
971 if ( !mcpart ) return 0;
972 pdgCode = mcpart->GetPdgCode();
978 //______________________________________________________________________________
979 void AliHFEsecVtx::GetESDPID(AliESDtrack *track, Int_t &recpid, Double_t &recprob)
982 // calculate likehood for esd pid
991 Double_t probability[5];
993 // get probability of the diffenrent particle types
994 track->GetESDpid(probability);
996 // find most probable particle in ESD pid
997 // 0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
998 ipid = TMath::LocMax(5,probability);
999 max = TMath::MaxElement(5,probability);
1005 //_____________________________________________________________________________
1006 void AliHFEsecVtx::AddHFEpairToArray(const AliHFEpairs* const pair)
1009 // Add a HFE pair to the array
1012 Int_t n = HFEpairs()->GetEntriesFast();
1013 if(n!=fNoOfHFEpairs)AliError(Form("fNoOfHFEpairs != HFEpairs()->GetEntriesFast %i != %i \n", fNoOfHFEpairs, n));
1014 new((*HFEpairs())[n]) AliHFEpairs(*pair);
1017 //_____________________________________________________________________________
1018 TClonesArray *AliHFEsecVtx::HFEpairs()
1021 // Returns the list of HFE pairs
1025 fHFEpairs = new TClonesArray("AliHFEpairs", 200);
1030 //_____________________________________________________________________________
1031 void AliHFEsecVtx::DeleteHFEpairs()
1034 // Delete the list of HFE pairs
1038 fHFEpairs->Delete();
1043 //_____________________________________________________________________________
1044 void AliHFEsecVtx::InitHFEpairs()
1047 // Initialization should be done before make all possible pairs for a given electron candidate
1053 //_____________________________________________________________________________
1054 void AliHFEsecVtx::AddHFEsecvtxToArray(const AliHFEsecVtxs* const secvtx)
1057 // Add a HFE secondary vertex to the array
1060 Int_t n = HFEsecvtxs()->GetEntriesFast();
1061 if(n!=fNoOfHFEsecvtxs)AliError(Form("fNoOfHFEsecvtxs != HFEsecvtxs()->GetEntriesFast %i != %i \n", fNoOfHFEsecvtxs, n));
1062 new((*HFEsecvtxs())[n]) AliHFEsecVtxs(*secvtx);
1065 //_____________________________________________________________________________
1066 TClonesArray *AliHFEsecVtx::HFEsecvtxs()
1069 // Returns the list of HFE secvtx
1073 fHFEsecvtxs = new TClonesArray("AliHFEsecVtxs", 200);
1078 //_____________________________________________________________________________
1079 void AliHFEsecVtx::DeleteHFEsecvtxs()
1082 // Delete the list of HFE pairs
1086 fHFEsecvtxs->Delete();
1087 //delete fHFEsecvtx;
1091 //_____________________________________________________________________________
1092 void AliHFEsecVtx::InitHFEsecvtxs()
1095 // Initialization should be done
1098 fNoOfHFEsecvtxs = 0;
1101 //_______________________________________________________________________________________________
1102 Bool_t AliHFEsecVtx::SingleTrackCut(AliESDtrack* track) const
1104 //if (track->Pt() < 1.0) return kFALSE;
1105 //if (TMath::Abs(track->Eta()) > 0.9) return kFALSE;
1106 //if (!(track->GetStatus() & AliESDtrack::kITSrefit)) return kFALSE;
1107 //if (!(track->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
1108 if (!(TESTBIT(track->GetITSClusterMap(),0))) return kFALSE; // ask hit on the first pixel layer
1109 //if (!(TESTBIT(track->GetITSClusterMap(),0) | TESTBIT(track->GetITSClusterMap(),1))) return kFALSE;
1114 track->GetImpactParameters(dcaR,dcaZ);
1115 if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) < 0.005) return kFALSE;
1116 if (TMath::Abs(TMath::Sqrt(dcaR*dcaR + dcaZ*dcaZ)) > 0.3) return kFALSE;
1120 //____________________________________________________________
1121 void AliHFEsecVtx::MakeContainer(){
1127 const Int_t nDimPair=6;
1128 Int_t nBinPair[nDimPair] = {200, 500, 314, 2000, 2000, 11};
1129 //Int_t nBinPair[nDimPair] = {200, 500, 314, 2000, 2000, 11, 1000, 1000, 60, 60};
1130 const Double_t kInvmassmin = 0., kInvmassmax = 20.;
1131 const Double_t kKFChi2min = 0, kKFChi2max= 50;
1132 const Double_t kOpenanglemin = 0, kOpenanglemax = 3.14;
1133 const Double_t kSignedLxymin = -10, kSignedLxymax= 10;
1134 const Double_t kKFIPmin = -10, kKFIPmax= 10;
1135 const Double_t kPairCodemin = -1, kPairCodemax= 10;
1136 //const Double_t kDCAsigmin = 0, kDCAsigmax= 5;
1137 //const Double_t kPtmin = 0, kPtmax= 30;
1139 Double_t* binEdgesPair[nDimPair];
1140 for(Int_t ivar = 0; ivar < nDimPair; ivar++)
1141 binEdgesPair[ivar] = new Double_t[nBinPair[ivar] + 1];
1143 for(Int_t i=0; i<=nBinPair[0]; i++) binEdgesPair[0][i]=(Double_t)kInvmassmin + (kInvmassmax - kInvmassmin)/nBinPair[0]*(Double_t)i;
1144 for(Int_t i=0; i<=nBinPair[1]; i++) binEdgesPair[1][i]=(Double_t)kKFChi2min + (kKFChi2max - kKFChi2min)/nBinPair[1]*(Double_t)i;
1145 for(Int_t i=0; i<=nBinPair[2]; i++) binEdgesPair[2][i]=(Double_t)kOpenanglemin + (kOpenanglemax - kOpenanglemin)/nBinPair[2]*(Double_t)i;
1146 for(Int_t i=0; i<=nBinPair[3]; i++) binEdgesPair[3][i]=(Double_t)kSignedLxymin + (kSignedLxymax - kSignedLxymin)/nBinPair[3]*(Double_t)i;
1147 for(Int_t i=0; i<=nBinPair[4]; i++) binEdgesPair[4][i]=(Double_t)kKFIPmin + (kKFIPmax - kKFIPmin)/nBinPair[4]*(Double_t)i;
1148 for(Int_t i=0; i<=nBinPair[5]; i++) binEdgesPair[5][i]=(Double_t)kPairCodemin + (kPairCodemax - kPairCodemin)/nBinPair[5]*(Double_t)i;
1149 /*for(Int_t i=0; i<=nBinPair[6]; i++) binEdgesPair[6][i]=(Double_t)kDCAsigmin + (kDCAsigmax - kDCAsigmin)/nBinPair[6]*(Double_t)i;
1150 for(Int_t i=0; i<=nBinPair[7]; i++) binEdgesPair[7][i]=binEdgesPair[6][i];
1151 for(Int_t i=0; i<=nBinPair[8]; i++) binEdgesPair[8][i]=(Double_t)kPtmin + (kPtmax - kPtmin)/nBinPair[8]*(Double_t)i;
1152 for(Int_t i=0; i<=nBinPair[9]; i++) binEdgesPair[9][i]=binEdgesPair[8][i];*/
1154 fPairQA = new THnSparseF("pairQA", "QA for Pair; invmass[GeV/c^2]; KF chi2; opening angle; signed Lxy; KF ip; pair code", nDimPair, nBinPair);
1155 //fPairQA = new THnSparseF("pairQA", "QA for Pair; invmass[GeV/c^2]; KF chi2; opening angle; signed Lxy; KF ip; pair code; dca sig trk1; dca sig trk2; pt trk1; pt trk2 ", nDimPair, nBinPair);
1156 for(Int_t idim = 0; idim < nDimPair; idim++){
1157 fPairQA->SetBinEdges(idim, binEdgesPair[idim]);
1160 fSecVtxList->AddAt(fPairQA,0);
1162 const Int_t nDimSecvtx=6;
1163 Double_t* binEdgesSecvtx[nDimSecvtx];
1164 Int_t nBinSecvtx[nDimSecvtx] = {200, 500, 2000, 2000, 11, 3};
1165 const Double_t kNtrksmin = 0, kNtrksmax= 3;
1166 for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
1167 binEdgesSecvtx[ivar] = new Double_t[nBinSecvtx[ivar] + 1];
1169 for(Int_t i=0; i<=nBinSecvtx[0]; i++) binEdgesSecvtx[0][i]=binEdgesPair[0][i];
1170 for(Int_t i=0; i<=nBinSecvtx[1]; i++) binEdgesSecvtx[1][i]=binEdgesPair[1][i];
1171 for(Int_t i=0; i<=nBinSecvtx[2]; i++) binEdgesSecvtx[2][i]=binEdgesPair[3][i];
1172 for(Int_t i=0; i<=nBinSecvtx[3]; i++) binEdgesSecvtx[3][i]=binEdgesPair[4][i];
1173 for(Int_t i=0; i<=nBinSecvtx[4]; i++) binEdgesSecvtx[4][i]=binEdgesPair[5][i];
1174 for(Int_t i=0; i<=nBinSecvtx[5]; i++) binEdgesSecvtx[5][i]=(Double_t)kNtrksmin + (kNtrksmax - kNtrksmin)/nBinSecvtx[5]*(Double_t)i;
1176 fSecvtxQA = new THnSparseF("secvtxQA", "QA for Secvtx; invmass[GeV/c^2]; KF chi2; signed Lxy; KF ip; pair code; n tracks ", nDimSecvtx, nBinSecvtx);
1177 for(Int_t idim = 0; idim < nDimSecvtx; idim++){
1178 fSecvtxQA->SetBinEdges(idim, binEdgesSecvtx[idim]);
1181 fSecVtxList->AddAt(fSecvtxQA,1);
1182 for(Int_t ivar = 0; ivar < nDimPair; ivar++)
1183 delete binEdgesPair[ivar];
1184 for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
1185 delete binEdgesSecvtx[ivar];