New task to do D2H QA checks on MC productions
authorprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 3 Apr 2012 22:14:14 +0000 (22:14 +0000)
committerprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 3 Apr 2012 22:14:14 +0000 (22:14 +0000)
PWGHF/CMakelibPWGHFvertexingHF.pkg
PWGHF/PWGHFvertexingHFLinkDef.h
PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx [new file with mode: 0644]
PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.h [new file with mode: 0644]
PWGHF/vertexingHF/macros/AddHFMCCheck.C [new file with mode: 0644]
PWGHF/vertexingHF/macros/PlotOutputMCCheck.C [new file with mode: 0644]

index ec43124..a8d1a88 100644 (file)
@@ -55,7 +55,8 @@ set ( SRCS
     vertexingHF/AliAnalysisTaskSELambdac.cxx 
     vertexingHF/AliAnalysisTaskSED0Mass.cxx 
     vertexingHF/AliAnalysisTaskSECharmFraction.cxx 
-    vertexingHF/AliAnalysisTaskSEDvsMultiplicity.cxx 
+    vertexingHF/AliAnalysisTaskSEDvsMultiplicity.cxx
+    vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx
     vertexingHF/AliCFVertexingHF.cxx 
     vertexingHF/AliCFVertexingHF2Prong.cxx 
     vertexingHF/AliCFVertexingHF3Prong.cxx 
index e5fdbc6..1177852 100644 (file)
@@ -34,6 +34,7 @@
 #pragma link C++ class AliAnalysisTaskSED0Mass+;
 #pragma link C++ class AliAnalysisTaskSECharmFraction+;
 #pragma link C++ class AliAnalysisTaskSEDvsMultiplicity+;
+#pragma link C++ class AliAnalysisTaskCheckHFMCProd+;
 #pragma link C++ class AliCFVertexingHF+;
 #pragma link C++ class AliCFVertexingHF2Prong+;
 #pragma link C++ class AliCFVertexingHF3Prong+;
diff --git a/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx b/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.cxx
new file mode 100644 (file)
index 0000000..6ccca48
--- /dev/null
@@ -0,0 +1,756 @@
+#include "AliAnalysisTaskSE.h"
+#include "AliAnalysisManager.h"
+#include "AliAnalysisDataContainer.h"
+#include "AliESDEvent.h"
+#include "AliStack.h"
+#include "AliCentrality.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliMultiplicity.h"
+#include <TParticle.h>
+#include <TSystem.h>
+#include <TTree.h>
+#include <TNtuple.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TChain.h>
+#include "AliESDInputHandlerRP.h"
+#include "AliAnalysisTaskCheckHFMCProd.h"
+
+/**************************************************************************
+ * Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id$ */ 
+
+//*************************************************************************
+// Implementation of class AliAnalysisTaskCheckHFMCProd
+// AliAnalysisTask to check MC production at ESD+Kine level
+// 
+//
+// Authors: F. Prino, prino@to.infn.it
+//          
+//*************************************************************************
+
+ClassImp(AliAnalysisTaskCheckHFMCProd)
+//______________________________________________________________________________
+AliAnalysisTaskCheckHFMCProd::AliAnalysisTaskCheckHFMCProd() : AliAnalysisTaskSE("HFMCChecks"), 
+  fOutput(0),
+  fHistoNEvents(0),
+  fHistoTracks(0),
+  fHistoSelTracks(0),
+  fHistoTracklets(0),
+  fHistoSPD3DVtxX(0),
+  fHistoSPD3DVtxY(0),
+  fHistoSPD3DVtxZ(0),
+  fHistoSPDZVtxX(0),
+  fHistoSPDZVtxY(0),
+  fHistoSPDZVtxZ(0),
+  fHistoTRKVtxX(0),
+  fHistoTRKVtxY(0),
+  fHistoTRKVtxZ(0),
+  fHistoNcharmed(0),
+  fHistoNbVsNc(0),
+  fPbPb(kFALSE),
+  fReadMC(kTRUE)
+{
+  //
+  DefineInput(0, TChain::Class());
+  DefineOutput(1, TList::Class());
+}
+
+
+//___________________________________________________________________________
+AliAnalysisTaskCheckHFMCProd::~AliAnalysisTaskCheckHFMCProd(){
+  //
+  if (AliAnalysisManager::GetAnalysisManager()->IsProofMode()) return;
+  if (fOutput) {
+    delete fOutput;
+    fOutput = 0;
+  }
+}
+   
+//___________________________________________________________________________
+void AliAnalysisTaskCheckHFMCProd::UserCreateOutputObjects() {
+  // create output histos
+
+  fOutput = new TList();
+  fOutput->SetOwner();
+  fOutput->SetName("OutputHistos");
+
+  fHistoNEvents = new TH1F("hNEvents", "Number of processed events",3,-0.5,2.5);
+  fHistoNEvents->Sumw2();
+  fHistoNEvents->SetMinimum(0);
+  fOutput->Add(fHistoNEvents);
+
+  Double_t maxMult=100.;
+  if(fPbPb) maxMult=10000.;
+  fHistoTracks = new TH1F("hTracks","",100,0.,maxMult*2);
+  fHistoTracks->Sumw2();
+  fOutput->Add(fHistoTracks);
+  fHistoSelTracks = new TH1F("hSelTracks","",100,0.,maxMult);
+  fHistoSelTracks->Sumw2();
+  fOutput->Add(fHistoSelTracks);
+  fHistoTracklets = new TH1F("hTracklets","",100,0.,maxMult);
+  fHistoTracklets->Sumw2();
+  fOutput->Add(fHistoTracklets);
+
+  fHistoSPD3DVtxX = new TH1F("hSPD3DvX","",100,-1.,1.);
+  fHistoSPD3DVtxX->Sumw2();
+  fOutput->Add(fHistoSPD3DVtxX);
+  fHistoSPD3DVtxY = new TH1F("hSPD3DvY","",100,-1.,1.);
+  fHistoSPD3DVtxY->Sumw2();
+  fOutput->Add(fHistoSPD3DVtxY);
+  fHistoSPD3DVtxZ = new TH1F("hSPD3DvZ","",100,-15.,15.);
+  fHistoSPD3DVtxZ->Sumw2();
+  fOutput->Add(fHistoSPD3DVtxZ);
+
+  fHistoSPDZVtxX = new TH1F("hSPDZvX","",100,-1.,1.);
+  fHistoSPDZVtxX->Sumw2();
+  fOutput->Add(fHistoSPDZVtxX);
+  fHistoSPDZVtxY = new TH1F("hSPDZvY","",100,-1.,1.);
+  fHistoSPDZVtxY->Sumw2();
+  fOutput->Add(fHistoSPDZVtxY);
+  fHistoSPDZVtxZ = new TH1F("hSPDZvZ","",100,-15.,15.);
+  fHistoSPDZVtxZ->Sumw2();
+  fOutput->Add(fHistoSPDZVtxZ);
+
+
+  fHistoTRKVtxX = new TH1F("hTRKvX","",100,-1.,1.);
+  fHistoTRKVtxX->Sumw2();
+  fOutput->Add(fHistoTRKVtxX);
+  fHistoTRKVtxY = new TH1F("hTRKvY","",100,-1.,1.);
+  fHistoTRKVtxY->Sumw2();
+  fOutput->Add(fHistoTRKVtxY);
+  fHistoTRKVtxZ = new TH1F("hTRKvZ","",100,-15.,15.);
+  fHistoTRKVtxZ->Sumw2();
+  fOutput->Add(fHistoTRKVtxZ);
+
+  Int_t nBinscb=11;
+  if(fPbPb) nBinscb=200;
+  Double_t maxncn=nBinscb-0.5;
+  fHistoNcharmed = new TH2F("hncharmed","",100,0.,maxMult,nBinscb,-0.5,maxncn);
+  fHistoNcharmed->Sumw2();
+  fOutput->Add(fHistoNcharmed);
+  fHistoNbVsNc = new TH2F("hnbvsnc","",nBinscb,-0.5,maxncn,nBinscb,-0.5,maxncn);
+  fHistoNbVsNc->Sumw2();
+  fOutput->Add(fHistoNbVsNc);
+
+  fHistYPtPrompt[0] = new TH2F("hyptd0prompt","D0 - Prompt",20,0.,20.,20,-2.,2.);
+  fHistYPtPrompt[1] = new TH2F("hyptdplusprompt","Dplus - Prompt",20,0.,20.,20,-2.,2.);
+  fHistYPtPrompt[2] = new TH2F("hyptdstarprompt","Dstar - Prompt",20,0.,20.,20,-2.,2.);
+  fHistYPtPrompt[3] = new TH2F("hyptdsprompt","Ds - Prompt",20,0.,20.,20,-2.,2.);
+  fHistYPtPrompt[4] = new TH2F("hyptlcprompt","Lc - Prompt",20,0.,20.,20,-2.,2.);
+
+  fHistYPtFeeddown[0] = new TH2F("hyptd0feeddown","D0 - Feeddown",20,0.,20.,20,-2.,2.);
+  fHistYPtFeeddown[1] = new TH2F("hyptdplusfeeddown","Dplus - Feeddown",20,0.,20.,20,-2.,2.);
+  fHistYPtFeeddown[2] = new TH2F("hyptdstarfeedown","Dstar - Feeddown",20,0.,20.,20,-2.,2.);
+  fHistYPtFeeddown[3] = new TH2F("hyptdsfeedown","Ds - Feeddown",20,0.,20.,20,-2.,2.);
+  fHistYPtFeeddown[4] = new TH2F("hyptlcfeedown","Lc - Feeddown",20,0.,20.,20,-2.,2.);
+
+  for(Int_t ih=0; ih<5; ih++){
+    fHistYPtPrompt[ih]->Sumw2();
+    fHistYPtPrompt[ih]->SetMinimum(0);
+    fOutput->Add(fHistYPtPrompt[ih]);
+    fHistYPtFeeddown[ih]->Sumw2();
+    fHistYPtFeeddown[ih]->SetMinimum(0);
+    fOutput->Add(fHistYPtFeeddown[ih]);
+  }
+
+  fHistYPtD0byDecChannel[0] = new TH2F("hyptD02","D0 - 2prong",20,0.,20.,20,-2.,2.);
+  fHistYPtD0byDecChannel[1] = new TH2F("hyptD04","D0 - 4prong",20,0.,20.,20,-2.,2.);
+  fHistYPtDplusbyDecChannel[0] = new TH2F("hyptDplusnonreson","Dplus - non reson",20,0.,20.,20,-2.,2.);
+  fHistYPtDplusbyDecChannel[1] = new TH2F("hyptDplusreson","Dplus - reson via K0*",20,0.,20.,20,-2.,2.);
+  fHistYPtDsbyDecChannel[0] = new TH2F("hyptdsphi","Ds - vis Phi",20,0.,20.,20,-2.,2.);
+  fHistYPtDsbyDecChannel[1] = new TH2F("hyptdsk0st","Ds - via k0*",20,0.,20.,20,-2.,2.);
+
+  for(Int_t ih=0; ih<2; ih++){
+
+    fHistYPtD0byDecChannel[ih]->Sumw2();
+    fHistYPtD0byDecChannel[ih]->SetMinimum(0);
+    fOutput->Add(fHistYPtD0byDecChannel[ih]);
+    fHistYPtDplusbyDecChannel[ih]->Sumw2();
+    fHistYPtDplusbyDecChannel[ih]->SetMinimum(0);
+    fOutput->Add(fHistYPtDplusbyDecChannel[ih]);
+    fHistYPtDsbyDecChannel[ih]->Sumw2();
+    fHistYPtDsbyDecChannel[ih]->SetMinimum(0);
+    fOutput->Add(fHistYPtDsbyDecChannel[ih]);
+  }
+    
+
+  PostData(1,fOutput);
+
+}
+//______________________________________________________________________________
+void AliAnalysisTaskCheckHFMCProd::UserExec(Option_t *)
+{
+  //
+
+  AliESDEvent *esd = (AliESDEvent*) (InputEvent());
+
+
+  if(!esd) {
+    printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
+    return;
+  } 
+
+  fHistoNEvents->Fill(0);
+
+  Int_t nTracks=esd->GetNumberOfTracks();
+  fHistoTracks->Fill(nTracks);
+  Int_t nSelTracks=0;
+  for(Int_t it=0; it<nTracks; it++){
+    AliESDtrack* tr=esd->GetTrack(it);
+    UInt_t status=tr->GetStatus();
+    if(!(status&AliESDtrack::kITSrefit)) continue;
+    if(!(status&AliESDtrack::kTPCin)) continue;
+    nSelTracks++;
+  }
+  fHistoSelTracks->Fill(nSelTracks);
+
+  const AliMultiplicity* mult=esd->GetMultiplicity();
+  Int_t nTracklets=mult->GetNumberOfTracklets();
+  fHistoTracklets->Fill(nTracklets);
+
+  const AliESDVertex *spdv=esd->GetVertex();
+  if(spdv && spdv->IsFromVertexer3D()){
+    fHistoSPD3DVtxX->Fill(spdv->GetXv());
+    fHistoSPD3DVtxY->Fill(spdv->GetYv());
+    fHistoSPD3DVtxZ->Fill(spdv->GetZv());
+  }
+  if(spdv && spdv->IsFromVertexerZ()){
+    fHistoSPDZVtxX->Fill(spdv->GetXv());
+    fHistoSPDZVtxY->Fill(spdv->GetYv());
+    fHistoSPDZVtxZ->Fill(spdv->GetZv());
+  }
+  const AliESDVertex *trkv=esd->GetPrimaryVertex();
+  if(trkv && trkv->GetNContributors()>1){
+    fHistoTRKVtxX->Fill(trkv->GetXv());
+    fHistoTRKVtxY->Fill(trkv->GetYv());
+    fHistoTRKVtxZ->Fill(trkv->GetZv());
+  }
+
+  AliStack* stack=0;
+
+  if(fReadMC){
+    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
+    if (!eventHandler) {
+      Printf("ERROR: Could not retrieve MC event handler");
+      return;
+    }
+    AliMCEvent* mcEvent = eventHandler->MCEvent();
+    if (!mcEvent) {
+      Printf("ERROR: Could not retrieve MC event");
+      return;
+    }
+    stack = mcEvent->Stack();
+    if (!stack) {
+      Printf("ERROR: stack not available");
+      return;
+    }
+  
+
+    Int_t nParticles=stack->GetNtrack();
+    Double_t dNchdy = 0.;
+    Int_t nb = 0, nc=0;
+    Int_t nCharmed=0.;
+    for (Int_t i=0;i<nParticles;i++){
+      TParticle* part = (TParticle*)stack->Particle(i);
+      Int_t absPdg=TMath::Abs(part->GetPdgCode());
+      if(absPdg==4) nc++;
+      if(absPdg==5) nb++;
+      if(stack->IsPhysicalPrimary(i)){
+       Double_t eta=part->Eta();
+       if(TMath::Abs(eta)<0.5) dNchdy+=0.6666;   // 2/3 for the ratio charged/all    
+      }
+      
+      Int_t iPart=-1;
+      Int_t iType=0;
+      if(absPdg==421){  
+       iType=CheckD0Decay(i,stack);
+       if(iType>=0) iPart=0;   
+      }
+      else if(absPdg==411){
+       iType=CheckDplusDecay(i,stack);
+       if(iType>=0) iPart=1;
+      }
+      else if(absPdg==413){
+       iType=CheckDstarDecay(i,stack);
+       if(iType>=0) iPart=2;
+      }
+      else if(absPdg==431){
+       iType=CheckDsDecay(i,stack);
+       if(iType==0 || iType==1) iPart=3;
+      }
+      else if(absPdg==4122){
+       iType=CheckLcDecay(i,stack);
+       if(iType>=0) iPart=4;
+      }
+      if(iPart<0) continue;
+      if(iType<0) continue;
+      nCharmed++;
+      Float_t rapid=-999.;
+      if (part->Energy() != TMath::Abs(part->Pz())){
+       rapid=0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
+      }
+      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);
+      }
+      
+      TParticle* runningpart=part;
+      Int_t iFromB=-1;
+      while(1){
+       Int_t labmoth=runningpart->GetFirstMother();
+       if(labmoth==-1) break;
+       TParticle *mot=(TParticle*)stack->Particle(labmoth);
+       Int_t pdgmoth=TMath::Abs(mot->GetPdgCode());
+       if(pdgmoth==5){ 
+         iFromB=1;
+         break;
+       }else if(pdgmoth==4){
+         iFromB=0;
+       break;
+       }
+       runningpart=mot;
+      }
+      if(iFromB<0) continue;
+      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);
+      
+    }
+    printf("  ---> %f %d %d %d\n",dNchdy,nCharmed,nc,nb);
+    fHistoNcharmed->Fill(dNchdy,nCharmed);
+    fHistoNbVsNc->Fill(nc,nb);
+  }
+
+  PostData(1,fOutput);
+  
+}
+//______________________________________________________________________________
+void AliAnalysisTaskCheckHFMCProd::Terminate(Option_t */*option*/)
+{
+  // Terminate analysis
+  fOutput = dynamic_cast<TList*> (GetOutputData(1));
+  if (!fOutput) {     
+    printf("ERROR: fOutput not available\n");
+    return;
+  }
+
+  return;
+}
+
+
+
+
+//______________________________________________________________________________
+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>4 || nDau<2) return -1;
+
+  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 && nDau!=2) return -1;
+
+  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 && nDau!=2) return -1;
+
+  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<0) 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);
+  Int_t pdgdp=dp->GetPdgCode();
+  Int_t nDau=dp->GetNDaughters();
+  if(nDau!=3 && nDau!=2) return -1;
+
+  if(nDau==3){
+    Int_t nKaons=0;
+    Int_t nPions=0;
+    Int_t nProtons=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 if(TMath::Abs(pdgdau)==2212){
+       if(pdgdp<0 && pdgdau>0) return -1;
+       if(pdgdp>0 && pdgdau<0) return -1;
+       nProtons++;
+      }else{
+       return -1;
+      }
+    }
+    if(nPions!=1) return -1;
+    if(nKaons!=1) return -1;
+    if(nProtons!=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;
+    Int_t nProtons=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 || TMath::Abs(pdgdau)==3124 || TMath::Abs(pdgdau)==2224
+        || TMath::Abs(pdgdau)==3122 || TMath::Abs(pdgdau)==311){
+       Int_t nDauRes=dau->GetNDaughters();
+       if(nDauRes!=2)  return -1;
+       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++;
+         }
+         if(TMath::Abs(pdgresdau)==2212){
+           if(pdgdp<0 && pdgresdau>0) return -1;
+           if(pdgdp>0 && pdgresdau<0) return -1;
+           nProtons++;
+         }
+       }
+      }else 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 if(TMath::Abs(pdgdau)==2212){
+       if(pdgdp<0 && pdgdau>0) return -1;
+       if(pdgdp>0 && pdgdau<0) return -1;
+       nProtons++;
+      }else{
+       return -1;
+      }
+    }
+    if(nPions!=1) return -1;
+    if(nKaons!=1) return -1;
+    if(nProtons!=1) return -1;
+    for(Int_t iDau=dp->GetFirstDaughter(); iDau<=dp->GetLastDaughter(); iDau++){
+      if(iDau<0) return -1;
+    }
+    return 1;
+  } 
+  return -1;
+}
+
diff --git a/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.h b/PWGHF/vertexingHF/AliAnalysisTaskCheckHFMCProd.h
new file mode 100644 (file)
index 0000000..4e04cbf
--- /dev/null
@@ -0,0 +1,91 @@
+#ifndef ALIANALYSISTASKCHECKHFMCPROD
+#define ALIANALYSISTASKCHECKHFMCPROD
+
+/* Copyright(c) 1998-2012, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */ 
+
+//*************************************************************************
+// Class AliAnalysisTaskCheckHFMCProd
+// AliAnalysisTask to check MC production at ESD+Kine level
+// 
+//
+// Author: F. Prino, prino@to.infn.it
+//          
+//*************************************************************************
+
+class TList;
+class TNtuple;
+class TH1F;
+class TH2F;
+class TTree;
+class TString;
+class AliESDEvent;
+class AliESDfriend;
+class AliStack;
+
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskCheckHFMCProd : public AliAnalysisTaskSE {
+
+ public:
+  
+  AliAnalysisTaskCheckHFMCProd();
+  virtual ~AliAnalysisTaskCheckHFMCProd();
+
+  virtual void   UserExec(Option_t *option);
+  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 SetReadMC(Bool_t opt) {fReadMC=opt;}
+  void SetPbPb() {fPbPb=kTRUE;}
+  void Setpp() {fPbPb=kFALSE;}
+
+ private:
+
+  AliAnalysisTaskCheckHFMCProd(const AliAnalysisTaskCheckHFMCProd &source);
+  AliAnalysisTaskCheckHFMCProd& operator=(const AliAnalysisTaskCheckHFMCProd &source);
+  
+  TList* fOutput;          //! list of output histos
+  TH1F* fHistoNEvents;     //! histo with N of events  
+  
+  TH1F* fHistoTracks;      //! histo with number of ESD tracks
+  TH1F* fHistoSelTracks;   //! histo with number of TPC+ITS refit ESD tracks
+  TH1F* fHistoTracklets;   //! histo with number of SPD tracklets
+  
+  TH1F* fHistoSPD3DVtxX;     //! histo with distr. of x coord. of SPD3D vertex
+  TH1F* fHistoSPD3DVtxY;     //! histo with distr. of y coord. of SPD3D vertex
+  TH1F* fHistoSPD3DVtxZ;     //! histo with distr. of z coord. of SPD3D vertex
+
+  TH1F* fHistoSPDZVtxX;     //! histo with distr. of x coord. of SPDZ vertex
+  TH1F* fHistoSPDZVtxY;     //! histo with distr. of y coord. of SPDZ vertex
+  TH1F* fHistoSPDZVtxZ;     //! histo with distr. of z coord. of SPDZ vertex
+
+  TH1F* fHistoTRKVtxX;     //! histo with distr. of x coord. of TRK vertex
+  TH1F* fHistoTRKVtxY;     //! histo with distr. of y coord. of TRK vertex
+  TH1F* fHistoTRKVtxZ;     //! histo with distr. of z coord. of TRK vertex
+
+  TH2F* fHistoNcharmed;   //! histo of D mesons vs. dN/dy
+  TH2F* fHistoNbVsNc;     //! histo of n. b quarks vs. n c. quarks
+
+  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* fHistYPtDsbyDecChannel[2]; //! histo of y vs. pt for Ds->phipi and Ds->K0*K
+
+  Bool_t fPbPb;
+  Bool_t fReadMC;
+
+  ClassDef(AliAnalysisTaskCheckHFMCProd,1);  
+};
+
+
+#endif
diff --git a/PWGHF/vertexingHF/macros/AddHFMCCheck.C b/PWGHF/vertexingHF/macros/AddHFMCCheck.C
new file mode 100644 (file)
index 0000000..8b10475
--- /dev/null
@@ -0,0 +1,54 @@
+AliAnalysisTaskCheckHFMCProd *AddHFMCCheck(Bool_t isPbPb=kTRUE, Bool_t readMC=kTRUE){
+
+  // Creates, configures and attaches to the train the task for QA of ITS standalone tracks
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AliAnalysisTaskCheckHFMCProd", "No analysis manager to connect to.");
+    return NULL;
+  }   
+  
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler()) {
+    ::Error("AliAnalysisTaskCheckHFMCProd", "This task requires an input event handler");
+    return NULL;
+  }   
+  
+  TString type = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
+  if(type.Contains("AOD")){
+    ::Error("AliAnalysisTaskCheckHFMCProd", "This task requires to run on ESD");
+    return NULL;
+  }
+  
+  //Bool_t isMC=kFALSE;
+  //if (mgr->GetMCtruthEventHandler()) isMC=kTRUE;
+  
+  // Add MC handler (for kinematics)
+  if(readMC){
+    AliMCEventHandler* handler = new AliMCEventHandler;
+    handler->SetReadTR(kFALSE);
+    mgr->SetMCtruthEventHandler(handler);
+  }
+  // Create and configure the task
+  AliAnalysisTaskCheckHFMCProd *task = new AliAnalysisTaskCheckHFMCProd();
+  if(isPbPb) task->SetPbPb();
+  task->SetReadMC(readMC);
+  mgr->AddTask(task);
+  
+  // Create ONLY the output containers for the data produced by the task.
+  // Get and connect other common input/output containers via the manager as below
+  //==============================================================================
+  TString outputFileName = AliAnalysisManager::GetCommonFileName();
+  outputFileName += ":HFMCCheck";
+  
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer("clistHFMCCheck",
+                                                           TList::Class(),
+                                                           AliAnalysisManager::kOutputContainer,
+                                                           outputFileName );
+  
+  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+  mgr->ConnectOutput(task, 1, coutput1);
+  return task;
+}   
diff --git a/PWGHF/vertexingHF/macros/PlotOutputMCCheck.C b/PWGHF/vertexingHF/macros/PlotOutputMCCheck.C
new file mode 100644 (file)
index 0000000..2e7aac3
--- /dev/null
@@ -0,0 +1,146 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TCanvas.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TStyle.h>
+#include <TFile.h>
+#endif
+
+/* $Id$ */ 
+
+// Macro to plot the output of AliAnalysisTaskCheckHFMCProd
+// Author: F. Prino, prino@to.infn.it
+
+
+void PlotOutputMCCheck(){
+  TFile *fil=new TFile("AnalysisResultsMerged.root");
+  TDirectoryFile* df=(TDirectoryFile*)fil->Get("HFMCCheck");
+  TList* l=(TList*)df->Get("clistHFMCCheck");
+  l->ls();
+
+  TH1F* hNEvents=(TH1F*)l->FindObject("hNEvents");
+  Int_t nAnalEv=hNEvents->GetBinContent(1);
+  printf("Number of events= %d\n",nAnalEv);
+
+
+  TCanvas* cv=new TCanvas("cv","Vertex");
+  cv->Divide(3,3);
+  cv->cd(1);
+  TH1F* hSPD3DvX=(TH1F*)l->FindObject("hSPD3DvX");
+  hSPD3DvX->Draw();
+  cv->cd(2);
+  TH1F* hSPD3DvY=(TH1F*)l->FindObject("hSPD3DvY");
+  hSPD3DvY->Draw();
+  cv->cd(3);
+  TH1F* hSPD3DvZ=(TH1F*)l->FindObject("hSPD3DvZ");
+  hSPD3DvZ->Draw();
+  cv->cd(4);
+  TH1F* hSPDZvX=(TH1F*)l->FindObject("hSPDZvX");
+  hSPDZvX->Draw();
+  cv->cd(5);
+  TH1F* hSPDZvY=(TH1F*)l->FindObject("hSPDZvY");
+  hSPDZvY->Draw();
+  cv->cd(6);
+  TH1F* hSPDZvZ=(TH1F*)l->FindObject("hSPDZvZ");
+  hSPDZvZ->Draw();
+  cv->cd(7);
+  TH1F* hTRKvX=(TH1F*)l->FindObject("hTRKvX");
+  hTRKvX->Draw();
+  cv->cd(8);
+  TH1F* hTRKvY=(TH1F*)l->FindObject("hTRKvY");
+  hTRKvY->Draw();
+  cv->cd(9);
+  TH1F* hTRKvZ=(TH1F*)l->FindObject("hTRKvZ");
+  hTRKvZ->Draw();
+
+  TCanvas* c1=new TCanvas("c1","Multipl");
+  c1->Divide(3,1);
+  c1->cd(1);
+  TH1F* hTracklets=(TH1F*)l->FindObject("hTracklets");
+  hTracklets->Draw();
+  c1->cd(2);
+  TH1F* hTracks=(TH1F*)l->FindObject("hTracks");
+  hTracks->Draw();
+  c1->cd(3);
+  TH1F* hSelTracks=(TH1F*)l->FindObject("hSelTracks");
+  hSelTracks->Draw();
+
+  TH1F* hncharmed=(TH1F*)l->FindObject("hncharmed");
+  TCanvas* cn=new TCanvas("cn","ncharm");
+  hncharmed->Draw("box");
+  hncharmed->GetXaxis()->SetTitle("dNch/dy");
+  hncharmed->GetYaxis()->SetTitle("N Charm hadrons in golden channels");
+  cn->Update();
+
+  TH1F* hnbvsnc=(TH1F*)l->FindObject("hnbvsnc");
+  TCanvas* cnhf=new TCanvas("cnhf","nb/c");
+  hnbvsnc->Draw("box");
+  hnbvsnc->GetXaxis()->SetTitle("Nc");
+  hnbvsnc->GetYaxis()->SetTitle("Nb");
+  cnhf->Update();
+
+  TH2F* hyptd0prompt=(TH2F*)l->FindObject("hyptd0prompt");
+  TH2F* hyptd0feeddown=(TH2F*)l->FindObject("hyptd0feeddown");
+  TH2F* hyptD02=(TH2F*)l->FindObject("hyptD02");
+  TH2F* hyptD04=(TH2F*)l->FindObject("hyptD04");
+
+  TCanvas* cd0=new TCanvas("cd0","D0");
+  cd0->Divide(2,2);
+  cd0->cd(1);
+  hyptd0prompt->Draw("colz");
+  cd0->cd(2);
+  hyptd0feeddown->Draw("colz");
+  cd0->cd(3);
+  hyptD02->Draw("colz");
+  cd0->cd(4);
+  hyptD04->Draw("colz");
+
+  TH2F* hyptdplusprompt=(TH2F*)l->FindObject("hyptdplusprompt");
+  TH2F* hyptdplusfeeddown=(TH2F*)l->FindObject("hyptdplusfeeddown");
+  TH2F* hyptDplusnonreson=(TH2F*)l->FindObject("hyptDplusnonreson");
+  TH2F* hyptDplusreson=(TH2F*)l->FindObject("hyptDplusreson");
+
+  TCanvas* cdplus=new TCanvas("cdplus","Dplus");
+  cdplus->Divide(2,2);
+  cdplus->cd(1);
+  hyptdplusprompt->Draw("colz");
+  cdplus->cd(2);
+  hyptdplusfeeddown->Draw("colz");
+  cdplus->cd(3);
+  hyptDplusnonreson->Draw("colz");
+  cdplus->cd(4);
+  hyptDplusreson->Draw("colz");
+
+  TH2F* hyptdsprompt=(TH2F*)l->FindObject("hyptdsprompt");
+  TH2F* hyptdsfeeddown=(TH2F*)l->FindObject("hyptdsfeedown");
+  TH2F* hyptdsphi=(TH2F*)l->FindObject("hyptdsphi");
+  TH2F* hyptdsK0st=(TH2F*)l->FindObject("hyptdsk0st");
+  
+  TCanvas* cds=new TCanvas("cds","Ds");
+  cds->Divide(2,2);
+  cds->cd(1);
+  hyptdsprompt->Draw("colz");
+  cds->cd(2);
+  hyptdsfeeddown->Draw("colz");
+  cds->cd(3);
+  hyptdsphi->Draw("colz");
+  cds->cd(4);
+  hyptdsK0st->Draw("colz");
+
+  TH2F* hyptdstarprompt=(TH2F*)l->FindObject("hyptdstarprompt");
+  TH2F* hyptdstarfeedown=(TH2F*)l->FindObject("hyptdstarfeedown");
+  TH2F* hyptlcprompt=(TH2F*)l->FindObject("hyptlcprompt");
+  TH2F* hyptlcfeedown=(TH2F*)l->FindObject("hyptlcfeedown");
+
+  TCanvas* cdstlc=new TCanvas("cdstls","Dstar LambdaC");
+  cdstlc->Divide(2,2);
+  cdstlc->cd(1);
+  hyptdstarprompt->Draw("colz");
+  cdstlc->cd(2);
+  hyptdstarfeedown->Draw("colz");
+  cdstlc->cd(3);
+  hyptlcprompt->Draw("colz");
+  cdstlc->cd(4);
+  hyptlcfeedown->Draw("colz");
+
+}