From ee9c9d9bb3584dedb8e24ca699cde818a99e0fd8 Mon Sep 17 00:00:00 2001 From: fprino Date: Fri, 8 Aug 2014 15:46:58 +0200 Subject: [PATCH] New methods to tag hadronic decays in the Kinematics.root --- .../AliAnalysisTaskCheckHFMCProd.cxx | 370 +----------- .../AliAnalysisTaskCheckHFMCProd.h | 4 - PWGHF/vertexingHF/AliVertexingHFUtils.cxx | 557 +++++++++++++++++- PWGHF/vertexingHF/AliVertexingHFUtils.h | 9 + 4 files changed, 586 insertions(+), 354 deletions(-) diff --git a/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx b/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx index af8b47a71f6..e69b9a832a0 100644 --- a/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx +++ b/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx @@ -11,6 +11,7 @@ #include "AliGenCocktailEventHeader.h" #include "AliGenEventHeader.h" #include "AliGenerator.h" +#include "AliVertexingHFUtils.h" #include "AliMultiplicity.h" #include #include @@ -490,25 +491,26 @@ void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *) Int_t iPart=-1; Int_t iType=0; Int_t iSpecies=-1; + Int_t dummy[4]; if(absPdg==421){ iSpecies=0; - iType=CheckD0Decay(i,stack); - if(iType>=0) iPart=0; + iType=AliVertexingHFUtils::CheckD0Decay(stack,i,dummy); + if(iType>0) iPart=0; } else if(absPdg==411){ iSpecies=1; - iType=CheckDplusDecay(i,stack); - if(iType>=0) iPart=1; + iType=AliVertexingHFUtils::CheckDplusDecay(stack,i,dummy); + if(iType>0) iPart=1; } else if(absPdg==413){ iSpecies=2; - iType=CheckDstarDecay(i,stack); - if(iType>=0) iPart=2; + iType=AliVertexingHFUtils::CheckDstarDecay(stack,i,dummy); + if(iType>0) iPart=2; } else if(absPdg==431){ iSpecies=3; - iType=CheckDsDecay(i,stack); - if(iType==0 || iType==1) iPart=3; + iType=AliVertexingHFUtils::CheckDsDecay(stack,i,dummy); + if(iType==1 || iType==2) iPart=3; } else if(absPdg==4122){ iSpecies=4; @@ -553,48 +555,12 @@ void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *) Double_t distz=part->Vz()-mcVert->GetZ(); Double_t distToVert=TMath::Sqrt(distx*distx+disty*disty+distz*distz); fHistMotherID->Fill(part->GetFirstMother()); - TParticle* runningpart=part; - Int_t iFromB=-1; - Int_t pdgmoth=-1; - if(fSearchUpToQuark){ - while(1){ - Int_t labmoth=runningpart->GetFirstMother(); - if(labmoth==-1) break; - TParticle *mot=(TParticle*)stack->Particle(labmoth); - pdgmoth=TMath::Abs(mot->GetPdgCode()); - if(pdgmoth==5){ - iFromB=1; - break; - }else if(pdgmoth==4){ - iFromB=0; - break; - } - runningpart=mot; - } - }else{ - iFromB=0; - while(1){ - Int_t labmoth=runningpart->GetFirstMother(); - if(labmoth==-1) break; - TParticle *mot=(TParticle*)stack->Particle(labmoth); - pdgmoth=TMath::Abs(mot->GetPdgCode()); - if(pdgmoth>=500 && pdgmoth<=599){ - iFromB=1; - break; - } - if(pdgmoth>=5000 && pdgmoth<=5999){ - iFromB=1; - break; - } - runningpart=mot; - } - } - - if(iFromB==0){ + Int_t iFromB=AliVertexingHFUtils::CheckOrigin(stack,part,fSearchUpToQuark); + if(iFromB==4){ fHistYPtPromptAllDecay[iSpecies]->Fill(part->Pt(),rapid); fHistOriginPrompt->Fill(distToVert); } - else if(iFromB==1){ + else if(iFromB==5){ fHistYPtFeeddownAllDecay[iSpecies]->Fill(part->Pt(),rapid); fHistOriginFeeddown->Fill(distToVert); } @@ -602,16 +568,16 @@ void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *) if(iPart<0) continue; if(iType<0) continue; nCharmed++; - if(iPart==0 && iType<=1){ - fHistYPtD0byDecChannel[iType]->Fill(part->Pt(),rapid); - }else if(iPart==1 && iType<=1){ - fHistYPtDplusbyDecChannel[iType]->Fill(part->Pt(),rapid); - }else if(iPart==3 && iType<=1){ - fHistYPtDsbyDecChannel[iType]->Fill(part->Pt(),rapid); + if(iPart==0 && iType>0 && iType<=2){ + fHistYPtD0byDecChannel[iType-1]->Fill(part->Pt(),rapid); + }else if(iPart==1 && iType>0 && iType<=2){ + fHistYPtDplusbyDecChannel[iType-1]->Fill(part->Pt(),rapid); + }else if(iPart==3 && iType>0 && iType<=2){ + fHistYPtDsbyDecChannel[iType-1]->Fill(part->Pt(),rapid); } - if(iFromB==0 && iPart>=0 && iPart<5) fHistYPtPrompt[iPart]->Fill(part->Pt(),rapid); - else if(iFromB==1 && iPart>=0 && iPart<5) fHistYPtFeeddown[iPart]->Fill(part->Pt(),rapid); + if(iFromB==4 && iPart>=0 && iPart<5) fHistYPtPrompt[iPart]->Fill(part->Pt(),rapid); + else if(iFromB==5 && iPart>=0 && iPart<5) fHistYPtFeeddown[iPart]->Fill(part->Pt(),rapid); } for(Int_t i=0; iParticle(labD0); - Int_t pdgdp=dp->GetPdgCode(); - Int_t nDau=dp->GetNDaughters(); - - if(nDau==2){ - Int_t nKaons=0; - Int_t nPions=0; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - TParticle* dau=(TParticle*)stack->Particle(iDau); - Int_t pdgdau=dau->GetPdgCode(); - if(TMath::Abs(pdgdau)==321){ - if(pdgdp>0 && pdgdau>0) return -1; - if(pdgdp<0 && pdgdau<0) return -1; - nKaons++; - }else if(TMath::Abs(pdgdau)==211){ - if(pdgdp<0 && pdgdau>0) return -1; - if(pdgdp>0 && pdgdau<0) return -1; - nPions++; - }else{ - return -1; - } - } - if(nPions!=1) return -1; - if(nKaons!=1) return -1; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - } - return 0; - } - - if(nDau==3 || nDau==4){ - Int_t nKaons=0; - Int_t nPions=0; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - TParticle* dau=(TParticle*)stack->Particle(iDau); - Int_t pdgdau=dau->GetPdgCode(); - if(TMath::Abs(pdgdau)==321){ - if(pdgdp>0 && pdgdau>0) return -1; - if(pdgdp<0 && pdgdau<0) return -1; - nKaons++; - }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){ - for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){ - if(resDau<0) return -1; - TParticle* resdau=(TParticle*)stack->Particle(resDau); - Int_t pdgresdau=resdau->GetPdgCode(); - if(TMath::Abs(pdgresdau)==321){ - if(pdgdp>0 && pdgresdau>0) return -1; - if(pdgdp<0 && pdgresdau<0) return -1; - nKaons++; - } - if(TMath::Abs(pdgresdau)==211){ - nPions++; - } - } - }else if(TMath::Abs(pdgdau)==211){ - nPions++; - }else{ - return -1; - } - } - if(nPions!=3) return -1; - if(nKaons!=1) return -1; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - } - return 1; - } - - return -1; -} - - -//______________________________________________________________________________ -Int_t AliAnalysisTaskCheckHFMCProd::CheckDplusDecay(Int_t labDplus, AliStack* stack) const{ - - if(labDplus<0) return -1; - TParticle* dp = (TParticle*)stack->Particle(labDplus); - Int_t pdgdp=dp->GetPdgCode(); - Int_t nDau=dp->GetNDaughters(); - - if(nDau==3){ - Int_t nKaons=0; - Int_t nPions=0; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - TParticle* dau=(TParticle*)stack->Particle(iDau); - Int_t pdgdau=dau->GetPdgCode(); - if(TMath::Abs(pdgdau)==321){ - if(pdgdp>0 && pdgdau>0) return -1; - if(pdgdp<0 && pdgdau<0) return -1; - nKaons++; - }else if(TMath::Abs(pdgdau)==211){ - if(pdgdp<0 && pdgdau>0) return -1; - if(pdgdp>0 && pdgdau<0) return -1; - nPions++; - }else{ - return -1; - } - } - if(nPions!=2) return -1; - if(nKaons!=1) return -1; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - } - return 0; - } - - if(nDau==2){ - Int_t nKaons=0; - Int_t nPions=0; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - TParticle* dau=(TParticle*)stack->Particle(iDau); - Int_t pdgdau=dau->GetPdgCode(); - if(TMath::Abs(pdgdau)==313){ - for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){ - if(resDau<0) return -1; - TParticle* resdau=(TParticle*)stack->Particle(resDau); - Int_t pdgresdau=resdau->GetPdgCode(); - if(TMath::Abs(pdgresdau)==321){ - if(pdgdp>0 && pdgresdau>0) return -1; - if(pdgdp<0 && pdgresdau<0) return -1; - nKaons++; - } - if(TMath::Abs(pdgresdau)==211){ - if(pdgdp<0 && pdgresdau>0) return -1; - if(pdgdp>0 && pdgresdau<0) return -1; - nPions++; - } - } - }else if(TMath::Abs(pdgdau)==211){ - if(pdgdp<0 && pdgdau>0) return -1; - if(pdgdp>0 && pdgdau<0) return -1; - nPions++; - }else{ - return -1; - } - } - if(nPions!=2) return -1; - if(nKaons!=1) return -1; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - } - return 1; - } - return -1; -} - -//______________________________________________________________________________ -Int_t AliAnalysisTaskCheckHFMCProd::CheckDsDecay(Int_t labDs, AliStack* stack) const{ - // Ds decay - if(labDs<0) return -1; - TParticle* dp = (TParticle*)stack->Particle(labDs); - Int_t pdgdp=dp->GetPdgCode(); - Int_t nDau=dp->GetNDaughters(); - - if(nDau==3){ - Int_t nKaons=0; - Int_t nPions=0; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - TParticle* dau=(TParticle*)stack->Particle(iDau); - Int_t pdgdau=dau->GetPdgCode(); - if(TMath::Abs(pdgdau)==321){ - nKaons++; - }else if(TMath::Abs(pdgdau)==211){ - if(pdgdp<0 && pdgdau>0) return -1; - if(pdgdp>0 && pdgdau<0) return -1; - nPions++; - }else{ - return -1; - } - } - if(nPions!=1) return -1; - if(nKaons!=2) return -1; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - } - return 2; - } - - if(nDau==2){ - Int_t nKaons=0; - Int_t nPions=0; - Bool_t isPhi=kFALSE; - Bool_t isk0st=kFALSE; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - TParticle* dau=(TParticle*)stack->Particle(iDau); - Int_t pdgdau=dau->GetPdgCode(); - if(TMath::Abs(pdgdau)==313){ - isk0st=kTRUE; - for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){ - if(resDau<0) return -1; - TParticle* resdau=(TParticle*)stack->Particle(resDau); - Int_t pdgresdau=resdau->GetPdgCode(); - if(TMath::Abs(pdgresdau)==321){ - nKaons++; - } - if(TMath::Abs(pdgresdau)==211){ - if(pdgdp<0 && pdgresdau>0) return -1; - if(pdgdp>0 && pdgresdau<0) return -1; - nPions++; - } - } - }else if(TMath::Abs(pdgdau)==333){ - isPhi=kTRUE; - for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){ - if(resDau<0) return -1; - TParticle* resdau=(TParticle*)stack->Particle(resDau); - if(!resdau) return -1; - Int_t pdgresdau=resdau->GetPdgCode(); - if(TMath::Abs(pdgresdau)==321){ - nKaons++; - }else{ - return -1; - } - } - }else if(TMath::Abs(pdgdau)==211){ - if(pdgdp<0 && pdgdau>0) return -1; - if(pdgdp>0 && pdgdau<0) return -1; - nPions++; - }else if(TMath::Abs(pdgdau)==321){ - nKaons++; - }else{ - return -1; - } - } - if(nPions!=1) return -1; - if(nKaons!=2) return -1; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - } - if(isk0st) return 1; - else if(isPhi) return 0; - else return 3; - } - return -1; -} - -//______________________________________________________________________________ -Int_t AliAnalysisTaskCheckHFMCProd::CheckDstarDecay(Int_t labDstar, AliStack* stack) const{ - - if(labDstar<0) return -1; - TParticle* dp = (TParticle*)stack->Particle(labDstar); - Int_t pdgdp=dp->GetPdgCode(); - Int_t nDau=dp->GetNDaughters(); - if(nDau!=2) return -1; - - Int_t nKaons=0; - Int_t nPions=0; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - TParticle* dau=(TParticle*)stack->Particle(iDau); - Int_t pdgdau=dau->GetPdgCode(); - if(TMath::Abs(pdgdau)==421){ - for(Int_t resDau=dau->GetFirstDaughter(); resDau<=dau->GetLastDaughter(); resDau++){ - if(resDau<0) return -1; - TParticle* resdau=(TParticle*)stack->Particle(resDau); - Int_t pdgresdau=resdau->GetPdgCode(); - if(TMath::Abs(pdgresdau)==321){ - if(pdgdp>0 && pdgresdau>0) return -1; - if(pdgdp<0 && pdgresdau<0) return -1; - nKaons++; - } - if(TMath::Abs(pdgresdau)==211){ - if(pdgdp<0 && pdgresdau>0) return -1; - if(pdgdp>0 && pdgresdau<0) return -1; - nPions++; - } - } - }else if(TMath::Abs(pdgdau)==211){ - if(pdgdp<0 && pdgdau>0) return -1; - if(pdgdp>0 && pdgdau<0) return -1; - nPions++; - }else{ - return -1; - } - } - if(nPions!=2) return -1; - if(nKaons!=1) return -1; - for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){ - if(iDau<0) return -1; - } - return 0; - -} - //______________________________________________________________________________ Int_t AliAnalysisTaskCheckHFMCProd::CheckLcDecay(Int_t labLc, AliStack* stack) const{ if(labLc<0) return -1; diff --git a/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.h b/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.h index 66758b4af4f..085e1b81007 100644 --- a/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.h +++ b/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.h @@ -40,10 +40,6 @@ class AliAnalysisTaskCheckHFMCProd : public AliAnalysisTaskSE { virtual void UserCreateOutputObjects(); virtual void Terminate(Option_t *option); - Int_t CheckD0Decay(Int_t labD0, AliStack* stack) const; - Int_t CheckDplusDecay(Int_t labDplus, AliStack* stack) const; - Int_t CheckDsDecay(Int_t labDs, AliStack* stack) const; - Int_t CheckDstarDecay(Int_t labDstar, AliStack* stack) const; Int_t CheckLcDecay(Int_t labLc, AliStack* stack) const; void SetSearchUpToQuark(Bool_t opt){fSearchUpToQuark=opt;}; diff --git a/PWGHF/vertexingHF/AliVertexingHFUtils.cxx b/PWGHF/vertexingHF/AliVertexingHFUtils.cxx index d35aeb598a4..78726f87eeb 100644 --- a/PWGHF/vertexingHF/AliVertexingHFUtils.cxx +++ b/PWGHF/vertexingHF/AliVertexingHFUtils.cxx @@ -20,6 +20,8 @@ #include #include #include +#include +#include "AliStack.h" #include "AliAODEvent.h" #include "AliAODMCHeader.h" #include "AliGenEventHeader.h" @@ -559,6 +561,38 @@ Bool_t AliVertexingHFUtils::HasCascadeCandidateAnyDaughInjected(AliAODRecoCascad return kFALSE; } //____________________________________________________________________________ +Int_t AliVertexingHFUtils::CheckOrigin(AliStack* stack, TParticle *mcPart, Bool_t searchUpToQuark){ + // checking whether the mother of the particles come from a charm or a bottom quark + + Int_t pdgGranma = 0; + Int_t mother = 0; + mother = mcPart->GetFirstMother(); + Int_t istep = 0; + Int_t abspdgGranma =0; + Bool_t isFromB=kFALSE; + Bool_t isQuarkFound=kFALSE; + while (mother >0 ){ + istep++; + TParticle* mcGranma = stack->Particle(mother); + if (mcGranma){ + pdgGranma = mcGranma->GetPdgCode(); + abspdgGranma = TMath::Abs(pdgGranma); + if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){ + isFromB=kTRUE; + } + if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE; + mother = mcGranma->GetFirstMother(); + }else{ + printf("CheckOrigin: Failed casting the mother particle!"); + break; + } + } + if(searchUpToQuark && !isQuarkFound) return 0; + if(isFromB) return 5; + else return 4; + +} +//____________________________________________________________________________ Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark){ // checking whether the mother of the particles come from a charm or a bottom quark @@ -590,10 +624,120 @@ Int_t AliVertexingHFUtils::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle * else return 4; } + //____________________________________________________________________________ -Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){ +Int_t AliVertexingHFUtils::CheckD0Decay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases + if(label<0) return -1; + TParticle* mcPart = stack->Particle(label); + Int_t pdgD=mcPart->GetPdgCode(); + if(TMath::Abs(pdgD)!=421) return -1; + + Int_t nDau=mcPart->GetNDaughters(); + + if(nDau==2){ + Int_t daughter0 = mcPart->GetDaughter(0); + Int_t daughter1 = mcPart->GetDaughter(1); + TParticle* mcPartDaughter0 = stack->Particle(daughter0); + TParticle* mcPartDaughter1 = stack->Particle(daughter1); + if(!mcPartDaughter0 || !mcPartDaughter1) return -1; + arrayDauLab[0]=daughter0; + arrayDauLab[1]=daughter1; + Int_t pdgCode0=mcPartDaughter0->GetPdgCode(); + Int_t pdgCode1=mcPartDaughter1->GetPdgCode(); + if(!(TMath::Abs(pdgCode0)==321 && TMath::Abs(pdgCode1)==211) && + !(TMath::Abs(pdgCode0)==211 && TMath::Abs(pdgCode1)==321)){ + return -1; + } + if(TMath::Abs(pdgCode0)==321 && (pdgD*pdgCode0)>0) return -1; + if(TMath::Abs(pdgCode1)==321 && (pdgD*pdgCode1)>0) return -1; + if((pdgCode0*pdgCode1)>0) return -1; + Double_t sumPxDau=mcPartDaughter0->Px()+mcPartDaughter1->Px(); + Double_t sumPyDau=mcPartDaughter0->Py()+mcPartDaughter1->Py(); + Double_t sumPzDau=mcPartDaughter0->Pz()+mcPartDaughter1->Pz(); + 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; + return 1; + } + + if(nDau==3 || nDau==4){ + Int_t nKaons=0; + Int_t nPions=0; + Double_t sumPxDau=0.; + Double_t sumPyDau=0.; + Double_t sumPzDau=0.; + Int_t nFoundKpi=0; + Int_t labelFirstDau = mcPart->GetDaughter(0); + for(Int_t iDau=0; iDauParticle(indDau); + if(!dau) return -1; + Int_t pdgdau=dau->GetPdgCode(); + if(TMath::Abs(pdgdau)==321){ + if(pdgD>0 && pdgdau>0) return -1; + if(pdgD<0 && pdgdau<0) return -1; + nKaons++; + sumPxDau+=dau->Px(); + sumPyDau+=dau->Py(); + sumPzDau+=dau->Pz(); + arrayDauLab[nFoundKpi++]=indDau; + if(nFoundKpi>4) return -1; + }else if(TMath::Abs(pdgdau)==211){ + nPions++; + sumPxDau+=dau->Px(); + sumPyDau+=dau->Py(); + sumPzDau+=dau->Pz(); + arrayDauLab[nFoundKpi++]=indDau; + if(nFoundKpi>4) return -1; + }else if(TMath::Abs(pdgdau)==113 || TMath::Abs(pdgdau)==313){ + 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){ + if(pdgD>0 && pdgresdau>0) return -1; + if(pdgD<0 && pdgresdau<0) return -1; + nKaons++; + sumPxDau+=resdau->Px(); + sumPyDau+=resdau->Py(); + sumPzDau+=resdau->Pz(); + arrayDauLab[nFoundKpi++]=indResDau; + if(nFoundKpi>4) return -1; + } + if(TMath::Abs(pdgresdau)==211){ + nPions++; + sumPxDau+=resdau->Px(); + sumPyDau+=resdau->Py(); + sumPzDau+=resdau->Pz(); + arrayDauLab[nFoundKpi++]=indResDau; + if(nFoundKpi>4) return -1; + } + } + }else{ + return -1; + } + } + if(nPions!=3) return -1; + if(nKaons!=1) return -1; + if(TMath::Abs(mcPart->Px()-sumPxDau)>0.001) return -1; + if(TMath::Abs(mcPart->Py()-sumPyDau)>0.001) return -1; + if(TMath::Abs(mcPart->Pz()-sumPzDau)>0.001) return -1; + return 2; + } + return -1; +} + +//____________________________________________________________________________ +Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){ + // Checks the D0 decay channel. Returns 1 for the D0->Kpi case, 2 for the D0->Kpipipi case, -1 in other cases Int_t pdgD=mcPart->GetPdgCode(); if(TMath::Abs(pdgD)!=421) return -1; @@ -698,6 +842,93 @@ Int_t AliVertexingHFUtils::CheckD0Decay(TClonesArray* arrayMC, AliAODMCParticle } return -1; } +//____________________________________________________________________________ +Int_t AliVertexingHFUtils::CheckDplusDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ + // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases + + 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; + + 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){ + if(pdgD*pdgdau>0) return -1; + nKaons++; + sumPxDau+=dau->Px(); + sumPyDau+=dau->Py(); + sumPzDau+=dau->Pz(); + arrayDauLab[nFoundKpi++]=indDau; + if(nFoundKpi>3) return -1; + }else if(TMath::Abs(pdgdau)==211){ + if(pdgD*pdgdau<0) return -1; + nPions++; + sumPxDau+=dau->Px(); + sumPyDau+=dau->Py(); + sumPzDau+=dau->Pz(); + arrayDauLab[nFoundKpi++]=indDau; + if(nFoundKpi>3) return -1; + }else if(TMath::Abs(pdgdau)==313){ + 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){ + if(pdgD*pdgresdau>0) return -1; + sumPxDau+=resdau->Px(); + sumPyDau+=resdau->Py(); + sumPzDau+=resdau->Pz(); + nKaons++; + arrayDauLab[nFoundKpi++]=indResDau; + if(nFoundKpi>3) return -1; + } + if(TMath::Abs(pdgresdau)==211){ + if(pdgD*pdgresdau<0) return -1; + sumPxDau+=resdau->Px(); + sumPyDau+=resdau->Py(); + sumPzDau+=resdau->Pz(); + nPions++; + arrayDauLab[nFoundKpi++]=indResDau; + if(nFoundKpi>3) return -1; + } + } + }else{ + return -1; + } + } + if(nPions!=2) return -1; + if(nKaons!=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) return 2; + } + + return -1; + +} + //____________________________________________________________________________ Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){ // Checks the Dplus decay channel. Returns 1 for the non-resonant case, 2 for the resonant case, -1 in other cases @@ -782,3 +1013,327 @@ Int_t AliVertexingHFUtils::CheckDplusDecay(TClonesArray* arrayMC, AliAODMCPartic return -1; } +//____________________________________________________________________________ +Int_t AliVertexingHFUtils::CheckDsDecay(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)!=431) 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(TClonesArray* arrayMC, AliAODMCParticle *mcPart, 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 + + Int_t pdgD=mcPart->GetPdgCode(); + if(TMath::Abs(pdgD)!=431) 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; iDau(arrayMC->At(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; + AliAODMCParticle* resdau=dynamic_cast(arrayMC->At(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::CheckDstarDecay(AliStack* stack, Int_t label, Int_t* arrayDauLab){ + // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases + + if(label<0) return -1; + TParticle* mcPart = stack->Particle(label); + Int_t pdgD=mcPart->GetPdgCode(); + if(TMath::Abs(pdgD)!=413) return -1; + + Int_t nDau=mcPart->GetNDaughters(); + if(nDau!=2) return -1; + + 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; + + for(Int_t iDau=0; iDauParticle(indDau); + if(!dau) return -1; + Int_t pdgdau=dau->GetPdgCode(); + if(TMath::Abs(pdgdau)==421){ + 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){ + if(pdgD*pdgresdau>0) return -1; + sumPxDau+=resdau->Px(); + sumPyDau+=resdau->Py(); + sumPzDau+=resdau->Pz(); + nKaons++; + arrayDauLab[nFoundKpi++]=indResDau; + if(nFoundKpi>3) return -1; + } + if(TMath::Abs(pdgresdau)==211){ + if(pdgD*pdgresdau<0) return -1; + sumPxDau+=resdau->Px(); + sumPyDau+=resdau->Py(); + sumPzDau+=resdau->Pz(); + nPions++; + arrayDauLab[nFoundKpi++]=indResDau; + if(nFoundKpi>3) return -1; + } + } + }else if(TMath::Abs(pdgdau)==211){ + if(pdgD*pdgdau<0) return -1; + nPions++; + sumPxDau+=dau->Px(); + sumPyDau+=dau->Py(); + sumPzDau+=dau->Pz(); + arrayDauLab[nFoundKpi++]=indDau; + if(nFoundKpi>3) return -1; + } + } + + if(nPions!=2) return -1; + if(nKaons!=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; + return 1; + +} + +//____________________________________________________________________________ +Int_t AliVertexingHFUtils::CheckDstarDecay(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Int_t* arrayDauLab){ + // Checks the Dstar decay channel. Returns 1 for D*->D0pi->Kpipi, -1 in other cases + + Int_t pdgD=mcPart->GetPdgCode(); + if(TMath::Abs(pdgD)!=413) return -1; + + Int_t nDau=mcPart->GetNDaughters(); + if(nDau!=2) return -1; + + 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; + + for(Int_t iDau=0; iDau(arrayMC->At(indDau)); + if(!dau) return -1; + Int_t pdgdau=dau->GetPdgCode(); + if(TMath::Abs(pdgdau)==421){ + 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; + AliAODMCParticle* resdau=dynamic_cast(arrayMC->At(indResDau)); + if(!resdau) return -1; + Int_t pdgresdau=resdau->GetPdgCode(); + if(TMath::Abs(pdgresdau)==321){ + if(pdgD*pdgresdau>0) return -1; + sumPxDau+=resdau->Px(); + sumPyDau+=resdau->Py(); + sumPzDau+=resdau->Pz(); + nKaons++; + arrayDauLab[nFoundKpi++]=indResDau; + if(nFoundKpi>3) return -1; + } + if(TMath::Abs(pdgresdau)==211){ + if(pdgD*pdgresdau<0) return -1; + sumPxDau+=resdau->Px(); + sumPyDau+=resdau->Py(); + sumPzDau+=resdau->Pz(); + nPions++; + arrayDauLab[nFoundKpi++]=indResDau; + if(nFoundKpi>3) return -1; + } + } + }else if(TMath::Abs(pdgdau)==211){ + if(pdgD*pdgdau<0) return -1; + nPions++; + sumPxDau+=dau->Px(); + sumPyDau+=dau->Py(); + sumPzDau+=dau->Pz(); + arrayDauLab[nFoundKpi++]=indDau; + if(nFoundKpi>3) return -1; + } + } + + if(nPions!=2) return -1; + if(nKaons!=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; + return 1; + +} + diff --git a/PWGHF/vertexingHF/AliVertexingHFUtils.h b/PWGHF/vertexingHF/AliVertexingHFUtils.h index 772330025ad..c68797625a9 100644 --- a/PWGHF/vertexingHF/AliVertexingHFUtils.h +++ b/PWGHF/vertexingHF/AliVertexingHFUtils.h @@ -20,11 +20,13 @@ #include "AliAODRecoDecayHF.h" #include "AliAODRecoCascadeHF.h" +class AliStack; class AliAODMCParticle; class AliAODMCHeader; class AliGenEventHeader; class AliAODEvent; class TProfile; +class TParticle; class TClonesArray; class TH1F; class TH2F; @@ -97,8 +99,15 @@ class AliVertexingHFUtils : public TObject{ static Double_t GetCorrectedNtracklets(TProfile* estimatorAvg, Double_t uncorrectedNacc, Double_t vtxZ, Double_t refMult); static Int_t CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE); + static Int_t CheckOrigin(AliStack* stack, TParticle *mcPart, Bool_t searchUpToQuark=kTRUE); + static Int_t CheckD0Decay(AliStack* stack, Int_t label, Int_t* arrayDauLab); 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 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); private: -- 2.43.0