New task for prompt charm fraction analysis (Andrea R)
authordainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 19 Jun 2009 09:43:16 +0000 (09:43 +0000)
committerdainese <dainese@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 19 Jun 2009 09:43:16 +0000 (09:43 +0000)
PWG3/CMake_libPWG3vertexingHF.txt
PWG3/PWG3vertexingHFLinkDef.h
PWG3/libPWG3vertexingHF.pkg
PWG3/vertexingHF/AddTaskCharmFraction.C [new file with mode: 0644]
PWG3/vertexingHF/AliAnalysisTaskCharmFraction.cxx [new file with mode: 0644]
PWG3/vertexingHF/AliAnalysisTaskCharmFraction.h [new file with mode: 0644]

index 3f76467953e9173764e46d019f013c22160ac41c..53de11d74c3ef6ebc9ac87311a612efa6aa8ba65 100644 (file)
@@ -10,6 +10,7 @@ set(SRCS
        vertexingHF/AliAnalysisTaskSESelectHF.cxx
        vertexingHF/AliAnalysisTaskSECompareHF.cxx
        vertexingHF/AliAnalysisTaskSEDplus.cxx
+       vertexingHF/AliAnalysisTaskCharmFraction.cxx
        vertexingHF/AliAnalysisTaskSEBkgLikeSignJPSI.cxx
        vertexingHF/AliAnalysisTaskSEBkgLikeSignD0.cxx
        vertexingHF/AliCFHeavyFlavourTask.cxx
index aefa12d40451ec97986ec494061372f776e70054..b8e86ccac41ad2106c17c99d8dc379d203a29f52 100644 (file)
@@ -15,6 +15,7 @@
 #pragma link C++ class AliAnalysisTaskSESelectHF+;
 #pragma link C++ class AliAnalysisTaskSECompareHF+;
 #pragma link C++ class AliAnalysisTaskSEDplus+;
+#pragma link C++ class AliAnalysisTaskCharmFraction+;
 #pragma link C++ class AliCFHeavyFlavourTask+;
 #pragma link C++ class AliCFHeavyFlavourTaskMultiVar+;
 #pragma link C++ class AliCFHeavyFlavourTaskMultiVarMultiStep+;
index 07fc50c93e19004be324bbb39b2945feb1858a04..14a37f51a2cded6c73a4fd1cc5af8a89acd1d16e 100644 (file)
@@ -9,6 +9,7 @@ SRCS:=   vertexingHF/AliAODRecoDecayHF.cxx \
        vertexingHF/AliAnalysisTaskSESelectHF.cxx \
        vertexingHF/AliAnalysisTaskSECompareHF.cxx \
        vertexingHF/AliAnalysisTaskSEDplus.cxx \
+       vertexingHF/AliAnalysisTaskCharmFraction.cxx \
        vertexingHF/AliCFHeavyFlavourTask.cxx \
        vertexingHF/AliCFHeavyFlavourTaskMultiVar.cxx \
        vertexingHF/AliCFHeavyFlavourTaskMultiVarMultiStep.cxx \
diff --git a/PWG3/vertexingHF/AddTaskCharmFraction.C b/PWG3/vertexingHF/AddTaskCharmFraction.C
new file mode 100644 (file)
index 0000000..9d77c86
--- /dev/null
@@ -0,0 +1,326 @@
+AliAnalysisTaskCharmFraction* AddTaskCharmFraction(
+         const char* fileout="d0D0.root",
+        Bool_t sideband=kFALSE,
+        Bool_t setD0usecuts=kTRUE,
+        Bool_t setcheckMC=kTRUE,
+        Bool_t setcheckMC_prompt=kTRUE,
+        Bool_t setcheckMC_fromB=kFALSE,
+        Bool_t setcheckMC_D0=kTRUE,
+        Bool_t setcheckMC_2prongs=kTRUE,
+        Bool_t setSkipD0star=kTRUE,
+        Bool_t setStudyPureBack=kFALSE)
+{  
+  //
+  // Configuration macro for the task to analyze the fraction of prompt charm
+  // using the D0 impact parameter
+  // andrea.rossi@ts.infn.it
+  //
+  //==========================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddTaskCharmFraction", "No analysis manager to connect to.");
+    return NULL;
+  }   
+
+  TString str=fileout,containername;
+  str.ReplaceAll(".root","");
+  str.Prepend("_");
+
+  AliAnalysisTaskCharmFraction *hfTask;
+  if(!sideband) {
+    hfTask = new AliAnalysisTaskCharmFraction("CharmFraction",10);
+  } else {
+    hfTask= new AliAnalysisTaskCharmFraction("CharmFractionSideB",10);
+    hfTask->SetSideBands(-2.);
+  }
+
+  hfTask->SetUseCuts(setD0usecuts);
+  hfTask->SetCheckMC(setcheckMC);
+  hfTask->SetCheckMC_D0(setcheckMC_D0);
+  hfTask->SetCheckMC_2prongs(setcheckMC_2prongs);
+  hfTask->SetCheckMC_prompt(setcheckMC_prompt);
+  hfTask->SetCheckMC_fromB(setcheckMC_fromB);
+  hfTask->SetCheckMC_fromDstar(setSkipD0star);
+  hfTask->SetStudyPureBackground(setStudyPureBack);
+  //  hfTask->SetSideBands(0);
+  //  hfTask->SetDebugLevel(2);
+  mgr->AddTask(hfTask);
+  //Now the same for sidebands
+  /*AliAnalysisTaskCharmFraction *hfTaskSideB 
+
+  
+  hfTaskSideB->SetUseCuts(fD0usecuts);
+  hfTaskSideB->SetCheckMC(fcheckMC);
+  hfTaskSideB->SetCheckMC_D0(fcheckMC_D0);
+  hfTaskSideB->SetCheckMC_2prongs(fcheckMC_2prongs);
+  hfTaskSideB->SetCheckMC_prompt(fcheckMC_prompt);
+  hfTaskSideB->SetCheckMC_fromB(fcheckMC_fromB);
+  hfTaskSideB->SetCheckMC_fromDstar(fSkipD0star);
+  hfTaskSideB->SetStudyPureBackground(fStudyPureBack);
+
+  //  hfTaskSideB->SetDebugLevel(2);
+  mgr->AddTask(hfTaskSideB);
+  */
+
+  // Create containers for input/output
+  AliAnalysisDataContainer *cinput =   mgr->GetCommonInputContainer();//mgr->CreateContainer("cinput",TChain::Class(),AliAnalysisManager::kInputContainer);
+  mgr->ConnectInput(hfTask,0,cinput);
+  //  mgr->ConnectInput(hfTaskSideB,0,cinput);
+
+  //Now container for general properties histograms
+  containername="coutputCptd0d0";
+  containername.Append(str.Data());
+  AliAnalysisDataContainer *coutputCptd0d0 = mgr->CreateContainer(containername.Data(),TH2::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          fileout);
+  mgr->ConnectOutput(hfTask,0,coutputCptd0d0);
+
+  containername="coutputSecVtxXY";
+  containername.Append(str.Data());
+  AliAnalysisDataContainer *coutputSecVtxXY = mgr->CreateContainer(containername.Data(),TH2::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          fileout);
+  mgr->ConnectOutput(hfTask,1,coutputSecVtxXY);
+
+
+  containername="coutputd0d0";
+  containername.Append(str.Data());
+  AliAnalysisDataContainer *coutputd0d0 = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          fileout);
+  mgr->ConnectOutput(hfTask,2,coutputd0d0);
+
+  containername="coutputCpt";
+  containername.Append(str.Data());
+  AliAnalysisDataContainer *coutputCpt = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          fileout);
+  mgr->ConnectOutput(hfTask,3,coutputCpt);
+
+  containername="coutputSecVtxZ";
+  containername.Append(str.Data());
+  AliAnalysisDataContainer *coutputSecVtxZ = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          fileout);
+  mgr->ConnectOutput(hfTask,4,coutputSecVtxZ);
+
+  containername="coutputSecVtxX";
+  containername.Append(str.Data());
+  AliAnalysisDataContainer *coutputSecVtxX = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          fileout);
+  mgr->ConnectOutput(hfTask,5,coutputSecVtxX);
+
+  containername="coutputSecVtxY";
+  containername.Append(str.Data());
+  AliAnalysisDataContainer *coutputSecVtxY = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          fileout);
+  mgr->ConnectOutput(hfTask,6,coutputSecVtxY);
+
+  containername="coutputSecVtxPhi";
+  containername.Append(str.Data());
+ AliAnalysisDataContainer *coutputSecVtxPhi = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          fileout);
+  mgr->ConnectOutput(hfTask,7,coutputSecVtxPhi);
+
+
+  //Now container for d0D0  
+  AliAnalysisDataContainer **coutput=new AliAnalysisDataContainer*[10];
+  containername="coutputAll";
+  containername.Append(str.Data());
+  AliAnalysisDataContainer *coutputAll = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          fileout);
+
+
+  TString name="coutput";
+  for(Int_t j=0;j<10;j++){
+    containername=name;
+    containername+=j;
+    containername.Append(str.Data());
+    coutput[j] = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                     AliAnalysisManager::kOutputContainer, 
+                                     fileout);
+    
+    mgr->ConnectOutput(hfTask,j+8,coutput[j]);
+  }
+  mgr->ConnectOutput(hfTask,18,coutputAll);
+  //Now container for MC d0D0  
+  AliAnalysisDataContainer **coutputMC=new AliAnalysisDataContainer*[10];
+  containername="coutputAllMC";
+  containername.Append(str.Data());
+  AliAnalysisDataContainer *coutputAllMC = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          fileout);
+
+
+  name="coutputMC";
+  for(Int_t j=0;j<10;j++){
+    containername=name;
+    containername+=j;
+    containername.Append(str.Data());
+    coutputMC[j] = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                     AliAnalysisManager::kOutputContainer, 
+                                     fileout);
+    
+    mgr->ConnectOutput(hfTask,j+19,coutputMC[j]);
+  }
+  mgr->ConnectOutput(hfTask,29,coutputAllMC);
+
+  //Now container for histo with d0 with respect to True Vtx
+  AliAnalysisDataContainer **coutputd0VtxTrue=new AliAnalysisDataContainer*[10];
+  containername="coutputd0VtxTrueAll";
+  containername.Append(str.Data());
+  AliAnalysisDataContainer *coutputd0VtxTrueAll = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          fileout);
+
+
+  name="coutputd0VtxTrue";
+  for(Int_t j=0;j<10;j++){
+    containername=name;
+    containername+=j;
+    containername.Append(str.Data());
+    coutputd0VtxTrue[j] = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                     AliAnalysisManager::kOutputContainer, 
+                                     fileout);
+    
+    mgr->ConnectOutput(hfTask,j+30,coutputd0VtxTrue[j]);
+  }
+  mgr->ConnectOutput(hfTask,40,coutputd0VtxTrueAll);
+  //INV MASS
+  containername="coutputD0InvMassAll";
+  containername.Append(str.Data());
+  AliAnalysisDataContainer *coutputD0InvMassAll = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                                                      AliAnalysisManager::kOutputContainer, 
+                                                                      fileout);
+  mgr->ConnectOutput(hfTask,41,coutputD0InvMassAll);
+  containername="coutputD0MCInvMassAll";
+  containername.Append(str.Data());
+  AliAnalysisDataContainer *coutputD0MCInvMassAll = mgr->CreateContainer(containername.Data(),TH1::Class(),
+                                                                        AliAnalysisManager::kOutputContainer, 
+                                                                        fileout);
+  mgr->ConnectOutput(hfTask,42,coutputD0MCInvMassAll);
+  
+ ////////
+ //NOW THE SAME FOR SIDE BANDS
+ /*
+ //Now container for general properties histograms
+
+ AliAnalysisDataContainer *coutputSBCptd0d0 = mgr->CreateContainer("coutputSBCptd0d0",TH2::Class(),
+                                                                  AliAnalysisManager::kOutputContainer, 
+                                                          "d0D0SideB.root");
+ mgr->ConnectOutput(hfTaskSideB,0,coutputSBCptd0d0);
+
+  AliAnalysisDataContainer *coutputSBSecVtxXY = mgr->CreateContainer("coutputSBSecVtxXY",TH2::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          "d0D0SideB.root");
+  mgr->ConnectOutput(hfTaskSideB,1,coutputSBSecVtxXY);
+
+
+  AliAnalysisDataContainer *coutputSBd0d0 = mgr->CreateContainer("coutputSBd0d0",TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          "d0D0SideB.root");
+  mgr->ConnectOutput(hfTaskSideB,2,coutputSBd0d0);
+
+  AliAnalysisDataContainer *coutputSBCpt = mgr->CreateContainer("coutputSBCpt",TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          "d0D0SideB.root");
+  mgr->ConnectOutput(hfTaskSideB,3,coutputSBCpt);
+
+  AliAnalysisDataContainer *coutputSBSecVtxZ = mgr->CreateContainer("coutputSBSecVtxZ",TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          "d0D0SideB.root");
+  mgr->ConnectOutput(hfTaskSideB,4,coutputSBSecVtxZ);
+
+  AliAnalysisDataContainer *coutputSBSecVtxX = mgr->CreateContainer("coutputSBSecVtxX",TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          "d0D0SideB.root");
+  mgr->ConnectOutput(hfTaskSideB,5,coutputSBSecVtxX);
+
+  AliAnalysisDataContainer *coutputSBSecVtxY = mgr->CreateContainer("coutputSBSecVtxY",TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          "d0D0SideB.root");
+  mgr->ConnectOutput(hfTaskSideB,6,coutputSBSecVtxY);
+
+ AliAnalysisDataContainer *coutputSBSecVtxPhi = mgr->CreateContainer("coutputSBSecVtxPhi",TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          "d0D0SideB.root");
+  mgr->ConnectOutput(hfTaskSideB,7,coutputSBSecVtxPhi);
+
+
+  //Now container for d0D0SideB  
+  AliAnalysisDataContainer **coutputSB=new AliAnalysisDataContainer*[10];
+  AliAnalysisDataContainer *coutputSBAll = mgr->CreateContainer("coutputSBAll",TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          "d0D0SideB.root");
+
+
+  TString name="coutputSB",strname;
+  for(Int_t j=0;j<10;j++){
+    strname=name;
+    strname+=j;
+    coutputSB[j] = mgr->CreateContainer(strname.Data(),TH1::Class(),
+                                     AliAnalysisManager::kOutputContainer, 
+                                     "d0D0SideB.root");
+    
+    mgr->ConnectOutput(hfTaskSideB,j+8,coutputSB[j]);
+  }
+  mgr->ConnectOutput(hfTaskSideB,18,coutputSBAll);
+  //Now container for MC d0D0SideB  
+  AliAnalysisDataContainer **coutputSBMC=new AliAnalysisDataContainer*[10];
+  AliAnalysisDataContainer *coutputSBAllMC = mgr->CreateContainer("coutputSBAllMC",TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          "d0D0SideB.root");
+
+
+  name="coutputSBMC";
+  for(Int_t j=0;j<10;j++){
+    strname=name;
+    strname+=j;
+    coutputSBMC[j] = mgr->CreateContainer(strname.Data(),TH1::Class(),
+                                     AliAnalysisManager::kOutputContainer, 
+                                     "d0D0SideB.root");
+    
+    mgr->ConnectOutput(hfTaskSideB,j+19,coutputSBMC[j]);
+  }
+  mgr->ConnectOutput(hfTaskSideB,29,coutputSBAllMC);
+
+  //Now container for histo with d0 with respect to True Vtx
+  AliAnalysisDataContainer **coutputSBd0VtxTrue=new AliAnalysisDataContainer*[10];
+  AliAnalysisDataContainer *coutputSBd0VtxTrueAll = mgr->CreateContainer("coutputSBd0VtxTrueAll",TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          "d0D0SideB.root");
+
+
+  name="coutputSBd0VtxTrue";
+  for(Int_t j=0;j<10;j++){
+    strname=name;
+    strname+=j;
+    coutputSBd0VtxTrue[j] = mgr->CreateContainer(strname.Data(),TH1::Class(),
+                                     AliAnalysisManager::kOutputContainer, 
+                                     "d0D0SideB.root");
+    
+    mgr->ConnectOutput(hfTaskSideB,j+30,coutputSBd0VtxTrue[j]);
+  }
+  mgr->ConnectOutput(hfTaskSideB,40,coutputSBd0VtxTrueAll);
+
+//INV MASS
+ AliAnalysisDataContainer *coutputSBD0InvMassAll = mgr->CreateContainer("coutputSBD0InvMassAll",TH1::Class(),
+                                                          AliAnalysisManager::kOutputContainer, 
+                                                          "d0D0SideB.root");
+ mgr->ConnectOutput(hfTaskSideB,41,coutputSBD0InvMassAll);
+
+ AliAnalysisDataContainer *coutputSBD0MCInvMassAll = mgr->CreateContainer("coutputSBD0MCInvMassAll",TH1::Class(),
+                                                                     AliAnalysisManager::kOutputContainer, 
+                                                          "d0D0SideB.root");
+ mgr->ConnectOutput(hfTaskSideB,42,coutputSBD0MCInvMassAll);
+
+ */
+
+  return hfTask;
+}
diff --git a/PWG3/vertexingHF/AliAnalysisTaskCharmFraction.cxx b/PWG3/vertexingHF/AliAnalysisTaskCharmFraction.cxx
new file mode 100644 (file)
index 0000000..cbc100c
--- /dev/null
@@ -0,0 +1,691 @@
+/**************************************************************************
+ * Copyright(c) 1998-2009, 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.                  *
+ **************************************************************************/
+
+/////////////////////////////////////////////////////////////
+//
+// Class AliAnalysisTaskCharmFraction
+// AliAnalysisTask for the extraction of the fraction of prompt charm
+// using the charm hadron impact parameter to the primary vertex
+//
+// Author: Andrea Rossi, andrea.rossi@ts.infn.it
+/////////////////////////////////////////////////////////////
+
+#include <TChain.h>
+#include <TTree.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TCanvas.h>
+#include <TDatabasePDG.h>
+#include <TMath.h>
+
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliAnalysisVertexingHF.h"
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODRecoDecayHF.h"
+#include "AliAODRecoDecay.h"
+#include "AliAODTrack.h"
+#include "AliAODVertex.h"
+#include "AliAODMCParticle.h"
+#include "AliAODMCHeader.h"
+#include "AliAnalysisTaskCharmFraction.h"
+
+ClassImp(AliAnalysisTaskCharmFraction)
+  //________________________________________________________________________
+  AliAnalysisTaskCharmFraction::AliAnalysisTaskCharmFraction(const char *name) 
+    : AliAnalysisTask(name, ""), fAOD(0), fVHF(0),
+      fArrayD0toKpi(0),
+      fArrayMC(0),
+      fAODmcHeader(0),
+      fhMass(0),
+      fhMassTrue(0),
+      fhCPtaVSd0d0(0),
+      fhd0xd0(0),
+      fhCPta(0),
+      fhSecVtxZ(0),
+      fhSecVtxX(0),
+      fhSecVtxY(0),
+      fhSecVtxXY(0),
+      fhSecVtxPhi(0),
+      fhd0D0(0),
+      fhd0D0pt(0x0),
+      fhd0D0VtxTrue(0),
+      fhd0D0VtxTruept(0x0),
+      fhMCd0D0(0),
+      fhMCd0D0pt(0x0),
+      fnbins(),
+      fD0usecuts(kTRUE),
+      fcheckMC(kFALSE),
+      fcheckMCD0(kFALSE),
+      fcheckMC2prongs(kFALSE),
+      fcheckMCprompt(kFALSE),
+      fcheckMCfromB(kFALSE),
+      fSkipD0star(kFALSE),
+      fStudyd0fromBTrue(kTRUE),
+      fStudyPureBackground(kFALSE),
+      fSideBands(0)
+{
+  // Constructor
+  //Side bands selection (see cxx): 999 -> no inv mass selection at all,
+  //                                >0 ->  both D0 and D0bar hypothesis are required to be kFALSE, 
+  //                                <0 -> only 1 inv mass hypothesis is required to be KFALSE
+  SetNPtBins();
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
+  DefineInput(0, TChain::Class());
+  //  DefineInput(1, TChain::Class());
+  // Output slot #0 writes into a TH1 container
+  DefineOutput(0, TH2F::Class());// hCPtaVsd0d0
+  DefineOutput(1, TH2F::Class());// hVtxXY
+  for(Int_t j=2;j<3*fnbins+11;j++){ // hCPta, hd0xd0, hVtZ, hVtX hVtY hVtPhi + nptbin hd0D0 +hd0D0 integreated +nptbin hd0VtxTrueD0 +hd0VtxTrueD0 integrated +nptbin hMCd0D0 +hMCd0D0 integrated
+    DefineOutput(j, TH1F::Class());
+  }
+  //  DefineOutput(fnbins+8, TH1F::Class());// hd0D0 pt integrated
+}
+
+//________________________________________________________________________
+AliAnalysisTaskCharmFraction::AliAnalysisTaskCharmFraction(const char *name,Int_t nptbins) 
+  : AliAnalysisTask(name, ""), fAOD(0), fVHF(0),
+    fArrayD0toKpi(0),
+    fArrayMC(0),
+    fAODmcHeader(0),
+    fhMass(0),
+    fhMassTrue(0),
+    fhCPtaVSd0d0(0),
+    fhd0xd0(0),
+    fhCPta(0),
+    fhSecVtxZ(0),
+    fhSecVtxX(0),
+    fhSecVtxY(0),
+    fhSecVtxXY(0),
+    fhSecVtxPhi(0),
+    fhd0D0(0),
+    fhd0D0pt(0x0),
+    fhd0D0VtxTrue(0),
+    fhd0D0VtxTruept(0x0),
+    fhMCd0D0(0),
+    fhMCd0D0pt(0x0),
+    fnbins(),
+    fD0usecuts(kTRUE),
+    fcheckMC(kFALSE),
+    fcheckMCD0(kFALSE),
+    fcheckMC2prongs(kFALSE),
+    fcheckMCprompt(kFALSE),
+    fcheckMCfromB(kFALSE),
+    fSkipD0star(kFALSE),
+    fStudyd0fromBTrue(kTRUE),
+    fStudyPureBackground(kFALSE),
+    fSideBands(0)        
+{
+  // Constructor
+  //Side bands selection (see cxx): 999 -> no inv mass selection at all,
+  //                                >0 ->  both D0 and D0bar hypothesis are required to be kFALSE, 
+  //                                <0 -> only 1 inv mass hypothesis is required to be KFALSE
+  SetNPtBins(nptbins);
+  // Define input and output slots here
+  // Input slot #0 works with a TChain
+  DefineInput(0, TChain::Class());
+  //  DefineInput(1, TChain::Class());
+  // Output slot #0 writes into a TH1 container
+  
+  DefineOutput(0, TH2F::Class());// hCPtaVsd0d0
+  DefineOutput(1, TH2F::Class());// hVtxXY
+  for(Int_t j=2;j<3*fnbins+13;j++){ // hCPta, hd0xd0, hVtZ, hVtX hVtY hVtPhi + nptbin hd0D0 +hd0D0 integreated +nptbin hd0VtxTrueD0 +hd0VtxTrueD0 integrated +nptbin hMCd0D0 +hMCd0D0 integrated
+    DefineOutput(j, TH1F::Class());
+  }
+  //  DefineOutput(fnbins+8, TH1F::Class());// hd0D0 pt integrated 
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskCharmFraction::ConnectInputData(Option_t *) 
+{
+  // Connect ESD or AOD here
+  // Called once
+
+  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
+  //  TTree* treeFriend = dynamic_cast<TTree*> (GetInputData(1));
+  if (!tree) {
+    Printf("ERROR: Could not read chain from input slot 0");
+  } else {
+    //  tree->AddFriend(treeFriend);
+    
+    // Disable all branches and enable only the needed ones
+    // The next two lines are different when data produced as AliESDEvent is read
+    
+    AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+    fAOD=new AliAODEvent();
+    fAOD->ReadFromTree(tree);
+    
+    fArrayD0toKpi = (TClonesArray*)fAOD->GetList()->FindObject("D0toKpi"); 
+    fArrayMC =      (TClonesArray*)fAOD->GetList()->FindObject("mcparticles"); 
+    
+    fAODmcHeader=(AliAODMCHeader*)fAOD->GetList()->FindObject("mcHeader");
+    
+    if (!aodH) {
+      Printf("ERROR: Could not get AODInputHandler");
+    } else
+      fAOD = aodH->GetEvent();
+    
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskCharmFraction::CreateOutputObjects()
+{
+  // Create histograms
+  // Called once
+  fhd0D0 = new TH1F("hd0D0","D^{0} impact par. plot  (All momenta)",1000,-1000.,1000.);
+  fhd0D0->SetXTitle("Impact parameter [#mum]");
+  fhd0D0->SetYTitle("Entries");
+
+  fhd0D0VtxTrue = new TH1F("hd0D0VtxTrue","D^{0} impact par. w.r.t. True Vtx (All momenta)",1000,-1000.,1000.);
+  fhd0D0VtxTrue->SetXTitle("Impact parameter [#mum]");
+  fhd0D0VtxTrue->SetYTitle("Entries");
+
+  fhMCd0D0 = new TH1F("hMCd0D0","D^{0} impact par. plot  (All momenta)",1000,-1000.,1000.);
+  fhMCd0D0->SetXTitle("MC Impact parameter [#mum]");
+  fhMCd0D0->SetYTitle("Entries");
+
+  fhCPtaVSd0d0=new TH2F("hCPtaVSd0d0","hCPtaVSd0d0",1000,-100000.,100000.,100,0.,1.);
+  fhSecVtxZ=new TH1F("hSecVtxZ","hSecVtxZ",1000,-8.,8.);
+  fhSecVtxX=new TH1F("hSecVtxX","hSecVtxX",1000,-3000.,3000.);
+  fhSecVtxY=new TH1F("hSecVtxY","hSecVtxY",1000,-3000.,3000.);
+  fhSecVtxXY=new TH2F("hSecVtxXY","hSecVtxXY",1000,-3000.,3000.,1000,-3000.,3000.);
+  fhSecVtxPhi=new TH1F("hSecVtxPhi","hSecVtxPhi",180,-180.1,180.1);
+  fhCPta=new TH1F("hCPta","hCPta",100,0.,1.);
+  fhd0xd0=new TH1F("hd0xd0","hd0xd0",1000,-100000.,100000.);
+  fhMassTrue=new TH1F("hMassTrue","D^{0} MC inv. Mass (All momenta)",600,1.600,2.200);
+  fhMass=new TH1F("hMass","D^{0} inv. Mass (All momenta)",600,1.600,2.200);
+  Double_t standardptbins[14]={0.,0.5,1.,2.,3.,5.,7.,10.,15.,20.,30.,40.,50.,100.};
+  
+  fhd0D0pt=new TH1F*[fnbins];
+  fhMCd0D0pt=new TH1F*[fnbins];
+  fhd0D0VtxTruept=new TH1F*[fnbins];
+  TString namehist="hd0D0pt";
+  TString titlehist="D^{0} impact par. plot ";
+  TString strnamept,strtitlept;
+  for(Int_t i=0;i<fnbins;i++){
+    strnamept=namehist;
+    strnamept+=i;// strnamept+=standardptbins[i+1];
+    //    strnamept.Append("_");
+    //strnamept+=standardptbins[i+1];
+
+    strtitlept=namehist;
+    strtitlept+=standardptbins[i];
+    strtitlept.Append("<= pt <");
+    strtitlept+=standardptbins[i+1];
+    strtitlept.Append(" [GeV/c]");
+    
+    fhd0D0pt[i] = new TH1F(strnamept.Data(),strtitlept.Data(),1000,-1000.,1000.);
+    fhd0D0pt[i]->SetXTitle("Impact parameter [#mum] ");
+    fhd0D0pt[i]->SetYTitle("Entries");
+
+    strnamept.ReplaceAll("hd0D0","hMCd0D0");
+    fhMCd0D0pt[i] = new TH1F(strnamept.Data(),strtitlept.Data(),1000,-1000.,1000.);
+    fhMCd0D0pt[i]->SetXTitle("MC Impact parameter [#mum] ");
+    fhMCd0D0pt[i]->SetYTitle("Entries");
+
+    strnamept.ReplaceAll("hMCd0D0","hd0D0VtxTrue");
+    fhd0D0VtxTruept[i] = new TH1F(strnamept.Data(),strtitlept.Data(),1000,-1000.,1000.);
+    fhd0D0VtxTruept[i]->SetXTitle("Impact parameter w.r.t. True Vtx [#mum] ");
+    fhd0D0VtxTruept[i]->SetYTitle("Entries");
+    
+  }
+
+
+  /*  fHistPt = new TH1F("fHistPt", "P_{T} distribution", 15, 0.1, 3.1);
+  fHistPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
+  fHistPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
+  fHistPt->SetMarkerStyle(kFullCircle);*/
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskCharmFraction::Exec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+
+  if (!fAOD) {
+    Printf("ERROR: fAOD not available");
+    return;
+  }
+
+  //Printf("There are %d tracks in this event", fAOD->GetNumberOfTracks());
+  //  Int_t nTotHF=0,nTotDstar=0,nTot3Prong=0;
+  Int_t nTotD0toKpi=0;
+  //  const  Int_t nptbins=10;
+  Double_t standardptbins[14]={0.,0.5,1.,2.,3.,5.,7.,10.,15.,20.,30.,40.,50.,100.};
+  Double_t pD[3],xD[3],pXtrTrue[2],pYtrTrue[2],pZtrTrue[2],xtr1[3],xtr2[3];
+  Double_t cutsD0[9]=
+       // cutsD0[0] = inv. mass half width [GeV]
+       // cutsD0[1] = dca [cm]
+       // cutsD0[2] = cosThetaStar
+       // cutsD0[3] = pTK [GeV/c]
+       // cutsD0[4] = pTPi [GeV/c]
+       // cutsD0[5] = d0K [cm]   upper limit!
+       // cutsD0[6] = d0Pi [cm]  upper limit!
+       // cutsD0[7] = d0d0 [cm^2]
+       // cutsD0[8] = cosThetaPoint
+    {1000.,
+     100000.,
+     1.1,
+     0.,
+     0.,
+     100000.,
+     100000.,
+     100000000.,
+     -1.1};   
+  if(fD0usecuts){
+    cutsD0[0]=0.033;
+    cutsD0[1]=0.03;
+    cutsD0[2]=0.8;
+    cutsD0[3]=0.;
+    cutsD0[4]=0.;
+    cutsD0[5]=100000.;
+    cutsD0[6]=100000.;
+    cutsD0[7]=-0.00005;
+    cutsD0[8]=0.8; 
+  }
+  const Double_t fixedInvMassCut=cutsD0[0];
+  AliAODVertex *vtx1=0;
+  // primary vertex
+  vtx1 = (AliAODVertex*)fAOD->GetPrimaryVertex();
+  Double_t vtxTrue[3];
+  fAODmcHeader->GetVertex(vtxTrue);
+  //   vtx1->Print();
+  AliAODRecoDecayHF *aodDMC=0x0;
+  Bool_t nomum=kFALSE,background=kFALSE;
+  AliAODMCParticle *mum1=0x0,*mum2=0x0;
+  AliAODMCParticle *b1=0x0,*b2=0x0;
+  AliAODMCParticle *sonpart=0x0,*grandmoth1=0x0;
+  
+  // make trkIDtoEntry register (temporary)
+  Int_t trkIDtoEntry[100000];
+  for(Int_t it=0;it<fAOD->GetNumberOfTracks();it++) {
+    AliAODTrack *track = fAOD->GetTrack(it);
+    if(track->GetID()<0) {
+      //printf("Track ID <0, id= %d\n",track->GetID());
+      return;
+    }
+    trkIDtoEntry[track->GetID()]=it;
+  }
+  
+  
+  // loop over D0->Kpi candidates
+  Int_t nD0toKpi = fArrayD0toKpi->GetEntriesFast();
+  nTotD0toKpi += nD0toKpi;
+  //   cout<<"Number of D0->Kpi: "<<nD0toKpi<<endl;
+  
+  for (Int_t iD0toKpi = 0; iD0toKpi < nD0toKpi; iD0toKpi++) {
+    if(aodDMC!=0x0)delete aodDMC;
+    mum1=0x0;
+    mum2=0x0;
+    b1=0x0;
+    b2=0x0;
+    sonpart=0x0;
+    grandmoth1=0x0;
+    
+    AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)fArrayD0toKpi->UncheckedAt(iD0toKpi);
+    Bool_t unsetvtx=kFALSE;
+    if(!d->GetOwnPrimaryVtx()) {
+      d->SetOwnPrimaryVtx(vtx1); // needed to compute all variables
+      unsetvtx=kTRUE;
+    }
+    
+    //    d->SetOwnPrimaryVtx(vtx1); 
+    //unsetvtx=kTRUE;
+    
+    Int_t okD0=0,okD0bar=0; 
+    Double_t invMassD0,invMassD0bar,cutmassvalue;
+    Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
+    nomum=kFALSE;
+    background=kFALSE;
+
+    cutsD0[0]=fixedInvMassCut;
+    if(fSideBands){
+      cutmassvalue=cutsD0[0];
+      if(fSideBands==999.)cutsD0[0]=1000.;
+      else cutsD0[0]=2.*TMath::Abs(fSideBands)*cutsD0[0];
+    }
+    if(d->SelectD0(cutsD0,okD0,okD0bar)) {
+      if(fSideBands){
+       if(fSideBands==999.)cutmassvalue=-1.;
+       d->InvMassD0(invMassD0,invMassD0bar);
+       
+       if(TMath::Abs(invMassD0-mD0PDG)    > TMath::Abs(fSideBands)*cutmassvalue) okD0 = okD0 && kTRUE;
+       else okD0=kFALSE;
+       if(TMath::Abs(invMassD0bar-mD0PDG) > TMath::Abs(fSideBands)*cutmassvalue) okD0bar = okD0bar && kTRUE;    // SIDE BANDS SELECTION
+       else okD0bar=kFALSE;
+
+       if(fSideBands<0.){
+         if(!okD0 && !okD0bar) continue;
+       }
+       else if (!okD0 || !okD0bar)continue;
+      }
+      
+      
+      // get daughter AOD tracks
+      AliAODTrack *trk0 = (AliAODTrack*)d->GetDaughter(0);
+      AliAODTrack *trk1 = (AliAODTrack*)d->GetDaughter(1);
+      
+      //Check if it's a True D0
+      if(fcheckMC){
+       while(!background){
+         
+         if(trk0->GetLabel()==-1||trk1->GetLabel()==-1){
+           //fake tracks!!
+           background=kTRUE;  //FAKE TRACK
+           break;
+         }
+         //      printf("Before entering the MC checks \n");
+         
+         b1=(AliAODMCParticle*)fArrayMC->At(trk0->GetLabel());
+         b2=(AliAODMCParticle*)fArrayMC->At(trk1->GetLabel());
+         
+         
+         if(b1->GetMother()==-1||b2->GetMother()==-1){
+           nomum=kTRUE; //Tracks with no mother 
+           //printf("Inside problem mothering \n");// FAKE DECAY VERTEX
+           background=kTRUE;
+           break;
+           // if(!fStudyPureBackground)continue;
+         }
+       
+         if(b1->GetMother()!=b2->GetMother()){
+           //Check the label of the mother is the same
+           background=kTRUE; // NOT SAME MOTHER
+           break;
+         }
+         //      printf("After first mothering \n");
+         mum1=(AliAODMCParticle*)fArrayMC->At(b1->GetMother());
+         mum2=(AliAODMCParticle*)fArrayMC->At(b2->GetMother());
+         //    if((mum1->GetPdgCode()!=mum2->GetPdgCode()))continue; //Check the mother is the same particle
+         // printf("After second mothering \n");
+         //      printf("Particle codes: tr1: %d, tr2: %d, mum1: %d, mum 2: %d \n",b1->GetPdgCode(),b2->GetPdgCode(),mum1->GetPdgCode(),mum2->GetPdgCode());
+         if(fcheckMCD0)if(TMath::Abs(mum1->GetPdgCode())!=421){
+           //^fStudyPureBackground)continue;//NPRONG
+           background=kTRUE;// NOT A D0
+           break;
+         }
+         if(fcheckMC2prongs)if(mum1->GetDaughter(1)-mum1->GetDaughter(0)+1!=2){
+           //^fStudyPureBackground)continue;// NPRONG
+           background=kTRUE; // NOT A 2 PRONG DECAY
+           break;
+         }
+         
+         if(mum1->GetMother()==-1){
+           background=kTRUE; // A particle coming from nothing
+           break;
+         }
+             
+         // ^fStudyPureBackground)continue;
+         //      printf("After second mothering \n");
+         grandmoth1=(AliAODMCParticle*)fArrayMC->At(mum1->GetMother());
+         if(fSkipD0star)if(TMath::Abs(grandmoth1->GetPdgCode())==413||TMath::Abs(grandmoth1->GetPdgCode())==423){
+           //^fStudyPureBackground)continue;
+           background=kTRUE;// D0 COMING FROM A D0*
+           break;
+         }
+         if(fcheckMCD0){//THIS CHECK IS NEEDE TO AVOID POSSIBLE FAILURE IN THE SECOND WHILE
+           while(TMath::Abs(grandmoth1->GetPdgCode())!=4&&TMath::Abs(grandmoth1->GetPdgCode())!=5){
+             if(grandmoth1->GetMother()==-1){
+               //            printf("Inside grandmothering \n");
+               //printf("mother=-1, pdgcode: %d \n",grandmoth1->GetPdgCode());
+               Int_t son=grandmoth1->GetDaughter(0);
+               sonpart=(AliAODMCParticle*)fArrayMC->At(son);
+               while(TMath::Abs(sonpart->GetPdgCode())!=421){
+                 //printf("mother=-1, pdgcode: %d \n",sonpart->GetPdgCode());
+                 son++;
+                 sonpart=(AliAODMCParticle*)fArrayMC->At(son);
+               }
+               break;
+             }
+             grandmoth1=(AliAODMCParticle*)fArrayMC->At(grandmoth1->GetMother());
+           }
+         }
+         if(fcheckMCprompt)if(TMath::Abs(grandmoth1->GetPdgCode())!=4){
+           //^fStudyPureBackground)continue;
+           background=kTRUE;
+           break;
+         }
+         if(fcheckMCfromB){
+           if(TMath::Abs(grandmoth1->GetPdgCode())!=5){
+             //^fStudyPureBackground)continue;
+             background=kTRUE;
+             break;
+           }
+         }
+         //else if(!fStudyPureBackground)continue;
+         break;
+       }
+       if(background^fStudyPureBackground)continue;
+       if(fStudyd0fromBTrue){
+         if(b1!=0x0&&b2!=0x0){
+           Int_t charge[2]={0,0};
+           if(b1->Charge()==-1)charge[0]=1;
+           else {
+             if(b2->Charge()==-1){
+               //printf("Same charges for prongs \n");
+               continue;
+             }
+             charge[1]=1;
+           }
+           /* if(!pXtrTrue[charge[0]]=b1->Px()){printf("!first MCtrack:Get momentum X failed \n");continue;}
+              if(!b1->Py(pYtrTrue[charge[0]])){printf("!first MCtrack:Get momentum Y failed \n");continue;}
+              if(!b1->Pz(pZtrTrue[charge[0]])){printf("!first MCtrack:Get momentum Z failed \n");continue;}
+           */
+           
+           pXtrTrue[charge[0]]=b1->Px();
+           pYtrTrue[charge[0]]=b1->Py();
+           pZtrTrue[charge[0]]=b1->Pz();
+           if(!b1->XvYvZv(xtr1)){
+             //printf("!first MCtrack:Get position failed \n");
+             continue;
+           }
+           /*if(!b2->Px(pXtrTrue[charge[1]])){printf("!second MCtrack:Get momentum X failed \n");continue;}
+             if(!b2->Py(pYtrTrue[charge[1]])){printf("!second MCtrack:Get momentum Y failed \n");continue;}
+             if(!b2->Pz(pZtrTrue[charge[1]])){printf("!second MCtrack:Get momentum Z failed \n");continue;}*/
+           pXtrTrue[charge[1]]=b2->Px();
+           pYtrTrue[charge[1]]=b2->Py();
+           pZtrTrue[charge[1]]=b2->Pz();
+           if(!nomum){
+             if(!(mum1==0x0||mum2==0x0)){
+               if(!mum1->PxPyPz(pD)){
+                 //printf("!D from B:Get momentum failed \n");
+                 continue;
+               }
+               if(!mum1->XvYvZv(xD)){
+                 //printf("!D from B:Get position failed \n");
+                 continue;
+               }
+               if(pXtrTrue[0]+pXtrTrue[1]!=pD[0]){
+                 //printf("Warning! MCinfo: Pt of the mother doesn't match the sum of the daughters pt: X component  %f + %f = %f Vs %f \n",pXtrTrue[0],pXtrTrue[1],pXtrTrue[0]+pXtrTrue[1],pD[0]);
+                 //printf("Warning! MCinfo: Pt of the mother doesn't match the sum of the daughters pt: Y componente %f + %f = %f Vs %f \n",pYtrTrue[0],pYtrTrue[1],pYtrTrue[0]+pYtrTrue[1],pD[1]);
+                 //printf("Warning!: Pt comparison from the sum of the daughters:%f Vs mother pT: %f \n",TMath::Sqrt((pXtrTrue[0]+pXtrTrue[1])*(pXtrTrue[0]+pXtrTrue[1])+(pYtrTrue[0]+pYtrTrue[1])*(pYtrTrue[0]+pYtrTrue[1])),TMath::Sqrt(pD[0]*pD[0]+pD[1]*pD[1]));
+                 //          return;
+               }
+               
+               //          if(!b2->PxPyPz(ptr2True)){printf("!second MCtrack:Get momentum failed \n");continue;}
+               if(!b2->XvYvZv(xtr2)){
+                 //printf("!second MCtrack:Get position failed \n");
+                 continue;
+               }
+               Double_t d0dummy[2]={0.,0.};//TEMPORARY : dummy d0 for AliAODRecoDeay constructor
+               aodDMC=new AliAODRecoDecayHF(vtxTrue,xD,2,0,pXtrTrue,pYtrTrue,pZtrTrue,d0dummy);
+               fhMassTrue->Fill(mum1->GetCalcMass());
+             }
+             //else printf("Warning! The candidate comes from two track with a different mother, \n -> the true mother d0 doesn't exist ! -> There will be different entries between observed and MC d0distributions \n");
+           }
+         }
+       }
+       //printf("End of MC checks \n");
+      }
+     
+      //printf("A real D0 was found!! \n");
+      
+    
+      //cout<<1e8*d->Prodd0d0()<<endl;
+      /*      if(fSideBands){  // The integral of the Hists will in general differ 
+                       // from the number of candidates which passed the selection
+       if(okD0)fhMass->Fill(d->InvMassD0(),1.);
+       if(okD0bar)fhMass->Fill(d->InvMassD0bar(),1.);
+      }
+      else { // The integral of the Hists will in general differ 
+       // from the number of candidates which passed the selection
+       if(okD0)fhMass->Fill(d->InvMassD0(),1.);
+       if(okD0bar)fhMass->Fill(d->InvMassD0bar(),1.);
+      }
+      */
+      
+      if(okD0)fhMass->Fill(d->InvMassD0(),1.);
+      if(okD0bar)fhMass->Fill(d->InvMassD0bar(),1.);
+   
+      fhCPtaVSd0d0->Fill(1e8*d->Prodd0d0(),d->CosPointingAngle());
+      fhSecVtxZ->Fill(d->GetSecVtxZ());
+      fhSecVtxX->Fill(d->GetSecVtxX()*10000.);
+      fhSecVtxY->Fill(d->GetSecVtxY()*10000.);
+      fhSecVtxXY->Fill(d->GetSecVtxX()*10000.,d->GetSecVtxY()*10000.);
+      fhSecVtxPhi->Fill(TMath::ATan2(d->GetSecVtxY()*10000.,d->GetSecVtxX()*10000.)*TMath::RadToDeg());
+      fhd0xd0->Fill(1e8*d->Prodd0d0());
+      fhCPta->Fill(d->CosPointingAngle());
+      
+    
+      //cout<<d->GetSecVtxX() <<endl;
+      
+      if(!trk0 || !trk1) {
+       /*      if(d->GetProngID(0)<0||d->GetProngID(1)<0){
+               printf("Wrong ProngID,%d or %d\n",d->GetProngID(0),d->GetProngID(1));continue;
+               }
+       */
+       trk0=fAOD->GetTrack(trkIDtoEntry[d->GetProngID(0)]);
+       trk1=fAOD->GetTrack(trkIDtoEntry[d->GetProngID(1)]);
+       //            cout<<"references to standard AOD not available"<<endl;
+      }
+      //           cout<<"pt of positive track: "<<trk0->Pt()<<endl;
+      
+      // make a AliExternalTrackParam from the D0 
+      // and calculate impact parameters to primary vertex
+      
+      
+      //       AliExternalTrackParam trackD0(d);
+      //cout<<"pt of D0: "<<d->Pt()<<" (AliAODRecoDecay); "<<trackD0.Pt()<<" (track)"<<endl;
+      //trackD0.Print();
+      Double_t ptD0=d->Pt();
+      for(Int_t j=0;j<fnbins;j++){
+       if(ptD0<standardptbins[j+1]){
+         fhd0D0VtxTruept[j]->Fill(d->AliAODRecoDecay::ImpParXY(vtxTrue)*10000.);
+         fhd0D0pt[j]->Fill(d->ImpParXY()*10000.);
+         // printf("HERE in tha class \n");
+         if(aodDMC!=0x0)fhMCd0D0pt[j]->Fill(aodDMC->ImpParXY()*10000.);
+         // printf("HERE2 in tha class \n");
+         break;
+       }
+      }
+      fhd0D0->Fill(d->ImpParXY()*10000.);         
+      fhd0D0VtxTrue->Fill(d->AliAODRecoDecay::ImpParXY(vtxTrue)*10000.);          
+      if(aodDMC!=0x0)fhMCd0D0->Fill(aodDMC->ImpParXY()*10000.);
+      //  if((1.864-3.*0.011.<d->InvMassD0()&&d->InvMassD0()<1.864+3.*0.011)||(1.864-3.*0.011<d->InvMassD0bar()&&d->InvMassD0bar()<1.864+3.*0.011)){
+      if(aodDMC!=0x0){
+       delete aodDMC;
+       aodDMC=0x0;
+      }
+      mum1=0x0;
+      mum2=0x0;
+      b1=0x0;
+      b2=0x0;
+      sonpart=0x0;
+      grandmoth1=0x0;
+    
+      //           }
+      //       Double_t dz[2],covdz[3];
+      //trackD0.PropagateToDCA(vtx1,fAOD->GetMagneticField(),1000.,dz,covdz);
+      //       cout<<"D0 impact parameter rphi: "<<dz[0]<<" +- "<<TMath::Sqrt(covdz[0])<<endl;
+    }
+    // printf("Before end of the event \n");
+    if(unsetvtx) d->UnsetOwnPrimaryVtx();
+         
+  }
+
+  
+  
+  // Post output data.
+  
+  PostData(0,fhCPtaVSd0d0);
+  PostData(1,fhSecVtxXY);
+  PostData(2,fhd0xd0);
+  PostData(3,fhCPta);
+  PostData(4,fhSecVtxZ);
+  PostData(5,fhSecVtxX);
+  PostData(6,fhSecVtxY);
+  PostData(7,fhSecVtxPhi);
+
+  
+  for(Int_t j=0;j<fnbins;j++){
+    
+    PostData(j+8, fhd0D0pt[j]);
+  }
+  
+  PostData(fnbins+8, fhd0D0);
+
+  for(Int_t j=0;j<fnbins;j++){
+    
+    PostData(j+fnbins+9,fhMCd0D0pt[j]);
+  }
+  PostData(2*fnbins+9, fhMCd0D0);
+
+  for(Int_t j=0;j<fnbins;j++){
+    
+    PostData(j+2*fnbins+10,fhd0D0VtxTruept[j]);
+  }
+  PostData(3*fnbins+10,fhd0D0VtxTrue);
+  PostData(3*fnbins+11,fhMassTrue);
+  PostData(3*fnbins+12,fhMass);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskCharmFraction::Terminate(Option_t *) 
+{
+  // Draw result to the screen
+  // Called once at the end of the query
+
+  fhd0D0 = dynamic_cast<TH1F*> (GetOutputData(fnbins+8));
+  if (!fhd0D0) {
+    Printf("ERROR: fhd0D0 not available");
+    return;
+  }
+   
+  TCanvas *c1 = new TCanvas("AliAnalysisTaskCharmFraction","D0 impact par",10,10,510,510);
+  c1->cd(1)->SetLogy();
+  fhd0D0->DrawCopy("E");
+
+  TCanvas *cd0D0pt=new TCanvas("cd0D0pt","cd0D0pt",0,0,800,700);
+  cd0D0pt->Divide(fnbins/2+fnbins%2,2);
+  for(Int_t j=0;j<fnbins;j++){
+    fhd0D0pt[j] = dynamic_cast<TH1F*> (GetOutputData(8+j));
+    if (!fhd0D0pt[j]) {
+      Printf("ERROR: fhd0D0pt not available");
+      return;
+    }
+    cd0D0pt->cd(j+1);
+    fhd0D0pt[j]->Draw();
+  }
+  
+}
diff --git a/PWG3/vertexingHF/AliAnalysisTaskCharmFraction.h b/PWG3/vertexingHF/AliAnalysisTaskCharmFraction.h
new file mode 100644 (file)
index 0000000..224a4f6
--- /dev/null
@@ -0,0 +1,81 @@
+#ifndef ALIANALYSISTASKCHARMFRACTION_H
+#define ALIANALYSISTASKCHARMFRACTION_H
+
+/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//*************************************************************************
+// Class AliAnalysisTaskCharmFraction
+// AliAnalysisTask for the extraction of the fraction of prompt charm
+// using the charm hadron impact parameter to the primary vertex
+// Author: Andrea Rossi andrea.rossi@ts.infn.it
+//*************************************************************************
+
+class TH1F;
+class TH2F;
+class AliAODDEvent;
+class AliAODMCHeader;
+
+#include "AliAnalysisTask.h"
+
+class AliAnalysisTaskCharmFraction : public AliAnalysisTask {
+ public:
+  AliAnalysisTaskCharmFraction(const char *name="AliAnalysisTaskCharmFraction");
+  AliAnalysisTaskCharmFraction(const char *name,Int_t nptbins);
+  
+  virtual ~AliAnalysisTaskCharmFraction() {}
+  
+  virtual void   ConnectInputData(Option_t *);
+  virtual void   CreateOutputObjects();
+  virtual void   Exec(Option_t *option);
+  virtual void   Terminate(Option_t *);
+  void SetNPtBins(Int_t nbins=10){fnbins=nbins;}
+  void SetCheckMC(Bool_t checkMC){fcheckMC=checkMC;}
+  void SetCheckMC_D0(Bool_t check_D0){fcheckMCD0=check_D0;}
+  void SetCheckMC_2prongs(Bool_t check2prongs){fcheckMC2prongs=check2prongs;}
+  void SetCheckMC_prompt(Bool_t checkprompt){fcheckMCprompt=checkprompt;}
+  void SetCheckMC_fromB(Bool_t checkfromB){fcheckMCfromB=checkfromB;}
+  void SetCheckMC_fromDstar(Bool_t skipD0star){fSkipD0star=skipD0star;}
+  void SetUseCuts(Bool_t usecuts){fD0usecuts=usecuts;}
+  void SetSideBands(Double_t sideband){fSideBands=sideband;}
+  void SetStudyPureBackground(Bool_t back){fStudyPureBackground=back;}
+ private:
+  AliAODEvent *fAOD;    //AOD object
+  AliAnalysisVertexingHF *fVHF;        // Vertexer heavy flavour
+  TClonesArray *fArrayD0toKpi;        // Array of D0->Kpi
+  TClonesArray *fArrayMC;             // Array of MC info
+  AliAODMCHeader *fAODmcHeader;        //AOD header
+  TH1F *fhMass;                         //!Inv Mass
+  TH1F *fhMassTrue;                     //!MC inv Mass
+  TH2F *fhCPtaVSd0d0;                  //! histo of CosPtAngle Vs d0xd0
+  TH1F *fhd0xd0;                         //! histo of d0d0
+  TH1F *fhCPta;                        //! histo of cosptangle
+  TH1F *fhSecVtxZ;                     //! histo of z coord of sec. vtx 
+  TH1F *fhSecVtxX;                      //! histo of x coord of sec. vtx
+  TH1F *fhSecVtxY;                      //! histo of x coord of sec. vtx
+  TH2F *fhSecVtxXY;                      //! histo of x coord of sec. vtx
+  TH1F *fhSecVtxPhi;                      //! histo of x coord of sec. vtx
+  TH1F        *fhd0D0;                 //!D0 impact par histo all bins
+  TH1F        **fhd0D0pt;                 //!D0 impact par histos per pt bin 
+  TH1F *fhd0D0VtxTrue;                    //! histo for D0 impact par w.r.t. true vertex, all bins integrated
+  TH1F      **fhd0D0VtxTruept;            //! histos for D0 impact par w.r.t. true vertex
+  TH1F   *fhMCd0D0;                      //!D0 MC impact par histo all bins
+  TH1F      **fhMCd0D0pt;              //!D0 MC impact par histos per pt bin 
+  Int_t        fnbins;                //Number of pt bins
+  Bool_t       fD0usecuts;            // Switch in the use of the cuts
+  Bool_t       fcheckMC;              //  Switch on MC check: minimum check is same mother
+  Bool_t       fcheckMCD0;           //  check the mother is a D0
+  Bool_t       fcheckMC2prongs;         //  check the decay is in two prongs
+  Bool_t       fcheckMCprompt;       //  check the D0 comes from a c quark
+  Bool_t       fcheckMCfromB;        //  check the D0 comes from a b quark
+  Bool_t       fSkipD0star;           // skip if D0 comes from a D*  
+  Bool_t  fStudyd0fromBTrue;         // Flag for analyze true impact par of D0 from B
+  Bool_t  fStudyPureBackground;      // Flag to study the background (reverse the selection on the signal)
+  Double_t  fSideBands;                //Side bands selection (see cxx)
+  AliAnalysisTaskCharmFraction(const AliAnalysisTaskCharmFraction&); // not implemented
+  AliAnalysisTaskCharmFraction& operator=(const AliAnalysisTaskCharmFraction&); // not implemented
+  
+  ClassDef(AliAnalysisTaskCharmFraction,1); // analysis task for prompt charm fraction
+};
+
+#endif