New methods to tag hadronic decays in the Kinematics.root
authorfprino <prino@to.infn.it>
Fri, 8 Aug 2014 13:46:58 +0000 (15:46 +0200)
committerfprino <prino@to.infn.it>
Fri, 8 Aug 2014 13:47:25 +0000 (15:47 +0200)
PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx
PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.h
PWGHF/vertexingHF/AliVertexingHFUtils.cxx
PWGHF/vertexingHF/AliVertexingHFUtils.h

index af8b47a..e69b9a8 100644 (file)
@@ -11,6 +11,7 @@
 #include "AliGenCocktailEventHeader.h"
 #include "AliGenEventHeader.h"
 #include "AliGenerator.h"
+#include "AliVertexingHFUtils.h"
 #include "AliMultiplicity.h"
 #include <TParticle.h>
 #include <TSystem.h>
@@ -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; i<nTracks; i++){
@@ -654,300 +620,6 @@ void AliAnalysisTaskCheckHFMCProd::Terminate(Option_t */*option*/)
 
 
 //______________________________________________________________________________
-Int_t AliAnalysisTaskCheckHFMCProd::CheckD0Decay(Int_t labD0, AliStack* stack) const{
-
-  if(labD0<0) return -1;
-  TParticle* dp = (TParticle*)stack->Particle(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;
   TParticle* dp = (TParticle*)stack->Particle(labLc);
index 66758b4..085e1b8 100644 (file)
@@ -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;};
 
index d35aeb5..78726f8 100644 (file)
@@ -20,6 +20,8 @@
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TF1.h>
+#include <TParticle.h>
+#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; iDau<nDau; iDau++){
+      Int_t indDau = labelFirstDau+iDau;
+      if(indDau<0) return -1;
+      TParticle* dau=stack->Particle(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;
@@ -699,6 +843,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; iDau<nDau; iDau++){
+      Int_t indDau = labelFirstDau+iDau;
+      if(indDau<0) return -1;
+      TParticle* dau=stack->Particle(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; iDau<nDau; iDau++){
+      Int_t indDau = labelFirstDau+iDau;
+      if(indDau<0) return -1;
+      TParticle* dau=stack->Particle(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<nDau; iDau++){
+      Int_t indDau = labelFirstDau+iDau;
+      if(indDau<0) return -1;
+      AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(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<AliAODMCParticle*>(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; iDau<nDau; iDau++){
+    Int_t indDau = labelFirstDau+iDau;
+    if(indDau<0) return -1;
+    TParticle* dau=stack->Particle(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<nDau; iDau++){
+    Int_t indDau = labelFirstDau+iDau;
+    if(indDau<0) return -1;
+    AliAODMCParticle* dau=dynamic_cast<AliAODMCParticle*>(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<AliAODMCParticle*>(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;
+  
+}
+
index 7723300..c687976 100644 (file)
 #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: