From 5dc814ec1d20c5d32b173b9b9a4e8e6e37de4035 Mon Sep 17 00:00:00 2001 From: fprino Date: Thu, 4 Dec 2014 01:27:35 +0100 Subject: [PATCH] Methods to check Lc decays in the Kinematics --- .../AliAnalysisTaskCheckHFMCProd.cxx | 40 ++- .../AliAnalysisTaskCheckHFMCProd.h | 5 +- PWGHF/vertexingHF/AliVertexingHFUtils.cxx | 291 +++++++++++++++++- PWGHF/vertexingHF/AliVertexingHFUtils.h | 4 + 4 files changed, 334 insertions(+), 6 deletions(-) diff --git a/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx b/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx index eaa62ed742b..d66ea5c7efd 100644 --- a/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx +++ b/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx @@ -78,6 +78,7 @@ AliAnalysisTaskCheckHFMCProd::AliAnalysisTaskCheckHFMCProd() : AliAnalysisTaskSE fHistMotherID(0), fHistDSpecies(0), fHistBSpecies(0), + fHistLcDecayChan(0), fHistNcollHFtype(0), fHistEtaPhiPtGenEle(0), fHistEtaPhiPtGenPi(0), @@ -252,6 +253,7 @@ void AliAnalysisTaskCheckHFMCProd::UserCreateOutputObjects() { fHistYPtD0byDecChannel[1] = new TH2F("hyptD04","D0 - 4prong",40,0.,40.,20,-2.,2.); fHistYPtDplusbyDecChannel[0] = new TH2F("hyptDplusnonreson","Dplus - non reson",40,0.,40.,20,-2.,2.); fHistYPtDplusbyDecChannel[1] = new TH2F("hyptDplusreson","Dplus - reson via K0*",40,0.,40.,20,-2.,2.); + fHistYPtDplusbyDecChannel[2] = new TH2F("hyptDplusKKpi","Dplus -> KKpi",40,0.,40.,20,-2.,2.); fHistYPtDsbyDecChannel[0] = new TH2F("hyptDsphi","Ds - vis Phi",40,0.,40.,20,-2.,2.); fHistYPtDsbyDecChannel[1] = new TH2F("hyptDsk0st","Ds - via k0*",40,0.,40.,20,-2.,2.); @@ -267,7 +269,10 @@ void AliAnalysisTaskCheckHFMCProd::UserCreateOutputObjects() { fHistYPtDsbyDecChannel[ih]->SetMinimum(0); fOutput->Add(fHistYPtDsbyDecChannel[ih]); } - + fHistYPtDplusbyDecChannel[2]->Sumw2(); + fHistYPtDplusbyDecChannel[2]->SetMinimum(0); + fOutput->Add(fHistYPtDplusbyDecChannel[2]); + fHistOriginPrompt=new TH1F("hOriginPrompt","",100,0.,0.5); fHistOriginPrompt->Sumw2(); fHistOriginPrompt->SetMinimum(0); @@ -305,6 +310,18 @@ void AliAnalysisTaskCheckHFMCProd::UserCreateOutputObjects() { fHistBSpecies->GetXaxis()->SetBinLabel(10,"Lb-"); fHistBSpecies->SetMinimum(0); fOutput->Add(fHistBSpecies); + fHistLcDecayChan=new TH1F("hLcDecayChan","",9,-2.5,6.5); + fHistLcDecayChan->GetXaxis()->SetBinLabel(1,"Violates p cons"); + fHistLcDecayChan->GetXaxis()->SetBinLabel(2,"Other decay"); + fHistLcDecayChan->GetXaxis()->SetBinLabel(3,"Error"); + fHistLcDecayChan->GetXaxis()->SetBinLabel(4,"pK#pi non res"); + fHistLcDecayChan->GetXaxis()->SetBinLabel(5,"pK#pi via K*0"); + fHistLcDecayChan->GetXaxis()->SetBinLabel(6,"pK#pi via #Delta++"); + fHistLcDecayChan->GetXaxis()->SetBinLabel(7,"pK#pi via #Lambda1520"); + fHistLcDecayChan->GetXaxis()->SetBinLabel(8,"pK0s#rightarrowp#pi#pi"); + fHistLcDecayChan->GetXaxis()->SetBinLabel(9,"#pi#Lambda#rightarrowp#pi#pi"); + fHistLcDecayChan->SetMinimum(0); + fOutput->Add(fHistLcDecayChan); fHistNcollHFtype=new TH2F("hNcollHFtype","",5,-1.5,3.5,30,-0.5,29.5); fOutput->Add(fHistNcollHFtype); @@ -454,6 +471,13 @@ void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *) } nColl=lgen->GetEntries(); fHistNcollHFtype->Fill(typeHF,nColl); + }else{ + TString genTitle=mcEvent->GenEventHeader()->GetTitle(); + if(genTitle.Contains("bchadr")) typeHF=1; + else if(genTitle.Contains("chadr")) typeHF=0; + else if(genTitle.Contains("bele")) typeHF=3; + else if(genTitle.Contains("cele")) typeHF=2; + fHistNcollHFtype->Fill(typeHF,1.); } Int_t nParticles=stack->GetNtrack(); Double_t dNchdy = 0.; @@ -500,6 +524,10 @@ void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *) else if(absPdg==411){ iSpecies=1; iType=AliVertexingHFUtils::CheckDplusDecay(stack,i,dummy); + if(iType<0){ + Int_t iTypeKKpi=AliVertexingHFUtils::CheckDplusKKpiDecay(stack,i,dummy); + if(iTypeKKpi>0) iType=3; + } if(iType>0) iPart=1; } else if(absPdg==413){ @@ -514,7 +542,13 @@ void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *) } else if(absPdg==4122){ iSpecies=4; - iType=CheckLcDecay(i,stack); + iType=AliVertexingHFUtils::CheckLcpKpiDecay(stack,i,dummy); + if(iType<0){ + Int_t iTypeV0=AliVertexingHFUtils::CheckLcV0bachelorDecay(stack,i,dummy); + if(iTypeV0==1) iType=5; + if(iTypeV0==2) iType=6; + } + fHistLcDecayChan->Fill(iType); if(iType>=0) iPart=4; } if(iSpecies>=0) fHistYPtAllDecay[iSpecies]->Fill(part->Pt(),rapid); @@ -570,7 +604,7 @@ void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *) nCharmed++; if(iPart==0 && iType>0 && iType<=2){ fHistYPtD0byDecChannel[iType-1]->Fill(part->Pt(),rapid); - }else if(iPart==1 && iType>0 && iType<=2){ + }else if(iPart==1 && iType>0 && iType<=3){ fHistYPtDplusbyDecChannel[iType-1]->Fill(part->Pt(),rapid); }else if(iPart==3 && iType>0 && iType<=2){ fHistYPtDsbyDecChannel[iType-1]->Fill(part->Pt(),rapid); diff --git a/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.h b/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.h index 085e1b81007..6bc2eb18d29 100644 --- a/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.h +++ b/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.h @@ -87,13 +87,14 @@ class AliAnalysisTaskCheckHFMCProd : public AliAnalysisTaskSE { TH2F* fHistYPtPrompt[5]; //! histo of y vs. pt from prompt D0, D+, D*, Ds, Lc TH2F* fHistYPtFeeddown[5]; //! histo of y vs. pt from feeddown D0, D+, D*, Ds, Lc TH2F* fHistYPtD0byDecChannel[2]; //! histo of y vs. pt for D0->Kpi and D0->Kpipipi - TH2F* fHistYPtDplusbyDecChannel[2]; //! histo of y vs. pt for D+->Kpipi and D+->K0*pi + TH2F* fHistYPtDplusbyDecChannel[3]; //! histo of y vs. pt for D+->Kpipi and D+->K0*pi TH2F* fHistYPtDsbyDecChannel[2]; //! histo of y vs. pt for Ds->phipi and Ds->K0*K TH1F* fHistOriginPrompt; //! histo of D production point (prompt) TH1F* fHistOriginFeeddown; //! histo of D production point (feeddown) TH1F* fHistMotherID; //! histo of mother ID TH1F* fHistDSpecies; //! histo of D hadron species TH1F* fHistBSpecies; //! histo of B hadron species + TH1F* fHistLcDecayChan; //! histo of Lc decay modes TH2F* fHistNcollHFtype; //! histo of B hadron species TH3F* fHistEtaPhiPtGenEle; //! histo of generated electrons TH3F* fHistEtaPhiPtGenPi; //! histo of generated pions @@ -108,7 +109,7 @@ class AliAnalysisTaskCheckHFMCProd : public AliAnalysisTaskSE { AliESDtrackCuts *fESDtrackCuts; // track selection Bool_t fReadMC; - ClassDef(AliAnalysisTaskCheckHFMCProd,6); + ClassDef(AliAnalysisTaskCheckHFMCProd,7); }; diff --git a/PWGHF/vertexingHF/AliVertexingHFUtils.cxx b/PWGHF/vertexingHF/AliVertexingHFUtils.cxx index a8b9bd0a398..8480850de22 100644 --- a/PWGHF/vertexingHF/AliVertexingHFUtils.cxx +++ b/PWGHF/vertexingHF/AliVertexingHFUtils.cxx @@ -1050,6 +1050,95 @@ Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCPartic return -1; +} +//____________________________________________________________________________ +Int_t AliVertexingHFUtils::CheckDplusKKpiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ + // Checks the Ds decay channel. Returns 1 for Ds->phipi->KKpi, 2 for Ds->K0*K->KKpi, 3 for the non-resonant case + + if(label<0) return -1; + TParticle* mcPart = stack->Particle(label); + Int_t pdgD=mcPart->GetPdgCode(); + if(TMath::Abs(pdgD)!=411) return -1; + + Int_t nDau=mcPart->GetNDaughters(); + Int_t labelFirstDau = mcPart->GetDaughter(0); + Int_t nKaons=0; + Int_t nPions=0; + Double_t sumPxDau=0.; + Double_t sumPyDau=0.; + Double_t sumPzDau=0.; + Int_t nFoundKpi=0; + Bool_t isPhi=kFALSE; + Bool_t isk0st=kFALSE; + + if(nDau==3 || nDau==2){ + for(Int_t iDau=0; iDauParticle(indDau); + if(!dau) return -1; + Int_t pdgdau=dau->GetPdgCode(); + if(TMath::Abs(pdgdau)==321){ + nKaons++; + sumPxDau+=dau->Px(); + sumPyDau+=dau->Py(); + sumPzDau+=dau->Pz(); + arrayDauLab[nFoundKpi++]=indDau; + if(nFoundKpi>3) return -1; + }else if(TMath::Abs(pdgdau)==211){ + nPions++; + sumPxDau+=dau->Px(); + sumPyDau+=dau->Py(); + sumPzDau+=dau->Pz(); + arrayDauLab[nFoundKpi++]=indDau; + if(nFoundKpi>3) return -1; + }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==333){ + if(TMath::Abs(pdgdau)==313) isk0st=kTRUE; + if(TMath::Abs(pdgdau)==333) isPhi=kTRUE; + Int_t nResDau=dau->GetNDaughters(); + if(nResDau!=2) return -1; + Int_t indFirstResDau=dau->GetDaughter(0); + for(Int_t resDau=0; resDau<2; resDau++){ + Int_t indResDau=indFirstResDau+resDau; + if(indResDau<0) return -1; + TParticle* resdau=stack->Particle(indResDau); + if(!resdau) return -1; + Int_t pdgresdau=resdau->GetPdgCode(); + if(TMath::Abs(pdgresdau)==321){ + sumPxDau+=resdau->Px(); + sumPyDau+=resdau->Py(); + sumPzDau+=resdau->Pz(); + nKaons++; + arrayDauLab[nFoundKpi++]=indResDau; + if(nFoundKpi>3) return -1; + } + if(TMath::Abs(pdgresdau)==211){ + sumPxDau+=resdau->Px(); + sumPyDau+=resdau->Py(); + sumPzDau+=resdau->Pz(); + nPions++; + arrayDauLab[nFoundKpi++]=indResDau; + if(nFoundKpi>3) return -1; + } + } + }else{ + return -1; + } + } + if(nPions!=1) return -1; + if(nKaons!=2) return -1; + if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; + if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; + if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; + if(nDau==3) return 3; + else if(nDau==2){ + if(isk0st) return 2; + if(isPhi) return 1; + } + } + + return -1; + } //____________________________________________________________________________ Int_t AliVertexingHFUtils::CheckDsDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ @@ -1207,7 +1296,7 @@ Int_t AliVertexingHFUtils::CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle arrayDauLab[nFoundKpi++]=indResDau; if(nFoundKpi>3) return -1; } - } + } }else{ return -1; } @@ -1375,3 +1464,203 @@ Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCPartic } +//____________________________________________________________________________ +Int_t AliVertexingHFUtils::CheckLcpKpiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ + // Checks the Lc->pKpi decay channel. Returns 1 for non-resonant decays and 2, 3 or 4 for resonant ones, -1 in other cases + + if(label<0) return -1; + TParticle* mcPart = stack->Particle(label); + Int_t pdgD=mcPart->GetPdgCode(); + if(TMath::Abs(pdgD)!=4122) return -1; + + Int_t nDau=mcPart->GetNDaughters(); + Int_t labelFirstDau = mcPart->GetDaughter(0); + Int_t nKaons=0; + Int_t nPions=0; + Int_t nProtons=0; + Double_t sumPxDau=0.; + Double_t sumPyDau=0.; + Double_t sumPzDau=0.; + Int_t nFoundpKpi=0; + + Int_t codeRes=-1; + if(nDau==3 || nDau==2){ + for(Int_t iDau=0; iDauParticle(indDau); + if(!dau) return -1; + Int_t pdgdau=dau->GetPdgCode(); + if(TMath::Abs(pdgdau)==321){ + nKaons++; + sumPxDau+=dau->Px(); + sumPyDau+=dau->Py(); + sumPzDau+=dau->Pz(); + arrayDauLab[nFoundpKpi++]=indDau; + if(nFoundpKpi>3) return -1; + }else if(TMath::Abs(pdgdau)==211){ + nPions++; + sumPxDau+=dau->Px(); + sumPyDau+=dau->Py(); + sumPzDau+=dau->Pz(); + arrayDauLab[nFoundpKpi++]=indDau; + if(nFoundpKpi>3) return -1; + }else if(TMath::Abs(pdgdau)==2212){ + nProtons++; + sumPxDau+=dau->Px(); + sumPyDau+=dau->Py(); + sumPzDau+=dau->Pz(); + arrayDauLab[nFoundpKpi++]=indDau; + if(nFoundpKpi>3) return -1; + }else if(TMath::Abs(pdgdau)==313 || TMath::Abs(pdgdau)==3124 || + TMath::Abs(pdgdau)==2224){ + codeRes=TMath::Abs(pdgdau); + Int_t nResDau=dau->GetNDaughters(); + if(nResDau!=2) return -1; + Int_t indFirstResDau=dau->GetDaughter(0); + for(Int_t resDau=0; resDau<2; resDau++){ + Int_t indResDau=indFirstResDau+resDau; + if(indResDau<0) return -1; + TParticle* resdau=stack->Particle(indResDau); + if(!resdau) return -1; + Int_t pdgresdau=resdau->GetPdgCode(); + if(TMath::Abs(pdgresdau)==321){ + sumPxDau+=resdau->Px(); + sumPyDau+=resdau->Py(); + sumPzDau+=resdau->Pz(); + nKaons++; + arrayDauLab[nFoundpKpi++]=indResDau; + if(nFoundpKpi>3) return -1; + }else if(TMath::Abs(pdgresdau)==211){ + sumPxDau+=resdau->Px(); + sumPyDau+=resdau->Py(); + sumPzDau+=resdau->Pz(); + nPions++; + arrayDauLab[nFoundpKpi++]=indResDau; + if(nFoundpKpi>3) return -1; + }else if(TMath::Abs(pdgresdau)==2212){ + sumPxDau+=resdau->Px(); + sumPyDau+=resdau->Py(); + sumPzDau+=resdau->Pz(); + nProtons++; + arrayDauLab[nFoundpKpi++]=indResDau; + if(nFoundpKpi>3) return -1; + } + } + }else{ + return -1; + } + } + if(nPions!=1) return -1; + if(nKaons!=1) return -1; + if(nProtons!=1) return -1; + if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; + if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; + if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; + if(nDau==3) return 1; + else if(nDau==2){ + if(codeRes==313) return 2; + else if(codeRes==2224) return 3; + else if(codeRes==3124) return 4; + } + } + return -1; + +} + + +//____________________________________________________________________________ +Int_t AliVertexingHFUtils::CheckLcV0bachelorDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ + // Checks the Lc->V0+bachelor decay channel. Returns 1 for pK0s, 2 for piLambda, 3 for pK0l -1 in other cases + + if(label<0) return -1; + TParticle* mcPart = stack->Particle(label); + Int_t pdgD=mcPart->GetPdgCode(); + if(TMath::Abs(pdgD)!=4122) return -1; + + Int_t nDau=mcPart->GetNDaughters(); + Int_t labelFirstDau = mcPart->GetDaughter(0); + Int_t nPions=0; + Int_t nProtons=0; + Double_t sumPxDau=0.; + Double_t sumPyDau=0.; + Double_t sumPzDau=0.; + Int_t nFoundppi=0; + + Int_t codeV0=-1; + if(nDau==2){ + for(Int_t iDau=0; iDauParticle(indDau); + if(!dau) return -1; + Int_t pdgdau=dau->GetPdgCode(); + if(TMath::Abs(pdgdau)==211){ + nPions++; + sumPxDau+=dau->Px(); + sumPyDau+=dau->Py(); + sumPzDau+=dau->Pz(); + arrayDauLab[nFoundppi++]=indDau; + if(nFoundppi>3) return -1; + }else if(TMath::Abs(pdgdau)==2212){ + nProtons++; + sumPxDau+=dau->Px(); + sumPyDau+=dau->Py(); + sumPzDau+=dau->Pz(); + arrayDauLab[nFoundppi++]=indDau; + if(nFoundppi>3) return -1; + }else if(TMath::Abs(pdgdau)==311 || TMath::Abs(pdgdau)==3122){ + codeV0=TMath::Abs(pdgdau); + TParticle* v0=dau; + if(codeV0==311){ + Int_t nK0Dau=dau->GetNDaughters(); + if(nK0Dau!=1) return -1; + Int_t indK0s=dau->GetDaughter(0); + if(indK0s<0) return -1; + v0=stack->Particle(indK0s); + if(!v0) return -1; + Int_t pdgK0sl=v0->GetPdgCode(); + codeV0=TMath::Abs(pdgK0sl); + } + Int_t nV0Dau=v0->GetNDaughters(); + if(nV0Dau!=2) return -1; + Int_t indFirstV0Dau=v0->GetDaughter(0); + for(Int_t v0Dau=0; v0Dau<2; v0Dau++){ + Int_t indV0Dau=indFirstV0Dau+v0Dau; + if(indV0Dau<0) return -1; + TParticle* v0dau=stack->Particle(indV0Dau); + if(!v0dau) return -1; + Int_t pdgv0dau=v0dau->GetPdgCode(); + if(TMath::Abs(pdgv0dau)==211){ + sumPxDau+=v0dau->Px(); + sumPyDau+=v0dau->Py(); + sumPzDau+=v0dau->Pz(); + nPions++; + arrayDauLab[nFoundppi++]=indV0Dau; + if(nFoundppi>3) return -1; + }else if(TMath::Abs(pdgv0dau)==2212){ + sumPxDau+=v0dau->Px(); + sumPyDau+=v0dau->Py(); + sumPzDau+=v0dau->Pz(); + nProtons++; + arrayDauLab[nFoundppi++]=indV0Dau; + if(nFoundppi>3) return -1; + } + } + }else{ + return -1; + } + } + if(nPions!=2) return -1; + if(nProtons!=1) return -1; + if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -2; + if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -2; + if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -2; + if(codeV0==310) return 1; + else if(codeV0==3122) return 2; + } + return -1; + +} + + diff --git a/PWGHF/vertexingHF/AliVertexingHFUtils.h b/PWGHF/vertexingHF/AliVertexingHFUtils.h index ed5580dc00f..bb074abf295 100644 --- a/PWGHF/vertexingHF/AliVertexingHFUtils.h +++ b/PWGHF/vertexingHF/AliVertexingHFUtils.h @@ -111,10 +111,14 @@ class AliVertexingHFUtils : public TObject{ static Int_t CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab); static Int_t CheckDplusDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab); static Int_t CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab); + static Int_t CheckDplusKKpiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab); static Int_t CheckDsDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab); static Int_t CheckDsDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab); static Int_t CheckDstarDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab); static Int_t CheckDstarDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab); + static Int_t CheckLcpKpiDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab); + static Int_t CheckLcV0bachelorDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab); + private: -- 2.43.0